src/jalview/analysis/AlignmentUtils.java

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