package jalview.io.vcf;
import static org.testng.Assert.assertEquals;
-import static org.testng.Assert.fail;
+import static org.testng.Assert.assertTrue;
import jalview.datamodel.AlignmentI;
import jalview.datamodel.DBRefEntry;
import jalview.datamodel.Sequence;
import jalview.datamodel.SequenceFeature;
import jalview.datamodel.SequenceI;
+import jalview.datamodel.features.SequenceFeatures;
import jalview.gui.AlignFrame;
import jalview.io.DataSourceType;
import jalview.io.FileLoader;
public class VCFLoaderTest
{
- // columns 9717- of gene P30419 from Ensembl (modified)
- private static final String FASTA =
- // forward strand 'gene'
- ">gene1/1-25 chromosome:GRCh38:17:45051610:45051634:1\n"
+ 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
+ */
+ ">gene1/1-25 chromosome:GRCh38:17:45051610:45051634:1\n"
+ "CAAGCTGGCGGACGAGAGTGTGACA\n"
- // and a 'made up' mini-transcript with two exons
+ ">transcript1/1-18\n--AGCTGGCG----AGAGTGTGAC-\n"
- +
- // 'reverse strand' gene (reverse complement)
- ">gene2/1-25 chromosome:GRCh38:17:45051610:45051634:-1\n"
+
+ /*
+ * reverse strand gene and transcript (reverse complement alleles!)
+ */
+ + ">gene2/1-25 chromosome:GRCh38:17:45051610:45051634:-1\n"
+ "TGTCACACTCTCGTCCGCCAGCTTG\n"
- // and its 'transcript'
- + ">transcript2/1-18\n"
- + "-GTCACACTCT----CGCCAGCT--\n";
+ + ">transcript2/1-18\n" + "-GTCACACTCT----CGCCAGCT--\n"
+
+ /*
+ * 'gene' on chromosome 5 with two transcripts
+ */
+ + ">gene3/1-25 chromosome:GRCh38:5:45051610:45051634:1\n"
+ + "CAAGCTGGCGGACGAGAGTGTGACA\n"
+ + ">transcript3/1-18\n--AGCTGGCG----AGAGTGTGAC-\n"
+ + ">transcript4/1-18\n-----TGG-GGACGAGAGTGTGA-A\n";
private static final String[] VCF = { "##fileformat=VCFv4.2",
"##INFO=<ID=AF,Number=A,Type=Float,Description=\"Allele Frequency, for each ALT allele, in the same order as listed\">",
"##reference=GRCh38",
"#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO",
- // SNP A/T in position 2 of gene sequence (precedes transcript)
- "17\t45051611\t.\tA\tT\t1666.64\tRF\tAC=15;AF=5.08130e-03",
+ // 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",
// SNP G/C in position 4 of gene sequence, position 2 of transcript
- // this is a mixed variant, the insertion G/GA is not transferred
- "17\t45051613\t.\tG\tGA,C\t1666.64\tRF\tAC=15;AF=3.08130e-03" };
+ // 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" };
@Test(groups = "Functional")
public void testDoLoad() throws IOException
/*
* verify variant feature(s) added to gene
+ * NB alleles at a locus may not be processed, and features added,
+ * in the order in which they appear in the VCF record as method
+ * VariantContext.getAlternateAlleles() does not guarantee order
+ * - order of assertions here matches what we find (is not important)
*/
List<SequenceFeature> geneFeatures = al.getSequenceAt(0)
.getSequenceFeatures();
- assertEquals(geneFeatures.size(), 2);
+ SequenceFeatures.sortFeatures(geneFeatures, true);
+ assertEquals(geneFeatures.size(), 4);
SequenceFeature sf = geneFeatures.get(0);
assertEquals(sf.getFeatureGroup(), "VCF");
assertEquals(sf.getBegin(), 2);
assertEquals(sf.getEnd(), 2);
+ assertEquals(sf.getScore(), 4.0e-03, DELTA);
+ assertEquals(sf.getValue(Gff3Helper.ALLELES), "A,C");
assertEquals(sf.getType(), SequenceOntologyI.SEQUENCE_VARIANT);
- assertEquals(sf.getScore(), 5.08130e-03, 0.000001f);
+ sf = geneFeatures.get(1);
+ assertEquals(sf.getFeatureGroup(), "VCF");
+ assertEquals(sf.getBegin(), 2);
+ assertEquals(sf.getEnd(), 2);
+ assertEquals(sf.getType(), SequenceOntologyI.SEQUENCE_VARIANT);
+ assertEquals(sf.getScore(), 5.0e-03, DELTA);
assertEquals(sf.getValue(Gff3Helper.ALLELES), "A,T");
- sf = geneFeatures.get(1);
+ sf = geneFeatures.get(2);
assertEquals(sf.getFeatureGroup(), "VCF");
assertEquals(sf.getBegin(), 4);
assertEquals(sf.getEnd(), 4);
assertEquals(sf.getType(), SequenceOntologyI.SEQUENCE_VARIANT);
- assertEquals(sf.getScore(), 3.08130e-03, 0.000001f);
+ assertEquals(sf.getScore(), 2.0e-03, DELTA);
assertEquals(sf.getValue(Gff3Helper.ALLELES), "G,C");
+ sf = geneFeatures.get(3);
+ assertEquals(sf.getFeatureGroup(), "VCF");
+ assertEquals(sf.getBegin(), 4);
+ assertEquals(sf.getEnd(), 4);
+ assertEquals(sf.getType(), SequenceOntologyI.SEQUENCE_VARIANT);
+ assertEquals(sf.getScore(), 3.0e-03, DELTA);
+ assertEquals(sf.getValue(Gff3Helper.ALLELES), "G,GA");
+
/*
* verify variant feature(s) added to transcript
*/
List<SequenceFeature> transcriptFeatures = al.getSequenceAt(1)
.getSequenceFeatures();
- assertEquals(transcriptFeatures.size(), 1);
+ assertEquals(transcriptFeatures.size(), 2);
sf = transcriptFeatures.get(0);
assertEquals(sf.getFeatureGroup(), "VCF");
assertEquals(sf.getBegin(), 2);
assertEquals(sf.getEnd(), 2);
assertEquals(sf.getType(), SequenceOntologyI.SEQUENCE_VARIANT);
- assertEquals(sf.getScore(), 3.08130e-03, 0.000001f);
+ assertEquals(sf.getScore(), 2.0e-03, DELTA);
assertEquals(sf.getValue(Gff3Helper.ALLELES), "G,C");
+ sf = transcriptFeatures.get(1);
+ assertEquals(sf.getFeatureGroup(), "VCF");
+ assertEquals(sf.getBegin(), 2);
+ assertEquals(sf.getEnd(), 2);
+ assertEquals(sf.getType(), SequenceOntologyI.SEQUENCE_VARIANT);
+ assertEquals(sf.getScore(), 3.0e-03, DELTA);
+ assertEquals(sf.getValue(Gff3Helper.ALLELES), "G,GA");
/*
- * verify variant feature(s) computed and added to protein
+ * 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();
* from Ensembl and transcripts computed)
*/
AlignmentI alignment = af.getViewport().getAlignment();
- SequenceI gene1 = alignment.getSequenceAt(0);
+ SequenceI gene1 = alignment.findName("gene1");
int[] to = new int[] { 45051610, 45051634 };
int[] from = new int[] { gene1.getStart(), gene1.getEnd() };
gene1.setGeneLoci("human", "GRCh38", "17", new MapList(from, to, 1, 1));
* which is chromosome 45051612-45051619,45051624-45051633
*/
to = new int[] { 45051612, 45051619, 45051624, 45051633 };
- SequenceI transcript1 = alignment.getSequenceAt(1);
+ SequenceI transcript1 = alignment.findName("transcript1");
from = new int[] { transcript1.getStart(), transcript1.getEnd() };
transcript1.setGeneLoci("human", "GRCh38", "17", new MapList(from, to,
1, 1));
/*
* map gene2 to chromosome reverse strand
*/
- SequenceI gene2 = alignment.getSequenceAt(2);
+ SequenceI gene2 = alignment.findName("gene2");
to = new int[] { 45051634, 45051610 };
from = new int[] { gene2.getStart(), gene2.getEnd() };
gene2.setGeneLoci("human", "GRCh38", "17", new MapList(from, to, 1, 1));
* which is chromosome 45051633-45051624,45051619-45051612
*/
to = new int[] { 45051633, 45051624, 45051619, 45051612 };
- SequenceI transcript2 = alignment.getSequenceAt(3);
+ SequenceI transcript2 = alignment.findName("transcript2");
from = new int[] { transcript2.getStart(), transcript2.getEnd() };
transcript2.setGeneLoci("human", "GRCh38", "17", new MapList(from, to,
1, 1));
product = new DBRefEntry("", "", "ENSP002", map);
transcript2.addDBRef(product);
+ /*
+ * map gene3 to chromosome
+ */
+ SequenceI gene3 = alignment.findName("gene3");
+ to = new int[] { 45051610, 45051634 };
+ from = new int[] { gene3.getStart(), gene3.getEnd() };
+ gene3.setGeneLoci("human", "GRCh38", "5", new MapList(from, to, 1, 1));
+
+ /*
+ * map 'transcript3' to chromosome
+ */
+ SequenceI transcript3 = alignment.findName("transcript3");
+ to = new int[] { 45051612, 45051619, 45051624, 45051633 };
+ from = new int[] { transcript3.getStart(), transcript3.getEnd() };
+ transcript3.setGeneLoci("human", "GRCh38", "5", new MapList(from, to,
+ 1, 1));
+
+ /*
+ * map 'transcript4' to chromosome
+ */
+ SequenceI transcript4 = alignment.findName("transcript4");
+ to = new int[] { 45051615, 45051617, 45051619, 45051632, 45051634,
+ 45051634 };
+ from = new int[] { transcript4.getStart(), transcript4.getEnd() };
+ transcript4.setGeneLoci("human", "GRCh38", "5", new MapList(from, to,
+ 1, 1));
+
+ /*
+ * add a protein product as a DBRef on transcript3
+ */
+ SequenceI peptide3 = new Sequence("ENSP003", "SWRECD");
+ mapList = new MapList(new int[] { 1, 18 }, new int[] { 1, 6 }, 3, 1);
+ map = new Mapping(peptide3, mapList);
+ product = new DBRefEntry("", "", "ENSP003", map);
+ transcript3.addDBRef(product);
+
return alignment;
}
/*
* verify variant feature(s) added to gene2
* gene/1-25 maps to chromosome 45051634- reverse strand
- * variants A/T at 45051611 and G/C at 45051613 map to
- * T/A and C/G at gene positions 24 and 22 respectively
+ * variants A/T, A/C at 45051611 and G/GA,G/C at 45051613 map to
+ * T/A, T/G and C/TC,C/G at gene positions 24 and 22 respectively
*/
List<SequenceFeature> geneFeatures = al.getSequenceAt(2)
.getSequenceFeatures();
- assertEquals(geneFeatures.size(), 2);
+ SequenceFeatures.sortFeatures(geneFeatures, true);
+ assertEquals(geneFeatures.size(), 4);
SequenceFeature sf = geneFeatures.get(0);
assertEquals(sf.getFeatureGroup(), "VCF");
assertEquals(sf.getBegin(), 22);
assertEquals(sf.getEnd(), 22);
assertEquals(sf.getType(), SequenceOntologyI.SEQUENCE_VARIANT);
- assertEquals(sf.getScore(), 3.08130e-03, 0.000001f);
+ assertEquals(sf.getScore(), 2.0e-03, DELTA);
assertEquals("C,G", sf.getValue(Gff3Helper.ALLELES));
sf = geneFeatures.get(1);
assertEquals(sf.getFeatureGroup(), "VCF");
+ assertEquals(sf.getBegin(), 22);
+ assertEquals(sf.getEnd(), 22);
+ assertEquals(sf.getType(), SequenceOntologyI.SEQUENCE_VARIANT);
+ assertEquals(sf.getScore(), 3.0e-03, DELTA);
+ assertEquals("C,TC", sf.getValue(Gff3Helper.ALLELES));
+
+ sf = geneFeatures.get(2);
+ assertEquals(sf.getFeatureGroup(), "VCF");
assertEquals(sf.getBegin(), 24);
assertEquals(sf.getEnd(), 24);
assertEquals(sf.getType(), SequenceOntologyI.SEQUENCE_VARIANT);
- assertEquals(sf.getScore(), 5.08130e-03, 0.000001f);
+ assertEquals(sf.getScore(), 4.0e-03, DELTA);
+ assertEquals("T,G", sf.getValue(Gff3Helper.ALLELES));
+
+ sf = geneFeatures.get(3);
+ assertEquals(sf.getFeatureGroup(), "VCF");
+ assertEquals(sf.getBegin(), 24);
+ assertEquals(sf.getEnd(), 24);
+ assertEquals(sf.getType(), SequenceOntologyI.SEQUENCE_VARIANT);
+ assertEquals(sf.getScore(), 5.0e-03, DELTA);
assertEquals("T,A", sf.getValue(Gff3Helper.ALLELES));
/*
* verify variant feature(s) added to transcript2
- * variant C/G at position 22 of gene overlaps and maps to
- * position 17 of transcript
+ * variants G/GA,G/C at position 22 of gene overlap and map to
+ * C/TC,C/G at position 17 of transcript
*/
List<SequenceFeature> transcriptFeatures = al.getSequenceAt(3)
.getSequenceFeatures();
- assertEquals(transcriptFeatures.size(), 1);
+ assertEquals(transcriptFeatures.size(), 2);
sf = transcriptFeatures.get(0);
assertEquals(sf.getFeatureGroup(), "VCF");
assertEquals(sf.getBegin(), 17);
assertEquals(sf.getEnd(), 17);
assertEquals(sf.getType(), SequenceOntologyI.SEQUENCE_VARIANT);
- assertEquals(sf.getScore(), 3.08130e-03, 0.000001f);
+ assertEquals(sf.getScore(), 2.0e-03, DELTA);
assertEquals("C,G", sf.getValue(Gff3Helper.ALLELES));
+ sf = transcriptFeatures.get(1);
+ assertEquals(sf.getFeatureGroup(), "VCF");
+ assertEquals(sf.getBegin(), 17);
+ assertEquals(sf.getEnd(), 17);
+ assertEquals(sf.getType(), SequenceOntologyI.SEQUENCE_VARIANT);
+ assertEquals(sf.getScore(), 3.0e-03, DELTA);
+ assertEquals("C,TC", sf.getValue(Gff3Helper.ALLELES));
+
/*
* verify variant feature(s) computed and added to protein
* last codon GCT varies to GGT giving A/G in the last peptide position
}
/**
- * Tests that where variant records have more than one SNP allele, a variant
- * feature is created for each, and the corresponding data values set on it
- *
- * @throws IOException
- */
- @Test(groups = "Functional")
- public void testDoLoad_multipleAlleles() throws IOException
- {
- fail("todo");
- }
-
- /**
* Tests that if VEP consequence (CSQ) data is present in the VCF data, then
* it is added to the variant feature, but restricted where possible to the
* consequences for a specific transcript
@Test(groups = "Functional")
public void testDoLoad_vepCsq() throws IOException
{
- fail("todo");
+ AlignmentI al = buildAlignment();
+
+ VCFLoader loader = new VCFLoader(al);
+
+ /*
+ * VCF data file with variants at gene3 positions
+ * 1 C/A
+ * 5 C/T
+ * 9 CGT/C (deletion)
+ * 13 C/G, C/T
+ * 17 A/AC (insertion), A/G
+ */
+ loader.doLoad("test/jalview/io/vcf/testVcf.dat", null);
+
+ /*
+ * verify variant feature(s) added to gene3
+ */
+ List<SequenceFeature> geneFeatures = al.findName("gene3")
+ .getSequenceFeatures();
+ SequenceFeatures.sortFeatures(geneFeatures, true);
+ assertEquals(geneFeatures.size(), 7);
+ SequenceFeature sf = geneFeatures.get(0);
+ assertEquals(sf.getBegin(), 1);
+ assertEquals(sf.getEnd(), 1);
+ assertEquals(sf.getScore(), 0.1f, DELTA);
+ assertEquals(sf.getValue("alleles"), "C,A");
+ // gene features include Consequence for all transcripts
+ assertEquals(((String) sf.getValue("CSQ")).split(",").length, 2);
+
+ sf = geneFeatures.get(1);
+ assertEquals(sf.getBegin(), 5);
+ assertEquals(sf.getEnd(), 5);
+ assertEquals(sf.getScore(), 0.2f, DELTA);
+ assertEquals(sf.getValue("alleles"), "C,T");
+ assertEquals(((String) sf.getValue("CSQ")).split(",").length, 2);
+
+ sf = geneFeatures.get(2);
+ assertEquals(sf.getBegin(), 9);
+ assertEquals(sf.getEnd(), 11); // deletion over 3 positions
+ assertEquals(sf.getScore(), 0.3f, DELTA);
+ assertEquals(sf.getValue("alleles"), "CGG,C");
+ assertEquals(((String) sf.getValue("CSQ")).split(",").length, 2);
+
+ sf = geneFeatures.get(3);
+ assertEquals(sf.getBegin(), 13);
+ assertEquals(sf.getEnd(), 13);
+ assertEquals(sf.getScore(), 0.5f, DELTA);
+ assertEquals(sf.getValue("alleles"), "C,T");
+ assertEquals(((String) sf.getValue("CSQ")).split(",").length, 2);
+
+ sf = geneFeatures.get(4);
+ assertEquals(sf.getBegin(), 13);
+ assertEquals(sf.getEnd(), 13);
+ assertEquals(sf.getScore(), 0.4f, DELTA);
+ assertEquals(sf.getValue("alleles"), "C,G");
+ assertEquals(((String) sf.getValue("CSQ")).split(",").length, 2);
+
+ sf = geneFeatures.get(5);
+ assertEquals(sf.getBegin(), 17);
+ assertEquals(sf.getEnd(), 17);
+ assertEquals(sf.getScore(), 0.7f, DELTA);
+ assertEquals(sf.getValue("alleles"), "A,G");
+ assertEquals(((String) sf.getValue("CSQ")).split(",").length, 2);
+
+ sf = geneFeatures.get(6);
+ assertEquals(sf.getBegin(), 17);
+ assertEquals(sf.getEnd(), 17); // insertion
+ assertEquals(sf.getScore(), 0.6f, DELTA);
+ assertEquals(sf.getValue("alleles"), "A,AC");
+ assertEquals(((String) sf.getValue("CSQ")).split(",").length, 2);
+
+ /*
+ * verify variant feature(s) added to transcript3
+ * at columns 5 (1), 17 (2), positions 3, 11
+ * note the deletion at columns 9-11 is not transferred since col 11
+ * has no mapping to transcript 3
+ */
+ List<SequenceFeature> transcriptFeatures = al.findName("transcript3")
+ .getSequenceFeatures();
+ SequenceFeatures.sortFeatures(transcriptFeatures, true);
+ assertEquals(transcriptFeatures.size(), 3);
+ sf = transcriptFeatures.get(0);
+ assertEquals(sf.getBegin(), 3);
+ assertEquals(sf.getEnd(), 3);
+ assertEquals(sf.getScore(), 0.2f, DELTA);
+ assertEquals(sf.getValue("alleles"), "C,T");
+ // transcript features only have Consequence for that transcripts
+ assertEquals(((String) sf.getValue("CSQ")).split(",").length, 1);
+ assertTrue(sf.getValue("CSQ").toString().contains("transcript3"));
+
+ sf = transcriptFeatures.get(1);
+ assertEquals(sf.getBegin(), 11);
+ assertEquals(sf.getEnd(), 11);
+ assertEquals(sf.getScore(), 0.7f, DELTA);
+ assertEquals(sf.getValue("alleles"), "A,G");
+ assertEquals(((String) sf.getValue("CSQ")).split(",").length, 1);
+ assertTrue(sf.getValue("CSQ").toString().contains("transcript3"));
+
+ sf = transcriptFeatures.get(2);
+ assertEquals(sf.getBegin(), 11);
+ assertEquals(sf.getEnd(), 11);
+ assertEquals(sf.getScore(), 0.6f, DELTA);
+ assertEquals(sf.getValue("alleles"), "A,AC");
+ assertEquals(((String) sf.getValue("CSQ")).split(",").length, 1);
+ assertTrue(sf.getValue("CSQ").toString().contains("transcript3"));
+
+ /*
+ * verify variant feature(s) added to transcript4
+ * at columns 13 (2) and 17 (2), positions 7 and 11
+ */
+ transcriptFeatures = al.findName("transcript4").getSequenceFeatures();
+ SequenceFeatures.sortFeatures(transcriptFeatures, true);
+ assertEquals(transcriptFeatures.size(), 4);
+ sf = transcriptFeatures.get(0);
+ assertEquals(sf.getBegin(), 7);
+ assertEquals(sf.getEnd(), 7);
+ assertEquals(sf.getScore(), 0.5f, DELTA);
+ assertEquals(sf.getValue("alleles"), "C,T");
+ assertEquals(((String) sf.getValue("CSQ")).split(",").length, 1);
+ assertTrue(sf.getValue("CSQ").toString().contains("transcript4"));
+
+ sf = transcriptFeatures.get(1);
+ assertEquals(sf.getBegin(), 7);
+ assertEquals(sf.getEnd(), 7);
+ assertEquals(sf.getScore(), 0.4f, DELTA);
+ assertEquals(sf.getValue("alleles"), "C,G");
+ assertEquals(((String) sf.getValue("CSQ")).split(",").length, 1);
+ assertTrue(sf.getValue("CSQ").toString().contains("transcript4"));
+
+ sf = transcriptFeatures.get(2);
+ assertEquals(sf.getBegin(), 11);
+ assertEquals(sf.getEnd(), 11);
+ assertEquals(sf.getScore(), 0.7f, DELTA);
+ assertEquals(sf.getValue("alleles"), "A,G");
+ assertEquals(((String) sf.getValue("CSQ")).split(",").length, 1);
+ assertTrue(sf.getValue("CSQ").toString().contains("transcript4"));
+
+ sf = transcriptFeatures.get(3);
+ assertEquals(sf.getBegin(), 11);
+ assertEquals(sf.getEnd(), 11);
+ assertEquals(sf.getScore(), 0.6f, DELTA);
+ assertEquals(sf.getValue("alleles"), "A,AC");
+ assertEquals(((String) sf.getValue("CSQ")).split(",").length, 1);
+ assertTrue(sf.getValue("CSQ").toString().contains("transcript4"));
}
}
git://source.jalview.org / jalview.git / commitdiff