src/jalview/analysis/Dna.java

   1 /*
   2  * Jalview - A Sequence Alignment Editor and Viewer ($$Version-Rel$$)
   3  * Copyright (C) $$Year-Rel$$ The Jalview Authors
   4  *
   5  * This file is part of Jalview.
   6  *
   7  * Jalview is free software: you can redistribute it and/or
   8  * modify it under the terms of the GNU General Public License
   9  * as published by the Free Software Foundation, either version 3
  10  * of the License, or (at your option) any later version.
  11  *
  12  * Jalview is distributed in the hope that it will be useful, but
  13  * WITHOUT ANY WARRANTY; without even the implied warranty
  14  * of MERCHANTABILITY or FITNESS FOR A PARTICULAR
  15  * PURPOSE.  See the GNU General Public License for more details.
  16  *
  17  * You should have received a copy of the GNU General Public License
  18  * along with Jalview.  If not, see <http://www.gnu.org/licenses/>.
  19  * The Jalview Authors are detailed in the 'AUTHORS' file.
  20  */
  21 package jalview.analysis;
  22
  23 import jalview.api.AlignViewportI;
  24 import jalview.datamodel.AlignedCodon;
  25 import jalview.datamodel.AlignedCodonFrame;
  26 import jalview.datamodel.Alignment;
  27 import jalview.datamodel.AlignmentAnnotation;
  28 import jalview.datamodel.AlignmentI;
  29 import jalview.datamodel.Annotation;
  30 import jalview.datamodel.DBRefEntry;
  31 import jalview.datamodel.DBRefSource;
  32 import jalview.datamodel.FeatureProperties;
  33 import jalview.datamodel.GraphLine;
  34 import jalview.datamodel.Mapping;
  35 import jalview.datamodel.Sequence;
  36 import jalview.datamodel.SequenceFeature;
  37 import jalview.datamodel.SequenceI;
  38 import jalview.schemes.ResidueProperties;
  39 import jalview.util.Comparison;
  40 import jalview.util.DBRefUtils;
  41 import jalview.util.MapList;
  42 import jalview.util.ShiftList;
  43
  44 import java.util.ArrayList;
  45 import java.util.Arrays;
  46 import java.util.Comparator;
  47 import java.util.Iterator;
  48 import java.util.List;
  49
  50 public class Dna
  51 {
  52   private static final String STOP_ASTERIX = "*";
  53
  54   private static final Comparator<AlignedCodon> comparator = new CodonComparator();
  55
  56   /*
  57    * 'final' variables describe the inputs to the translation, which should not
  58    * be modified.
  59    */
  60   private final List<SequenceI> selection;
  61
  62   private final String[] seqstring;
  63
  64   private final Iterator<int[]> contigs;
  65
  66   private final char gapChar;
  67
  68   private final AlignmentAnnotation[] annotations;
  69
  70   private final int dnaWidth;
  71
  72   private final AlignmentI dataset;
  73
  74   private ShiftList vismapping;
  75
  76   private int[] startcontigs;
  77
  78   /*
  79    * Working variables for the translation.
  80    *
  81    * The width of the translation-in-progress protein alignment.
  82    */
  83   private int aaWidth = 0;
  84
  85   /*
  86    * This array will be built up so that position i holds the codon positions
  87    * e.g. [7, 9, 10] that match column i (base 0) in the aligned translation.
  88    * Note this implies a contract that if two codons do not align exactly, their
  89    * translated products must occupy different column positions.
  90    */
  91   private AlignedCodon[] alignedCodons;
  92
  93   /**
  94    * Constructor given a viewport and the visible contigs.
  95    *
  96    * @param viewport
  97    * @param visibleContigs
  98    */
  99   public Dna(AlignViewportI viewport, Iterator<int[]> visibleContigs)
 100   {
 101     this.selection = Arrays.asList(viewport.getSequenceSelection());
 102     this.seqstring = viewport.getViewAsString(true);
 103     this.contigs = visibleContigs;
 104     this.gapChar = viewport.getGapCharacter();
 105     this.annotations = viewport.getAlignment().getAlignmentAnnotation();
 106     this.dnaWidth = viewport.getAlignment().getWidth();
 107     this.dataset = viewport.getAlignment().getDataset();
 108     initContigs();
 109   }
 110
 111   /**
 112    * Initialise contigs used as starting point for translateCodingRegion
 113    */
 114   private void initContigs()
 115   {
 116     vismapping = new ShiftList(); // map from viscontigs to seqstring
 117     // intervals
 118
 119     int npos = 0;
 120     int[] lastregion = null;
 121     ArrayList<Integer> tempcontigs = new ArrayList<>();
 122     while (contigs.hasNext())
 123     {
 124       int[] region = contigs.next();
 125       if (lastregion == null)
 126       {
 127         vismapping.addShift(npos, region[0]);
 128       }
 129       else
 130       {
 131         // hidden region
 132         vismapping.addShift(npos, region[0] - lastregion[1] + 1);
 133       }
 134       lastregion = region;
 135       tempcontigs.add(region[0]);
 136       tempcontigs.add(region[1]);
 137     }
 138
 139     startcontigs = new int[tempcontigs.size()];
 140     int i = 0;
 141     for (Integer val : tempcontigs)
 142     {
 143       startcontigs[i] = val;
 144       i++;
 145     }
 146     tempcontigs = null;
 147   }
 148
 149   /**
 150    * Test whether codon positions cdp1 should align before, with, or after cdp2.
 151    * Returns zero if all positions match (or either argument is null). Returns
 152    * -1 if any position in the first codon precedes the corresponding position
 153    * in the second codon. Else returns +1 (some position in the second codon
 154    * precedes the corresponding position in the first).
 155    *
 156    * Note this is not necessarily symmetric, for example:
 157    * <ul>
 158    * <li>compareCodonPos([2,5,6], [3,4,5]) returns -1</li>
 159    * <li>compareCodonPos([3,4,5], [2,5,6]) also returns -1</li>
 160    * </ul>
 161    *
 162    * @param ac1
 163    * @param ac2
 164    * @return
 165    */
 166   public static final int compareCodonPos(AlignedCodon ac1, AlignedCodon ac2)
 167   {
 168     return comparator.compare(ac1, ac2);
 169     // return jalview_2_8_2compare(ac1, ac2);
 170   }
 171
 172   /**
 173    * Codon comparison up to Jalview 2.8.2. This rule is sequence order dependent
 174    * - see http://issues.jalview.org/browse/JAL-1635
 175    *
 176    * @param ac1
 177    * @param ac2
 178    * @return
 179    */
 180   private static int jalview_2_8_2compare(AlignedCodon ac1,
 181           AlignedCodon ac2)
 182   {
 183     if (ac1 == null || ac2 == null || (ac1.equals(ac2)))
 184     {
 185       return 0;
 186     }
 187     if (ac1.pos1 < ac2.pos1 || ac1.pos2 < ac2.pos2 || ac1.pos3 < ac2.pos3)
 188     {
 189       // one base in cdp1 precedes the corresponding base in the other codon
 190       return -1;
 191     }
 192     // one base in cdp1 appears after the corresponding base in the other codon.
 193     return 1;
 194   }
 195
 196   /**
 197    *
 198    * @return
 199    */
 200   public AlignmentI translateCdna()
 201   {
 202     AlignedCodonFrame acf = new AlignedCodonFrame();
 203
 204     alignedCodons = new AlignedCodon[dnaWidth];
 205
 206     int s;
 207     int sSize = selection.size();
 208     List<SequenceI> pepseqs = new ArrayList<>();
 209     for (s = 0; s < sSize; s++)
 210     {
 211       SequenceI newseq = translateCodingRegion(selection.get(s),
 212               seqstring[s], acf, pepseqs);
 213
 214       if (newseq != null)
 215       {
 216         pepseqs.add(newseq);
 217         SequenceI ds = newseq;
 218         if (dataset != null)
 219         {
 220           while (ds.getDatasetSequence() != null)
 221           {
 222             ds = ds.getDatasetSequence();
 223           }
 224           dataset.addSequence(ds);
 225         }
 226       }
 227     }
 228
 229     SequenceI[] newseqs = pepseqs.toArray(new SequenceI[pepseqs.size()]);
 230     AlignmentI al = new Alignment(newseqs);
 231     // ensure we look aligned.
 232     al.padGaps();
 233     // link the protein translation to the DNA dataset
 234     al.setDataset(dataset);
 235     translateAlignedAnnotations(al, acf);
 236     al.addCodonFrame(acf);
 237     return al;
 238   }
 239
 240   /**
 241    * fake the collection of DbRefs with associated exon mappings to identify if
 242    * a translation would generate distinct product in the currently selected
 243    * region.
 244    *
 245    * @param selection
 246    * @param viscontigs
 247    * @return
 248    */
 249   public static boolean canTranslate(SequenceI[] selection,
 250           int viscontigs[])
 251   {
 252     for (int gd = 0; gd < selection.length; gd++)
 253     {
 254       SequenceI dna = selection[gd];
 255       DBRefEntry[] dnarefs = DBRefUtils.selectRefs(dna.getDBRefs(),
 256               jalview.datamodel.DBRefSource.DNACODINGDBS);
 257       if (dnarefs != null)
 258       {
 259         // intersect with pep
 260         List<DBRefEntry> mappedrefs = new ArrayList<>();
 261         DBRefEntry[] refs = dna.getDBRefs();
 262         for (int d = 0; d < refs.length; d++)
 263         {
 264           if (refs[d].getMap() != null && refs[d].getMap().getMap() != null
 265                   && refs[d].getMap().getMap().getFromRatio() == 3
 266                   && refs[d].getMap().getMap().getToRatio() == 1)
 267           {
 268             mappedrefs.add(refs[d]); // add translated protein maps
 269           }
 270         }
 271         dnarefs = mappedrefs.toArray(new DBRefEntry[mappedrefs.size()]);
 272         for (int d = 0; d < dnarefs.length; d++)
 273         {
 274           Mapping mp = dnarefs[d].getMap();
 275           if (mp != null)
 276           {
 277             for (int vc = 0; vc < viscontigs.length; vc += 2)
 278             {
 279               int[] mpr = mp.locateMappedRange(viscontigs[vc],
 280                       viscontigs[vc + 1]);
 281               if (mpr != null)
 282               {
 283                 return true;
 284               }
 285             }
 286           }
 287         }
 288       }
 289     }
 290     return false;
 291   }
 292
 293   /**
 294    * Translate nucleotide alignment annotations onto translated amino acid
 295    * alignment using codon mapping codons
 296    *
 297    * @param al
 298    *          the translated protein alignment
 299    */
 300   protected void translateAlignedAnnotations(AlignmentI al,
 301           AlignedCodonFrame acf)
 302   {
 303     // Can only do this for columns with consecutive codons, or where
 304     // annotation is sequence associated.
 305
 306     if (annotations != null)
 307     {
 308       for (AlignmentAnnotation annotation : annotations)
 309       {
 310         /*
 311          * Skip hidden or autogenerated annotation. Also (for now), RNA
 312          * secondary structure annotation. If we want to show this against
 313          * protein we need a smarter way to 'translate' without generating
 314          * invalid (unbalanced) structure annotation.
 315          */
 316         if (annotation.autoCalculated || !annotation.visible
 317                 || annotation.isRNA())
 318         {
 319           continue;
 320         }
 321
 322         int aSize = aaWidth;
 323         Annotation[] anots = (annotation.annotations == null) ? null
 324                 : new Annotation[aSize];
 325         if (anots != null)
 326         {
 327           for (int a = 0; a < aSize; a++)
 328           {
 329             // process through codon map.
 330             if (a < alignedCodons.length && alignedCodons[a] != null
 331                     && alignedCodons[a].pos1 == (alignedCodons[a].pos3 - 2))
 332             {
 333               anots[a] = getCodonAnnotation(alignedCodons[a],
 334                       annotation.annotations);
 335             }
 336           }
 337         }
 338
 339         AlignmentAnnotation aa = new AlignmentAnnotation(annotation.label,
 340                 annotation.description, anots);
 341         aa.graph = annotation.graph;
 342         aa.graphGroup = annotation.graphGroup;
 343         aa.graphHeight = annotation.graphHeight;
 344         if (annotation.getThreshold() != null)
 345         {
 346           aa.setThreshold(new GraphLine(annotation.getThreshold()));
 347         }
 348         if (annotation.hasScore)
 349         {
 350           aa.setScore(annotation.getScore());
 351         }
 352
 353         final SequenceI seqRef = annotation.sequenceRef;
 354         if (seqRef != null)
 355         {
 356           SequenceI aaSeq = acf.getAaForDnaSeq(seqRef);
 357           if (aaSeq != null)
 358           {
 359             // aa.compactAnnotationArray(); // throw away alignment annotation
 360             // positioning
 361             aa.setSequenceRef(aaSeq);
 362             // rebuild mapping
 363             aa.createSequenceMapping(aaSeq, aaSeq.getStart(), true);
 364             aa.adjustForAlignment();
 365             aaSeq.addAlignmentAnnotation(aa);
 366           }
 367         }
 368         al.addAnnotation(aa);
 369       }
 370     }
 371   }
 372
 373   private static Annotation getCodonAnnotation(AlignedCodon is,
 374           Annotation[] annotations)
 375   {
 376     // Have a look at all the codon positions for annotation and put the first
 377     // one found into the translated annotation pos.
 378     int contrib = 0;
 379     Annotation annot = null;
 380     for (int p = 1; p <= 3; p++)
 381     {
 382       int dnaCol = is.getBaseColumn(p);
 383       if (annotations[dnaCol] != null)
 384       {
 385         if (annot == null)
 386         {
 387           annot = new Annotation(annotations[dnaCol]);
 388           contrib = 1;
 389         }
 390         else
 391         {
 392           // merge with last
 393           Annotation cpy = new Annotation(annotations[dnaCol]);
 394           if (annot.colour == null)
 395           {
 396             annot.colour = cpy.colour;
 397           }
 398           if (annot.description == null || annot.description.length() == 0)
 399           {
 400             annot.description = cpy.description;
 401           }
 402           if (annot.displayCharacter == null)
 403           {
 404             annot.displayCharacter = cpy.displayCharacter;
 405           }
 406           if (annot.secondaryStructure == 0)
 407           {
 408             annot.secondaryStructure = cpy.secondaryStructure;
 409           }
 410           annot.value += cpy.value;
 411           contrib++;
 412         }
 413       }
 414     }
 415     if (contrib > 1)
 416     {
 417       annot.value /= contrib;
 418     }
 419     return annot;
 420   }
 421
 422   /**
 423    * Translate a na sequence
 424    *
 425    * @param selection
 426    *          sequence displayed under viscontigs visible columns
 427    * @param seqstring
 428    *          ORF read in some global alignment reference frame
 429    * @param acf
 430    *          Definition of global ORF alignment reference frame
 431    * @param proteinSeqs
 432    * @return sequence ready to be added to alignment.
 433    */
 434   protected SequenceI translateCodingRegion(SequenceI selection,
 435           String seqstring, AlignedCodonFrame acf,
 436           List<SequenceI> proteinSeqs)
 437   {
 438     List<int[]> skip = new ArrayList<>();
 439     int[] skipint = null;
 440     int npos = 0;
 441     int vc = 0;
 442
 443     int[] scontigs = new int[startcontigs.length];
 444     System.arraycopy(startcontigs, 0, scontigs, 0, startcontigs.length);
 445
 446     // allocate a roughly sized buffer for the protein sequence
 447     StringBuilder protein = new StringBuilder(seqstring.length() / 2);
 448     String seq = seqstring.replace('U', 'T').replace('u', 'T');
 449     char codon[] = new char[3];
 450     int cdp[] = new int[3];
 451     int rf = 0;
 452     int lastnpos = 0;
 453     int nend;
 454     int aspos = 0;
 455     int resSize = 0;
 456     for (npos = 0, nend = seq.length(); npos < nend; npos++)
 457     {
 458       if (!Comparison.isGap(seq.charAt(npos)))
 459       {
 460         cdp[rf] = npos; // store position
 461         codon[rf++] = seq.charAt(npos); // store base
 462       }
 463       if (rf == 3)
 464       {
 465         /*
 466          * Filled up a reading frame...
 467          */
 468         AlignedCodon alignedCodon = new AlignedCodon(cdp[0], cdp[1],
 469                 cdp[2]);
 470         String aa = ResidueProperties.codonTranslate(new String(codon));
 471         rf = 0;
 472         final String gapString = String.valueOf(gapChar);
 473         if (aa == null)
 474         {
 475           aa = gapString;
 476           if (skipint == null)
 477           {
 478             skipint = new int[] { alignedCodon.pos1,
 479                 alignedCodon.pos3 /*
 480                                    * cdp[0],
 481                                    * cdp[2]
 482                                    */ };
 483           }
 484           skipint[1] = alignedCodon.pos3; // cdp[2];
 485         }
 486         else
 487         {
 488           if (skipint != null)
 489           {
 490             // edit scontigs
 491             skipint[0] = vismapping.shift(skipint[0]);
 492             skipint[1] = vismapping.shift(skipint[1]);
 493             for (vc = 0; vc < scontigs.length;)
 494             {
 495               if (scontigs[vc + 1] < skipint[0])
 496               {
 497                 // before skipint starts
 498                 vc += 2;
 499                 continue;
 500               }
 501               if (scontigs[vc] > skipint[1])
 502               {
 503                 // finished editing so
 504                 break;
 505               }
 506               // Edit the contig list to include the skipped region which did
 507               // not translate
 508               int[] t;
 509               // from : s1 e1 s2 e2 s3 e3
 510               // to s: s1 e1 s2 k0 k1 e2 s3 e3
 511               // list increases by one unless one boundary (s2==k0 or e2==k1)
 512               // matches, and decreases by one if skipint intersects whole
 513               // visible contig
 514               if (scontigs[vc] <= skipint[0])
 515               {
 516                 if (skipint[0] == scontigs[vc])
 517                 {
 518                   // skipint at start of contig
 519                   // shift the start of this contig
 520                   if (scontigs[vc + 1] > skipint[1])
 521                   {
 522                     scontigs[vc] = skipint[1];
 523                     vc += 2;
 524                   }
 525                   else
 526                   {
 527                     if (scontigs[vc + 1] == skipint[1])
 528                     {
 529                       // remove the contig
 530                       t = new int[scontigs.length - 2];
 531                       if (vc > 0)
 532                       {
 533                         System.arraycopy(scontigs, 0, t, 0, vc - 1);
 534                       }
 535                       if (vc + 2 < t.length)
 536                       {
 537                         System.arraycopy(scontigs, vc + 2, t, vc,
 538                                 t.length - vc + 2);
 539                       }
 540                       scontigs = t;
 541                     }
 542                     else
 543                     {
 544                       // truncate contig to before the skipint region
 545                       scontigs[vc + 1] = skipint[0] - 1;
 546                       vc += 2;
 547                     }
 548                   }
 549                 }
 550                 else
 551                 {
 552                   // scontig starts before start of skipint
 553                   if (scontigs[vc + 1] < skipint[1])
 554                   {
 555                     // skipint truncates end of scontig
 556                     scontigs[vc + 1] = skipint[0] - 1;
 557                     vc += 2;
 558                   }
 559                   else
 560                   {
 561                     // divide region to new contigs
 562                     t = new int[scontigs.length + 2];
 563                     System.arraycopy(scontigs, 0, t, 0, vc + 1);
 564                     t[vc + 1] = skipint[0];
 565                     t[vc + 2] = skipint[1];
 566                     System.arraycopy(scontigs, vc + 1, t, vc + 3,
 567                             scontigs.length - (vc + 1));
 568                     scontigs = t;
 569                     vc += 4;
 570                   }
 571                 }
 572               }
 573             }
 574             skip.add(skipint);
 575             skipint = null;
 576           }
 577           if (aa.equals("STOP"))
 578           {
 579             aa = STOP_ASTERIX;
 580           }
 581           resSize++;
 582         }
 583         boolean findpos = true;
 584         while (findpos)
 585         {
 586           /*
 587            * Compare this codon's base positions with those currently aligned to
 588            * this column in the translation.
 589            */
 590           final int compareCodonPos = compareCodonPos(alignedCodon,
 591                   alignedCodons[aspos]);
 592           switch (compareCodonPos)
 593           {
 594           case -1:
 595
 596             /*
 597              * This codon should precede the mapped positions - need to insert a
 598              * gap in all prior sequences.
 599              */
 600             insertAAGap(aspos, proteinSeqs);
 601             findpos = false;
 602             break;
 603
 604           case +1:
 605
 606             /*
 607              * This codon belongs after the aligned codons at aspos. Prefix it
 608              * with a gap and try the next position.
 609              */
 610             aa = gapString + aa;
 611             aspos++;
 612             break;
 613
 614           case 0:
 615
 616             /*
 617              * Exact match - codon 'belongs' at this translated position.
 618              */
 619             findpos = false;
 620           }
 621         }
 622         protein.append(aa);
 623         lastnpos = npos;
 624         if (alignedCodons[aspos] == null)
 625         {
 626           // mark this column as aligning to this aligned reading frame
 627           alignedCodons[aspos] = alignedCodon;
 628         }
 629         else if (!alignedCodons[aspos].equals(alignedCodon))
 630         {
 631           throw new IllegalStateException(
 632                   "Tried to coalign " + alignedCodons[aspos].toString()
 633                           + " with " + alignedCodon.toString());
 634         }
 635         if (aspos >= aaWidth)
 636         {
 637           // update maximum alignment width
 638           aaWidth = aspos;
 639         }
 640         // ready for next translated reading frame alignment position (if any)
 641         aspos++;
 642       }
 643     }
 644     if (resSize > 0)
 645     {
 646       SequenceI newseq = new Sequence(selection.getName(),
 647               protein.toString());
 648       if (rf != 0)
 649       {
 650         final String errMsg = "trimming contigs for incomplete terminal codon.";
 651         System.err.println(errMsg);
 652         // map and trim contigs to ORF region
 653         vc = scontigs.length - 1;
 654         lastnpos = vismapping.shift(lastnpos); // place npos in context of
 655         // whole dna alignment (rather
 656         // than visible contigs)
 657         // incomplete ORF could be broken over one or two visible contig
 658         // intervals.
 659         while (vc >= 0 && scontigs[vc] > lastnpos)
 660         {
 661           if (vc > 0 && scontigs[vc - 1] > lastnpos)
 662           {
 663             vc -= 2;
 664           }
 665           else
 666           {
 667             // correct last interval in list.
 668             scontigs[vc] = lastnpos;
 669           }
 670         }
 671
 672         if (vc > 0 && (vc + 1) < scontigs.length)
 673         {
 674           // truncate map list to just vc elements
 675           int t[] = new int[vc + 1];
 676           System.arraycopy(scontigs, 0, t, 0, vc + 1);
 677           scontigs = t;
 678         }
 679         if (vc <= 0)
 680         {
 681           scontigs = null;
 682         }
 683       }
 684       if (scontigs != null)
 685       {
 686         npos = 0;
 687         // map scontigs to actual sequence positions on selection
 688         for (vc = 0; vc < scontigs.length; vc += 2)
 689         {
 690           scontigs[vc] = selection.findPosition(scontigs[vc]); // not from 1!
 691           scontigs[vc + 1] = selection.findPosition(scontigs[vc + 1]); // exclusive
 692           if (scontigs[vc + 1] == selection.getEnd())
 693           {
 694             break;
 695           }
 696         }
 697         // trim trailing empty intervals.
 698         if ((vc + 2) < scontigs.length)
 699         {
 700           int t[] = new int[vc + 2];
 701           System.arraycopy(scontigs, 0, t, 0, vc + 2);
 702           scontigs = t;
 703         }
 704         /*
 705          * delete intervals in scontigs which are not translated. 1. map skip
 706          * into sequence position intervals 2. truncate existing ranges and add
 707          * new ranges to exclude untranslated regions. if (skip.size()>0) {
 708          * Vector narange = new Vector(); for (vc=0; vc<scontigs.length; vc++) {
 709          * narange.addElement(new int[] {scontigs[vc]}); } int sint=0,iv[]; vc =
 710          * 0; while (sint<skip.size()) { skipint = (int[]) skip.elementAt(sint);
 711          * do { iv = (int[]) narange.elementAt(vc); if (iv[0]>=skipint[0] &&
 712          * iv[0]<=skipint[1]) { if (iv[0]==skipint[0]) { // delete beginning of
 713          * range } else { // truncate range and create new one if necessary iv =
 714          * (int[]) narange.elementAt(vc+1); if (iv[0]<=skipint[1]) { // truncate
 715          * range iv[0] = skipint[1]; } else { } } } else if (iv[0]<skipint[0]) {
 716          * iv = (int[]) narange.elementAt(vc+1); } } while (iv[0]) } }
 717          */
 718         MapList map = new MapList(scontigs, new int[] { 1, resSize }, 3, 1);
 719
 720         transferCodedFeatures(selection, newseq, map);
 721
 722         /*
 723          * Construct a dataset sequence for our new peptide.
 724          */
 725         SequenceI rseq = newseq.deriveSequence();
 726
 727         /*
 728          * Store a mapping (between the dataset sequences for the two
 729          * sequences).
 730          */
 731         // SIDE-EFFECT: acf stores the aligned sequence reseq; to remove!
 732         acf.addMap(selection, rseq, map);
 733         return rseq;
 734       }
 735     }
 736     // register the mapping somehow
 737     //
 738     return null;
 739   }
 740
 741   /**
 742    * Insert a gap into the aligned proteins and the codon mapping array.
 743    *
 744    * @param pos
 745    * @param proteinSeqs
 746    * @return
 747    */
 748   protected void insertAAGap(int pos, List<SequenceI> proteinSeqs)
 749   {
 750     aaWidth++;
 751     for (SequenceI seq : proteinSeqs)
 752     {
 753       seq.insertCharAt(pos, gapChar);
 754     }
 755
 756     checkCodonFrameWidth();
 757     if (pos < aaWidth)
 758     {
 759       aaWidth++;
 760
 761       /*
 762        * Shift from [pos] to the end one to the right, and null out [pos]
 763        */
 764       System.arraycopy(alignedCodons, pos, alignedCodons, pos + 1,
 765               alignedCodons.length - pos - 1);
 766       alignedCodons[pos] = null;
 767     }
 768   }
 769
 770   /**
 771    * Check the codons array can accommodate a single insertion, if not resize
 772    * it.
 773    */
 774   protected void checkCodonFrameWidth()
 775   {
 776     if (alignedCodons[alignedCodons.length - 1] != null)
 777     {
 778       /*
 779        * arraycopy insertion would bump a filled slot off the end, so expand.
 780        */
 781       AlignedCodon[] c = new AlignedCodon[alignedCodons.length + 10];
 782       System.arraycopy(alignedCodons, 0, c, 0, alignedCodons.length);
 783       alignedCodons = c;
 784     }
 785   }
 786
 787   /**
 788    * Given a peptide newly translated from a dna sequence, copy over and set any
 789    * features on the peptide from the DNA.
 790    *
 791    * @param dna
 792    * @param pep
 793    * @param map
 794    */
 795   private static void transferCodedFeatures(SequenceI dna, SequenceI pep,
 796           MapList map)
 797   {
 798     DBRefEntry[] dnarefs = DBRefUtils.selectRefs(dna.getDBRefs(),
 799             DBRefSource.DNACODINGDBS);
 800     if (dnarefs != null)
 801     {
 802       // intersect with pep
 803       for (int d = 0; d < dnarefs.length; d++)
 804       {
 805         Mapping mp = dnarefs[d].getMap();
 806         if (mp != null)
 807         {
 808         }
 809       }
 810     }
 811     for (SequenceFeature sf : dna.getFeatures().getAllFeatures())
 812     {
 813         if (FeatureProperties.isCodingFeature(null, sf.getType()))
 814         {
 815           // if (map.intersectsFrom(sf[f].begin, sf[f].end))
 816           {
 817
 818           }
 819         }
 820     }
 821   }
 822
 823   /**
 824    * Returns an alignment consisting of the reversed (and optionally
 825    * complemented) sequences set in this object's constructor
 826    *
 827    * @param complement
 828    * @return
 829    */
 830   public AlignmentI reverseCdna(boolean complement)
 831   {
 832     int sSize = selection.size();
 833     List<SequenceI> reversed = new ArrayList<>();
 834     for (int s = 0; s < sSize; s++)
 835     {
 836       SequenceI newseq = reverseSequence(selection.get(s).getName(),
 837               seqstring[s], complement);
 838
 839       if (newseq != null)
 840       {
 841         reversed.add(newseq);
 842       }
 843     }
 844
 845     SequenceI[] newseqs = reversed.toArray(new SequenceI[reversed.size()]);
 846     AlignmentI al = new Alignment(newseqs);
 847     ((Alignment) al).createDatasetAlignment();
 848     return al;
 849   }
 850
 851   /**
 852    * Returns a reversed, and optionally complemented, sequence. The new
 853    * sequence's name is the original name with "|rev" or "|revcomp" appended.
 854    * aAcCgGtT and DNA ambiguity codes are complemented, any other characters are
 855    * left unchanged.
 856    *
 857    * @param seq
 858    * @param complement
 859    * @return
 860    */
 861   public static SequenceI reverseSequence(String seqName, String sequence,
 862           boolean complement)
 863   {
 864     String newName = seqName + "|rev" + (complement ? "comp" : "");
 865     char[] originalSequence = sequence.toCharArray();
 866     int length = originalSequence.length;
 867     char[] reversedSequence = new char[length];
 868     int bases = 0;
 869     for (int i = 0; i < length; i++)
 870     {
 871       char c = complement ? getComplement(originalSequence[i])
 872               : originalSequence[i];
 873       reversedSequence[length - i - 1] = c;
 874       if (!Comparison.isGap(c))
 875       {
 876         bases++;
 877       }
 878     }
 879     SequenceI reversed = new Sequence(newName, reversedSequence, 1, bases);
 880     return reversed;
 881   }
 882
 883   /**
 884    * Returns dna complement (preserving case) for aAcCgGtTuU. Ambiguity codes
 885    * are treated as on http://reverse-complement.com/. Anything else is left
 886    * unchanged.
 887    *
 888    * @param c
 889    * @return
 890    */
 891   public static char getComplement(char c)
 892   {
 893     char result = c;
 894     switch (c)
 895     {
 896     case '-':
 897     case '.':
 898     case ' ':
 899       break;
 900     case 'a':
 901       result = 't';
 902       break;
 903     case 'A':
 904       result = 'T';
 905       break;
 906     case 'c':
 907       result = 'g';
 908       break;
 909     case 'C':
 910       result = 'G';
 911       break;
 912     case 'g':
 913       result = 'c';
 914       break;
 915     case 'G':
 916       result = 'C';
 917       break;
 918     case 't':
 919       result = 'a';
 920       break;
 921     case 'T':
 922       result = 'A';
 923       break;
 924     case 'u':
 925       result = 'a';
 926       break;
 927     case 'U':
 928       result = 'A';
 929       break;
 930     case 'r':
 931       result = 'y';
 932       break;
 933     case 'R':
 934       result = 'Y';
 935       break;
 936     case 'y':
 937       result = 'r';
 938       break;
 939     case 'Y':
 940       result = 'R';
 941       break;
 942     case 'k':
 943       result = 'm';
 944       break;
 945     case 'K':
 946       result = 'M';
 947       break;
 948     case 'm':
 949       result = 'k';
 950       break;
 951     case 'M':
 952       result = 'K';
 953       break;
 954     case 'b':
 955       result = 'v';
 956       break;
 957     case 'B':
 958       result = 'V';
 959       break;
 960     case 'v':
 961       result = 'b';
 962       break;
 963     case 'V':
 964       result = 'B';
 965       break;
 966     case 'd':
 967       result = 'h';
 968       break;
 969     case 'D':
 970       result = 'H';
 971       break;
 972     case 'h':
 973       result = 'd';
 974       break;
 975     case 'H':
 976       result = 'D';
 977       break;
 978     }
 979
 980     return result;
 981   }
 982 }