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[] { "##fileformat=VCFv4.2",
  20       "#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO",
  21       "20\t3\t.\tC\tG\t.\tPASS\tDP=100", // SNP C/G
  22       "20\t7\t.\tG\tGA\t.\tPASS\tDP=100", // insertion G/GA
  23       "18\t2\t.\tACG\tA\t.\tPASS\tDP=100" }; // deletion ACG/A
  24
  25   // gnomAD exome variant dataset
  26   private static final String VCF_PATH = "/Volumes/gjb/smacgowan/NOBACK/resources/gnomad/gnomad.exomes.r2.0.1.sites.vcf.gz";
  27
  28   // "https://storage.cloud.google.com/gnomad-public/release/2.0.1/vcf/exomes/gnomad.exomes.r2.0.1.sites.vcf.gz";
  29
  30   /**
  31    * A test to exercise some basic functionality of the htsjdk VCF reader,
  32    * reading from a non-index VCF file
  33    *
  34    * @throws IOException
  35    */
  36   @Test(groups = "Functional")
  37   public void testReadVcf_plain() throws IOException
  38   {
  39     File f = writeVcfFile();
  40     VCFReader reader = new VCFReader(f.getAbsolutePath());
  41     CloseableIterator<VariantContext> variants = reader.iterator();
  42
  43     /*
  44      * SNP C/G variant
  45      */
  46     VariantContext vc = variants.next();
  47     assertTrue(vc.isSNP());
  48     Allele ref = vc.getReference();
  49     assertEquals(ref.getBaseString(), "C");
  50     List<Allele> alleles = vc.getAlleles();
  51     assertEquals(alleles.size(), 2);
  52     assertTrue(alleles.get(0).isReference());
  53     assertEquals(alleles.get(0).getBaseString(), "C");
  54     assertFalse(alleles.get(1).isReference());
  55     assertEquals(alleles.get(1).getBaseString(), "G");
  56
  57     /*
  58      * Insertion G -> GA
  59      */
  60     vc = variants.next();
  61     assertFalse(vc.isSNP());
  62     assertTrue(vc.isSimpleInsertion());
  63     ref = vc.getReference();
  64     assertEquals(ref.getBaseString(), "G");
  65     alleles = vc.getAlleles();
  66     assertEquals(alleles.size(), 2);
  67     assertTrue(alleles.get(0).isReference());
  68     assertEquals(alleles.get(0).getBaseString(), "G");
  69     assertFalse(alleles.get(1).isReference());
  70     assertEquals(alleles.get(1).getBaseString(), "GA");
  71
  72     /*
  73      * Deletion ACG -> A
  74      */
  75     vc = variants.next();
  76     assertFalse(vc.isSNP());
  77     assertTrue(vc.isSimpleDeletion());
  78     ref = vc.getReference();
  79     assertEquals(ref.getBaseString(), "ACG");
  80     alleles = vc.getAlleles();
  81     assertEquals(alleles.size(), 2);
  82     assertTrue(alleles.get(0).isReference());
  83     assertEquals(alleles.get(0).getBaseString(), "ACG");
  84     assertFalse(alleles.get(1).isReference());
  85     assertEquals(alleles.get(1).getBaseString(), "A");
  86
  87     assertFalse(variants.hasNext());
  88
  89     variants.close();
  90     reader.close();
  91   }
  92
  93   /**
  94    * Creates a temporary file to be read by the htsjdk VCF reader
  95    *
  96    * @return
  97    * @throws IOException
  98    */
  99   protected File writeVcfFile() throws IOException
 100   {
 101     File f = File.createTempFile("Test", "vcf");
 102     f.deleteOnExit();
 103     PrintWriter pw = new PrintWriter(f);
 104     for (String vcfLine : VCF)
 105     {
 106       pw.println(vcfLine);
 107     }
 108     pw.close();
 109     return f;
 110   }
 111
 112   /**
 113    * A 'test' that demonstrates querying an indexed VCF file for features in a
 114    * specified interval
 115    *
 116    * @throws IOException
 117    */
 118   @Test
 119   public void testQuery_indexed() throws IOException
 120   {
 121     /*
 122      * if not specified, assumes index file is filename.tbi
 123      */
 124     VCFReader reader = new VCFReader(VCF_PATH);
 125
 126     /*
 127      * gene NMT1 (human) is on chromosome 17
 128      * GCHR38 (Ensembl): 45051610-45109016
 129      * GCHR37 (gnoMAD): 43128978-43186384
 130      * CDS begins at offset 9720, first CDS variant at offset 9724
 131      */
 132     CloseableIterator<VariantContext> features = reader.query("17",
 133             43128978 + 9724, 43128978 + 9734); // first 11 CDS positions
 134
 135     assertEquals(printNext(features), 43138702);
 136     assertEquals(printNext(features), 43138704);
 137     assertEquals(printNext(features), 43138707);
 138     assertEquals(printNext(features), 43138708);
 139     assertEquals(printNext(features), 43138710);
 140     assertEquals(printNext(features), 43138711);
 141     assertFalse(features.hasNext());
 142
 143     features.close();
 144     reader.close();
 145   }
 146
 147   /**
 148    * Prints the toString value of the next variant, and returns its start
 149    * location
 150    *
 151    * @param features
 152    * @return
 153    */
 154   protected int printNext(CloseableIterator<VariantContext> features)
 155   {
 156     VariantContext next = features.next();
 157     System.out.println(next.toString());
 158     return next.getStart();
 159   }
 160
 161   // "https://storage.cloud.google.com/gnomad-public/release/2.0.1/vcf/exomes/gnomad.exomes.r2.0.1.sites.vcf.gz";
 162
 163   /**
 164    * Test the query method that wraps a non-indexed VCF file
 165    *
 166    * @throws IOException
 167    */
 168   @Test(groups = "Functional")
 169   public void testQuery_plain() throws IOException
 170   {
 171     File f = writeVcfFile();
 172     VCFReader reader = new VCFReader(f.getAbsolutePath());
 173
 174     /*
 175      * query for overlap of 5-8 - should find variant at 7
 176      */
 177     CloseableIterator<VariantContext> variants = reader.query("20", 5, 8);
 178
 179     /*
 180      * INDEL G/GA variant
 181      */
 182     VariantContext vc = variants.next();
 183     assertTrue(vc.isIndel());
 184     assertEquals(vc.getStart(), 7);
 185     assertEquals(vc.getEnd(), 7);
 186     Allele ref = vc.getReference();
 187     assertEquals(ref.getBaseString(), "G");
 188     List<Allele> alleles = vc.getAlleles();
 189     assertEquals(alleles.size(), 2);
 190     assertTrue(alleles.get(0).isReference());
 191     assertEquals(alleles.get(0).getBaseString(), "G");
 192     assertFalse(alleles.get(1).isReference());
 193     assertEquals(alleles.get(1).getBaseString(), "GA");
 194
 195     assertFalse(variants.hasNext());
 196
 197     variants.close();
 198     reader.close();
 199   }
 200 }