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