X-Git-Url: http://source.jalview.org/gitweb/?a=blobdiff_plain;f=src%2Fjalview%2Fio%2Fvcf%2FVCFLoader.java;h=ea5b8e06f723f3ac512a903c780508e2158184f5;hb=57738a1f3c19b1c3a00bd3ac5108f8cd0af32f99;hp=ac707d8f01cea046948f2809c07256b0739c0ae2;hpb=25da4ccb5679905ada8b21e4d21fd416c392ab76;p=jalview.git
diff --git a/src/jalview/io/vcf/VCFLoader.java b/src/jalview/io/vcf/VCFLoader.java
index ac707d8..ea5b8e0 100644
--- a/src/jalview/io/vcf/VCFLoader.java
+++ b/src/jalview/io/vcf/VCFLoader.java
@@ -1,29 +1,29 @@
+/*
+ * 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 .
+ * The Jalview Authors are detailed in the 'AUTHORS' file.
+ */
package jalview.io.vcf;
-import jalview.analysis.Dna;
-import jalview.api.AlignViewControllerGuiI;
-import jalview.bin.Cache;
-import jalview.datamodel.DBRefEntry;
-import jalview.datamodel.GeneLociI;
-import jalview.datamodel.Mapping;
-import jalview.datamodel.SequenceFeature;
-import jalview.datamodel.SequenceI;
-import jalview.datamodel.features.FeatureAttributeType;
-import jalview.datamodel.features.FeatureSource;
-import jalview.datamodel.features.FeatureSources;
-import jalview.ext.ensembl.EnsemblMap;
-import jalview.ext.htsjdk.HtsContigDb;
-import jalview.ext.htsjdk.VCFReader;
-import jalview.io.gff.Gff3Helper;
-import jalview.io.gff.SequenceOntologyI;
-import jalview.util.MapList;
-import jalview.util.MappingUtils;
-import jalview.util.MessageManager;
+import java.util.Locale;
import java.io.File;
import java.io.IOException;
-import java.io.UnsupportedEncodingException;
-import java.net.URLDecoder;
import java.util.ArrayList;
import java.util.HashMap;
import java.util.HashSet;
@@ -48,6 +48,27 @@ import htsjdk.variant.vcf.VCFHeaderLine;
import htsjdk.variant.vcf.VCFHeaderLineCount;
import htsjdk.variant.vcf.VCFHeaderLineType;
import htsjdk.variant.vcf.VCFInfoHeaderLine;
+import jalview.analysis.Dna;
+import jalview.api.AlignViewControllerGuiI;
+import jalview.bin.Cache;
+import jalview.bin.Console;
+import jalview.datamodel.DBRefEntry;
+import jalview.datamodel.GeneLociI;
+import jalview.datamodel.Mapping;
+import jalview.datamodel.SequenceFeature;
+import jalview.datamodel.SequenceI;
+import jalview.datamodel.features.FeatureAttributeType;
+import jalview.datamodel.features.FeatureSource;
+import jalview.datamodel.features.FeatureSources;
+import jalview.ext.ensembl.EnsemblMap;
+import jalview.ext.htsjdk.HtsContigDb;
+import jalview.ext.htsjdk.VCFReader;
+import jalview.io.gff.Gff3Helper;
+import jalview.io.gff.SequenceOntologyI;
+import jalview.util.MapList;
+import jalview.util.MappingUtils;
+import jalview.util.MessageManager;
+import jalview.util.StringUtils;
/**
* A class to read VCF data (using the htsjdk) and add variants as sequence
@@ -57,17 +78,7 @@ 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";
+ private static final String VCF_ENCODABLE = ":;=%,";
/*
* Jalview feature attributes for VCF fixed column data
@@ -85,13 +96,14 @@ public class VCFLoader
private static final String DEFAULT_SPECIES = "homo_sapiens";
/**
- * A class to model the mapping from sequence to VCF coordinates. Cases include
+ * A class to model the mapping from sequence to VCF coordinates. Cases
+ * include
*
* - a direct 1:1 mapping where the sequence is one of the VCF contigs
- * - a mapping of sequence to chromosomal coordinates, where sequence and VCF
- * use the same reference assembly
- * - a modified mapping of sequence to chromosomal coordinates, where sequence
- * and VCF use different reference assembles
+ * - a mapping of sequence to chromosomal coordinates, where sequence and
+ * VCF use the same reference assembly
+ * - a modified mapping of sequence to chromosomal coordinates, where
+ * sequence and VCF use different reference assembles
*
*/
class VCFMap
@@ -135,15 +147,20 @@ public class VCFLoader
private static final String VCF_SPECIES = "VCF_SPECIES"; // default is human
- private static final String DEFAULT_REFERENCE = "grch37"; // fallback default is human GRCh37
+ 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 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_ALLELE_NUM_KEY = "ALLELE_NUM"; // 0 (ref),
+ // 1...
+
private static final String CSQ_FEATURE_KEY = "Feature"; // Ensembl stable id
/*
@@ -214,8 +231,11 @@ public class VCFLoader
* 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
@@ -301,10 +321,11 @@ public class VCFLoader
*/
public SequenceI loadVCFContig(String contig)
{
- VCFHeaderLine headerLine = header.getOtherHeaderLine(VCFHeader.REFERENCE_KEY);
+ VCFHeaderLine headerLine = header
+ .getOtherHeaderLine(VCFHeader.REFERENCE_KEY);
if (headerLine == null)
{
- Cache.log.error("VCF reference header not found");
+ Console.error("VCF reference header not found");
return null;
}
String ref = headerLine.getValue();
@@ -326,7 +347,7 @@ public class VCFLoader
}
else
{
- Cache.log.error("VCF reference not found: " + ref);
+ Console.error("VCF reference not found: " + ref);
}
return seq;
@@ -402,9 +423,9 @@ 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.
+ * 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.
*
* 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
@@ -419,11 +440,11 @@ public class VCFLoader
{
if (reference == null)
{
- Cache.log.error("No VCF ##reference found, defaulting to "
+ Console.error("No VCF ##reference found, defaulting to "
+ DEFAULT_REFERENCE + ":" + DEFAULT_SPECIES);
reference = DEFAULT_REFERENCE; // default to GRCh37 if not specified
}
- reference = reference.toLowerCase();
+ reference = reference.toLowerCase(Locale.ROOT);
/*
* for a non-human species, or other assembly identifier,
@@ -438,7 +459,7 @@ public class VCFLoader
String[] tokens = token.split("=");
if (tokens.length == 2)
{
- if (reference.contains(tokens[0].trim().toLowerCase()))
+ if (reference.contains(tokens[0].trim().toLowerCase(Locale.ROOT)))
{
vcfAssembly = tokens[1].trim();
break;
@@ -455,7 +476,7 @@ public class VCFLoader
String[] tokens = token.split("=");
if (tokens.length == 2)
{
- if (reference.contains(tokens[0].trim().toLowerCase()))
+ if (reference.contains(tokens[0].trim().toLowerCase(Locale.ROOT)))
{
vcfSpecies = tokens[1].trim();
break;
@@ -560,7 +581,7 @@ public class VCFLoader
{
for (Pattern p : filters)
{
- if (p.matcher(id.toUpperCase()).matches())
+ if (p.matcher(id.toUpperCase(Locale.ROOT)).matches())
{
return true;
}
@@ -654,7 +675,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);
@@ -665,13 +686,12 @@ public class VCFLoader
/**
* Transfers VCF features to sequences to which this sequence has a mapping.
- * If the mapping is 3:1, computes peptide variants from nucleotide variants.
*
* @param seq
*/
protected void transferAddedFeatures(SequenceI seq)
{
- DBRefEntry[] dbrefs = seq.getDBRefs();
+ List dbrefs = seq.getDBRefs();
if (dbrefs == null)
{
return;
@@ -713,8 +733,8 @@ 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
+ * Tries to add overlapping variants read from a VCF file to the given
+ * sequence, and returns the number of variant features added
*
* @param seq
* @return
@@ -766,7 +786,7 @@ public class VCFLoader
GeneLociI seqCoords = seq.getGeneLoci();
if (seqCoords == null)
{
- Cache.log.warn(String.format(
+ Console.warn(String.format(
"Can't query VCF for %s as chromosome coordinates not known",
seq.getName()));
return null;
@@ -781,7 +801,7 @@ public class VCFLoader
// returned with the Ensembl sequence; todo: support aliases?
if (!vcfSpecies.equalsIgnoreCase(species))
{
- Cache.log.warn("No VCF loaded to " + seq.getName()
+ Console.warn("No VCF loaded to " + seq.getName()
+ " as species not matched");
return null;
}
@@ -811,10 +831,9 @@ public class VCFLoader
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],
- vcfAssembly));
+ Console.error(String.format("Failed to map %s:%s:%s:%d:%d to %s",
+ species, chromosome, seqRef, range[0], range[1],
+ vcfAssembly));
continue;
}
else
@@ -888,12 +907,19 @@ public class VCFLoader
int[] featureRange = map.map.locateInFrom(variant.getStart(),
variant.getEnd());
+ /*
+ * only take features whose range is fully mappable to sequence positions
+ */
if (featureRange != null)
{
int featureStart = Math.min(featureRange[0], featureRange[1]);
int featureEnd = Math.max(featureRange[0], featureRange[1]);
- count += addAlleleFeatures(seq, variant, featureStart,
- featureEnd, forwardStrand);
+ if (featureEnd - featureStart == variant.getEnd()
+ - variant.getStart())
+ {
+ count += addAlleleFeatures(seq, variant, featureStart,
+ featureEnd, forwardStrand);
+ }
}
}
variants.close();
@@ -903,8 +929,8 @@ public class VCFLoader
* RuntimeException throwable by htsjdk
*/
String msg = String.format("Error reading VCF for %s:%d-%d: %s ",
- map.chromosome, vcfStart, vcfEnd);
- Cache.log.error(msg);
+ map.chromosome, vcfStart, vcfEnd, e.getLocalizedMessage());
+ Console.error(msg);
}
}
@@ -968,8 +994,8 @@ public class VCFLoader
/**
* Inspects one allele and attempts to add a variant feature for it to the
* 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).
+ * store in the feature's key-value map. Answers the number of features added
+ * (0 or 1).
*
* @param seq
* @param variant
@@ -1017,10 +1043,10 @@ public class VCFLoader
* 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);
+ String consequence = getConsequenceForAlleleAndFeature(variant,
+ CSQ_FIELD, altAlleleIndex, csqAlleleFieldIndex,
+ csqAlleleNumberFieldIndex,
+ seq.getName().toLowerCase(Locale.ROOT), csqFeatureFieldIndex);
/*
* pick out the ontology term for the consequence type
@@ -1060,8 +1086,8 @@ 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.
+ * passed. htsjdk saves filters as a HashSet, so the order when reassembled
+ * into a list may be different.
*
* @param variant
* @return
@@ -1105,9 +1131,9 @@ public class VCFLoader
}
/**
- * 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:
+ * 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:
*
* - VEP (or SnpEff) Consequence annotation is included in the VCF
* - sequence id can be matched to VEP Feature (or SnpEff Feature_ID)
@@ -1185,8 +1211,7 @@ public class VCFLoader
*/
private String getConsequenceForAlleleAndFeature(VariantContext variant,
String vcfInfoId, int altAlleleIndex, int alleleFieldIndex,
- int alleleNumberFieldIndex,
- String seqName, int featureFieldIndex)
+ int alleleNumberFieldIndex, String seqName, int featureFieldIndex)
{
if (alleleFieldIndex == -1 || featureFieldIndex == -1)
{
@@ -1211,8 +1236,8 @@ public class VCFLoader
if (csqFields.length > featureFieldIndex)
{
String featureIdentifier = csqFields[featureFieldIndex];
- if (featureIdentifier.length() > 4
- && seqName.indexOf(featureIdentifier.toLowerCase()) > -1)
+ if (featureIdentifier.length() > 4 && seqName
+ .indexOf(featureIdentifier.toLowerCase(Locale.ROOT)) > -1)
{
/*
* feature (transcript) matched - now check for allele match
@@ -1336,45 +1361,20 @@ public class VCFLoader
String value = getAttributeValue(variant, key, index);
if (value != null && isValid(variant, key, value))
{
- value = decodeSpecialCharacters(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.
- *
- * @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.
+ * 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
@@ -1391,7 +1391,7 @@ public class VCFLoader
VCFInfoHeaderLine infoHeader = header.getInfoHeaderLine(infoId);
if (infoHeader == null)
{
- Cache.log.error("Field " + infoId + " has no INFO header");
+ Console.error("Field " + infoId + " has no INFO header");
return false;
}
VCFHeaderLineType infoType = infoHeader.getType();
@@ -1432,7 +1432,7 @@ public class VCFLoader
if (!badData.contains(token))
{
badData.add(token);
- Cache.log.error(String.format("Invalid VCF data at %s:%d %s=%s",
+ Console.error(String.format("Invalid VCF data at %s:%d %s=%s",
variant.getContig(), variant.getStart(), infoId, value));
}
}
@@ -1489,12 +1489,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);
}
}
@@ -1604,8 +1599,8 @@ public class VCFLoader
* @param toRef
* @return
*/
- protected int[] findSubsumedRangeMapping(int[] queryRange, String chromosome,
- String species, String fromRef, String toRef)
+ protected int[] findSubsumedRangeMapping(int[] queryRange,
+ String chromosome, String species, String fromRef, String toRef)
{
String key = makeRangesKey(chromosome, species, fromRef, toRef);
if (assemblyMappings.containsKey(key))
@@ -1627,7 +1622,8 @@ public class VCFLoader
*/
int offset = queryRange[0] - fromRange[0];
int mappedRangeFrom = toRange[0] + offset;
- int mappedRangeTo = mappedRangeFrom + (queryRange[1] - queryRange[0]);
+ int mappedRangeTo = mappedRangeFrom
+ + (queryRange[1] - queryRange[0]);
return new int[] { mappedRangeFrom, mappedRangeTo };
}
}
@@ -1650,7 +1646,7 @@ public class VCFLoader
SequenceI targetSequence, MapList mapping)
{
int[] mappedRange = mapping.locateInTo(sf.getBegin(), sf.getEnd());
-
+
if (mappedRange != null)
{
String group = sf.getFeatureGroup();
@@ -1674,7 +1670,6 @@ public class VCFLoader
protected static String makeRangesKey(String chromosome, String species,
String fromRef, String toRef)
{
- return species + EXCL + chromosome + EXCL + fromRef + EXCL
- + toRef;
+ return species + EXCL + chromosome + EXCL + fromRef + EXCL + toRef;
}
}