a02cc5a82c587253b55355bf28b1472d695b89a8
[jalview.git] / test / jalview / io / vcf / VCFLoaderTest.java
1 package jalview.io.vcf;
2
3 import static org.testng.Assert.assertEquals;
4
5 import jalview.bin.Cache;
6 import jalview.datamodel.AlignmentI;
7 import jalview.datamodel.DBRefEntry;
8 import jalview.datamodel.Mapping;
9 import jalview.datamodel.Sequence;
10 import jalview.datamodel.SequenceFeature;
11 import jalview.datamodel.SequenceI;
12 import jalview.datamodel.features.SequenceFeatures;
13 import jalview.gui.AlignFrame;
14 import jalview.io.DataSourceType;
15 import jalview.io.FileLoader;
16 import jalview.io.gff.Gff3Helper;
17 import jalview.io.gff.SequenceOntologyI;
18 import jalview.util.MapList;
19
20 import java.io.File;
21 import java.io.IOException;
22 import java.io.PrintWriter;
23 import java.util.List;
24 import java.util.Map;
25
26 import org.testng.annotations.BeforeClass;
27 import org.testng.annotations.Test;
28
29 public class VCFLoaderTest
30 {
31   private static final float DELTA = 0.00001f;
32
33   // columns 9717- of gene P30419 from Ensembl (much modified)
34   private static final String FASTA = ""
35           +
36           /*
37            * forward strand 'gene' and 'transcript' with two exons
38            */
39           ">gene1/1-25 chromosome:GRCh38:17:45051610:45051634:1\n"
40           + "CAAGCTGGCGGACGAGAGTGTGACA\n"
41           + ">transcript1/1-18\n--AGCTGGCG----AGAGTGTGAC-\n"
42
43           /*
44            * reverse strand gene and transcript (reverse complement alleles!)
45            */
46           + ">gene2/1-25 chromosome:GRCh38:17:45051610:45051634:-1\n"
47           + "TGTCACACTCTCGTCCGCCAGCTTG\n"
48           + ">transcript2/1-18\n" + "-GTCACACTCT----CGCCAGCT--\n"
49
50           /*
51            * 'gene' on chromosome 5 with two transcripts
52            */
53           + ">gene3/1-25 chromosome:GRCh38:5:45051610:45051634:1\n"
54           + "CAAGCTGGCGGACGAGAGTGTGACA\n"
55           + ">transcript3/1-18\n--AGCTGGCG----AGAGTGTGAC-\n"
56           + ">transcript4/1-18\n-----TGG-GGACGAGAGTGTGA-A\n";
57
58   private static final String[] VCF = { "##fileformat=VCFv4.2",
59       "##INFO=<ID=AF,Number=A,Type=Float,Description=\"Allele Frequency, for each ALT allele, in the same order as listed\">",
60       "##reference=Homo_sapiens/GRCh38",
61       "#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO",
62       // A/T,C variants in position 2 of gene sequence (precedes transcript)
63       // should create 2 variant features with respective scores
64       "17\t45051611\t.\tA\tT,C\t1666.64\tRF\tAC=15;AF=5.0e-03,4.0e-03",
65       // SNP G/C in position 4 of gene sequence, position 2 of transcript
66       // insertion G/GA is transferred to nucleotide but not to peptide
67       "17\t45051613\t.\tG\tGA,C\t1666.64\tRF\tAC=15;AF=3.0e-03,2.0e-03" };
68
69   @BeforeClass
70   public void setUp()
71   {
72     /*
73      * configure to capture all available VCF and VEP (CSQ) fields
74      */
75     Cache.loadProperties("test/jalview/io/testProps.jvprops");
76     Cache.setProperty("VCF_FIELDS", ".*");
77     Cache.setProperty("VEP_FIELDS", ".*");
78     Cache.initLogger();
79   }
80
81   @Test(groups = "Functional")
82   public void testDoLoad() throws IOException
83   {
84     AlignmentI al = buildAlignment();
85
86     File f = makeVcf();
87     VCFLoader loader = new VCFLoader(f.getPath());
88
89     loader.doLoad(al.getSequencesArray(), null);
90
91     /*
92      * verify variant feature(s) added to gene
93      * NB alleles at a locus may not be processed, and features added,
94      * in the order in which they appear in the VCF record as method
95      * VariantContext.getAlternateAlleles() does not guarantee order
96      * - order of assertions here matches what we find (is not important) 
97      */
98     List<SequenceFeature> geneFeatures = al.getSequenceAt(0)
99             .getSequenceFeatures();
100     SequenceFeatures.sortFeatures(geneFeatures, true);
101     assertEquals(geneFeatures.size(), 4);
102     SequenceFeature sf = geneFeatures.get(0);
103     assertEquals(sf.getFeatureGroup(), "VCF");
104     assertEquals(sf.getBegin(), 2);
105     assertEquals(sf.getEnd(), 2);
106     assertEquals(sf.getScore(), 4.0e-03, DELTA);
107     assertEquals(sf.getValue(Gff3Helper.ALLELES), "A,C");
108     assertEquals(sf.getType(), SequenceOntologyI.SEQUENCE_VARIANT);
109     sf = geneFeatures.get(1);
110     assertEquals(sf.getFeatureGroup(), "VCF");
111     assertEquals(sf.getBegin(), 2);
112     assertEquals(sf.getEnd(), 2);
113     assertEquals(sf.getType(), SequenceOntologyI.SEQUENCE_VARIANT);
114     assertEquals(sf.getScore(), 5.0e-03, DELTA);
115     assertEquals(sf.getValue(Gff3Helper.ALLELES), "A,T");
116
117     sf = geneFeatures.get(2);
118     assertEquals(sf.getFeatureGroup(), "VCF");
119     assertEquals(sf.getBegin(), 4);
120     assertEquals(sf.getEnd(), 4);
121     assertEquals(sf.getType(), SequenceOntologyI.SEQUENCE_VARIANT);
122     assertEquals(sf.getScore(), 2.0e-03, DELTA);
123     assertEquals(sf.getValue(Gff3Helper.ALLELES), "G,C");
124
125     sf = geneFeatures.get(3);
126     assertEquals(sf.getFeatureGroup(), "VCF");
127     assertEquals(sf.getBegin(), 4);
128     assertEquals(sf.getEnd(), 4);
129     assertEquals(sf.getType(), SequenceOntologyI.SEQUENCE_VARIANT);
130     assertEquals(sf.getScore(), 3.0e-03, DELTA);
131     assertEquals(sf.getValue(Gff3Helper.ALLELES), "G,GA");
132
133     /*
134      * verify variant feature(s) added to transcript
135      */
136     List<SequenceFeature> transcriptFeatures = al.getSequenceAt(1)
137             .getSequenceFeatures();
138     assertEquals(transcriptFeatures.size(), 2);
139     sf = transcriptFeatures.get(0);
140     assertEquals(sf.getFeatureGroup(), "VCF");
141     assertEquals(sf.getBegin(), 2);
142     assertEquals(sf.getEnd(), 2);
143     assertEquals(sf.getType(), SequenceOntologyI.SEQUENCE_VARIANT);
144     assertEquals(sf.getScore(), 2.0e-03, DELTA);
145     assertEquals(sf.getValue(Gff3Helper.ALLELES), "G,C");
146     sf = transcriptFeatures.get(1);
147     assertEquals(sf.getFeatureGroup(), "VCF");
148     assertEquals(sf.getBegin(), 2);
149     assertEquals(sf.getEnd(), 2);
150     assertEquals(sf.getType(), SequenceOntologyI.SEQUENCE_VARIANT);
151     assertEquals(sf.getScore(), 3.0e-03, DELTA);
152     assertEquals(sf.getValue(Gff3Helper.ALLELES), "G,GA");
153
154     /*
155      * verify SNP variant feature(s) computed and added to protein
156      * first codon AGC varies to ACC giving S/T
157      */
158     DBRefEntry[] dbRefs = al.getSequenceAt(1).getDBRefs();
159     SequenceI peptide = null;
160     for (DBRefEntry dbref : dbRefs)
161     {
162       if (dbref.getMap().getMap().getFromRatio() == 3)
163       {
164         peptide = dbref.getMap().getTo();
165       }
166     }
167     List<SequenceFeature> proteinFeatures = peptide.getSequenceFeatures();
168     assertEquals(proteinFeatures.size(), 1);
169     sf = proteinFeatures.get(0);
170     assertEquals(sf.getFeatureGroup(), "VCF");
171     assertEquals(sf.getBegin(), 1);
172     assertEquals(sf.getEnd(), 1);
173     assertEquals(sf.getType(), SequenceOntologyI.SEQUENCE_VARIANT);
174     assertEquals(sf.getDescription(), "p.Ser1Thr");
175   }
176
177   private File makeVcf() throws IOException
178   {
179     File f = File.createTempFile("Test", ".vcf");
180     f.deleteOnExit();
181     PrintWriter pw = new PrintWriter(f);
182     for (String vcfLine : VCF)
183     {
184       pw.println(vcfLine);
185     }
186     pw.close();
187     return f;
188   }
189
190   /**
191    * Make a simple alignment with one 'gene' and one 'transcript'
192    * 
193    * @return
194    */
195   private AlignmentI buildAlignment()
196   {
197     AlignFrame af = new FileLoader().LoadFileWaitTillLoaded(FASTA,
198             DataSourceType.PASTE);
199
200     /*
201      * map gene1 sequence to chromosome (normally done when the sequence is fetched
202      * from Ensembl and transcripts computed)
203      */
204     AlignmentI alignment = af.getViewport().getAlignment();
205     SequenceI gene1 = alignment.findName("gene1");
206     int[] to = new int[] { 45051610, 45051634 };
207     int[] from = new int[] { gene1.getStart(), gene1.getEnd() };
208     gene1.setGeneLoci("homo_sapiens", "GRCh38", "17", new MapList(from, to,
209             1, 1));
210
211     /*
212      * map 'transcript1' to chromosome via 'gene1'
213      * transcript1/1-18 is gene1/3-10,15-24
214      * which is chromosome 45051612-45051619,45051624-45051633
215      */
216     to = new int[] { 45051612, 45051619, 45051624, 45051633 };
217     SequenceI transcript1 = alignment.findName("transcript1");
218     from = new int[] { transcript1.getStart(), transcript1.getEnd() };
219     transcript1.setGeneLoci("homo_sapiens", "GRCh38", "17", new MapList(
220             from, to,
221             1, 1));
222
223     /*
224      * map gene2 to chromosome reverse strand
225      */
226     SequenceI gene2 = alignment.findName("gene2");
227     to = new int[] { 45051634, 45051610 };
228     from = new int[] { gene2.getStart(), gene2.getEnd() };
229     gene2.setGeneLoci("homo_sapiens", "GRCh38", "17", new MapList(from, to,
230             1, 1));
231
232     /*
233      * map 'transcript2' to chromosome via 'gene2'
234      * transcript2/1-18 is gene2/2-11,16-23
235      * which is chromosome 45051633-45051624,45051619-45051612
236      */
237     to = new int[] { 45051633, 45051624, 45051619, 45051612 };
238     SequenceI transcript2 = alignment.findName("transcript2");
239     from = new int[] { transcript2.getStart(), transcript2.getEnd() };
240     transcript2.setGeneLoci("homo_sapiens", "GRCh38", "17", new MapList(
241             from, to,
242             1, 1));
243
244     /*
245      * add a protein product as a DBRef on transcript1
246      */
247     SequenceI peptide1 = new Sequence("ENSP001", "SWRECD");
248     MapList mapList = new MapList(new int[] { 1, 18 }, new int[] { 1, 6 },
249             3, 1);
250     Mapping map = new Mapping(peptide1, mapList);
251     DBRefEntry product = new DBRefEntry("", "", "ENSP001", map);
252     transcript1.addDBRef(product);
253
254     /*
255      * add a protein product as a DBRef on transcript2
256      */
257     SequenceI peptide2 = new Sequence("ENSP002", "VTLSPA");
258     mapList = new MapList(new int[] { 1, 18 }, new int[] { 1, 6 }, 3, 1);
259     map = new Mapping(peptide2, mapList);
260     product = new DBRefEntry("", "", "ENSP002", map);
261     transcript2.addDBRef(product);
262
263     /*
264      * map gene3 to chromosome 
265      */
266     SequenceI gene3 = alignment.findName("gene3");
267     to = new int[] { 45051610, 45051634 };
268     from = new int[] { gene3.getStart(), gene3.getEnd() };
269     gene3.setGeneLoci("homo_sapiens", "GRCh38", "5", new MapList(from, to,
270             1, 1));
271
272     /*
273      * map 'transcript3' to chromosome
274      */
275     SequenceI transcript3 = alignment.findName("transcript3");
276     to = new int[] { 45051612, 45051619, 45051624, 45051633 };
277     from = new int[] { transcript3.getStart(), transcript3.getEnd() };
278     transcript3.setGeneLoci("homo_sapiens", "GRCh38", "5", new MapList(
279             from, to,
280             1, 1));
281
282     /*
283      * map 'transcript4' to chromosome
284      */
285     SequenceI transcript4 = alignment.findName("transcript4");
286     to = new int[] { 45051615, 45051617, 45051619, 45051632, 45051634,
287         45051634 };
288     from = new int[] { transcript4.getStart(), transcript4.getEnd() };
289     transcript4.setGeneLoci("homo_sapiens", "GRCh38", "5", new MapList(
290             from, to,
291             1, 1));
292
293     /*
294      * add a protein product as a DBRef on transcript3
295      */
296     SequenceI peptide3 = new Sequence("ENSP003", "SWRECD");
297     mapList = new MapList(new int[] { 1, 18 }, new int[] { 1, 6 }, 3, 1);
298     map = new Mapping(peptide3, mapList);
299     product = new DBRefEntry("", "", "ENSP003", map);
300     transcript3.addDBRef(product);
301
302     return alignment;
303   }
304
305   /**
306    * Test with 'gene' and 'transcript' mapped to the reverse strand of the
307    * chromosome. The VCF variant positions (in forward coordinates) should get
308    * correctly located on sequence positions.
309    * 
310    * @throws IOException
311    */
312   @Test(groups = "Functional")
313   public void testDoLoad_reverseStrand() throws IOException
314   {
315     AlignmentI al = buildAlignment();
316
317     File f = makeVcf();
318
319     VCFLoader loader = new VCFLoader(f.getPath());
320
321     loader.doLoad(al.getSequencesArray(), null);
322
323     /*
324      * verify variant feature(s) added to gene2
325      * gene2/1-25 maps to chromosome 45051634- reverse strand
326      */
327     List<SequenceFeature> geneFeatures = al.getSequenceAt(2)
328             .getSequenceFeatures();
329     SequenceFeatures.sortFeatures(geneFeatures, true);
330     assertEquals(geneFeatures.size(), 4);
331
332     /*
333      * variant A/T at 45051611 maps to T/A at gene position 24
334      */
335     SequenceFeature sf = geneFeatures.get(3);
336     assertEquals(sf.getFeatureGroup(), "VCF");
337     assertEquals(sf.getBegin(), 24);
338     assertEquals(sf.getEnd(), 24);
339     assertEquals(sf.getType(), SequenceOntologyI.SEQUENCE_VARIANT);
340     assertEquals(sf.getScore(), 5.0e-03, DELTA);
341     assertEquals(sf.getValue(Gff3Helper.ALLELES), "T,A");
342
343     /*
344      * variant A/C at 45051611 maps to T/G at gene position 24
345      */
346     sf = geneFeatures.get(2);
347     assertEquals(sf.getFeatureGroup(), "VCF");
348     assertEquals(sf.getBegin(), 24);
349     assertEquals(sf.getEnd(), 24);
350     assertEquals(sf.getType(), SequenceOntologyI.SEQUENCE_VARIANT);
351     assertEquals(sf.getScore(), 4.0e-03, DELTA);
352     assertEquals(sf.getValue(Gff3Helper.ALLELES), "T,G");
353
354     /*
355      * variant G/C at 45051613 maps to C/G at gene position 22
356      */
357     sf = geneFeatures.get(1);
358     assertEquals(sf.getFeatureGroup(), "VCF");
359     assertEquals(sf.getBegin(), 22);
360     assertEquals(sf.getEnd(), 22);
361     assertEquals(sf.getType(), SequenceOntologyI.SEQUENCE_VARIANT);
362     assertEquals(sf.getScore(), 2.0e-03, DELTA);
363     assertEquals(sf.getValue(Gff3Helper.ALLELES), "C,G");
364
365     /*
366      * insertion G/GA at 45051613 maps to an insertion at
367      * the preceding position (21) on reverse strand gene
368      * reference: CAAGC -> GCTTG/21-25
369      * genomic variant: CAAGAC (G/GA)
370      * gene variant: GTCTTG (G/GT at 21)
371      */
372     sf = geneFeatures.get(0);
373     assertEquals(sf.getFeatureGroup(), "VCF");
374     assertEquals(sf.getBegin(), 21);
375     assertEquals(sf.getEnd(), 21);
376     assertEquals(sf.getType(), SequenceOntologyI.SEQUENCE_VARIANT);
377     assertEquals(sf.getScore(), 3.0e-03, DELTA);
378     assertEquals(sf.getValue(Gff3Helper.ALLELES), "G,GT");
379
380     /*
381      * verify 2 variant features added to transcript2
382      */
383     List<SequenceFeature> transcriptFeatures = al.getSequenceAt(3)
384             .getSequenceFeatures();
385     assertEquals(transcriptFeatures.size(), 2);
386
387     /*
388      * insertion G/GT at position 21 of gene maps to position 16 of transcript
389      */
390     sf = transcriptFeatures.get(0);
391     assertEquals(sf.getFeatureGroup(), "VCF");
392     assertEquals(sf.getBegin(), 16);
393     assertEquals(sf.getEnd(), 16);
394     assertEquals(sf.getType(), SequenceOntologyI.SEQUENCE_VARIANT);
395     assertEquals(sf.getScore(), 3.0e-03, DELTA);
396     assertEquals(sf.getValue(Gff3Helper.ALLELES), "G,GT");
397
398     /*
399      * SNP C/G at position 22 of gene maps to position 17 of transcript
400      */
401     sf = transcriptFeatures.get(1);
402     assertEquals(sf.getFeatureGroup(), "VCF");
403     assertEquals(sf.getBegin(), 17);
404     assertEquals(sf.getEnd(), 17);
405     assertEquals(sf.getType(), SequenceOntologyI.SEQUENCE_VARIANT);
406     assertEquals(sf.getScore(), 2.0e-03, DELTA);
407     assertEquals(sf.getValue(Gff3Helper.ALLELES), "C,G");
408
409     /*
410      * verify variant feature(s) computed and added to protein
411      * last codon GCT varies to GGT giving A/G in the last peptide position
412      */
413     DBRefEntry[] dbRefs = al.getSequenceAt(3).getDBRefs();
414     SequenceI peptide = null;
415     for (DBRefEntry dbref : dbRefs)
416     {
417       if (dbref.getMap().getMap().getFromRatio() == 3)
418       {
419         peptide = dbref.getMap().getTo();
420       }
421     }
422     List<SequenceFeature> proteinFeatures = peptide.getSequenceFeatures();
423     assertEquals(proteinFeatures.size(), 1);
424     sf = proteinFeatures.get(0);
425     assertEquals(sf.getFeatureGroup(), "VCF");
426     assertEquals(sf.getBegin(), 6);
427     assertEquals(sf.getEnd(), 6);
428     assertEquals(sf.getType(), SequenceOntologyI.SEQUENCE_VARIANT);
429     assertEquals(sf.getDescription(), "p.Ala6Gly");
430   }
431
432   /**
433    * Tests that if VEP consequence (CSQ) data is present in the VCF data, then
434    * it is added to the variant feature, but restricted where possible to the
435    * consequences for a specific transcript
436    * 
437    * @throws IOException
438    */
439   @Test(groups = "Functional")
440   public void testDoLoad_vepCsq() throws IOException
441   {
442     AlignmentI al = buildAlignment();
443
444     VCFLoader loader = new VCFLoader("test/jalview/io/vcf/testVcf.vcf");
445
446     /*
447      * VCF data file with variants at gene3 positions
448      * 1 C/A
449      * 5 C/T
450      * 9 CGT/C (deletion)
451      * 13 C/G, C/T
452      * 17 A/AC (insertion), A/G
453      */
454     loader.doLoad(al.getSequencesArray(), null);
455
456     /*
457      * verify variant feature(s) added to gene3
458      */
459     List<SequenceFeature> geneFeatures = al.findName("gene3")
460             .getSequenceFeatures();
461     SequenceFeatures.sortFeatures(geneFeatures, true);
462     assertEquals(geneFeatures.size(), 7);
463     SequenceFeature sf = geneFeatures.get(0);
464     assertEquals(sf.getBegin(), 1);
465     assertEquals(sf.getEnd(), 1);
466     assertEquals(sf.getScore(), 0.1f, DELTA);
467     assertEquals(sf.getValue("alleles"), "C,A");
468     // gene features include Consequence for all transcripts
469     Map map = (Map) sf.getValue("CSQ");
470     assertEquals(map.size(), 9);
471
472     sf = geneFeatures.get(1);
473     assertEquals(sf.getBegin(), 5);
474     assertEquals(sf.getEnd(), 5);
475     assertEquals(sf.getScore(), 0.2f, DELTA);
476     assertEquals(sf.getValue("alleles"), "C,T");
477     map = (Map) sf.getValue("CSQ");
478     assertEquals(map.size(), 9);
479
480     sf = geneFeatures.get(2);
481     assertEquals(sf.getBegin(), 9);
482     assertEquals(sf.getEnd(), 11); // deletion over 3 positions
483     assertEquals(sf.getScore(), 0.3f, DELTA);
484     assertEquals(sf.getValue("alleles"), "CGG,C");
485     map = (Map) sf.getValue("CSQ");
486     assertEquals(map.size(), 9);
487
488     sf = geneFeatures.get(3);
489     assertEquals(sf.getBegin(), 13);
490     assertEquals(sf.getEnd(), 13);
491     assertEquals(sf.getScore(), 0.5f, DELTA);
492     assertEquals(sf.getValue("alleles"), "C,T");
493     map = (Map) sf.getValue("CSQ");
494     assertEquals(map.size(), 9);
495
496     sf = geneFeatures.get(4);
497     assertEquals(sf.getBegin(), 13);
498     assertEquals(sf.getEnd(), 13);
499     assertEquals(sf.getScore(), 0.4f, DELTA);
500     assertEquals(sf.getValue("alleles"), "C,G");
501     map = (Map) sf.getValue("CSQ");
502     assertEquals(map.size(), 9);
503
504     sf = geneFeatures.get(5);
505     assertEquals(sf.getBegin(), 17);
506     assertEquals(sf.getEnd(), 17);
507     assertEquals(sf.getScore(), 0.7f, DELTA);
508     assertEquals(sf.getValue("alleles"), "A,G");
509     map = (Map) sf.getValue("CSQ");
510     assertEquals(map.size(), 9);
511
512     sf = geneFeatures.get(6);
513     assertEquals(sf.getBegin(), 17);
514     assertEquals(sf.getEnd(), 17); // insertion
515     assertEquals(sf.getScore(), 0.6f, DELTA);
516     assertEquals(sf.getValue("alleles"), "A,AC");
517     map = (Map) sf.getValue("CSQ");
518     assertEquals(map.size(), 9);
519
520     /*
521      * verify variant feature(s) added to transcript3
522      * at columns 5 (1), 17 (2), positions 3, 11
523      * note the deletion at columns 9-11 is not transferred since col 11
524      * has no mapping to transcript 3 
525      */
526     List<SequenceFeature> transcriptFeatures = al.findName("transcript3")
527             .getSequenceFeatures();
528     SequenceFeatures.sortFeatures(transcriptFeatures, true);
529     assertEquals(transcriptFeatures.size(), 3);
530     sf = transcriptFeatures.get(0);
531     assertEquals(sf.getBegin(), 3);
532     assertEquals(sf.getEnd(), 3);
533     assertEquals(sf.getScore(), 0.2f, DELTA);
534     assertEquals(sf.getValue("alleles"), "C,T");
535     // transcript features only have Consequence for that transcripts
536     map = (Map) sf.getValue("CSQ");
537     assertEquals(map.size(), 9);
538     assertEquals(sf.getValueAsString("CSQ", "Feature"), "transcript3");
539
540     sf = transcriptFeatures.get(1);
541     assertEquals(sf.getBegin(), 11);
542     assertEquals(sf.getEnd(), 11);
543     assertEquals(sf.getScore(), 0.7f, DELTA);
544     assertEquals(sf.getValue("alleles"), "A,G");
545     assertEquals(map.size(), 9);
546     assertEquals(sf.getValueAsString("CSQ", "Feature"), "transcript3");
547
548     sf = transcriptFeatures.get(2);
549     assertEquals(sf.getBegin(), 11);
550     assertEquals(sf.getEnd(), 11);
551     assertEquals(sf.getScore(), 0.6f, DELTA);
552     assertEquals(sf.getValue("alleles"), "A,AC");
553     assertEquals(map.size(), 9);
554     assertEquals(sf.getValueAsString("CSQ", "Feature"), "transcript3");
555
556     /*
557      * verify variants computed on protein product for transcript3
558      * peptide is SWRECD
559      * codon variants are AGC/AGT position 1 which is synonymous
560      * and GAG/GGG which is E/G in position 4
561      * the insertion variant is not transferred to the peptide
562      */
563     DBRefEntry[] dbRefs = al.findName("transcript3").getDBRefs();
564     SequenceI peptide = null;
565     for (DBRefEntry dbref : dbRefs)
566     {
567       if (dbref.getMap().getMap().getFromRatio() == 3)
568       {
569         peptide = dbref.getMap().getTo();
570       }
571     }
572     List<SequenceFeature> proteinFeatures = peptide.getSequenceFeatures();
573     assertEquals(proteinFeatures.size(), 1);
574     sf = proteinFeatures.get(0);
575     assertEquals(sf.getFeatureGroup(), "VCF");
576     assertEquals(sf.getBegin(), 4);
577     assertEquals(sf.getEnd(), 4);
578     assertEquals(sf.getType(), SequenceOntologyI.SEQUENCE_VARIANT);
579     assertEquals(sf.getDescription(), "p.Glu4Gly");
580
581     /*
582      * verify variant feature(s) added to transcript4
583      * at columns 13 (2) and 17 (2), positions 7 and 11
584      */
585     transcriptFeatures = al.findName("transcript4").getSequenceFeatures();
586     SequenceFeatures.sortFeatures(transcriptFeatures, true);
587     assertEquals(transcriptFeatures.size(), 4);
588     sf = transcriptFeatures.get(0);
589     assertEquals(sf.getBegin(), 7);
590     assertEquals(sf.getEnd(), 7);
591     assertEquals(sf.getScore(), 0.5f, DELTA);
592     assertEquals(sf.getValue("alleles"), "C,T");
593     assertEquals(map.size(), 9);
594     assertEquals(sf.getValueAsString("CSQ", "Feature"), "transcript4");
595
596     sf = transcriptFeatures.get(1);
597     assertEquals(sf.getBegin(), 7);
598     assertEquals(sf.getEnd(), 7);
599     assertEquals(sf.getScore(), 0.4f, DELTA);
600     assertEquals(sf.getValue("alleles"), "C,G");
601     assertEquals(map.size(), 9);
602     assertEquals(sf.getValueAsString("CSQ", "Feature"), "transcript4");
603
604     sf = transcriptFeatures.get(2);
605     assertEquals(sf.getBegin(), 11);
606     assertEquals(sf.getEnd(), 11);
607     assertEquals(sf.getScore(), 0.7f, DELTA);
608     assertEquals(sf.getValue("alleles"), "A,G");
609     assertEquals(map.size(), 9);
610     assertEquals(sf.getValueAsString("CSQ", "Feature"), "transcript4");
611
612     sf = transcriptFeatures.get(3);
613     assertEquals(sf.getBegin(), 11);
614     assertEquals(sf.getEnd(), 11);
615     assertEquals(sf.getScore(), 0.6f, DELTA);
616     assertEquals(sf.getValue("alleles"), "A,AC");
617     assertEquals(map.size(), 9);
618     assertEquals(sf.getValueAsString("CSQ", "Feature"), "transcript4");
619   }
620
621   /**
622    * A test that demonstrates loading a contig sequence from an indexed sequence
623    * database which is the reference for a VCF file
624    * 
625    * @throws IOException
626    */
627   @Test(groups = "Functional")
628   public void testLoadVCFContig() throws IOException
629   {
630     VCFLoader loader = new VCFLoader(
631             "test/jalview/io/vcf/testVcf2.vcf");
632
633     SequenceI seq = loader.loadVCFContig("contig123");
634     assertEquals(seq.getLength(), 15);
635     assertEquals(seq.getSequenceAsString(), "AAAAACCCCCGGGGG");
636     List<SequenceFeature> features = seq.getSequenceFeatures();
637     SequenceFeatures.sortFeatures(features, true);
638     assertEquals(features.size(), 2);
639     SequenceFeature sf = features.get(0);
640     assertEquals(sf.getBegin(), 8);
641     assertEquals(sf.getEnd(), 8);
642     assertEquals(sf.getDescription(), "C,A");
643     sf = features.get(1);
644     assertEquals(sf.getBegin(), 12);
645     assertEquals(sf.getEnd(), 12);
646     assertEquals(sf.getDescription(), "G,T");
647
648     seq = loader.loadVCFContig("contig789");
649     assertEquals(seq.getLength(), 25);
650     assertEquals(seq.getSequenceAsString(), "GGGGGTTTTTAAAAACCCCCGGGGG");
651     features = seq.getSequenceFeatures();
652     SequenceFeatures.sortFeatures(features, true);
653     assertEquals(features.size(), 2);
654     sf = features.get(0);
655     assertEquals(sf.getBegin(), 2);
656     assertEquals(sf.getEnd(), 2);
657     assertEquals(sf.getDescription(), "G,T");
658     sf = features.get(1);
659     assertEquals(sf.getBegin(), 21);
660     assertEquals(sf.getEnd(), 21);
661     assertEquals(sf.getDescription(), "G,A");
662
663     seq = loader.loadVCFContig("contig456");
664     assertEquals(seq.getLength(), 20);
665     assertEquals(seq.getSequenceAsString(), "CCCCCGGGGGTTTTTAAAAA");
666     features = seq.getSequenceFeatures();
667     SequenceFeatures.sortFeatures(features, true);
668     assertEquals(features.size(), 1);
669     sf = features.get(0);
670     assertEquals(sf.getBegin(), 15);
671     assertEquals(sf.getEnd(), 15);
672     assertEquals(sf.getDescription(), "T,C");
673   }
674 }