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