JAL-4366 test reconstruction of peptide alignment given 3di alignment and peptide...
[jalview.git] / test / jalview / ext / htsjdk / VCFReaderTest.java
1 /*
2  * Jalview - A Sequence Alignment Editor and Viewer ($$Version-Rel$$)
3  * Copyright (C) $$Year-Rel$$ The Jalview Authors
4  * 
5  * This file is part of Jalview.
6  * 
7  * Jalview is free software: you can redistribute it and/or
8  * modify it under the terms of the GNU General Public License 
9  * as published by the Free Software Foundation, either version 3
10  * of the License, or (at your option) any later version.
11  *  
12  * Jalview is distributed in the hope that it will be useful, but 
13  * WITHOUT ANY WARRANTY; without even the implied warranty 
14  * of MERCHANTABILITY or FITNESS FOR A PARTICULAR 
15  * PURPOSE.  See the GNU General Public License for more details.
16  * 
17  * You should have received a copy of the GNU General Public License
18  * along with Jalview.  If not, see <http://www.gnu.org/licenses/>.
19  * The Jalview Authors are detailed in the 'AUTHORS' file.
20  */
21 package jalview.ext.htsjdk;
22
23 import static org.testng.Assert.assertEquals;
24 import static org.testng.Assert.assertFalse;
25 import static org.testng.Assert.assertTrue;
26 import htsjdk.samtools.util.CloseableIterator;
27 import htsjdk.variant.variantcontext.Allele;
28 import htsjdk.variant.variantcontext.VariantContext;
29
30 import java.io.File;
31 import java.io.IOException;
32 import java.io.PrintWriter;
33 import java.util.List;
34
35 import org.testng.annotations.Test;
36
37 public class VCFReaderTest
38 {
39   private static final String[] VCF = new String[] { "##fileformat=VCFv4.2",
40       "#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO",
41       "20\t3\t.\tC\tG\t.\tPASS\tDP=100", // SNP C/G
42       "20\t7\t.\tG\tGA\t.\tPASS\tDP=100", // insertion G/GA
43       "18\t2\t.\tACG\tA\t.\tPASS\tDP=100" }; // deletion ACG/A
44
45   // gnomAD exome variant dataset
46   private static final String VCF_PATH = "/Volumes/gjb/smacgowan/NOBACK/resources/gnomad/gnomad.exomes.r2.0.1.sites.vcf.gz";
47
48   // "https://storage.cloud.google.com/gnomad-public/release/2.0.1/vcf/exomes/gnomad.exomes.r2.0.1.sites.vcf.gz";
49
50   /**
51    * A test to exercise some basic functionality of the htsjdk VCF reader,
52    * reading from a non-index VCF file
53    * 
54    * @throws IOException
55    */
56   @Test(groups = "Functional")
57   public void testReadVcf_plain() throws IOException
58   {
59     File f = writeVcfFile();
60     VCFReader reader = new VCFReader(f.getAbsolutePath());
61     CloseableIterator<VariantContext> variants = reader.iterator();
62
63     /*
64      * SNP C/G variant
65      */
66     VariantContext vc = variants.next();
67     assertTrue(vc.isSNP());
68     Allele ref = vc.getReference();
69     assertEquals(ref.getBaseString(), "C");
70     List<Allele> alleles = vc.getAlleles();
71     assertEquals(alleles.size(), 2);
72     assertTrue(alleles.get(0).isReference());
73     assertEquals(alleles.get(0).getBaseString(), "C");
74     assertFalse(alleles.get(1).isReference());
75     assertEquals(alleles.get(1).getBaseString(), "G");
76
77     /*
78      * Insertion G -> GA
79      */
80     vc = variants.next();
81     assertFalse(vc.isSNP());
82     assertTrue(vc.isSimpleInsertion());
83     ref = vc.getReference();
84     assertEquals(ref.getBaseString(), "G");
85     alleles = vc.getAlleles();
86     assertEquals(alleles.size(), 2);
87     assertTrue(alleles.get(0).isReference());
88     assertEquals(alleles.get(0).getBaseString(), "G");
89     assertFalse(alleles.get(1).isReference());
90     assertEquals(alleles.get(1).getBaseString(), "GA");
91
92     /*
93      * Deletion ACG -> A
94      */
95     vc = variants.next();
96     assertFalse(vc.isSNP());
97     assertTrue(vc.isSimpleDeletion());
98     ref = vc.getReference();
99     assertEquals(ref.getBaseString(), "ACG");
100     alleles = vc.getAlleles();
101     assertEquals(alleles.size(), 2);
102     assertTrue(alleles.get(0).isReference());
103     assertEquals(alleles.get(0).getBaseString(), "ACG");
104     assertFalse(alleles.get(1).isReference());
105     assertEquals(alleles.get(1).getBaseString(), "A");
106
107     assertFalse(variants.hasNext());
108
109     variants.close();
110     reader.close();
111   }
112
113   /**
114    * Creates a temporary file to be read by the htsjdk VCF reader
115    * 
116    * @return
117    * @throws IOException
118    */
119   protected File writeVcfFile() throws IOException
120   {
121     File f = File.createTempFile("Test", "vcf");
122     f.deleteOnExit();
123     PrintWriter pw = new PrintWriter(f);
124     for (String vcfLine : VCF)
125     {
126       pw.println(vcfLine);
127     }
128     pw.close();
129     return f;
130   }
131
132   /**
133    * A 'test' that demonstrates querying an indexed VCF file for features in a
134    * specified interval
135    * 
136    * @throws IOException
137    */
138   @Test
139   public void testQuery_indexed() throws IOException
140   {
141     /*
142      * if not specified, assumes index file is filename.tbi
143      */
144     VCFReader reader = new VCFReader(VCF_PATH);
145
146     /*
147      * gene NMT1 (human) is on chromosome 17
148      * GCHR38 (Ensembl): 45051610-45109016
149      * GCHR37 (gnoMAD): 43128978-43186384
150      * CDS begins at offset 9720, first CDS variant at offset 9724
151      */
152     CloseableIterator<VariantContext> features = reader.query("17",
153             43128978 + 9724, 43128978 + 9734); // first 11 CDS positions
154
155     assertEquals(printNext(features), 43138702);
156     assertEquals(printNext(features), 43138704);
157     assertEquals(printNext(features), 43138707);
158     assertEquals(printNext(features), 43138708);
159     assertEquals(printNext(features), 43138710);
160     assertEquals(printNext(features), 43138711);
161     assertFalse(features.hasNext());
162
163     features.close();
164     reader.close();
165   }
166
167   /**
168    * Prints the toString value of the next variant, and returns its start
169    * location
170    * 
171    * @param features
172    * @return
173    */
174   protected int printNext(CloseableIterator<VariantContext> features)
175   {
176     VariantContext next = features.next();
177     System.out.println(next.toString());
178     return next.getStart();
179   }
180
181   // "https://storage.cloud.google.com/gnomad-public/release/2.0.1/vcf/exomes/gnomad.exomes.r2.0.1.sites.vcf.gz";
182
183   /**
184    * Test the query method that wraps a non-indexed VCF file
185    * 
186    * @throws IOException
187    */
188   @Test(groups = "Functional")
189   public void testQuery_plain() throws IOException
190   {
191     File f = writeVcfFile();
192     VCFReader reader = new VCFReader(f.getAbsolutePath());
193
194     /*
195      * query for overlap of 5-8 - should find variant at 7
196      */
197     CloseableIterator<VariantContext> variants = reader.query("20", 5, 8);
198
199     /*
200      * INDEL G/GA variant
201      */
202     VariantContext vc = variants.next();
203     assertTrue(vc.isIndel());
204     assertEquals(vc.getStart(), 7);
205     assertEquals(vc.getEnd(), 7);
206     Allele ref = vc.getReference();
207     assertEquals(ref.getBaseString(), "G");
208     List<Allele> alleles = vc.getAlleles();
209     assertEquals(alleles.size(), 2);
210     assertTrue(alleles.get(0).isReference());
211     assertEquals(alleles.get(0).getBaseString(), "G");
212     assertFalse(alleles.get(1).isReference());
213     assertEquals(alleles.get(1).getBaseString(), "GA");
214
215     assertFalse(variants.hasNext());
216
217     variants.close();
218     reader.close();
219   }
220 }