JAL-3251 SequenceToSequenceMapping now SequenceMapping, MappingType++
[jalview.git] / src / jalview / analysis / AlignmentUtils.java
index 9278361..f1084a7 100644 (file)
@@ -24,7 +24,6 @@ import static jalview.io.gff.GffConstants.CLINICAL_SIGNIFICANCE;
 
 import jalview.datamodel.AlignedCodon;
 import jalview.datamodel.AlignedCodonFrame;
-import jalview.datamodel.AlignedCodonFrame.SequenceToSequenceMapping;
 import jalview.datamodel.Alignment;
 import jalview.datamodel.AlignmentAnnotation;
 import jalview.datamodel.AlignmentI;
@@ -36,6 +35,7 @@ import jalview.datamodel.Sequence;
 import jalview.datamodel.SequenceFeature;
 import jalview.datamodel.SequenceGroup;
 import jalview.datamodel.SequenceI;
+import jalview.datamodel.SequenceMapping;
 import jalview.datamodel.features.SequenceFeatures;
 import jalview.io.gff.Gff3Helper;
 import jalview.io.gff.SequenceOntologyI;
@@ -74,12 +74,15 @@ import java.util.TreeMap;
  */
 public class AlignmentUtils
 {
-
   private static final int CODON_LENGTH = 3;
 
   private static final String SEQUENCE_VARIANT = "sequence_variant:";
 
-  private static final String ID = "ID";
+  /*
+   * the 'id' attribute is provided for variant features fetched from
+   * Ensembl using its REST service with JSON format 
+   */
+  public static final String VARIANT_ID = "id";
 
   /**
    * A data model to hold the 'normal' base value at a position, and an optional
@@ -465,7 +468,7 @@ public class AlignmentUtils
     {
       String lastCodon = String.valueOf(cdnaSeqChars,
               cdnaLength - CODON_LENGTH, CODON_LENGTH).toUpperCase();
-      for (String stop : ResidueProperties.STOP)
+      for (String stop : ResidueProperties.STOP_CODONS)
       {
         if (lastCodon.equals(stop))
         {
@@ -536,7 +539,8 @@ public class AlignmentUtils
        * allow * in protein to match untranslatable in dna
        */
       final char aaRes = aaSeqChars[aaPos];
-      if ((translated == null || "STOP".equals(translated)) && aaRes == '*')
+      if ((translated == null || ResidueProperties.STOP.equals(translated))
+              && aaRes == '*')
       {
         continue;
       }
@@ -568,7 +572,8 @@ public class AlignmentUtils
     if (dnaPos == cdnaSeqChars.length - CODON_LENGTH)
     {
       String codon = String.valueOf(cdnaSeqChars, dnaPos, CODON_LENGTH);
-      if ("STOP".equals(ResidueProperties.codonTranslate(codon)))
+      if (ResidueProperties.STOP
+              .equals(ResidueProperties.codonTranslate(codon)))
       {
         return true;
       }
@@ -1102,7 +1107,7 @@ public class AlignmentUtils
         SequenceI prot = mapping.findAlignedSequence(dnaSeq, protein);
         if (prot != null)
         {
-          Mapping seqMap = mapping.getMappingForSequence(dnaSeq);
+          Mapping seqMap = mapping.getMappingForSequence(dnaSeq, false);
           addCodonPositions(dnaSeq, prot, protein.getGapCharacter(), seqMap,
                   alignedCodons);
           unmappedProtein.remove(prot);
@@ -1886,7 +1891,7 @@ public class AlignmentUtils
    * @param seqMappings
    *          the set of mappings involving dnaSeq
    * @param aMapping
-   *          an initial candidate from seqMappings
+   *          a transcript-to-peptide mapping
    * @return
    */
   static SequenceI findCdsForProtein(List<AlignedCodonFrame> mappings,
@@ -1911,7 +1916,15 @@ public class AlignmentUtils
     if (mappedFromLength == dnaLength
             || mappedFromLength == dnaLength - CODON_LENGTH)
     {
-      return seqDss;
+      /*
+       * if sequence has CDS features, this is a transcript with no UTR
+       * - do not take this as the CDS sequence! (JAL-2789)
+       */
+      if (seqDss.getFeatures().getFeaturesByOntology(SequenceOntologyI.CDS)
+              .isEmpty())
+      {
+        return seqDss;
+      }
     }
 
     /*
@@ -1922,7 +1935,7 @@ public class AlignmentUtils
             .findMappingsForSequence(proteinProduct, mappings);
     for (AlignedCodonFrame acf : mappingsToPeptide)
     {
-      for (SequenceToSequenceMapping map : acf.getMappings())
+      for (SequenceMapping map : acf.getMappings())
       {
         Mapping mapping = map.getMapping();
         if (mapping != aMapping
@@ -1936,10 +1949,12 @@ public class AlignmentUtils
           {
             /*
             * found a 3:1 mapping to the protein product which covers
-            * the whole dna sequence i.e. is from CDS; finally check it
-            * is from the dna start sequence
+            * the whole dna sequence i.e. is from CDS; finally check the CDS
+            * is mapped from the given dna start sequence
             */
             SequenceI cdsSeq = map.getFromSeq();
+            // todo this test is weak if seqMappings contains multiple mappings;
+            // we get away with it if transcript:cds relationship is 1:1
             List<AlignedCodonFrame> dnaToCdsMaps = MappingUtils
                     .findMappingsForSequence(cdsSeq, seqMappings);
             if (!dnaToCdsMaps.isEmpty())
@@ -2054,9 +2069,11 @@ public class AlignmentUtils
   protected static List<DBRefEntry> propagateDBRefsToCDS(SequenceI cdsSeq,
           SequenceI contig, SequenceI proteinProduct, Mapping mapping)
   {
+
     // gather direct refs from contig congruent with mapping
     List<DBRefEntry> direct = new ArrayList<>();
     HashSet<String> directSources = new HashSet<>();
+
     if (contig.getDBRefs() != null)
     {
       for (DBRefEntry dbr : contig.getDBRefs())
@@ -2233,12 +2250,13 @@ public class AlignmentUtils
     int mappedDnaLength = MappingUtils.getLength(ranges);
 
     /*
-     * if not a whole number of codons, something is wrong,
-     * abort mapping
+     * if not a whole number of codons, truncate mapping
      */
-    if (mappedDnaLength % CODON_LENGTH > 0)
+    int codonRemainder = mappedDnaLength % CODON_LENGTH;
+    if (codonRemainder > 0)
     {
-      return null;
+      mappedDnaLength -= codonRemainder;
+      MappingUtils.removeEndPositions(codonRemainder, ranges);
     }
 
     int proteinLength = proteinSeq.getLength();
@@ -2410,7 +2428,8 @@ public class AlignmentUtils
   static int computePeptideVariants(SequenceI peptide, int peptidePos,
           List<DnaVariant>[] codonVariants)
   {
-    String residue = String.valueOf(peptide.getCharAt(peptidePos - 1));
+    String residue = String
+            .valueOf(peptide.getCharAt(peptidePos - peptide.getStart()));
     int count = 0;
     String base1 = codonVariants[0].get(0).base;
     String base2 = codonVariants[1].get(0).base;
@@ -2428,10 +2447,17 @@ public class AlignmentUtils
         {
           for (String base : alleles.split(","))
           {
-            String codon = base + base2 + base3;
-            if (addPeptideVariant(peptide, peptidePos, residue, var, codon))
+            if (!base1.equalsIgnoreCase(base))
             {
-              count++;
+              String codon = base.toUpperCase() + base2.toLowerCase()
+                      + base3.toLowerCase();
+              String canonical = base1.toUpperCase() + base2.toLowerCase()
+                      + base3.toLowerCase();
+              if (addPeptideVariant(peptide, peptidePos, residue, var,
+                      codon, canonical))
+              {
+                count++;
+              }
             }
           }
         }
@@ -2450,10 +2476,17 @@ public class AlignmentUtils
         {
           for (String base : alleles.split(","))
           {
-            String codon = base1 + base + base3;
-            if (addPeptideVariant(peptide, peptidePos, residue, var, codon))
+            if (!base2.equalsIgnoreCase(base))
             {
-              count++;
+              String codon = base1.toLowerCase() + base.toUpperCase()
+                      + base3.toLowerCase();
+              String canonical = base1.toLowerCase() + base2.toUpperCase()
+                      + base3.toLowerCase();
+              if (addPeptideVariant(peptide, peptidePos, residue, var,
+                      codon, canonical))
+              {
+                count++;
+              }
             }
           }
         }
@@ -2472,10 +2505,17 @@ public class AlignmentUtils
         {
           for (String base : alleles.split(","))
           {
-            String codon = base1 + base2 + base;
-            if (addPeptideVariant(peptide, peptidePos, residue, var, codon))
+            if (!base3.equalsIgnoreCase(base))
             {
-              count++;
+              String codon = base1.toLowerCase() + base2.toLowerCase()
+                      + base.toUpperCase();
+              String canonical = base1.toLowerCase() + base2.toLowerCase()
+                      + base3.toUpperCase();
+              if (addPeptideVariant(peptide, peptidePos, residue, var,
+                      codon, canonical))
+              {
+                count++;
+              }
             }
           }
         }
@@ -2486,20 +2526,22 @@ public class AlignmentUtils
   }
 
   /**
-   * Helper method that adds a peptide variant feature, provided the given codon
-   * translates to a value different to the current residue (is a non-synonymous
-   * variant). ID and clinical_significance attributes of the dna variant (if
-   * present) are copied to the new feature.
+   * Helper method that adds a peptide variant feature. ID and
+   * clinical_significance attributes of the dna variant (if present) are copied
+   * to the new feature.
    * 
    * @param peptide
    * @param peptidePos
    * @param residue
    * @param var
    * @param codon
+   *          the variant codon e.g. aCg
+   * @param canonical
+   *          the 'normal' codon e.g. aTg
    * @return true if a feature was added, else false
    */
   static boolean addPeptideVariant(SequenceI peptide, int peptidePos,
-          String residue, DnaVariant var, String codon)
+          String residue, DnaVariant var, String codon, String canonical)
   {
     /*
      * get peptide translation of codon e.g. GAT -> D
@@ -2507,57 +2549,71 @@ public class AlignmentUtils
      * e.g. multibase variants or HGMD_MUTATION etc
      * are currently ignored here
      */
-    String trans = codon.contains("-") ? "-"
+    String trans = codon.contains("-") ? null
             : (codon.length() > CODON_LENGTH ? null
                     : ResidueProperties.codonTranslate(codon));
-    if (trans != null && !trans.equals(residue))
+    if (trans == null)
+    {
+      return false;
+    }
+    String desc = canonical + "/" + codon;
+    String featureType = "";
+    if (trans.equals(residue))
+    {
+      featureType = SequenceOntologyI.SYNONYMOUS_VARIANT;
+    }
+    else if (ResidueProperties.STOP.equals(trans))
+    {
+      featureType = SequenceOntologyI.STOP_GAINED;
+    }
+    else
     {
       String residue3Char = StringUtils
               .toSentenceCase(ResidueProperties.aa2Triplet.get(residue));
       String trans3Char = StringUtils
               .toSentenceCase(ResidueProperties.aa2Triplet.get(trans));
-      String desc = "p." + residue3Char + peptidePos + trans3Char;
-      SequenceFeature sf = new SequenceFeature(
-              SequenceOntologyI.SEQUENCE_VARIANT, desc, peptidePos,
-              peptidePos, var.getSource());
-      StringBuilder attributes = new StringBuilder(32);
-      String id = (String) var.variant.getValue(ID);
-      if (id != null)
-      {
-        if (id.startsWith(SEQUENCE_VARIANT))
-        {
-          id = id.substring(SEQUENCE_VARIANT.length());
-        }
-        sf.setValue(ID, id);
-        attributes.append(ID).append("=").append(id);
-        // TODO handle other species variants JAL-2064
-        StringBuilder link = new StringBuilder(32);
-        try
-        {
-          link.append(desc).append(" ").append(id).append(
-                  "|http://www.ensembl.org/Homo_sapiens/Variation/Summary?v=")
-                  .append(URLEncoder.encode(id, "UTF-8"));
-          sf.addLink(link.toString());
-        } catch (UnsupportedEncodingException e)
-        {
-          // as if
-        }
-      }
-      String clinSig = (String) var.variant.getValue(CLINICAL_SIGNIFICANCE);
-      if (clinSig != null)
+      desc = "p." + residue3Char + peptidePos + trans3Char;
+      featureType = SequenceOntologyI.NONSYNONYMOUS_VARIANT;
+    }
+    SequenceFeature sf = new SequenceFeature(featureType, desc, peptidePos,
+            peptidePos, var.getSource());
+
+    StringBuilder attributes = new StringBuilder(32);
+    String id = (String) var.variant.getValue(VARIANT_ID);
+    if (id != null)
+    {
+      if (id.startsWith(SEQUENCE_VARIANT))
       {
-        sf.setValue(CLINICAL_SIGNIFICANCE, clinSig);
-        attributes.append(";").append(CLINICAL_SIGNIFICANCE).append("=")
-                .append(clinSig);
+        id = id.substring(SEQUENCE_VARIANT.length());
       }
-      peptide.addSequenceFeature(sf);
-      if (attributes.length() > 0)
+      sf.setValue(VARIANT_ID, id);
+      attributes.append(VARIANT_ID).append("=").append(id);
+      // TODO handle other species variants JAL-2064
+      StringBuilder link = new StringBuilder(32);
+      try
+      {
+        link.append(desc).append(" ").append(id).append(
+                "|http://www.ensembl.org/Homo_sapiens/Variation/Summary?v=")
+                .append(URLEncoder.encode(id, "UTF-8"));
+        sf.addLink(link.toString());
+      } catch (UnsupportedEncodingException e)
       {
-        sf.setAttributes(attributes.toString());
+        // as if
       }
-      return true;
     }
-    return false;
+    String clinSig = (String) var.variant.getValue(CLINICAL_SIGNIFICANCE);
+    if (clinSig != null)
+    {
+      sf.setValue(CLINICAL_SIGNIFICANCE, clinSig);
+      attributes.append(";").append(CLINICAL_SIGNIFICANCE).append("=")
+              .append(clinSig);
+    }
+    peptide.addSequenceFeature(sf);
+    if (attributes.length() > 0)
+    {
+      sf.setAttributes(attributes.toString());
+    }
+    return true;
   }
 
   /**