X-Git-Url: http://source.jalview.org/gitweb/?a=blobdiff_plain;f=src%2Fjalview%2Fio%2Fvcf%2FVCFLoader.java;h=29b300486332c456068cfad098ccfd620f828122;hb=cb8e52fbbc5f725e3f7f48c672cdddb0690bd978;hp=7bf7791f2b5d98a916313b804098aa3d088e10e0;hpb=cbeb7ad59d51b468c54ca3db2a2a7693060a2509;p=jalview.git

diff --git a/src/jalview/io/vcf/VCFLoader.java b/src/jalview/io/vcf/VCFLoader.java
index 7bf7791..29b3004 100644
--- a/src/jalview/io/vcf/VCFLoader.java
+++ b/src/jalview/io/vcf/VCFLoader.java
@@ -1,33 +1,37 @@
+/*
+ * 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.util.Locale;
 
 import java.io.File;
 import java.io.IOException;
 import java.util.ArrayList;
 import java.util.HashMap;
+import java.util.HashSet;
+import java.util.Iterator;
 import java.util.List;
 import java.util.Map;
 import java.util.Map.Entry;
+import java.util.Set;
 import java.util.regex.Pattern;
 import java.util.regex.PatternSyntaxException;
 
@@ -35,13 +39,36 @@ import htsjdk.samtools.SAMException;
 import htsjdk.samtools.SAMSequenceDictionary;
 import htsjdk.samtools.SAMSequenceRecord;
 import htsjdk.samtools.util.CloseableIterator;
+import htsjdk.tribble.TribbleException;
 import htsjdk.variant.variantcontext.Allele;
 import htsjdk.variant.variantcontext.VariantContext;
+import htsjdk.variant.vcf.VCFConstants;
 import htsjdk.variant.vcf.VCFHeader;
 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.bin.Console;
+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
@@ -51,6 +78,21 @@ import htsjdk.variant.vcf.VCFInfoHeaderLine;
  */
 public class VCFLoader
 {
+  private static final String VCF_ENCODABLE = ":;=%,";
+
+  /*
+   * Jalview feature attributes for VCF fixed column data
+   */
+  private static final String VCF_POS = "POS";
+
+  private static final String VCF_ID = "ID";
+
+  private static final String VCF_QUAL = "QUAL";
+
+  private static final String VCF_FILTER = "FILTER";
+
+  private static final String NO_VALUE = VCFConstants.MISSING_VALUE_v4; // '.'
+
   private static final String DEFAULT_SPECIES = "homo_sapiens";
 
   /**
@@ -209,6 +251,12 @@ public class VCFLoader
    */
   Map<Integer, String> vepFieldsOfInterest;
 
+  /*
+   * key:value for which rejected data has been seen
+   * (the error is logged only once for each combination)
+   */
+  private Set<String> badData;
+
   /**
    * Constructor given a VCF file
    * 
@@ -267,7 +315,7 @@ public class VCFLoader
     VCFHeaderLine headerLine = header.getOtherHeaderLine(VCFHeader.REFERENCE_KEY);
     if (headerLine == null)
     {
-      Cache.log.error("VCF reference header not found");
+      Console.error("VCF reference header not found");
       return null;
     }
     String ref = headerLine.getValue();
@@ -289,7 +337,7 @@ public class VCFLoader
     }
     else
     {
-      Cache.log.error("VCF reference not found: " + ref);
+      Console.error("VCF reference not found: " + ref);
     }
 
     return seq;
@@ -382,11 +430,11 @@ public class VCFLoader
   {
     if (reference == null)
     {
-      Cache.log.error("No VCF ##reference found, defaulting to "
+      Console.error("No VCF ##reference found, defaulting to "
               + DEFAULT_REFERENCE + ":" + DEFAULT_SPECIES);
       reference = DEFAULT_REFERENCE; // default to GRCh37 if not specified
     }
-    reference = reference.toLowerCase();
+    reference = reference.toLowerCase(Locale.ROOT);
 
     /*
      * for a non-human species, or other assembly identifier,
@@ -401,7 +449,7 @@ public class VCFLoader
       String[] tokens = token.split("=");
       if (tokens.length == 2)
       {
-        if (reference.contains(tokens[0].trim().toLowerCase()))
+        if (reference.contains(tokens[0].trim().toLowerCase(Locale.ROOT)))
         {
           vcfAssembly = tokens[1].trim();
           break;
@@ -418,7 +466,7 @@ public class VCFLoader
         String[] tokens = token.split("=");
         if (tokens.length == 2)
         {
-          if (reference.contains(tokens[0].trim().toLowerCase()))
+          if (reference.contains(tokens[0].trim().toLowerCase(Locale.ROOT)))
           {
             vcfSpecies = tokens[1].trim();
             break;
@@ -523,7 +571,7 @@ public class VCFLoader
   {
     for (Pattern p : filters)
     {
-      if (p.matcher(id.toUpperCase()).matches())
+      if (p.matcher(id.toUpperCase(Locale.ROOT)).matches())
       {
         return true;
       }
@@ -617,7 +665,7 @@ public class VCFLoader
     {
       try
       {
-      patterns.add(Pattern.compile(token.toUpperCase()));
+      patterns.add(Pattern.compile(token.toUpperCase(Locale.ROOT)));
       } catch (PatternSyntaxException e)
       {
         System.err.println("Invalid pattern ignored: " + token);
@@ -628,13 +676,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;
@@ -654,7 +701,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
       {
@@ -728,7 +776,7 @@ public class VCFLoader
     GeneLociI seqCoords = seq.getGeneLoci();
     if (seqCoords == null)
     {
-      Cache.log.warn(String.format(
+      Console.warn(String.format(
               "Can't query VCF for %s as chromosome coordinates not known",
               seq.getName()));
       return null;
@@ -743,7 +791,7 @@ public class VCFLoader
     // returned with the Ensembl sequence; todo: support aliases?
     if (!vcfSpecies.equalsIgnoreCase(species))
     {
-      Cache.log.warn("No VCF loaded to " + seq.getName()
+      Console.warn("No VCF loaded to " + seq.getName()
               + " as species not matched");
       return null;
     }
@@ -773,7 +821,7 @@ public class VCFLoader
               vcfAssembly);
       if (newRange == null)
       {
-        Cache.log.error(
+        Console.error(
                 String.format("Failed to map %s:%s:%s:%d:%d to %s", species,
                         chromosome, seqRef, range[0], range[1],
                         vcfAssembly));
@@ -839,24 +887,42 @@ public class VCFLoader
     {
       int vcfStart = Math.min(range[0], range[1]);
       int vcfEnd = Math.max(range[0], range[1]);
-      CloseableIterator<VariantContext> variants = reader
-              .query(map.chromosome, vcfStart, vcfEnd);
-      while (variants.hasNext())
+      try
       {
-        VariantContext variant = variants.next();
+        CloseableIterator<VariantContext> variants = reader
+                .query(map.chromosome, vcfStart, vcfEnd);
+        while (variants.hasNext())
+        {
+          VariantContext variant = variants.next();
 
-        int[] featureRange = map.map.locateInFrom(variant.getStart(),
-                variant.getEnd());
+          int[] featureRange = map.map.locateInFrom(variant.getStart(),
+                  variant.getEnd());
 
-        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);
+          /*
+           * 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]);
+            if (featureEnd - featureStart == variant.getEnd()
+                    - variant.getStart())
+            {
+              count += addAlleleFeatures(seq, variant, featureStart,
+                      featureEnd, forwardStrand);
+            }
+          }
         }
+        variants.close();
+      } catch (TribbleException e)
+      {
+        /*
+         * RuntimeException throwable by htsjdk
+         */
+        String msg = String.format("Error reading VCF for %s:%d-%d: %s ",
+                map.chromosome, vcfStart, vcfEnd,e.getLocalizedMessage());
+        Console.error(msg);
       }
