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     for (AlignmentAnnotation aa : al.getAlignmentAnnotation())
1337     {
1338       if (anyType || types.contains(aa.label))
1339       {
1340         if ((aa.sequenceRef != null)
1341                 && (forSequences == null || forSequences
1342                         .contains(aa.sequenceRef)))
1343         {
1344           aa.visible = doShow;
1345         }
1346       }
1347     }
1348   }
1349
1350   /**
1351    * Returns true if either sequence has a cross-reference to the other
1352    *
1353    * @param seq1
1354    * @param seq2
1355    * @return
1356    */
1357   public static boolean haveCrossRef(SequenceI seq1, SequenceI seq2)
1358   {
1359     // Note: moved here from class CrossRef as the latter class has dependencies
1360     // not availability to the applet's classpath
1361     return hasCrossRef(seq1, seq2) || hasCrossRef(seq2, seq1);
1362   }
1363
1364   /**
1365    * Returns true if seq1 has a cross-reference to seq2. Currently this assumes
1366    * that sequence name is structured as Source|AccessionId.
1367    *
1368    * @param seq1
1369    * @param seq2
1370    * @return
1371    */
1372   public static boolean hasCrossRef(SequenceI seq1, SequenceI seq2)
1373   {
1374     if (seq1 == null || seq2 == null)
1375     {
1376       return false;
1377     }
1378     String name = seq2.getName();
1379     final DBRefEntry[] xrefs = seq1.getDBRefs();
1380     if (xrefs != null)
1381     {
1382       for (DBRefEntry xref : xrefs)
1383       {
1384         String xrefName = xref.getSource() + "|" + xref.getAccessionId();
1385         // case-insensitive test, consistent with DBRefEntry.equalRef()
1386         if (xrefName.equalsIgnoreCase(name))
1387         {
1388           return true;
1389         }
1390       }
1391     }
1392     return false;
1393   }
1394
1395   /**
1396    * Constructs an alignment consisting of the mapped (CDS) regions in the given
1397    * nucleotide sequences, and updates mappings to match. The CDS sequences are
1398    * added to the original alignment's dataset, which is shared by the new
1399    * alignment. Mappings from nucleotide to CDS, and from CDS to protein, are
1400    * added to the alignment dataset.
1401    *
1402    * @param dna
1403    *          aligned dna sequences
1404    * @param dataset
1405    *          - throws error if not given a dataset
1406    * @return an alignment whose sequences are the cds-only parts of the dna
1407    *         sequences (or null if no mappings are found)
1408    */
1409   public static AlignmentI makeCdsAlignment(SequenceI[] dna,
1410           AlignmentI dataset)
1411   {
1412     if (dataset.getDataset() != null)
1413     {
1414       throw new Error(
1415               "IMPLEMENTATION ERROR: dataset.getDataset() must be null!");
1416     }
1417     List<SequenceI> cdsSeqs = new ArrayList<SequenceI>();
1418     List<AlignedCodonFrame> mappings = dataset.getCodonFrames();
1419
1420
1421     /*
1422      * construct CDS sequences from the (cds-to-protein) mappings made earlier;
1423      * this makes it possible to model multiple products from dna (e.g. EMBL);
1424      * however it does mean we don't have the EMBL protein_id (a property on
1425      * the CDS features) in order to make the CDS sequence name :-(
1426      */
1427     for (SequenceI seq : dna)
1428     {
1429       SequenceI seqDss = seq.getDatasetSequence() == null ? seq : seq
1430               .getDatasetSequence();
1431       List<AlignedCodonFrame> seqMappings = MappingUtils
1432               .findMappingsForSequence(seq, mappings);
1433       for (AlignedCodonFrame mapping : seqMappings)
1434       {
1435         List<Mapping> mappingsFromSequence = mapping.getMappingsFromSequence(seq);
1436
1437         for (Mapping aMapping : mappingsFromSequence)
1438         {
1439           if (aMapping.getMap().getFromRatio() == 1)
1440           {
1441             /*
1442              * not a dna-to-protein mapping (likely dna-to-cds)
1443              */
1444             continue;
1445           }
1446
1447           /*
1448            * check for an existing CDS sequence i.e. a 3:1 mapping to
1449            * the dna mapping's product
1450            */
1451           SequenceI cdsSeq = null;
1452           // TODO better mappings collection data model so we can do
1453           // a table lookup instead of double loops to find mappings
1454           SequenceI proteinProduct = aMapping.getTo();
1455           for (AlignedCodonFrame acf : MappingUtils
1456                   .findMappingsForSequence(proteinProduct, mappings))
1457           {
1458             for (SequenceToSequenceMapping map : acf.getMappings())
1459             {
1460               if (map.getMapping().getMap().getFromRatio() == 3
1461                       && proteinProduct == map.getMapping().getTo()
1462                       && seqDss != map.getFromSeq())
1463               {
1464                 /*
1465                  * found a 3:1 mapping to the protein product which is not
1466                  * from the dna sequence...assume it is from the CDS sequence
1467                  * TODO mappings data model that brings together related
1468                  * dna-cds-protein mappings in one object
1469                  */
1470                 cdsSeq = map.getFromSeq();
1471               }
1472             }
1473           }
1474           if (cdsSeq != null)
1475           {
1476             /*
1477              * mappings are always to dataset sequences so create an aligned
1478              * sequence to own it; add the dataset sequence to the dataset
1479              */
1480             SequenceI derivedSequence = cdsSeq.deriveSequence();
1481             cdsSeqs.add(derivedSequence);
1482             if (!dataset.getSequences().contains(cdsSeq))
1483             {
1484               dataset.addSequence(cdsSeq);
1485             }
1486             continue;
1487           }
1488
1489           /*
1490            * didn't find mapped CDS sequence - construct it and add
1491            * its dataset sequence to the dataset
1492            */
1493           cdsSeq = makeCdsSequence(seq.getDatasetSequence(), aMapping);
1494           SequenceI cdsSeqDss = cdsSeq.createDatasetSequence();
1495           cdsSeqs.add(cdsSeq);
1496           if (!dataset.getSequences().contains(cdsSeqDss))
1497           {
1498             dataset.addSequence(cdsSeqDss);
1499           }
1500
1501           /*
1502            * add a mapping from CDS to the (unchanged) mapped to range
1503            */
1504           List<int[]> cdsRange = Collections.singletonList(new int[] { 1,
1505               cdsSeq.getLength() });
1506           MapList map = new MapList(cdsRange, aMapping.getMap()
1507                   .getToRanges(), aMapping.getMap().getFromRatio(),
1508                   aMapping.getMap().getToRatio());
1509           AlignedCodonFrame cdsToProteinMapping = new AlignedCodonFrame();
1510           cdsToProteinMapping.addMap(cdsSeq, proteinProduct, map);
1511
1512           /*
1513            * guard against duplicating the mapping if repeating this action
1514            */
1515           if (!mappings.contains(cdsToProteinMapping))
1516           {
1517             mappings.add(cdsToProteinMapping);
1518           }
1519
1520           /*
1521            * add another mapping from original 'from' range to CDS
1522            */
1523           AlignedCodonFrame dnaToProteinMapping = new AlignedCodonFrame();
1524           map = new MapList(aMapping.getMap().getFromRanges(), cdsRange, 1,
1525                   1);
1526           dnaToProteinMapping.addMap(seq.getDatasetSequence(), cdsSeq, map);
1527           if (!mappings.contains(dnaToProteinMapping))
1528           {
1529             mappings.add(dnaToProteinMapping);
1530           }
1531
1532
1533           /*
1534            * transfer any features on dna that overlap the CDS
1535            */
1536           transferFeatures(seq, cdsSeq, map, null, SequenceOntologyI.CDS);
1537         }
1538       }
1539     }
1540
1541     AlignmentI cds = new Alignment(cdsSeqs.toArray(new SequenceI[cdsSeqs
1542             .size()]));
1543     cds.setDataset((Alignment) dataset);
1544
1545     return cds;
1546   }
1547
1548   /**
1549    * Helper method that makes a CDS sequence as defined by the mappings from the
1550    * given sequence i.e. extracts the 'mapped from' ranges (which may be on
1551    * forward or reverse strand).
1552    *
1553    * @param seq
1554    * @param mapping
1555    * @return CDS sequence (as a dataset sequence)
1556    */
1557   static SequenceI makeCdsSequence(SequenceI seq, Mapping mapping)
1558   {
1559     char[] seqChars = seq.getSequence();
1560     List<int[]> fromRanges = mapping.getMap().getFromRanges();
1561     int cdsWidth = MappingUtils.getLength(fromRanges);
1562     char[] newSeqChars = new char[cdsWidth];
1563
1564     int newPos = 0;
1565     for (int[] range : fromRanges)
1566     {
1567       if (range[0] <= range[1])
1568       {
1569         // forward strand mapping - just copy the range
1570         int length = range[1] - range[0] + 1;
1571         System.arraycopy(seqChars, range[0] - 1, newSeqChars, newPos,
1572                 length);
1573         newPos += length;
1574       }
1575       else
1576       {
1577         // reverse strand mapping - copy and complement one by one
1578         for (int i = range[0]; i >= range[1]; i--)
1579         {
1580           newSeqChars[newPos++] = Dna.getComplement(seqChars[i - 1]);
1581         }
1582       }
1583     }
1584
1585     SequenceI newSeq = new Sequence(seq.getName() + "|"
1586             + mapping.getTo().getName(), newSeqChars, 1, newPos);
1587     return newSeq;
1588   }
1589
1590   /**
1591    * Transfers co-located features on 'fromSeq' to 'toSeq', adjusting the
1592    * feature start/end ranges, optionally omitting specified feature types.
1593    * Returns the number of features copied.
1594    *
1595    * @param fromSeq
1596    * @param toSeq
1597    * @param select
1598    *          if not null, only features of this type are copied (including
1599    *          subtypes in the Sequence Ontology)
1600    * @param mapping
1601    *          the mapping from 'fromSeq' to 'toSeq'
1602    * @param omitting
1603    */
1604   public static int transferFeatures(SequenceI fromSeq, SequenceI toSeq,
1605           MapList mapping, String select, String... omitting)
1606   {
1607     SequenceI copyTo = toSeq;
1608     while (copyTo.getDatasetSequence() != null)
1609     {
1610       copyTo = copyTo.getDatasetSequence();
1611     }
1612
1613     SequenceOntologyI so = SequenceOntologyFactory.getInstance();
1614     int count = 0;
1615     SequenceFeature[] sfs = fromSeq.getSequenceFeatures();
1616     if (sfs != null)
1617     {
1618       for (SequenceFeature sf : sfs)
1619       {
1620         String type = sf.getType();
1621         if (select != null && !so.isA(type, select))
1622         {
1623           continue;
1624         }
1625         boolean omit = false;
1626         for (String toOmit : omitting)
1627         {
1628           if (type.equals(toOmit))
1629           {
1630             omit = true;
1631           }
1632         }
1633         if (omit)
1634         {
1635           continue;
1636         }
1637
1638         /*
1639          * locate the mapped range - null if either start or end is
1640          * not mapped (no partial overlaps are calculated)
1641          */
1642         int start = sf.getBegin();
1643         int end = sf.getEnd();
1644         int[] mappedTo = mapping.locateInTo(start, end);
1645         /*
1646          * if whole exon range doesn't map, try interpreting it
1647          * as 5' or 3' exon overlapping the CDS range
1648          */
1649         if (mappedTo == null)
1650         {
1651           mappedTo = mapping.locateInTo(end, end);
1652           if (mappedTo != null)
1653           {
1654             /*
1655              * end of exon is in CDS range - 5' overlap
1656              * to a range from the start of the peptide
1657              */
1658             mappedTo[0] = 1;
1659           }
1660         }
1661         if (mappedTo == null)
1662         {
1663           mappedTo = mapping.locateInTo(start, start);
1664           if (mappedTo != null)
1665           {
1666             /*
1667              * start of exon is in CDS range - 3' overlap
1668              * to a range up to the end of the peptide
1669              */
1670             mappedTo[1] = toSeq.getLength();
1671           }
1672         }
1673         if (mappedTo != null)
1674         {
1675           SequenceFeature copy = new SequenceFeature(sf);
1676           copy.setBegin(Math.min(mappedTo[0], mappedTo[1]));
1677           copy.setEnd(Math.max(mappedTo[0], mappedTo[1]));
1678           copyTo.addSequenceFeature(copy);
1679           count++;
1680         }
1681       }
1682     }
1683     return count;
1684   }
1685
1686   /**
1687    * Returns a mapping from dna to protein by inspecting sequence features of
1688    * type "CDS" on the dna.
1689    *
1690    * @param dnaSeq
1691    * @param proteinSeq
1692    * @return
1693    */
1694   public static MapList mapCdsToProtein(SequenceI dnaSeq,
1695           SequenceI proteinSeq)
1696   {
1697     List<int[]> ranges = findCdsPositions(dnaSeq);
1698     int mappedDnaLength = MappingUtils.getLength(ranges);
1699
1700     int proteinLength = proteinSeq.getLength();
1701     int proteinStart = proteinSeq.getStart();
1702     int proteinEnd = proteinSeq.getEnd();
1703
1704     /*
1705      * incomplete start codon may mean X at start of peptide
1706      * we ignore both for mapping purposes
1707      */
1708     if (proteinSeq.getCharAt(0) == 'X')
1709     {
1710       // todo JAL-2022 support startPhase > 0
1711       proteinStart++;
1712       proteinLength--;
1713     }
1714     List<int[]> proteinRange = new ArrayList<int[]>();
1715
1716     /*
1717      * dna length should map to protein (or protein plus stop codon)
1718      */
1719     int codesForResidues = mappedDnaLength / 3;
1720     if (codesForResidues == (proteinLength + 1))
1721     {
1722       // assuming extra codon is for STOP and not in peptide
1723       codesForResidues--;
1724     }
1725     if (codesForResidues == proteinLength)
1726     {
1727       proteinRange.add(new int[] { proteinStart, proteinEnd });
1728       return new MapList(ranges, proteinRange, 3, 1);
1729     }
1730     return null;
1731   }
1732
1733   /**
1734    * Returns a list of CDS ranges found (as sequence positions base 1), i.e. of
1735    * start/end positions of sequence features of type "CDS" (or a sub-type of
1736    * CDS in the Sequence Ontology). The ranges are sorted into ascending start
1737    * position order, so this method is only valid for linear CDS in the same
1738    * sense as the protein product.
1739    *
1740    * @param dnaSeq
1741    * @return
1742    */
1743   public static List<int[]> findCdsPositions(SequenceI dnaSeq)
1744   {
1745     List<int[]> result = new ArrayList<int[]>();
1746     SequenceFeature[] sfs = dnaSeq.getSequenceFeatures();
1747     if (sfs == null)
1748     {
1749       return result;
1750     }
1751
1752     SequenceOntologyI so = SequenceOntologyFactory.getInstance();
1753     int startPhase = 0;
1754
1755     for (SequenceFeature sf : sfs)
1756     {
1757       /*
1758        * process a CDS feature (or a sub-type of CDS)
1759        */
1760       if (so.isA(sf.getType(), SequenceOntologyI.CDS))
1761       {
1762         int phase = 0;
1763         try
1764         {
1765           phase = Integer.parseInt(sf.getPhase());
1766         } catch (NumberFormatException e)
1767         {
1768           // ignore
1769         }
1770         /*
1771          * phase > 0 on first codon means 5' incomplete - skip to the start
1772          * of the next codon; example ENST00000496384
1773          */
1774         int begin = sf.getBegin();
1775         int end = sf.getEnd();
1776         if (result.isEmpty())
1777         {
1778           begin += phase;
1779           if (begin > end)
1780           {
1781             // shouldn't happen!
1782             System.err
1783                     .println("Error: start phase extends beyond start CDS in "
1784                             + dnaSeq.getName());
1785           }
1786         }
1787         result.add(new int[] { begin, end });
1788       }
1789     }
1790
1791     /*
1792      * remove 'startPhase' positions (usually 0) from the first range
1793      * so we begin at the start of a complete codon
1794      */
1795     if (!result.isEmpty())
1796     {
1797       // TODO JAL-2022 correctly model start phase > 0
1798       result.get(0)[0] += startPhase;
1799     }
1800
1801     /*
1802      * Finally sort ranges by start position. This avoids a dependency on
1803      * keeping features in order on the sequence (if they are in order anyway,
1804      * the sort will have almost no work to do). The implicit assumption is CDS
1805      * ranges are assembled in order. Other cases should not use this method,
1806      * but instead construct an explicit mapping for CDS (e.g. EMBL parsing).
1807      */
1808     Collections.sort(result, new Comparator<int[]>()
1809     {
1810       @Override
1811       public int compare(int[] o1, int[] o2)
1812       {
1813         return Integer.compare(o1[0], o2[0]);
1814       }
1815     });
1816     return result;
1817   }
1818
1819   /**
1820    * Maps exon features from dna to protein, and computes variants in peptide
1821    * product generated by variants in dna, and adds them as sequence_variant
1822    * features on the protein sequence. Returns the number of variant features
1823    * added.
1824    *
1825    * @param dnaSeq
1826    * @param peptide
1827    * @param dnaToProtein
1828    */
1829   public static int computeProteinFeatures(SequenceI dnaSeq,
1830           SequenceI peptide, MapList dnaToProtein)
1831   {
1832     while (dnaSeq.getDatasetSequence() != null)
1833     {
1834       dnaSeq = dnaSeq.getDatasetSequence();
1835     }
1836     while (peptide.getDatasetSequence() != null)
1837     {
1838       peptide = peptide.getDatasetSequence();
1839     }
1840
1841     transferFeatures(dnaSeq, peptide, dnaToProtein, SequenceOntologyI.EXON);
1842
1843     /*
1844      * compute protein variants from dna variants and codon mappings;
1845      * NB - alternatively we could retrieve this using the REST service e.g.
1846      * http://rest.ensembl.org/overlap/translation
1847      * /ENSP00000288602?feature=transcript_variation;content-type=text/xml
1848      * which would be a bit slower but possibly more reliable
1849      */
1850
1851     /*
1852      * build a map with codon variations for each potentially varying peptide
1853      */
1854     LinkedHashMap<Integer, List<DnaVariant>[]> variants = buildDnaVariantsMap(
1855             dnaSeq, dnaToProtein);
1856
1857     /*
1858      * scan codon variations, compute peptide variants and add to peptide sequence
1859      */
1860     int count = 0;
1861     for (Entry<Integer, List<DnaVariant>[]> variant : variants.entrySet())
1862     {
1863       int peptidePos = variant.getKey();
1864       List<DnaVariant>[] codonVariants = variant.getValue();
1865       count += computePeptideVariants(peptide, peptidePos, codonVariants);
1866     }
1867
1868     /*
1869      * sort to get sequence features in start position order
1870      * - would be better to store in Sequence as a TreeSet or NCList?
1871      */
1872     if (peptide.getSequenceFeatures() != null)
1873     {
1874       Arrays.sort(peptide.getSequenceFeatures(),
1875               new Comparator<SequenceFeature>()
1876               {
1877                 @Override
1878                 public int compare(SequenceFeature o1, SequenceFeature o2)
1879                 {
1880                   int c = Integer.compare(o1.getBegin(), o2.getBegin());
1881                   return c == 0 ? Integer.compare(o1.getEnd(), o2.getEnd())
1882                           : c;
1883                 }
1884               });
1885     }
1886     return count;
1887   }
1888
1889   /**
1890    * Computes non-synonymous peptide variants from codon variants and adds them
1891    * as sequence_variant features on the protein sequence (one feature per
1892    * allele variant). Selected attributes (variant id, clinical significance)
1893    * are copied over to the new features.
1894    *
1895    * @param peptide
1896    *          the protein sequence
1897    * @param peptidePos
1898    *          the position to compute peptide variants for
1899    * @param codonVariants
1900    *          a list of dna variants per codon position
1901    * @return the number of features added
1902    */
1903   static int computePeptideVariants(SequenceI peptide, int peptidePos,
1904           List<DnaVariant>[] codonVariants)
1905   {
1906     String residue = String.valueOf(peptide.getCharAt(peptidePos - 1));
1907     int count = 0;
1908     String base1 = codonVariants[0].get(0).base;
1909     String base2 = codonVariants[1].get(0).base;
1910     String base3 = codonVariants[2].get(0).base;
1911
1912     /*
1913      * variants in first codon base
1914      */
1915     for (DnaVariant var : codonVariants[0])
1916     {
1917       if (var.variant != null)
1918       {
1919         String alleles = (String) var.variant.getValue("alleles");
1920         if (alleles != null)
1921         {
1922           for (String base : alleles.split(","))
1923           {
1924             String codon = base + base2 + base3;
1925             if (addPeptideVariant(peptide, peptidePos, residue, var, codon))
1926             {
1927               count++;
1928             }
1929           }
1930         }
1931       }
1932     }
1933
1934     /*
1935      * variants in second codon base
1936      */
1937     for (DnaVariant var : codonVariants[1])
1938     {
1939       if (var.variant != null)
1940       {
1941         String alleles = (String) var.variant.getValue("alleles");
1942         if (alleles != null)
1943         {
1944           for (String base : alleles.split(","))
1945           {
1946             String codon = base1 + base + base3;
1947             if (addPeptideVariant(peptide, peptidePos, residue, var, codon))
1948             {
1949               count++;
1950             }
1951           }
1952         }
1953       }
1954     }
1955
1956     /*
1957      * variants in third codon base
1958      */
1959     for (DnaVariant var : codonVariants[2])
1960     {
1961       if (var.variant != null)
1962       {
1963         String alleles = (String) var.variant.getValue("alleles");
1964         if (alleles != null)
1965         {
1966           for (String base : alleles.split(","))
1967           {
1968             String codon = base1 + base2 + base;
1969             if (addPeptideVariant(peptide, peptidePos, residue, var, codon))
1970             {
1971               count++;
1972             }
1973           }
1974         }
1975       }
1976     }
1977
1978     return count;
1979   }
1980
1981   /**
1982    * Helper method that adds a peptide variant feature, provided the given codon
1983    * translates to a value different to the current residue (is a non-synonymous
1984    * variant). ID and clinical_significance attributes of the dna variant (if
1985    * present) are copied to the new feature.
1986    *
1987    * @param peptide
1988    * @param peptidePos
1989    * @param residue
1990    * @param var
1991    * @param codon
1992    * @return true if a feature was added, else false
1993    */
1994   static boolean addPeptideVariant(SequenceI peptide, int peptidePos,
1995           String residue, DnaVariant var, String codon)
1996   {
1997     /*
1998      * get peptide translation of codon e.g. GAT -> D
1999      * note that variants which are not single alleles,
2000      * e.g. multibase variants or HGMD_MUTATION etc
2001      * are currently ignored here
2002      */
2003     String trans = codon.contains("-") ? "-"
2004             : (codon.length() > 3 ? null : ResidueProperties
2005                     .codonTranslate(codon));
2006     if (trans != null && !trans.equals(residue))
2007     {
2008       String residue3Char = StringUtils
2009               .toSentenceCase(ResidueProperties.aa2Triplet.get(residue));
2010       String trans3Char = StringUtils
2011               .toSentenceCase(ResidueProperties.aa2Triplet.get(trans));
2012       String desc = "p." + residue3Char + peptidePos + trans3Char;
2013       // set score to 0f so 'graduated colour' option is offered! JAL-2060
2014       SequenceFeature sf = new SequenceFeature(
2015               SequenceOntologyI.SEQUENCE_VARIANT, desc, peptidePos,
2016               peptidePos, 0f, "Jalview");
2017       StringBuilder attributes = new StringBuilder(32);
2018       String id = (String) var.variant.getValue(ID);
2019       if (id != null)
2020       {
2021         if (id.startsWith(SEQUENCE_VARIANT))
2022         {
2023           id = id.substring(SEQUENCE_VARIANT.length());
2024         }
2025         sf.setValue(ID, id);
2026         attributes.append(ID).append("=").append(id);
2027         // TODO handle other species variants
2028         StringBuilder link = new StringBuilder(32);
2029         try
2030         {
2031           link.append(desc).append(" ").append(id)
2032                   .append("|http://www.ensembl.org/Homo_sapiens/Variation/Summary?v=")
2033                   .append(URLEncoder.encode(id, "UTF-8"));
2034           sf.addLink(link.toString());
2035         } catch (UnsupportedEncodingException e)
2036         {
2037           // as if
2038         }
2039       }
2040       String clinSig = (String) var.variant
2041               .getValue(CLINICAL_SIGNIFICANCE);
2042       if (clinSig != null)
2043       {
2044         sf.setValue(CLINICAL_SIGNIFICANCE, clinSig);
2045         attributes.append(";").append(CLINICAL_SIGNIFICANCE).append("=")
2046                 .append(clinSig);
2047       }
2048       peptide.addSequenceFeature(sf);
2049       if (attributes.length() > 0)
2050       {
2051         sf.setAttributes(attributes.toString());
2052       }
2053       return true;
2054     }
2055     return false;
2056   }
2057
2058   /**
2059    * Builds a map whose key is position in the protein sequence, and value is a
2060    * list of the base and all variants for each corresponding codon position
2061    *
2062    * @param dnaSeq
2063    * @param dnaToProtein
2064    * @return
2065    */
2066   static LinkedHashMap<Integer, List<DnaVariant>[]> buildDnaVariantsMap(
2067           SequenceI dnaSeq, MapList dnaToProtein)
2068   {
2069     /*
2070      * map from peptide position to all variants of the codon which codes for it
2071      * LinkedHashMap ensures we keep the peptide features in sequence order
2072      */
2073     LinkedHashMap<Integer, List<DnaVariant>[]> variants = new LinkedHashMap<Integer, List<DnaVariant>[]>();
2074     SequenceOntologyI so = SequenceOntologyFactory.getInstance();
2075
2076     SequenceFeature[] dnaFeatures = dnaSeq.getSequenceFeatures();
2077     if (dnaFeatures == null)
2078     {
2079       return variants;
2080     }
2081
2082     int dnaStart = dnaSeq.getStart();
2083     int[] lastCodon = null;
2084     int lastPeptidePostion = 0;
2085
2086     /*
2087      * build a map of codon variations for peptides
2088      */
2089     for (SequenceFeature sf : dnaFeatures)
2090     {
2091       int dnaCol = sf.getBegin();
2092       if (dnaCol != sf.getEnd())
2093       {
2094         // not handling multi-locus variant features
2095         continue;
2096       }
2097       if (so.isA(sf.getType(), SequenceOntologyI.SEQUENCE_VARIANT))
2098       {
2099         int[] mapsTo = dnaToProtein.locateInTo(dnaCol, dnaCol);
2100         if (mapsTo == null)
2101         {
2102           // feature doesn't lie within coding region
2103           continue;
2104         }
2105         int peptidePosition = mapsTo[0];
2106         List<DnaVariant>[] codonVariants = variants.get(peptidePosition);
2107         if (codonVariants == null)
2108         {
2109           codonVariants = new ArrayList[3];
2110           codonVariants[0] = new ArrayList<DnaVariant>();
2111           codonVariants[1] = new ArrayList<DnaVariant>();
2112           codonVariants[2] = new ArrayList<DnaVariant>();
2113           variants.put(peptidePosition, codonVariants);
2114         }
2115
2116         /*
2117          * extract dna variants to a string array
2118          */
2119         String alls = (String) sf.getValue("alleles");
2120         if (alls == null)
2121         {
2122           continue;
2123         }
2124         String[] alleles = alls.toUpperCase().split(",");
2125         int i = 0;
2126         for (String allele : alleles)
2127         {
2128           alleles[i++] = allele.trim(); // lose any space characters "A, G"
2129         }
2130
2131         /*
2132          * get this peptide's codon positions e.g. [3, 4, 5] or [4, 7, 10]
2133          */
2134         int[] codon = peptidePosition == lastPeptidePostion ? lastCodon
2135                 : MappingUtils.flattenRanges(dnaToProtein.locateInFrom(
2136                         peptidePosition, peptidePosition));
2137         lastPeptidePostion = peptidePosition;
2138         lastCodon = codon;
2139
2140         /*
2141          * save nucleotide (and any variant) for each codon position
2142          */
2143         for (int codonPos = 0; codonPos < 3; codonPos++)
2144         {
2145           String nucleotide = String.valueOf(
2146                   dnaSeq.getCharAt(codon[codonPos] - dnaStart))
2147                   .toUpperCase();
2148           List<DnaVariant> codonVariant = codonVariants[codonPos];
2149           if (codon[codonPos] == dnaCol)
2150           {
2151             if (!codonVariant.isEmpty()
2152                     && codonVariant.get(0).variant == null)
2153             {
2154               /*
2155                * already recorded base value, add this variant
2156                */
2157               codonVariant.get(0).variant = sf;
2158             }
2159             else
2160             {
2161               /*
2162                * add variant with base value
2163                */
2164               codonVariant.add(new DnaVariant(nucleotide, sf));
2165             }
2166           }
2167           else if (codonVariant.isEmpty())
2168           {
2169             /*
2170              * record (possibly non-varying) base value
2171              */
2172             codonVariant.add(new DnaVariant(nucleotide));
2173           }
2174         }
2175       }
2176     }
2177     return variants;
2178   }
2179
2180   /**
2181    * Makes an alignment with a copy of the given sequences, adding in any
2182    * non-redundant sequences which are mapped to by the cross-referenced
2183    * sequences.
2184    *
2185    * @param seqs
2186    * @param xrefs
2187    * @return
2188    */
2189   public static AlignmentI makeCopyAlignment(SequenceI[] seqs,
2190           SequenceI[] xrefs)
2191   {
2192     AlignmentI copy = new Alignment(new Alignment(seqs));
2193
2194     /*
2195      * add mappings between sequences to the new alignment
2196      */
2197     AlignedCodonFrame mappings = new AlignedCodonFrame();
2198     copy.addCodonFrame(mappings);
2199     for (int i = 0; i < copy.getHeight(); i++)
2200     {
2201       SequenceI from = seqs[i];
2202       SequenceI to = copy.getSequenceAt(i);
2203       if (to.getDatasetSequence() != null)
2204       {
2205         to = to.getDatasetSequence();
2206       }
2207       int start = from.getStart();
2208       int end = from.getEnd();
2209       MapList map = new MapList(new int[] { start, end }, new int[] {
2210           start, end }, 1, 1);
2211       mappings.addMap(to, from, map);
2212     }
2213
2214     SequenceIdMatcher matcher = new SequenceIdMatcher(seqs);
2215     if (xrefs != null)
2216     {
2217       for (SequenceI xref : xrefs)
2218       {
2219         DBRefEntry[] dbrefs = xref.getDBRefs();
2220         if (dbrefs != null)
2221         {
2222           for (DBRefEntry dbref : dbrefs)
2223           {
2224             if (dbref.getMap() == null || dbref.getMap().getTo() == null)
2225             {
2226               continue;
2227             }
2228             SequenceI mappedTo = dbref.getMap().getTo();
2229             SequenceI match = matcher.findIdMatch(mappedTo);
2230             if (match == null)
2231             {
2232               matcher.add(mappedTo);
2233               copy.addSequence(mappedTo);
2234             }
2235           }
2236         }
2237       }
2238     }
2239     return copy;
2240   }
2241
2242   /**
2243    * Try to align sequences in 'unaligned' to match the alignment of their
2244    * mapped regions in 'aligned'. For example, could use this to align CDS
2245    * sequences which are mapped to their parent cDNA sequences.
2246    *
2247    * This method handles 1:1 mappings (dna-to-dna or protein-to-protein). For
2248    * dna-to-protein or protein-to-dna use alternative methods.
2249    *
2250    * @param unaligned
2251    *          sequences to be aligned
2252    * @param aligned
2253    *          holds aligned sequences and their mappings
2254    * @return
2255    */
2256   public static int alignAs(AlignmentI unaligned, AlignmentI aligned)
2257   {
2258     List<SequenceI> unmapped = new ArrayList<SequenceI>();
2259     Map<Integer, Map<SequenceI, Character>> columnMap = buildMappedColumnsMap(
2260             unaligned, aligned, unmapped);
2261     int width = columnMap.size();
2262     char gap = unaligned.getGapCharacter();
2263     int realignedCount = 0;
2264
2265     for (SequenceI seq : unaligned.getSequences())
2266     {
2267       if (!unmapped.contains(seq))
2268       {
2269         char[] newSeq = new char[width];
2270         Arrays.fill(newSeq, gap);
2271         int newCol = 0;
2272         int lastCol = 0;
2273
2274         /*
2275          * traverse the map to find columns populated
2276          * by our sequence
2277          */
2278         for (Integer column : columnMap.keySet())
2279         {
2280           Character c = columnMap.get(column).get(seq);
2281           if (c != null)
2282           {
2283             /*
2284              * sequence has a character at this position
2285              *
2286              */
2287             newSeq[newCol] = c;
2288             lastCol = newCol;
2289           }
2290           newCol++;
2291         }
2292
2293         /*
2294          * trim trailing gaps
2295          */
2296         if (lastCol < width)
2297         {
2298           char[] tmp = new char[lastCol + 1];
2299           System.arraycopy(newSeq, 0, tmp, 0, lastCol + 1);
2300           newSeq = tmp;
2301         }
2302         seq.setSequence(String.valueOf(newSeq));
2303         realignedCount++;
2304       }
2305     }
2306     return realignedCount;
2307   }
2308
2309   /**
2310    * Returns a map whose key is alignment column number (base 1), and whose
2311    * values are a map of sequence characters in that column.
2312    *
2313    * @param unaligned
2314    * @param aligned
2315    * @param unmapped
2316    * @return
2317    */
2318   static Map<Integer, Map<SequenceI, Character>> buildMappedColumnsMap(
2319           AlignmentI unaligned, AlignmentI aligned, List<SequenceI> unmapped)
2320   {
2321     /*
2322      * Map will hold, for each aligned column position, a map of
2323      * {unalignedSequence, sequenceCharacter} at that position.
2324      * TreeMap keeps the entries in ascending column order.
2325      */
2326     Map<Integer, Map<SequenceI, Character>> map = new TreeMap<Integer, Map<SequenceI, Character>>();
2327
2328     /*
2329      * r any sequences that have no mapping so can't be realigned
2330      */
2331     unmapped.addAll(unaligned.getSequences());
2332
2333     List<AlignedCodonFrame> mappings = aligned.getCodonFrames();
2334
2335     for (SequenceI seq : unaligned.getSequences())
2336     {
2337       for (AlignedCodonFrame mapping : mappings)
2338       {
2339         SequenceI fromSeq = mapping.findAlignedSequence(seq, aligned);
2340         if (fromSeq != null)
2341         {
2342           Mapping seqMap = mapping.getMappingBetween(fromSeq, seq);
2343           if (addMappedPositions(seq, fromSeq, seqMap, map))
2344           {
2345             unmapped.remove(seq);
2346           }
2347         }
2348       }
2349     }
2350     return map;
2351   }
2352
2353   /**
2354    * Helper method that adds to a map the mapped column positions of a sequence. <br>
2355    * For example if aaTT-Tg-gAAA is mapped to TTTAAA then the map should record
2356    * that columns 3,4,6,10,11,12 map to characters T,T,T,A,A,A of the mapped to
2357    * sequence.
2358    *
2359    * @param seq
2360    *          the sequence whose column positions we are recording
2361    * @param fromSeq
2362    *          a sequence that is mapped to the first sequence
2363    * @param seqMap
2364    *          the mapping from 'fromSeq' to 'seq'
2365    * @param map
2366    *          a map to add the column positions (in fromSeq) of the mapped
2367    *          positions of seq
2368    * @return
2369    */
2370   static boolean addMappedPositions(SequenceI seq, SequenceI fromSeq,
2371           Mapping seqMap, Map<Integer, Map<SequenceI, Character>> map)
2372   {
2373     if (seqMap == null)
2374     {
2375       return false;
2376     }
2377
2378     char[] fromChars = fromSeq.getSequence();
2379     int toStart = seq.getStart();
2380     char[] toChars = seq.getSequence();
2381
2382     /*
2383      * traverse [start, end, start, end...] ranges in fromSeq
2384      */
2385     for (int[] fromRange : seqMap.getMap().getFromRanges())
2386     {
2387       for (int i = 0; i < fromRange.length - 1; i += 2)
2388       {
2389         boolean forward = fromRange[i + 1] >= fromRange[i];
2390
2391         /*
2392          * find the range mapped to (sequence positions base 1)
2393          */
2394         int[] range = seqMap.locateMappedRange(fromRange[i],
2395                 fromRange[i + 1]);
2396         if (range == null)
2397         {
2398           System.err.println("Error in mapping " + seqMap + " from "
2399                   + fromSeq.getName());
2400           return false;
2401         }
2402         int fromCol = fromSeq.findIndex(fromRange[i]);
2403         int mappedCharPos = range[0];
2404
2405         /*
2406          * walk over the 'from' aligned sequence in forward or reverse
2407          * direction; when a non-gap is found, record the column position
2408          * of the next character of the mapped-to sequence; stop when all
2409          * the characters of the range have been counted
2410          */
2411         while (mappedCharPos <= range[1])
2412         {
2413           if (!Comparison.isGap(fromChars[fromCol - 1]))
2414           {
2415             /*
2416              * mapped from sequence has a character in this column
2417              * record the column position for the mapped to character
2418              */
2419             Map<SequenceI, Character> seqsMap = map.get(fromCol);
2420             if (seqsMap == null)
2421             {
2422               seqsMap = new HashMap<SequenceI, Character>();
2423               map.put(fromCol, seqsMap);
2424             }
2425             seqsMap.put(seq, toChars[mappedCharPos - toStart]);
2426             mappedCharPos++;
2427           }
2428           fromCol += (forward ? 1 : -1);
2429         }
2430       }
2431     }
2432     return true;
2433   }
2434
2435   // strictly temporary hack until proper criteria for aligning protein to cds
2436   // are in place; this is so Ensembl -> fetch xrefs Uniprot aligns the Uniprot
2437   public static boolean looksLikeEnsembl(AlignmentI alignment)
2438   {
2439     for (SequenceI seq : alignment.getSequences())
2440     {
2441       String name = seq.getName();
2442       if (!name.startsWith("ENSG") && !name.startsWith("ENST"))
2443       {
2444         return false;
2445       }
2446     }
2447     return true;
2448   }
2449 }