test/jalview/io/vcf/VCFLoaderTest.java

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