-      variants.close();
     }
 
     return count;
@@ -970,7 +1036,7 @@ public class VCFLoader
      */
     String consequence = getConsequenceForAlleleAndFeature(variant, CSQ_FIELD,
             altAlleleIndex, csqAlleleFieldIndex,
-            csqAlleleNumberFieldIndex, seq.getName().toLowerCase(),
+            csqAlleleNumberFieldIndex, seq.getName().toLowerCase(Locale.ROOT),
             csqFeatureFieldIndex);
 
     /*
@@ -986,7 +1052,20 @@ public class VCFLoader
             featureEnd, FEATURE_GROUP_VCF);
     sf.setSource(sourceId);
 
-    sf.setValue(Gff3Helper.ALLELES, alleles);
+    /*
+     * save the derived alleles as a named attribute; this will be
+     * needed when Jalview computes derived peptide variants
+     */
+    addFeatureAttribute(sf, Gff3Helper.ALLELES, alleles);
+
+    /*
+     * add selected VCF fixed column data as feature attributes
+     */
+    addFeatureAttribute(sf, VCF_POS, String.valueOf(variant.getStart()));
+    addFeatureAttribute(sf, VCF_ID, variant.getID());
+    addFeatureAttribute(sf, VCF_QUAL,
+            String.valueOf(variant.getPhredScaledQual()));
+    addFeatureAttribute(sf, VCF_FILTER, getFilter(variant));
 
     addAlleleProperties(variant, sf, altAlleleIndex, consequence);
 
