X-Git-Url: http://source.jalview.org/gitweb/?a=blobdiff_plain;ds=sidebyside;f=src%2Fjalview%2Fio%2Fvcf%2FVCFLoader.java;h=ac707d8f01cea046948f2809c07256b0739c0ae2;hb=456e5c765ac1e85336fd9f9b1a35453069bb2298;hp=9d98b7e6c675c2967424790c4f22e77c86161df1;hpb=b87ae5ac68939a1b964682046e8b07958fae219a;p=jalview.git

diff --git a/src/jalview/io/vcf/VCFLoader.java b/src/jalview/io/vcf/VCFLoader.java
index 9d98b7e..ac707d8 100644
--- a/src/jalview/io/vcf/VCFLoader.java
+++ b/src/jalview/io/vcf/VCFLoader.java
@@ -1,10 +1,8 @@
 package jalview.io.vcf;
 
-import jalview.analysis.AlignmentUtils;
 import jalview.analysis.Dna;
 import jalview.api.AlignViewControllerGuiI;
 import jalview.bin.Cache;
-import jalview.datamodel.AlignmentI;
 import jalview.datamodel.DBRefEntry;
 import jalview.datamodel.GeneLociI;
 import jalview.datamodel.Mapping;
@@ -14,6 +12,7 @@ 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;
@@ -21,12 +20,18 @@ import jalview.util.MapList;
 import jalview.util.MappingUtils;
 import jalview.util.MessageManager;
 
+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.Map;
 import java.util.Map.Entry;
+import java.util.Set;
 import java.util.regex.Pattern;
 import java.util.regex.PatternSyntaxException;
 
