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