JAL-2759 Moved propagateInsertions out of HiddenColumns
[jalview.git] / 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 }