X-Git-Url: http://source.jalview.org/gitweb/?a=blobdiff_plain;f=src%2Fjalview%2Fio%2Fvcf%2FVCFLoader.java;h=ac707d8f01cea046948f2809c07256b0739c0ae2;hb=ba2c0286e9a47af8ffd46057f91a22df0d4850ea;hp=a85802ae2117001c7a1377c1284334dffbdc43dc;hpb=946a0619168de3fbe9858401d13bcd1981d0a0df;p=jalview.git

diff --git a/src/jalview/io/vcf/VCFLoader.java b/src/jalview/io/vcf/VCFLoader.java
index a85802a..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,14 +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.SortedMap;
-import java.util.TreeMap;
+import java.util.Set;
 import java.util.regex.Pattern;
 import java.util.regex.PatternSyntaxException;
 
@@ -36,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;
@@ -52,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>
@@ -73,11 +105,17 @@ public class VCFLoader
       chromosome = chr;
       map = m;
     }
+
+    @Override
+    public String toString()
+    {
+      return chromosome + ":" + map.toString();
+    }
   }
 
   /*
    * 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";
 
@@ -88,13 +126,25 @@ 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
    */
-  private static final String ALLELE_KEY = "Allele";
-
-  private static final String ALLELE_NUM_KEY = "ALLELE_NUM"; // 0 (ref), 1...
-  private static final String FEATURE_KEY = "Feature"; // Ensembl stable id
+  private static final String CSQ_CONSEQUENCE_KEY = "Consequence";
+  private static final String CSQ_ALLELE_KEY = "Allele";
+  private static final String CSQ_ALLELE_NUM_KEY = "ALLELE_NUM"; // 0 (ref), 1...
+  private static final String CSQ_FEATURE_KEY = "Feature"; // Ensembl stable id
 
   /*
    * default VCF INFO key for VEP consequence data
@@ -109,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 = ",";
@@ -131,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 
@@ -142,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;
@@ -157,10 +213,14 @@ public class VCFLoader
    * CSQ (consequence) data (if declared in the VCF INFO header for CSQ)
    * see http://www.ensembl.org/info/docs/tools/vep/vep_formats.html
    */
+  private int csqConsequenceFieldIndex = -1;
   private int csqAlleleFieldIndex = -1;
   private int csqAlleleNumberFieldIndex = -1;
   private int csqFeatureFieldIndex = -1;
 
+  // todo the same fields for SnpEff ANN data if wanted
+  // see http://snpeff.sourceforge.net/SnpEff_manual.html#input
+
   /*
    * a unique identifier under which to save metadata about feature
    * attributes (selected INFO field data)
@@ -180,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)
     {
@@ -212,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;
@@ -267,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++;
@@ -279,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);
@@ -314,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
    * 
@@ -419,15 +603,19 @@ public class VCFLoader
       int index = 0;
       for (String field : format)
       {
-        if (ALLELE_NUM_KEY.equals(field))
+        if (CSQ_CONSEQUENCE_KEY.equals(field))
+        {
+          csqConsequenceFieldIndex = index;
+        }
+        if (CSQ_ALLELE_NUM_KEY.equals(field))
         {
           csqAlleleNumberFieldIndex = index;
         }
-        if (ALLELE_KEY.equals(field))
+        if (CSQ_ALLELE_KEY.equals(field))
         {
           csqAlleleFieldIndex = index;
         }
-        if (FEATURE_KEY.equals(field))
+        if (CSQ_FEATURE_KEY.equals(field))
         {
           csqFeatureFieldIndex = index;
         }
@@ -503,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
       {
@@ -524,37 +713,38 @@ 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;
     }
 
-    return addVcfVariants(seq, reader, vcfMap, vcfAssembly);
+    /*
+     * work with the dataset sequence here
+     */
+    SequenceI dss = seq.getDatasetSequence();
+    if (dss == null)
+    {
+      dss = seq;
+    }
+    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
@@ -585,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())
     {
@@ -624,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
@@ -670,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();
 
@@ -752,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;
   }
 
   /**
@@ -858,9 +967,9 @@ public class VCFLoader
 
   /**
    * Inspects one allele and attempts to add a variant feature for it to the
-   * sequence. We extract as much as possible of the additional data associated
-   * with this allele to store in the feature's key-value map. Answers the
-   * number of features added (0 or 1).
+   * sequence. The additional data associated with this allele is extracted to
+   * store in the feature's key-value map. Answers the number of features added (0
+   * or 1).
    * 
    * @param seq
    * @param variant
@@ -880,26 +989,68 @@ public class VCFLoader
     String allele = alt.getBaseString();
 
     /*
+     * insertion after a genomic base, if on reverse strand, has to be 
+     * converted to insertion of complement after the preceding position 
+     */
+    int referenceLength = reference.length();
+    if (!forwardStrand && allele.length() > referenceLength
+            && allele.startsWith(reference))
+    {
+      featureStart -= referenceLength;
+      featureEnd = featureStart;
+      char insertAfter = seq.getCharAt(featureStart - seq.getStart());
+      reference = Dna.reverseComplement(String.valueOf(insertAfter));
+      allele = allele.substring(referenceLength) + reference;
+    }
+
+    /*
      * build the ref,alt allele description e.g. "G,A", using the base
      * complement if the sequence is on the reverse strand
      */
