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