JAL-2738 match consequence "Allele" to VCF allele; add non-SNP variants
[jalview.git] / src / jalview / analysis / AlignmentUtils.java
index 2b9b9f9..9a46a91 100644 (file)
@@ -36,6 +36,7 @@ import jalview.datamodel.SequenceFeature;
 import jalview.datamodel.SequenceGroup;
 import jalview.datamodel.SequenceI;
 import jalview.datamodel.features.SequenceFeatures;
+import jalview.io.gff.Gff3Helper;
 import jalview.io.gff.SequenceOntologyI;
 import jalview.schemes.ResidueProperties;
 import jalview.util.Comparison;
@@ -105,6 +106,15 @@ public class AlignmentUtils
     {
       return variant == null ? null : variant.getFeatureGroup();
     }
+
+    /**
+     * toString for aid in the debugger only
+     */
+    @Override
+    public String toString()
+    {
+      return base + ":" + (variant == null ? "" : variant.getDescription());
+    }
   }
 
   /**
@@ -1623,7 +1633,7 @@ public class AlignmentUtils
           AlignmentI dataset, SequenceI[] products)
   {
     if (dataset == null || dataset.getDataset() != null)
-    {
+  {
       throw new IllegalArgumentException(
               "IMPLEMENTATION ERROR: dataset.getDataset() must be null!");
     }
@@ -1635,10 +1645,10 @@ public class AlignmentUtils
     {
       productSeqs = new HashSet<SequenceI>();
       for (SequenceI seq : products)
-      {
-        productSeqs.add(seq.getDatasetSequence() == null ? seq
-                : seq.getDatasetSequence());
-      }
+    {
+        productSeqs.add(seq.getDatasetSequence() == null ? seq : seq
+                .getDatasetSequence());
+    }
     }
 
     /*
@@ -1660,15 +1670,15 @@ public class AlignmentUtils
       List<AlignedCodonFrame> seqMappings = MappingUtils
               .findMappingsForSequence(dnaSeq, mappings);
       for (AlignedCodonFrame mapping : seqMappings)
-      {
+    {
         List<Mapping> mappingsFromSequence = mapping
                 .getMappingsFromSequence(dnaSeq);
 
         for (Mapping aMapping : mappingsFromSequence)
-        {
+      {
           MapList mapList = aMapping.getMap();
           if (mapList.getFromRatio() == 1)
-          {
+        {
             /*
              * not a dna-to-protein mapping (likely dna-to-cds)
              */
@@ -1694,15 +1704,15 @@ public class AlignmentUtils
           if (cdsSeq != null)
           {
             if (!foundSeqs.contains(cdsSeq))
-            {
+          {
               foundSeqs.add(cdsSeq);
               SequenceI derivedSequence = cdsSeq.deriveSequence();
               cdsSeqs.add(derivedSequence);
               if (!dataset.getSequences().contains(cdsSeq))
-              {
+            {
                 dataset.addSequence(cdsSeq);
-              }
             }
+          }
             continue;
           }
 
@@ -1730,9 +1740,8 @@ public class AlignmentUtils
           /*
            * add a mapping from CDS to the (unchanged) mapped to range
            */
-          List<int[]> cdsRange = Collections
-                  .singletonList(new int[]
-                  { 1, cdsSeq.getLength() });
+          List<int[]> cdsRange = Collections.singletonList(new int[] { 1,
+              cdsSeq.getLength() });
           MapList cdsToProteinMap = new MapList(cdsRange,
                   mapList.getToRanges(), mapList.getFromRatio(),
                   mapList.getToRatio());
@@ -1782,40 +1791,46 @@ public class AlignmentUtils
 
           for (DBRefEntry primRef : dnaDss.getPrimaryDBRefs())
           {
-            // creates a complementary cross-reference to the source sequence's
-            // primary reference.
-
-            DBRefEntry cdsCrossRef = new DBRefEntry(primRef.getSource(),
-                    primRef.getSource() + ":" + primRef.getVersion(),
-                    primRef.getAccessionId());
-            cdsCrossRef
-                    .setMap(new Mapping(dnaDss, new MapList(dnaToCdsMap)));
+            /*
+             * create a cross-reference from CDS to the source sequence's
+             * primary reference and vice versa
+             */
+
+            String source = primRef.getSource();
+            String version = primRef.getVersion();
+            DBRefEntry cdsCrossRef = new DBRefEntry(source, source + ":"
+                    + version, primRef.getAccessionId());
+            cdsCrossRef.setMap(new Mapping(dnaDss, new MapList(dnaToCdsMap
+                    .getInverse())));
             cdsSeqDss.addDBRef(cdsCrossRef);
 
+            dnaSeq.addDBRef(new DBRefEntry(source, version, cdsSeq
+                    .getName(), new Mapping(cdsSeqDss, dnaToCdsMap)));
+
             // problem here is that the cross-reference is synthesized -
             // cdsSeq.getName() may be like 'CDS|dnaaccession' or
             // 'CDS|emblcdsacc'
             // assuming cds version same as dna ?!?
 
-            DBRefEntry proteinToCdsRef = new DBRefEntry(primRef.getSource(),
-                    primRef.getVersion(), cdsSeq.getName());
+            DBRefEntry proteinToCdsRef = new DBRefEntry(source, version,
+                    cdsSeq.getName());
             //
-            proteinToCdsRef.setMap(
-                    new Mapping(cdsSeqDss, cdsToProteinMap.getInverse()));
+            proteinToCdsRef.setMap(new Mapping(cdsSeqDss, cdsToProteinMap
+                    .getInverse()));
             proteinProduct.addDBRef(proteinToCdsRef);
-          }
+        }
 
           /*
            * transfer any features on dna that overlap the CDS
            */
           transferFeatures(dnaSeq, cdsSeq, dnaToCdsMap, null,
                   SequenceOntologyI.CDS);
