src/jalview/analysis/AlignmentUtils.java

   1 /*
   2  * Jalview - A Sequence Alignment Editor and Viewer ($$Version-Rel$$)
   3  * Copyright (C) $$Year-Rel$$ The Jalview Authors
   4  *
   5  * This file is part of Jalview.
   6  *
   7  * Jalview is free software: you can redistribute it and/or
   8  * modify it under the terms of the GNU General Public License
   9  * as published by the Free Software Foundation, either version 3
  10  * of the License, or (at your option) any later version.
  11  *
  12  * Jalview is distributed in the hope that it will be useful, but
  13  * WITHOUT ANY WARRANTY; without even the implied warranty
  14  * of MERCHANTABILITY or FITNESS FOR A PARTICULAR
  15  * PURPOSE.  See the GNU General Public License for more details.
  16  *
  17  * You should have received a copy of the GNU General Public License
  18  * along with Jalview.  If not, see <http://www.gnu.org/licenses/>.
  19  * The Jalview Authors are detailed in the 'AUTHORS' file.
  20  */
  21 package jalview.analysis;
  22
  23 import static jalview.io.gff.GffConstants.CLINICAL_SIGNIFICANCE;
  24
  25 import jalview.datamodel.AlignedCodon;
  26 import jalview.datamodel.AlignedCodonFrame;
  27 import jalview.datamodel.Alignment;
  28 import jalview.datamodel.AlignmentAnnotation;
  29 import jalview.datamodel.AlignmentI;
  30 import jalview.datamodel.DBRefEntry;
  31 import jalview.datamodel.IncompleteCodonException;
  32 import jalview.datamodel.Mapping;
  33 import jalview.datamodel.Sequence;
  34 import jalview.datamodel.SequenceFeature;
  35 import jalview.datamodel.SequenceGroup;
  36 import jalview.datamodel.SequenceI;
  37 import jalview.io.gff.SequenceOntologyFactory;
  38 import jalview.io.gff.SequenceOntologyI;
  39 import jalview.schemes.ResidueProperties;
  40 import jalview.util.Comparison;
  41 import jalview.util.MapList;
  42 import jalview.util.MappingUtils;
  43
  44 import java.io.UnsupportedEncodingException;
  45 import java.net.URLEncoder;
  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   private static final String SEQUENCE_VARIANT = "sequence_variant:";
  73   private static final String ID = "ID";
  74
  75   /**
  76    * A data model to hold the 'normal' base value at a position, and an optional
  77    * sequence variant feature
  78    */
  79   static class DnaVariant
  80   {
  81     String base;
  82
  83     SequenceFeature variant;
  84
  85     DnaVariant(String nuc)
  86     {
  87       base = nuc;
  88     }
  89
  90     DnaVariant(String nuc, SequenceFeature var)
  91     {
  92       base = nuc;
  93       variant = var;
  94     }
  95   }
  96
  97   /**
  98    * given an existing alignment, create a new alignment including all, or up to
  99    * flankSize additional symbols from each sequence's dataset sequence
 100    *
 101    * @param core
 102    * @param flankSize
 103    * @return AlignmentI
 104    */
 105   public static AlignmentI expandContext(AlignmentI core, int flankSize)
 106   {
 107     List<SequenceI> sq = new ArrayList<SequenceI>();
 108     int maxoffset = 0;
 109     for (SequenceI s : core.getSequences())
 110     {
 111       SequenceI newSeq = s.deriveSequence();
 112       final int newSeqStart = newSeq.getStart() - 1;
 113       if (newSeqStart > maxoffset
 114               && newSeq.getDatasetSequence().getStart() < s.getStart())
 115       {
 116         maxoffset = newSeqStart;
 117       }
 118       sq.add(newSeq);
 119     }
 120     if (flankSize > -1)
 121     {
 122       maxoffset = Math.min(maxoffset, flankSize);
 123     }
 124
 125     /*
 126      * now add offset left and right to create an expanded alignment
 127      */
 128     for (SequenceI s : sq)
 129     {
 130       SequenceI ds = s;
 131       while (ds.getDatasetSequence() != null)
 132       {
 133         ds = ds.getDatasetSequence();
 134       }
 135       int s_end = s.findPosition(s.getStart() + s.getLength());
 136       // find available flanking residues for sequence
 137       int ustream_ds = s.getStart() - ds.getStart();
 138       int dstream_ds = ds.getEnd() - s_end;
 139
 140       // build new flanked sequence
 141
 142       // compute gap padding to start of flanking sequence
 143       int offset = maxoffset - ustream_ds;
 144
 145       // padding is gapChar x ( maxoffset - min(ustream_ds, flank)
 146       if (flankSize >= 0)
 147       {
 148         if (flankSize < ustream_ds)
 149         {
 150           // take up to flankSize residues
 151           offset = maxoffset - flankSize;
 152           ustream_ds = flankSize;
 153         }
 154         if (flankSize <= dstream_ds)
 155         {
 156           dstream_ds = flankSize - 1;
 157         }
 158       }
 159       // TODO use Character.toLowerCase to avoid creating String objects?
 160       char[] upstream = new String(ds.getSequence(s.getStart() - 1
 161               - ustream_ds, s.getStart() - 1)).toLowerCase().toCharArray();
 162       char[] downstream = new String(ds.getSequence(s_end - 1, s_end
 163               + dstream_ds)).toLowerCase().toCharArray();
 164       char[] coreseq = s.getSequence();
 165       char[] nseq = new char[offset + upstream.length + downstream.length
 166               + coreseq.length];
 167       char c = core.getGapCharacter();
 168
 169       int p = 0;
 170       for (; p < offset; p++)
 171       {
 172         nseq[p] = c;
 173       }
 174
 175       System.arraycopy(upstream, 0, nseq, p, upstream.length);
 176       System.arraycopy(coreseq, 0, nseq, p + upstream.length,
 177               coreseq.length);
 178       System.arraycopy(downstream, 0, nseq, p + coreseq.length
 179               + upstream.length, downstream.length);
 180       s.setSequence(new String(nseq));
 181       s.setStart(s.getStart() - ustream_ds);
 182       s.setEnd(s_end + downstream.length);
 183     }
 184     AlignmentI newAl = new jalview.datamodel.Alignment(
 185             sq.toArray(new SequenceI[0]));
 186     for (SequenceI s : sq)
 187     {
 188       if (s.getAnnotation() != null)
 189       {
 190         for (AlignmentAnnotation aa : s.getAnnotation())
 191         {
 192           aa.adjustForAlignment(); // JAL-1712 fix
 193           newAl.addAnnotation(aa);
 194         }
 195       }
 196     }
 197     newAl.setDataset(core.getDataset());
 198     return newAl;
 199   }
 200
 201   /**
 202    * Returns the index (zero-based position) of a sequence in an alignment, or
 203    * -1 if not found.
 204    *
 205    * @param al
 206    * @param seq
 207    * @return
 208    */
 209   public static int getSequenceIndex(AlignmentI al, SequenceI seq)
 210   {
 211     int result = -1;
 212     int pos = 0;
 213     for (SequenceI alSeq : al.getSequences())
 214     {
 215       if (alSeq == seq)
 216       {
 217         result = pos;
 218         break;
 219       }
 220       pos++;
 221     }
 222     return result;
 223   }
 224
 225   /**
 226    * Returns a map of lists of sequences in the alignment, keyed by sequence
 227    * name. For use in mapping between different alignment views of the same
 228    * sequences.
 229    *
 230    * @see jalview.datamodel.AlignmentI#getSequencesByName()
 231    */
 232   public static Map<String, List<SequenceI>> getSequencesByName(
 233           AlignmentI al)
 234   {
 235     Map<String, List<SequenceI>> theMap = new LinkedHashMap<String, List<SequenceI>>();
 236     for (SequenceI seq : al.getSequences())
 237     {
 238       String name = seq.getName();
 239       if (name != null)
 240       {
 241         List<SequenceI> seqs = theMap.get(name);
 242         if (seqs == null)
 243         {
 244           seqs = new ArrayList<SequenceI>();
 245           theMap.put(name, seqs);
 246         }
 247         seqs.add(seq);
 248       }
 249     }
 250     return theMap;
 251   }
 252
 253   /**
 254    * Build mapping of protein to cDNA alignment. Mappings are made between
 255    * sequences where the cDNA translates to the protein sequence. Any new
 256    * mappings are added to the protein alignment. Returns true if any mappings
 257    * either already exist or were added, else false.
 258    *
 259    * @param proteinAlignment
 260    * @param cdnaAlignment
 261    * @return
 262    */
 263   public static boolean mapProteinAlignmentToCdna(
 264           final AlignmentI proteinAlignment, final AlignmentI cdnaAlignment)
 265   {
 266     if (proteinAlignment == null || cdnaAlignment == null)
 267     {
 268       return false;
 269     }
 270
 271     Set<SequenceI> mappedDna = new HashSet<SequenceI>();
 272     Set<SequenceI> mappedProtein = new HashSet<SequenceI>();
 273
 274     /*
 275      * First pass - map sequences where cross-references exist. This include
 276      * 1-to-many mappings to support, for example, variant cDNA.
 277      */
 278     boolean mappingPerformed = mapProteinToCdna(proteinAlignment,
 279             cdnaAlignment, mappedDna, mappedProtein, true);
 280
 281     /*
 282      * Second pass - map sequences where no cross-references exist. This only
 283      * does 1-to-1 mappings and assumes corresponding sequences are in the same
 284      * order in the alignments.
 285      */
 286     mappingPerformed |= mapProteinToCdna(proteinAlignment, cdnaAlignment,
 287             mappedDna, mappedProtein, false);
 288     return mappingPerformed;
 289   }
 290
 291   /**
 292    * Make mappings between compatible sequences (where the cDNA translation
 293    * matches the protein).
 294    *
 295    * @param proteinAlignment
 296    * @param cdnaAlignment
 297    * @param mappedDna
 298    *          a set of mapped DNA sequences (to add to)
 299    * @param mappedProtein
 300    *          a set of mapped Protein sequences (to add to)
 301    * @param xrefsOnly
 302    *          if true, only map sequences where xrefs exist
 303    * @return
 304    */
 305   protected static boolean mapProteinToCdna(
 306           final AlignmentI proteinAlignment,
 307           final AlignmentI cdnaAlignment, Set<SequenceI> mappedDna,
 308           Set<SequenceI> mappedProtein, boolean xrefsOnly)
 309   {
 310     boolean mappingExistsOrAdded = false;
 311     List<SequenceI> thisSeqs = proteinAlignment.getSequences();
 312     for (SequenceI aaSeq : thisSeqs)
 313     {
 314       boolean proteinMapped = false;
 315       AlignedCodonFrame acf = new AlignedCodonFrame();
 316
 317       for (SequenceI cdnaSeq : cdnaAlignment.getSequences())
 318       {
 319         /*
 320          * Always try to map if sequences have xref to each other; this supports
 321          * variant cDNA or alternative splicing for a protein sequence.
 322          *
 323          * If no xrefs, try to map progressively, assuming that alignments have
 324          * mappable sequences in corresponding order. These are not
 325          * many-to-many, as that would risk mixing species with similar cDNA
 326          * sequences.
 327          */
 328         if (xrefsOnly && !AlignmentUtils.haveCrossRef(aaSeq, cdnaSeq))
 329         {
 330           continue;
 331         }
 332
 333         /*
 334          * Don't map non-xrefd sequences more than once each. This heuristic
 335          * allows us to pair up similar sequences in ordered alignments.
 336          */
 337         if (!xrefsOnly
 338                 && (mappedProtein.contains(aaSeq) || mappedDna
 339                         .contains(cdnaSeq)))
 340         {
 341           continue;
 342         }
 343         if (mappingExists(proteinAlignment.getCodonFrames(),
 344                 aaSeq.getDatasetSequence(), cdnaSeq.getDatasetSequence()))
 345         {
 346           mappingExistsOrAdded = true;
 347         }
 348         else
 349         {
 350           MapList map = mapCdnaToProtein(aaSeq, cdnaSeq);
 351           if (map != null)
 352           {
 353             acf.addMap(cdnaSeq, aaSeq, map);
 354             mappingExistsOrAdded = true;
 355             proteinMapped = true;
 356             mappedDna.add(cdnaSeq);
 357             mappedProtein.add(aaSeq);
 358           }
 359         }
 360       }
 361       if (proteinMapped)
 362       {
 363         proteinAlignment.addCodonFrame(acf);
 364       }
 365     }
 366     return mappingExistsOrAdded;
 367   }
 368
 369   /**
 370    * Answers true if the mappings include one between the given (dataset)
 371    * sequences.
 372    */
 373   public static boolean mappingExists(List<AlignedCodonFrame> mappings,
 374           SequenceI aaSeq, SequenceI cdnaSeq)
 375   {
 376     if (mappings != null)
 377     {
 378       for (AlignedCodonFrame acf : mappings)
 379       {
 380         if (cdnaSeq == acf.getDnaForAaSeq(aaSeq))
 381         {
 382           return true;
 383         }
 384       }
 385     }
 386     return false;
 387   }
 388
 389   /**
 390    * Builds a mapping (if possible) of a cDNA to a protein sequence.
 391    * <ul>
 392    * <li>first checks if the cdna translates exactly to the protein sequence</li>
 393    * <li>else checks for translation after removing a STOP codon</li>
 394    * <li>else checks for translation after removing a START codon</li>
 395    * <li>if that fails, inspect CDS features on the cDNA sequence</li>
 396    * </ul>
 397    * Returns null if no mapping is determined.
 398    *
 399    * @param proteinSeq
 400    *          the aligned protein sequence
 401    * @param cdnaSeq
 402    *          the aligned cdna sequence
 403    * @return
 404    */
 405   public static MapList mapCdnaToProtein(SequenceI proteinSeq,
 406           SequenceI cdnaSeq)
 407   {
 408     /*
 409      * Here we handle either dataset sequence set (desktop) or absent (applet).
 410      * Use only the char[] form of the sequence to avoid creating possibly large
 411      * String objects.
 412      */
 413     final SequenceI proteinDataset = proteinSeq.getDatasetSequence();
 414     char[] aaSeqChars = proteinDataset != null ? proteinDataset
 415             .getSequence() : proteinSeq.getSequence();
 416     final SequenceI cdnaDataset = cdnaSeq.getDatasetSequence();
 417     char[] cdnaSeqChars = cdnaDataset != null ? cdnaDataset.getSequence()
 418             : cdnaSeq.getSequence();
 419     if (aaSeqChars == null || cdnaSeqChars == null)
 420     {
 421       return null;
 422     }
 423
 424     /*
 425      * cdnaStart/End, proteinStartEnd are base 1 (for dataset sequence mapping)
 426      */
 427     final int mappedLength = 3 * aaSeqChars.length;
 428     int cdnaLength = cdnaSeqChars.length;
 429     int cdnaStart = cdnaSeq.getStart();
 430     int cdnaEnd = cdnaSeq.getEnd();
 431     final int proteinStart = proteinSeq.getStart();
 432     final int proteinEnd = proteinSeq.getEnd();
 433
 434     /*
 435      * If lengths don't match, try ignoring stop codon (if present)
 436      */
 437     if (cdnaLength != mappedLength && cdnaLength > 2)
 438     {
 439       String lastCodon = String.valueOf(cdnaSeqChars, cdnaLength - 3, 3)
 440               .toUpperCase();
 441       for (String stop : ResidueProperties.STOP)
 442       {
 443         if (lastCodon.equals(stop))
 444         {
 445           cdnaEnd -= 3;
 446           cdnaLength -= 3;
 447           break;
 448         }
 449       }
 450     }
 451
 452     /*
 453      * If lengths still don't match, try ignoring start codon.
 454      */
 455     int startOffset = 0;
 456     if (cdnaLength != mappedLength
 457             && cdnaLength > 2
 458             && String.valueOf(cdnaSeqChars, 0, 3).toUpperCase()
 459                     .equals(ResidueProperties.START))
 460     {
 461       startOffset += 3;
 462       cdnaStart += 3;
 463       cdnaLength -= 3;
 464     }
 465
 466     if (translatesAs(cdnaSeqChars, startOffset, aaSeqChars))
 467     {
 468       /*
 469        * protein is translation of dna (+/- start/stop codons)
 470        */
 471       MapList map = new MapList(new int[] { cdnaStart, cdnaEnd }, new int[]
 472       { proteinStart, proteinEnd }, 3, 1);
 473       return map;
 474     }
 475
 476     /*
 477      * translation failed - try mapping CDS annotated regions of dna
 478      */
 479     return mapCdsToProtein(cdnaSeq, proteinSeq);
 480   }
 481
 482   /**
 483    * Test whether the given cdna sequence, starting at the given offset,
 484    * translates to the given amino acid sequence, using the standard translation
 485    * table. Designed to fail fast i.e. as soon as a mismatch position is found.
 486    *
 487    * @param cdnaSeqChars
 488    * @param cdnaStart
 489    * @param aaSeqChars
 490    * @return
 491    */
 492   protected static boolean translatesAs(char[] cdnaSeqChars, int cdnaStart,
 493           char[] aaSeqChars)
 494   {
 495     if (cdnaSeqChars == null || aaSeqChars == null)
 496     {
 497       return false;
 498     }
 499
 500     int aaPos = 0;
 501     int dnaPos = cdnaStart;
 502     for (; dnaPos < cdnaSeqChars.length - 2
 503             && aaPos < aaSeqChars.length; dnaPos += 3, aaPos++)
 504     {
 505       String codon = String.valueOf(cdnaSeqChars, dnaPos, 3);
 506       final String translated = ResidueProperties.codonTranslate(codon);
 507
 508       /*
 509        * allow * in protein to match untranslatable in dna
 510        */
 511       final char aaRes = aaSeqChars[aaPos];
 512       if ((translated == null || "STOP".equals(translated)) && aaRes == '*')
 513       {
 514         continue;
 515       }
 516       if (translated == null || !(aaRes == translated.charAt(0)))
 517       {
 518         // debug
 519         // System.out.println(("Mismatch at " + i + "/" + aaResidue + ": "
 520         // + codon + "(" + translated + ") != " + aaRes));
 521         return false;
 522       }
 523     }
 524
 525     /*
 526      * check we matched all of the protein sequence
 527      */
 528     if (aaPos != aaSeqChars.length)
 529     {
 530       return false;
 531     }
 532
 533     /*
 534      * check we matched all of the dna except
 535      * for optional trailing STOP codon
 536      */
 537     if (dnaPos == cdnaSeqChars.length)
 538     {
 539       return true;
 540     }
 541     if (dnaPos == cdnaSeqChars.length - 3)
 542     {
 543       String codon = String.valueOf(cdnaSeqChars, dnaPos, 3);
 544       if ("STOP".equals(ResidueProperties.codonTranslate(codon)))
 545       {
 546         return true;
 547       }
 548     }
 549     return false;
 550   }
 551
 552   /**
 553    * Align sequence 'seq' to match the alignment of a mapped sequence. Note this
 554    * currently assumes that we are aligning cDNA to match protein.
 555    *
 556    * @param seq
 557    *          the sequence to be realigned
 558    * @param al
 559    *          the alignment whose sequence alignment is to be 'copied'
 560    * @param gap
 561    *          character string represent a gap in the realigned sequence
 562    * @param preserveUnmappedGaps
 563    * @param preserveMappedGaps
 564    * @return true if the sequence was realigned, false if it could not be
 565    */
 566   public static boolean alignSequenceAs(SequenceI seq, AlignmentI al,
 567           String gap, boolean preserveMappedGaps,
 568           boolean preserveUnmappedGaps)
 569   {
 570     /*
 571      * Get any mappings from the source alignment to the target (dataset)
 572      * sequence.
 573      */
 574     // TODO there may be one AlignedCodonFrame per dataset sequence, or one with
 575     // all mappings. Would it help to constrain this?
 576     List<AlignedCodonFrame> mappings = al.getCodonFrame(seq);
 577     if (mappings == null || mappings.isEmpty())
 578     {
 579       return false;
 580     }
 581
 582     /*
 583      * Locate the aligned source sequence whose dataset sequence is mapped. We
 584      * just take the first match here (as we can't align like more than one
 585      * sequence).
 586      */
 587     SequenceI alignFrom = null;
 588     AlignedCodonFrame mapping = null;
 589     for (AlignedCodonFrame mp : mappings)
 590     {
 591       alignFrom = mp.findAlignedSequence(seq, al);
 592       if (alignFrom != null)
 593       {
 594         mapping = mp;
 595         break;
 596       }
 597     }
 598
 599     if (alignFrom == null)
 600     {
 601       return false;
 602     }
 603     alignSequenceAs(seq, alignFrom, mapping, gap, al.getGapCharacter(),
 604             preserveMappedGaps, preserveUnmappedGaps);
 605     return true;
 606   }
 607
 608   /**
 609    * Align sequence 'alignTo' the same way as 'alignFrom', using the mapping to
 610    * match residues and codons. Flags control whether existing gaps in unmapped
 611    * (intron) and mapped (exon) regions are preserved or not. Gaps between
 612    * intron and exon are only retained if both flags are set.
 613    *
 614    * @param alignTo
 615    * @param alignFrom
 616    * @param mapping
 617    * @param myGap
 618    * @param sourceGap
 619    * @param preserveUnmappedGaps
 620    * @param preserveMappedGaps
 621    */
 622   public static void alignSequenceAs(SequenceI alignTo,
 623           SequenceI alignFrom, AlignedCodonFrame mapping, String myGap,
 624           char sourceGap, boolean preserveMappedGaps,
 625           boolean preserveUnmappedGaps)
 626   {
 627     // TODO generalise to work for Protein-Protein, dna-dna, dna-protein
 628
 629     // aligned and dataset sequence positions, all base zero
 630     int thisSeqPos = 0;
 631     int sourceDsPos = 0;
 632
 633     int basesWritten = 0;
 634     char myGapChar = myGap.charAt(0);
 635     int ratio = myGap.length();
 636
 637     int fromOffset = alignFrom.getStart() - 1;
 638     int toOffset = alignTo.getStart() - 1;
 639     int sourceGapMappedLength = 0;
 640     boolean inExon = false;
 641     final char[] thisSeq = alignTo.getSequence();
 642     final char[] thatAligned = alignFrom.getSequence();
 643     StringBuilder thisAligned = new StringBuilder(2 * thisSeq.length);
 644
 645     /*
 646      * Traverse the 'model' aligned sequence
 647      */
 648     for (char sourceChar : thatAligned)
 649     {
 650       if (sourceChar == sourceGap)
 651       {
 652         sourceGapMappedLength += ratio;
 653         continue;
 654       }
 655
 656       /*
 657        * Found a non-gap character. Locate its mapped region if any.
 658        */
 659       sourceDsPos++;
 660       // Note mapping positions are base 1, our sequence positions base 0
 661       int[] mappedPos = mapping.getMappedRegion(alignTo, alignFrom,
 662               sourceDsPos + fromOffset);
 663       if (mappedPos == null)
 664       {
 665         /*
 666          * unmapped position; treat like a gap
 667          */
 668         sourceGapMappedLength += ratio;
 669         // System.err.println("Can't align: no codon mapping to residue "
 670         // + sourceDsPos + "(" + sourceChar + ")");
 671         // return;
 672         continue;
 673       }
 674
 675       int mappedCodonStart = mappedPos[0]; // position (1...) of codon start
 676       int mappedCodonEnd = mappedPos[mappedPos.length - 1]; // codon end pos
 677       StringBuilder trailingCopiedGap = new StringBuilder();
 678
 679       /*
 680        * Copy dna sequence up to and including this codon. Optionally, include
 681        * gaps before the codon starts (in introns) and/or after the codon starts
 682        * (in exons).
 683        *
 684        * Note this only works for 'linear' splicing, not reverse or interleaved.
 685        * But then 'align dna as protein' doesn't make much sense otherwise.
 686        */
 687       int intronLength = 0;
 688       while (basesWritten + toOffset < mappedCodonEnd
 689               && thisSeqPos < thisSeq.length)
 690       {
 691         final char c = thisSeq[thisSeqPos++];
 692         if (c != myGapChar)
 693         {
 694           basesWritten++;
 695           int sourcePosition = basesWritten + toOffset;
 696           if (sourcePosition < mappedCodonStart)
 697           {
 698             /*
 699              * Found an unmapped (intron) base. First add in any preceding gaps
 700              * (if wanted).
 701              */
 702             if (preserveUnmappedGaps && trailingCopiedGap.length() > 0)
 703             {
 704               thisAligned.append(trailingCopiedGap.toString());
 705               intronLength += trailingCopiedGap.length();
 706               trailingCopiedGap = new StringBuilder();
 707             }
 708             intronLength++;
 709             inExon = false;
 710           }
 711           else
 712           {
 713             final boolean startOfCodon = sourcePosition == mappedCodonStart;
 714             int gapsToAdd = calculateGapsToInsert(preserveMappedGaps,
 715                     preserveUnmappedGaps, sourceGapMappedLength, inExon,
 716                     trailingCopiedGap.length(), intronLength, startOfCodon);
 717             for (int i = 0; i < gapsToAdd; i++)
 718             {
 719               thisAligned.append(myGapChar);
 720             }
 721             sourceGapMappedLength = 0;
 722             inExon = true;
 723           }
 724           thisAligned.append(c);
 725           trailingCopiedGap = new StringBuilder();
 726         }
 727         else
 728         {
 729           if (inExon && preserveMappedGaps)
 730           {
 731             trailingCopiedGap.append(myGapChar);
 732           }
 733           else if (!inExon && preserveUnmappedGaps)
 734           {
 735             trailingCopiedGap.append(myGapChar);
 736           }
 737         }
 738       }
 739     }
 740
 741     /*
 742      * At end of model aligned sequence. Copy any remaining target sequence, optionally
 743      * including (intron) gaps.
 744      */
 745     while (thisSeqPos < thisSeq.length)
 746     {
 747       final char c = thisSeq[thisSeqPos++];
 748       if (c != myGapChar || preserveUnmappedGaps)
 749       {
 750         thisAligned.append(c);
 751       }
 752       sourceGapMappedLength--;
 753     }
 754
 755     /*
 756      * finally add gaps to pad for any trailing source gaps or
 757      * unmapped characters
 758      */
 759     if (preserveUnmappedGaps)
 760     {
 761       while (sourceGapMappedLength > 0)
 762       {
 763         thisAligned.append(myGapChar);
 764         sourceGapMappedLength--;
 765       }
 766     }
 767
 768     /*
 769      * All done aligning, set the aligned sequence.
 770      */
 771     alignTo.setSequence(new String(thisAligned));
 772   }
 773
 774   /**
 775    * Helper method to work out how many gaps to insert when realigning.
 776    *
 777    * @param preserveMappedGaps
 778    * @param preserveUnmappedGaps
 779    * @param sourceGapMappedLength
 780    * @param inExon
 781    * @param trailingCopiedGap
 782    * @param intronLength
 783    * @param startOfCodon
 784    * @return
 785    */
 786   protected static int calculateGapsToInsert(boolean preserveMappedGaps,
 787           boolean preserveUnmappedGaps, int sourceGapMappedLength,
 788           boolean inExon, int trailingGapLength, int intronLength,
 789           final boolean startOfCodon)
 790   {
 791     int gapsToAdd = 0;
 792     if (startOfCodon)
 793     {
 794       /*
 795        * Reached start of codon. Ignore trailing gaps in intron unless we are
 796        * preserving gaps in both exon and intron. Ignore them anyway if the
 797        * protein alignment introduces a gap at least as large as the intronic
 798        * region.
 799        */
 800       if (inExon && !preserveMappedGaps)
 801       {
 802         trailingGapLength = 0;
 803       }
 804       if (!inExon && !(preserveMappedGaps && preserveUnmappedGaps))
 805       {
 806         trailingGapLength = 0;
 807       }
 808       if (inExon)
 809       {
 810         gapsToAdd = Math.max(sourceGapMappedLength, trailingGapLength);
 811       }
 812       else
 813       {
 814         if (intronLength + trailingGapLength <= sourceGapMappedLength)
 815         {
 816           gapsToAdd = sourceGapMappedLength - intronLength;
 817         }
 818         else
 819         {
 820           gapsToAdd = Math.min(intronLength + trailingGapLength
 821                   - sourceGapMappedLength, trailingGapLength);
 822         }
 823       }
 824     }
 825     else
 826     {
 827       /*
 828        * second or third base of codon; check for any gaps in dna
 829        */
 830       if (!preserveMappedGaps)
 831       {
 832         trailingGapLength = 0;
 833       }
 834       gapsToAdd = Math.max(sourceGapMappedLength, trailingGapLength);
 835     }
 836     return gapsToAdd;
 837   }
 838
 839   /**
 840    * Realigns the given protein to match the alignment of the dna, using codon
 841    * mappings to translate aligned codon positions to protein residues.
 842    *
 843    * @param protein
 844    *          the alignment whose sequences are realigned by this method
 845    * @param dna
 846    *          the dna alignment whose alignment we are 'copying'
 847    * @return the number of sequences that were realigned
 848    */
 849   public static int alignProteinAsDna(AlignmentI protein, AlignmentI dna)
 850   {
 851     List<SequenceI> unmappedProtein = new ArrayList<SequenceI>();
 852     Map<AlignedCodon, Map<SequenceI, AlignedCodon>> alignedCodons = buildCodonColumnsMap(
 853             protein, dna, unmappedProtein);
 854     return alignProteinAs(protein, alignedCodons, unmappedProtein);
 855   }
 856
 857   /**
 858    * Builds a map whose key is an aligned codon position (3 alignment column
 859    * numbers base 0), and whose value is a map from protein sequence to each
 860    * protein's peptide residue for that codon. The map generates an ordering of
 861    * the codons, and allows us to read off the peptides at each position in
 862    * order to assemble 'aligned' protein sequences.
 863    *
 864    * @param protein
 865    *          the protein alignment
 866    * @param dna
 867    *          the coding dna alignment
 868    * @param unmappedProtein
 869    *          any unmapped proteins are added to this list
 870    * @return
 871    */
 872   protected static Map<AlignedCodon, Map<SequenceI, AlignedCodon>> buildCodonColumnsMap(
 873           AlignmentI protein, AlignmentI dna,
 874           List<SequenceI> unmappedProtein)
 875   {
 876     /*
 877      * maintain a list of any proteins with no mappings - these will be
 878      * rendered 'as is' in the protein alignment as we can't align them
 879      */
 880     unmappedProtein.addAll(protein.getSequences());
 881
 882     List<AlignedCodonFrame> mappings = protein.getCodonFrames();
 883
 884     /*
 885      * Map will hold, for each aligned codon position e.g. [3, 5, 6], a map of
 886      * {dnaSequence, {proteinSequence, codonProduct}} at that position. The
 887      * comparator keeps the codon positions ordered.
 888      */
 889     Map<AlignedCodon, Map<SequenceI, AlignedCodon>> alignedCodons = new TreeMap<AlignedCodon, Map<SequenceI, AlignedCodon>>(
 890             new CodonComparator());
 891
 892     for (SequenceI dnaSeq : dna.getSequences())
 893     {
 894       for (AlignedCodonFrame mapping : mappings)
 895       {
 896         SequenceI prot = mapping.findAlignedSequence(dnaSeq, protein);
 897         if (prot != null)
 898         {
 899           Mapping seqMap = mapping.getMappingForSequence(dnaSeq);
 900           addCodonPositions(dnaSeq, prot, protein.getGapCharacter(),
 901                   seqMap, alignedCodons);
 902           unmappedProtein.remove(prot);
 903         }
 904       }
 905     }
 906
 907     /*
 908      * Finally add any unmapped peptide start residues (e.g. for incomplete
 909      * codons) as if at the codon position before the second residue
 910      */
 911     // TODO resolve JAL-2022 so this fudge can be removed
 912     int mappedSequenceCount = protein.getHeight() - unmappedProtein.size();
 913     addUnmappedPeptideStarts(alignedCodons, mappedSequenceCount);
 914
 915     return alignedCodons;
 916   }
 917
 918   /**
 919    * Scans for any protein mapped from position 2 (meaning unmapped start
 920    * position e.g. an incomplete codon), and synthesizes a 'codon' for it at the
 921    * preceding position in the alignment
 922    *
 923    * @param alignedCodons
 924    *          the codon-to-peptide map
 925    * @param mappedSequenceCount
 926    *          the number of distinct sequences in the map
 927    */
 928   protected static void addUnmappedPeptideStarts(
 929           Map<AlignedCodon, Map<SequenceI, AlignedCodon>> alignedCodons,
 930           int mappedSequenceCount)
 931   {
 932     // TODO delete this ugly hack once JAL-2022 is resolved
 933     // i.e. we can model startPhase > 0 (incomplete start codon)
 934
 935     List<SequenceI> sequencesChecked = new ArrayList<SequenceI>();
 936     AlignedCodon lastCodon = null;
 937     Map<SequenceI, AlignedCodon> toAdd = new HashMap<SequenceI, AlignedCodon>();
 938
 939     for (Entry<AlignedCodon, Map<SequenceI, AlignedCodon>> entry : alignedCodons
 940             .entrySet())
 941     {
 942       for (Entry<SequenceI, AlignedCodon> sequenceCodon : entry.getValue()
 943               .entrySet())
 944       {
 945         SequenceI seq = sequenceCodon.getKey();
 946         if (sequencesChecked.contains(seq))
 947         {
 948           continue;
 949         }
 950         sequencesChecked.add(seq);
 951         AlignedCodon codon = sequenceCodon.getValue();
 952         if (codon.peptideCol > 1)
 953         {
 954           System.err
 955                   .println("Problem mapping protein with >1 unmapped start positions: "
 956                           + seq.getName());
 957         }
 958         else if (codon.peptideCol == 1)
 959         {
 960           /*
 961            * first position (peptideCol == 0) was unmapped - add it
 962            */
 963           if (lastCodon != null)
 964           {
 965             AlignedCodon firstPeptide = new AlignedCodon(lastCodon.pos1,
 966                     lastCodon.pos2, lastCodon.pos3, String.valueOf(seq
 967                             .getCharAt(0)), 0);
 968             toAdd.put(seq, firstPeptide);
 969           }
 970           else
 971           {
 972             /*
 973              * unmapped residue at start of alignment (no prior column) -
 974              * 'insert' at nominal codon [0, 0, 0]
 975              */
 976             AlignedCodon firstPeptide = new AlignedCodon(0, 0, 0,
 977                     String.valueOf(seq.getCharAt(0)), 0);
 978             toAdd.put(seq, firstPeptide);
 979           }
 980         }
 981         if (sequencesChecked.size() == mappedSequenceCount)
 982         {
 983           // no need to check past first mapped position in all sequences
 984           break;
 985         }
 986       }
 987       lastCodon = entry.getKey();
 988     }
 989
 990     /*
 991      * add any new codons safely after iterating over the map
 992      */
 993     for (Entry<SequenceI, AlignedCodon> startCodon : toAdd.entrySet())
 994     {
 995       addCodonToMap(alignedCodons, startCodon.getValue(),
 996               startCodon.getKey());
 997     }
 998   }
 999