@@ -34,8 +39,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;
@@ -50,6 +57,33 @@ import htsjdk.variant.vcf.VCFInfoHeaderLine;
  */
 public class VCFLoader
 {
+  private static final String ENCODED_COMMA = "%2C";
+
+  private static final String ENCODED_PERCENT = "%25";
+
+  private static final String ENCODED_EQUALS = "%3D";
+
+  private static final String ENCODED_SEMICOLON = "%3B";
+
+  private static final String ENCODED_COLON = "%3A";
+
+  private static final String UTF_8 = "UTF-8";
+
+  /*
+   * 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";
+
   /**
    * A class to model the mapping from sequence to VCF coordinates. Cases include
    * <ul>
@@ -81,7 +115,7 @@ public class VCFLoader
 
   /*
    * Lookup keys, and default values, for Preference entries that describe
-   * patterns for VCF and VEP fields to capture 
+   * patterns for VCF and VEP fields to capture
    */
   private static final String VEP_FIELDS_PREF = "VEP_FIELDS";
 
@@ -92,6 +126,18 @@ public class VCFLoader
   private static final String DEFAULT_VEP_FIELDS = ".*";// "Allele,Consequence,IMPACT,SWISSPROT,SIFT,PolyPhen,CLIN_SIG";
 
   /*
+   * Lookup keys, and default values, for Preference entries that give
+   * mappings from tokens in the 'reference' header to species or assembly
+   */
+  private static final String VCF_ASSEMBLY = "VCF_ASSEMBLY";
+
+  private static final String DEFAULT_VCF_ASSEMBLY = "assembly19=GRCh37,hs37=GRCh37,grch37=GRCh37,grch38=GRCh38";
+
+  private static final String VCF_SPECIES = "VCF_SPECIES"; // default is human
+
+  private static final String DEFAULT_REFERENCE = "grch37"; // fallback default is human GRCh37
+
+  /*
    * keys to fields of VEP CSQ consequence data
    * see https://www.ensembl.org/info/docs/tools/vep/vep_formats.html
    */
@@ -113,12 +159,6 @@ public class VCFLoader
   private static final String PIPE_REGEX = "\\|";
 
   /*
-   * key for Allele Frequency output by VEP
-   * see http://www.ensembl.org/info/docs/tools/vep/vep_formats.html
-   */
-  private static final String ALLELE_FREQUENCY_KEY = "AF";
-
-  /*
    * delimiter that separates multiple consequence data blocks
    */
   private static final String COMMA = ",";
@@ -135,9 +175,9 @@ public class VCFLoader
   private static final String EXCL = "!";
 
   /*
-   * the alignment we are associating VCF data with
+   * the VCF file we are processing
    */
-  private AlignmentI al;
+  protected String vcfFilePath;
 
   /*
    * mappings between VCF and sequence reference assembly regions, as 
@@ -146,12 +186,24 @@ public class VCFLoader
    */
   private Map<String, Map<int[], int[]>> assemblyMappings;
 
+  private VCFReader reader;
+
   /*
    * holds details of the VCF header lines (metadata)
    */
   private VCFHeader header;
 
   /*
+   * species (as a valid Ensembl term) the VCF is for 
+   */
+  private String vcfSpecies;
+
+  /*
+   * genome assembly version (as a valid Ensembl identifier) the VCF is for 
+   */
+  private String vcfAssembly;
+
+  /*
    * a Dictionary of contigs (if present) referenced in the VCF file
    */
   private SAMSequenceDictionary dictionary;
@@ -188,30 +240,42 @@ 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 an alignment context
+   * Constructor given a VCF file
    * 
    * @param alignment
    */
-  public VCFLoader(AlignmentI alignment)
+  public VCFLoader(String vcfFile)
   {
-    al = alignment;
+    try
+    {
+      initialise(vcfFile);
+    } catch (IOException e)
+    {
+      System.err.println("Error opening VCF file: " + e.getMessage());
+    }
 
     // map of species!chromosome!fromAssembly!toAssembly to {fromRange, toRange}
     assemblyMappings = new HashMap<>();
   }
 
   /**
-   * Starts a new thread to query and load VCF variant data on to the alignment
+   * Starts a new thread to query and load VCF variant data on to the given
+   * sequences
    * <p>
    * This method is not thread safe - concurrent threads should use separate
    * instances of this class.
    * 
-   * @param filePath
+   * @param seqs
    * @param gui
    */
-  public void loadVCF(final String filePath,
-          final AlignViewControllerGuiI gui)
+  public void loadVCF(SequenceI[] seqs, final AlignViewControllerGuiI gui)
   {
     if (gui != null)
     {
@@ -220,54 +284,70 @@ public class VCFLoader
 
     new Thread()
     {
-
       @Override
       public void run()
       {
-        VCFLoader.this.doLoad(filePath, gui);
+        VCFLoader.this.doLoad(seqs, gui);
       }
-
     }.start();
   }
 
   /**
-   * Loads VCF on to an alignment - provided it can be related to one or more
-   * sequence's chromosomal coordinates
+   * Reads the specified contig sequence and adds its VCF variants to it
    * 
-   * @param filePath
-   * @param gui
-   *          optional callback handler for messages
+   * @param contig
+   *          the id of a single sequence (contig) to load
+   * @return
    */
-  protected void doLoad(String filePath, AlignViewControllerGuiI gui)
+  public SequenceI loadVCFContig(String contig)
   {
-    VCFReader reader = null;
-    try
+    VCFHeaderLine headerLine = header.getOtherHeaderLine(VCFHeader.REFERENCE_KEY);
+    if (headerLine == null)
     {
-      // long start = System.currentTimeMillis();
-      reader = new VCFReader(filePath);
-
-      header = reader.getFileHeader();
-
-      try
-      {
-        dictionary = header.getSequenceDictionary();
-      } catch (SAMException e)
-      {
-        // ignore - thrown if any contig line lacks length info
-      }
+      Cache.log.error("VCF reference header not found");
+      return null;
+    }
+    String ref = headerLine.getValue();
+    if (ref.startsWith("file://"))
+    {
+      ref = ref.substring(7);
+    }
+    setSpeciesAndAssembly(ref);
 
-      sourceId = filePath;
+    SequenceI seq = null;
+    File dbFile = new File(ref);
 
-      saveMetadata(sourceId);
+    if (dbFile.exists())
+    {
+      HtsContigDb db = new HtsContigDb("", dbFile);
+      seq = db.getSequenceProxy(contig);
+      loadSequenceVCF(seq);
+      db.close();
+    }
+    else
+    {
+      Cache.log.error("VCF reference not found: " + ref);
+    }
 
-      /*
-       * get offset of CSQ ALLELE_NUM and Feature if declared
-       */
-      parseCsqHeader();
+    return seq;
+  }
 
+  /**
+   * Loads VCF on to one or more sequences
+   * 
+   * @param seqs
+   * @param gui
+   *          optional callback handler for messages
+   */
+  protected void doLoad(SequenceI[] seqs, AlignViewControllerGuiI gui)
+  {
+    try
+    {
       VCFHeaderLine ref = header
               .getOtherHeaderLine(VCFHeader.REFERENCE_KEY);
-      String vcfAssembly = ref.getValue();
+      String reference = ref == null ? null : ref.getValue();
+
+      setSpeciesAndAssembly(reference);
 
       int varCount = 0;
       int seqCount = 0;
@@ -275,9 +355,9 @@ public class VCFLoader
       /*
        * query for VCF overlapping each sequence in turn
        */
-      for (SequenceI seq : al.getSequences())
+      for (SequenceI seq : seqs)
       {
-        int added = loadSequenceVCF(seq, reader, vcfAssembly);
+        int added = loadSequenceVCF(seq);
         if (added > 0)
         {
           seqCount++;
@@ -287,7 +367,6 @@ public class VCFLoader
       }
       if (gui != null)
       {
-        // long elapsed = System.currentTimeMillis() - start;
         String msg = MessageManager.formatMessage("label.added_vcf",
                 varCount, seqCount);
         gui.setStatus(msg);
@@ -322,6 +401,103 @@ public class VCFLoader
   }
 
   /**
+   * Attempts to determine and save the species and genome assembly version to
+   * which the VCF data applies. This may be done by parsing the {@code reference}
+   * header line, configured in a property file, or (potentially) confirmed
+   * interactively by the user.
+   * <p>
+   * The saved values should be identifiers valid for Ensembl's REST service
+   * {@code map} endpoint, so they can be used (if necessary) to retrieve the
+   * mapping between VCF coordinates and sequence coordinates.
+   * 
+   * @param reference
+   * @see https://rest.ensembl.org/documentation/info/assembly_map
+   * @see https://rest.ensembl.org/info/assembly/human?content-type=text/xml
+   * @see https://rest.ensembl.org/info/species?content-type=text/xml
+   */
+  protected void setSpeciesAndAssembly(String reference)
+  {
+    if (reference == null)
+    {
+      Cache.log.error("No VCF ##reference found, defaulting to "
+              + DEFAULT_REFERENCE + ":" + DEFAULT_SPECIES);
+      reference = DEFAULT_REFERENCE; // default to GRCh37 if not specified
+    }
+    reference = reference.toLowerCase();
+
+    /*
+     * for a non-human species, or other assembly identifier,
+     * specify as a Jalview property file entry e.g.
+     * VCF_ASSEMBLY = hs37=GRCh37,assembly19=GRCh37
+     * VCF_SPECIES = c_elegans=celegans
+     * to map a token in the reference header to a value
+     */
+    String prop = Cache.getDefault(VCF_ASSEMBLY, DEFAULT_VCF_ASSEMBLY);
+    for (String token : prop.split(","))
+    {
+      String[] tokens = token.split("=");
+      if (tokens.length == 2)
+      {
+        if (reference.contains(tokens[0].trim().toLowerCase()))
+        {
+          vcfAssembly = tokens[1].trim();
+          break;
+        }
+      }
+    }
+
+    vcfSpecies = DEFAULT_SPECIES;
+    prop = Cache.getProperty(VCF_SPECIES);
+    if (prop != null)
+    {
+      for (String token : prop.split(","))
+      {
+        String[] tokens = token.split("=");
+        if (tokens.length == 2)
+        {
+          if (reference.contains(tokens[0].trim().toLowerCase()))
+          {
+            vcfSpecies = tokens[1].trim();
+            break;
+          }
+        }
+      }
+    }
+  }
+
+  /**
+   * Opens the VCF file and parses header data
+   * 
+   * @param filePath
+   * @throws IOException
+   */
+  private void initialise(String filePath) throws IOException
+  {
+    vcfFilePath = filePath;
+
+    reader = new VCFReader(filePath);
+
+    header = reader.getFileHeader();
+
+    try
+    {
+      dictionary = header.getSequenceDictionary();
+    } catch (SAMException e)
+    {
+      // ignore - thrown if any contig line lacks length info
+    }
+
+    sourceId = filePath;
+
+    saveMetadata(sourceId);
+
+    /*
+     * get offset of CSQ ALLELE_NUM and Feature if declared
+     */
+    parseCsqHeader();
+  }
+
+  /**
    * Reads metadata (such as INFO field descriptions and datatypes) and saves
    * them for future reference
    * 
@@ -515,7 +691,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
       {
@@ -536,21 +713,15 @@ public class VCFLoader
   }
 
   /**
-   * Tries to add overlapping variants read from a VCF file to the given
-   * sequence, and returns the number of variant features added. Note that this
-   * requires the sequence to hold information as to its species, chromosomal
-   * positions and reference assembly, in order to be able to map the VCF
-   * variants to the sequence (or not)
+   * Tries to add overlapping variants read from a VCF file to the given sequence,
+   * and returns the number of variant features added
    * 
    * @param seq
-   * @param reader
-   * @param vcfAssembly
    * @return
    */
-  protected int loadSequenceVCF(SequenceI seq, VCFReader reader,
-          String vcfAssembly)
+  protected int loadSequenceVCF(SequenceI seq)
   {
-    VCFMap vcfMap = getVcfMap(seq, vcfAssembly);
+    VCFMap vcfMap = getVcfMap(seq);
     if (vcfMap == null)
     {
       return 0;
@@ -564,17 +735,16 @@ public class VCFLoader
     {
       dss = seq;
     }
-    return addVcfVariants(dss, reader, vcfMap, vcfAssembly);
+    return addVcfVariants(dss, vcfMap);
   }
 
   /**
    * Answers a map from sequence coordinates to VCF chromosome ranges
    * 
    * @param seq
-   * @param vcfAssembly
    * @return
    */
-  private VCFMap getVcfMap(SequenceI seq, String vcfAssembly)
+  private VCFMap getVcfMap(SequenceI seq)
   {
     /*
      * simplest case: sequence has id and length matching a VCF contig
@@ -605,34 +775,28 @@ public class VCFLoader
     String species = seqCoords.getSpeciesId();
     String chromosome = seqCoords.getChromosomeId();
     String seqRef = seqCoords.getAssemblyId();
-    MapList map = seqCoords.getMap();
+    MapList map = seqCoords.getMapping();
 
-    if (!vcfSpeciesMatchesSequence(vcfAssembly, species))
+    // note this requires the configured species to match that
+    // returned with the Ensembl sequence; todo: support aliases?
+    if (!vcfSpecies.equalsIgnoreCase(species))
     {
+      Cache.log.warn("No VCF loaded to " + seq.getName()
+              + " as species not matched");
       return null;
     }
 
-    if (vcfAssemblyMatchesSequence(vcfAssembly, seqRef))
+    if (seqRef.equalsIgnoreCase(vcfAssembly))
     {
       return new VCFMap(chromosome, map);
     }
 
-    if (!"GRCh38".equalsIgnoreCase(seqRef) // Ensembl
-            || !vcfAssembly.contains("Homo_sapiens_assembly19")) // gnomAD
-    {
-      return null;
-    }
-
     /*
-     * map chromosomal coordinates from sequence to VCF if the VCF
-     * data has a different reference assembly to the sequence
+     * VCF data has a different reference assembly to the sequence:
+     * query Ensembl to map chromosomal coordinates from sequence to VCF
      */
-    // TODO generalise for cases other than GRCh38 -> GRCh37 !
-    // - or get the user to choose in a dialog
-
     List<int[]> toVcfRanges = new ArrayList<>();
     List<int[]> fromSequenceRanges = new ArrayList<>();
-    String toRef = "GRCh37";
 
     for (int[] range : map.getToRanges())
     {
@@ -644,12 +808,13 @@ public class VCFLoader
       }
 
       int[] newRange = mapReferenceRange(range, chromosome, "human", seqRef,
-              toRef);
+              vcfAssembly);
       if (newRange == null)
       {
         Cache.log.error(
                 String.format("Failed to map %s:%s:%s:%d:%d to %s", species,
-                        chromosome, seqRef, range[0], range[1], toRef));
+                        chromosome, seqRef, range[0], range[1],
+                        vcfAssembly));
         continue;
       }
       else
@@ -690,76 +855,16 @@ public class VCFLoader
   }
 
   /**
-   * Answers true if we determine that the VCF data uses the same reference
-   * assembly as the sequence, else false
-   * 
-   * @param vcfAssembly
-   * @param seqRef
-   * @return
-   */
-  private boolean vcfAssemblyMatchesSequence(String vcfAssembly,
-          String seqRef)
-  {
-    // TODO improve on this stub, which handles gnomAD and
-    // hopes for the best for other cases
-
-    if ("GRCh38".equalsIgnoreCase(seqRef) // Ensembl
-            && vcfAssembly.contains("Homo_sapiens_assembly19")) // gnomAD
-    {
-      return false;
-    }
-    return true;
-  }
-
-  /**
-   * Answers true if the species inferred from the VCF reference identifier
-   * matches that for the sequence
-   * 
-   * @param vcfAssembly
-   * @param speciesId
-   * @return
-   */
-  boolean vcfSpeciesMatchesSequence(String vcfAssembly, String speciesId)
-  {
-    // PROBLEM 1
-    // there are many aliases for species - how to equate one with another?
-    // PROBLEM 2
-    // VCF ##reference header is an unstructured URI - how to extract species?
-    // perhaps check if ref includes any (Ensembl) alias of speciesId??
-    // TODO ask the user to confirm this??
-
-    if (vcfAssembly.contains("Homo_sapiens") // gnomAD exome data example
-            && "HOMO_SAPIENS".equals(speciesId)) // Ensembl species id
-    {
-      return true;
-    }
-
-    if (vcfAssembly.contains("c_elegans") // VEP VCF response example
-            && "CAENORHABDITIS_ELEGANS".equals(speciesId)) // Ensembl
-    {
-      return true;
-    }
-
-    // this is not a sustainable solution...
-
-    return false;
-  }
-
-  /**
    * Queries the VCF reader for any variants that overlap the mapped chromosome
    * ranges of the sequence, and adds as variant features. Returns the number of
    * overlapping variants found.
    * 
    * @param seq
-   * @param reader
    * @param map
    *          mapping from sequence to VCF coordinates
-   * @param vcfAssembly
-   *          the '##reference' identifier for the VCF reference assembly
    * @return
    */
-  protected int addVcfVariants(SequenceI seq, VCFReader reader,
-          VCFMap map, String vcfAssembly)
+  protected int addVcfVariants(SequenceI seq, VCFMap map)
   {
     boolean forwardStrand = map.map.isToForwardStrand();
 
@@ -772,54 +877,38 @@ 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();
-
-        int[] featureRange = map.map.locateInFrom(variant.getStart(),
-                variant.getEnd());
-
-        if (featureRange != null)
+        CloseableIterator<VariantContext> variants = reader
+                .query(map.chromosome, vcfStart, vcfEnd);
+        while (variants.hasNext())
         {
-          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();
-    }
+          VariantContext variant = variants.next();
 
-    return count;
-  }
+          int[] featureRange = map.map.locateInFrom(variant.getStart(),
+                  variant.getEnd());
 
-  /**
-   * A convenience method to get the AF value for the given alternate allele
-   * index
-   * 
-   * @param variant
-   * @param alleleIndex
-   * @return
-   */
-  protected float getAlleleFrequency(VariantContext variant, int alleleIndex)
-  {
-    float score = 0f;
-    String attributeValue = getAttributeValue(variant,
-            ALLELE_FREQUENCY_KEY, alleleIndex);
-    if (attributeValue != null)
-    {
-      try
-      {
-        score = Float.parseFloat(attributeValue);
-      } catch (NumberFormatException e)
+          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)
       {
-        // leave as 0
+        /*
+         * RuntimeException throwable by htsjdk
+         */
+        String msg = String.format("Error reading VCF for %s:%d-%d: %s ",
+                map.chromosome, vcfStart, vcfEnd);
+        Cache.log.error(msg);
       }
     }
 
-    return score;
+    return count;
   }
 
   /**
@@ -939,19 +1028,29 @@ public class VCFLoader
     String type = SequenceOntologyI.SEQUENCE_VARIANT;
     if (consequence != null)
     {
-      type = getOntologyTerm(seq, variant, altAlleleIndex,
-            consequence);
+      type = getOntologyTerm(consequence);
     }
 
-    float score = getAlleleFrequency(variant, altAlleleIndex);
-
     SequenceFeature sf = new SequenceFeature(type, alleles, featureStart,
-            featureEnd, score, FEATURE_GROUP_VCF);
+            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, seq, sf, altAlleleIndex, consequence);
+    addAlleleProperties(variant, sf, altAlleleIndex, consequence);
 
     seq.addSequenceFeature(sf);
 
@@ -959,6 +1058,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:
@@ -967,18 +1113,18 @@ public class VCFLoader
    * <li>sequence id can be matched to VEP Feature (or SnpEff Feature_ID)</li>
    * </ul>
    * 
-   * @param seq
-   * @param variant
-   * @param altAlleleIndex
    * @param consequence
    * @return
    * @see http://www.sequenceontology.org/browser/current_svn/term/SO:0001060
    */
-  String getOntologyTerm(SequenceI seq, VariantContext variant,
-          int altAlleleIndex, String consequence)
+  String getOntologyTerm(String consequence)
   {
     String type = SequenceOntologyI.SEQUENCE_VARIANT;
 
+    /*
+     * could we associate Consequence data with this allele and feature (transcript)?
+     * if so, prefer the consequence term from that data
+     */
     if (csqAlleleFieldIndex == -1) // && snpEffAlleleFieldIndex == -1
     {
       /*
@@ -987,10 +1133,6 @@ public class VCFLoader
       return type;
     }
 
-    /*
-     * can we associate Consequence data with this allele and feature (transcript)?
-     * if so, prefer the consequence term from that data
-     */
     if (consequence != null)
     {
       String[] csqFields = consequence.split(PIPE_REGEX);
@@ -1121,7 +1263,6 @@ public class VCFLoader
    * Add any allele-specific VCF key-value data to the sequence feature
    * 
    * @param variant
-   * @param seq
    * @param sf
    * @param altAlelleIndex
    *          (0, 1..)
@@ -1129,7 +1270,7 @@ public class VCFLoader
    *          if not null, the consequence specific to this sequence (transcript
    *          feature) and allele
    */
-  protected void addAlleleProperties(VariantContext variant, SequenceI seq,
+  protected void addAlleleProperties(VariantContext variant,
           SequenceFeature sf, final int altAlelleIndex, String consequence)
   {
     Map<String, Object> atts = variant.getAttributes();
@@ -1144,7 +1285,7 @@ public class VCFLoader
        */
       if (CSQ_FIELD.equals(key))
       {
-        addConsequences(variant, seq, sf, consequence);
+        addConsequences(variant, sf, consequence);
         continue;
       }
 
@@ -1193,14 +1334,110 @@ public class VCFLoader
        * take the index'th value
        */
       String value = getAttributeValue(variant, key, index);
-      if (value != null)
+      if (value != null && isValid(variant, key, value))
       {
-        sf.setValue(key, value);
+        value = decodeSpecialCharacters(value);
+        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.
+   * 
+   * @param value
+   * @return
+   */
+  protected static String decodeSpecialCharacters(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))
+    {
+      return value;
+    }
+
+    value = value.replace(ENCODED_COLON, ":")
+            .replace(ENCODED_SEMICOLON, ";").replace(ENCODED_EQUALS, "=")
+            .replace(ENCODED_PERCENT, "%").replace(ENCODED_COMMA, ",");
+    return 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)
+    {
+      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;
+  }
+
+  /**
+   * 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));
+    }
+  }
+
+  /**
    * Inspects CSQ data blocks (consequences) and adds attributes on the sequence
    * feature.
    * <p>
@@ -1209,15 +1446,13 @@ public class VCFLoader
    * transcript (sequence) being processed)
    * 
    * @param variant
-   * @param seq
    * @param sf
    * @param myConsequence
    */
-  protected void addConsequences(VariantContext variant, SequenceI seq,
-          SequenceFeature sf, String myConsequence)
+  protected void addConsequences(VariantContext variant, SequenceFeature sf,
+          String myConsequence)
   {
     Object value = variant.getAttribute(CSQ_FIELD);
-    // TODO if CSQ not present, try ANN (for SnpEff consequence data)?
 
     if (value == null || !(value instanceof List<?>))
     {
@@ -1250,6 +1485,16 @@ 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
+               */
+              try
+              {
+                field = URLDecoder.decode(field, UTF_8);
+              } catch (UnsupportedEncodingException e)
+              {
+              }
               csqValues.put(id, field);
             }
           }