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 jalview.datamodel.AlignedCodon;
  24 import jalview.datamodel.AlignedCodonFrame;
  25 import jalview.datamodel.Alignment;
  26 import jalview.datamodel.AlignmentAnnotation;
  27 import jalview.datamodel.AlignmentI;
  28 import jalview.datamodel.DBRefEntry;
  29 import jalview.datamodel.DBRefSource;
  30 import jalview.datamodel.FeatureProperties;
  31 import jalview.datamodel.IncompleteCodonException;
  32 import jalview.datamodel.Mapping;
  33 import jalview.datamodel.SearchResults;
  34 import jalview.datamodel.Sequence;
  35 import jalview.datamodel.SequenceFeature;
  36 import jalview.datamodel.SequenceGroup;
  37 import jalview.datamodel.SequenceI;
  38 import jalview.io.gff.SequenceOntologyFactory;
  39 import jalview.io.gff.SequenceOntologyI;
  40 import jalview.schemes.ResidueProperties;
  41 import jalview.util.Comparison;
  42 import jalview.util.DBRefUtils;
  43 import jalview.util.MapList;
  44 import jalview.util.MappingUtils;
  45
  46 import java.util.ArrayList;
  47 import java.util.Arrays;
  48 import java.util.Collection;
  49 import java.util.Collections;
  50 import java.util.Comparator;
  51 import java.util.HashMap;
  52 import java.util.HashSet;
  53 import java.util.Iterator;
  54 import java.util.LinkedHashMap;
  55 import java.util.List;
  56 import java.util.Map;
  57 import java.util.Map.Entry;
  58 import java.util.Set;
  59 import java.util.TreeMap;
  60
  61 /**
  62  * grab bag of useful alignment manipulation operations Expect these to be
  63  * refactored elsewhere at some point.
  64  *
  65  * @author jimp
  66  *
  67  */
  68 public class AlignmentUtils
  69 {
  70
  71   /**
  72    * given an existing alignment, create a new alignment including all, or up to
  73    * flankSize additional symbols from each sequence's dataset sequence
  74    *
  75    * @param core
  76    * @param flankSize
  77    * @return AlignmentI
  78    */
  79   public static AlignmentI expandContext(AlignmentI core, int flankSize)
  80   {
  81     List<SequenceI> sq = new ArrayList<SequenceI>();
  82     int maxoffset = 0;
  83     for (SequenceI s : core.getSequences())
  84     {
  85       SequenceI newSeq = s.deriveSequence();
  86       final int newSeqStart = newSeq.getStart() - 1;
  87       if (newSeqStart > maxoffset
  88               && newSeq.getDatasetSequence().getStart() < s.getStart())
  89       {
  90         maxoffset = newSeqStart;
  91       }
  92       sq.add(newSeq);
  93     }
  94     if (flankSize > -1)
  95     {
  96       maxoffset = Math.min(maxoffset, flankSize);
  97     }
  98
  99     /*
 100      * now add offset left and right to create an expanded alignment
 101      */
 102     for (SequenceI s : sq)
 103     {
 104       SequenceI ds = s;
 105       while (ds.getDatasetSequence() != null)
 106       {
 107         ds = ds.getDatasetSequence();
 108       }
 109       int s_end = s.findPosition(s.getStart() + s.getLength());
 110       // find available flanking residues for sequence
 111       int ustream_ds = s.getStart() - ds.getStart();
 112       int dstream_ds = ds.getEnd() - s_end;
 113
 114       // build new flanked sequence
 115
 116       // compute gap padding to start of flanking sequence
 117       int offset = maxoffset - ustream_ds;
 118
 119       // padding is gapChar x ( maxoffset - min(ustream_ds, flank)
 120       if (flankSize >= 0)
 121       {
 122         if (flankSize < ustream_ds)
 123         {
 124           // take up to flankSize residues
 125           offset = maxoffset - flankSize;
 126           ustream_ds = flankSize;
 127         }
 128         if (flankSize <= dstream_ds)
 129         {
 130           dstream_ds = flankSize - 1;
 131         }
 132       }
 133       // TODO use Character.toLowerCase to avoid creating String objects?
 134       char[] upstream = new String(ds.getSequence(s.getStart() - 1
 135               - ustream_ds, s.getStart() - 1)).toLowerCase().toCharArray();
 136       char[] downstream = new String(ds.getSequence(s_end - 1, s_end
 137               + dstream_ds)).toLowerCase().toCharArray();
 138       char[] coreseq = s.getSequence();
 139       char[] nseq = new char[offset + upstream.length + downstream.length
 140               + coreseq.length];
 141       char c = core.getGapCharacter();
 142
 143       int p = 0;
 144       for (; p < offset; p++)
 145       {
 146         nseq[p] = c;
 147       }
 148
 149       System.arraycopy(upstream, 0, nseq, p, upstream.length);
 150       System.arraycopy(coreseq, 0, nseq, p + upstream.length,
 151               coreseq.length);
 152       System.arraycopy(downstream, 0, nseq, p + coreseq.length
 153               + upstream.length, downstream.length);
 154       s.setSequence(new String(nseq));
 155       s.setStart(s.getStart() - ustream_ds);
 156       s.setEnd(s_end + downstream.length);
 157     }
 158     AlignmentI newAl = new jalview.datamodel.Alignment(
 159             sq.toArray(new SequenceI[0]));
 160     for (SequenceI s : sq)
 161     {
 162       if (s.getAnnotation() != null)
 163       {
 164         for (AlignmentAnnotation aa : s.getAnnotation())
 165         {
 166           aa.adjustForAlignment(); // JAL-1712 fix
 167           newAl.addAnnotation(aa);
 168         }
 169       }
 170     }
 171     newAl.setDataset(core.getDataset());
 172     return newAl;
 173   }
 174
 175   /**
 176    * Returns the index (zero-based position) of a sequence in an alignment, or
 177    * -1 if not found.
 178    *
 179    * @param al
 180    * @param seq
 181    * @return
 182    */
 183   public static int getSequenceIndex(AlignmentI al, SequenceI seq)
 184   {
 185     int result = -1;
 186     int pos = 0;
 187     for (SequenceI alSeq : al.getSequences())
 188     {
 189       if (alSeq == seq)
 190       {
 191         result = pos;
 192         break;
 193       }
 194       pos++;
 195     }
 196     return result;
 197   }
 198
 199   /**
 200    * Returns a map of lists of sequences in the alignment, keyed by sequence
 201    * name. For use in mapping between different alignment views of the same
 202    * sequences.
 203    *
 204    * @see jalview.datamodel.AlignmentI#getSequencesByName()
 205    */
 206   public static Map<String, List<SequenceI>> getSequencesByName(
 207           AlignmentI al)
 208   {
 209     Map<String, List<SequenceI>> theMap = new LinkedHashMap<String, List<SequenceI>>();
 210     for (SequenceI seq : al.getSequences())
 211     {
 212       String name = seq.getName();
 213       if (name != null)
 214       {
 215         List<SequenceI> seqs = theMap.get(name);
 216         if (seqs == null)
 217         {
 218           seqs = new ArrayList<SequenceI>();
 219           theMap.put(name, seqs);
 220         }
 221         seqs.add(seq);
 222       }
 223     }
 224     return theMap;
 225   }
 226
 227   /**
 228    * Build mapping of protein to cDNA alignment. Mappings are made between
 229    * sequences where the cDNA translates to the protein sequence. Any new
 230    * mappings are added to the protein alignment. Returns true if any mappings
 231    * either already exist or were added, else false.
 232    *
 233    * @param proteinAlignment
 234    * @param cdnaAlignment
 235    * @return
 236    */
 237   public static boolean mapProteinAlignmentToCdna(
 238           final AlignmentI proteinAlignment, final AlignmentI cdnaAlignment)
 239   {
 240     if (proteinAlignment == null || cdnaAlignment == null)
 241     {
 242       return false;
 243     }
 244
 245     Set<SequenceI> mappedDna = new HashSet<SequenceI>();
 246     Set<SequenceI> mappedProtein = new HashSet<SequenceI>();
 247
 248     /*
 249      * First pass - map sequences where cross-references exist. This include
 250      * 1-to-many mappings to support, for example, variant cDNA.
 251      */
 252     boolean mappingPerformed = mapProteinToCdna(proteinAlignment,
 253             cdnaAlignment, mappedDna, mappedProtein, true);
 254
 255     /*
 256      * Second pass - map sequences where no cross-references exist. This only
 257      * does 1-to-1 mappings and assumes corresponding sequences are in the same
 258      * order in the alignments.
 259      */
 260     mappingPerformed |= mapProteinToCdna(proteinAlignment, cdnaAlignment,
 261             mappedDna, mappedProtein, false);
 262     return mappingPerformed;
 263   }
 264
 265   /**
 266    * Make mappings between compatible sequences (where the cDNA translation
 267    * matches the protein).
 268    *
 269    * @param proteinAlignment
 270    * @param cdnaAlignment
 271    * @param mappedDna
 272    *          a set of mapped DNA sequences (to add to)
 273    * @param mappedProtein
 274    *          a set of mapped Protein sequences (to add to)
 275    * @param xrefsOnly
 276    *          if true, only map sequences where xrefs exist
 277    * @return
 278    */
 279   protected static boolean mapProteinToCdna(
 280           final AlignmentI proteinAlignment,
 281           final AlignmentI cdnaAlignment, Set<SequenceI> mappedDna,
 282           Set<SequenceI> mappedProtein, boolean xrefsOnly)
 283   {
 284     boolean mappingExistsOrAdded = false;
 285     List<SequenceI> thisSeqs = proteinAlignment.getSequences();
 286     for (SequenceI aaSeq : thisSeqs)
 287     {
 288       boolean proteinMapped = false;
 289       AlignedCodonFrame acf = new AlignedCodonFrame();
 290
 291       for (SequenceI cdnaSeq : cdnaAlignment.getSequences())
 292       {
 293         /*
 294          * Always try to map if sequences have xref to each other; this supports
 295          * variant cDNA or alternative splicing for a protein sequence.
 296          *
 297          * If no xrefs, try to map progressively, assuming that alignments have
 298          * mappable sequences in corresponding order. These are not
 299          * many-to-many, as that would risk mixing species with similar cDNA
 300          * sequences.
 301          */
 302         if (xrefsOnly && !AlignmentUtils.haveCrossRef(aaSeq, cdnaSeq))
 303         {
 304           continue;
 305         }
 306
 307         /*
 308          * Don't map non-xrefd sequences more than once each. This heuristic
 309          * allows us to pair up similar sequences in ordered alignments.
 310          */
 311         if (!xrefsOnly
 312                 && (mappedProtein.contains(aaSeq) || mappedDna
 313                         .contains(cdnaSeq)))
 314         {
 315           continue;
 316         }
 317         if (mappingExists(proteinAlignment.getCodonFrames(),
 318                 aaSeq.getDatasetSequence(), cdnaSeq.getDatasetSequence()))
 319         {
 320           mappingExistsOrAdded = true;
 321         }
 322         else
 323         {
 324           MapList map = mapProteinSequenceToCdna(aaSeq, cdnaSeq);
 325           if (map != null)
 326           {
 327             acf.addMap(cdnaSeq, aaSeq, map);
 328             mappingExistsOrAdded = true;
 329             proteinMapped = true;
 330             mappedDna.add(cdnaSeq);
 331             mappedProtein.add(aaSeq);
 332           }
 333         }
 334       }
 335       if (proteinMapped)
 336       {
 337         proteinAlignment.addCodonFrame(acf);
 338       }
 339     }
 340     return mappingExistsOrAdded;
 341   }
 342
 343   /**
 344    * Answers true if the mappings include one between the given (dataset)
 345    * sequences.
 346    */
 347   public static boolean mappingExists(List<AlignedCodonFrame> mappings,
 348           SequenceI aaSeq, SequenceI cdnaSeq)
 349   {
 350     if (mappings != null)
 351     {
 352       for (AlignedCodonFrame acf : mappings)
 353       {
 354         if (cdnaSeq == acf.getDnaForAaSeq(aaSeq))
 355         {
 356           return true;
 357         }
 358       }
 359     }
 360     return false;
 361   }
 362
 363   /**
 364    * Build a mapping (if possible) of a protein to a cDNA sequence. The cDNA
 365    * must be three times the length of the protein, possibly after ignoring
 366    * start and/or stop codons, and must translate to the protein. Returns null
 367    * if no mapping is determined.
 368    *
 369    * @param proteinSeqs
 370    * @param cdnaSeq
 371    * @return
 372    */
 373   public static MapList mapProteinSequenceToCdna(SequenceI proteinSeq,
 374           SequenceI cdnaSeq)
 375   {
 376     /*
 377      * Here we handle either dataset sequence set (desktop) or absent (applet).
 378      * Use only the char[] form of the sequence to avoid creating possibly large
 379      * String objects.
 380      */
 381     final SequenceI proteinDataset = proteinSeq.getDatasetSequence();
 382     char[] aaSeqChars = proteinDataset != null ? proteinDataset
 383             .getSequence() : proteinSeq.getSequence();
 384     final SequenceI cdnaDataset = cdnaSeq.getDatasetSequence();
 385     char[] cdnaSeqChars = cdnaDataset != null ? cdnaDataset.getSequence()
 386             : cdnaSeq.getSequence();
 387     if (aaSeqChars == null || cdnaSeqChars == null)
 388     {
 389       return null;
 390     }
 391
 392     /*
 393      * cdnaStart/End, proteinStartEnd are base 1 (for dataset sequence mapping)
 394      */
 395     final int mappedLength = 3 * aaSeqChars.length;
 396     int cdnaLength = cdnaSeqChars.length;
 397     int cdnaStart = cdnaSeq.getStart();
 398     int cdnaEnd = cdnaSeq.getEnd();
 399     final int proteinStart = proteinSeq.getStart();
 400     final int proteinEnd = proteinSeq.getEnd();
 401
 402     /*
 403      * If lengths don't match, try ignoring stop codon.
 404      */
 405     if (cdnaLength != mappedLength && cdnaLength > 2)
 406     {
 407       String lastCodon = String.valueOf(cdnaSeqChars, cdnaLength - 3, 3)
 408               .toUpperCase();
 409       for (String stop : ResidueProperties.STOP)
 410       {
 411         if (lastCodon.equals(stop))
 412         {
 413           cdnaEnd -= 3;
 414           cdnaLength -= 3;
 415           break;
 416         }
 417       }
 418     }
 419
 420     /*
 421      * If lengths still don't match, try ignoring start codon.
 422      */
 423     int startOffset = 0;
 424     if (cdnaLength != mappedLength
 425             && cdnaLength > 2
 426             && String.valueOf(cdnaSeqChars, 0, 3).toUpperCase()
 427                     .equals(ResidueProperties.START))
 428     {
 429       startOffset += 3;
 430       cdnaStart += 3;
 431       cdnaLength -= 3;
 432     }
 433
 434     if (cdnaLength != mappedLength)
 435     {
 436       return null;
 437     }
 438     if (!translatesAs(cdnaSeqChars, startOffset, aaSeqChars))
 439     {
 440       return null;
 441     }
 442     MapList map = new MapList(new int[] { cdnaStart, cdnaEnd }, new int[] {
 443         proteinStart, proteinEnd }, 3, 1);
 444     return map;
 445   }
 446
 447   /**
 448    * Test whether the given cdna sequence, starting at the given offset,
 449    * translates to the given amino acid sequence, using the standard translation
 450    * table. Designed to fail fast i.e. as soon as a mismatch position is found.
 451    *
 452    * @param cdnaSeqChars
 453    * @param cdnaStart
 454    * @param aaSeqChars
 455    * @return
 456    */
 457   protected static boolean translatesAs(char[] cdnaSeqChars, int cdnaStart,
 458           char[] aaSeqChars)
 459   {
 460     if (cdnaSeqChars == null || aaSeqChars == null)
 461     {
 462       return false;
 463     }
 464
 465     int aaResidue = 0;
 466     for (int i = cdnaStart; i < cdnaSeqChars.length - 2
 467             && aaResidue < aaSeqChars.length; i += 3, aaResidue++)
 468     {
 469       String codon = String.valueOf(cdnaSeqChars, i, 3);
 470       final String translated = ResidueProperties.codonTranslate(codon);
 471       /*
 472        * allow * in protein to match untranslatable in dna
 473        */
 474       final char aaRes = aaSeqChars[aaResidue];
 475       if ((translated == null || "STOP".equals(translated)) && aaRes == '*')
 476       {
 477         continue;
 478       }
 479       if (translated == null || !(aaRes == translated.charAt(0)))
 480       {
 481         // debug
 482         // System.out.println(("Mismatch at " + i + "/" + aaResidue + ": "
 483         // + codon + "(" + translated + ") != " + aaRes));
 484         return false;
 485       }
 486     }
 487     // fail if we didn't match all of the aa sequence
 488     return (aaResidue == aaSeqChars.length);
 489   }
 490
 491   /**
 492    * Align sequence 'seq' to match the alignment of a mapped sequence. Note this
 493    * currently assumes that we are aligning cDNA to match protein.
 494    *
 495    * @param seq
 496    *          the sequence to be realigned
 497    * @param al
 498    *          the alignment whose sequence alignment is to be 'copied'
 499    * @param gap
 500    *          character string represent a gap in the realigned sequence
 501    * @param preserveUnmappedGaps
 502    * @param preserveMappedGaps
 503    * @return true if the sequence was realigned, false if it could not be
 504    */
 505   public static boolean alignSequenceAs(SequenceI seq, AlignmentI al,
 506           String gap, boolean preserveMappedGaps,
 507           boolean preserveUnmappedGaps)
 508   {
 509     /*
 510      * Get any mappings from the source alignment to the target (dataset)
 511      * sequence.
 512      */
 513     // TODO there may be one AlignedCodonFrame per dataset sequence, or one with
 514     // all mappings. Would it help to constrain this?
 515     List<AlignedCodonFrame> mappings = al.getCodonFrame(seq);
 516     if (mappings == null || mappings.isEmpty())
 517     {
 518       return false;
 519     }
 520
 521     /*
 522      * Locate the aligned source sequence whose dataset sequence is mapped. We
 523      * just take the first match here (as we can't align like more than one
 524      * sequence).
 525      */
 526     SequenceI alignFrom = null;
 527     AlignedCodonFrame mapping = null;
 528     for (AlignedCodonFrame mp : mappings)
 529     {
 530       alignFrom = mp.findAlignedSequence(seq.getDatasetSequence(), al);
 531       if (alignFrom != null)
 532       {
 533         mapping = mp;
 534         break;
 535       }
 536     }
 537
 538     if (alignFrom == null)
 539     {
 540       return false;
 541     }
 542     alignSequenceAs(seq, alignFrom, mapping, gap, al.getGapCharacter(),
 543             preserveMappedGaps, preserveUnmappedGaps);
 544     return true;
 545   }
 546
 547   /**
 548    * Align sequence 'alignTo' the same way as 'alignFrom', using the mapping to
 549    * match residues and codons. Flags control whether existing gaps in unmapped
 550    * (intron) and mapped (exon) regions are preserved or not. Gaps between
 551    * intron and exon are only retained if both flags are set.
 552    *
 553    * @param alignTo
 554    * @param alignFrom
 555    * @param mapping
 556    * @param myGap
 557    * @param sourceGap
 558    * @param preserveUnmappedGaps
 559    * @param preserveMappedGaps
 560    */
 561   public static void alignSequenceAs(SequenceI alignTo,
 562           SequenceI alignFrom, AlignedCodonFrame mapping, String myGap,
 563           char sourceGap, boolean preserveMappedGaps,
 564           boolean preserveUnmappedGaps)
 565   {
 566     // TODO generalise to work for Protein-Protein, dna-dna, dna-protein
 567
 568     // aligned and dataset sequence positions, all base zero
 569     int thisSeqPos = 0;
 570     int sourceDsPos = 0;
 571
 572     int basesWritten = 0;
 573     char myGapChar = myGap.charAt(0);
 574     int ratio = myGap.length();
 575
 576     int fromOffset = alignFrom.getStart() - 1;
 577     int toOffset = alignTo.getStart() - 1;
 578     int sourceGapMappedLength = 0;
 579     boolean inExon = false;
 580     final char[] thisSeq = alignTo.getSequence();
 581     final char[] thatAligned = alignFrom.getSequence();
 582     StringBuilder thisAligned = new StringBuilder(2 * thisSeq.length);
 583
 584     /*
 585      * Traverse the 'model' aligned sequence
 586      */
 587     for (char sourceChar : thatAligned)
 588     {
 589       if (sourceChar == sourceGap)
 590       {
 591         sourceGapMappedLength += ratio;
 592         continue;
 593       }
 594
 595       /*
 596        * Found a non-gap character. Locate its mapped region if any.
 597        */
 598       sourceDsPos++;
 599       // Note mapping positions are base 1, our sequence positions base 0
 600       int[] mappedPos = mapping.getMappedRegion(alignTo, alignFrom,
 601               sourceDsPos + fromOffset);
 602       if (mappedPos == null)
 603       {
 604         /*
 605          * unmapped position; treat like a gap
 606          */
 607         sourceGapMappedLength += ratio;
 608         // System.err.println("Can't align: no codon mapping to residue "
 609         // + sourceDsPos + "(" + sourceChar + ")");
 610         // return;
 611         continue;
 612       }
 613
 614       int mappedCodonStart = mappedPos[0]; // position (1...) of codon start
 615       int mappedCodonEnd = mappedPos[mappedPos.length - 1]; // codon end pos
 616       StringBuilder trailingCopiedGap = new StringBuilder();
 617
 618       /*
 619        * Copy dna sequence up to and including this codon. Optionally, include
 620        * gaps before the codon starts (in introns) and/or after the codon starts
 621        * (in exons).
 622        *
 623        * Note this only works for 'linear' splicing, not reverse or interleaved.
 624        * But then 'align dna as protein' doesn't make much sense otherwise.
 625        */
 626       int intronLength = 0;
 627       while (basesWritten + toOffset < mappedCodonEnd
 628               && thisSeqPos < thisSeq.length)
 629       {
 630         final char c = thisSeq[thisSeqPos++];
 631         if (c != myGapChar)
 632         {
 633           basesWritten++;
 634           int sourcePosition = basesWritten + toOffset;
 635           if (sourcePosition < mappedCodonStart)
 636           {
 637             /*
 638              * Found an unmapped (intron) base. First add in any preceding gaps
 639              * (if wanted).
 640              */
 641             if (preserveUnmappedGaps && trailingCopiedGap.length() > 0)
 642             {
 643               thisAligned.append(trailingCopiedGap.toString());
 644               intronLength += trailingCopiedGap.length();
 645               trailingCopiedGap = new StringBuilder();
 646             }
 647             intronLength++;
 648             inExon = false;
 649           }
 650           else
 651           {
 652             final boolean startOfCodon = sourcePosition == mappedCodonStart;
 653             int gapsToAdd = calculateGapsToInsert(preserveMappedGaps,
 654                     preserveUnmappedGaps, sourceGapMappedLength, inExon,
 655                     trailingCopiedGap.length(), intronLength, startOfCodon);
 656             for (int i = 0; i < gapsToAdd; i++)
 657             {
 658               thisAligned.append(myGapChar);
 659             }
 660             sourceGapMappedLength = 0;
 661             inExon = true;
 662           }
 663           thisAligned.append(c);
 664           trailingCopiedGap = new StringBuilder();
 665         }
 666         else
 667         {
 668           if (inExon && preserveMappedGaps)
 669           {
 670             trailingCopiedGap.append(myGapChar);
 671           }
 672           else if (!inExon && preserveUnmappedGaps)
 673           {
 674             trailingCopiedGap.append(myGapChar);
 675           }
 676         }
 677       }
 678     }
 679
 680     /*
 681      * At end of model aligned sequence. Copy any remaining target sequence, optionally
 682      * including (intron) gaps.
 683      */
 684     while (thisSeqPos < thisSeq.length)
 685     {
 686       final char c = thisSeq[thisSeqPos++];
 687       if (c != myGapChar || preserveUnmappedGaps)
 688       {
 689         thisAligned.append(c);
 690       }
 691       sourceGapMappedLength--;
 692     }
 693
 694     /*
 695      * finally add gaps to pad for any trailing source gaps or
 696      * unmapped characters
 697      */
 698     if (preserveUnmappedGaps)
 699     {
 700       while (sourceGapMappedLength > 0)
 701       {
 702         thisAligned.append(myGapChar);
 703         sourceGapMappedLength--;
 704       }
 705     }
 706
 707     /*
 708      * All done aligning, set the aligned sequence.
 709      */
 710     alignTo.setSequence(new String(thisAligned));
 711   }
 712
 713   /**
 714    * Helper method to work out how many gaps to insert when realigning.
 715    *
 716    * @param preserveMappedGaps
 717    * @param preserveUnmappedGaps
 718    * @param sourceGapMappedLength
 719    * @param inExon
 720    * @param trailingCopiedGap
 721    * @param intronLength
 722    * @param startOfCodon
 723    * @return
 724    */
 725   protected static int calculateGapsToInsert(boolean preserveMappedGaps,
 726           boolean preserveUnmappedGaps, int sourceGapMappedLength,
 727           boolean inExon, int trailingGapLength, int intronLength,
 728           final boolean startOfCodon)
 729   {
 730     int gapsToAdd = 0;
 731     if (startOfCodon)
 732     {
 733       /*
 734        * Reached start of codon. Ignore trailing gaps in intron unless we are
 735        * preserving gaps in both exon and intron. Ignore them anyway if the
 736        * protein alignment introduces a gap at least as large as the intronic
 737        * region.
 738        */
 739       if (inExon && !preserveMappedGaps)
 740       {
 741         trailingGapLength = 0;
 742       }
 743       if (!inExon && !(preserveMappedGaps && preserveUnmappedGaps))
 744       {
 745         trailingGapLength = 0;
 746       }
 747       if (inExon)
 748       {
 749         gapsToAdd = Math.max(sourceGapMappedLength, trailingGapLength);
 750       }
 751       else
 752       {
 753         if (intronLength + trailingGapLength <= sourceGapMappedLength)
 754         {
 755           gapsToAdd = sourceGapMappedLength - intronLength;
 756         }
 757         else
 758         {
 759           gapsToAdd = Math.min(intronLength + trailingGapLength
 760                   - sourceGapMappedLength, trailingGapLength);
 761         }
 762       }
 763     }
 764     else
 765     {
 766       /*
 767        * second or third base of codon; check for any gaps in dna
 768        */
 769       if (!preserveMappedGaps)
 770       {
 771         trailingGapLength = 0;
 772       }
 773       gapsToAdd = Math.max(sourceGapMappedLength, trailingGapLength);
 774     }
 775     return gapsToAdd;
 776   }
 777
 778   /**
 779    * Returns a list of sequences mapped from the given sequences and aligned
 780    * (gapped) in the same way. For example, the cDNA for aligned protein, where
 781    * a single gap in protein generates three gaps in cDNA.
 782    *
 783    * @param sequences
 784    * @param gapCharacter
 785    * @param mappings
 786    * @return
 787    */
 788   public static List<SequenceI> getAlignedTranslation(
 789           List<SequenceI> sequences, char gapCharacter,
 790           Set<AlignedCodonFrame> mappings)
 791   {
 792     List<SequenceI> alignedSeqs = new ArrayList<SequenceI>();
 793
 794     for (SequenceI seq : sequences)
 795     {
 796       List<SequenceI> mapped = getAlignedTranslation(seq, gapCharacter,
 797               mappings);
 798       alignedSeqs.addAll(mapped);
 799     }
 800     return alignedSeqs;
 801   }
 802
 803   /**
 804    * Returns sequences aligned 'like' the source sequence, as mapped by the
 805    * given mappings. Normally we expect zero or one 'mapped' sequences, but this
 806    * will support 1-to-many as well.
 807    *
 808    * @param seq
 809    * @param gapCharacter
 810    * @param mappings
 811    * @return
 812    */
 813   protected static List<SequenceI> getAlignedTranslation(SequenceI seq,
 814           char gapCharacter, Set<AlignedCodonFrame> mappings)
 815   {
 816     List<SequenceI> result = new ArrayList<SequenceI>();
 817     for (AlignedCodonFrame mapping : mappings)
 818     {
 819       if (mapping.involvesSequence(seq))
 820       {
 821         SequenceI mapped = getAlignedTranslation(seq, gapCharacter, mapping);
 822         if (mapped != null)
 823         {
 824           result.add(mapped);
 825         }
 826       }
 827     }
 828     return result;
 829   }
 830
 831   /**
 832    * Returns the translation of 'seq' (as held in the mapping) with
 833    * corresponding alignment (gaps).
 834    *
 835    * @param seq
 836    * @param gapCharacter
 837    * @param mapping
 838    * @return
 839    */
 840   protected static SequenceI getAlignedTranslation(SequenceI seq,
 841           char gapCharacter, AlignedCodonFrame mapping)
 842   {
 843     String gap = String.valueOf(gapCharacter);
 844     boolean toDna = false;
 845     int fromRatio = 1;
 846     SequenceI mapTo = mapping.getDnaForAaSeq(seq);
 847     if (mapTo != null)
 848     {
 849       // mapping is from protein to nucleotide
 850       toDna = true;
 851       // should ideally get gap count ratio from mapping
 852       gap = String.valueOf(new char[] { gapCharacter, gapCharacter,
 853           gapCharacter });
 854     }
 855     else
 856     {
 857       // mapping is from nucleotide to protein
 858       mapTo = mapping.getAaForDnaSeq(seq);
 859       fromRatio = 3;
 860     }
 861     StringBuilder newseq = new StringBuilder(seq.getLength()
 862             * (toDna ? 3 : 1));
 863
 864     int residueNo = 0; // in seq, base 1
 865     int[] phrase = new int[fromRatio];
 866     int phraseOffset = 0;
 867     int gapWidth = 0;
 868     boolean first = true;
 869     final Sequence alignedSeq = new Sequence("", "");
 870
 871     for (char c : seq.getSequence())
 872     {
 873       if (c == gapCharacter)
 874       {
 875         gapWidth++;
 876         if (gapWidth >= fromRatio)
 877         {
 878           newseq.append(gap);
 879           gapWidth = 0;
 880         }
 881       }
 882       else
 883       {
 884         phrase[phraseOffset++] = residueNo + 1;
 885         if (phraseOffset == fromRatio)
 886         {
 887           /*
 888            * Have read a whole codon (or protein residue), now translate: map
 889            * source phrase to positions in target sequence add characters at
 890            * these positions to newseq Note mapping positions are base 1, our
 891            * sequence positions base 0.
 892            */
 893           SearchResults sr = new SearchResults();
 894           for (int pos : phrase)
 895           {
 896             mapping.markMappedRegion(seq, pos, sr);
 897           }
 898           newseq.append(sr.getCharacters());
 899           if (first)
 900           {
 901             first = false;
 902             // Hack: Copy sequence dataset, name and description from
 903             // SearchResults.match[0].sequence
 904             // TODO? carry over sequence names from original 'complement'
 905             // alignment
 906             SequenceI mappedTo = sr.getResultSequence(0);
 907             alignedSeq.setName(mappedTo.getName());
 908             alignedSeq.setDescription(mappedTo.getDescription());
 909             alignedSeq.setDatasetSequence(mappedTo);
 910           }
 911           phraseOffset = 0;
 912         }
 913         residueNo++;
 914       }
 915     }
 916     alignedSeq.setSequence(newseq.toString());
 917     return alignedSeq;
 918   }
 919
 920   /**
 921    * Realigns the given protein to match the alignment of the dna, using codon
 922    * mappings to translate aligned codon positions to protein residues.
 923    *
 924    * @param protein
 925    *          the alignment whose sequences are realigned by this method
 926    * @param dna
 927    *          the dna alignment whose alignment we are 'copying'
 928    * @return the number of sequences that were realigned
 929    */
 930   public static int alignProteinAsDna(AlignmentI protein, AlignmentI dna)
 931   {
 932     List<SequenceI> unmappedProtein = new ArrayList<SequenceI>();
 933     Map<AlignedCodon, Map<SequenceI, AlignedCodon>> alignedCodons = buildCodonColumnsMap(
 934             protein, dna, unmappedProtein);
 935     return alignProteinAs(protein, alignedCodons, unmappedProtein);
 936   }
 937
 938   /**
 939    * Builds a map whose key is an aligned codon position (3 alignment column
 940    * numbers base 0), and whose value is a map from protein sequence to each
 941    * protein's peptide residue for that codon. The map generates an ordering of
 942    * the codons, and allows us to read off the peptides at each position in
 943    * order to assemble 'aligned' protein sequences.
 944    *
 945    * @param protein
 946    *          the protein alignment
 947    * @param dna
 948    *          the coding dna alignment
 949    * @param unmappedProtein
 950    *          any unmapped proteins are added to this list
 951    * @return
 952    */
 953   protected static Map<AlignedCodon, Map<SequenceI, AlignedCodon>> buildCodonColumnsMap(
 954           AlignmentI protein, AlignmentI dna,
 955           List<SequenceI> unmappedProtein)
 956   {
 957     /*
 958      * maintain a list of any proteins with no mappings - these will be
 959      * rendered 'as is' in the protein alignment as we can't align them
 960      */
 961     unmappedProtein.addAll(protein.getSequences());
 962
 963     List<AlignedCodonFrame> mappings = protein.getCodonFrames();
 964
 965     /*
 966      * Map will hold, for each aligned codon position e.g. [3, 5, 6], a map of
 967      * {dnaSequence, {proteinSequence, codonProduct}} at that position. The
 968      * comparator keeps the codon positions ordered.
 969      */
 970     Map<AlignedCodon, Map<SequenceI, AlignedCodon>> alignedCodons = new TreeMap<AlignedCodon, Map<SequenceI, AlignedCodon>>(
 971             new CodonComparator());
 972
 973     for (SequenceI dnaSeq : dna.getSequences())
 974     {
 975       for (AlignedCodonFrame mapping : mappings)
 976       {
 977         SequenceI prot = mapping.findAlignedSequence(
 978                 dnaSeq.getDatasetSequence(), protein);
 979         if (prot != null)
 980         {
 981           Mapping seqMap = mapping.getMappingForSequence(dnaSeq);
 982           addCodonPositions(dnaSeq, prot, protein.getGapCharacter(),
 983                   seqMap, alignedCodons);
 984           unmappedProtein.remove(prot);
 985         }
 986       }
 987     }
 988
 989     /*
 990      * Finally add any unmapped peptide start residues (e.g. for incomplete
 991      * codons) as if at the codon position before the second residue
 992      */
 993     int mappedSequenceCount = protein.getHeight() - unmappedProtein.size();
 994     addUnmappedPeptideStarts(alignedCodons, mappedSequenceCount);
 995
 996     return alignedCodons;
 997   }
 998
 999   /**
1000    * Scans for any protein mapped from position 2 (meaning unmapped start
1001    * position e.g. an incomplete codon), and synthesizes a 'codon' for it at the
1002    * preceding position in the alignment
1003    *
1004    * @param alignedCodons
1005    *          the codon-to-peptide map
1006    * @param mappedSequenceCount
1007    *          the number of distinct sequences in the map
1008    */
1009   protected static void addUnmappedPeptideStarts(
1010           Map<AlignedCodon, Map<SequenceI, AlignedCodon>> alignedCodons,
1011           int mappedSequenceCount)
1012   {
1013     // TODO there must be an easier way! root problem is that our mapping data
1014     // model does not include phase so can't map part of a codon to a peptide
1015     List<SequenceI> sequencesChecked = new ArrayList<SequenceI>();
1016     AlignedCodon lastCodon = null;
1017     Map<SequenceI, AlignedCodon> toAdd = new HashMap<SequenceI, AlignedCodon>();
1018
1019     for (Entry<AlignedCodon, Map<SequenceI, AlignedCodon>> entry : alignedCodons
1020             .entrySet())
1021     {
1022       for (Entry<SequenceI, AlignedCodon> sequenceCodon : entry.getValue()
1023               .entrySet())
1024       {
1025         SequenceI seq = sequenceCodon.getKey();
1026         if (sequencesChecked.contains(seq))
1027         {
1028           continue;
1029         }
1030         sequencesChecked.add(seq);
1031         AlignedCodon codon = sequenceCodon.getValue();
1032         if (codon.peptideCol > 1)
1033         {
1034           System.err
1035                   .println("Problem mapping protein with >1 unmapped start positions: "
1036                           + seq.getName());
1037         }
1038         else if (codon.peptideCol == 1)
1039         {
1040           /*
1041            * first position (peptideCol == 0) was unmapped - add it
1042            */
1043           if (lastCodon != null)
1044           {
1045             AlignedCodon firstPeptide = new AlignedCodon(lastCodon.pos1,
1046                     lastCodon.pos2, lastCodon.pos3, String.valueOf(seq
1047                             .getCharAt(0)), 0);
1048             toAdd.put(seq, firstPeptide);
1049           }
1050           else
1051           {
1052             /*
1053              * unmapped residue at start of alignment (no prior column) -
1054              * 'insert' at nominal codon [0, 0, 0]
1055              */
1056             AlignedCodon firstPeptide = new AlignedCodon(0, 0, 0,
1057                     String.valueOf(seq.getCharAt(0)), 0);
1058             toAdd.put(seq, firstPeptide);
1059           }
1060         }
1061         if (sequencesChecked.size() == mappedSequenceCount)
1062         {
1063           // no need to check past first mapped position in all sequences
1064           break;
1065         }
1066       }
1067       lastCodon = entry.getKey();
1068     }
1069
1070     /*
1071      * add any new codons safely after iterating over the map
1072      */
1073     for (Entry<SequenceI, AlignedCodon> startCodon : toAdd.entrySet())
1074     {
1075       addCodonToMap(alignedCodons, startCodon.getValue(),
1076               startCodon.getKey());
1077     }
1078   }
1079
1080   /**
1081    * Update the aligned protein sequences to match the codon alignments given in
1082    * the map.
1083    *
1084    * @param protein
1085    * @param alignedCodons
1086    *          an ordered map of codon positions (columns), with sequence/peptide
1087    *          values present in each column
1088    * @param unmappedProtein
1089    * @return
1090    */
1091   protected static int alignProteinAs(AlignmentI protein,
1092           Map<AlignedCodon, Map<SequenceI, AlignedCodon>> alignedCodons,
1093           List<SequenceI> unmappedProtein)
1094   {
1095     /*
1096      * Prefill aligned sequences with gaps before inserting aligned protein
1097      * residues.
1098      */
1099     int alignedWidth = alignedCodons.size();
1100     char[] gaps = new char[alignedWidth];
1101     Arrays.fill(gaps, protein.getGapCharacter());
1102     String allGaps = String.valueOf(gaps);
1103     for (SequenceI seq : protein.getSequences())
1104     {
1105       if (!unmappedProtein.contains(seq))
1106       {
1107         seq.setSequence(allGaps);
1108       }
1109     }
1110
1111     int column = 0;
1112     for (AlignedCodon codon : alignedCodons.keySet())
1113     {
1114       final Map<SequenceI, AlignedCodon> columnResidues = alignedCodons
1115               .get(codon);
1116       for (Entry<SequenceI, AlignedCodon> entry : columnResidues.entrySet())
1117       {
1118         // place translated codon at its column position in sequence
1119         entry.getKey().getSequence()[column] = entry.getValue().product
1120                 .charAt(0);
1121       }
1122       column++;
1123     }
1124     return 0;
1125   }
1126
1127   /**
1128    * Populate the map of aligned codons by traversing the given sequence
1129    * mapping, locating the aligned positions of mapped codons, and adding those
1130    * positions and their translation products to the map.
1131    *
1132    * @param dna
1133    *          the aligned sequence we are mapping from
1134    * @param protein
1135    *          the sequence to be aligned to the codons
1136    * @param gapChar
1137    *          the gap character in the dna sequence
1138    * @param seqMap
1139    *          a mapping to a sequence translation
1140    * @param alignedCodons
1141    *          the map we are building up
1142    */
1143   static void addCodonPositions(SequenceI dna, SequenceI protein,
1144           char gapChar, Mapping seqMap,
1145           Map<AlignedCodon, Map<SequenceI, AlignedCodon>> alignedCodons)
1146   {
1147     Iterator<AlignedCodon> codons = seqMap.getCodonIterator(dna, gapChar);
1148
1149     /*
1150      * add codon positions, and their peptide translations, to the alignment
1151      * map, while remembering the first codon mapped
1152      */
1153     while (codons.hasNext())
1154     {
1155       try
1156       {
1157         AlignedCodon codon = codons.next();
1158         addCodonToMap(alignedCodons, codon, protein);
1159       } catch (IncompleteCodonException e)
1160       {
1161         // possible incomplete trailing codon - ignore
1162       }
1163     }
1164   }
1165
1166   /**
1167    * Helper method to add a codon-to-peptide entry to the aligned codons map
1168    *
1169    * @param alignedCodons
1170    * @param codon
1171    * @param protein
1172    */
1173   protected static void addCodonToMap(
1174           Map<AlignedCodon, Map<SequenceI, AlignedCodon>> alignedCodons,
1175           AlignedCodon codon, SequenceI protein)
1176   {
1177     Map<SequenceI, AlignedCodon> seqProduct = alignedCodons.get(codon);
1178     if (seqProduct == null)
1179     {
1180       seqProduct = new HashMap<SequenceI, AlignedCodon>();
1181       alignedCodons.put(codon, seqProduct);
1182     }
1183     seqProduct.put(protein, codon);
1184   }
1185
1186   /**
1187    * Returns true if a cDNA/Protein mapping either exists, or could be made,
1188    * between at least one pair of sequences in the two alignments. Currently,
1189    * the logic is:
1190    * <ul>
1191    * <li>One alignment must be nucleotide, and the other protein</li>
1192    * <li>At least one pair of sequences must be already mapped, or mappable</li>
1193    * <li>Mappable means the nucleotide translation matches the protein sequence</li>
1194    * <li>The translation may ignore start and stop codons if present in the
1195    * nucleotide</li>
1196    * </ul>
1197    *
1198    * @param al1
1199    * @param al2
1200    * @return
1201    */
1202   public static boolean isMappable(AlignmentI al1, AlignmentI al2)
1203   {
1204     if (al1 == null || al2 == null)
1205     {
1206       return false;
1207     }
1208
1209     /*
1210      * Require one nucleotide and one protein
1211      */
1212     if (al1.isNucleotide() == al2.isNucleotide())
1213     {
1214       return false;
1215     }
1216     AlignmentI dna = al1.isNucleotide() ? al1 : al2;
1217     AlignmentI protein = dna == al1 ? al2 : al1;
1218     List<AlignedCodonFrame> mappings = protein.getCodonFrames();
1219     for (SequenceI dnaSeq : dna.getSequences())
1220     {
1221       for (SequenceI proteinSeq : protein.getSequences())
1222       {
1223         if (isMappable(dnaSeq, proteinSeq, mappings))
1224         {
1225           return true;
1226         }
1227       }
1228     }
1229     return false;
1230   }
1231
1232   /**
1233    * Returns true if the dna sequence is mapped, or could be mapped, to the
1234    * protein sequence.
1235    *
1236    * @param dnaSeq
1237    * @param proteinSeq
1238    * @param mappings
1239    * @return
1240    */
1241   protected static boolean isMappable(SequenceI dnaSeq,
1242           SequenceI proteinSeq, List<AlignedCodonFrame> mappings)
1243   {
1244     if (dnaSeq == null || proteinSeq == null)
1245     {
1246       return false;
1247     }
1248
1249     SequenceI dnaDs = dnaSeq.getDatasetSequence() == null ? dnaSeq : dnaSeq
1250             .getDatasetSequence();
1251     SequenceI proteinDs = proteinSeq.getDatasetSequence() == null ? proteinSeq
1252             : proteinSeq.getDatasetSequence();
1253
1254     for (AlignedCodonFrame mapping : mappings)
1255     {
1256       if (proteinDs == mapping.getAaForDnaSeq(dnaDs))
1257       {
1258         /*
1259          * already mapped
1260          */
1261         return true;
1262       }
1263     }
1264
1265     /*
1266      * Just try to make a mapping (it is not yet stored), test whether
1267      * successful.
1268      */
1269     return mapProteinSequenceToCdna(proteinDs, dnaDs) != null;
1270   }
1271
1272   /**
1273    * Finds any reference annotations associated with the sequences in
1274    * sequenceScope, that are not already added to the alignment, and adds them
1275    * to the 'candidates' map. Also populates a lookup table of annotation
1276    * labels, keyed by calcId, for use in constructing tooltips or the like.
1277    *
1278    * @param sequenceScope
1279    *          the sequences to scan for reference annotations
1280    * @param labelForCalcId
1281    *          (optional) map to populate with label for calcId
1282    * @param candidates
1283    *          map to populate with annotations for sequence
1284    * @param al
1285    *          the alignment to check for presence of annotations
1286    */
1287   public static void findAddableReferenceAnnotations(
1288           List<SequenceI> sequenceScope,
1289           Map<String, String> labelForCalcId,
1290           final Map<SequenceI, List<AlignmentAnnotation>> candidates,
1291           AlignmentI al)
1292   {
1293     if (sequenceScope == null)
1294     {
1295       return;
1296     }
1297
1298     /*
1299      * For each sequence in scope, make a list of any annotations on the
1300      * underlying dataset sequence which are not already on the alignment.
1301      *
1302      * Add to a map of { alignmentSequence, <List of annotations to add> }
1303      */
1304     for (SequenceI seq : sequenceScope)
1305     {
1306       SequenceI dataset = seq.getDatasetSequence();
1307       if (dataset == null)
1308       {
1309         continue;
1310       }
1311       AlignmentAnnotation[] datasetAnnotations = dataset.getAnnotation();
1312       if (datasetAnnotations == null)
1313       {
1314         continue;
1315       }
1316       final List<AlignmentAnnotation> result = new ArrayList<AlignmentAnnotation>();
1317       for (AlignmentAnnotation dsann : datasetAnnotations)
1318       {
1319         /*
1320          * Find matching annotations on the alignment. If none is found, then
1321          * add this annotation to the list of 'addable' annotations for this
1322          * sequence.
1323          */
1324         final Iterable<AlignmentAnnotation> matchedAlignmentAnnotations = al
1325                 .findAnnotations(seq, dsann.getCalcId(), dsann.label);
1326         if (!matchedAlignmentAnnotations.iterator().hasNext())
1327         {
1328           result.add(dsann);
1329           if (labelForCalcId != null)
1330           {
1331             labelForCalcId.put(dsann.getCalcId(), dsann.label);
1332           }
1333         }
1334       }
1335       /*
1336        * Save any addable annotations for this sequence
1337        */
1338       if (!result.isEmpty())
1339       {
1340         candidates.put(seq, result);
1341       }
1342     }
1343   }
1344
1345   /**
1346    * Adds annotations to the top of the alignment annotations, in the same order
1347    * as their related sequences.
1348    *
1349    * @param annotations
1350    *          the annotations to add
1351    * @param alignment
1352    *          the alignment to add them to
1353    * @param selectionGroup
1354    *          current selection group (or null if none)
1355    */
1356   public static void addReferenceAnnotations(
1357           Map<SequenceI, List<AlignmentAnnotation>> annotations,
1358           final AlignmentI alignment, final SequenceGroup selectionGroup)
1359   {
1360     for (SequenceI seq : annotations.keySet())
1361     {
1362       for (AlignmentAnnotation ann : annotations.get(seq))
1363       {
1364         AlignmentAnnotation copyAnn = new AlignmentAnnotation(ann);
1365         int startRes = 0;
1366         int endRes = ann.annotations.length;
1367         if (selectionGroup != null)
1368         {
1369           startRes = selectionGroup.getStartRes();
1370           endRes = selectionGroup.getEndRes();
1371         }
1372         copyAnn.restrict(startRes, endRes);
1373
1374         /*
1375          * Add to the sequence (sets copyAnn.datasetSequence), unless the
1376          * original annotation is already on the sequence.
1377          */
1378         if (!seq.hasAnnotation(ann))
1379         {
1380           seq.addAlignmentAnnotation(copyAnn);
1381         }
1382         // adjust for gaps
1383         copyAnn.adjustForAlignment();
1384         // add to the alignment and set visible
1385         alignment.addAnnotation(copyAnn);
1386         copyAnn.visible = true;
1387       }
1388     }
1389   }
1390
1391   /**
1392    * Set visibility of alignment annotations of specified types (labels), for
1393    * specified sequences. This supports controls like
1394    * "Show all secondary structure", "Hide all Temp factor", etc.
1395    *
1396    * @al the alignment to scan for annotations
1397    * @param types
1398    *          the types (labels) of annotations to be updated
1399    * @param forSequences
1400    *          if not null, only annotations linked to one of these sequences are
1401    *          in scope for update; if null, acts on all sequence annotations
1402    * @param anyType
1403    *          if this flag is true, 'types' is ignored (label not checked)
1404    * @param doShow
1405    *          if true, set visibility on, else set off
1406    */
1407   public static void showOrHideSequenceAnnotations(AlignmentI al,
1408           Collection<String> types, List<SequenceI> forSequences,
1409           boolean anyType, boolean doShow)
1410   {
1411     for (AlignmentAnnotation aa : al.getAlignmentAnnotation())
1412     {
1413       if (anyType || types.contains(aa.label))
1414       {
1415         if ((aa.sequenceRef != null)
1416                 && (forSequences == null || forSequences
1417                         .contains(aa.sequenceRef)))
1418         {
1419           aa.visible = doShow;
1420         }
1421       }
1422     }
1423   }
1424
1425   /**
1426    * Returns true if either sequence has a cross-reference to the other
1427    *
1428    * @param seq1
1429    * @param seq2
1430    * @return
1431    */
1432   public static boolean haveCrossRef(SequenceI seq1, SequenceI seq2)
1433   {
1434     // Note: moved here from class CrossRef as the latter class has dependencies
1435     // not availability to the applet's classpath
1436     return hasCrossRef(seq1, seq2) || hasCrossRef(seq2, seq1);
1437   }
1438
1439   /**
1440    * Returns true if seq1 has a cross-reference to seq2. Currently this assumes
1441    * that sequence name is structured as Source|AccessionId.
1442    *
1443    * @param seq1
1444    * @param seq2
1445    * @return
1446    */
1447   public static boolean hasCrossRef(SequenceI seq1, SequenceI seq2)
1448   {
1449     if (seq1 == null || seq2 == null)
1450     {
1451       return false;
1452     }
1453     String name = seq2.getName();
1454     final DBRefEntry[] xrefs = seq1.getDBRefs();
1455     if (xrefs != null)
1456     {
1457       for (DBRefEntry xref : xrefs)
1458       {
1459         String xrefName = xref.getSource() + "|" + xref.getAccessionId();
1460         // case-insensitive test, consistent with DBRefEntry.equalRef()
1461         if (xrefName.equalsIgnoreCase(name))
1462         {
1463           return true;
1464         }
1465       }
1466     }
1467     return false;
1468   }
1469
1470   /**
1471    * Constructs an alignment consisting of the mapped (CDS) regions in the given
1472    * nucleotide sequences, and updates mappings to match. The new sequences are
1473    * aligned as per the original sequences (with gapped columns omitted).
1474    *
1475    * @param dna
1476    *          aligned dna sequences
1477    * @param mappings
1478    *          from dna to protein; these are replaced with new mappings
1479    * @param gapChar
1480    * @return an alignment whose sequences are the cds-only parts of the dna
1481    *         sequences (or null if no mappings are found)
1482    */
1483   public static AlignmentI makeCdsAlignment(SequenceI[] dna,
1484           List<AlignedCodonFrame> mappings, char gapChar)
1485   {
1486     List<int[]> cdsColumns = findCdsColumns(dna);
1487
1488     /*
1489      * create CDS sequences and new mappings
1490      * (from cdna to cds, and cds to peptide)
1491      */
1492     List<AlignedCodonFrame> newMappings = new ArrayList<AlignedCodonFrame>();
1493     List<SequenceI> cdsSequences = new ArrayList<SequenceI>();
1494
1495     for (SequenceI dnaSeq : dna)
1496     {
1497       final SequenceI ds = dnaSeq.getDatasetSequence();
1498       List<AlignedCodonFrame> seqMappings = MappingUtils
1499               .findMappingsForSequence(ds, mappings);
1500       for (AlignedCodonFrame acf : seqMappings)
1501       {
1502         AlignedCodonFrame newMapping = new AlignedCodonFrame();
1503         final List<SequenceI> mappedCds = makeCdsSequences(dnaSeq, acf,
1504                 cdsColumns, newMapping, gapChar);
1505         if (!mappedCds.isEmpty())
1506         {
1507           cdsSequences.addAll(mappedCds);
1508           newMappings.add(newMapping);
1509         }
1510       }
1511     }
1512     AlignmentI al = new Alignment(
1513             cdsSequences.toArray(new SequenceI[cdsSequences.size()]));
1514     al.setGapCharacter(gapChar);
1515     al.setDataset(null);
1516
1517     /*
1518      * Replace the old mappings with the new ones
1519      */
1520     mappings.clear();
1521     mappings.addAll(newMappings);
1522
1523     return al;
1524   }
1525
1526   /**
1527    * Returns a consolidated list of column ranges where at least one sequence
1528    * has a CDS feature. This assumes CDS features are on genomic sequence i.e.
1529    * are for contiguous CDS ranges (no gaps).
1530    *
1531    * @param seqs
1532    * @return
1533    */
1534   public static List<int[]> findCdsColumns(SequenceI[] seqs)
1535   {
1536     // TODO use refactored code from AlignViewController
1537     // markColumnsContainingFeatures, not reinvent the wheel!
1538
1539     List<int[]> result = new ArrayList<int[]>();
1540     for (SequenceI seq : seqs)
1541     {
1542       result.addAll(findCdsColumns(seq));
1543     }
1544
1545     /*
1546      * sort and compact the list into ascending, non-overlapping ranges
1547      */
1548     Collections.sort(result, new Comparator<int[]>()
1549     {
1550       @Override
1551       public int compare(int[] o1, int[] o2)
1552       {
1553         return Integer.compare(o1[0], o2[0]);
1554       }
1555     });
1556     result = MapList.coalesceRanges(result);
1557
1558     return result;
1559   }
1560
1561   public static List<int[]> findCdsColumns(SequenceI seq)
1562   {
1563     List<int[]> result = new ArrayList<int[]>();
1564     SequenceOntologyI so = SequenceOntologyFactory.getInstance();
1565     SequenceFeature[] sfs = seq.getSequenceFeatures();
1566     if (sfs != null)
1567     {
1568       for (SequenceFeature sf : sfs)
1569       {
1570         if (so.isA(sf.getType(), SequenceOntologyI.CDS))
1571         {
1572           int colStart = seq.findIndex(sf.getBegin());
1573           int colEnd = seq.findIndex(sf.getEnd());
1574           result.add(new int[] { colStart, colEnd });
1575         }
1576       }
1577     }
1578     return result;
1579   }
1580
1581   /**
1582    * Answers true if all sequences have a gap at (or do not extend to) the
1583    * specified column position (base 1)
1584    *
1585    * @param seqs
1586    * @param col
1587    * @return
1588    */
1589   public static boolean isGappedColumn(List<SequenceI> seqs, int col)
1590   {
1591     if (seqs != null)
1592     {
1593       for (SequenceI seq : seqs)
1594       {
1595         if (!Comparison.isGap(seq.getCharAt(col - 1)))
1596         {
1597           return false;
1598         }
1599       }
1600     }
1601     return true;
1602   }
1603
1604   /**
1605    * Returns the column ranges (base 1) of each aligned sequence that are
1606    * involved in any mapping. This is a helper method for aligning protein
1607    * products of aligned transcripts.
1608    *
1609    * @param mappedSequences
1610    *          (possibly gapped) dna sequences
1611    * @param mappings
1612    * @return
1613    */
1614   protected static List<List<int[]>> getMappedColumns(
1615           List<SequenceI> mappedSequences, List<AlignedCodonFrame> mappings)
1616   {
1617     List<List<int[]>> result = new ArrayList<List<int[]>>();
1618     for (SequenceI seq : mappedSequences)
1619     {
1620       List<int[]> columns = new ArrayList<int[]>();
1621       List<AlignedCodonFrame> seqMappings = MappingUtils
1622               .findMappingsForSequence(seq, mappings);
1623       for (AlignedCodonFrame mapping : seqMappings)
1624       {
1625         List<Mapping> maps = mapping.getMappingsForSequence(seq);
1626         for (Mapping map : maps)
1627         {
1628           /*
1629            * Get the codon regions as { [2, 5], [7, 12], [14, 14] etc }
1630            * Find and add the overall aligned column range for each
1631            */
1632           for (int[] cdsRange : map.getMap().getFromRanges())
1633           {
1634             int startPos = cdsRange[0];
1635             int endPos = cdsRange[1];
1636             int startCol = seq.findIndex(startPos);
1637             int endCol = seq.findIndex(endPos);
1638             columns.add(new int[] { startCol, endCol });
1639           }
1640         }
1641       }
1642       result.add(columns);
1643     }
1644     return result;
1645   }
1646
1647   /**
1648    * Helper method to make cds-only sequences and populate their mappings to
1649    * protein products
1650    * <p>
1651    * For example, if ggCCaTTcGAg has mappings [3, 4, 6, 7, 9, 10] to protein
1652    * then generate a sequence CCTTGA with mapping [1, 6] to the same protein
1653    * residues
1654    * <p>
1655    * Typically eukaryotic dna will include cds encoding for a single peptide
1656    * sequence i.e. return a single result. Bacterial dna may have overlapping
1657    * cds mappings coding for multiple peptides so return multiple results
1658    * (example EMBL KF591215).
1659    *
1660    * @param dnaSeq
1661    *          a dna aligned sequence
1662    * @param mapping
1663    *          containing one or more mappings of the sequence to protein
1664    * @param ungappedCdsColumns
1665    * @param newMappings
1666    *          the new mapping to populate, from the cds-only sequences to their
1667    *          mapped protein sequences
1668    * @return
1669    */
1670   protected static List<SequenceI> makeCdsSequences(SequenceI dnaSeq,
1671           AlignedCodonFrame mapping, List<int[]> ungappedCdsColumns,
1672           AlignedCodonFrame newMappings, char gapChar)
1673   {
1674     List<SequenceI> cdsSequences = new ArrayList<SequenceI>();
1675     List<Mapping> seqMappings = mapping.getMappingsForSequence(dnaSeq);
1676
1677     for (Mapping seqMapping : seqMappings)
1678     {
1679       SequenceI cds = makeCdsSequence(dnaSeq, seqMapping,
1680               ungappedCdsColumns, gapChar);
1681       cdsSequences.add(cds);
1682
1683       /*
1684        * add new mappings, from dna to cds, and from cds to peptide
1685        */
1686       MapList dnaToCds = addCdsMappings(dnaSeq.getDatasetSequence(), cds,
1687               seqMapping, newMappings);
1688
1689       /*
1690        * transfer any features on dna that overlap the CDS
1691        */
1692       transferFeatures(dnaSeq, cds, dnaToCds, null, SequenceOntologyI.CDS);
1693     }
1694     return cdsSequences;
1695   }
1696
1697   /**
1698    * Transfers co-located features on 'fromSeq' to 'toSeq', adjusting the
1699    * feature start/end ranges, optionally omitting specified feature types.
1700    * Returns the number of features copied.
1701    *
1702    * @param fromSeq
1703    * @param toSeq
1704    * @param select
1705    *          if not null, only features of this type are copied (including
1706    *          subtypes in the Sequence Ontology)
1707    * @param mapping
1708    *          the mapping from 'fromSeq' to 'toSeq'
1709    * @param omitting
1710    */
1711   public static int transferFeatures(SequenceI fromSeq, SequenceI toSeq,
1712           MapList mapping, String select, String... omitting)
1713   {
1714     SequenceI copyTo = toSeq;
1715     while (copyTo.getDatasetSequence() != null)
1716     {
1717       copyTo = copyTo.getDatasetSequence();
1718     }
1719
1720     SequenceOntologyI so = SequenceOntologyFactory.getInstance();
1721     int count = 0;
1722     SequenceFeature[] sfs = fromSeq.getSequenceFeatures();
1723     if (sfs != null)
1724     {
1725       for (SequenceFeature sf : sfs)
1726       {
1727         String type = sf.getType();
1728         if (select != null && !so.isA(type, select))
1729         {
1730           continue;
1731         }
1732         boolean omit = false;
1733         for (String toOmit : omitting)
1734         {
1735           if (type.equals(toOmit))
1736           {
1737             omit = true;
1738           }
1739         }
1740         if (omit)
1741         {
1742           continue;
1743         }
1744
1745         /*
1746          * locate the mapped range - null if either start or end is
1747          * not mapped (no partial overlaps are calculated)
1748          */
1749         int start = sf.getBegin();
1750         int end = sf.getEnd();
1751         int[] mappedTo = mapping.locateInTo(start, end);
1752         /*
1753          * if whole exon range doesn't map, try interpreting it
1754          * as 5' or 3' exon overlapping the CDS range
1755          */
1756         if (mappedTo == null)
1757         {
1758           mappedTo = mapping.locateInTo(end, end);
1759           if (mappedTo != null)
1760           {
1761             /*
1762              * end of exon is in CDS range - 5' overlap
1763              * to a range from the start of the peptide
1764              */
1765             mappedTo[0] = 1;
1766           }
1767         }
1768         if (mappedTo == null)
1769         {
1770           mappedTo = mapping.locateInTo(start, start);
1771           if (mappedTo != null)
1772           {
1773             /*
1774              * start of exon is in CDS range - 3' overlap
1775              * to a range up to the end of the peptide
1776              */
1777             mappedTo[1] = toSeq.getLength();
1778           }
1779         }
1780         if (mappedTo != null)
1781         {
1782           SequenceFeature copy = new SequenceFeature(sf);
1783           copy.setBegin(Math.min(mappedTo[0], mappedTo[1]));
1784           copy.setEnd(Math.max(mappedTo[0], mappedTo[1]));
1785           copyTo.addSequenceFeature(copy);
1786           count++;
1787         }
1788       }
1789     }
1790     return count;
1791   }
1792
1793   /**
1794    * Creates and adds mappings
1795    * <ul>
1796    * <li>from cds to peptide</li>
1797    * <li>from dna to cds</li>
1798    * </ul>
1799    * and returns the dna-to-cds mapping
1800    *
1801    * @param dnaSeq
1802    * @param cdsSeq
1803    * @param dnaMapping
1804    * @param newMappings
1805    * @return
1806    */
1807   protected static MapList addCdsMappings(SequenceI dnaSeq,
1808           SequenceI cdsSeq, Mapping dnaMapping,
1809           AlignedCodonFrame newMappings)
1810   {
1811     cdsSeq.createDatasetSequence();
1812
1813     /*
1814      * CDS to peptide is just a contiguous 3:1 mapping, with
1815      * the peptide ranges taken unchanged from the dna mapping
1816      */
1817     List<int[]> cdsRanges = new ArrayList<int[]>();
1818     SequenceI cdsDataset = cdsSeq.getDatasetSequence();
1819     cdsRanges.add(new int[] { 1, cdsDataset.getLength() });
1820     MapList cdsToPeptide = new MapList(cdsRanges, dnaMapping.getMap()
1821             .getToRanges(), 3, 1);
1822     newMappings.addMap(cdsDataset, dnaMapping.getTo(), cdsToPeptide);
1823
1824     /*
1825      * dna 'from' ranges map 1:1 to the contiguous extracted CDS
1826      */
1827     MapList dnaToCds = new MapList(dnaMapping.getMap().getFromRanges(),
1828             cdsRanges, 1, 1);
1829     newMappings.addMap(dnaSeq, cdsDataset, dnaToCds);
1830     return dnaToCds;
1831   }
1832
1833   /**
1834    * Makes and returns a CDS-only sequence, where the CDS regions are identified
1835    * as the 'from' ranges of the mapping on the dna.
1836    *
1837    * @param dnaSeq
1838    *          nucleotide sequence
1839    * @param seqMapping
1840    *          mappings from CDS regions of nucleotide
1841    * @param ungappedCdsColumns
1842    * @return
1843    */
1844   protected static SequenceI makeCdsSequence(SequenceI dnaSeq,
1845           Mapping seqMapping, List<int[]> ungappedCdsColumns, char gapChar)
1846   {
1847     int cdsWidth = MappingUtils.getLength(ungappedCdsColumns);
1848
1849     /*
1850      * populate CDS columns with the aligned
1851      * column character if that column is mapped (which may be a gap
1852      * if an intron interrupts a codon), else with a gap
1853      */
1854     List<int[]> fromRanges = seqMapping.getMap().getFromRanges();
1855     char[] cdsChars = new char[cdsWidth];
1856     int pos = 0;
1857     for (int[] columns : ungappedCdsColumns)
1858     {
1859       for (int i = columns[0]; i <= columns[1]; i++)
1860       {
1861         char dnaChar = dnaSeq.getCharAt(i - 1);
1862         if (Comparison.isGap(dnaChar))
1863         {
1864           cdsChars[pos] = gapChar;
1865         }
1866         else
1867         {
1868           int seqPos = dnaSeq.findPosition(i - 1);
1869           if (MappingUtils.contains(fromRanges, seqPos))
1870           {
1871             cdsChars[pos] = dnaChar;
1872           }
1873           else
1874           {
1875             cdsChars[pos] = gapChar;
1876           }
1877         }
1878         pos++;
1879       }
1880     }
1881     SequenceI cdsSequence = new Sequence(dnaSeq.getName(),
1882             String.valueOf(cdsChars));
1883
1884     transferDbRefs(seqMapping.getTo(), cdsSequence);
1885
1886     return cdsSequence;
1887   }
1888
1889   /**
1890    * Locate any xrefs to CDS databases on the protein product and attach to the
1891    * CDS sequence. Also add as a sub-token of the sequence name.
1892    *
1893    * @param from
1894    * @param to
1895    */
1896   protected static void transferDbRefs(SequenceI from, SequenceI to)
1897   {
1898     String cdsAccId = FeatureProperties.getCodingFeature(DBRefSource.EMBL);
1899     DBRefEntry[] cdsRefs = DBRefUtils.selectRefs(from.getDBRefs(),
1900             DBRefSource.CODINGDBS);
1901     if (cdsRefs != null)
1902     {
1903       for (DBRefEntry cdsRef : cdsRefs)
1904       {
1905         to.addDBRef(new DBRefEntry(cdsRef));
1906         cdsAccId = cdsRef.getAccessionId();
1907       }
1908     }
1909     if (!to.getName().contains(cdsAccId))
1910     {
1911       to.setName(to.getName() + "|" + cdsAccId);
1912     }
1913   }
1914 }