JAL-2759 Moved propagateInsertions out of HiddenColumns
[jalview.git] / src / jalview / analysis / AlignmentUtils.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 static jalview.io.gff.GffConstants.CLINICAL_SIGNIFICANCE;
24
25 import jalview.datamodel.AlignedCodon;
26 import jalview.datamodel.AlignedCodonFrame;
27 import jalview.datamodel.AlignedCodonFrame.SequenceToSequenceMapping;
28 import jalview.datamodel.Alignment;
29 import jalview.datamodel.AlignmentAnnotation;
30 import jalview.datamodel.AlignmentI;
31 import jalview.datamodel.DBRefEntry;
32 import jalview.datamodel.IncompleteCodonException;
33 import jalview.datamodel.Mapping;
34 import jalview.datamodel.Sequence;
35 import jalview.datamodel.SequenceFeature;
36 import jalview.datamodel.SequenceGroup;
37 import jalview.datamodel.SequenceI;
38 import jalview.datamodel.features.SequenceFeatures;
39 import jalview.io.gff.SequenceOntologyI;
40 import jalview.schemes.ResidueProperties;
41 import jalview.util.Comparison;
42 import jalview.util.DBRefUtils;
43 import jalview.util.IntRangeComparator;
44 import jalview.util.MapList;
45 import jalview.util.MappingUtils;
46 import jalview.util.StringUtils;
47
48 import java.io.UnsupportedEncodingException;
49 import java.net.URLEncoder;
50 import java.util.ArrayList;
51 import java.util.Arrays;
52 import java.util.Collection;
53 import java.util.Collections;
54 import java.util.HashMap;
55 import java.util.HashSet;
56 import java.util.Iterator;
57 import java.util.LinkedHashMap;
58 import java.util.List;
59 import java.util.Map;
60 import java.util.Map.Entry;
61 import java.util.NoSuchElementException;
62 import java.util.Set;
63 import java.util.SortedMap;
64 import java.util.TreeMap;
65
66 /**
67  * grab bag of useful alignment manipulation operations Expect these to be
68  * refactored elsewhere at some point.
69  * 
70  * @author jimp
71  * 
72  */
73 public class AlignmentUtils
74 {
75
76   private static final int CODON_LENGTH = 3;
77
78   private static final String SEQUENCE_VARIANT = "sequence_variant:";
79
80   private static final String ID = "ID";
81
82   /**
83    * A data model to hold the 'normal' base value at a position, and an optional
84    * sequence variant feature
85    */
86   static final class DnaVariant
87   {
88     final String base;
89
90     SequenceFeature variant;
91
92     DnaVariant(String nuc)
93     {
94       base = nuc;
95       variant = null;
96     }
97
98     DnaVariant(String nuc, SequenceFeature var)
99     {
100       base = nuc;
101       variant = var;
102     }
103
104     public String getSource()
105     {
106       return variant == null ? null : variant.getFeatureGroup();
107     }
108   }
109
110   /**
111    * given an existing alignment, create a new alignment including all, or up to
112    * flankSize additional symbols from each sequence's dataset sequence
113    * 
114    * @param core
115    * @param flankSize
116    * @return AlignmentI
117    */
118   public static AlignmentI expandContext(AlignmentI core, int flankSize)
119   {
120     List<SequenceI> sq = new ArrayList<SequenceI>();
121     int maxoffset = 0;
122     for (SequenceI s : core.getSequences())
123     {
124       SequenceI newSeq = s.deriveSequence();
125       final int newSeqStart = newSeq.getStart() - 1;
126       if (newSeqStart > maxoffset
127               && newSeq.getDatasetSequence().getStart() < s.getStart())
128       {
129         maxoffset = newSeqStart;
130       }
131       sq.add(newSeq);
132     }
133     if (flankSize > -1)
134     {
135       maxoffset = Math.min(maxoffset, flankSize);
136     }
137
138     /*
139      * now add offset left and right to create an expanded alignment
140      */
141     for (SequenceI s : sq)
142     {
143       SequenceI ds = s;
144       while (ds.getDatasetSequence() != null)
145       {
146         ds = ds.getDatasetSequence();
147       }
148       int s_end = s.findPosition(s.getStart() + s.getLength());
149       // find available flanking residues for sequence
150       int ustream_ds = s.getStart() - ds.getStart();
151       int dstream_ds = ds.getEnd() - s_end;
152
153       // build new flanked sequence
154
155       // compute gap padding to start of flanking sequence
156       int offset = maxoffset - ustream_ds;
157
158       // padding is gapChar x ( maxoffset - min(ustream_ds, flank)
159       if (flankSize >= 0)
160       {
161         if (flankSize < ustream_ds)
162         {
163           // take up to flankSize residues
164           offset = maxoffset - flankSize;
165           ustream_ds = flankSize;
166         }
167         if (flankSize <= dstream_ds)
168         {
169           dstream_ds = flankSize - 1;
170         }
171       }
172       // TODO use Character.toLowerCase to avoid creating String objects?
173       char[] upstream = new String(ds
174               .getSequence(s.getStart() - 1 - ustream_ds, s.getStart() - 1))
175                       .toLowerCase().toCharArray();
176       char[] downstream = new String(
177               ds.getSequence(s_end - 1, s_end + dstream_ds)).toLowerCase()
178                       .toCharArray();
179       char[] coreseq = s.getSequence();
180       char[] nseq = new char[offset + upstream.length + downstream.length
181               + coreseq.length];
182       char c = core.getGapCharacter();
183
184       int p = 0;
185       for (; p < offset; p++)
186       {
187         nseq[p] = c;
188       }
189
190       System.arraycopy(upstream, 0, nseq, p, upstream.length);
191       System.arraycopy(coreseq, 0, nseq, p + upstream.length,
192               coreseq.length);
193       System.arraycopy(downstream, 0, nseq,
194               p + coreseq.length + upstream.length, downstream.length);
195       s.setSequence(new String(nseq));
196       s.setStart(s.getStart() - ustream_ds);
197       s.setEnd(s_end + downstream.length);
198     }
199     AlignmentI newAl = new jalview.datamodel.Alignment(
200             sq.toArray(new SequenceI[0]));
201     for (SequenceI s : sq)
202     {
203       if (s.getAnnotation() != null)
204       {
205         for (AlignmentAnnotation aa : s.getAnnotation())
206         {
207           aa.adjustForAlignment(); // JAL-1712 fix
208           newAl.addAnnotation(aa);
209         }
210       }
211     }
212     newAl.setDataset(core.getDataset());
213     return newAl;
214   }
215
216   /**
217    * Returns the index (zero-based position) of a sequence in an alignment, or
218    * -1 if not found.
219    * 
220    * @param al
221    * @param seq
222    * @return
223    */
224   public static int getSequenceIndex(AlignmentI al, SequenceI seq)
225   {
226     int result = -1;
227     int pos = 0;
228     for (SequenceI alSeq : al.getSequences())
229     {
230       if (alSeq == seq)
231       {
232         result = pos;
233         break;
234       }
235       pos++;
236     }
237     return result;
238   }
239
240   /**
241    * Returns a map of lists of sequences in the alignment, keyed by sequence
242    * name. For use in mapping between different alignment views of the same
243    * sequences.
244    * 
245    * @see jalview.datamodel.AlignmentI#getSequencesByName()
246    */
247   public static Map<String, List<SequenceI>> getSequencesByName(
248           AlignmentI al)
249   {
250     Map<String, List<SequenceI>> theMap = new LinkedHashMap<String, List<SequenceI>>();
251     for (SequenceI seq : al.getSequences())
252     {
253       String name = seq.getName();
254       if (name != null)
255       {
256         List<SequenceI> seqs = theMap.get(name);
257         if (seqs == null)
258         {
259           seqs = new ArrayList<SequenceI>();
260           theMap.put(name, seqs);
261         }
262         seqs.add(seq);
263       }
264     }
265     return theMap;
266   }
267
268   /**
269    * Build mapping of protein to cDNA alignment. Mappings are made between
270    * sequences where the cDNA translates to the protein sequence. Any new
271    * mappings are added to the protein alignment. Returns true if any mappings
272    * either already exist or were added, else false.
273    * 
274    * @param proteinAlignment
275    * @param cdnaAlignment
276    * @return
277    */
278   public static boolean mapProteinAlignmentToCdna(
279           final AlignmentI proteinAlignment, final AlignmentI cdnaAlignment)
280   {
281     if (proteinAlignment == null || cdnaAlignment == null)
282     {
283       return false;
284     }
285
286     Set<SequenceI> mappedDna = new HashSet<SequenceI>();
287     Set<SequenceI> mappedProtein = new HashSet<SequenceI>();
288
289     /*
290      * First pass - map sequences where cross-references exist. This include
291      * 1-to-many mappings to support, for example, variant cDNA.
292      */
293     boolean mappingPerformed = mapProteinToCdna(proteinAlignment,
294             cdnaAlignment, mappedDna, mappedProtein, true);
295
296     /*
297      * Second pass - map sequences where no cross-references exist. This only
298      * does 1-to-1 mappings and assumes corresponding sequences are in the same
299      * order in the alignments.
300      */
301     mappingPerformed |= mapProteinToCdna(proteinAlignment, cdnaAlignment,
302             mappedDna, mappedProtein, false);
303     return mappingPerformed;
304   }
305
306   /**
307    * Make mappings between compatible sequences (where the cDNA translation
308    * matches the protein).
309    * 
310    * @param proteinAlignment
311    * @param cdnaAlignment
312    * @param mappedDna
313    *          a set of mapped DNA sequences (to add to)
314    * @param mappedProtein
315    *          a set of mapped Protein sequences (to add to)
316    * @param xrefsOnly
317    *          if true, only map sequences where xrefs exist
318    * @return
319    */
320   protected static boolean mapProteinToCdna(
321           final AlignmentI proteinAlignment, final AlignmentI cdnaAlignment,
322           Set<SequenceI> mappedDna, Set<SequenceI> mappedProtein,
323           boolean xrefsOnly)
324   {
325     boolean mappingExistsOrAdded = false;
326     List<SequenceI> thisSeqs = proteinAlignment.getSequences();
327     for (SequenceI aaSeq : thisSeqs)
328     {
329       boolean proteinMapped = false;
330       AlignedCodonFrame acf = new AlignedCodonFrame();
331
332       for (SequenceI cdnaSeq : cdnaAlignment.getSequences())
333       {
334         /*
335          * Always try to map if sequences have xref to each other; this supports
336          * variant cDNA or alternative splicing for a protein sequence.
337          * 
338          * If no xrefs, try to map progressively, assuming that alignments have
339          * mappable sequences in corresponding order. These are not
340          * many-to-many, as that would risk mixing species with similar cDNA
341          * sequences.
342          */
343         if (xrefsOnly && !AlignmentUtils.haveCrossRef(aaSeq, cdnaSeq))
344         {
345           continue;
346         }
347
348         /*
349          * Don't map non-xrefd sequences more than once each. This heuristic
350          * allows us to pair up similar sequences in ordered alignments.
351          */
352         if (!xrefsOnly && (mappedProtein.contains(aaSeq)
353                 || mappedDna.contains(cdnaSeq)))
354         {
355           continue;
356         }
357         if (mappingExists(proteinAlignment.getCodonFrames(),
358                 aaSeq.getDatasetSequence(), cdnaSeq.getDatasetSequence()))
359         {
360           mappingExistsOrAdded = true;
361         }
362         else
363         {
364           MapList map = mapCdnaToProtein(aaSeq, cdnaSeq);
365           if (map != null)
366           {
367             acf.addMap(cdnaSeq, aaSeq, map);
368             mappingExistsOrAdded = true;
369             proteinMapped = true;
370             mappedDna.add(cdnaSeq);
371             mappedProtein.add(aaSeq);
372           }
373         }
374       }
375       if (proteinMapped)
376       {
377         proteinAlignment.addCodonFrame(acf);
378       }
379     }
380     return mappingExistsOrAdded;
381   }
382
383   /**
384    * Answers true if the mappings include one between the given (dataset)
385    * sequences.
386    */
387   public static boolean mappingExists(List<AlignedCodonFrame> mappings,
388           SequenceI aaSeq, SequenceI cdnaSeq)
389   {
390     if (mappings != null)
391     {
392       for (AlignedCodonFrame acf : mappings)
393       {
394         if (cdnaSeq == acf.getDnaForAaSeq(aaSeq))
395         {
396           return true;
397         }
398       }
399     }
400     return false;
401   }
402
403   /**
404    * Builds a mapping (if possible) of a cDNA to a protein sequence.
405    * <ul>
406    * <li>first checks if the cdna translates exactly to the protein
407    * sequence</li>
408    * <li>else checks for translation after removing a STOP codon</li>
409    * <li>else checks for translation after removing a START codon</li>
410    * <li>if that fails, inspect CDS features on the cDNA sequence</li>
411    * </ul>
412    * Returns null if no mapping is determined.
413    * 
414    * @param proteinSeq
415    *          the aligned protein sequence
416    * @param cdnaSeq
417    *          the aligned cdna sequence
418    * @return
419    */
420   public static MapList mapCdnaToProtein(SequenceI proteinSeq,
421           SequenceI cdnaSeq)
422   {
423     /*
424      * Here we handle either dataset sequence set (desktop) or absent (applet).
425      * Use only the char[] form of the sequence to avoid creating possibly large
426      * String objects.
427      */
428     final SequenceI proteinDataset = proteinSeq.getDatasetSequence();
429     char[] aaSeqChars = proteinDataset != null
430             ? proteinDataset.getSequence()
431             : proteinSeq.getSequence();
432     final SequenceI cdnaDataset = cdnaSeq.getDatasetSequence();
433     char[] cdnaSeqChars = cdnaDataset != null ? cdnaDataset.getSequence()
434             : cdnaSeq.getSequence();
435     if (aaSeqChars == null || cdnaSeqChars == null)
436     {
437       return null;
438     }
439
440     /*
441      * cdnaStart/End, proteinStartEnd are base 1 (for dataset sequence mapping)
442      */
443     final int mappedLength = CODON_LENGTH * aaSeqChars.length;
444     int cdnaLength = cdnaSeqChars.length;
445     int cdnaStart = cdnaSeq.getStart();
446     int cdnaEnd = cdnaSeq.getEnd();
447     final int proteinStart = proteinSeq.getStart();
448     final int proteinEnd = proteinSeq.getEnd();
449
450     /*
451      * If lengths don't match, try ignoring stop codon (if present)
452      */
453     if (cdnaLength != mappedLength && cdnaLength > 2)
454     {
455       String lastCodon = String.valueOf(cdnaSeqChars,
456               cdnaLength - CODON_LENGTH, CODON_LENGTH).toUpperCase();
457       for (String stop : ResidueProperties.STOP)
458       {
459         if (lastCodon.equals(stop))
460         {
461           cdnaEnd -= CODON_LENGTH;
462           cdnaLength -= CODON_LENGTH;
463           break;
464         }
465       }
466     }
467
468     /*
469      * If lengths still don't match, try ignoring start codon.
470      */
471     int startOffset = 0;
472     if (cdnaLength != mappedLength && cdnaLength > 2
473             && String.valueOf(cdnaSeqChars, 0, CODON_LENGTH).toUpperCase()
474                     .equals(ResidueProperties.START))
475     {
476       startOffset += CODON_LENGTH;
477       cdnaStart += CODON_LENGTH;
478       cdnaLength -= CODON_LENGTH;
479     }
480
481     if (translatesAs(cdnaSeqChars, startOffset, aaSeqChars))
482     {
483       /*
484        * protein is translation of dna (+/- start/stop codons)
485        */
486       MapList map = new MapList(new int[] { cdnaStart, cdnaEnd },
487               new int[]
488               { proteinStart, proteinEnd }, CODON_LENGTH, 1);
489       return map;
490     }
491
492     /*
493      * translation failed - try mapping CDS annotated regions of dna
494      */
495     return mapCdsToProtein(cdnaSeq, proteinSeq);
496   }
497
498   /**
499    * Test whether the given cdna sequence, starting at the given offset,
500    * translates to the given amino acid sequence, using the standard translation
501    * table. Designed to fail fast i.e. as soon as a mismatch position is found.
502    * 
503    * @param cdnaSeqChars
504    * @param cdnaStart
505    * @param aaSeqChars
506    * @return
507    */
508   protected static boolean translatesAs(char[] cdnaSeqChars, int cdnaStart,
509           char[] aaSeqChars)
510   {
511     if (cdnaSeqChars == null || aaSeqChars == null)
512     {
513       return false;
514     }
515
516     int aaPos = 0;
517     int dnaPos = cdnaStart;
518     for (; dnaPos < cdnaSeqChars.length - 2
519             && aaPos < aaSeqChars.length; dnaPos += CODON_LENGTH, aaPos++)
520     {
521       String codon = String.valueOf(cdnaSeqChars, dnaPos, CODON_LENGTH);
522       final String translated = ResidueProperties.codonTranslate(codon);
523
524       /*
525        * allow * in protein to match untranslatable in dna
526        */
527       final char aaRes = aaSeqChars[aaPos];
528       if ((translated == null || "STOP".equals(translated)) && aaRes == '*')
529       {
530         continue;
531       }
532       if (translated == null || !(aaRes == translated.charAt(0)))
533       {
534         // debug
535         // System.out.println(("Mismatch at " + i + "/" + aaResidue + ": "
536         // + codon + "(" + translated + ") != " + aaRes));
537         return false;
538       }
539     }
540
541     /*
542      * check we matched all of the protein sequence
543      */
544     if (aaPos != aaSeqChars.length)
545     {
546       return false;
547     }
548
549     /*
550      * check we matched all of the dna except
551      * for optional trailing STOP codon
552      */
553     if (dnaPos == cdnaSeqChars.length)
554     {
555       return true;
556     }
557     if (dnaPos == cdnaSeqChars.length - CODON_LENGTH)
558     {
559       String codon = String.valueOf(cdnaSeqChars, dnaPos, CODON_LENGTH);
560       if ("STOP".equals(ResidueProperties.codonTranslate(codon)))
561       {
562         return true;
563       }
564     }
565     return false;
566   }
567
568   /**
569    * Align sequence 'seq' to match the alignment of a mapped sequence. Note this
570    * currently assumes that we are aligning cDNA to match protein.
571    * 
572    * @param seq
573    *          the sequence to be realigned
574    * @param al
575    *          the alignment whose sequence alignment is to be 'copied'
576    * @param gap
577    *          character string represent a gap in the realigned sequence
578    * @param preserveUnmappedGaps
579    * @param preserveMappedGaps
580    * @return true if the sequence was realigned, false if it could not be
581    */
582   public static boolean alignSequenceAs(SequenceI seq, AlignmentI al,
583           String gap, boolean preserveMappedGaps,
584           boolean preserveUnmappedGaps)
585   {
586     /*
587      * Get any mappings from the source alignment to the target (dataset)
588      * sequence.
589      */
590     // TODO there may be one AlignedCodonFrame per dataset sequence, or one with
591     // all mappings. Would it help to constrain this?
592     List<AlignedCodonFrame> mappings = al.getCodonFrame(seq);
593     if (mappings == null || mappings.isEmpty())
594     {
595       return false;
596     }
597
598     /*
599      * Locate the aligned source sequence whose dataset sequence is mapped. We
600      * just take the first match here (as we can't align like more than one
601      * sequence).
602      */
603     SequenceI alignFrom = null;
604     AlignedCodonFrame mapping = null;
605     for (AlignedCodonFrame mp : mappings)
606     {
607       alignFrom = mp.findAlignedSequence(seq, al);
608       if (alignFrom != null)
609       {
610         mapping = mp;
611         break;
612       }
613     }
614
615     if (alignFrom == null)
616     {
617       return false;
618     }
619     alignSequenceAs(seq, alignFrom, mapping, gap, al.getGapCharacter(),
620             preserveMappedGaps, preserveUnmappedGaps);
621     return true;
622   }
623
624   /**
625    * Align sequence 'alignTo' the same way as 'alignFrom', using the mapping to
626    * match residues and codons. Flags control whether existing gaps in unmapped
627    * (intron) and mapped (exon) regions are preserved or not. Gaps between
628    * intron and exon are only retained if both flags are set.
629    * 
630    * @param alignTo
631    * @param alignFrom
632    * @param mapping
633    * @param myGap
634    * @param sourceGap
635    * @param preserveUnmappedGaps
636    * @param preserveMappedGaps
637    */
638   public static void alignSequenceAs(SequenceI alignTo, SequenceI alignFrom,
639           AlignedCodonFrame mapping, String myGap, char sourceGap,
640           boolean preserveMappedGaps, boolean preserveUnmappedGaps)
641   {
642     // TODO generalise to work for Protein-Protein, dna-dna, dna-protein
643
644     // aligned and dataset sequence positions, all base zero
645     int thisSeqPos = 0;
646     int sourceDsPos = 0;
647
648     int basesWritten = 0;
649     char myGapChar = myGap.charAt(0);
650     int ratio = myGap.length();
651
652     int fromOffset = alignFrom.getStart() - 1;
653     int toOffset = alignTo.getStart() - 1;
654     int sourceGapMappedLength = 0;
655     boolean inExon = false;
656     final int toLength = alignTo.getLength();
657     final int fromLength = alignFrom.getLength();
658     StringBuilder thisAligned = new StringBuilder(2 * toLength);
659
660     /*
661      * Traverse the 'model' aligned sequence
662      */
663     for (int i = 0; i < fromLength; i++)
664     {
665       char sourceChar = alignFrom.getCharAt(i);
666       if (sourceChar == sourceGap)
667       {
668         sourceGapMappedLength += ratio;
669         continue;
670       }
671
672       /*
673        * Found a non-gap character. Locate its mapped region if any.
674        */
675       sourceDsPos++;
676       // Note mapping positions are base 1, our sequence positions base 0
677       int[] mappedPos = mapping.getMappedRegion(alignTo, alignFrom,
678               sourceDsPos + fromOffset);
679       if (mappedPos == null)
680       {
681         /*
682          * unmapped position; treat like a gap
683          */
684         sourceGapMappedLength += ratio;
685         // System.err.println("Can't align: no codon mapping to residue "
686         // + sourceDsPos + "(" + sourceChar + ")");
687         // return;
688         continue;
689       }
690
691       int mappedCodonStart = mappedPos[0]; // position (1...) of codon start
692       int mappedCodonEnd = mappedPos[mappedPos.length - 1]; // codon end pos
693       StringBuilder trailingCopiedGap = new StringBuilder();
694
695       /*
696        * Copy dna sequence up to and including this codon. Optionally, include
697        * gaps before the codon starts (in introns) and/or after the codon starts
698        * (in exons).
699        * 
700        * Note this only works for 'linear' splicing, not reverse or interleaved.
701        * But then 'align dna as protein' doesn't make much sense otherwise.
702        */
703       int intronLength = 0;
704       while (basesWritten + toOffset < mappedCodonEnd
705               && thisSeqPos < toLength)
706       {
707         final char c = alignTo.getCharAt(thisSeqPos++);
708         if (c != myGapChar)
709         {
710           basesWritten++;
711           int sourcePosition = basesWritten + toOffset;
712           if (sourcePosition < mappedCodonStart)
713           {
714             /*
715              * Found an unmapped (intron) base. First add in any preceding gaps
716              * (if wanted).
717              */
718             if (preserveUnmappedGaps && trailingCopiedGap.length() > 0)
719             {
720               thisAligned.append(trailingCopiedGap.toString());
721               intronLength += trailingCopiedGap.length();
722               trailingCopiedGap = new StringBuilder();
723             }
724             intronLength++;
725             inExon = false;
726           }
727           else
728           {
729             final boolean startOfCodon = sourcePosition == mappedCodonStart;
730             int gapsToAdd = calculateGapsToInsert(preserveMappedGaps,
731                     preserveUnmappedGaps, sourceGapMappedLength, inExon,
732                     trailingCopiedGap.length(), intronLength, startOfCodon);
733             for (int k = 0; k < gapsToAdd; k++)
734             {
735               thisAligned.append(myGapChar);
736             }
737             sourceGapMappedLength = 0;
738             inExon = true;
739           }
740           thisAligned.append(c);
741           trailingCopiedGap = new StringBuilder();
742         }
743         else
744         {
745           if (inExon && preserveMappedGaps)
746           {
747             trailingCopiedGap.append(myGapChar);
748           }
749           else if (!inExon && preserveUnmappedGaps)
750           {
751             trailingCopiedGap.append(myGapChar);
752           }
753         }
754       }
755     }
756
757     /*
758      * At end of model aligned sequence. Copy any remaining target sequence, optionally
759      * including (intron) gaps.
760      */
761     while (thisSeqPos < toLength)
762     {
763       final char c = alignTo.getCharAt(thisSeqPos++);
764       if (c != myGapChar || preserveUnmappedGaps)
765       {
766         thisAligned.append(c);
767       }
768       sourceGapMappedLength--;
769     }
770
771     /*
772      * finally add gaps to pad for any trailing source gaps or
773      * unmapped characters
774      */
775     if (preserveUnmappedGaps)
776     {
777       while (sourceGapMappedLength > 0)
778       {
779         thisAligned.append(myGapChar);
780         sourceGapMappedLength--;
781       }
782     }
783
784     /*
785      * All done aligning, set the aligned sequence.
786      */
787     alignTo.setSequence(new String(thisAligned));
788   }
789
790   /**
791    * Helper method to work out how many gaps to insert when realigning.
792    * 
793    * @param preserveMappedGaps
794    * @param preserveUnmappedGaps
795    * @param sourceGapMappedLength
796    * @param inExon
797    * @param trailingCopiedGap
798    * @param intronLength
799    * @param startOfCodon
800    * @return
801    */
802   protected static int calculateGapsToInsert(boolean preserveMappedGaps,
803           boolean preserveUnmappedGaps, int sourceGapMappedLength,
804           boolean inExon, int trailingGapLength, int intronLength,
805           final boolean startOfCodon)
806   {
807     int gapsToAdd = 0;
808     if (startOfCodon)
809     {
810       /*
811        * Reached start of codon. Ignore trailing gaps in intron unless we are
812        * preserving gaps in both exon and intron. Ignore them anyway if the
813        * protein alignment introduces a gap at least as large as the intronic
814        * region.
815        */
816       if (inExon && !preserveMappedGaps)
817       {
818         trailingGapLength = 0;
819       }
820       if (!inExon && !(preserveMappedGaps && preserveUnmappedGaps))
821       {
822         trailingGapLength = 0;
823       }
824       if (inExon)
825       {
826         gapsToAdd = Math.max(sourceGapMappedLength, trailingGapLength);
827       }
828       else
829       {
830         if (intronLength + trailingGapLength <= sourceGapMappedLength)
831         {
832           gapsToAdd = sourceGapMappedLength - intronLength;
833         }
834         else
835         {
836           gapsToAdd = Math.min(
837                   intronLength + trailingGapLength - sourceGapMappedLength,
838                   trailingGapLength);
839         }
840       }
841     }
842     else
843     {
844       /*
845        * second or third base of codon; check for any gaps in dna
846        */
847       if (!preserveMappedGaps)
848       {
849         trailingGapLength = 0;
850       }
851       gapsToAdd = Math.max(sourceGapMappedLength, trailingGapLength);
852     }
853     return gapsToAdd;
854   }
855
856   /**
857    * Realigns the given protein to match the alignment of the dna, using codon
858    * mappings to translate aligned codon positions to protein residues.
859    * 
860    * @param protein
861    *          the alignment whose sequences are realigned by this method
862    * @param dna
863    *          the dna alignment whose alignment we are 'copying'
864    * @return the number of sequences that were realigned
865    */
866   public static int alignProteinAsDna(AlignmentI protein, AlignmentI dna)
867   {
868     if (protein.isNucleotide() || !dna.isNucleotide())
869     {
870       System.err.println("Wrong alignment type in alignProteinAsDna");
871       return 0;
872     }
873     List<SequenceI> unmappedProtein = new ArrayList<SequenceI>();
874     Map<AlignedCodon, Map<SequenceI, AlignedCodon>> alignedCodons = buildCodonColumnsMap(
875             protein, dna, unmappedProtein);
876     return alignProteinAs(protein, alignedCodons, unmappedProtein);
877   }
878
879   /**
880    * Realigns the given dna to match the alignment of the protein, using codon
881    * mappings to translate aligned peptide positions to codons.
882    * 
883    * Always produces a padded CDS alignment.
884    * 
885    * @param dna
886    *          the alignment whose sequences are realigned by this method
887    * @param protein
888    *          the protein alignment whose alignment we are 'copying'
889    * @return the number of sequences that were realigned
890    */
891   public static int alignCdsAsProtein(AlignmentI dna, AlignmentI protein)
892   {
893     if (protein.isNucleotide() || !dna.isNucleotide())
894     {
895       System.err.println("Wrong alignment type in alignProteinAsDna");
896       return 0;
897     }
898     // todo: implement this
899     List<AlignedCodonFrame> mappings = protein.getCodonFrames();
900     int alignedCount = 0;
901     int width = 0; // alignment width for padding CDS
902     for (SequenceI dnaSeq : dna.getSequences())
903     {
904       if (alignCdsSequenceAsProtein(dnaSeq, protein, mappings,
905               dna.getGapCharacter()))
906       {
907         alignedCount++;
908       }
909       width = Math.max(dnaSeq.getLength(), width);
910     }
911     int oldwidth;
912     int diff;
913     for (SequenceI dnaSeq : dna.getSequences())
914     {
915       oldwidth = dnaSeq.getLength();
916       diff = width - oldwidth;
917       if (diff > 0)
918       {
919         dnaSeq.insertCharAt(oldwidth, diff, dna.getGapCharacter());
920       }
921     }
922     return alignedCount;
923   }
924
925   /**
926    * Helper method to align (if possible) the dna sequence to match the
927    * alignment of a mapped protein sequence. This is currently limited to
928    * handling coding sequence only.
929    * 
930    * @param cdsSeq
931    * @param protein
932    * @param mappings
933    * @param gapChar
934    * @return
935    */
936   static boolean alignCdsSequenceAsProtein(SequenceI cdsSeq,
937           AlignmentI protein, List<AlignedCodonFrame> mappings,
938           char gapChar)
939   {
940     SequenceI cdsDss = cdsSeq.getDatasetSequence();
941     if (cdsDss == null)
942     {
943       System.err
944               .println("alignCdsSequenceAsProtein needs aligned sequence!");
945       return false;
946     }
947
948     List<AlignedCodonFrame> dnaMappings = MappingUtils
949             .findMappingsForSequence(cdsSeq, mappings);
950     for (AlignedCodonFrame mapping : dnaMappings)
951     {
952       SequenceI peptide = mapping.findAlignedSequence(cdsSeq, protein);
953       if (peptide != null)
954       {
955         final int peptideLength = peptide.getLength();
956         Mapping map = mapping.getMappingBetween(cdsSeq, peptide);
957         if (map != null)
958         {
959           MapList mapList = map.getMap();
960           if (map.getTo() == peptide.getDatasetSequence())
961           {
962             mapList = mapList.getInverse();
963           }
964           final int cdsLength = cdsDss.getLength();
965           int mappedFromLength = MappingUtils.getLength(mapList
966                   .getFromRanges());
967           int mappedToLength = MappingUtils
968                   .getLength(mapList.getToRanges());
969           boolean addStopCodon = (cdsLength == mappedFromLength
970                   * CODON_LENGTH + CODON_LENGTH)
971                   || (peptide.getDatasetSequence()
972                           .getLength() == mappedFromLength - 1);
973           if (cdsLength != mappedToLength && !addStopCodon)
974           {
975             System.err.println(String.format(
976                     "Can't align cds as protein (length mismatch %d/%d): %s",
977                     cdsLength, mappedToLength, cdsSeq.getName()));
978           }
979
980           /*
981            * pre-fill the aligned cds sequence with gaps
982            */
983           char[] alignedCds = new char[peptideLength * CODON_LENGTH
984                   + (addStopCodon ? CODON_LENGTH : 0)];
985           Arrays.fill(alignedCds, gapChar);
986
987           /*
988            * walk over the aligned peptide sequence and insert mapped 
989            * codons for residues in the aligned cds sequence 
990            */
991           int copiedBases = 0;
992           int cdsStart = cdsDss.getStart();
993           int proteinPos = peptide.getStart() - 1;
994           int cdsCol = 0;
995
996           for (int col = 0; col < peptideLength; col++)
997           {
998             char residue = peptide.getCharAt(col);
999
1000             if (Comparison.isGap(residue))
1001             {
1002               cdsCol += CODON_LENGTH;
1003             }
1004             else
1005             {
1006               proteinPos++;
1007               int[] codon = mapList.locateInTo(proteinPos, proteinPos);
1008               if (codon == null)
1009               {
1010                 // e.g. incomplete start codon, X in peptide
1011                 cdsCol += CODON_LENGTH;
1012               }
1013               else
1014               {
1015                 for (int j = codon[0]; j <= codon[1]; j++)
1016                 {
1017                   char mappedBase = cdsDss.getCharAt(j - cdsStart);
1018                   alignedCds[cdsCol++] = mappedBase;
1019                   copiedBases++;
1020                 }
1021               }
1022             }
1023           }
1024
1025           /*
1026            * append stop codon if not mapped from protein,
1027            * closing it up to the end of the mapped sequence
1028            */
1029           if (copiedBases == cdsLength - CODON_LENGTH)
1030           {
1031             for (int i = alignedCds.length - 1; i >= 0; i--)
1032             {
1033               if (!Comparison.isGap(alignedCds[i]))
1034               {
1035                 cdsCol = i + 1; // gap just after end of sequence
1036                 break;
1037               }
1038             }
1039             for (int i = cdsLength - CODON_LENGTH; i < cdsLength; i++)
1040             {
1041               alignedCds[cdsCol++] = cdsDss.getCharAt(i);
1042             }
1043           }
1044           cdsSeq.setSequence(new String(alignedCds));
1045           return true;
1046         }
1047       }
1048     }
1049     return false;
1050   }
1051
1052   /**
1053    * Builds a map whose key is an aligned codon position (3 alignment column
1054    * numbers base 0), and whose value is a map from protein sequence to each
1055    * protein's peptide residue for that codon. The map generates an ordering of
1056    * the codons, and allows us to read off the peptides at each position in
1057    * order to assemble 'aligned' protein sequences.
1058    * 
1059    * @param protein
1060    *          the protein alignment
1061    * @param dna
1062    *          the coding dna alignment
1063    * @param unmappedProtein
1064    *          any unmapped proteins are added to this list
1065    * @return
1066    */
1067   protected static Map<AlignedCodon, Map<SequenceI, AlignedCodon>> buildCodonColumnsMap(
1068           AlignmentI protein, AlignmentI dna,
1069           List<SequenceI> unmappedProtein)
1070   {
1071     /*
1072      * maintain a list of any proteins with no mappings - these will be
1073      * rendered 'as is' in the protein alignment as we can't align them
1074      */
1075     unmappedProtein.addAll(protein.getSequences());
1076
1077     List<AlignedCodonFrame> mappings = protein.getCodonFrames();
1078
1079     /*
1080      * Map will hold, for each aligned codon position e.g. [3, 5, 6], a map of
1081      * {dnaSequence, {proteinSequence, codonProduct}} at that position. The
1082      * comparator keeps the codon positions ordered.
1083      */
1084     Map<AlignedCodon, Map<SequenceI, AlignedCodon>> alignedCodons = new TreeMap<AlignedCodon, Map<SequenceI, AlignedCodon>>(
1085             new CodonComparator());
1086
1087     for (SequenceI dnaSeq : dna.getSequences())
1088     {
1089       for (AlignedCodonFrame mapping : mappings)
1090       {
1091         SequenceI prot = mapping.findAlignedSequence(dnaSeq, protein);
1092         if (prot != null)
1093         {
1094           Mapping seqMap = mapping.getMappingForSequence(dnaSeq);
1095           addCodonPositions(dnaSeq, prot, protein.getGapCharacter(), seqMap,
1096                   alignedCodons);
1097           unmappedProtein.remove(prot);
1098         }
1099       }
1100     }
1101
1102     /*
1103      * Finally add any unmapped peptide start residues (e.g. for incomplete
1104      * codons) as if at the codon position before the second residue
1105      */
1106     // TODO resolve JAL-2022 so this fudge can be removed
1107     int mappedSequenceCount = protein.getHeight() - unmappedProtein.size();
1108     addUnmappedPeptideStarts(alignedCodons, mappedSequenceCount);
1109
1110     return alignedCodons;
1111   }
1112
1113   /**
1114    * Scans for any protein mapped from position 2 (meaning unmapped start
1115    * position e.g. an incomplete codon), and synthesizes a 'codon' for it at the
1116    * preceding position in the alignment
1117    * 
1118    * @param alignedCodons
1119    *          the codon-to-peptide map
1120    * @param mappedSequenceCount
1121    *          the number of distinct sequences in the map
1122    */
1123   protected static void addUnmappedPeptideStarts(
1124           Map<AlignedCodon, Map<SequenceI, AlignedCodon>> alignedCodons,
1125           int mappedSequenceCount)
1126   {
1127     // TODO delete this ugly hack once JAL-2022 is resolved
1128     // i.e. we can model startPhase > 0 (incomplete start codon)
1129
1130     List<SequenceI> sequencesChecked = new ArrayList<SequenceI>();
1131     AlignedCodon lastCodon = null;
1132     Map<SequenceI, AlignedCodon> toAdd = new HashMap<SequenceI, AlignedCodon>();
1133
1134     for (Entry<AlignedCodon, Map<SequenceI, AlignedCodon>> entry : alignedCodons
1135             .entrySet())
1136     {
1137       for (Entry<SequenceI, AlignedCodon> sequenceCodon : entry.getValue()
1138               .entrySet())
1139       {
1140         SequenceI seq = sequenceCodon.getKey();
1141         if (sequencesChecked.contains(seq))
1142         {
1143           continue;
1144         }
1145         sequencesChecked.add(seq);
1146         AlignedCodon codon = sequenceCodon.getValue();
1147         if (codon.peptideCol > 1)
1148         {
1149           System.err.println(
1150                   "Problem mapping protein with >1 unmapped start positions: "
1151                           + seq.getName());
1152         }
1153         else if (codon.peptideCol == 1)
1154         {
1155           /*
1156            * first position (peptideCol == 0) was unmapped - add it
1157            */
1158           if (lastCodon != null)
1159           {
1160             AlignedCodon firstPeptide = new AlignedCodon(lastCodon.pos1,
1161                     lastCodon.pos2, lastCodon.pos3,
1162                     String.valueOf(seq.getCharAt(0)), 0);
1163             toAdd.put(seq, firstPeptide);
1164           }
1165           else
1166           {
1167             /*
1168              * unmapped residue at start of alignment (no prior column) -
1169              * 'insert' at nominal codon [0, 0, 0]
1170              */
1171             AlignedCodon firstPeptide = new AlignedCodon(0, 0, 0,
1172                     String.valueOf(seq.getCharAt(0)), 0);
1173             toAdd.put(seq, firstPeptide);
1174           }
1175         }
1176         if (sequencesChecked.size() == mappedSequenceCount)
1177         {
1178           // no need to check past first mapped position in all sequences
1179           break;
1180         }
1181       }
1182       lastCodon = entry.getKey();
1183     }
1184
1185     /*
1186      * add any new codons safely after iterating over the map
1187      */
1188     for (Entry<SequenceI, AlignedCodon> startCodon : toAdd.entrySet())
1189     {
1190       addCodonToMap(alignedCodons, startCodon.getValue(),
1191               startCodon.getKey());
1192     }
1193   }
1194
1195   /**
1196    * Update the aligned protein sequences to match the codon alignments given in
1197    * the map.
1198    * 
1199    * @param protein
1200    * @param alignedCodons
1201    *          an ordered map of codon positions (columns), with sequence/peptide
1202    *          values present in each column
1203    * @param unmappedProtein
1204    * @return
1205    */
1206   protected static int alignProteinAs(AlignmentI protein,
1207           Map<AlignedCodon, Map<SequenceI, AlignedCodon>> alignedCodons,
1208           List<SequenceI> unmappedProtein)
1209   {
1210     /*
1211      * prefill peptide sequences with gaps 
1212      */
1213     int alignedWidth = alignedCodons.size();
1214     char[] gaps = new char[alignedWidth];
1215     Arrays.fill(gaps, protein.getGapCharacter());
1216     Map<SequenceI, char[]> peptides = new HashMap<>();
1217     for (SequenceI seq : protein.getSequences())
1218     {
1219       if (!unmappedProtein.contains(seq))
1220       {
1221         peptides.put(seq, Arrays.copyOf(gaps, gaps.length));
1222       }
1223     }
1224
1225     /*
1226      * Traverse the codons left to right (as defined by CodonComparator)
1227      * and insert peptides in each column where the sequence is mapped.
1228      * This gives a peptide 'alignment' where residues are aligned if their
1229      * corresponding codons occupy the same columns in the cdna alignment.
1230      */
1231     int column = 0;
1232     for (AlignedCodon codon : alignedCodons.keySet())
1233     {
1234       final Map<SequenceI, AlignedCodon> columnResidues = alignedCodons
1235               .get(codon);
1236       for (Entry<SequenceI, AlignedCodon> entry : columnResidues.entrySet())
1237       {
1238         char residue = entry.getValue().product.charAt(0);
1239         peptides.get(entry.getKey())[column] = residue;
1240       }
1241       column++;
1242     }
1243
1244     /*
1245      * and finally set the constructed sequences
1246      */
1247     for (Entry<SequenceI, char[]> entry : peptides.entrySet())
1248     {
1249       entry.getKey().setSequence(new String(entry.getValue()));
1250     }
1251
1252     return 0;
1253   }
1254
1255   /**
1256    * Populate the map of aligned codons by traversing the given sequence
1257    * mapping, locating the aligned positions of mapped codons, and adding those
1258    * positions and their translation products to the map.
1259    * 
1260    * @param dna
1261    *          the aligned sequence we are mapping from
1262    * @param protein
1263    *          the sequence to be aligned to the codons
1264    * @param gapChar
1265    *          the gap character in the dna sequence
1266    * @param seqMap
1267    *          a mapping to a sequence translation
1268    * @param alignedCodons
1269    *          the map we are building up
1270    */
1271   static void addCodonPositions(SequenceI dna, SequenceI protein,
1272           char gapChar, Mapping seqMap,
1273           Map<AlignedCodon, Map<SequenceI, AlignedCodon>> alignedCodons)
1274   {
1275     Iterator<AlignedCodon> codons = seqMap.getCodonIterator(dna, gapChar);
1276
1277     /*
1278      * add codon positions, and their peptide translations, to the alignment
1279      * map, while remembering the first codon mapped
1280      */
1281     while (codons.hasNext())
1282     {
1283       try
1284       {
1285         AlignedCodon codon = codons.next();
1286         addCodonToMap(alignedCodons, codon, protein);
1287       } catch (IncompleteCodonException e)
1288       {
1289         // possible incomplete trailing codon - ignore
1290       } catch (NoSuchElementException e)
1291       {
1292         // possibly peptide lacking STOP
1293       }
1294     }
1295   }
1296
1297   /**
1298    * Helper method to add a codon-to-peptide entry to the aligned codons map
1299    * 
1300    * @param alignedCodons
1301    * @param codon
1302    * @param protein
1303    */
1304   protected static void addCodonToMap(
1305           Map<AlignedCodon, Map<SequenceI, AlignedCodon>> alignedCodons,
1306           AlignedCodon codon, SequenceI protein)
1307   {
1308     Map<SequenceI, AlignedCodon> seqProduct = alignedCodons.get(codon);
1309     if (seqProduct == null)
1310     {
1311       seqProduct = new HashMap<SequenceI, AlignedCodon>();
1312       alignedCodons.put(codon, seqProduct);
1313     }
1314     seqProduct.put(protein, codon);
1315   }
1316
1317   /**
1318    * Returns true if a cDNA/Protein mapping either exists, or could be made,
1319    * between at least one pair of sequences in the two alignments. Currently,
1320    * the logic is:
1321    * <ul>
1322    * <li>One alignment must be nucleotide, and the other protein</li>
1323    * <li>At least one pair of sequences must be already mapped, or mappable</li>
1324    * <li>Mappable means the nucleotide translation matches the protein
1325    * sequence</li>
1326    * <li>The translation may ignore start and stop codons if present in the
1327    * nucleotide</li>
1328    * </ul>
1329    * 
1330    * @param al1
1331    * @param al2
1332    * @return
1333    */
1334   public static boolean isMappable(AlignmentI al1, AlignmentI al2)
1335   {
1336     if (al1 == null || al2 == null)
1337     {
1338       return false;
1339     }
1340
1341     /*
1342      * Require one nucleotide and one protein
1343      */
1344     if (al1.isNucleotide() == al2.isNucleotide())
1345     {
1346       return false;
1347     }
1348     AlignmentI dna = al1.isNucleotide() ? al1 : al2;
1349     AlignmentI protein = dna == al1 ? al2 : al1;
1350     List<AlignedCodonFrame> mappings = protein.getCodonFrames();
1351     for (SequenceI dnaSeq : dna.getSequences())
1352     {
1353       for (SequenceI proteinSeq : protein.getSequences())
1354       {
1355         if (isMappable(dnaSeq, proteinSeq, mappings))
1356         {
1357           return true;
1358         }
1359       }
1360     }
1361     return false;
1362   }
1363
1364   /**
1365    * Returns true if the dna sequence is mapped, or could be mapped, to the
1366    * protein sequence.
1367    * 
1368    * @param dnaSeq
1369    * @param proteinSeq
1370    * @param mappings
1371    * @return
1372    */
1373   protected static boolean isMappable(SequenceI dnaSeq,
1374           SequenceI proteinSeq, List<AlignedCodonFrame> mappings)
1375   {
1376     if (dnaSeq == null || proteinSeq == null)
1377     {
1378       return false;
1379     }
1380
1381     SequenceI dnaDs = dnaSeq.getDatasetSequence() == null ? dnaSeq
1382             : dnaSeq.getDatasetSequence();
1383     SequenceI proteinDs = proteinSeq.getDatasetSequence() == null
1384             ? proteinSeq
1385             : proteinSeq.getDatasetSequence();
1386
1387     for (AlignedCodonFrame mapping : mappings)
1388     {
1389       if (proteinDs == mapping.getAaForDnaSeq(dnaDs))
1390       {
1391         /*
1392          * already mapped
1393          */
1394         return true;
1395       }
1396     }
1397
1398     /*
1399      * Just try to make a mapping (it is not yet stored), test whether
1400      * successful.
1401      */
1402     return mapCdnaToProtein(proteinDs, dnaDs) != null;
1403   }
1404
1405   /**
1406    * Finds any reference annotations associated with the sequences in
1407    * sequenceScope, that are not already added to the alignment, and adds them
1408    * to the 'candidates' map. Also populates a lookup table of annotation
1409    * labels, keyed by calcId, for use in constructing tooltips or the like.
1410    * 
1411    * @param sequenceScope
1412    *          the sequences to scan for reference annotations
1413    * @param labelForCalcId
1414    *          (optional) map to populate with label for calcId
1415    * @param candidates
1416    *          map to populate with annotations for sequence
1417    * @param al
1418    *          the alignment to check for presence of annotations
1419    */
1420   public static void findAddableReferenceAnnotations(
1421           List<SequenceI> sequenceScope, Map<String, String> labelForCalcId,
1422           final Map<SequenceI, List<AlignmentAnnotation>> candidates,
1423           AlignmentI al)
1424   {
1425     if (sequenceScope == null)
1426     {
1427       return;
1428     }
1429
1430     /*
1431      * For each sequence in scope, make a list of any annotations on the
1432      * underlying dataset sequence which are not already on the alignment.
1433      * 
1434      * Add to a map of { alignmentSequence, <List of annotations to add> }
1435      */
1436     for (SequenceI seq : sequenceScope)
1437     {
1438       SequenceI dataset = seq.getDatasetSequence();
1439       if (dataset == null)
1440       {
1441         continue;
1442       }
1443       AlignmentAnnotation[] datasetAnnotations = dataset.getAnnotation();
1444       if (datasetAnnotations == null)
1445       {
1446         continue;
1447       }
1448       final List<AlignmentAnnotation> result = new ArrayList<AlignmentAnnotation>();
1449       for (AlignmentAnnotation dsann : datasetAnnotations)
1450       {
1451         /*
1452          * Find matching annotations on the alignment. If none is found, then
1453          * add this annotation to the list of 'addable' annotations for this
1454          * sequence.
1455          */
1456         final Iterable<AlignmentAnnotation> matchedAlignmentAnnotations = al
1457                 .findAnnotations(seq, dsann.getCalcId(), dsann.label);
1458         if (!matchedAlignmentAnnotations.iterator().hasNext())
1459         {
1460           result.add(dsann);
1461           if (labelForCalcId != null)
1462           {
1463             labelForCalcId.put(dsann.getCalcId(), dsann.label);
1464           }
1465         }
1466       }
1467       /*
1468        * Save any addable annotations for this sequence
1469        */
1470       if (!result.isEmpty())
1471       {
1472         candidates.put(seq, result);
1473       }
1474     }
1475   }
1476
1477   /**
1478    * Adds annotations to the top of the alignment annotations, in the same order
1479    * as their related sequences.
1480    * 
1481    * @param annotations
1482    *          the annotations to add
1483    * @param alignment
1484    *          the alignment to add them to
1485    * @param selectionGroup
1486    *          current selection group (or null if none)
1487    */
1488   public static void addReferenceAnnotations(
1489           Map<SequenceI, List<AlignmentAnnotation>> annotations,
1490           final AlignmentI alignment, final SequenceGroup selectionGroup)
1491   {
1492     for (SequenceI seq : annotations.keySet())
1493     {
1494       for (AlignmentAnnotation ann : annotations.get(seq))
1495       {
1496         AlignmentAnnotation copyAnn = new AlignmentAnnotation(ann);
1497         int startRes = 0;
1498         int endRes = ann.annotations.length;
1499         if (selectionGroup != null)
1500         {
1501           startRes = selectionGroup.getStartRes();
1502           endRes = selectionGroup.getEndRes();
1503         }
1504         copyAnn.restrict(startRes, endRes);
1505
1506         /*
1507          * Add to the sequence (sets copyAnn.datasetSequence), unless the
1508          * original annotation is already on the sequence.
1509          */
1510         if (!seq.hasAnnotation(ann))
1511         {
1512           seq.addAlignmentAnnotation(copyAnn);
1513         }
1514         // adjust for gaps
1515         copyAnn.adjustForAlignment();
1516         // add to the alignment and set visible
1517         alignment.addAnnotation(copyAnn);
1518         copyAnn.visible = true;
1519       }
1520     }
1521   }
1522
1523   /**
1524    * Set visibility of alignment annotations of specified types (labels), for
1525    * specified sequences. This supports controls like "Show all secondary
1526    * structure", "Hide all Temp factor", etc.
1527    * 
1528    * @al the alignment to scan for annotations
1529    * @param types
1530    *          the types (labels) of annotations to be updated
1531    * @param forSequences
1532    *          if not null, only annotations linked to one of these sequences are
1533    *          in scope for update; if null, acts on all sequence annotations
1534    * @param anyType
1535    *          if this flag is true, 'types' is ignored (label not checked)
1536    * @param doShow
1537    *          if true, set visibility on, else set off
1538    */
1539   public static void showOrHideSequenceAnnotations(AlignmentI al,
1540           Collection<String> types, List<SequenceI> forSequences,
1541           boolean anyType, boolean doShow)
1542   {
1543     AlignmentAnnotation[] anns = al.getAlignmentAnnotation();
1544     if (anns != null)
1545     {
1546       for (AlignmentAnnotation aa : anns)
1547       {
1548         if (anyType || types.contains(aa.label))
1549         {
1550           if ((aa.sequenceRef != null) && (forSequences == null
1551                   || forSequences.contains(aa.sequenceRef)))
1552           {
1553             aa.visible = doShow;
1554           }
1555         }
1556       }
1557     }
1558   }
1559
1560   /**
1561    * Returns true if either sequence has a cross-reference to the other
1562    * 
1563    * @param seq1
1564    * @param seq2
1565    * @return
1566    */
1567   public static boolean haveCrossRef(SequenceI seq1, SequenceI seq2)
1568   {
1569     // Note: moved here from class CrossRef as the latter class has dependencies
1570     // not availability to the applet's classpath
1571     return hasCrossRef(seq1, seq2) || hasCrossRef(seq2, seq1);
1572   }
1573
1574   /**
1575    * Returns true if seq1 has a cross-reference to seq2. Currently this assumes
1576    * that sequence name is structured as Source|AccessionId.
1577    * 
1578    * @param seq1
1579    * @param seq2
1580    * @return
1581    */
1582   public static boolean hasCrossRef(SequenceI seq1, SequenceI seq2)
1583   {
1584     if (seq1 == null || seq2 == null)
1585     {
1586       return false;
1587     }
1588     String name = seq2.getName();
1589     final DBRefEntry[] xrefs = seq1.getDBRefs();
1590     if (xrefs != null)
1591     {
1592       for (DBRefEntry xref : xrefs)
1593       {
1594         String xrefName = xref.getSource() + "|" + xref.getAccessionId();
1595         // case-insensitive test, consistent with DBRefEntry.equalRef()
1596         if (xrefName.equalsIgnoreCase(name))
1597         {
1598           return true;
1599         }
1600       }
1601     }
1602     return false;
1603   }
1604
1605   /**
1606    * Constructs an alignment consisting of the mapped (CDS) regions in the given
1607    * nucleotide sequences, and updates mappings to match. The CDS sequences are
1608    * added to the original alignment's dataset, which is shared by the new
1609    * alignment. Mappings from nucleotide to CDS, and from CDS to protein, are
1610    * added to the alignment dataset.
1611    * 
1612    * @param dna
1613    *          aligned nucleotide (dna or cds) sequences
1614    * @param dataset
1615    *          the alignment dataset the sequences belong to
1616    * @param products
1617    *          (optional) to restrict results to CDS that map to specified
1618    *          protein products
1619    * @return an alignment whose sequences are the cds-only parts of the dna
1620    *         sequences (or null if no mappings are found)
1621    */
1622   public static AlignmentI makeCdsAlignment(SequenceI[] dna,
1623           AlignmentI dataset, SequenceI[] products)
1624   {
1625     if (dataset == null || dataset.getDataset() != null)
1626     {
1627       throw new IllegalArgumentException(
1628               "IMPLEMENTATION ERROR: dataset.getDataset() must be null!");
1629     }
1630     List<SequenceI> foundSeqs = new ArrayList<SequenceI>();
1631     List<SequenceI> cdsSeqs = new ArrayList<SequenceI>();
1632     List<AlignedCodonFrame> mappings = dataset.getCodonFrames();
1633     HashSet<SequenceI> productSeqs = null;
1634     if (products != null)
1635     {
1636       productSeqs = new HashSet<SequenceI>();
1637       for (SequenceI seq : products)
1638       {
1639         productSeqs.add(seq.getDatasetSequence() == null ? seq
1640                 : seq.getDatasetSequence());
1641       }
1642     }
1643
1644     /*
1645      * Construct CDS sequences from mappings on the alignment dataset.
1646      * The logic is:
1647      * - find the protein product(s) mapped to from each dna sequence
1648      * - if the mapping covers the whole dna sequence (give or take start/stop
1649      *   codon), take the dna as the CDS sequence
1650      * - else search dataset mappings for a suitable dna sequence, i.e. one
1651      *   whose whole sequence is mapped to the protein 
1652      * - if no sequence found, construct one from the dna sequence and mapping
1653      *   (and add it to dataset so it is found if this is repeated)
1654      */
1655     for (SequenceI dnaSeq : dna)
1656     {
1657       SequenceI dnaDss = dnaSeq.getDatasetSequence() == null ? dnaSeq
1658               : dnaSeq.getDatasetSequence();
1659
1660       List<AlignedCodonFrame> seqMappings = MappingUtils
1661               .findMappingsForSequence(dnaSeq, mappings);
1662       for (AlignedCodonFrame mapping : seqMappings)
1663       {
1664         List<Mapping> mappingsFromSequence = mapping
1665                 .getMappingsFromSequence(dnaSeq);
1666
1667         for (Mapping aMapping : mappingsFromSequence)
1668         {
1669           MapList mapList = aMapping.getMap();
1670           if (mapList.getFromRatio() == 1)
1671           {
1672             /*
1673              * not a dna-to-protein mapping (likely dna-to-cds)
1674              */
1675             continue;
1676           }
1677
1678           /*
1679            * skip if mapping is not to one of the target set of proteins
1680            */
1681           SequenceI proteinProduct = aMapping.getTo();
1682           if (productSeqs != null && !productSeqs.contains(proteinProduct))
1683           {
1684             continue;
1685           }
1686
1687           /*
1688            * try to locate the CDS from the dataset mappings;
1689            * guard against duplicate results (for the case that protein has
1690            * dbrefs to both dna and cds sequences)
1691            */
1692           SequenceI cdsSeq = findCdsForProtein(mappings, dnaSeq,
1693                   seqMappings, aMapping);
1694           if (cdsSeq != null)
1695           {
1696             if (!foundSeqs.contains(cdsSeq))
1697             {
1698               foundSeqs.add(cdsSeq);
1699               SequenceI derivedSequence = cdsSeq.deriveSequence();
1700               cdsSeqs.add(derivedSequence);
1701               if (!dataset.getSequences().contains(cdsSeq))
1702               {
1703                 dataset.addSequence(cdsSeq);
1704               }
1705             }
1706             continue;
1707           }
1708
1709           /*
1710            * didn't find mapped CDS sequence - construct it and add
1711            * its dataset sequence to the dataset
1712            */
1713           cdsSeq = makeCdsSequence(dnaSeq.getDatasetSequence(), aMapping,
1714                   dataset).deriveSequence();
1715           // cdsSeq has a name constructed as CDS|<dbref>
1716           // <dbref> will be either the accession for the coding sequence,
1717           // marked in the /via/ dbref to the protein product accession
1718           // or it will be the original nucleotide accession.
1719           SequenceI cdsSeqDss = cdsSeq.getDatasetSequence();
1720
1721           cdsSeqs.add(cdsSeq);
1722
1723           if (!dataset.getSequences().contains(cdsSeqDss))
1724           {
1725             // check if this sequence is a newly created one
1726             // so needs adding to the dataset
1727             dataset.addSequence(cdsSeqDss);
1728           }
1729
1730           /*
1731            * add a mapping from CDS to the (unchanged) mapped to range
1732            */
1733           List<int[]> cdsRange = Collections
1734                   .singletonList(new int[]
1735                   { 1, cdsSeq.getLength() });
1736           MapList cdsToProteinMap = new MapList(cdsRange,
1737                   mapList.getToRanges(), mapList.getFromRatio(),
1738                   mapList.getToRatio());
1739           AlignedCodonFrame cdsToProteinMapping = new AlignedCodonFrame();
1740           cdsToProteinMapping.addMap(cdsSeqDss, proteinProduct,
1741                   cdsToProteinMap);
1742
1743           /*
1744            * guard against duplicating the mapping if repeating this action
1745            */
1746           if (!mappings.contains(cdsToProteinMapping))
1747           {
1748             mappings.add(cdsToProteinMapping);
1749           }
1750
1751           propagateDBRefsToCDS(cdsSeqDss, dnaSeq.getDatasetSequence(),
1752                   proteinProduct, aMapping);
1753           /*
1754            * add another mapping from original 'from' range to CDS
1755            */
1756           AlignedCodonFrame dnaToCdsMapping = new AlignedCodonFrame();
1757           MapList dnaToCdsMap = new MapList(mapList.getFromRanges(),
1758                   cdsRange, 1, 1);
1759           dnaToCdsMapping.addMap(dnaSeq.getDatasetSequence(), cdsSeqDss,
1760                   dnaToCdsMap);
1761           if (!mappings.contains(dnaToCdsMapping))
1762           {
1763             mappings.add(dnaToCdsMapping);
1764           }
1765
1766           /*
1767            * add DBRef with mapping from protein to CDS
1768            * (this enables Get Cross-References from protein alignment)
1769            * This is tricky because we can't have two DBRefs with the
1770            * same source and accession, so need a different accession for
1771            * the CDS from the dna sequence
1772            */
1773
1774           // specific use case:
1775           // Genomic contig ENSCHR:1, contains coding regions for ENSG01,
1776           // ENSG02, ENSG03, with transcripts and products similarly named.
1777           // cannot add distinct dbrefs mapping location on ENSCHR:1 to ENSG01
1778
1779           // JBPNote: ?? can't actually create an example that demonstrates we
1780           // need to
1781           // synthesize an xref.
1782
1783           for (DBRefEntry primRef : dnaDss.getPrimaryDBRefs())
1784           {
1785             // creates a complementary cross-reference to the source sequence's
1786             // primary reference.
1787
1788             DBRefEntry cdsCrossRef = new DBRefEntry(primRef.getSource(),
1789                     primRef.getSource() + ":" + primRef.getVersion(),
1790                     primRef.getAccessionId());
1791             cdsCrossRef
1792                     .setMap(new Mapping(dnaDss, new MapList(dnaToCdsMap)));
1793             cdsSeqDss.addDBRef(cdsCrossRef);
1794
1795             // problem here is that the cross-reference is synthesized -
1796             // cdsSeq.getName() may be like 'CDS|dnaaccession' or
1797             // 'CDS|emblcdsacc'
1798             // assuming cds version same as dna ?!?
1799
1800             DBRefEntry proteinToCdsRef = new DBRefEntry(primRef.getSource(),
1801                     primRef.getVersion(), cdsSeq.getName());
1802             //
1803             proteinToCdsRef.setMap(
1804                     new Mapping(cdsSeqDss, cdsToProteinMap.getInverse()));
1805             proteinProduct.addDBRef(proteinToCdsRef);
1806           }
1807
1808           /*
1809            * transfer any features on dna that overlap the CDS
1810            */
1811           transferFeatures(dnaSeq, cdsSeq, dnaToCdsMap, null,
1812                   SequenceOntologyI.CDS);
1813         }
1814       }
1815     }
1816
1817     AlignmentI cds = new Alignment(
1818             cdsSeqs.toArray(new SequenceI[cdsSeqs.size()]));
1819     cds.setDataset(dataset);
1820
1821     return cds;
1822   }
1823
1824   /**
1825    * A helper method that finds a CDS sequence in the alignment dataset that is
1826    * mapped to the given protein sequence, and either is, or has a mapping from,
1827    * the given dna sequence.
1828    * 
1829    * @param mappings
1830    *          set of all mappings on the dataset
1831    * @param dnaSeq
1832    *          a dna (or cds) sequence we are searching from
1833    * @param seqMappings
1834    *          the set of mappings involving dnaSeq
1835    * @param aMapping
1836    *          an initial candidate from seqMappings
1837    * @return
1838    */
1839   static SequenceI findCdsForProtein(List<AlignedCodonFrame> mappings,
1840           SequenceI dnaSeq, List<AlignedCodonFrame> seqMappings,
1841           Mapping aMapping)
1842   {
1843     /*
1844      * TODO a better dna-cds-protein mapping data representation to allow easy
1845      * navigation; until then this clunky looping around lists of mappings
1846      */
1847     SequenceI seqDss = dnaSeq.getDatasetSequence() == null ? dnaSeq
1848             : dnaSeq.getDatasetSequence();
1849     SequenceI proteinProduct = aMapping.getTo();
1850
1851     /*
1852      * is this mapping from the whole dna sequence (i.e. CDS)?
1853      * allowing for possible stop codon on dna but not peptide
1854      */
1855     int mappedFromLength = MappingUtils
1856             .getLength(aMapping.getMap().getFromRanges());
1857     int dnaLength = seqDss.getLength();
1858     if (mappedFromLength == dnaLength
1859             || mappedFromLength == dnaLength - CODON_LENGTH)
1860     {
1861       return seqDss;
1862     }
1863
1864     /*
1865      * looks like we found the dna-to-protein mapping; search for the
1866      * corresponding cds-to-protein mapping
1867      */
1868     List<AlignedCodonFrame> mappingsToPeptide = MappingUtils
1869             .findMappingsForSequence(proteinProduct, mappings);
1870     for (AlignedCodonFrame acf : mappingsToPeptide)
1871     {
1872       for (SequenceToSequenceMapping map : acf.getMappings())
1873       {
1874         Mapping mapping = map.getMapping();
1875         if (mapping != aMapping
1876                 && mapping.getMap().getFromRatio() == CODON_LENGTH
1877                 && proteinProduct == mapping.getTo()
1878                 && seqDss != map.getFromSeq())
1879         {
1880           mappedFromLength = MappingUtils
1881                   .getLength(mapping.getMap().getFromRanges());
1882           if (mappedFromLength == map.getFromSeq().getLength())
1883           {
1884             /*
1885             * found a 3:1 mapping to the protein product which covers
1886             * the whole dna sequence i.e. is from CDS; finally check it
1887             * is from the dna start sequence
1888             */
1889             SequenceI cdsSeq = map.getFromSeq();
1890             List<AlignedCodonFrame> dnaToCdsMaps = MappingUtils
1891                     .findMappingsForSequence(cdsSeq, seqMappings);
1892             if (!dnaToCdsMaps.isEmpty())
1893             {
1894               return cdsSeq;
1895             }
1896           }
1897         }
1898       }
1899     }
1900     return null;
1901   }
1902
1903   /**
1904    * Helper method that makes a CDS sequence as defined by the mappings from the
1905    * given sequence i.e. extracts the 'mapped from' ranges (which may be on
1906    * forward or reverse strand).
1907    * 
1908    * @param seq
1909    * @param mapping
1910    * @param dataset
1911    *          - existing dataset. We check for sequences that look like the CDS
1912    *          we are about to construct, if one exists already, then we will
1913    *          just return that one.
1914    * @return CDS sequence (as a dataset sequence)
1915    */
1916   static SequenceI makeCdsSequence(SequenceI seq, Mapping mapping,
1917           AlignmentI dataset)
1918   {
1919     char[] seqChars = seq.getSequence();
1920     List<int[]> fromRanges = mapping.getMap().getFromRanges();
1921     int cdsWidth = MappingUtils.getLength(fromRanges);
1922     char[] newSeqChars = new char[cdsWidth];
1923
1924     int newPos = 0;
1925     for (int[] range : fromRanges)
1926     {
1927       if (range[0] <= range[1])
1928       {
1929         // forward strand mapping - just copy the range
1930         int length = range[1] - range[0] + 1;
1931         System.arraycopy(seqChars, range[0] - 1, newSeqChars, newPos,
1932                 length);
1933         newPos += length;
1934       }
1935       else
1936       {
1937         // reverse strand mapping - copy and complement one by one
1938         for (int i = range[0]; i >= range[1]; i--)
1939         {
1940           newSeqChars[newPos++] = Dna.getComplement(seqChars[i - 1]);
1941         }
1942       }
1943     }
1944
1945     /*
1946      * assign 'from id' held in the mapping if set (e.g. EMBL protein_id),
1947      * else generate a sequence name
1948      */
1949     String mapFromId = mapping.getMappedFromId();
1950     String seqId = "CDS|" + (mapFromId != null ? mapFromId : seq.getName());
1951     SequenceI newSeq = new Sequence(seqId, newSeqChars, 1, newPos);
1952     if (dataset != null)
1953     {
1954       SequenceI[] matches = dataset.findSequenceMatch(newSeq.getName());
1955       if (matches != null)
1956       {
1957         boolean matched = false;
1958         for (SequenceI mtch : matches)
1959         {
1960           if (mtch.getStart() != newSeq.getStart())
1961           {
1962             continue;
1963           }
1964           if (mtch.getEnd() != newSeq.getEnd())
1965           {
1966             continue;
1967           }
1968           if (!Arrays.equals(mtch.getSequence(), newSeq.getSequence()))
1969           {
1970             continue;
1971           }
1972           if (!matched)
1973           {
1974             matched = true;
1975             newSeq = mtch;
1976           }
1977           else
1978           {
1979             System.err.println(
1980                     "JAL-2154 regression: warning - found (and ignnored a duplicate CDS sequence):"
1981                             + mtch.toString());
1982           }
1983         }
1984       }
1985     }
1986     // newSeq.setDescription(mapFromId);
1987
1988     return newSeq;
1989   }
1990
1991   /**
1992    * add any DBRefEntrys to cdsSeq from contig that have a Mapping congruent to
1993    * the given mapping.
1994    * 
1995    * @param cdsSeq
1996    * @param contig
1997    * @param mapping
1998    * @return list of DBRefEntrys added.
1999    */
2000   public static List<DBRefEntry> propagateDBRefsToCDS(SequenceI cdsSeq,
2001           SequenceI contig, SequenceI proteinProduct, Mapping mapping)
2002   {
2003
2004     // gather direct refs from contig congrent with mapping
2005     List<DBRefEntry> direct = new ArrayList<DBRefEntry>();
2006     HashSet<String> directSources = new HashSet<String>();
2007     if (contig.getDBRefs() != null)
2008     {
2009       for (DBRefEntry dbr : contig.getDBRefs())
2010       {
2011         if (dbr.hasMap() && dbr.getMap().getMap().isTripletMap())
2012         {
2013           MapList map = dbr.getMap().getMap();
2014           // check if map is the CDS mapping
2015           if (mapping.getMap().equals(map))
2016           {
2017             direct.add(dbr);
2018             directSources.add(dbr.getSource());
2019           }
2020         }
2021       }
2022     }
2023     DBRefEntry[] onSource = DBRefUtils.selectRefs(
2024             proteinProduct.getDBRefs(),
2025             directSources.toArray(new String[0]));
2026     List<DBRefEntry> propagated = new ArrayList<DBRefEntry>();
2027
2028     // and generate appropriate mappings
2029     for (DBRefEntry cdsref : direct)
2030     {
2031       // clone maplist and mapping
2032       MapList cdsposmap = new MapList(
2033               Arrays.asList(new int[][]
2034               { new int[] { cdsSeq.getStart(), cdsSeq.getEnd() } }),
2035               cdsref.getMap().getMap().getToRanges(), 3, 1);
2036       Mapping cdsmap = new Mapping(cdsref.getMap().getTo(),
2037               cdsref.getMap().getMap());
2038
2039       // create dbref
2040       DBRefEntry newref = new DBRefEntry(cdsref.getSource(),
2041               cdsref.getVersion(), cdsref.getAccessionId(),
2042               new Mapping(cdsmap.getTo(), cdsposmap));
2043
2044       // and see if we can map to the protein product for this mapping.
2045       // onSource is the filtered set of accessions on protein that we are
2046       // tranferring, so we assume accession is the same.
2047       if (cdsmap.getTo() == null && onSource != null)
2048       {
2049         List<DBRefEntry> sourceRefs = DBRefUtils.searchRefs(onSource,
2050                 cdsref.getAccessionId());
2051         if (sourceRefs != null)
2052         {
2053           for (DBRefEntry srcref : sourceRefs)
2054           {
2055             if (srcref.getSource().equalsIgnoreCase(cdsref.getSource()))
2056             {
2057               // we have found a complementary dbref on the protein product, so
2058               // update mapping's getTo
2059               newref.getMap().setTo(proteinProduct);
2060             }
2061           }
2062         }
2063       }
2064       cdsSeq.addDBRef(newref);
2065       propagated.add(newref);
2066     }
2067     return propagated;
2068   }
2069
2070   /**
2071    * Transfers co-located features on 'fromSeq' to 'toSeq', adjusting the
2072    * feature start/end ranges, optionally omitting specified feature types.
2073    * Returns the number of features copied.
2074    * 
2075    * @param fromSeq
2076    * @param toSeq
2077    * @param mapping
2078    *          the mapping from 'fromSeq' to 'toSeq'
2079    * @param select
2080    *          if not null, only features of this type are copied (including
2081    *          subtypes in the Sequence Ontology)
2082    * @param omitting
2083    */
2084   public static int transferFeatures(SequenceI fromSeq, SequenceI toSeq,
2085           MapList mapping, String select, String... omitting)
2086   {
2087     SequenceI copyTo = toSeq;
2088     while (copyTo.getDatasetSequence() != null)
2089     {
2090       copyTo = copyTo.getDatasetSequence();
2091     }
2092
2093     /*
2094      * get features, optionally restricted by an ontology term
2095      */
2096     List<SequenceFeature> sfs = select == null ? fromSeq.getFeatures()
2097             .getPositionalFeatures() : fromSeq.getFeatures()
2098             .getFeaturesByOntology(select);
2099
2100     int count = 0;
2101     for (SequenceFeature sf : sfs)
2102     {
2103       String type = sf.getType();
2104       boolean omit = false;
2105       for (String toOmit : omitting)
2106       {
2107         if (type.equals(toOmit))
2108         {
2109           omit = true;
2110         }
2111       }
2112       if (omit)
2113       {
2114         continue;
2115       }
2116
2117       /*
2118        * locate the mapped range - null if either start or end is
2119        * not mapped (no partial overlaps are calculated)
2120        */
2121       int start = sf.getBegin();
2122       int end = sf.getEnd();
2123       int[] mappedTo = mapping.locateInTo(start, end);
2124       /*
2125        * if whole exon range doesn't map, try interpreting it
2126        * as 5' or 3' exon overlapping the CDS range
2127        */
2128       if (mappedTo == null)
2129       {
2130         mappedTo = mapping.locateInTo(end, end);
2131         if (mappedTo != null)
2132         {
2133           /*
2134            * end of exon is in CDS range - 5' overlap
2135            * to a range from the start of the peptide
2136            */
2137           mappedTo[0] = 1;
2138         }
2139       }
2140       if (mappedTo == null)
2141       {
2142         mappedTo = mapping.locateInTo(start, start);
2143         if (mappedTo != null)
2144         {
2145           /*
2146            * start of exon is in CDS range - 3' overlap
2147            * to a range up to the end of the peptide
2148            */
2149           mappedTo[1] = toSeq.getLength();
2150         }
2151       }
2152       if (mappedTo != null)
2153       {
2154         int newBegin = Math.min(mappedTo[0], mappedTo[1]);
2155         int newEnd = Math.max(mappedTo[0], mappedTo[1]);
2156         SequenceFeature copy = new SequenceFeature(sf, newBegin, newEnd,
2157                 sf.getFeatureGroup(), sf.getScore());
2158         copyTo.addSequenceFeature(copy);
2159         count++;
2160       }
2161     }
2162     return count;
2163   }
2164
2165   /**
2166    * Returns a mapping from dna to protein by inspecting sequence features of
2167    * type "CDS" on the dna. A mapping is constructed if the total CDS feature
2168    * length is 3 times the peptide length (optionally after dropping a trailing
2169    * stop codon). This method does not check whether the CDS nucleotide sequence
2170    * translates to the peptide sequence.
2171    * 
2172    * @param dnaSeq
2173    * @param proteinSeq
2174    * @return
2175    */
2176   public static MapList mapCdsToProtein(SequenceI dnaSeq,
2177           SequenceI proteinSeq)
2178   {
2179     List<int[]> ranges = findCdsPositions(dnaSeq);
2180     int mappedDnaLength = MappingUtils.getLength(ranges);
2181
2182     /*
2183      * if not a whole number of codons, something is wrong,
2184      * abort mapping
2185      */
2186     if (mappedDnaLength % CODON_LENGTH > 0)
2187     {
2188       return null;
2189     }
2190
2191     int proteinLength = proteinSeq.getLength();
2192     int proteinStart = proteinSeq.getStart();
2193     int proteinEnd = proteinSeq.getEnd();
2194
2195     /*
2196      * incomplete start codon may mean X at start of peptide
2197      * we ignore both for mapping purposes
2198      */
2199     if (proteinSeq.getCharAt(0) == 'X')
2200     {
2201       // todo JAL-2022 support startPhase > 0
2202       proteinStart++;
2203       proteinLength--;
2204     }
2205     List<int[]> proteinRange = new ArrayList<int[]>();
2206
2207     /*
2208      * dna length should map to protein (or protein plus stop codon)
2209      */
2210     int codesForResidues = mappedDnaLength / CODON_LENGTH;
2211     if (codesForResidues == (proteinLength + 1))
2212     {
2213       // assuming extra codon is for STOP and not in peptide
2214       // todo: check trailing codon is indeed a STOP codon
2215       codesForResidues--;
2216       mappedDnaLength -= CODON_LENGTH;
2217       MappingUtils.removeEndPositions(CODON_LENGTH, ranges);
2218     }
2219
2220     if (codesForResidues == proteinLength)
2221     {
2222       proteinRange.add(new int[] { proteinStart, proteinEnd });
2223       return new MapList(ranges, proteinRange, CODON_LENGTH, 1);
2224     }
2225     return null;
2226   }
2227
2228   /**
2229    * Returns a list of CDS ranges found (as sequence positions base 1), i.e. of
2230    * [start, end] positions of sequence features of type "CDS" (or a sub-type of
2231    * CDS in the Sequence Ontology). The ranges are sorted into ascending start
2232    * position order, so this method is only valid for linear CDS in the same
2233    * sense as the protein product.
2234    * 
2235    * @param dnaSeq
2236    * @return
2237    */
2238   public static List<int[]> findCdsPositions(SequenceI dnaSeq)
2239   {
2240     List<int[]> result = new ArrayList<int[]>();
2241
2242     List<SequenceFeature> sfs = dnaSeq.getFeatures().getFeaturesByOntology(
2243             SequenceOntologyI.CDS);
2244     if (sfs.isEmpty())
2245     {
2246       return result;
2247     }
2248     SequenceFeatures.sortFeatures(sfs, true);
2249
2250     for (SequenceFeature sf : sfs)
2251     {
2252       int phase = 0;
2253       try
2254       {
2255         phase = Integer.parseInt(sf.getPhase());
2256       } catch (NumberFormatException e)
2257       {
2258         // ignore
2259       }
2260       /*
2261        * phase > 0 on first codon means 5' incomplete - skip to the start
2262        * of the next codon; example ENST00000496384
2263        */
2264       int begin = sf.getBegin();
2265       int end = sf.getEnd();
2266       if (result.isEmpty() && phase > 0)
2267       {
2268         begin += phase;
2269         if (begin > end)
2270         {
2271           // shouldn't happen!
2272           System.err
2273                   .println("Error: start phase extends beyond start CDS in "
2274                           + dnaSeq.getName());
2275         }
2276       }
2277       result.add(new int[] { begin, end });
2278     }
2279
2280     /*
2281      * Finally sort ranges by start position. This avoids a dependency on 
2282      * keeping features in order on the sequence (if they are in order anyway,
2283      * the sort will have almost no work to do). The implicit assumption is CDS
2284      * ranges are assembled in order. Other cases should not use this method,
2285      * but instead construct an explicit mapping for CDS (e.g. EMBL parsing).
2286      */
2287     Collections.sort(result, IntRangeComparator.ASCENDING);
2288     return result;
2289   }
2290
2291   /**
2292    * Maps exon features from dna to protein, and computes variants in peptide
2293    * product generated by variants in dna, and adds them as sequence_variant
2294    * features on the protein sequence. Returns the number of variant features
2295    * added.
2296    * 
2297    * @param dnaSeq
2298    * @param peptide
2299    * @param dnaToProtein
2300    */
2301   public static int computeProteinFeatures(SequenceI dnaSeq,
2302           SequenceI peptide, MapList dnaToProtein)
2303   {
2304     while (dnaSeq.getDatasetSequence() != null)
2305     {
2306       dnaSeq = dnaSeq.getDatasetSequence();
2307     }
2308     while (peptide.getDatasetSequence() != null)
2309     {
2310       peptide = peptide.getDatasetSequence();
2311     }
2312
2313     transferFeatures(dnaSeq, peptide, dnaToProtein, SequenceOntologyI.EXON);
2314
2315     /*
2316      * compute protein variants from dna variants and codon mappings;
2317      * NB - alternatively we could retrieve this using the REST service e.g.
2318      * http://rest.ensembl.org/overlap/translation
2319      * /ENSP00000288602?feature=transcript_variation;content-type=text/xml
2320      * which would be a bit slower but possibly more reliable
2321      */
2322
2323     /*
2324      * build a map with codon variations for each potentially varying peptide
2325      */
2326     LinkedHashMap<Integer, List<DnaVariant>[]> variants = buildDnaVariantsMap(
2327             dnaSeq, dnaToProtein);
2328
2329     /*
2330      * scan codon variations, compute peptide variants and add to peptide sequence
2331      */
2332     int count = 0;
2333     for (Entry<Integer, List<DnaVariant>[]> variant : variants.entrySet())
2334     {
2335       int peptidePos = variant.getKey();
2336       List<DnaVariant>[] codonVariants = variant.getValue();
2337       count += computePeptideVariants(peptide, peptidePos, codonVariants);
2338     }
2339
2340     return count;
2341   }
2342
2343   /**
2344    * Computes non-synonymous peptide variants from codon variants and adds them
2345    * as sequence_variant features on the protein sequence (one feature per
2346    * allele variant). Selected attributes (variant id, clinical significance)
2347    * are copied over to the new features.
2348    * 
2349    * @param peptide
2350    *          the protein sequence
2351    * @param peptidePos
2352    *          the position to compute peptide variants for
2353    * @param codonVariants
2354    *          a list of dna variants per codon position
2355    * @return the number of features added
2356    */
2357   static int computePeptideVariants(SequenceI peptide, int peptidePos,
2358           List<DnaVariant>[] codonVariants)
2359   {
2360     String residue = String.valueOf(peptide.getCharAt(peptidePos - 1));
2361     int count = 0;
2362     String base1 = codonVariants[0].get(0).base;
2363     String base2 = codonVariants[1].get(0).base;
2364     String base3 = codonVariants[2].get(0).base;
2365
2366     /*
2367      * variants in first codon base
2368      */
2369     for (DnaVariant var : codonVariants[0])
2370     {
2371       if (var.variant != null)
2372       {
2373         String alleles = (String) var.variant.getValue("alleles");
2374         if (alleles != null)
2375         {
2376           for (String base : alleles.split(","))
2377           {
2378             String codon = base + base2 + base3;
2379             if (addPeptideVariant(peptide, peptidePos, residue, var, codon))
2380             {
2381               count++;
2382             }
2383           }
2384         }
2385       }
2386     }
2387
2388     /*
2389      * variants in second codon base
2390      */
2391     for (DnaVariant var : codonVariants[1])
2392     {
2393       if (var.variant != null)
2394       {
2395         String alleles = (String) var.variant.getValue("alleles");
2396         if (alleles != null)
2397         {
2398           for (String base : alleles.split(","))
2399           {
2400             String codon = base1 + base + base3;
2401             if (addPeptideVariant(peptide, peptidePos, residue, var, codon))
2402             {
2403               count++;
2404             }
2405           }
2406         }
2407       }
2408     }
2409
2410     /*
2411      * variants in third codon base
2412      */
2413     for (DnaVariant var : codonVariants[2])
2414     {
2415       if (var.variant != null)
2416       {
2417         String alleles = (String) var.variant.getValue("alleles");
2418         if (alleles != null)
2419         {
2420           for (String base : alleles.split(","))
2421           {
2422             String codon = base1 + base2 + base;
2423             if (addPeptideVariant(peptide, peptidePos, residue, var, codon))
2424             {
2425               count++;
2426             }
2427           }
2428         }
2429       }
2430     }
2431
2432     return count;
2433   }
2434
2435   /**
2436    * Helper method that adds a peptide variant feature, provided the given codon
2437    * translates to a value different to the current residue (is a non-synonymous
2438    * variant). ID and clinical_significance attributes of the dna variant (if
2439    * present) are copied to the new feature.
2440    * 
2441    * @param peptide
2442    * @param peptidePos
2443    * @param residue
2444    * @param var
2445    * @param codon
2446    * @return true if a feature was added, else false
2447    */
2448   static boolean addPeptideVariant(SequenceI peptide, int peptidePos,
2449           String residue, DnaVariant var, String codon)
2450   {
2451     /*
2452      * get peptide translation of codon e.g. GAT -> D
2453      * note that variants which are not single alleles,
2454      * e.g. multibase variants or HGMD_MUTATION etc
2455      * are currently ignored here
2456      */
2457     String trans = codon.contains("-") ? "-"
2458             : (codon.length() > CODON_LENGTH ? null
2459                     : ResidueProperties.codonTranslate(codon));
2460     if (trans != null && !trans.equals(residue))
2461     {
2462       String residue3Char = StringUtils
2463               .toSentenceCase(ResidueProperties.aa2Triplet.get(residue));
2464       String trans3Char = StringUtils
2465               .toSentenceCase(ResidueProperties.aa2Triplet.get(trans));
2466       String desc = "p." + residue3Char + peptidePos + trans3Char;
2467       SequenceFeature sf = new SequenceFeature(
2468               SequenceOntologyI.SEQUENCE_VARIANT, desc, peptidePos,
2469               peptidePos, var.getSource());
2470       StringBuilder attributes = new StringBuilder(32);
2471       String id = (String) var.variant.getValue(ID);
2472       if (id != null)
2473       {
2474         if (id.startsWith(SEQUENCE_VARIANT))
2475         {
2476           id = id.substring(SEQUENCE_VARIANT.length());
2477         }
2478         sf.setValue(ID, id);
2479         attributes.append(ID).append("=").append(id);
2480         // TODO handle other species variants JAL-2064
2481         StringBuilder link = new StringBuilder(32);
2482         try
2483         {
2484           link.append(desc).append(" ").append(id).append(
2485                   "|http://www.ensembl.org/Homo_sapiens/Variation/Summary?v=")
2486                   .append(URLEncoder.encode(id, "UTF-8"));
2487           sf.addLink(link.toString());
2488         } catch (UnsupportedEncodingException e)
2489         {
2490           // as if
2491         }
2492       }
2493       String clinSig = (String) var.variant.getValue(CLINICAL_SIGNIFICANCE);
2494       if (clinSig != null)
2495       {
2496         sf.setValue(CLINICAL_SIGNIFICANCE, clinSig);
2497         attributes.append(";").append(CLINICAL_SIGNIFICANCE).append("=")
2498                 .append(clinSig);
2499       }
2500       peptide.addSequenceFeature(sf);
2501       if (attributes.length() > 0)
2502       {
2503         sf.setAttributes(attributes.toString());
2504       }
2505       return true;
2506     }
2507     return false;
2508   }
2509
2510   /**
2511    * Builds a map whose key is position in the protein sequence, and value is a
2512    * list of the base and all variants for each corresponding codon position
2513    * 
2514    * @param dnaSeq
2515    * @param dnaToProtein
2516    * @return
2517    */
2518   @SuppressWarnings("unchecked")
2519   static LinkedHashMap<Integer, List<DnaVariant>[]> buildDnaVariantsMap(
2520           SequenceI dnaSeq, MapList dnaToProtein)
2521   {
2522     /*
2523      * map from peptide position to all variants of the codon which codes for it
2524      * LinkedHashMap ensures we keep the peptide features in sequence order
2525      */
2526     LinkedHashMap<Integer, List<DnaVariant>[]> variants = new LinkedHashMap<Integer, List<DnaVariant>[]>();
2527
2528     List<SequenceFeature> dnaFeatures = dnaSeq.getFeatures()
2529             .getFeaturesByOntology(SequenceOntologyI.SEQUENCE_VARIANT);
2530     if (dnaFeatures.isEmpty())
2531     {
2532       return variants;
2533     }
2534
2535     int dnaStart = dnaSeq.getStart();
2536     int[] lastCodon = null;
2537     int lastPeptidePostion = 0;
2538
2539     /*
2540      * build a map of codon variations for peptides
2541      */
2542     for (SequenceFeature sf : dnaFeatures)
2543     {
2544       int dnaCol = sf.getBegin();
2545       if (dnaCol != sf.getEnd())
2546       {
2547         // not handling multi-locus variant features
2548         continue;
2549       }
2550       int[] mapsTo = dnaToProtein.locateInTo(dnaCol, dnaCol);
2551       if (mapsTo == null)
2552       {
2553         // feature doesn't lie within coding region
2554         continue;
2555       }
2556       int peptidePosition = mapsTo[0];
2557       List<DnaVariant>[] codonVariants = variants.get(peptidePosition);
2558       if (codonVariants == null)
2559       {
2560         codonVariants = new ArrayList[CODON_LENGTH];
2561         codonVariants[0] = new ArrayList<DnaVariant>();
2562         codonVariants[1] = new ArrayList<DnaVariant>();
2563         codonVariants[2] = new ArrayList<DnaVariant>();
2564         variants.put(peptidePosition, codonVariants);
2565       }
2566
2567       /*
2568        * extract dna variants to a string array
2569        */
2570       String alls = (String) sf.getValue("alleles");
2571       if (alls == null)
2572       {
2573         continue;
2574       }
2575       String[] alleles = alls.toUpperCase().split(",");
2576       int i = 0;
2577       for (String allele : alleles)
2578       {
2579         alleles[i++] = allele.trim(); // lose any space characters "A, G"
2580       }
2581
2582       /*
2583        * get this peptide's codon positions e.g. [3, 4, 5] or [4, 7, 10]
2584        */
2585       int[] codon = peptidePosition == lastPeptidePostion ? lastCodon
2586               : MappingUtils.flattenRanges(dnaToProtein.locateInFrom(
2587                       peptidePosition, peptidePosition));
2588       lastPeptidePostion = peptidePosition;
2589       lastCodon = codon;
2590
2591       /*
2592        * save nucleotide (and any variant) for each codon position
2593        */
2594       for (int codonPos = 0; codonPos < CODON_LENGTH; codonPos++)
2595       {
2596         String nucleotide = String.valueOf(
2597                 dnaSeq.getCharAt(codon[codonPos] - dnaStart)).toUpperCase();
2598         List<DnaVariant> codonVariant = codonVariants[codonPos];
2599         if (codon[codonPos] == dnaCol)
2600         {
2601           if (!codonVariant.isEmpty()
2602                   && codonVariant.get(0).variant == null)
2603           {
2604             /*
2605              * already recorded base value, add this variant
2606              */
2607             codonVariant.get(0).variant = sf;
2608           }
2609           else
2610           {
2611             /*
2612              * add variant with base value
2613              */
2614             codonVariant.add(new DnaVariant(nucleotide, sf));
2615           }
2616         }
2617         else if (codonVariant.isEmpty())
2618         {
2619           /*
2620            * record (possibly non-varying) base value
2621            */
2622           codonVariant.add(new DnaVariant(nucleotide));
2623         }
2624       }
2625     }
2626     return variants;
2627   }
2628
2629   /**
2630    * Makes an alignment with a copy of the given sequences, adding in any
2631    * non-redundant sequences which are mapped to by the cross-referenced
2632    * sequences.
2633    * 
2634    * @param seqs
2635    * @param xrefs
2636    * @param dataset
2637    *          the alignment dataset shared by the new copy
2638    * @return
2639    */
2640   public static AlignmentI makeCopyAlignment(SequenceI[] seqs,
2641           SequenceI[] xrefs, AlignmentI dataset)
2642   {
2643     AlignmentI copy = new Alignment(new Alignment(seqs));
2644     copy.setDataset(dataset);
2645     boolean isProtein = !copy.isNucleotide();
2646     SequenceIdMatcher matcher = new SequenceIdMatcher(seqs);
2647     if (xrefs != null)
2648     {
2649       for (SequenceI xref : xrefs)
2650       {
2651         DBRefEntry[] dbrefs = xref.getDBRefs();
2652         if (dbrefs != null)
2653         {
2654           for (DBRefEntry dbref : dbrefs)
2655           {
2656             if (dbref.getMap() == null || dbref.getMap().getTo() == null
2657                     || dbref.getMap().getTo().isProtein() != isProtein)
2658             {
2659               continue;
2660             }
2661             SequenceI mappedTo = dbref.getMap().getTo();
2662             SequenceI match = matcher.findIdMatch(mappedTo);
2663             if (match == null)
2664             {
2665               matcher.add(mappedTo);
2666               copy.addSequence(mappedTo);
2667             }
2668           }
2669         }
2670       }
2671     }
2672     return copy;
2673   }
2674
2675   /**
2676    * Try to align sequences in 'unaligned' to match the alignment of their
2677    * mapped regions in 'aligned'. For example, could use this to align CDS
2678    * sequences which are mapped to their parent cDNA sequences.
2679    * 
2680    * This method handles 1:1 mappings (dna-to-dna or protein-to-protein). For
2681    * dna-to-protein or protein-to-dna use alternative methods.
2682    * 
2683    * @param unaligned
2684    *          sequences to be aligned
2685    * @param aligned
2686    *          holds aligned sequences and their mappings
2687    * @return
2688    */
2689   public static int alignAs(AlignmentI unaligned, AlignmentI aligned)
2690   {
2691     /*
2692      * easy case - aligning a copy of aligned sequences
2693      */
2694     if (alignAsSameSequences(unaligned, aligned))
2695     {
2696       return unaligned.getHeight();
2697     }
2698
2699     /*
2700      * fancy case - aligning via mappings between sequences
2701      */
2702     List<SequenceI> unmapped = new ArrayList<SequenceI>();
2703     Map<Integer, Map<SequenceI, Character>> columnMap = buildMappedColumnsMap(
2704             unaligned, aligned, unmapped);
2705     int width = columnMap.size();
2706     char gap = unaligned.getGapCharacter();
2707     int realignedCount = 0;
2708     // TODO: verify this loop scales sensibly for very wide/high alignments
2709
2710     for (SequenceI seq : unaligned.getSequences())
2711     {
2712       if (!unmapped.contains(seq))
2713       {
2714         char[] newSeq = new char[width];
2715         Arrays.fill(newSeq, gap); // JBPComment - doubt this is faster than the
2716                                   // Integer iteration below
2717         int newCol = 0;
2718         int lastCol = 0;
2719
2720         /*
2721          * traverse the map to find columns populated
2722          * by our sequence
2723          */
2724         for (Integer column : columnMap.keySet())
2725         {
2726           Character c = columnMap.get(column).get(seq);
2727           if (c != null)
2728           {
2729             /*
2730              * sequence has a character at this position
2731              * 
2732              */
2733             newSeq[newCol] = c;
2734             lastCol = newCol;
2735           }
2736           newCol++;
2737         }
2738
2739         /*
2740          * trim trailing gaps
2741          */
2742         if (lastCol < width)
2743         {
2744           char[] tmp = new char[lastCol + 1];
2745           System.arraycopy(newSeq, 0, tmp, 0, lastCol + 1);
2746           newSeq = tmp;
2747         }
2748         // TODO: optimise SequenceI to avoid char[]->String->char[]
2749         seq.setSequence(String.valueOf(newSeq));
2750         realignedCount++;
2751       }
2752     }
2753     return realignedCount;
2754   }
2755
2756   /**
2757    * If unaligned and aligned sequences share the same dataset sequences, then
2758    * simply copies the aligned sequences to the unaligned sequences and returns
2759    * true; else returns false
2760    * 
2761    * @param unaligned
2762    *          - sequences to be aligned based on aligned
2763    * @param aligned
2764    *          - 'guide' alignment containing sequences derived from same dataset
2765    *          as unaligned
2766    * @return
2767    */
2768   static boolean alignAsSameSequences(AlignmentI unaligned,
2769           AlignmentI aligned)
2770   {
2771     if (aligned.getDataset() == null || unaligned.getDataset() == null)
2772     {
2773       return false; // should only pass alignments with datasets here
2774     }
2775
2776     // map from dataset sequence to alignment sequence(s)
2777     Map<SequenceI, List<SequenceI>> alignedDatasets = new HashMap<SequenceI, List<SequenceI>>();
2778     for (SequenceI seq : aligned.getSequences())
2779     {
2780       SequenceI ds = seq.getDatasetSequence();
2781       if (alignedDatasets.get(ds) == null)
2782       {
2783         alignedDatasets.put(ds, new ArrayList<SequenceI>());
2784       }
2785       alignedDatasets.get(ds).add(seq);
2786     }
2787
2788     /*
2789      * first pass - check whether all sequences to be aligned share a dataset
2790      * sequence with an aligned sequence
2791      */
2792     for (SequenceI seq : unaligned.getSequences())
2793     {
2794       if (!alignedDatasets.containsKey(seq.getDatasetSequence()))
2795       {
2796         return false;
2797       }
2798     }
2799
2800     /*
2801      * second pass - copy aligned sequences;
2802      * heuristic rule: pair off sequences in order for the case where 
2803      * more than one shares the same dataset sequence 
2804      */
2805     for (SequenceI seq : unaligned.getSequences())
2806     {
2807       List<SequenceI> alignedSequences = alignedDatasets
2808               .get(seq.getDatasetSequence());
2809       // TODO: getSequenceAsString() will be deprecated in the future
2810       // TODO: need to leave to SequenceI implementor to update gaps
2811       seq.setSequence(alignedSequences.get(0).getSequenceAsString());
2812       if (alignedSequences.size() > 0)
2813       {
2814         // pop off aligned sequences (except the last one)
2815         alignedSequences.remove(0);
2816       }
2817     }
2818
2819     return true;
2820   }
2821
2822   /**
2823    * Returns a map whose key is alignment column number (base 1), and whose
2824    * values are a map of sequence characters in that column.
2825    * 
2826    * @param unaligned
2827    * @param aligned
2828    * @param unmapped
2829    * @return
2830    */
2831   static SortedMap<Integer, Map<SequenceI, Character>> buildMappedColumnsMap(
2832           AlignmentI unaligned, AlignmentI aligned,
2833           List<SequenceI> unmapped)
2834   {
2835     /*
2836      * Map will hold, for each aligned column position, a map of
2837      * {unalignedSequence, characterPerSequence} at that position.
2838      * TreeMap keeps the entries in ascending column order. 
2839      */
2840     SortedMap<Integer, Map<SequenceI, Character>> map = new TreeMap<Integer, Map<SequenceI, Character>>();
2841
2842     /*
2843      * record any sequences that have no mapping so can't be realigned
2844      */
2845     unmapped.addAll(unaligned.getSequences());
2846
2847     List<AlignedCodonFrame> mappings = aligned.getCodonFrames();
2848
2849     for (SequenceI seq : unaligned.getSequences())
2850     {
2851       for (AlignedCodonFrame mapping : mappings)
2852       {
2853         SequenceI fromSeq = mapping.findAlignedSequence(seq, aligned);
2854         if (fromSeq != null)
2855         {
2856           Mapping seqMap = mapping.getMappingBetween(fromSeq, seq);
2857           if (addMappedPositions(seq, fromSeq, seqMap, map))
2858           {
2859             unmapped.remove(seq);
2860           }
2861         }
2862       }
2863     }
2864     return map;
2865   }
2866
2867   /**
2868    * Helper method that adds to a map the mapped column positions of a sequence.
2869    * <br>
2870    * For example if aaTT-Tg-gAAA is mapped to TTTAAA then the map should record
2871    * that columns 3,4,6,10,11,12 map to characters T,T,T,A,A,A of the mapped to
2872    * sequence.
2873    * 
2874    * @param seq
2875    *          the sequence whose column positions we are recording
2876    * @param fromSeq
2877    *          a sequence that is mapped to the first sequence
2878    * @param seqMap
2879    *          the mapping from 'fromSeq' to 'seq'
2880    * @param map
2881    *          a map to add the column positions (in fromSeq) of the mapped
2882    *          positions of seq
2883    * @return
2884    */
2885   static boolean addMappedPositions(SequenceI seq, SequenceI fromSeq,
2886           Mapping seqMap, Map<Integer, Map<SequenceI, Character>> map)
2887   {
2888     if (seqMap == null)
2889     {
2890       return false;
2891     }
2892
2893     /*
2894      * invert mapping if it is from unaligned to aligned sequence
2895      */
2896     if (seqMap.getTo() == fromSeq.getDatasetSequence())
2897     {
2898       seqMap = new Mapping(seq.getDatasetSequence(),
2899               seqMap.getMap().getInverse());
2900     }
2901
2902     int toStart = seq.getStart();
2903
2904     /*
2905      * traverse [start, end, start, end...] ranges in fromSeq
2906      */
2907     for (int[] fromRange : seqMap.getMap().getFromRanges())
2908     {
2909       for (int i = 0; i < fromRange.length - 1; i += 2)
2910       {
2911         boolean forward = fromRange[i + 1] >= fromRange[i];
2912
2913         /*
2914          * find the range mapped to (sequence positions base 1)
2915          */
2916         int[] range = seqMap.locateMappedRange(fromRange[i],
2917                 fromRange[i + 1]);
2918         if (range == null)
2919         {
2920           System.err.println("Error in mapping " + seqMap + " from "
2921                   + fromSeq.getName());
2922           return false;
2923         }
2924         int fromCol = fromSeq.findIndex(fromRange[i]);
2925         int mappedCharPos = range[0];
2926
2927         /*
2928          * walk over the 'from' aligned sequence in forward or reverse
2929          * direction; when a non-gap is found, record the column position
2930          * of the next character of the mapped-to sequence; stop when all
2931          * the characters of the range have been counted
2932          */
2933         while (mappedCharPos <= range[1] && fromCol <= fromSeq.getLength()
2934                 && fromCol >= 0)
2935         {
2936           if (!Comparison.isGap(fromSeq.getCharAt(fromCol - 1)))
2937           {
2938             /*
2939              * mapped from sequence has a character in this column
2940              * record the column position for the mapped to character
2941              */
2942             Map<SequenceI, Character> seqsMap = map.get(fromCol);
2943             if (seqsMap == null)
2944             {
2945               seqsMap = new HashMap<SequenceI, Character>();
2946               map.put(fromCol, seqsMap);
2947             }
2948             seqsMap.put(seq, seq.getCharAt(mappedCharPos - toStart));
2949             mappedCharPos++;
2950           }
2951           fromCol += (forward ? 1 : -1);
2952         }
2953       }
2954     }
2955     return true;
2956   }
2957
2958   // strictly temporary hack until proper criteria for aligning protein to cds
2959   // are in place; this is so Ensembl -> fetch xrefs Uniprot aligns the Uniprot
2960   public static boolean looksLikeEnsembl(AlignmentI alignment)
2961   {
2962     for (SequenceI seq : alignment.getSequences())
2963     {
2964       String name = seq.getName();
2965       if (!name.startsWith("ENSG") && !name.startsWith("ENST"))
2966       {
2967         return false;
2968       }
2969     }
2970     return true;
2971   }
2972 }