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