test/jalview/ext/htsjdk/VCFReaderTest.java

   1 package jalview.ext.htsjdk;
   2
   3 import static org.testng.Assert.assertEquals;
   4 import static org.testng.Assert.assertFalse;
   5 import static org.testng.Assert.assertTrue;
   6 import htsjdk.samtools.util.CloseableIterator;
   7 import htsjdk.variant.variantcontext.Allele;
   8 import htsjdk.variant.variantcontext.VariantContext;
   9
  10 import java.io.File;
  11 import java.io.IOException;
  12 import java.io.PrintWriter;
  13 import java.util.List;
  14
  15 import org.testng.annotations.Test;
  16
  17 public class VCFReaderTest
  18 {
  19   private static final String[] VCF = new String[] {
  20       "##fileformat=VCFv4.2",
  21       "#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO",
  22       "20\t3\t.\tC\tG\t.\tPASS\tDP=100", // SNP C/G
  23       "20\t7\t.\tG\tGA\t.\tPASS\tDP=100", // insertion G/GA
  24       "18\t2\t.\tACG\tA\t.\tPASS\tDP=100" }; // deletion ACG/A
  25
  26   // gnomAD exome variant dataset
  27   private static final String VCF_PATH = "/Volumes/gjb/smacgowan/NOBACK/resources/gnomad/gnomad.exomes.r2.0.1.sites.vcf.gz";
  28
  29   /**
  30    * A test to exercise some basic functionality of the htsjdk VCF reader
  31    *
  32    * @throws IOException
  33    */
  34   @Test(groups = "Functional")
  35   public void testReadVcf_plain() throws IOException
  36   {
  37     File f = writeVcfFile();
  38     VCFReader reader = new VCFReader(f.getAbsolutePath());
  39     CloseableIterator<VariantContext> variants = reader.iterator();
  40
  41     /*
  42      * SNP C/G variant
  43      */
  44     VariantContext vc = variants.next();
  45     assertTrue(vc.isSNP());
  46     Allele ref = vc.getReference();
  47     assertEquals(ref.getBaseString(), "C");
  48     List<Allele> alleles = vc.getAlleles();
  49     assertEquals(alleles.size(), 2);
  50     assertTrue(alleles.get(0).isReference());
  51     assertEquals(alleles.get(0).getBaseString(), "C");
  52     assertFalse(alleles.get(1).isReference());
  53     assertEquals(alleles.get(1).getBaseString(), "G");
  54
  55     /*
  56      * Insertion G -> GA
  57      */
  58     vc = variants.next();
  59     assertFalse(vc.isSNP());
  60     assertTrue(vc.isSimpleInsertion());
  61     ref = vc.getReference();
  62     assertEquals(ref.getBaseString(), "G");
  63     alleles = vc.getAlleles();
  64     assertEquals(alleles.size(), 2);
  65     assertTrue(alleles.get(0).isReference());
  66     assertEquals(alleles.get(0).getBaseString(), "G");
  67     assertFalse(alleles.get(1).isReference());
  68     assertEquals(alleles.get(1).getBaseString(), "GA");
  69
  70     /*
  71      * Deletion ACG -> A
  72      */
  73     vc = variants.next();
  74     assertFalse(vc.isSNP());
  75     assertTrue(vc.isSimpleDeletion());
  76     ref = vc.getReference();
  77     assertEquals(ref.getBaseString(), "ACG");
  78     alleles = vc.getAlleles();
  79     assertEquals(alleles.size(), 2);
  80     assertTrue(alleles.get(0).isReference());
  81     assertEquals(alleles.get(0).getBaseString(), "ACG");
  82     assertFalse(alleles.get(1).isReference());
  83     assertEquals(alleles.get(1).getBaseString(), "A");
  84
  85     assertFalse(variants.hasNext());
  86
  87     variants.close();
  88     reader.close();
  89   }
  90
  91   /**
  92    * Creates a temporary file to be read by the htsjdk VCF reader
  93    *
  94    * @return
  95    * @throws IOException
  96    */
  97   protected File writeVcfFile() throws IOException
  98   {
  99     File f = File.createTempFile("Test", "vcf");
 100     f.deleteOnExit();
 101     PrintWriter pw = new PrintWriter(f);
 102     for (String vcfLine : VCF) {
 103       pw.println(vcfLine);
 104     }
 105     pw.close();
 106     return f;
 107   }
 108
 109   // "https://storage.cloud.google.com/gnomad-public/release/2.0.1/vcf/exomes/gnomad.exomes.r2.0.1.sites.vcf.gz";
 110
 111   /**
 112    * A 'test' that demonstrates querying an indexed VCF file for features in a
 113    * specified interval
 114    *
 115    * @throws IOException
 116    */
 117   @Test
 118   public void testQuery_indexed() throws IOException
 119   {
 120     /*
 121      * if not specified, assumes index file is filename.tbi
 122      */
 123     VCFReader reader = new VCFReader(VCF_PATH);
 124
 125     /*
 126      * gene NMT1 (human) is on chromosome 17
 127      * GCHR38 (Ensembl): 45051610-45109016
 128      * GCHR37 (gnoMAD): 43128978-43186384
 129      * CDS begins at offset 9720, first CDS variant at offset 9724
 130      */
 131     CloseableIterator<VariantContext> features = reader.query("17",
 132             43128978 + 9724, 43128978 + 9734); // first 11 CDS positions
 133
 134     assertEquals(features.next().getStart(), 43138702);
 135     assertEquals(features.next().getStart(), 43138704);
 136     assertEquals(features.next().getStart(), 43138707);
 137     assertEquals(features.next().getStart(), 43138708);
 138     assertEquals(features.next().getStart(), 43138710);
 139     assertEquals(features.next().getStart(), 43138711);
 140     assertFalse(features.hasNext());
 141
 142     features.close();
 143     reader.close();
 144   }
 145 }