Merge branch 'feature/JAL-3469clinvarVCF' into spike/clinvar
[jalview.git] / src / jalview / io / vcf / VCFLoader.java
index 5544bd6..1abe638 100644 (file)
@@ -1,6 +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;
@@ -20,13 +39,16 @@ 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;
 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;
@@ -55,6 +77,19 @@ 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";
@@ -222,9 +257,9 @@ public class VCFLoader
   private Set<String> badData;
 
   /**
-   * Constructor given a VCF file
+   * Constructor given a path to a VCF file
    * 
-   * @param alignment
+   * @param vcfFile
    */
   public VCFLoader(String vcfFile)
   {
@@ -535,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;
       }
@@ -629,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);
@@ -666,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
       {
@@ -877,7 +913,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);
       }
     }
@@ -886,7 +922,11 @@ public class VCFLoader
   }
 
   /**
-   * A convenience method to get an attribute value for an alternate allele
+   * A convenience method to get an attribute value for an alternate allele.
+   * {@code alleleIndex} is the position in the list of values for the allele.
+   * If {@alleleIndex == -1} then all values are concatenated (comma-separated).
+   * This is the case for fields declared with "Number=." i.e. values are not
+   * related to specific alleles.
    * 
    * @param variant
    * @param attributeName
@@ -898,16 +938,25 @@ public class VCFLoader
   {
     Object att = variant.getAttribute(attributeName);
 
+    String result = null;
     if (att instanceof String)
     {
-      return NO_VALUE.equals(att) ? null : (String) att;
+      result = (String) att;
     }
-    else if (att instanceof ArrayList)
+    else if (att instanceof List<?>)
     {
-      return ((List<String>) att).get(alleleIndex);
+      List<String> theList = (List<String>) att;
+      if (alleleIndex == -1)
+      {
+        result = StringUtils.listToDelimitedString(theList, ",");
+      }
+      else
+      {
+        result = theList.get(alleleIndex);
+      }
     }
 
-    return null;
+    return result;
   }
 
   /**
@@ -1009,7 +1058,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);
 
@@ -1019,6 +1081,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:
@@ -1212,14 +1321,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') 
@@ -1244,6 +1345,10 @@ public class VCFLoader
          */
         index++;
       }
+      else if (number == VCFHeaderLineCount.UNBOUNDED) // .
+      {
+        index = -1;
+      }
       else if (number != VCFHeaderLineCount.A)
       {
         /*
@@ -1258,7 +1363,12 @@ public class VCFLoader
       String value = getAttributeValue(variant, key, index);
       if (value != null && isValid(variant, key, 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);
       }
     }
   }
@@ -1377,6 +1487,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);
             }
           }