X-Git-Url: http://source.jalview.org/gitweb/?a=blobdiff_plain;f=test%2Fjalview%2Fio%2Fvcf%2FVCFLoaderTest.java;h=f512d0b297c8fea4d543974433a98e3a88c0f0fc;hb=57738a1f3c19b1c3a00bd3ac5108f8cd0af32f99;hp=a87c160c458d3466337f489c836cf191c29c1119;hpb=1b0014839c94c711c09a94641ecb5036072bb5b3;p=jalview.git diff --git a/test/jalview/io/vcf/VCFLoaderTest.java b/test/jalview/io/vcf/VCFLoaderTest.java index a87c160..f512d0b 100644 --- a/test/jalview/io/vcf/VCFLoaderTest.java +++ b/test/jalview/io/vcf/VCFLoaderTest.java @@ -1,20 +1,24 @@ 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.assertTrue; import jalview.bin.Cache; +import jalview.bin.Console; import jalview.datamodel.AlignmentI; import jalview.datamodel.DBRefEntry; 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.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; @@ -24,6 +28,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 @@ -31,11 +36,10 @@ public class VCFLoaderTest private static final float DELTA = 0.00001f; // columns 9717- of gene P30419 from Ensembl (much modified) - private static final String FASTA = "" - + - /* - * forward strand 'gene' and 'transcript' with two exons - */ + private static final String FASTA = "" + + /* + * forward strand 'gene' and 'transcript' with two exons + */ ">gene1/1-25 chromosome:GRCh38:17:45051610:45051634:1\n" + "CAAGCTGGCGGACGAGAGTGTGACA\n" + ">transcript1/1-18\n--AGCTGGCG----AGAGTGTGAC-\n" @@ -44,8 +48,8 @@ public class VCFLoaderTest * reverse strand gene and transcript (reverse complement alleles!) */ + ">gene2/1-25 chromosome:GRCh38:17:45051610:45051634:-1\n" - + "TGTCACACTCTCGTCCGCCAGCTTG\n" - + ">transcript2/1-18\n" + "-GTCACACTCT----CGCCAGCT--\n" + + "TGTCACACTCTCGTCCGCCAGCTTG\n" + ">transcript2/1-18\n" + + "-GTCACACTCT----CGCCAGCT--\n" /* * 'gene' on chromosome 5 with two transcripts @@ -56,15 +60,22 @@ 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=", + "##INFO= 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); @@ -133,23 +163,35 @@ 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"); + 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 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 +200,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); @@ -168,7 +210,7 @@ public class VCFLoaderTest * verify SNP variant feature(s) computed and added to protein * first codon AGC varies to ACC giving S/T */ - DBRefEntry[] dbRefs = al.getSequenceAt(1).getDBRefs(); + List dbRefs = al.getSequenceAt(1).getDBRefs(); SequenceI peptide = null; for (DBRefEntry dbref : dbRefs) { @@ -178,16 +220,14 @@ public class VCFLoaderTest } } List proteinFeatures = peptide.getSequenceFeatures(); - 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"); + + /* + * JAL-3187 don't precompute protein features, do dynamically instead + */ + assertTrue(proteinFeatures.isEmpty()); } - private File makeVcf() throws IOException + private File makeVcfFile() throws IOException { File f = File.createTempFile("Test", ".vcf"); f.deleteOnExit(); @@ -218,8 +258,8 @@ public class VCFLoaderTest SequenceI gene1 = alignment.findName("gene1"); int[] to = new int[] { 45051610, 45051634 }; int[] from = new int[] { gene1.getStart(), gene1.getEnd() }; - gene1.setGeneLoci("homo_sapiens", "GRCh38", "17", new MapList(from, to, - 1, 1)); + gene1.setGeneLoci("homo_sapiens", "GRCh38", "17", + new MapList(from, to, 1, 1)); /* * map 'transcript1' to chromosome via 'gene1' @@ -229,9 +269,8 @@ public class VCFLoaderTest to = new int[] { 45051612, 45051619, 45051624, 45051633 }; SequenceI transcript1 = alignment.findName("transcript1"); from = new int[] { transcript1.getStart(), transcript1.getEnd() }; - transcript1.setGeneLoci("homo_sapiens", "GRCh38", "17", new MapList( - from, to, - 1, 1)); + transcript1.setGeneLoci("homo_sapiens", "GRCh38", "17", + new MapList(from, to, 1, 1)); /* * map gene2 to chromosome reverse strand @@ -239,8 +278,8 @@ public class VCFLoaderTest SequenceI gene2 = alignment.findName("gene2"); to = new int[] { 45051634, 45051610 }; from = new int[] { gene2.getStart(), gene2.getEnd() }; - gene2.setGeneLoci("homo_sapiens", "GRCh38", "17", new MapList(from, to, - 1, 1)); + gene2.setGeneLoci("homo_sapiens", "GRCh38", "17", + new MapList(from, to, 1, 1)); /* * map 'transcript2' to chromosome via 'gene2' @@ -250,9 +289,8 @@ public class VCFLoaderTest to = new int[] { 45051633, 45051624, 45051619, 45051612 }; SequenceI transcript2 = alignment.findName("transcript2"); from = new int[] { transcript2.getStart(), transcript2.getEnd() }; - transcript2.setGeneLoci("homo_sapiens", "GRCh38", "17", new MapList( - from, to, - 1, 1)); + transcript2.setGeneLoci("homo_sapiens", "GRCh38", "17", + new MapList(from, to, 1, 1)); /* * add a protein product as a DBRef on transcript1 @@ -279,8 +317,8 @@ public class VCFLoaderTest SequenceI gene3 = alignment.findName("gene3"); to = new int[] { 45051610, 45051634 }; from = new int[] { gene3.getStart(), gene3.getEnd() }; - gene3.setGeneLoci("homo_sapiens", "GRCh38", "5", new MapList(from, to, - 1, 1)); + gene3.setGeneLoci("homo_sapiens", "GRCh38", "5", + new MapList(from, to, 1, 1)); /* * map 'transcript3' to chromosome @@ -288,9 +326,8 @@ public class VCFLoaderTest SequenceI transcript3 = alignment.findName("transcript3"); to = new int[] { 45051612, 45051619, 45051624, 45051633 }; from = new int[] { transcript3.getStart(), transcript3.getEnd() }; - transcript3.setGeneLoci("homo_sapiens", "GRCh38", "5", new MapList( - from, to, - 1, 1)); + transcript3.setGeneLoci("homo_sapiens", "GRCh38", "5", + new MapList(from, to, 1, 1)); /* * map 'transcript4' to chromosome @@ -299,9 +336,8 @@ public class VCFLoaderTest to = new int[] { 45051615, 45051617, 45051619, 45051632, 45051634, 45051634 }; from = new int[] { transcript4.getStart(), transcript4.getEnd() }; - transcript4.setGeneLoci("homo_sapiens", "GRCh38", "5", new MapList( - from, to, - 1, 1)); + transcript4.setGeneLoci("homo_sapiens", "GRCh38", "5", + new MapList(from, to, 1, 1)); /* * add a protein product as a DBRef on transcript3 @@ -327,7 +363,7 @@ public class VCFLoaderTest { AlignmentI al = buildAlignment(); - File f = makeVcf(); + File f = makeVcfFile(); VCFLoader loader = new VCFLoader(f.getPath()); @@ -340,102 +376,103 @@ 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 * last codon GCT varies to GGT giving A/G in the last peptide position */ - DBRefEntry[] dbRefs = al.getSequenceAt(3).getDBRefs(); + List dbRefs = al.getSequenceAt(3).getDBRefs(); SequenceI peptide = null; for (DBRefEntry dbref : dbRefs) { @@ -445,13 +482,11 @@ public class VCFLoaderTest } } List proteinFeatures = peptide.getSequenceFeatures(); - 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"); + + /* + * JAL-3187 don't precompute protein features, do dynamically instead + */ + assertTrue(proteinFeatures.isEmpty()); } /** @@ -494,6 +529,7 @@ public class VCFLoaderTest // gene features include Consequence for all transcripts Map map = (Map) sf.getValue("CSQ"); assertEquals(map.size(), 9); + assertEquals(map.get("PolyPhen"), "Bad"); sf = geneFeatures.get(1); assertEquals(sf.getBegin(), 5); @@ -503,6 +539,7 @@ public class VCFLoaderTest assertEquals(sf.getValue("alleles"), "C,T"); map = (Map) sf.getValue("CSQ"); assertEquals(map.size(), 9); + assertEquals(map.get("PolyPhen"), "Bad;;"); // %3B%3B decoded sf = geneFeatures.get(2); assertEquals(sf.getBegin(), 9); @@ -595,7 +632,7 @@ public class VCFLoaderTest * and GAG/GGG which is E/G in position 4 * the insertion variant is not transferred to the peptide */ - DBRefEntry[] dbRefs = al.findName("transcript3").getDBRefs(); + List dbRefs = al.findName("transcript3").getDBRefs(); SequenceI peptide = null; for (DBRefEntry dbref : dbRefs) { @@ -605,20 +642,24 @@ public class VCFLoaderTest } } List proteinFeatures = peptide.getSequenceFeatures(); - SequenceFeatures.sortFeatures(proteinFeatures, true); - assertEquals(proteinFeatures.size(), 2); - sf = proteinFeatures.get(0); - assertEquals(sf.getFeatureGroup(), "VCF"); - assertEquals(sf.getBegin(), 1); - assertEquals(sf.getEnd(), 1); - assertEquals(sf.getType(), SequenceOntologyI.SYNONYMOUS_VARIANT); - assertEquals(sf.getDescription(), "agC/agT"); - sf = proteinFeatures.get(1); - assertEquals(sf.getFeatureGroup(), "VCF"); - assertEquals(sf.getBegin(), 4); - assertEquals(sf.getEnd(), 4); - assertEquals(sf.getType(), SequenceOntologyI.NONSYNONYMOUS_VARIANT); - assertEquals(sf.getDescription(), "p.Glu4Gly"); + /* + * JAL-3187 don't precompute protein features, do dynamically instead + */ + assertTrue(proteinFeatures.isEmpty()); + // SequenceFeatures.sortFeatures(proteinFeatures, true); + // assertEquals(proteinFeatures.size(), 2); + // sf = proteinFeatures.get(0); + // assertEquals(sf.getFeatureGroup(), "VCF"); + // assertEquals(sf.getBegin(), 1); + // assertEquals(sf.getEnd(), 1); + // assertEquals(sf.getType(), SequenceOntologyI.SYNONYMOUS_VARIANT); + // assertEquals(sf.getDescription(), "agC/agT"); + // sf = proteinFeatures.get(1); + // assertEquals(sf.getFeatureGroup(), "VCF"); + // assertEquals(sf.getBegin(), 4); + // assertEquals(sf.getEnd(), 4); + // assertEquals(sf.getType(), SequenceOntologyI.NONSYNONYMOUS_VARIANT); + // assertEquals(sf.getDescription(), "p.Glu4Gly"); /* * verify variant feature(s) added to transcript4 @@ -673,8 +714,7 @@ public class VCFLoaderTest @Test(groups = "Functional") public void testLoadVCFContig() throws IOException { - VCFLoader loader = new VCFLoader( - "test/jalview/io/vcf/testVcf2.vcf"); + VCFLoader loader = new VCFLoader("test/jalview/io/vcf/testVcf2.vcf"); SequenceI seq = loader.loadVCFContig("contig123"); assertEquals(seq.getLength(), 15);