1000   /**
1001    * Update the aligned protein sequences to match the codon alignments given in
1002    * the map.
1003    *
1004    * @param protein
1005    * @param alignedCodons
1006    *          an ordered map of codon positions (columns), with sequence/peptide
1007    *          values present in each column
1008    * @param unmappedProtein
1009    * @return
1010    */
1011   protected static int alignProteinAs(AlignmentI protein,
1012           Map<AlignedCodon, Map<SequenceI, AlignedCodon>> alignedCodons,
1013           List<SequenceI> unmappedProtein)
1014   {
1015     /*
1016      * Prefill aligned sequences with gaps before inserting aligned protein
1017      * residues.
1018      */
1019     int alignedWidth = alignedCodons.size();
1020     char[] gaps = new char[alignedWidth];
1021     Arrays.fill(gaps, protein.getGapCharacter());
1022     String allGaps = String.valueOf(gaps);
1023     for (SequenceI seq : protein.getSequences())
1024     {
1025       if (!unmappedProtein.contains(seq))
1026       {
1027         seq.setSequence(allGaps);
1028       }
1029     }
1030
1031     int column = 0;
1032     for (AlignedCodon codon : alignedCodons.keySet())
1033     {
1034       final Map<SequenceI, AlignedCodon> columnResidues = alignedCodons
1035               .get(codon);
1036       for (Entry<SequenceI, AlignedCodon> entry : columnResidues.entrySet())
1037       {
1038         // place translated codon at its column position in sequence
1039         entry.getKey().getSequence()[column] = entry.getValue().product
1040                 .charAt(0);
1041       }
1042       column++;
1043     }
1044     return 0;
1045   }
1046
1047   /**
1048    * Populate the map of aligned codons by traversing the given sequence
1049    * mapping, locating the aligned positions of mapped codons, and adding those
1050    * positions and their translation products to the map.
1051    *
1052    * @param dna
1053    *          the aligned sequence we are mapping from
1054    * @param protein
1055    *          the sequence to be aligned to the codons
1056    * @param gapChar
1057    *          the gap character in the dna sequence
1058    * @param seqMap
1059    *          a mapping to a sequence translation
1060    * @param alignedCodons
1061    *          the map we are building up
1062    */
1063   static void addCodonPositions(SequenceI dna, SequenceI protein,
1064           char gapChar, Mapping seqMap,
1065           Map<AlignedCodon, Map<SequenceI, AlignedCodon>> alignedCodons)
1066   {
1067     Iterator<AlignedCodon> codons = seqMap.getCodonIterator(dna, gapChar);
1068
1069     /*
1070      * add codon positions, and their peptide translations, to the alignment
1071      * map, while remembering the first codon mapped
1072      */
1073     while (codons.hasNext())
1074     {
1075       try
1076       {
1077         AlignedCodon codon = codons.next();
1078         addCodonToMap(alignedCodons, codon, protein);
1079       } catch (IncompleteCodonException e)
1080       {
1081         // possible incomplete trailing codon - ignore
1082       } catch (NoSuchElementException e)
1083       {
1084         // possibly peptide lacking STOP
1085       }
1086     }
1087   }
1088
1089   /**
1090    * Helper method to add a codon-to-peptide entry to the aligned codons map
1091    *
1092    * @param alignedCodons
1093    * @param codon
1094    * @param protein
1095    */
1096   protected static void addCodonToMap(
1097           Map<AlignedCodon, Map<SequenceI, AlignedCodon>> alignedCodons,
1098           AlignedCodon codon, SequenceI protein)
1099   {
1100     Map<SequenceI, AlignedCodon> seqProduct = alignedCodons.get(codon);
1101     if (seqProduct == null)
1102     {
1103       seqProduct = new HashMap<SequenceI, AlignedCodon>();
1104       alignedCodons.put(codon, seqProduct);
1105     }
1106     seqProduct.put(protein, codon);
1107   }
1108
1109   /**
1110    * Returns true if a cDNA/Protein mapping either exists, or could be made,
1111    * between at least one pair of sequences in the two alignments. Currently,
1112    * the logic is:
1113    * <ul>
1114    * <li>One alignment must be nucleotide, and the other protein</li>
1115    * <li>At least one pair of sequences must be already mapped, or mappable</li>
1116    * <li>Mappable means the nucleotide translation matches the protein sequence</li>
1117    * <li>The translation may ignore start and stop codons if present in the
1118    * nucleotide</li>
1119    * </ul>
1120    *
1121    * @param al1
1122    * @param al2
1123    * @return
1124    */
1125   public static boolean isMappable(AlignmentI al1, AlignmentI al2)
1126   {
1127     if (al1 == null || al2 == null)
1128     {
1129       return false;
1130     }
1131
1132     /*
1133      * Require one nucleotide and one protein
1134      */
1135     if (al1.isNucleotide() == al2.isNucleotide())
1136     {
1137       return false;
1138     }
1139     AlignmentI dna = al1.isNucleotide() ? al1 : al2;
1140     AlignmentI protein = dna == al1 ? al2 : al1;
1141     List<AlignedCodonFrame> mappings = protein.getCodonFrames();
1142     for (SequenceI dnaSeq : dna.getSequences())
1143     {
1144       for (SequenceI proteinSeq : protein.getSequences())
1145       {
1146         if (isMappable(dnaSeq, proteinSeq, mappings))
1147         {
1148           return true;
1149         }
1150       }
1151     }
1152     return false;
1153   }
1154
1155   /**
1156    * Returns true if the dna sequence is mapped, or could be mapped, to the
1157    * protein sequence.
1158    *
1159    * @param dnaSeq
1160    * @param proteinSeq
1161    * @param mappings
1162    * @return
1163    */
1164   protected static boolean isMappable(SequenceI dnaSeq,
1165           SequenceI proteinSeq, List<AlignedCodonFrame> mappings)
1166   {
1167     if (dnaSeq == null || proteinSeq == null)
1168     {
1169       return false;
1170     }
1171
1172     SequenceI dnaDs = dnaSeq.getDatasetSequence() == null ? dnaSeq : dnaSeq
1173             .getDatasetSequence();
1174     SequenceI proteinDs = proteinSeq.getDatasetSequence() == null ? proteinSeq
1175             : proteinSeq.getDatasetSequence();
1176
1177     for (AlignedCodonFrame mapping : mappings)
1178     {
1179       if (proteinDs == mapping.getAaForDnaSeq(dnaDs))
1180       {
1181         /*
1182          * already mapped
1183          */
1184         return true;
1185       }
1186     }
1187
1188     /*
1189      * Just try to make a mapping (it is not yet stored), test whether
1190      * successful.
1191      */
1192     return mapCdnaToProtein(proteinDs, dnaDs) != null;
1193   }
1194
1195   /**
1196    * Finds any reference annotations associated with the sequences in
1197    * sequenceScope, that are not already added to the alignment, and adds them
1198    * to the 'candidates' map. Also populates a lookup table of annotation
1199    * labels, keyed by calcId, for use in constructing tooltips or the like.
1200    *
1201    * @param sequenceScope
1202    *          the sequences to scan for reference annotations
1203    * @param labelForCalcId
1204    *          (optional) map to populate with label for calcId
1205    * @param candidates
1206    *          map to populate with annotations for sequence
1207    * @param al
1208    *          the alignment to check for presence of annotations
1209    */
1210   public static void findAddableReferenceAnnotations(
1211           List<SequenceI> sequenceScope,
1212           Map<String, String> labelForCalcId,
1213           final Map<SequenceI, List<AlignmentAnnotation>> candidates,
1214           AlignmentI al)
1215   {
1216     if (sequenceScope == null)
1217     {
1218       return;
1219     }
1220
1221     /*
1222      * For each sequence in scope, make a list of any annotations on the
1223      * underlying dataset sequence which are not already on the alignment.
1224      *
1225      * Add to a map of { alignmentSequence, <List of annotations to add> }
1226      */
1227     for (SequenceI seq : sequenceScope)
1228     {
1229       SequenceI dataset = seq.getDatasetSequence();
1230       if (dataset == null)
1231       {
1232         continue;
1233       }
1234       AlignmentAnnotation[] datasetAnnotations = dataset.getAnnotation();
1235       if (datasetAnnotations == null)
1236       {
1237         continue;
1238       }
1239       final List<AlignmentAnnotation> result = new ArrayList<AlignmentAnnotation>();
1240       for (AlignmentAnnotation dsann : datasetAnnotations)
1241       {
1242         /*
1243          * Find matching annotations on the alignment. If none is found, then
1244          * add this annotation to the list of 'addable' annotations for this
1245          * sequence.
1246          */
1247         final Iterable<AlignmentAnnotation> matchedAlignmentAnnotations = al
1248                 .findAnnotations(seq, dsann.getCalcId(), dsann.label);
1249         if (!matchedAlignmentAnnotations.iterator().hasNext())
1250         {
1251           result.add(dsann);
1252           if (labelForCalcId != null)
1253           {
1254             labelForCalcId.put(dsann.getCalcId(), dsann.label);
1255           }
1256         }
1257       }
1258       /*
1259        * Save any addable annotations for this sequence
1260        */
1261       if (!result.isEmpty())
1262       {
1263         candidates.put(seq, result);
1264       }
1265     }
1266   }
1267
1268   /**
1269    * Adds annotations to the top of the alignment annotations, in the same order
1270    * as their related sequences.
1271    *
1272    * @param annotations
1273    *          the annotations to add
1274    * @param alignment
1275    *          the alignment to add them to
1276    * @param selectionGroup
1277    *          current selection group (or null if none)
1278    */
1279   public static void addReferenceAnnotations(
1280           Map<SequenceI, List<AlignmentAnnotation>> annotations,
1281           final AlignmentI alignment, final SequenceGroup selectionGroup)
1282   {
1283     for (SequenceI seq : annotations.keySet())
1284     {
1285       for (AlignmentAnnotation ann : annotations.get(seq))
1286       {
1287         AlignmentAnnotation copyAnn = new AlignmentAnnotation(ann);
1288         int startRes = 0;
1289         int endRes = ann.annotations.length;
1290         if (selectionGroup != null)
1291         {
1292           startRes = selectionGroup.getStartRes();
1293           endRes = selectionGroup.getEndRes();
1294         }
1295         copyAnn.restrict(startRes, endRes);
1296
1297         /*
1298          * Add to the sequence (sets copyAnn.datasetSequence), unless the
1299          * original annotation is already on the sequence.
1300          */
1301         if (!seq.hasAnnotation(ann))
1302         {
1303           seq.addAlignmentAnnotation(copyAnn);
1304         }
1305         // adjust for gaps
1306         copyAnn.adjustForAlignment();
1307         // add to the alignment and set visible
1308         alignment.addAnnotation(copyAnn);
1309         copyAnn.visible = true;
1310       }
1311     }
1312   }
1313
1314   /**
1315    * Set visibility of alignment annotations of specified types (labels), for
1316    * specified sequences. This supports controls like
1317    * "Show all secondary structure", "Hide all Temp factor", etc.
1318    *
1319    * @al the alignment to scan for annotations
1320    * @param types
1321    *          the types (labels) of annotations to be updated
1322    * @param forSequences
1323    *          if not null, only annotations linked to one of these sequences are
1324    *          in scope for update; if null, acts on all sequence annotations
1325    * @param anyType
1326    *          if this flag is true, 'types' is ignored (label not checked)
1327    * @param doShow
1328    *          if true, set visibility on, else set off
1329    */
1330   public static void showOrHideSequenceAnnotations(AlignmentI al,
1331           Collection<String> types, List<SequenceI> forSequences,
1332           boolean anyType, boolean doShow)
1333   {
1334     for (AlignmentAnnotation aa : al.getAlignmentAnnotation())
1335     {
1336       if (anyType || types.contains(aa.label))
1337       {
1338         if ((aa.sequenceRef != null)
1339                 && (forSequences == null || forSequences
1340                         .contains(aa.sequenceRef)))
1341         {
1342           aa.visible = doShow;
1343         }
1344       }
1345     }
1346   }
1347
1348   /**
1349    * Returns true if either sequence has a cross-reference to the other
1350    *
1351    * @param seq1
1352    * @param seq2
1353    * @return
1354    */
1355   public static boolean haveCrossRef(SequenceI seq1, SequenceI seq2)
1356   {
1357     // Note: moved here from class CrossRef as the latter class has dependencies
1358     // not availability to the applet's classpath
1359     return hasCrossRef(seq1, seq2) || hasCrossRef(seq2, seq1);
1360   }
1361
1362   /**
1363    * Returns true if seq1 has a cross-reference to seq2. Currently this assumes
1364    * that sequence name is structured as Source|AccessionId.
1365    *
1366    * @param seq1
1367    * @param seq2
1368    * @return
1369    */
1370   public static boolean hasCrossRef(SequenceI seq1, SequenceI seq2)
1371   {
1372     if (seq1 == null || seq2 == null)
1373     {
1374       return false;
1375     }
1376     String name = seq2.getName();
1377     final DBRefEntry[] xrefs = seq1.getDBRefs();
1378     if (xrefs != null)
1379     {
1380       for (DBRefEntry xref : xrefs)
1381       {
1382         String xrefName = xref.getSource() + "|" + xref.getAccessionId();
1383         // case-insensitive test, consistent with DBRefEntry.equalRef()
1384         if (xrefName.equalsIgnoreCase(name))
1385         {
1386           return true;
1387         }
1388       }
1389     }
1390     return false;
1391   }
1392
1393   /**
1394    * Constructs an alignment consisting of the mapped (CDS) regions in the given
1395    * nucleotide sequences, and updates mappings to match. The CDS sequences are
1396    * added to the original alignment's dataset, which is shared by the new
1397    * alignment. Mappings from nucleotide to CDS, and from CDS to protein, are
1398    * added to the alignment dataset.
1399    *
1400    * @param dna
1401    *          aligned dna sequences
1402    * @param mappings
1403    *          from dna to protein
1404    * @param al
1405    * @return an alignment whose sequences are the cds-only parts of the dna
1406    *         sequences (or null if no mappings are found)
1407    */
1408   public static AlignmentI makeCdsAlignment(SequenceI[] dna,
1409           List<AlignedCodonFrame> mappings, AlignmentI al)
1410   {
1411     List<SequenceI> cdsSeqs = new ArrayList<SequenceI>();
1412
1413     for (SequenceI seq : dna)
1414     {
1415       AlignedCodonFrame cdsMappings = new AlignedCodonFrame();
1416       List<AlignedCodonFrame> seqMappings = MappingUtils
1417               .findMappingsForSequence(seq, mappings);
1418       List<AlignedCodonFrame> alignmentMappings = al.getCodonFrames();
1419       for (AlignedCodonFrame mapping : seqMappings)
1420       {
1421         for (Mapping aMapping : mapping.getMappingsFromSequence(seq))
1422         {
1423           SequenceI cdsSeq = makeCdsSequence(seq.getDatasetSequence(),
1424                   aMapping);
1425           cdsSeqs.add(cdsSeq);
1426
1427           /*
1428            * add a mapping from CDS to the (unchanged) mapped to range
1429            */
1430           List<int[]> cdsRange = Collections.singletonList(new int[] { 1,
1431               cdsSeq.getLength() });
1432           MapList map = new MapList(cdsRange, aMapping.getMap()
1433                   .getToRanges(), aMapping.getMap().getFromRatio(),
1434                   aMapping.getMap().getToRatio());
1435           cdsMappings.addMap(cdsSeq, aMapping.getTo(), map);
1436
1437           /*
1438            * add another mapping from original 'from' range to CDS
1439            */
1440           map = new MapList(aMapping.getMap().getFromRanges(), cdsRange, 1,
1441                   1);
1442           cdsMappings.addMap(seq.getDatasetSequence(), cdsSeq, map);
1443
1444           alignmentMappings.add(cdsMappings);
1445
1446           /*
1447            * transfer any features on dna that overlap the CDS
1448            */
1449           transferFeatures(seq, cdsSeq, map, null, SequenceOntologyI.CDS);
1450         }
1451       }
1452     }
1453
1454     /*
1455      * add CDS seqs to shared dataset
1456      */
1457     Alignment dataset = al.getDataset();
1458     for (SequenceI seq : cdsSeqs)
1459     {
1460       if (!dataset.getSequences().contains(seq.getDatasetSequence()))
1461       {
1462         dataset.addSequence(seq.getDatasetSequence());
1463       }
1464     }
1465     AlignmentI cds = new Alignment(cdsSeqs.toArray(new SequenceI[cdsSeqs
1466             .size()]));
1467     cds.setDataset(dataset);
1468
1469     return cds;
1470   }
1471
1472   /**
1473    * Helper method that makes a CDS sequence as defined by the mappings from the
1474    * given sequence i.e. extracts the 'mapped from' ranges (which may be on
1475    * forward or reverse strand).
1476    *
1477    * @param seq
1478    * @param mapping
1479    * @return
1480    */
1481   static SequenceI makeCdsSequence(SequenceI seq, Mapping mapping)
1482   {
1483     char[] seqChars = seq.getSequence();
1484     List<int[]> fromRanges = mapping.getMap().getFromRanges();
1485     int cdsWidth = MappingUtils.getLength(fromRanges);
1486     char[] newSeqChars = new char[cdsWidth];
1487
1488     int newPos = 0;
1489     for (int[] range : fromRanges)
1490     {
1491       if (range[0] <= range[1])
1492       {
1493         // forward strand mapping - just copy the range
1494         int length = range[1] - range[0] + 1;
1495         System.arraycopy(seqChars, range[0] - 1, newSeqChars, newPos,
1496                 length);
1497         newPos += length;
1498       }
1499       else
1500       {
1501         // reverse strand mapping - copy and complement one by one
1502         for (int i = range[0]; i >= range[1]; i--)
1503         {
1504           newSeqChars[newPos++] = Dna.getComplement(seqChars[i - 1]);
1505         }
1506       }
1507     }
1508
1509     SequenceI newSeq = new Sequence(seq.getName() + "|"
1510             + mapping.getTo().getName(), newSeqChars, 1, newPos);
1511     newSeq.createDatasetSequence();
1512     return newSeq;
1513   }
1514
1515   /**
1516    * Transfers co-located features on 'fromSeq' to 'toSeq', adjusting the
1517    * feature start/end ranges, optionally omitting specified feature types.
1518    * Returns the number of features copied.
1519    *
1520    * @param fromSeq
1521    * @param toSeq
1522    * @param select
1523    *          if not null, only features of this type are copied (including
1524    *          subtypes in the Sequence Ontology)
1525    * @param mapping
1526    *          the mapping from 'fromSeq' to 'toSeq'
1527    * @param omitting
1528    */
1529   public static int transferFeatures(SequenceI fromSeq, SequenceI toSeq,
1530           MapList mapping, String select, String... omitting)
1531   {
1532     SequenceI copyTo = toSeq;
1533     while (copyTo.getDatasetSequence() != null)
1534     {
1535       copyTo = copyTo.getDatasetSequence();
1536     }
1537
1538     SequenceOntologyI so = SequenceOntologyFactory.getInstance();
1539     int count = 0;
1540     SequenceFeature[] sfs = fromSeq.getSequenceFeatures();
1541     if (sfs != null)
1542     {
1543       for (SequenceFeature sf : sfs)
1544       {
1545         String type = sf.getType();
1546         if (select != null && !so.isA(type, select))
1547         {
1548           continue;
1549         }
1550         boolean omit = false;
1551         for (String toOmit : omitting)
1552         {
1553           if (type.equals(toOmit))
1554           {
1555             omit = true;
1556           }
1557         }
1558         if (omit)
1559         {
1560           continue;
1561         }
1562
1563         /*
1564          * locate the mapped range - null if either start or end is
1565          * not mapped (no partial overlaps are calculated)
1566          */
1567         int start = sf.getBegin();
1568         int end = sf.getEnd();
1569         int[] mappedTo = mapping.locateInTo(start, end);
1570         /*
1571          * if whole exon range doesn't map, try interpreting it
1572          * as 5' or 3' exon overlapping the CDS range
1573          */
1574         if (mappedTo == null)
1575         {
1576           mappedTo = mapping.locateInTo(end, end);
1577           if (mappedTo != null)
1578           {
1579             /*
1580              * end of exon is in CDS range - 5' overlap
1581              * to a range from the start of the peptide
1582              */
1583             mappedTo[0] = 1;
1584           }
1585         }
1586         if (mappedTo == null)
1587         {
1588           mappedTo = mapping.locateInTo(start, start);
1589           if (mappedTo != null)
1590           {
1591             /*
1592              * start of exon is in CDS range - 3' overlap
1593              * to a range up to the end of the peptide
1594              */
1595             mappedTo[1] = toSeq.getLength();
1596           }
1597         }
1598         if (mappedTo != null)
1599         {
1600           SequenceFeature copy = new SequenceFeature(sf);
1601           copy.setBegin(Math.min(mappedTo[0], mappedTo[1]));
1602           copy.setEnd(Math.max(mappedTo[0], mappedTo[1]));
1603           copyTo.addSequenceFeature(copy);
1604           count++;
1605         }
1606       }
1607     }
1608     return count;
1609   }
1610
1611   /**
1612    * Returns a mapping from dna to protein by inspecting sequence features of
1613    * type "CDS" on the dna.
1614    *
1615    * @param dnaSeq
1616    * @param proteinSeq
1617    * @return
1618    */
1619   public static MapList mapCdsToProtein(SequenceI dnaSeq,
1620           SequenceI proteinSeq)
1621   {
1622     List<int[]> ranges = findCdsPositions(dnaSeq);
1623     int mappedDnaLength = MappingUtils.getLength(ranges);
1624
1625     int proteinLength = proteinSeq.getLength();
1626     int proteinStart = proteinSeq.getStart();
1627     int proteinEnd = proteinSeq.getEnd();
1628
1629     /*
1630      * incomplete start codon may mean X at start of peptide
1631      * we ignore both for mapping purposes
1632      */
1633     if (proteinSeq.getCharAt(0) == 'X')
1634     {
1635       // todo JAL-2022 support startPhase > 0
1636       proteinStart++;
1637       proteinLength--;
1638     }
1639     List<int[]> proteinRange = new ArrayList<int[]>();
1640
1641     /*
1642      * dna length should map to protein (or protein plus stop codon)
1643      */
1644     int codesForResidues = mappedDnaLength / 3;
1645     if (codesForResidues == (proteinLength + 1))
1646     {
1647       // assuming extra codon is for STOP and not in peptide
1648       codesForResidues--;
1649     }
1650     if (codesForResidues == proteinLength)
1651     {
1652       proteinRange.add(new int[] { proteinStart, proteinEnd });
1653       return new MapList(ranges, proteinRange, 3, 1);
1654     }
1655     return null;
1656   }
1657
1658   /**
1659    * Returns a list of CDS ranges found (as sequence positions base 1), i.e. of
1660    * start/end positions of sequence features of type "CDS" (or a sub-type of
1661    * CDS in the Sequence Ontology). The ranges are sorted into ascending start
1662    * position order, so this method is only valid for linear CDS in the same
1663    * sense as the protein product.
1664    *
1665    * @param dnaSeq
1666    * @return
1667    */
1668   public static List<int[]> findCdsPositions(SequenceI dnaSeq)
1669   {
1670     List<int[]> result = new ArrayList<int[]>();
1671     SequenceFeature[] sfs = dnaSeq.getSequenceFeatures();
1672     if (sfs == null)
1673     {
1674       return result;
1675     }
1676
1677     SequenceOntologyI so = SequenceOntologyFactory.getInstance();
1678     int startPhase = 0;
1679
1680     for (SequenceFeature sf : sfs)
1681     {
1682       /*
1683        * process a CDS feature (or a sub-type of CDS)
1684        */
1685       if (so.isA(sf.getType(), SequenceOntologyI.CDS))
1686       {
1687         int phase = 0;
1688         try
1689         {
1690           phase = Integer.parseInt(sf.getPhase());
1691         } catch (NumberFormatException e)
1692         {
1693           // ignore
1694         }
1695         /*
1696          * phase > 0 on first codon means 5' incomplete - skip to the start
1697          * of the next codon; example ENST00000496384
1698          */
1699         int begin = sf.getBegin();
1700         int end = sf.getEnd();
1701         if (result.isEmpty())
1702         {
1703           begin += phase;
1704           if (begin > end)
1705           {
1706             // shouldn't happen!
1707             System.err
1708                     .println("Error: start phase extends beyond start CDS in "
1709                             + dnaSeq.getName());
1710           }
1711         }
1712         result.add(new int[] { begin, end });
1713       }
1714     }
1715
1716     /*
1717      * remove 'startPhase' positions (usually 0) from the first range
1718      * so we begin at the start of a complete codon
1719      */
1720     if (!result.isEmpty())
1721     {
1722       // TODO JAL-2022 correctly model start phase > 0
1723       result.get(0)[0] += startPhase;
1724     }
1725
1726     /*
1727      * Finally sort ranges by start position. This avoids a dependency on
1728      * keeping features in order on the sequence (if they are in order anyway,
1729      * the sort will have almost no work to do). The implicit assumption is CDS
1730      * ranges are assembled in order. Other cases should not use this method,
1731      * but instead construct an explicit mapping for CDS (e.g. EMBL parsing).
1732      */
1733     Collections.sort(result, new Comparator<int[]>()
1734     {
1735       @Override
1736       public int compare(int[] o1, int[] o2)
1737       {
1738         return Integer.compare(o1[0], o2[0]);
1739       }
1740     });
1741     return result;
1742   }
1743
1744   /**
1745    * Maps exon features from dna to protein, and computes variants in peptide
1746    * product generated by variants in dna, and adds them as sequence_variant
1747    * features on the protein sequence. Returns the number of variant features
1748    * added.
1749    *
1750    * @param dnaSeq
1751    * @param peptide
1752    * @param dnaToProtein
1753    */
1754   public static int computeProteinFeatures(SequenceI dnaSeq,
1755           SequenceI peptide, MapList dnaToProtein)
1756   {
1757     while (dnaSeq.getDatasetSequence() != null)
1758     {
1759       dnaSeq = dnaSeq.getDatasetSequence();
1760     }
1761     while (peptide.getDatasetSequence() != null)
1762     {
1763       peptide = peptide.getDatasetSequence();
1764     }
1765
1766     transferFeatures(dnaSeq, peptide, dnaToProtein, SequenceOntologyI.EXON);
1767
1768     /*
1769      * compute protein variants from dna variants and codon mappings;
1770      * NB - alternatively we could retrieve this using the REST service e.g.
1771      * http://rest.ensembl.org/overlap/translation
1772      * /ENSP00000288602?feature=transcript_variation;content-type=text/xml
1773      * which would be a bit slower but possibly more reliable
1774      */
1775
1776     /*
1777      * build a map with codon variations for each potentially varying peptide
1778      */
1779     LinkedHashMap<Integer, List<DnaVariant>[]> variants = buildDnaVariantsMap(
1780             dnaSeq, dnaToProtein);
1781
1782     /*
1783      * scan codon variations, compute peptide variants and add to peptide sequence
1784      */
1785     int count = 0;
1786     for (Entry<Integer, List<DnaVariant>[]> variant : variants.entrySet())
1787     {
1788       int peptidePos = variant.getKey();
1789       List<DnaVariant>[] codonVariants = variant.getValue();
1790       count += computePeptideVariants(peptide, peptidePos, codonVariants);
1791     }
1792
1793     /*
1794      * sort to get sequence features in start position order
1795      * - would be better to store in Sequence as a TreeSet or NCList?
1796      */
1797     Arrays.sort(peptide.getSequenceFeatures(),
1798             new Comparator<SequenceFeature>()
1799             {
1800               @Override
1801               public int compare(SequenceFeature o1, SequenceFeature o2)
1802               {
1803                 int c = Integer.compare(o1.getBegin(), o2.getBegin());
1804                 return c == 0 ? Integer.compare(o1.getEnd(), o2.getEnd())
1805                         : c;
1806               }
1807             });
1808     return count;
1809   }
1810
1811   /**
1812    * Computes non-synonymous peptide variants from codon variants and adds them
1813    * as sequence_variant features on the protein sequence (one feature per
1814    * allele variant). Selected attributes (variant id, clinical significance)
1815    * are copied over to the new features.
1816    *
1817    * @param peptide
1818    *          the protein sequence
1819    * @param peptidePos
1820    *          the position to compute peptide variants for
1821    * @param codonVariants
1822    *          a list of dna variants per codon position
1823    * @return the number of features added
1824    */
1825   static int computePeptideVariants(SequenceI peptide, int peptidePos,
1826           List<DnaVariant>[] codonVariants)
1827   {
1828     String residue = String.valueOf(peptide.getCharAt(peptidePos - 1));
1829     int count = 0;
1830     String base1 = codonVariants[0].get(0).base;
1831     String base2 = codonVariants[1].get(0).base;
1832     String base3 = codonVariants[2].get(0).base;
1833
1834     /*
1835      * variants in first codon base
1836      */
1837     for (DnaVariant var : codonVariants[0])
1838     {
1839       if (var.variant != null)
1840       {
1841         String alleles = (String) var.variant.getValue("alleles");
1842         if (alleles != null)
1843         {
1844           for (String base : alleles.split(","))
1845           {
1846             String codon = base + base2 + base3;
1847             if (addPeptideVariant(peptide, peptidePos, residue, var, codon))
1848             {
1849               count++;
1850             }
1851           }
1852         }
1853       }
1854     }
1855
1856     /*
1857      * variants in second codon base
1858      */
1859     for (DnaVariant var : codonVariants[1])
1860     {
1861       if (var.variant != null)
1862       {
1863         String alleles = (String) var.variant.getValue("alleles");
1864         if (alleles != null)
1865         {
1866           for (String base : alleles.split(","))
1867           {
1868             String codon = base1 + base + base3;
1869             if (addPeptideVariant(peptide, peptidePos, residue, var, codon))
1870             {
1871               count++;
1872             }
1873           }
1874         }
1875       }
1876     }
1877
1878     /*
1879      * variants in third codon base
1880      */
1881     for (DnaVariant var : codonVariants[2])
1882     {
1883       if (var.variant != null)
1884       {
1885         String alleles = (String) var.variant.getValue("alleles");
1886         if (alleles != null)
1887         {
1888           for (String base : alleles.split(","))
1889           {
1890             String codon = base1 + base2 + base;
1891             if (addPeptideVariant(peptide, peptidePos, residue, var, codon))
1892             {
1893               count++;
1894             }
1895           }
1896         }
1897       }
1898     }
1899
1900     return count;
1901   }
1902
1903   /**
1904    * Helper method that adds a peptide variant feature, provided the given codon
1905    * translates to a value different to the current residue (is a non-synonymous
1906    * variant). ID and clinical_significance attributes of the dna variant (if
1907    * present) are copied to the new feature.
1908    *
1909    * @param peptide
1910    * @param peptidePos
1911    * @param residue
1912    * @param var
1913    * @param codon
1914    * @return true if a feature was added, else false
1915    */
1916   static boolean addPeptideVariant(SequenceI peptide, int peptidePos,
1917           String residue, DnaVariant var, String codon)
1918   {
1919     /*
1920      * get peptide translation of codon e.g. GAT -> D
1921      * note that variants which are not single alleles,
1922      * e.g. multibase variants or HGMD_MUTATION etc
1923      * are currently ignored here
1924      */
1925     String trans = codon.contains("-") ? "-"
1926             : (codon.length() > 3 ? null : ResidueProperties
1927                     .codonTranslate(codon));
1928     if (trans != null && !trans.equals(residue))
1929     {
1930       String desc = residue + "->" + trans;
1931       // set score to 0f so 'graduated colour' option is offered! JAL-2060
1932       SequenceFeature sf = new SequenceFeature(
1933               SequenceOntologyI.SEQUENCE_VARIANT, desc, peptidePos,
1934               peptidePos, 0f, null);
1935       StringBuilder attributes = new StringBuilder(32);
1936       String id = (String) var.variant.getValue(ID);
1937       if (id != null)
1938       {
1939         if (id.startsWith(SEQUENCE_VARIANT))
1940         {
1941           id = id.substring(SEQUENCE_VARIANT.length());
1942         }
1943         sf.setValue(ID, id);
1944         attributes.append(ID).append("=").append(id);
1945         // TODO handle other species variants
1946         StringBuilder link = new StringBuilder(32);
1947         try
1948         {
1949           link.append(desc).append(" ").append(id)
1950                   .append("|http://www.ensembl.org/Homo_sapiens/Variation/Summary?v=")
1951                   .append(URLEncoder.encode(id, "UTF-8"));
1952           sf.addLink(link.toString());
1953         } catch (UnsupportedEncodingException e)
1954         {
1955           // as if
1956         }
1957       }
1958       String clinSig = (String) var.variant
1959               .getValue(CLINICAL_SIGNIFICANCE);
1960       if (clinSig != null)
1961       {
1962         sf.setValue(CLINICAL_SIGNIFICANCE, clinSig);
1963         attributes.append(";").append(CLINICAL_SIGNIFICANCE).append("=")
1964                 .append(clinSig);
1965       }
1966       peptide.addSequenceFeature(sf);
1967       if (attributes.length() > 0)
1968       {
1969         sf.setAttributes(attributes.toString());
1970       }
1971       return true;
1972     }
1973     return false;
1974   }
1975
1976   /**
1977    * Builds a map whose key is position in the protein sequence, and value is a
1978    * list of the base and all variants for each corresponding codon position
1979    *
1980    * @param dnaSeq
1981    * @param dnaToProtein
1982    * @return
1983    */
1984   static LinkedHashMap<Integer, List<DnaVariant>[]> buildDnaVariantsMap(
1985           SequenceI dnaSeq, MapList dnaToProtein)
1986   {
1987     /*
1988      * map from peptide position to all variants of the codon which codes for it
1989      * LinkedHashMap ensures we keep the peptide features in sequence order
1990      */
1991     LinkedHashMap<Integer, List<DnaVariant>[]> variants = new LinkedHashMap<Integer, List<DnaVariant>[]>();
1992     SequenceOntologyI so = SequenceOntologyFactory.getInstance();
1993
1994     SequenceFeature[] dnaFeatures = dnaSeq.getSequenceFeatures();
1995     if (dnaFeatures == null)
1996     {
1997       return variants;
1998     }
1999
2000     int dnaStart = dnaSeq.getStart();
2001     int[] lastCodon = null;
2002     int lastPeptidePostion = 0;
2003
2004     /*
2005      * build a map of codon variations for peptides
2006      */
2007     for (SequenceFeature sf : dnaFeatures)
2008     {
2009       int dnaCol = sf.getBegin();
2010       if (dnaCol != sf.getEnd())
2011       {
2012         // not handling multi-locus variant features
2013         continue;
2014       }
2015       if (so.isA(sf.getType(), SequenceOntologyI.SEQUENCE_VARIANT))
2016       {
2017         int[] mapsTo = dnaToProtein.locateInTo(dnaCol, dnaCol);
2018         if (mapsTo == null)
2019         {
2020           // feature doesn't lie within coding region
2021           continue;
2022         }
2023         int peptidePosition = mapsTo[0];
2024         List<DnaVariant>[] codonVariants = variants.get(peptidePosition);
2025         if (codonVariants == null)
2026         {
2027           codonVariants = new ArrayList[3];
2028           codonVariants[0] = new ArrayList<DnaVariant>();
2029           codonVariants[1] = new ArrayList<DnaVariant>();
2030           codonVariants[2] = new ArrayList<DnaVariant>();
2031           variants.put(peptidePosition, codonVariants);
2032         }
2033
2034         /*
2035          * extract dna variants to a string array
2036          */
2037         String alls = (String) sf.getValue("alleles");
2038         if (alls == null)
2039         {
2040           continue;
2041         }
2042         String[] alleles = alls.toUpperCase().split(",");
2043         int i = 0;
2044         for (String allele : alleles)
2045         {
2046           alleles[i++] = allele.trim(); // lose any space characters "A, G"
2047         }
2048
2049         /*
2050          * get this peptide's codon positions e.g. [3, 4, 5] or [4, 7, 10]
2051          */
2052         int[] codon = peptidePosition == lastPeptidePostion ? lastCodon
2053                 : MappingUtils.flattenRanges(dnaToProtein.locateInFrom(
2054                         peptidePosition, peptidePosition));
2055         lastPeptidePostion = peptidePosition;
2056         lastCodon = codon;
2057
2058         /*
2059          * save nucleotide (and any variant) for each codon position
2060          */
2061         for (int codonPos = 0; codonPos < 3; codonPos++)
2062         {
2063           String nucleotide = String.valueOf(
2064                   dnaSeq.getCharAt(codon[codonPos] - dnaStart))
2065                   .toUpperCase();
2066           List<DnaVariant> codonVariant = codonVariants[codonPos];
2067           if (codon[codonPos] == dnaCol)
2068           {
2069             if (!codonVariant.isEmpty()
2070                     && codonVariant.get(0).variant == null)
2071             {
2072               /*
2073                * already recorded base value, add this variant
2074                */
2075               codonVariant.get(0).variant = sf;
2076             }
2077             else
2078             {
2079               /*
2080                * add variant with base value
2081                */
2082               codonVariant.add(new DnaVariant(nucleotide, sf));
2083             }
2084           }
2085           else if (codonVariant.isEmpty())
2086           {
2087             /*
2088              * record (possibly non-varying) base value
2089              */
2090             codonVariant.add(new DnaVariant(nucleotide));
2091           }
2092         }
2093       }
2094     }
2095     return variants;
2096   }
2097
2098   /**
2099    * Makes an alignment with a copy of the given sequences, adding in any
2100    * non-redundant sequences which are mapped to by the cross-referenced
2101    * sequences.
2102    *
2103    * @param seqs
2104    * @param xrefs
2105    * @return
2106    */
2107   public static AlignmentI makeCopyAlignment(SequenceI[] seqs,
2108           SequenceI[] xrefs)
2109   {
2110     AlignmentI copy = new Alignment(new Alignment(seqs));
2111
2112     SequenceIdMatcher matcher = new SequenceIdMatcher(seqs);
2113     if (xrefs != null)
2114     {
2115       for (SequenceI xref : xrefs)
2116       {
2117         DBRefEntry[] dbrefs = xref.getDBRefs();
2118         if (dbrefs != null)
2119         {
2120           for (DBRefEntry dbref : dbrefs)
2121           {
2122             if (dbref.getMap() == null || dbref.getMap().getTo() == null)
2123             {
2124               continue;
2125             }
2126             SequenceI mappedTo = dbref.getMap().getTo();
2127             SequenceI match = matcher.findIdMatch(mappedTo);
2128             if (match == null)
2129             {
2130               matcher.add(mappedTo);
2131               copy.addSequence(mappedTo);
2132             }
2133           }
2134         }
2135       }
2136     }
2137     return copy;
2138   }
2139
2140   /**
2141    * Try to align sequences in 'unaligned' to match the alignment of their
2142    * mapped regions in 'aligned'. For example, could use this to align CDS
2143    * sequences which are mapped to their parent cDNA sequences.
2144    *
2145    * This method handles 1:1 mappings (dna-to-dna or protein-to-protein). For
2146    * dna-to-protein or protein-to-dna use alternative methods.
2147    *
2148    * @param unaligned
2149    *          sequences to be aligned
2150    * @param aligned
2151    *          holds aligned sequences and their mappings
2152    * @return
2153    */
2154   public static int alignAs(AlignmentI unaligned, AlignmentI aligned)
2155   {
2156     List<SequenceI> unmapped = new ArrayList<SequenceI>();
2157     Map<Integer, Map<SequenceI, Character>> columnMap = buildMappedColumnsMap(
2158             unaligned, aligned, unmapped);
2159     int width = columnMap.size();
2160     char gap = unaligned.getGapCharacter();
2161     int realignedCount = 0;
2162
2163     for (SequenceI seq : unaligned.getSequences())
2164     {
2165       if (!unmapped.contains(seq))
2166       {
2167         char[] newSeq = new char[width];
2168         Arrays.fill(newSeq, gap);
2169         int newCol = 0;
2170         int lastCol = 0;
2171
2172         /*
2173          * traverse the map to find columns populated
2174          * by our sequence
2175          */
2176         for (Integer column : columnMap.keySet())
2177         {
2178           Character c = columnMap.get(column).get(seq);
2179           if (c != null)
2180           {
2181             /*
2182              * sequence has a character at this position
2183              *
2184              */
2185             newSeq[newCol] = c;
2186             lastCol = newCol;
2187           }
2188           newCol++;
2189         }
2190
2191         /*
2192          * trim trailing gaps
2193          */
2194         if (lastCol < width)
2195         {
2196           char[] tmp = new char[lastCol + 1];
2197           System.arraycopy(newSeq, 0, tmp, 0, lastCol + 1);
2198           newSeq = tmp;
2199         }
2200         seq.setSequence(String.valueOf(newSeq));
2201         realignedCount++;
2202       }
2203     }
2204     return realignedCount;
2205   }
2206
2207   /**
2208    * Returns a map whose key is alignment column number (base 1), and whose
2209    * values are a map of sequence characters in that column.
2210    *
2211    * @param unaligned
2212    * @param aligned
2213    * @param unmapped
2214    * @return
2215    */
2216   static Map<Integer, Map<SequenceI, Character>> buildMappedColumnsMap(
2217           AlignmentI unaligned, AlignmentI aligned, List<SequenceI> unmapped)
2218   {
2219     /*
2220      * Map will hold, for each aligned column position, a map of
2221      * {unalignedSequence, sequenceCharacter} at that position.
2222      * TreeMap keeps the entries in ascending column order.
2223      */
2224     Map<Integer, Map<SequenceI, Character>> map = new TreeMap<Integer, Map<SequenceI, Character>>();
2225
2226     /*
2227      * r any sequences that have no mapping so can't be realigned
2228      */
2229     unmapped.addAll(unaligned.getSequences());
2230
2231     List<AlignedCodonFrame> mappings = aligned.getCodonFrames();
2232
2233     for (SequenceI seq : unaligned.getSequences())
2234     {
2235       for (AlignedCodonFrame mapping : mappings)
2236       {
2237         SequenceI fromSeq = mapping.findAlignedSequence(seq, aligned);
2238         if (fromSeq != null)
2239         {
2240           Mapping seqMap = mapping.getMappingBetween(fromSeq, seq);
2241           if (addMappedPositions(seq, fromSeq, seqMap, map))
2242           {
2243             unmapped.remove(seq);
2244           }
2245         }
2246       }
2247     }
2248     return map;
2249   }
2250
2251   /**
2252    * Helper method that adds to a map the mapped column positions of a sequence. <br>
2253    * For example if aaTT-Tg-gAAA is mapped to TTTAAA then the map should record
2254    * that columns 3,4,6,10,11,12 map to characters T,T,T,A,A,A of the mapped to
2255    * sequence.
2256    *
2257    * @param seq
2258    *          the sequence whose column positions we are recording
2259    * @param fromSeq
2260    *          a sequence that is mapped to the first sequence
2261    * @param seqMap
2262    *          the mapping from 'fromSeq' to 'seq'
2263    * @param map
2264    *          a map to add the column positions (in fromSeq) of the mapped
2265    *          positions of seq
2266    * @return
2267    */
2268   static boolean addMappedPositions(SequenceI seq, SequenceI fromSeq,
2269           Mapping seqMap, Map<Integer, Map<SequenceI, Character>> map)
2270   {
2271     if (seqMap == null)
2272     {
2273       return false;
2274     }
2275
2276     char[] fromChars = fromSeq.getSequence();
2277     int toStart = seq.getStart();
2278     char[] toChars = seq.getSequence();
2279
2280     /*
2281      * traverse [start, end, start, end...] ranges in fromSeq
2282      */
2283     for (int[] fromRange : seqMap.getMap().getFromRanges())
2284     {
2285       for (int i = 0; i < fromRange.length - 1; i += 2)
2286       {
2287         boolean forward = fromRange[i + 1] >= fromRange[i];
2288
2289         /*
2290          * find the range mapped to (sequence positions base 1)
2291          */
2292         int[] range = seqMap.locateMappedRange(fromRange[i],
2293                 fromRange[i + 1]);
2294         if (range == null)
2295         {
2296           System.err.println("Error in mapping " + seqMap + " from "
2297                   + fromSeq.getName());
2298           return false;
2299         }
2300         int fromCol = fromSeq.findIndex(fromRange[i]);
2301         int mappedCharPos = range[0];
2302
2303         /*
2304          * walk over the 'from' aligned sequence in forward or reverse
2305          * direction; when a non-gap is found, record the column position
2306          * of the next character of the mapped-to sequence; stop when all
2307          * the characters of the range have been counted
2308          */
2309         while (mappedCharPos <= range[1])
2310         {
2311           if (!Comparison.isGap(fromChars[fromCol - 1]))
2312           {
2313             /*
2314              * mapped from sequence has a character in this column
2315              * record the column position for the mapped to character
2316              */
2317             Map<SequenceI, Character> seqsMap = map.get(fromCol);
2318             if (seqsMap == null)
2319             {
2320               seqsMap = new HashMap<SequenceI, Character>();
2321               map.put(fromCol, seqsMap);
2322             }
2323             seqsMap.put(seq, toChars[mappedCharPos - toStart]);
2324             mappedCharPos++;
2325           }
2326           fromCol += (forward ? 1 : -1);
2327         }
2328       }
2329     }
2330     return true;
2331   }
2332
2333   // strictly temporary hack until proper criteria for aligning protein to cds
2334   // are in place; this is so Ensembl -> fetch xrefs Uniprot aligns the Uniprot
2335   public static boolean looksLikeEnsembl(AlignmentI alignment)
2336   {
2337     for (SequenceI seq : alignment.getSequences())
2338     {
2339       String name = seq.getName();
2340       if (!name.startsWith("ENSG") && !name.startsWith("ENST"))
2341       {
2342         return false;
2343       }
2344     }
2345     return true;
2346   }
2347 }