-    // TODO check how structural variants are shown on reverse strand
     StringBuilder sb = new StringBuilder();
     sb.append(forwardStrand ? reference : Dna.reverseComplement(reference));
     sb.append(COMMA);
     sb.append(forwardStrand ? allele : Dna.reverseComplement(allele));
     String alleles = sb.toString(); // e.g. G,A
 
+    /*
+     * pick out the consequence data (if any) that is for the current allele
+     * and feature (transcript) that matches the current sequence
+     */
+    String consequence = getConsequenceForAlleleAndFeature(variant, CSQ_FIELD,
+            altAlleleIndex, csqAlleleFieldIndex,
+            csqAlleleNumberFieldIndex, seq.getName().toLowerCase(),
+            csqFeatureFieldIndex);
+
+    /*
+     * pick out the ontology term for the consequence type
+     */
     String type = SequenceOntologyI.SEQUENCE_VARIANT;
-    float score = getAlleleFrequency(variant, altAlleleIndex);
+    if (consequence != null)
+    {
+      type = getOntologyTerm(consequence);
+    }
 
     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);
 
-    addAlleleProperties(variant, seq, sf, altAlleleIndex);
+    /*
+     * 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);
 
     seq.addSequenceFeature(sf);
 
@@ -907,16 +1058,220 @@ 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:
+   * <ul>
+   * <li>VEP (or SnpEff) Consequence annotation is included in the VCF</li>
+   * <li>sequence id can be matched to VEP Feature (or SnpEff Feature_ID)</li>
+   * </ul>
+   * 
+   * @param consequence
+   * @return
+   * @see http://www.sequenceontology.org/browser/current_svn/term/SO:0001060
+   */
+  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
+    {
+      /*
+       * no Consequence data so we can't refine the ontology term
+       */
+      return type;
+    }
+
+    if (consequence != null)
+    {
+      String[] csqFields = consequence.split(PIPE_REGEX);
+      if (csqFields.length > csqConsequenceFieldIndex)
+      {
+        type = csqFields[csqConsequenceFieldIndex];
+      }
+    }
+    else
+    {
+      // todo the same for SnpEff consequence data matching if wanted
+    }
+
+    /*
+     * if of the form (e.g.) missense_variant&splice_region_variant,
+     * just take the first ('most severe') consequence
+     */
+    if (type != null)
+    {
+      int pos = type.indexOf('&');
+      if (pos > 0)
+      {
+        type = type.substring(0, pos);
+      }
+    }
+    return type;
+  }
+
+  /**
+   * Returns matched consequence data if it can be found, else null.
+   * <ul>
+   * <li>inspects the VCF data for key 'vcfInfoId'</li>
+   * <li>splits this on comma (to distinct consequences)</li>
+   * <li>returns the first consequence (if any) where</li>
+   * <ul>
+   * <li>the allele matches the altAlleleIndex'th allele of variant</li>
+   * <li>the feature matches the sequence name (e.g. transcript id)</li>
+   * </ul>
+   * </ul>
+   * If matched, the consequence is returned (as pipe-delimited fields).
+   * 
+   * @param variant
+   * @param vcfInfoId
+   * @param altAlleleIndex
+   * @param alleleFieldIndex
+   * @param alleleNumberFieldIndex
+   * @param seqName
+   * @param featureFieldIndex
+   * @return
+   */
+  private String getConsequenceForAlleleAndFeature(VariantContext variant,
+          String vcfInfoId, int altAlleleIndex, int alleleFieldIndex,
+          int alleleNumberFieldIndex,
+          String seqName, int featureFieldIndex)
+  {
+    if (alleleFieldIndex == -1 || featureFieldIndex == -1)
+    {
+      return null;
+    }
+    Object value = variant.getAttribute(vcfInfoId);
+
+    if (value == null || !(value instanceof List<?>))
+    {
+      return null;
+    }
+
+    /*
+     * inspect each consequence in turn (comma-separated blocks
+     * extracted by htsjdk)
+     */
+    List<String> consequences = (List<String>) value;
+
+    for (String consequence : consequences)
+    {
+      String[] csqFields = consequence.split(PIPE_REGEX);
+      if (csqFields.length > featureFieldIndex)
+      {
+        String featureIdentifier = csqFields[featureFieldIndex];
+        if (featureIdentifier.length() > 4
+                && seqName.indexOf(featureIdentifier.toLowerCase()) > -1)
+        {
+          /*
+           * feature (transcript) matched - now check for allele match
+           */
+          if (matchAllele(variant, altAlleleIndex, csqFields,
+                  alleleFieldIndex, alleleNumberFieldIndex))
+          {
+            return consequence;
+          }
+        }
+      }
+    }
+    return null;
+  }
+
+  private boolean matchAllele(VariantContext variant, int altAlleleIndex,
+          String[] csqFields, int alleleFieldIndex,
+          int alleleNumberFieldIndex)
+  {
+    /*
+     * if ALLELE_NUM is present, it must match altAlleleIndex
+     * NB first alternate allele is 1 for ALLELE_NUM, 0 for altAlleleIndex
+     */
+    if (alleleNumberFieldIndex > -1)
+    {
+      if (csqFields.length <= alleleNumberFieldIndex)
+      {
+        return false;
+      }
+      String alleleNum = csqFields[alleleNumberFieldIndex];
+      return String.valueOf(altAlleleIndex + 1).equals(alleleNum);
+    }
+
+    /*
+     * else consequence allele must match variant allele
+     */
+    if (alleleFieldIndex > -1 && csqFields.length > alleleFieldIndex)
+    {
+      String csqAllele = csqFields[alleleFieldIndex];
+      String vcfAllele = variant.getAlternateAllele(altAlleleIndex)
+              .getBaseString();
+      return csqAllele.equals(vcfAllele);
+    }
+    return false;
+  }
+
+  /**
    * Add any allele-specific VCF key-value data to the sequence feature
    * 
    * @param variant
-   * @param seq
    * @param sf
    * @param altAlelleIndex
    *          (0, 1..)
+   * @param consequence
+   *          if not null, the consequence specific to this sequence (transcript
+   *          feature) and allele
    */
