JAL-1925 update source version in license
[jalview.git] / src / jalview / analysis / AlignmentUtils.java
index a4aeac7..f2262fb 100644 (file)
@@ -1,6 +1,6 @@
 /*
- * Jalview - A Sequence Alignment Editor and Viewer ($$Version-Rel$$)
- * Copyright (C) $$Year-Rel$$ The Jalview Authors
+ * Jalview - A Sequence Alignment Editor and Viewer (Version 2.9.0b2)
+ * Copyright (C) 2015 The Jalview Authors
  * 
  * This file is part of Jalview.
  * 
  */
 package jalview.analysis;
 
-import java.util.ArrayList;
-import java.util.Arrays;
-import java.util.Collection;
-import java.util.HashMap;
-import java.util.HashSet;
-import java.util.Iterator;
-import java.util.LinkedHashMap;
-import java.util.LinkedHashSet;
-import java.util.List;
-import java.util.Map;
-import java.util.Map.Entry;
-import java.util.Set;
-import java.util.TreeMap;
-
 import jalview.datamodel.AlignedCodon;
 import jalview.datamodel.AlignedCodonFrame;
 import jalview.datamodel.Alignment;
@@ -52,6 +38,20 @@ import jalview.util.DBRefUtils;
 import jalview.util.MapList;
 import jalview.util.MappingUtils;
 
+import java.util.ArrayList;
+import java.util.Arrays;
+import java.util.Collection;
+import java.util.HashMap;
+import java.util.HashSet;
+import java.util.Iterator;
+import java.util.LinkedHashMap;
+import java.util.LinkedHashSet;
+import java.util.List;
+import java.util.Map;
+import java.util.Map.Entry;
+import java.util.Set;
+import java.util.TreeMap;
+
 /**
  * grab bag of useful alignment manipulation operations Expect these to be
  * refactored elsewhere at some point.
@@ -77,18 +77,22 @@ public class AlignmentUtils
     for (SequenceI s : core.getSequences())
     {
       SequenceI newSeq = s.deriveSequence();
-      if (newSeq.getStart() > maxoffset
+      final int newSeqStart = newSeq.getStart() - 1;
+      if (newSeqStart > maxoffset
               && newSeq.getDatasetSequence().getStart() < s.getStart())
       {
-        maxoffset = newSeq.getStart();
+        maxoffset = newSeqStart;
       }
       sq.add(newSeq);
     }
     if (flankSize > -1)
     {
-      maxoffset = flankSize;
+      maxoffset = Math.min(maxoffset, flankSize);
     }
-    // now add offset to create a new expanded alignment
+
+    /*
+     * now add offset left and right to create an expanded alignment
+     */
     for (SequenceI s : sq)
     {
       SequenceI ds = s;
@@ -98,8 +102,8 @@ public class AlignmentUtils
       }
       int s_end = s.findPosition(s.getStart() + s.getLength());
       // find available flanking residues for sequence
-      int ustream_ds = s.getStart() - ds.getStart(), dstream_ds = ds
-              .getEnd() - s_end;
+      int ustream_ds = s.getStart() - ds.getStart();
+      int dstream_ds = ds.getEnd() - s_end;
 
       // build new flanked sequence
 
@@ -115,27 +119,27 @@ public class AlignmentUtils
           offset = maxoffset - flankSize;
           ustream_ds = flankSize;
         }
-        if (flankSize < dstream_ds)
+        if (flankSize <= dstream_ds)
         {
-          dstream_ds = flankSize;
+          dstream_ds = flankSize - 1;
         }
       }
+      // TODO use Character.toLowerCase to avoid creating String objects?
       char[] upstream = new String(ds.getSequence(s.getStart() - 1
               - ustream_ds, s.getStart() - 1)).toLowerCase().toCharArray();
