From b53d89acfa678df63d6870176b4c7ec9285f52ee Mon Sep 17 00:00:00 2001 From: gmungoc Date: Thu, 25 Jul 2019 12:48:33 +0200 Subject: [PATCH] JAL-3375 ignore '.' values for VCF data --- .../datamodel/features/FeatureAttributes.java | 23 +++ src/jalview/io/vcf/VCFLoader.java | 4 +- test/jalview/io/vcf/VCFLoaderTest.java | 150 +++++++++++++------- 3 files changed, 121 insertions(+), 56 deletions(-) diff --git a/src/jalview/datamodel/features/FeatureAttributes.java b/src/jalview/datamodel/features/FeatureAttributes.java index 10249f3..26e6f8b 100644 --- a/src/jalview/datamodel/features/FeatureAttributes.java +++ b/src/jalview/datamodel/features/FeatureAttributes.java @@ -371,4 +371,27 @@ public class FeatureAttributes } return null; } + + /** + * Resets all attribute metadata + */ + public void clear() + { + attributes.clear(); + } + + /** + * Resets attribute metadata for one feature type + * + * @param featureType + */ + public void clear(String featureType) + { + Map map = attributes.get(featureType); + if (map != null) + { + map.clear(); + } + + } } diff --git a/src/jalview/io/vcf/VCFLoader.java b/src/jalview/io/vcf/VCFLoader.java index 7bf7791..053b52f 100644 --- a/src/jalview/io/vcf/VCFLoader.java +++ b/src/jalview/io/vcf/VCFLoader.java @@ -51,6 +51,8 @@ import htsjdk.variant.vcf.VCFInfoHeaderLine; */ public class VCFLoader { + private static final String NO_VALUE = "."; + private static final String DEFAULT_SPECIES = "homo_sapiens"; /** @@ -877,7 +879,7 @@ public class VCFLoader if (att instanceof String) { - return (String) att; + return NO_VALUE.equals(att) ? null : (String) att; } else if (att instanceof ArrayList) { diff --git a/test/jalview/io/vcf/VCFLoaderTest.java b/test/jalview/io/vcf/VCFLoaderTest.java index a87c160..6ba6fbe 100644 --- a/test/jalview/io/vcf/VCFLoaderTest.java +++ b/test/jalview/io/vcf/VCFLoaderTest.java @@ -1,6 +1,9 @@ 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 jalview.bin.Cache; import jalview.datamodel.AlignmentI; @@ -9,6 +12,8 @@ import jalview.datamodel.Mapping; import jalview.datamodel.Sequence; import jalview.datamodel.SequenceFeature; import jalview.datamodel.SequenceI; +import jalview.datamodel.features.FeatureAttributes; +import jalview.datamodel.features.FeatureAttributes.Datatype; import jalview.datamodel.features.SequenceFeatures; import jalview.gui.AlignFrame; import jalview.io.DataSourceType; @@ -24,6 +29,7 @@ import java.util.List; import java.util.Map; import org.testng.annotations.BeforeClass; +import org.testng.annotations.BeforeTest; import org.testng.annotations.Test; public class VCFLoaderTest @@ -56,15 +62,18 @@ 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=", "##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 + // should create 2 variant features with respective AF values "17\t45051611\t.\tA\tT,C\t1666.64\tRF\tAC=15;AF=5.0e-03,4.0e-03", // 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\tRF\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() @@ -79,12 +88,21 @@ public class VCFLoaderTest 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); @@ -99,7 +117,7 @@ public class VCFLoaderTest List 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); @@ -108,12 +126,12 @@ public class VCFLoaderTest assertEquals(Float.parseFloat((String) sf.getValue("AF")), 4.0e-03, DELTA); assertEquals(sf.getValue(Gff3Helper.ALLELES), "A,C"); - assertEquals(sf.getType(), SequenceOntologyI.SEQUENCE_VARIANT); + assertEquals(sf.getType(), SEQUENCE_VARIANT); 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); @@ -123,7 +141,7 @@ public class VCFLoaderTest 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); @@ -133,23 +151,33 @@ public class VCFLoaderTest 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"); + 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 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); @@ -158,7 +186,7 @@ public class VCFLoaderTest 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); @@ -178,16 +206,27 @@ public class VCFLoaderTest } } List proteinFeatures = peptide.getSequenceFeatures(); - assertEquals(proteinFeatures.size(), 1); + assertEquals(proteinFeatures.size(), 3); 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"); + + /* + * check that sequence_variant attribute AF has been clocked as + * numeric with correct min and max values + * (i.e. invalid values have been ignored - JAL-3375) + */ + FeatureAttributes fa = FeatureAttributes.getInstance(); + assertSame(fa.getDatatype(SEQUENCE_VARIANT, "AF"), Datatype.Number); + float[] minmax = fa.getMinMax(SEQUENCE_VARIANT, "AF"); + assertEquals(minmax[0], 0.002f); + assertEquals(minmax[1], 0.005f); } - private File makeVcf() throws IOException + private File makeVcfFile() throws IOException { File f = File.createTempFile("Test", ".vcf"); f.deleteOnExit(); @@ -327,7 +366,7 @@ public class VCFLoaderTest { AlignmentI al = buildAlignment(); - File f = makeVcf(); + File f = makeVcfFile(); VCFLoader loader = new VCFLoader(f.getPath()); @@ -340,96 +379,97 @@ public class VCFLoaderTest List 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 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 @@ -445,7 +485,7 @@ public class VCFLoaderTest } } List proteinFeatures = peptide.getSequenceFeatures(); - assertEquals(proteinFeatures.size(), 1); + assertEquals(proteinFeatures.size(), 3); sf = proteinFeatures.get(0); assertEquals(sf.getFeatureGroup(), "VCF"); assertEquals(sf.getBegin(), 6); -- 1.7.10.2