-  protected void addAlleleProperties(VariantContext variant, SequenceI seq,
-          SequenceFeature sf, final int altAlelleIndex)
+  protected void addAlleleProperties(VariantContext variant,
+          SequenceFeature sf, final int altAlelleIndex, String consequence)
   {
     Map<String, Object> atts = variant.getAttributes();
 
@@ -930,7 +1285,7 @@ public class VCFLoader
        */
       if (CSQ_FIELD.equals(key))
       {
-        addConsequences(variant, seq, sf, altAlelleIndex);
+        addConsequences(variant, sf, consequence);
         continue;
       }
 
@@ -979,39 +1334,127 @@ 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 for the current allele (and transcript if applicable)
-   * <p>
-   * Allele matching: if field ALLELE_NUM is present, it must match
-   * altAlleleIndex. If not present, then field Allele value must match the VCF
-   * Allele.
+   * feature.
    * <p>
-   * Transcript matching: if sequence name can be identified to at least one of
-   * the consequences' Feature values, then select only consequences that match
-   * the value (i.e. consequences for the current transcript sequence). If not,
-   * take all consequences (this is the case when adding features to the gene
-   * sequence).
+   * If <code>myConsequence</code> is not null, then this is the specific
+   * consequence data (pipe-delimited fields) that is for the current allele and
+   * transcript (sequence) being processed)
    * 
    * @param variant
-   * @param seq
    * @param sf
-   * @param altAlelleIndex
-   *          (0, 1..)
+   * @param myConsequence
    */