@@ -996,6 +1075,53 @@ public class VCFLoader
   }
 
   /**
+   * Answers the VCF FILTER value for the variant - or an approximation to it.
+   * This field is either PASS, or a semi-colon separated list of filters not
+   * passed. htsjdk saves filters as a HashSet, so the order when reassembled into
+   * a list may be different.
+   * 
+   * @param variant
+   * @return
+   */
+  String getFilter(VariantContext variant)
+  {
+    Set<String> filters = variant.getFilters();
+    if (filters.isEmpty())
+    {
+      return NO_VALUE;
+    }
+    Iterator<String> iterator = filters.iterator();
+    String first = iterator.next();
+    if (filters.size() == 1)
+    {
+      return first;
+    }
+
+    StringBuilder sb = new StringBuilder(first);
+    while (iterator.hasNext())
+    {
+      sb.append(";").append(iterator.next());
+    }
+
+    return sb.toString();
+  }
+
+  /**
+   * Adds one feature attribute unless the value is null, empty or '.'
+   * 
+   * @param sf
+   * @param key
+   * @param value
+   */
+  void addFeatureAttribute(SequenceFeature sf, String key, String value)
+  {
+    if (value != null && !value.isEmpty() && !NO_VALUE.equals(value))
+    {
+      sf.setValue(key, value);
+    }
+  }
+
+  /**
    * Determines the Sequence Ontology term to use for the variant feature type in
    * Jalview. The default is 'sequence_variant', but a more specific term is used
    * if:
@@ -1103,7 +1229,7 @@ public class VCFLoader
       {
         String featureIdentifier = csqFields[featureFieldIndex];
         if (featureIdentifier.length() > 4
-                && seqName.indexOf(featureIdentifier.toLowerCase()) > -1)
+                && seqName.indexOf(featureIdentifier.toLowerCase(Locale.ROOT)) > -1)
         {
           /*
            * feature (transcript) matched - now check for allele match
@@ -1189,14 +1315,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') 
@@ -1233,10 +1351,81 @@ public class VCFLoader
        * take the index'th value
        */
       String value = getAttributeValue(variant, key, index);
-      if (value != null)
+      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);
+      }
+    }
+  }
+
+  /**
+   * Answers true for '.', null, or an empty value, or if the INFO type is String.
+   * If the INFO type is Integer or Float, answers false if the value is not in
+   * valid format.
+   * 
+   * @param variant
+   * @param infoId
+   * @param value
+   * @return
+   */
+  protected boolean isValid(VariantContext variant, String infoId,
+          String value)
+  {
+    if (value == null || value.isEmpty() || NO_VALUE.equals(value))
+    {
+      return true;
+    }
+    VCFInfoHeaderLine infoHeader = header.getInfoHeaderLine(infoId);
+    if (infoHeader == null)
+    {
+      Console.error("Field " + infoId + " has no INFO header");
+      return false;
+    }
+    VCFHeaderLineType infoType = infoHeader.getType();
+    try
+    {
+      if (infoType == VCFHeaderLineType.Integer)
+      {
+        Integer.parseInt(value);
+      }
+      else if (infoType == VCFHeaderLineType.Float)
       {
-        sf.setValue(key, value);
+        Float.parseFloat(value);
       }
+    } catch (NumberFormatException e)
+    {
+      logInvalidValue(variant, infoId, value);
+      return false;
+    }
+    return true;
+  }
+
+  /**
+   * Logs an error message for malformed data; duplicate messages (same id and
+   * value) are not logged
+   * 
+   * @param variant
+   * @param infoId
+   * @param value
+   */
+  private void logInvalidValue(VariantContext variant, String infoId,
+          String value)
+  {
+    if (badData == null)
+    {
+      badData = new HashSet<>();
+    }
+    String token = infoId + ":" + value;
+    if (!badData.contains(token))
+    {
+      badData.add(token);
+      Console.error(String.format("Invalid VCF data at %s:%d %s=%s",
+              variant.getContig(), variant.getStart(), infoId, value));
     }
   }
 
@@ -1288,6 +1477,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);
             }
           }