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;
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;
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";
/**
*/
Map<Integer, String> vepFieldsOfInterest;
+ /*
+ * key:value for which rejected data has been seen
+ * (the error is logged only once for each combination)
+ */
+ private Set<String> badData;
+
/**
* Constructor given a VCF file
*
{
int vcfStart = Math.min(range[0], range[1]);
int vcfEnd = Math.max(range[0], range[1]);
- CloseableIterator<VariantContext> variants = reader
- .query(map.chromosome, vcfStart, vcfEnd);
- while (variants.hasNext())
+ try
{
- VariantContext variant = variants.next();
+ CloseableIterator<VariantContext> variants = reader
+ .query(map.chromosome, vcfStart, vcfEnd);
+ while (variants.hasNext())
+ {
+ VariantContext variant = variants.next();
- int[] featureRange = map.map.locateInFrom(variant.getStart(),
- variant.getEnd());
+ int[] featureRange = map.map.locateInFrom(variant.getStart(),
+ variant.getEnd());
- if (featureRange != null)
- {
- int featureStart = Math.min(featureRange[0], featureRange[1]);
- int featureEnd = Math.max(featureRange[0], featureRange[1]);
- count += addAlleleFeatures(seq, variant, featureStart, featureEnd,
- forwardStrand);
+ if (featureRange != null)
+ {
+ int featureStart = Math.min(featureRange[0], featureRange[1]);
+ int featureEnd = Math.max(featureRange[0], featureRange[1]);
+ count += addAlleleFeatures(seq, variant, featureStart,
+ featureEnd, forwardStrand);
+ }
}
+ variants.close();
+ } catch (TribbleException e)
+ {
+ /*
+ * RuntimeException throwable by htsjdk
+ */
+ String msg = String.format("Error reading VCF for %s:%d-%d: %s ",
+ map.chromosome, vcfStart, vcfEnd);
+ Cache.log.error(msg);
}
- variants.close();
}
return count;
featureEnd, FEATURE_GROUP_VCF);
sf.setSource(sourceId);
- sf.setValue(Gff3Helper.ALLELES, alleles);
+ /*
+ * save the derived alleles as a named attribute; this will be
+ * needed when Jalview computes derived peptide variants
+ */
+ addFeatureAttribute(sf, Gff3Helper.ALLELES, alleles);
+
+ /*
+ * add selected VCF fixed column data as feature attributes
+ */
+ addFeatureAttribute(sf, VCF_POS, String.valueOf(variant.getStart()));
+ addFeatureAttribute(sf, VCF_ID, variant.getID());
+ addFeatureAttribute(sf, VCF_QUAL,
+ String.valueOf(variant.getPhredScaledQual()));
+ addFeatureAttribute(sf, VCF_FILTER, getFilter(variant));
addAlleleProperties(variant, sf, altAlleleIndex, consequence);
}
/**
+ * 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:
* take the index'th value
*/
String value = getAttributeValue(variant, key, index);
- if (value != null)
+ if (value != null && isValid(variant, key, value))
{
value = decodeSpecialCharacters(value);
- sf.setValue(key, value);
+ addFeatureAttribute(sf, key, 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>
package jalview.io.vcf;
+import static jalview.io.gff.SequenceOntologyI.SEQUENCE_VARIANT;
import static org.testng.Assert.assertEquals;
+import static org.testng.Assert.assertNull;
import static org.testng.Assert.assertSame;
import static org.testng.Assert.assertTrue;
import jalview.datamodel.Sequence;
import jalview.datamodel.SequenceFeature;
import jalview.datamodel.SequenceI;
+import jalview.datamodel.features.FeatureAttributes;
import jalview.datamodel.features.SequenceFeatures;
import jalview.gui.AlignFrame;
import jalview.io.DataSourceType;
import jalview.io.FileLoader;
import jalview.io.gff.Gff3Helper;
-import jalview.io.gff.SequenceOntologyI;
import jalview.util.MapList;
import java.io.File;
import java.util.Map;
import org.testng.annotations.BeforeClass;
+import org.testng.annotations.BeforeTest;
import org.testng.annotations.Test;
public class VCFLoaderTest
+ ">transcript4/1-18\n-----TGG-GGACGAGAGTGTGA-A\n";
private static final String[] VCF = { "##fileformat=VCFv4.2",
+ // fields other than AF are ignored when parsing as they have no INFO definition
"##INFO=<ID=AF,Number=A,Type=Float,Description=\"Allele Frequency, for each ALT allele, in the same order as listed\">",
+ "##INFO=<ID=AC_Female,Number=A,Type=Integer,Description=\"Allele count in Female genotypes\"",
+ "##INFO=<ID=AF_AFR,Number=A,Type=Float,Description=\"Allele Frequency among African/African American genotypes\"",
"##reference=Homo_sapiens/GRCh38",
"#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO",
// A/T,C variants in position 2 of gene sequence (precedes transcript)
- // should create 2 variant features with respective scores
- "17\t45051611\t.\tA\tT,C\t1666.64\tRF\tAC=15;AF=5.0e-03,4.0e-03",
+ // should create 2 variant features with respective AF values
+ // malformed values for AC_Female and AF_AFR should be ignored
+ "17\t45051611\trs384765\tA\tT,C\t1666.64\tRF;XYZ\tAC=15;AF=5.0e-03,4.0e-03;AC_Female=12,3d;AF_AFR=low,2.3e-4",
// SNP G/C in position 4 of gene sequence, position 2 of transcript
// insertion G/GA is transferred to nucleotide but not to peptide
- "17\t45051613\t.\tG\tGA,C\t1666.64\tRF\tAC=15;AF=3.0e-03,2.0e-03" };
+ "17\t45051613\t.\tG\tGA,C\t1666.65\t.\tAC=15;AF=3.0e-03,2.0e-03",
+ // '.' in INFO field should be ignored
+ "17\t45051615\t.\tG\tC\t1666.66\tRF\tAC=16;AF=." };
@BeforeClass(alwaysRun = true)
public void setUp()
Cache.initLogger();
}
+ @BeforeTest(alwaysRun = true)
+ public void setUpBeforeTest()
+ {
+ /*
+ * clear down feature attributes metadata
+ */
+ FeatureAttributes.getInstance().clear();
+ }
+
@Test(groups = "Functional")
public void testDoLoad() throws IOException
{
AlignmentI al = buildAlignment();
- File f = makeVcf();
+ File f = makeVcfFile();
VCFLoader loader = new VCFLoader(f.getPath());
loader.doLoad(al.getSequencesArray(), null);
List<SequenceFeature> geneFeatures = al.getSequenceAt(0)
.getSequenceFeatures();
SequenceFeatures.sortFeatures(geneFeatures, true);
- assertEquals(geneFeatures.size(), 4);
+ assertEquals(geneFeatures.size(), 5);
SequenceFeature sf = geneFeatures.get(0);
assertEquals(sf.getFeatureGroup(), "VCF");
assertEquals(sf.getBegin(), 2);
assertEquals(sf.getEnd(), 2);
assertEquals(sf.getScore(), 0f);
- assertEquals(Float.parseFloat((String) sf.getValue("AF")), 4.0e-03,
- DELTA);
+ assertEquals(sf.getValue("AF"), "4.0e-03");
+ assertEquals(sf.getValue("AF_AFR"), "2.3e-4");
assertEquals(sf.getValue(Gff3Helper.ALLELES), "A,C");
- assertEquals(sf.getType(), SequenceOntologyI.SEQUENCE_VARIANT);
+ assertEquals(sf.getType(), SEQUENCE_VARIANT);
+ assertEquals(sf.getValue("POS"), "45051611");
+ assertEquals(sf.getValue("ID"), "rs384765");
+ assertEquals(sf.getValue("QUAL"), "1666.64");
+ assertEquals(sf.getValue("FILTER"), "RF;XYZ");
+ // malformed integer for AC_Female is ignored (JAL-3375)
+ assertNull(sf.getValue("AC_Female"));
+
sf = geneFeatures.get(1);
assertEquals(sf.getFeatureGroup(), "VCF");
assertEquals(sf.getBegin(), 2);
assertEquals(sf.getEnd(), 2);
- assertEquals(sf.getType(), SequenceOntologyI.SEQUENCE_VARIANT);
+ assertEquals(sf.getType(), SEQUENCE_VARIANT);
assertEquals(sf.getScore(), 0f);
assertEquals(Float.parseFloat((String) sf.getValue("AF")), 5.0e-03,
DELTA);
+ assertEquals(sf.getValue("AC_Female"), "12");
+ // malformed float for AF_AFR is ignored (JAL-3375)
+ assertNull(sf.getValue("AC_AFR"));
assertEquals(sf.getValue(Gff3Helper.ALLELES), "A,T");
sf = geneFeatures.get(2);
assertEquals(sf.getFeatureGroup(), "VCF");
assertEquals(sf.getBegin(), 4);
assertEquals(sf.getEnd(), 4);
- assertEquals(sf.getType(), SequenceOntologyI.SEQUENCE_VARIANT);
+ assertEquals(sf.getType(), SEQUENCE_VARIANT);
assertEquals(sf.getScore(), 0f);
assertEquals(Float.parseFloat((String) sf.getValue("AF")), 2.0e-03,
DELTA);
assertEquals(sf.getFeatureGroup(), "VCF");
assertEquals(sf.getBegin(), 4);
assertEquals(sf.getEnd(), 4);
- assertEquals(sf.getType(), SequenceOntologyI.SEQUENCE_VARIANT);
+ assertEquals(sf.getType(), SEQUENCE_VARIANT);
assertEquals(sf.getScore(), 0f);
assertEquals(Float.parseFloat((String) sf.getValue("AF")), 3.0e-03,
DELTA);
assertEquals(sf.getValue(Gff3Helper.ALLELES), "G,GA");
+ assertNull(sf.getValue("ID")); // '.' is ignored
+ assertNull(sf.getValue("FILTER")); // '.' is ignored
+
+ sf = geneFeatures.get(4);
+ assertEquals(sf.getFeatureGroup(), "VCF");
+ assertEquals(sf.getBegin(), 6);
+ assertEquals(sf.getEnd(), 6);
+ assertEquals(sf.getType(), SEQUENCE_VARIANT);
+ assertEquals(sf.getScore(), 0f);
+ // AF=. should not have been captured
+ assertNull(sf.getValue("AF"));
+ assertEquals(sf.getValue(Gff3Helper.ALLELES), "G,C");
/*
* verify variant feature(s) added to transcript
*/
List<SequenceFeature> transcriptFeatures = al.getSequenceAt(1)
.getSequenceFeatures();
- assertEquals(transcriptFeatures.size(), 2);
+ assertEquals(transcriptFeatures.size(), 3);
sf = transcriptFeatures.get(0);
assertEquals(sf.getFeatureGroup(), "VCF");
assertEquals(sf.getBegin(), 2);
assertEquals(sf.getEnd(), 2);
- assertEquals(sf.getType(), SequenceOntologyI.SEQUENCE_VARIANT);
+ assertEquals(sf.getType(), SEQUENCE_VARIANT);
assertEquals(sf.getScore(), 0f);
assertEquals(Float.parseFloat((String) sf.getValue("AF")), 2.0e-03,
DELTA);
assertEquals(sf.getFeatureGroup(), "VCF");
assertEquals(sf.getBegin(), 2);
assertEquals(sf.getEnd(), 2);
- assertEquals(sf.getType(), SequenceOntologyI.SEQUENCE_VARIANT);
+ assertEquals(sf.getType(), SEQUENCE_VARIANT);
assertEquals(sf.getScore(), 0f);
assertEquals(Float.parseFloat((String) sf.getValue("AF")), 3.0e-03,
DELTA);
* JAL-3187 don't precompute protein features, do dynamically instead
*/
assertTrue(proteinFeatures.isEmpty());
- // assertEquals(proteinFeatures.size(), 1);
- // sf = proteinFeatures.get(0);
- // assertEquals(sf.getFeatureGroup(), "VCF");
- // assertEquals(sf.getBegin(), 1);
- // assertEquals(sf.getEnd(), 1);
- // assertEquals(sf.getType(), SequenceOntologyI.NONSYNONYMOUS_VARIANT);
- // assertEquals(sf.getDescription(), "p.Ser1Thr");
}
- private File makeVcf() throws IOException
+ private File makeVcfFile() throws IOException
{
File f = File.createTempFile("Test", ".vcf");
f.deleteOnExit();
{
AlignmentI al = buildAlignment();
- File f = makeVcf();
+ File f = makeVcfFile();
VCFLoader loader = new VCFLoader(f.getPath());
List<SequenceFeature> geneFeatures = al.getSequenceAt(2)
.getSequenceFeatures();
SequenceFeatures.sortFeatures(geneFeatures, true);
- assertEquals(geneFeatures.size(), 4);
+ assertEquals(geneFeatures.size(), 5);
+ SequenceFeature sf;
/*
- * variant A/T at 45051611 maps to T/A at gene position 24
+ * insertion G/GA at 45051613 maps to an insertion at
+ * the preceding position (21) on reverse strand gene
+ * reference: CAAGC -> GCTTG/21-25
+ * genomic variant: CAAGAC (G/GA)
+ * gene variant: GTCTTG (G/GT at 21)
*/
- SequenceFeature sf = geneFeatures.get(3);
+ sf = geneFeatures.get(1);
assertEquals(sf.getFeatureGroup(), "VCF");
- assertEquals(sf.getBegin(), 24);
- assertEquals(sf.getEnd(), 24);
- assertEquals(sf.getType(), SequenceOntologyI.SEQUENCE_VARIANT);
+ assertEquals(sf.getBegin(), 21);
+ assertEquals(sf.getEnd(), 21);
+ assertEquals(sf.getType(), SEQUENCE_VARIANT);
assertEquals(sf.getScore(), 0f);
- assertEquals(Float.parseFloat((String) sf.getValue("AF")), 5.0e-03,
+ assertEquals(sf.getValue(Gff3Helper.ALLELES), "G,GT");
+ assertEquals(Float.parseFloat((String) sf.getValue("AF")), 3.0e-03,
DELTA);
- assertEquals(sf.getValue(Gff3Helper.ALLELES), "T,A");
/*
- * variant A/C at 45051611 maps to T/G at gene position 24
+ * variant G/C at 45051613 maps to C/G at gene position 22
*/
sf = geneFeatures.get(2);
assertEquals(sf.getFeatureGroup(), "VCF");
- assertEquals(sf.getBegin(), 24);
- assertEquals(sf.getEnd(), 24);
- assertEquals(sf.getType(), SequenceOntologyI.SEQUENCE_VARIANT);
+ assertEquals(sf.getBegin(), 22);
+ assertEquals(sf.getEnd(), 22);
+ assertEquals(sf.getType(), SEQUENCE_VARIANT);
assertEquals(sf.getScore(), 0f);
- assertEquals(Float.parseFloat((String) sf.getValue("AF")), 4.0e-03,
+ assertEquals(sf.getValue(Gff3Helper.ALLELES), "C,G");
+ assertEquals(Float.parseFloat((String) sf.getValue("AF")), 2.0e-03,
DELTA);
- assertEquals(sf.getValue(Gff3Helper.ALLELES), "T,G");
/*
- * variant G/C at 45051613 maps to C/G at gene position 22
+ * variant A/C at 45051611 maps to T/G at gene position 24
*/
- sf = geneFeatures.get(1);
+ sf = geneFeatures.get(3);
assertEquals(sf.getFeatureGroup(), "VCF");
- assertEquals(sf.getBegin(), 22);
- assertEquals(sf.getEnd(), 22);
- assertEquals(sf.getType(), SequenceOntologyI.SEQUENCE_VARIANT);
+ assertEquals(sf.getBegin(), 24);
+ assertEquals(sf.getEnd(), 24);
+ assertEquals(sf.getType(), SEQUENCE_VARIANT);
assertEquals(sf.getScore(), 0f);
- assertEquals(Float.parseFloat((String) sf.getValue("AF")), 2.0e-03,
+ assertEquals(sf.getValue(Gff3Helper.ALLELES), "T,G");
+ assertEquals(Float.parseFloat((String) sf.getValue("AF")), 4.0e-03,
DELTA);
- assertEquals(sf.getValue(Gff3Helper.ALLELES), "C,G");
/*
- * insertion G/GA at 45051613 maps to an insertion at
- * the preceding position (21) on reverse strand gene
- * reference: CAAGC -> GCTTG/21-25
- * genomic variant: CAAGAC (G/GA)
- * gene variant: GTCTTG (G/GT at 21)
+ * variant A/T at 45051611 maps to T/A at gene position 24
*/
- sf = geneFeatures.get(0);
+ sf = geneFeatures.get(4);
assertEquals(sf.getFeatureGroup(), "VCF");
- assertEquals(sf.getBegin(), 21);
- assertEquals(sf.getEnd(), 21);
- assertEquals(sf.getType(), SequenceOntologyI.SEQUENCE_VARIANT);
+ assertEquals(sf.getBegin(), 24);
+ assertEquals(sf.getEnd(), 24);
+ assertEquals(sf.getType(), SEQUENCE_VARIANT);
assertEquals(sf.getScore(), 0f);
- assertEquals(Float.parseFloat((String) sf.getValue("AF")), 3.0e-03,
+ assertEquals(sf.getValue(Gff3Helper.ALLELES), "T,A");
+ assertEquals(Float.parseFloat((String) sf.getValue("AF")), 5.0e-03,
DELTA);
- assertEquals(sf.getValue(Gff3Helper.ALLELES), "G,GT");
/*
- * verify 2 variant features added to transcript2
+ * verify 3 variant features added to transcript2
*/
List<SequenceFeature> transcriptFeatures = al.getSequenceAt(3)
.getSequenceFeatures();
- assertEquals(transcriptFeatures.size(), 2);
+ assertEquals(transcriptFeatures.size(), 3);
/*
* insertion G/GT at position 21 of gene maps to position 16 of transcript
*/
- sf = transcriptFeatures.get(0);
+ sf = transcriptFeatures.get(1);
assertEquals(sf.getFeatureGroup(), "VCF");
assertEquals(sf.getBegin(), 16);
assertEquals(sf.getEnd(), 16);
- assertEquals(sf.getType(), SequenceOntologyI.SEQUENCE_VARIANT);
+ assertEquals(sf.getType(), SEQUENCE_VARIANT);
assertEquals(sf.getScore(), 0f);
+ assertEquals(sf.getValue(Gff3Helper.ALLELES), "G,GT");
assertEquals(Float.parseFloat((String) sf.getValue("AF")), 3.0e-03,
DELTA);
- assertEquals(sf.getValue(Gff3Helper.ALLELES), "G,GT");
/*
* SNP C/G at position 22 of gene maps to position 17 of transcript
*/
- sf = transcriptFeatures.get(1);
+ sf = transcriptFeatures.get(2);
assertEquals(sf.getFeatureGroup(), "VCF");
assertEquals(sf.getBegin(), 17);
assertEquals(sf.getEnd(), 17);
- assertEquals(sf.getType(), SequenceOntologyI.SEQUENCE_VARIANT);
+ assertEquals(sf.getType(), SEQUENCE_VARIANT);
assertEquals(sf.getScore(), 0f);
+ assertEquals(sf.getValue(Gff3Helper.ALLELES), "C,G");
assertEquals(Float.parseFloat((String) sf.getValue("AF")), 2.0e-03,
DELTA);
- assertEquals(sf.getValue(Gff3Helper.ALLELES), "C,G");
/*
* verify variant feature(s) computed and added to protein
}
}
List<SequenceFeature> proteinFeatures = peptide.getSequenceFeatures();
+
/*
* JAL-3187 don't precompute protein features, do dynamically instead
*/
assertTrue(proteinFeatures.isEmpty());
- // assertEquals(proteinFeatures.size(), 1);
- // sf = proteinFeatures.get(0);
- // assertEquals(sf.getFeatureGroup(), "VCF");
- // assertEquals(sf.getBegin(), 6);
- // assertEquals(sf.getEnd(), 6);
- // assertEquals(sf.getType(), SequenceOntologyI.NONSYNONYMOUS_VARIANT);
- // assertEquals(sf.getDescription(), "p.Ala6Gly");
}
/**
git://source.jalview.org / jalview.git / commitdiff