JAL-3761 extra check that VCF feature is fully mappable to sequence
[jalview.git] / src / jalview / io / vcf / VCFLoader.java
index cbdd66c..f965c4e 100644 (file)
@@ -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.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;
-
 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
@@ -654,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;
@@ -877,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();
@@ -892,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);
       }
     }