-      char[] downstream = new String(ds.getSequence(s_end - 1, s_end + 1
+      char[] downstream = new String(ds.getSequence(s_end - 1, s_end
               + dstream_ds)).toLowerCase().toCharArray();
       char[] coreseq = s.getSequence();
       char[] nseq = new char[offset + upstream.length + downstream.length
               + coreseq.length];
       char c = core.getGapCharacter();
-      // TODO could lowercase the flanking regions
+
       int p = 0;
       for (; p < offset; p++)
       {
         nseq[p] = c;
       }
-      // s.setSequence(new String(upstream).toLowerCase()+new String(coreseq) +
-      // new String(downstream).toLowerCase());
+
       System.arraycopy(upstream, 0, nseq, p, upstream.length);
       System.arraycopy(coreseq, 0, nseq, p + upstream.length,
               coreseq.length);
@@ -153,6 +157,7 @@ public class AlignmentUtils
       {
         for (AlignmentAnnotation aa : s.getAnnotation())
         {
+          aa.adjustForAlignment(); // JAL-1712 fix
           newAl.addAnnotation(aa);
         }
       }
@@ -223,9 +228,8 @@ public class AlignmentUtils
    * @param cdnaAlignment
    * @return
    */
-  public static boolean mapProteinToCdna(
-          final AlignmentI proteinAlignment,
-          final AlignmentI cdnaAlignment)
+  public static boolean mapProteinAlignmentToCdna(
+          final AlignmentI proteinAlignment, final AlignmentI cdnaAlignment)
   {
     if (proteinAlignment == null || cdnaAlignment == null)
     {
@@ -271,7 +275,7 @@ public class AlignmentUtils
           final AlignmentI cdnaAlignment, Set<SequenceI> mappedDna,
           Set<SequenceI> mappedProtein, boolean xrefsOnly)
   {
-    boolean mappingPerformed = false;
+    boolean mappingExistsOrAdded = false;
     List<SequenceI> thisSeqs = proteinAlignment.getSequences();
     for (SequenceI aaSeq : thisSeqs)
     {
@@ -304,14 +308,18 @@ public class AlignmentUtils
         {
           continue;
         }
-        if (!mappingExists(proteinAlignment.getCodonFrames(),
+        if (mappingExists(proteinAlignment.getCodonFrames(),
                 aaSeq.getDatasetSequence(), cdnaSeq.getDatasetSequence()))
         {
-          MapList map = mapProteinToCdna(aaSeq, cdnaSeq);
+          mappingExistsOrAdded = true;
+        }
+        else
+        {
+          MapList map = mapProteinSequenceToCdna(aaSeq, cdnaSeq);
           if (map != null)
           {
             acf.addMap(cdnaSeq, aaSeq, map);
-            mappingPerformed = true;
+            mappingExistsOrAdded = true;
             proteinMapped = true;
             mappedDna.add(cdnaSeq);
             mappedProtein.add(aaSeq);
@@ -323,7 +331,7 @@ public class AlignmentUtils
         proteinAlignment.addCodonFrame(acf);
       }
     }
-    return mappingPerformed;
+    return mappingExistsOrAdded;
   }
 
   /**
@@ -356,7 +364,7 @@ public class AlignmentUtils
    * @param cdnaSeq
    * @return
    */
-  public static MapList mapProteinToCdna(SequenceI proteinSeq,
+  public static MapList mapProteinSequenceToCdna(SequenceI proteinSeq,
           SequenceI cdnaSeq)
   {
     /*
@@ -380,10 +388,10 @@ public class AlignmentUtils
      */
     final int mappedLength = 3 * aaSeqChars.length;
     int cdnaLength = cdnaSeqChars.length;
-    int cdnaStart = 1;
-    int cdnaEnd = cdnaLength;
-    final int proteinStart = 1;
-    final int proteinEnd = aaSeqChars.length;
+    int cdnaStart = cdnaSeq.getStart();
+    int cdnaEnd = cdnaSeq.getEnd();
+    final int proteinStart = proteinSeq.getStart();
+    final int proteinEnd = proteinSeq.getEnd();
 
     /*
      * If lengths don't match, try ignoring stop codon.
@@ -406,12 +414,13 @@ public class AlignmentUtils
     /*
      * If lengths still don't match, try ignoring start codon.
      */
+    int startOffset = 0;
     if (cdnaLength != mappedLength
             && cdnaLength > 2
             && String.valueOf(cdnaSeqChars, 0, 3).toUpperCase()
-                    .equals(
-                    ResidueProperties.START))
+                    .equals(ResidueProperties.START))
     {
+      startOffset += 3;
       cdnaStart += 3;
       cdnaLength -= 3;
     }
@@ -420,13 +429,12 @@ public class AlignmentUtils
     {
       return null;
     }
-    if (!translatesAs(cdnaSeqChars, cdnaStart - 1, aaSeqChars))
+    if (!translatesAs(cdnaSeqChars, startOffset, aaSeqChars))
     {
       return null;
     }
-    MapList map = new MapList(new int[]
-    { cdnaStart, cdnaEnd }, new int[]
-    { proteinStart, proteinEnd }, 3, 1);
+    MapList map = new MapList(new int[] { cdnaStart, cdnaEnd }, new int[] {
+        proteinStart, proteinEnd }, 3, 1);
     return map;
   }
 
@@ -443,23 +451,26 @@ public class AlignmentUtils
   protected static boolean translatesAs(char[] cdnaSeqChars, int cdnaStart,
           char[] aaSeqChars)
   {
+    if (cdnaSeqChars == null || aaSeqChars == null)
+    {
+      return false;
+    }
+
     int aaResidue = 0;
     for (int i = cdnaStart; i < cdnaSeqChars.length - 2
             && aaResidue < aaSeqChars.length; i += 3, aaResidue++)
     {
       String codon = String.valueOf(cdnaSeqChars, i, 3);
-      final String translated = ResidueProperties.codonTranslate(
-              codon);
+      final String translated = ResidueProperties.codonTranslate(codon);
       /*
-       * ? allow X in protein to match untranslatable in dna ?
+       * allow * in protein to match untranslatable in dna
        */
       final char aaRes = aaSeqChars[aaResidue];
-      if ((translated == null || "STOP".equals(translated)) && aaRes == 'X')
+      if ((translated == null || "STOP".equals(translated)) && aaRes == '*')
       {
         continue;
       }
-      if (translated == null
-              || !(aaRes == translated.charAt(0)))
+      if (translated == null || !(aaRes == translated.charAt(0)))
       {
         // debug
         // System.out.println(("Mismatch at " + i + "/" + aaResidue + ": "
@@ -490,7 +501,8 @@ public class AlignmentUtils
           boolean preserveUnmappedGaps)
   {
     /*
-     * Get any mappings from the source alignment to the target (dataset) sequence.
+     * Get any mappings from the source alignment to the target (dataset)
+     * sequence.
      */
     // TODO there may be one AlignedCodonFrame per dataset sequence, or one with
     // all mappings. Would it help to constrain this?
@@ -499,7 +511,7 @@ public class AlignmentUtils
     {
       return false;
     }
-  
+
     /*
      * Locate the aligned source sequence whose dataset sequence is mapped. We
      * just take the first match here (as we can't align cDNA like more than one
@@ -516,7 +528,7 @@ public class AlignmentUtils
         break;
       }
     }
-  
+
     if (alignFrom == null)
     {
       return false;
@@ -541,15 +553,15 @@ public class AlignmentUtils
    * @param preserveMappedGaps
    */
   public static void alignSequenceAs(SequenceI alignTo,
-          SequenceI alignFrom,
-          AlignedCodonFrame mapping, String myGap, char sourceGap,
-          boolean preserveMappedGaps, boolean preserveUnmappedGaps)
+          SequenceI alignFrom, AlignedCodonFrame mapping, String myGap,
+          char sourceGap, boolean preserveMappedGaps,
+          boolean preserveUnmappedGaps)
   {
     // TODO generalise to work for Protein-Protein, dna-dna, dna-protein
     final char[] thisSeq = alignTo.getSequence();
     final char[] thatAligned = alignFrom.getSequence();
     StringBuilder thisAligned = new StringBuilder(2 * thisSeq.length);
-  
+
     // aligned and dataset sequence positions, all base zero
     int thisSeqPos = 0;
     int sourceDsPos = 0;
@@ -561,6 +573,8 @@ public class AlignmentUtils
     /*
      * Traverse the aligned protein sequence.
      */
+    int fromOffset = alignFrom.getStart() - 1;
+    int toOffset = alignTo.getStart() - 1;
     int sourceGapMappedLength = 0;
     boolean inExon = false;
     for (char sourceChar : thatAligned)
@@ -577,7 +591,7 @@ public class AlignmentUtils
       sourceDsPos++;
       // Note mapping positions are base 1, our sequence positions base 0
       int[] mappedPos = mapping.getMappedRegion(alignTo, alignFrom,
-              sourceDsPos);
+              sourceDsPos + fromOffset);
       if (mappedPos == null)
       {
         /*
@@ -601,14 +615,15 @@ public class AlignmentUtils
        * But then 'align dna as protein' doesn't make much sense otherwise.
        */
       int intronLength = 0;
-      while (basesWritten < mappedCodonEnd && thisSeqPos < thisSeq.length)
+      while (basesWritten + toOffset < mappedCodonEnd
+              && thisSeqPos < thisSeq.length)
       {
         final char c = thisSeq[thisSeqPos++];
         if (c != myGapChar)
         {
           basesWritten++;
-
-          if (basesWritten < mappedCodonStart)
+          int sourcePosition = basesWritten + toOffset;
+          if (sourcePosition < mappedCodonStart)
           {
             /*
              * Found an unmapped (intron) base. First add in any preceding gaps
@@ -625,7 +640,7 @@ public class AlignmentUtils
           }
           else
           {
-            final boolean startOfCodon = basesWritten == mappedCodonStart;
+            final boolean startOfCodon = sourcePosition == mappedCodonStart;
             int gapsToAdd = calculateGapsToInsert(preserveMappedGaps,
                     preserveUnmappedGaps, sourceGapMappedLength, inExon,
                     trailingCopiedGap.length(), intronLength, startOfCodon);
@@ -686,8 +701,8 @@ public class AlignmentUtils
    */
   protected static int calculateGapsToInsert(boolean preserveMappedGaps,
           boolean preserveUnmappedGaps, int sourceGapMappedLength,
-          boolean inExon, int trailingGapLength,
-          int intronLength, final boolean startOfCodon)
+          boolean inExon, int trailingGapLength, int intronLength,
+          final boolean startOfCodon)
   {
     int gapsToAdd = 0;
     if (startOfCodon)
@@ -811,8 +826,8 @@ public class AlignmentUtils
       // mapping is from protein to nucleotide
       toDna = true;
       // should ideally get gap count ratio from mapping
-      gap = String.valueOf(new char[]
-      { gapCharacter, gapCharacter, gapCharacter });
+      gap = String.valueOf(new char[] { gapCharacter, gapCharacter,
+          gapCharacter });
     }
     else
     {
@@ -857,7 +872,7 @@ public class AlignmentUtils
           {
             mapping.markMappedRegion(seq, pos, sr);
           }
-          newseq.append(sr.toString());
+          newseq.append(sr.getCharacters());
           if (first)
           {
             first = false;
@@ -891,6 +906,9 @@ public class AlignmentUtils
    */
   public static int alignProteinAsDna(AlignmentI protein, AlignmentI dna)
   {
+    List<SequenceI> unmappedProtein = new ArrayList<SequenceI>();
+    unmappedProtein.addAll(protein.getSequences());
+
     Set<AlignedCodonFrame> mappings = protein.getCodonFrames();
 
     /*
@@ -911,10 +929,11 @@ public class AlignmentUtils
         {
           addCodonPositions(dnaSeq, prot, protein.getGapCharacter(),
                   seqMap, alignedCodons);
+          unmappedProtein.remove(prot);
         }
       }
     }
-    return alignProteinAs(protein, alignedCodons);
+    return alignProteinAs(protein, alignedCodons, unmappedProtein);
   }
 
   /**
@@ -925,10 +944,12 @@ public class AlignmentUtils
    * @param alignedCodons
    *          an ordered map of codon positions (columns), with sequence/peptide
    *          values present in each column
+   * @param unmappedProtein
    * @return
    */
   protected static int alignProteinAs(AlignmentI protein,
-          Map<AlignedCodon, Map<SequenceI, String>> alignedCodons)
+          Map<AlignedCodon, Map<SequenceI, String>> alignedCodons,
+          List<SequenceI> unmappedProtein)
   {
     /*
      * Prefill aligned sequences with gaps before inserting aligned protein
@@ -940,15 +961,18 @@ public class AlignmentUtils
     String allGaps = String.valueOf(gaps);
     for (SequenceI seq : protein.getSequences())
     {
-      seq.setSequence(allGaps);
+      if (!unmappedProtein.contains(seq))
+      {
+        seq.setSequence(allGaps);
+      }
     }
 
     int column = 0;
     for (AlignedCodon codon : alignedCodons.keySet())
     {
-      final Map<SequenceI, String> columnResidues = alignedCodons.get(codon);
-      for (Entry<SequenceI, String> entry : columnResidues
-              .entrySet())
+      final Map<SequenceI, String> columnResidues = alignedCodons
+              .get(codon);
+      for (Entry<SequenceI, String> entry : columnResidues.entrySet())
       {
         // place translated codon at its column position in sequence
         entry.getKey().getSequence()[column] = entry.getValue().charAt(0);
@@ -975,8 +999,7 @@ public class AlignmentUtils
    *          the map we are building up
    */
   static void addCodonPositions(SequenceI dna, SequenceI protein,
-          char gapChar,
-          Mapping seqMap,
+          char gapChar, Mapping seqMap,
           Map<AlignedCodon, Map<SequenceI, String>> alignedCodons)
   {
     Iterator<AlignedCodon> codons = seqMap.getCodonIterator(dna, gapChar);
@@ -1011,6 +1034,11 @@ public class AlignmentUtils
    */
   public static boolean isMappable(AlignmentI al1, AlignmentI al2)
   {
+    if (al1 == null || al2 == null)
+    {
+      return false;
+    }
+
     /*
      * Require one nucleotide and one protein
      */
@@ -1043,18 +1071,26 @@ public class AlignmentUtils
    * @param mappings
    * @return
    */
-  public static boolean isMappable(SequenceI dnaSeq, SequenceI proteinSeq,
-          Set<AlignedCodonFrame> mappings)
+  protected static boolean isMappable(SequenceI dnaSeq,
+          SequenceI proteinSeq, Set<AlignedCodonFrame> mappings)
   {
-    SequenceI dnaDs = dnaSeq.getDatasetSequence() == null ? dnaSeq : dnaSeq.getDatasetSequence();
+    if (dnaSeq == null || proteinSeq == null)
+    {
+      return false;
+    }
+
+    SequenceI dnaDs = dnaSeq.getDatasetSequence() == null ? dnaSeq : dnaSeq
+            .getDatasetSequence();
     SequenceI proteinDs = proteinSeq.getDatasetSequence() == null ? proteinSeq
             : proteinSeq.getDatasetSequence();
-    
+
     /*
      * Already mapped?
      */
-    for (AlignedCodonFrame mapping : mappings) {
-      if ( proteinDs == mapping.getAaForDnaSeq(dnaDs)) {
+    for (AlignedCodonFrame mapping : mappings)
+    {
+      if (proteinDs == mapping.getAaForDnaSeq(dnaDs))
+      {
         return true;
       }
     }
@@ -1063,7 +1099,7 @@ public class AlignmentUtils
      * Just try to make a mapping (it is not yet stored), test whether
      * successful.
      */
-    return mapProteinToCdna(proteinDs, dnaDs) != null;
+    return mapProteinSequenceToCdna(proteinDs, dnaDs) != null;
   }
 
   /**
@@ -1082,7 +1118,8 @@ public class AlignmentUtils
    *          the alignment to check for presence of annotations
    */
   public static void findAddableReferenceAnnotations(
-          List<SequenceI> sequenceScope, Map<String, String> labelForCalcId,
+          List<SequenceI> sequenceScope,
+          Map<String, String> labelForCalcId,
           final Map<SequenceI, List<AlignmentAnnotation>> candidates,
           AlignmentI al)
   {
@@ -1090,7 +1127,7 @@ public class AlignmentUtils
     {
       return;
     }
-  
+
     /*
      * For each sequence in scope, make a list of any annotations on the
      * underlying dataset sequence which are not already on the alignment.
@@ -1118,8 +1155,7 @@ public class AlignmentUtils
          * sequence.
          */
         final Iterable<AlignmentAnnotation> matchedAlignmentAnnotations = al
-                .findAnnotations(seq, dsann.getCalcId(),
-                        dsann.label);
+                .findAnnotations(seq, dsann.getCalcId(), dsann.label);
         if (!matchedAlignmentAnnotations.iterator().hasNext())
         {
           result.add(dsann);
@@ -1167,7 +1203,7 @@ public class AlignmentUtils
           endRes = selectionGroup.getEndRes();
         }
         copyAnn.restrict(startRes, endRes);
-  
+
         /*
          * Add to the sequence (sets copyAnn.datasetSequence), unless the
          * original annotation is already on the sequence.
@@ -1205,8 +1241,7 @@ public class AlignmentUtils
           Collection<String> types, List<SequenceI> forSequences,
           boolean anyType, boolean doShow)
   {
-    for (AlignmentAnnotation aa : al
-            .getAlignmentAnnotation())
+    for (AlignmentAnnotation aa : al.getAlignmentAnnotation())
     {
       if (anyType || types.contains(aa.label))
       {
@@ -1236,7 +1271,7 @@ public class AlignmentUtils
 
   /**
    * Returns true if seq1 has a cross-reference to seq2. Currently this assumes
-   * that sequence name is structured as Source|AccessId.
+   * that sequence name is structured as Source|AccessionId.
    * 
    * @param seq1
    * @param seq2
@@ -1281,7 +1316,7 @@ public class AlignmentUtils
   {
     Set<AlignedCodonFrame> newMappings = new LinkedHashSet<AlignedCodonFrame>();
     List<SequenceI> exonSequences = new ArrayList<SequenceI>();
-    
+
     for (SequenceI dnaSeq : dna)
     {
       final SequenceI ds = dnaSeq.getDatasetSequence();
@@ -1384,11 +1419,9 @@ public class AlignmentUtils
        * contiguous exons
        */
       List<int[]> exonRange = new ArrayList<int[]>();
-      exonRange.add(new int[]
-      { 1, newSequence.length() });
+      exonRange.add(new int[] { 1, newSequence.length() });
       MapList map = new MapList(exonRange, seqMapping.getMap()
-              .getToRanges(),
-              3, 1);
+              .getToRanges(), 3, 1);
       newMapping.addMap(exon.getDatasetSequence(), seqMapping.getTo(), map);
       MapList cdsToDnaMap = new MapList(dnaExonRanges, exonRange, 1, 1);
       newMapping.addMap(dnaSeq, exon.getDatasetSequence(), cdsToDnaMap);