JAL-1705 compute peptide variants from dna; add "consequence" feature on
[jalview.git] / src / jalview / analysis / CrossRef.java
index e96d9d7..68f6c93 100644 (file)
@@ -27,13 +27,22 @@ import jalview.datamodel.DBRefEntry;
 import jalview.datamodel.DBRefSource;
 import jalview.datamodel.Mapping;
 import jalview.datamodel.Sequence;
+import jalview.datamodel.SequenceFeature;
 import jalview.datamodel.SequenceI;
+import jalview.io.gff.SequenceOntology;
+import jalview.schemes.ResidueProperties;
 import jalview.util.DBRefUtils;
+import jalview.util.MapList;
+import jalview.util.MappingUtils;
+import jalview.util.StringUtils;
 import jalview.ws.SequenceFetcher;
 import jalview.ws.seqfetcher.ASequenceFetcher;
 
 import java.util.ArrayList;
+import java.util.Collections;
+import java.util.LinkedHashMap;
 import java.util.List;
+import java.util.Map.Entry;
 import java.util.Vector;
 
 /**
@@ -269,6 +278,12 @@ public class CrossRef
                 // map should be from protein seq to its coding dna
                 cf.addMap(rsq, dss, xref.getMap().getMap().getInverse());
               }
+
+              /*
+               * compute peptide variants from dna variants
+               */
+              rsq.createDatasetSequence();
+              computeProteinVariants(seq, rsq, xref.getMap().getMap());
             }
             found = true;
           }