-        }
       }
     }
+    }
 
-    AlignmentI cds = new Alignment(
-            cdsSeqs.toArray(new SequenceI[cdsSeqs.size()]));
+    AlignmentI cds = new Alignment(cdsSeqs.toArray(new SequenceI[cdsSeqs
+            .size()]));
     cds.setDataset(dataset);
 
     return cds;
@@ -2365,7 +2380,7 @@ public class AlignmentUtils
     {
       if (var.variant != null)
       {
-        String alleles = (String) var.variant.getValue("alleles");
+        String alleles = (String) var.variant.getValue(Gff3Helper.ALLELES);
         if (alleles != null)
         {
           for (String base : alleles.split(","))
@@ -2387,7 +2402,7 @@ public class AlignmentUtils
     {
       if (var.variant != null)
       {
-        String alleles = (String) var.variant.getValue("alleles");
+        String alleles = (String) var.variant.getValue(Gff3Helper.ALLELES);
         if (alleles != null)
         {
           for (String base : alleles.split(","))
@@ -2409,7 +2424,7 @@ public class AlignmentUtils
     {
       if (var.variant != null)
       {
-        String alleles = (String) var.variant.getValue("alleles");
+        String alleles = (String) var.variant.getValue(Gff3Helper.ALLELES);
         if (alleles != null)
         {
           for (String base : alleles.split(","))
@@ -2504,7 +2519,10 @@ public class AlignmentUtils
 
   /**
    * Builds a map whose key is position in the protein sequence, and value is a
-   * list of the base and all variants for each corresponding codon position
+   * list of the base and all variants for each corresponding codon position.
+   * <p>
+   * This depends on dna variants being held as a comma-separated list as
+   * property "alleles" on variant features.
    * 
    * @param dnaSeq
    * @param dnaToProtein
@@ -2542,6 +2560,30 @@ public class AlignmentUtils
         // not handling multi-locus variant features
         continue;
       }
+
+      /*
+       * ignore variant if not a SNP
+       */
+      String alls = (String) sf.getValue(Gff3Helper.ALLELES);
+      if (alls == null)
+      {
+        continue; // non-SNP VCF variant perhaps - can't process this
+      }
+
+      String[] alleles = alls.toUpperCase().split(",");
+      boolean isSnp = true;
+      for (String allele : alleles)
+      {
+        if (allele.trim().length() > 1)
+        {
+          isSnp = false;
+        }
+      }
+      if (!isSnp)
+      {
+        continue;
+      }
+
       int[] mapsTo = dnaToProtein.locateInTo(dnaCol, dnaCol);
       if (mapsTo == null)
       {
@@ -2560,21 +2602,6 @@ public class AlignmentUtils
       }
 
       /*
-       * extract dna variants to a string array
-       */
-      String alls = (String) sf.getValue("alleles");
-      if (alls == null)
-      {
-        continue;
-      }
-      String[] alleles = alls.toUpperCase().split(",");
-      int i = 0;
-      for (String allele : alleles)
-      {
-        alleles[i++] = allele.trim(); // lose any space characters "A, G"
-      }
-
-      /*
        * get this peptide's codon positions e.g. [3, 4, 5] or [4, 7, 10]
        */
       int[] codon = peptidePosition == lastPeptidePostion ? lastCodon