X-Git-Url: http://source.jalview.org/gitweb/?a=blobdiff_plain;f=src%2Fjalview%2Fio%2Fvcf%2FVCFLoader.java;h=d6a530b7956a086a11c8c8d406e65d0f8b088327;hb=2305d687bac5ed8a76ecb93c54c7b090a928e362;hp=bb2948d265e6a4b4745474d9b69c69775fa1b1e5;hpb=94fe1bfb6ca65c8a787a336ebffde44df6795803;p=jalview.git

diff --git a/src/jalview/io/vcf/VCFLoader.java b/src/jalview/io/vcf/VCFLoader.java
index bb2948d..d6a530b 100644
--- a/src/jalview/io/vcf/VCFLoader.java
+++ b/src/jalview/io/vcf/VCFLoader.java
@@ -1,3 +1,23 @@
+/*
+ * 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;
@@ -19,16 +39,19 @@ 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.io.UnsupportedEncodingException;
-import java.net.URLDecoder;
 import java.util.ArrayList;
 import java.util.HashMap;
+import java.util.HashSet;
+import java.util.Iterator;
 import java.util.List;
+import java.util.Locale;
 import java.util.Map;
 import java.util.Map.Entry;
+import java.util.Set;
 import java.util.regex.Pattern;
 import java.util.regex.PatternSyntaxException;
 
@@ -36,8 +59,10 @@ 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;
@@ -52,17 +77,20 @@ import htsjdk.variant.vcf.VCFInfoHeaderLine;
  */
 public class VCFLoader
 {
-  private static final String ENCODED_COMMA = "%2C";
+  private static final String VCF_ENCODABLE = ":;=%,";
 
-  private static final String ENCODED_PERCENT = "%25";
+  /*
+   * Jalview feature attributes for VCF fixed column data
+   */
+  private static final String VCF_POS = "POS";
 
-  private static final String ENCODED_EQUALS = "%3D";
+  private static final String VCF_ID = "ID";
 
-  private static final String ENCODED_SEMICOLON = "%3B";
+  private static final String VCF_QUAL = "QUAL";
 
-  private static final String ENCODED_COLON = "%3A";
+  private static final String VCF_FILTER = "FILTER";
 
-  private static final String UTF_8 = "UTF-8";
+  private static final String NO_VALUE = VCFConstants.MISSING_VALUE_v4; // '.'
 
   private static final String DEFAULT_SPECIES = "homo_sapiens";
 
@@ -222,6 +250,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
    * 
@@ -536,7 +570,7 @@ public class VCFLoader
   {
     for (Pattern p : filters)
     {
-      if (p.matcher(id.toUpperCase()).matches())
+      if (p.matcher(id.toUpperCase(Locale.ROOT)).matches())
       {
         return true;
       }
@@ -630,7 +664,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);
@@ -853,24 +887,35 @@ 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);
+          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);
+          }
         }
+        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());
+        Cache.log.error(msg);
       }
-      variants.close();
     }
 
     return count;
@@ -1000,7 +1045,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);
 
@@ -1010,6 +1068,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:
@@ -1239,41 +1344,82 @@ public class VCFLoader
        * take the index'th value
        */
       String value = getAttributeValue(variant, key, index);
-      if (value != null)
+      if (value != null && isValid(variant, key, value))
       {
-        value = decodeSpecialCharacters(value);
-        sf.setValue(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);
       }
     }
   }
 
   /**
-   * Decodes colon, semicolon, equals sign, percent sign, comma to their decoded
-   * form. The VCF specification (para 1.2) requires these to be encoded where not
-   * used with their special meaning in the VCF syntax. Note that general URL
-   * decoding should not be applied, since this would incorrectly decode (for
-   * example) a '+' sign.
+   * 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 static String decodeSpecialCharacters(String value)
+  protected boolean isValid(VariantContext variant, String infoId,
+          String value)
   {
-    /*
-     * avoid regex compilation if it is not needed!
-     */
-    if (!value.contains(ENCODED_COLON) && !value.contains(ENCODED_SEMICOLON)
-            && !value.contains(ENCODED_EQUALS)
-            && !value.contains(ENCODED_PERCENT)
-            && !value.contains(ENCODED_COMMA))
+    if (value == null || value.isEmpty() || NO_VALUE.equals(value))
+    {
+      return true;
+    }
+    VCFInfoHeaderLine infoHeader = header.getInfoHeaderLine(infoId);
+    if (infoHeader == null)
     {
-      return value;
+      Cache.log.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)
+      {
+        Float.parseFloat(value);
+      }
+    } catch (NumberFormatException e)
+    {
+      logInvalidValue(variant, infoId, value);
+      return false;
+    }
+    return true;
+  }
 
-    value = value.replace(ENCODED_COLON, ":")
-            .replace(ENCODED_SEMICOLON, ";").replace(ENCODED_EQUALS, "=")
-            .replace(ENCODED_PERCENT, "%").replace(ENCODED_COMMA, ",");
-    return value;
+  /**
+   * 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);
+      Cache.log.error(String.format("Invalid VCF data at %s:%d %s=%s",
+              variant.getContig(), variant.getStart(), infoId, value));
+    }
   }
 
   /**
@@ -1328,12 +1474,7 @@ public class VCFLoader
                * VCF spec requires encoding of special characters e.g. '='
                * so decode them here before storing
                */
-              try
-              {
-                field = URLDecoder.decode(field, UTF_8);
-              } catch (UnsupportedEncodingException e)
-              {
-              }
+              field = StringUtils.urlDecode(field, VCF_ENCODABLE);
               csqValues.put(id, field);
             }
           }