-  protected void addConsequences(VariantContext variant, SequenceI seq,
-          SequenceFeature sf, int altAlelleIndex)
+  protected void addConsequences(VariantContext variant, SequenceFeature sf,
+          String myConsequence)
   {
     Object value = variant.getAttribute(CSQ_FIELD);
 
-    if (value == null || !(value instanceof ArrayList<?>))
+    if (value == null || !(value instanceof List<?>))
     {
       return;
     }
@@ -1019,43 +1462,17 @@ public class VCFLoader
     List<String> consequences = (List<String>) value;
 
     /*
-     * if CSQ data includes 'Feature', and any value matches the sequence name,
-     * then restrict consequence data to only the matching value (transcript)
-     * i.e. just pick out consequences for the transcript the variant feature is on
-     */
-    String seqName = seq.getName()== null ? "" : seq.getName().toLowerCase();
-    String matchFeature = null;
-    if (csqFeatureFieldIndex > -1)
-    {
-      for (String consequence : consequences)
-      {
-        String[] csqFields = consequence.split(PIPE_REGEX);
-        if (csqFields.length > csqFeatureFieldIndex)
-        {
-          String featureIdentifier = csqFields[csqFeatureFieldIndex];
-          if (featureIdentifier.length() > 4
-                  && seqName.indexOf(featureIdentifier.toLowerCase()) > -1)
-          {
-            matchFeature = featureIdentifier;
-          }
-        }
-      }
-    }
-
-    /*
-     * inspect CSQ consequences; where possible restrict to the consequence
+     * inspect CSQ consequences; restrict to the consequence
      * associated with the current transcript (Feature)
      */
-    SortedMap<String, String> csqValues = new TreeMap<>(
-            String.CASE_INSENSITIVE_ORDER);
+    Map<String, String> csqValues = new HashMap<>();
 
     for (String consequence : consequences)
     {
-      String[] csqFields = consequence.split(PIPE_REGEX);
-
-      if (includeConsequence(csqFields, matchFeature, variant,
-              altAlelleIndex))
+      if (myConsequence == null || myConsequence.equals(consequence))
       {
+        String[] csqFields = consequence.split(PIPE_REGEX);
+
         /*
          * inspect individual fields of this consequence, copying non-null
          * values which are 'fields of interest'
@@ -1068,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);
             }
           }
@@ -1083,72 +1510,6 @@ public class VCFLoader
   }
 
   /**
-   * Answers true if we want to associate this block of consequence data with
-   * the specified alternate allele of the VCF variant.
-   * <p>
-   * If consequence data includes the ALLELE_NUM field, then this has to match
-   * altAlleleIndex. Otherwise the Allele field of the consequence data has to
-   * match the allele value.
-   * <p>
-   * Optionally (if matchFeature is not null), restrict to only include
-   * consequences whose Feature value matches. This allows us to attach
-   * consequences to their respective transcripts.
-   * 
-   * @param csqFields
-   * @param matchFeature
-   * @param variant
-   * @param altAlelleIndex
-   *          (0, 1..)
-   * @return
-   */
-  protected boolean includeConsequence(String[] csqFields,
-          String matchFeature, VariantContext variant, int altAlelleIndex)
-  {
-    /*
-     * check consequence is for the current transcript
-     */
-    if (matchFeature != null)
-    {
-      if (csqFields.length <= csqFeatureFieldIndex)
-      {
-        return false;
-      }
-      String featureIdentifier = csqFields[csqFeatureFieldIndex];
-      if (!featureIdentifier.equals(matchFeature))
-      {
-        return false; // consequence is for a different transcript
-      }
-    }
-
-    /*
-     * if ALLELE_NUM is present, it must match altAlleleIndex
-     * NB first alternate allele is 1 for ALLELE_NUM, 0 for altAlleleIndex
-     */
-    if (csqAlleleNumberFieldIndex > -1)
-    {
-      if (csqFields.length <= csqAlleleNumberFieldIndex)
-      {
-        return false;
-      }
-      String alleleNum = csqFields[csqAlleleNumberFieldIndex];
-      return String.valueOf(altAlelleIndex + 1).equals(alleleNum);
-    }
-
-    /*
-     * else consequence allele must match variant allele
-     */
-    if (csqAlleleFieldIndex > -1 && csqFields.length > csqAlleleFieldIndex)
-    {
-      String csqAllele = csqFields[csqAlleleFieldIndex];
-      String vcfAllele = variant.getAlternateAllele(altAlelleIndex)
-              .getBaseString();
-      return csqAllele.equals(vcfAllele);
-    }
-
-    return false;
-  }
-
-  /**
    * A convenience method to complement a dna base and return the string value
    * of its complement
    *