JAL-3761 extra check that VCF feature is fully mappable to sequence

[jalview.git] / src / jalview / io / vcf / VCFLoader.java

diff --git a/src/jalview/io/vcf/VCFLoader.java b/src/jalview/io/vcf/VCFLoader.java
index decff23..f965c4e 100644 (file)
--- a/src/jalview/io/vcf/VCFLoader.java
+++ b/src/jalview/io/vcf/VCFLoader.java
@@ -1,26 +1,25 @@
+/*
+ * Jalview - A Sequence Alignment Editor and Viewer ($$Version-Rel$$)
+ * Copyright (C) $$Year-Rel$$ The Jalview Authors
+ * 
+ * This file is part of Jalview.
+ * 
+ * Jalview is free software: you can redistribute it and/or
+ * modify it under the terms of the GNU General Public License 
+ * as published by the Free Software Foundation, either version 3
+ * of the License, or (at your option) any later version.
+ *  
+ * Jalview is distributed in the hope that it will be useful, but 
+ * WITHOUT ANY WARRANTY; without even the implied warranty 
+ * of MERCHANTABILITY or FITNESS FOR A PARTICULAR 
+ * PURPOSE.  See the GNU General Public License for more details.
+ * 
+ * You should have received a copy of the GNU General Public License
+ * along with Jalview.  If not, see <http://www.gnu.org/licenses/>.
+ * The Jalview Authors are detailed in the 'AUTHORS' file.
+ */
 package jalview.io.vcf;
 
-import jalview.analysis.AlignmentUtils;
-import jalview.analysis.Dna;
-import jalview.api.AlignViewControllerGuiI;
-import jalview.bin.Cache;
-import jalview.datamodel.DBRefEntry;
-import jalview.datamodel.GeneLociI;
-import jalview.datamodel.Mapping;
-import jalview.datamodel.SequenceFeature;
-import jalview.datamodel.SequenceI;
-import jalview.datamodel.features.FeatureAttributeType;
-import jalview.datamodel.features.FeatureSource;
-import jalview.datamodel.features.FeatureSources;
-import jalview.ext.ensembl.EnsemblMap;
-import jalview.ext.htsjdk.HtsContigDb;
-import jalview.ext.htsjdk.VCFReader;
-import jalview.io.gff.Gff3Helper;
-import jalview.io.gff.SequenceOntologyI;
-import jalview.util.MapList;
-import jalview.util.MappingUtils;
-import jalview.util.MessageManager;
-
 import java.io.File;
 import java.io.IOException;
 import java.util.ArrayList;
@@ -47,6 +46,26 @@ import htsjdk.variant.vcf.VCFHeaderLine;
 import htsjdk.variant.vcf.VCFHeaderLineCount;
 import htsjdk.variant.vcf.VCFHeaderLineType;
 import htsjdk.variant.vcf.VCFInfoHeaderLine;
+import jalview.analysis.Dna;
+import jalview.api.AlignViewControllerGuiI;
+import jalview.bin.Cache;
+import jalview.datamodel.DBRefEntry;
+import jalview.datamodel.GeneLociI;
+import jalview.datamodel.Mapping;
+import jalview.datamodel.SequenceFeature;
+import jalview.datamodel.SequenceI;
+import jalview.datamodel.features.FeatureAttributeType;
+import jalview.datamodel.features.FeatureSource;
+import jalview.datamodel.features.FeatureSources;
+import jalview.ext.ensembl.EnsemblMap;
+import jalview.ext.htsjdk.HtsContigDb;
+import jalview.ext.htsjdk.VCFReader;
+import jalview.io.gff.Gff3Helper;
+import jalview.io.gff.SequenceOntologyI;
+import jalview.util.MapList;
+import jalview.util.MappingUtils;
+import jalview.util.MessageManager;
+import jalview.util.StringUtils;
 
 /**
  * A class to read VCF data (using the htsjdk) and add variants as sequence
@@ -56,6 +75,8 @@ import htsjdk.variant.vcf.VCFInfoHeaderLine;
  */
 public class VCFLoader
 {
+  private static final String VCF_ENCODABLE = ":;=%,";
+
   /*
    * Jalview feature attributes for VCF fixed column data
    */
@@ -652,13 +673,12 @@ public class VCFLoader
 
   /**
    * Transfers VCF features to sequences to which this sequence has a mapping.
-   * If the mapping is 3:1, computes peptide variants from nucleotide variants.
    * 
    * @param seq
    */
   protected void transferAddedFeatures(SequenceI seq)
   {
-    DBRefEntry[] dbrefs = seq.getDBRefs();
+    List<DBRefEntry> dbrefs = seq.getDBRefs();
     if (dbrefs == null)
     {
       return;
@@ -678,7 +698,8 @@ public class VCFLoader
         /*
          * dna-to-peptide product mapping
          */
-        AlignmentUtils.computeProteinFeatures(seq, mapTo, map);
+        // JAL-3187 render on the fly instead
+        // AlignmentUtils.computeProteinFeatures(seq, mapTo, map);
       }
       else
       {
@@ -874,12 +895,19 @@ public class VCFLoader
           int[] featureRange = map.map.locateInFrom(variant.getStart(),
                   variant.getEnd());
 
+          /*
+           * only take features whose range is fully mappable to sequence positions
+           */
           if (featureRange != null)
           {
             int featureStart = Math.min(featureRange[0], featureRange[1]);
             int featureEnd = Math.max(featureRange[0], featureRange[1]);
-            count += addAlleleFeatures(seq, variant, featureStart,
-                    featureEnd, forwardStrand);
+            if (featureEnd - featureStart == variant.getEnd()
+                    - variant.getStart())
+            {
+              count += addAlleleFeatures(seq, variant, featureStart,
+                      featureEnd, forwardStrand);
+            }
           }
         }
         variants.close();
@@ -889,7 +917,7 @@ public class VCFLoader
          * RuntimeException throwable by htsjdk
          */
         String msg = String.format("Error reading VCF for %s:%d-%d: %s ",
-                map.chromosome, vcfStart, vcfEnd);
+                map.chromosome, vcfStart, vcfEnd,e.getLocalizedMessage());
         Cache.log.error(msg);
       }
     }
@@ -1284,14 +1312,6 @@ public class VCFLoader
       }
 
       /*
-       * filter out fields we don't want to capture
-       */
-      if (!vcfFieldsOfInterest.contains(key))
-      {
-        continue;
-      }
-
-      /*
        * we extract values for other data which are allele-specific; 
        * these may be per alternate allele (INFO[key].Number = 'A') 
        * or per allele including reference (INFO[key].Number = 'R') 
@@ -1330,6 +1350,11 @@ public class VCFLoader
       String value = getAttributeValue(variant, key, index);
       if (value != null && isValid(variant, key, value))
       {
+        /*
+         * decode colon, semicolon, equals sign, percent sign, comma (only)
+         * as required by the VCF specification (para 1.2)
+         */
+        value = StringUtils.urlDecode(value, VCF_ENCODABLE);
         addFeatureAttribute(sf, key, value);
       }
     }
@@ -1449,6 +1474,11 @@ public class VCFLoader
             String id = vepFieldsOfInterest.get(i);
             if (id != null)
             {
+              /*
+               * VCF spec requires encoding of special characters e.g. '='
+               * so decode them here before storing
+               */
+              field = StringUtils.urlDecode(field, VCF_ENCODABLE);
               csqValues.put(id, field);
             }
           }