@@ -570,6 +585,195 @@ public class CrossRef
   }
 
   /**
+   * Computes variants in peptide product generated by variants in dna, and adds
+   * them as sequence_variant features on the protein sequence. Returns the
+   * number of variant features added.
+   * 
+   * @param dnaSeq
+   * @param peptide
+   * @param dnaToProtein
+   */
+  protected static int computeProteinVariants(SequenceI dnaSeq,
+          SequenceI peptide, MapList dnaToProtein)
+  {
+    /*
+     * map from peptide position to all variant features of the codon for it
+     * LinkedHashMap ensures we add the peptide features in sequence order
+     */
+    LinkedHashMap<Integer, String[][]> variants = new LinkedHashMap<Integer, String[][]>();
+    SequenceOntology so = SequenceOntology.getInstance();
+  
+    SequenceFeature[] dnaFeatures = dnaSeq.getSequenceFeatures();
+    if (dnaFeatures == null)
+    {
+      return 0;
+    }
+  
+    int[] lastCodon = null;
+    int lastPeptidePostion = 0;
+
+    /*
+     * build a map of codon variations for peptides
+     */
+    for (SequenceFeature sf : dnaFeatures)
+    {
+      int dnaCol = sf.getBegin();
+      if (dnaCol != sf.getEnd())
+      {
+        // not handling multi-locus variant features
+        continue;
+      }
+      if (so.isSequenceVariant(sf.getType()))
+      {
+        int[] mapsTo = dnaToProtein.locateInTo(dnaCol, dnaCol);
+        if (mapsTo == null)
+        {
+          // feature doesn't lie within coding region
+          continue;
+        }
+        int peptidePosition = mapsTo[0];
+        String[][] codonVariants = variants.get(peptidePosition);
+        if (codonVariants == null)
+        {
+          codonVariants = new String[3][];
+          variants.put(peptidePosition, codonVariants);
+        }
+
+        /*
+         * extract dna variants to a string array
+         */
+        String alls = (String) sf.getValue("alleles");
+        if (alls == null)
+        {
+          continue;
+        }
+        String[] alleles = alls.split(",");
+
+        /*
+         * get this peptides codon positions e.g. [3, 4, 5] or [4, 7, 10]
+         */
+        int[] codon = peptidePosition == lastPeptidePostion ? lastCodon
+                : MappingUtils.flattenRanges(dnaToProtein.locateInFrom(
+                        peptidePosition, peptidePosition));
+        lastPeptidePostion = peptidePosition;
+        lastCodon = codon;
+
+        /*
+         * save nucleotide (and this variant) for each codon position
+         */
+        for (int codonPos = 0; codonPos < 3; codonPos++)
+        {
+          String nucleotide = String.valueOf(dnaSeq
+                  .getCharAt(codon[codonPos] - 1));
+          if (codon[codonPos] == dnaCol)
+          {
+            /*
+             * record current dna base and its alleles
+             */
+            String[] dnaVariants = new String[alleles.length + 1];
+            dnaVariants[0] = nucleotide;
+            System.arraycopy(alleles, 0, dnaVariants, 1, alleles.length);
+            codonVariants[codonPos] = dnaVariants;
+          }
+          else if (codonVariants[codonPos] == null)
+          {
+            /*
+             * record current dna base only 
+             * (at least until we find any variation and overwrite it)
+             */
+            codonVariants[codonPos] = new String[] { nucleotide };
+          }
+        }
+      }
+    }
+  
+    /*
+     * scan codon variations, compute peptide variants and add to peptide sequence
+     */
+    int count = 0;
+    for (Entry<Integer, String[][]> variant : variants.entrySet())
+    {
+      int peptidePos = variant.getKey();
+      String[][] codonVariants = variant.getValue();
+      String residue = String.valueOf(peptide.getCharAt(peptidePos - 1)); // 0-based
+      List<String> peptideVariants = computePeptideVariants(codonVariants,
+              residue);
+      if (!peptideVariants.isEmpty())
+      {
+        Collections.sort(peptideVariants);
+        String desc = StringUtils.listToDelimitedString(peptideVariants,
+                ", ");
+        SequenceFeature sf = new SequenceFeature(
+                SequenceOntology.SEQUENCE_VARIANT, desc, peptidePos,
+                peptidePos, Float.NaN, null);
+        peptide.getDatasetSequence().addSequenceFeature(sf);
+        count++;
+      }
+    }
+    return count;
+  }
+
+  /**
+   * Returns a non-redundant list of all peptide translations generated by the
+   * given dna variants, excluding the current residue value
+   * 
+   * @param codonVariants
+   *          an array of base values for codon positions 1, 2, 3
+   * @param residue
+   *          the current residue translation
+   * @return
+   */
+  protected static List<String> computePeptideVariants(
+          String[][] codonVariants, String residue)
+  {
+    List<String> result = new ArrayList<String>();
+    for (String base1 : codonVariants[0])
+    {
+      for (String base2 : codonVariants[1])
+      {
+        for (String base3 : codonVariants[2])
+        {
+          String codon = base1 + base2 + base3;
+          // TODO: report frameshift/insertion/deletion
+          // and multiple-base variants?!
+          String peptide = codon.contains("-") ? "-" : ResidueProperties
+                  .codonTranslate(codon);
+          if (peptide != null && !result.contains(peptide)
+                  && !peptide.equals(residue))
+          {
+            result.add(peptide);
+          }
+        }
+      }
+    }
+    return result;
+  }
+
+  /**
+   * Computes a list of all peptide variants given dna variants
+   * 
+   * @param dnaSeq
+   *          the coding dna sequence
+   * @param codonVariants
+   *          variant features for each codon position (null if no variant)
+   * @param residue
+   *          the canonical protein translation
+   * @return
+   */
+  protected static List<String> computePeptideVariants(SequenceI dnaSeq,
+          SequenceFeature[] codonVariants, String residue)
+  {
+    List<String> result = new ArrayList<String>();
+    int[][] dnaVariants = new int[3][];
+    for (int i = 0; i < 3; i++)
+    {
+
+    }
+    // TODO Auto-generated method stub
+    return null;
+  }
+
+  /**
    * precalculate different products that can be found for seqs in dataset and
    * return them.
    *