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   // "https://storage.cloud.google.com/gnomad-public/release/2.0.1/vcf/exomes/gnomad.exomes.r2.0.1.sites.vcf.gz";
  30
  31   /**
  32    * A test to exercise some basic functionality of the htsjdk VCF reader
  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       pw.println(vcfLine);
 106     }
 107     pw.close();
 108     return f;
 109   }
 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(printNext(features), 43138702);
 135     assertEquals(printNext(features), 43138704);
 136     assertEquals(printNext(features), 43138707);
 137     assertEquals(printNext(features), 43138708);
 138     assertEquals(printNext(features), 43138710);
 139     assertEquals(printNext(features), 43138711);
 140     assertFalse(features.hasNext());
 141
 142     features.close();
 143     reader.close();
 144   }
 145
 146   /**
 147    * Prints the toString value of the next variant, and returns its start
 148    * location
 149    *
 150    * @param features
 151    * @return
 152    */
 153   protected int printNext(CloseableIterator<VariantContext> features)
 154   {
 155     VariantContext next = features.next();
 156     System.out.println(next.toString());
 157     return next.getStart();
 158   }
 159 }