From 2958419af2db03ef23cfd96acb66296e080b978a Mon Sep 17 00:00:00 2001 From: gmungoc Date: Mon, 25 Sep 2017 14:42:40 +0100 Subject: [PATCH] JAL-2738 query unindexed file, unit test for VCFLoader --- src/jalview/ext/htsjdk/VCFReader.java | 120 +++++++++++++++++++- test/jalview/ext/htsjdk/VCFReaderTest.java | 43 +++++++- test/jalview/io/vcf/VCFLoaderTest.java | 164 ++++++++++++++++++++++++++++ 3 files changed, 322 insertions(+), 5 deletions(-) create mode 100644 test/jalview/io/vcf/VCFLoaderTest.java diff --git a/src/jalview/ext/htsjdk/VCFReader.java b/src/jalview/ext/htsjdk/VCFReader.java index c5e09e0..14c057f 100644 --- a/src/jalview/ext/htsjdk/VCFReader.java +++ b/src/jalview/ext/htsjdk/VCFReader.java @@ -72,9 +72,10 @@ public class VCFReader implements Closeable, Iterable /** * Queries for records overlapping the region specified. Note that this method - * requires a VCF file with an associated index. If no index exists a - * TribbleException will be thrown. Client code should call close() on the - * iterator when finished with it. + * is performant if the VCF file is indexed, and may be very slow if it is + * not. + *

+ * Client code should call close() on the iterator when finished with it. * * @param chrom * the chromosome to query @@ -87,7 +88,107 @@ public class VCFReader implements Closeable, Iterable public CloseableIterator query(final String chrom, final int start, final int end) { - return reader == null ? null : reader.query(chrom, start, end); + if (reader == null) { + return null; + } + if (indexed) + { + return reader.query(chrom, start, end); + } + else + { + return queryUnindexed(chrom, start, end); + } + } + + /** + * Returns an iterator over variant records read from a flat file which + * overlap the specified chromosomal positions. Call close() on the iterator + * when finished with it! + * + * @param chrom + * @param start + * @param end + * @return + */ + protected CloseableIterator queryUnindexed( + final String chrom, final int start, final int end) + { + final CloseableIterator it = reader.iterator(); + + return new CloseableIterator() + { + boolean atEnd = false; + + // prime look-ahead buffer with next matching record + private VariantContext next = findNext(); + + private VariantContext findNext() + { + if (atEnd) + { + return null; + } + VariantContext variant = null; + while (it.hasNext()) + { + variant = it.next(); + int vstart = variant.getStart(); + + if (vstart > end) + { + atEnd = true; + close(); + return null; + } + + int vend = variant.getEnd(); + // todo what is the undeprecated way to get + // the chromosome for the variant? + if (chrom.equals(variant.getChr()) && (vstart <= end) + && (vend >= start)) + { + return variant; + } + } + return null; + } + + @Override + public boolean hasNext() + { + boolean hasNext = !atEnd && (next != null); + if (!hasNext) + { + close(); + } + return hasNext; + } + + @Override + public VariantContext next() + { + /* + * return the next match, and then re-prime + * it with the following one (if any) + */ + VariantContext temp = next; + next = findNext(); + return temp; + } + + @Override + public void remove() + { + // not implemented + } + + @Override + public void close() + { + it.close(); + } + }; } /** @@ -99,4 +200,15 @@ public class VCFReader implements Closeable, Iterable { return reader == null ? null : reader.getFileHeader(); } + + /** + * Answers true if we are processing a tab-indexed VCF file, false if it is a + * plain text (uncompressed) file. + * + * @return + */ + public boolean isIndex() + { + return indexed; + } } diff --git a/test/jalview/ext/htsjdk/VCFReaderTest.java b/test/jalview/ext/htsjdk/VCFReaderTest.java index 29d752c..bf617ae 100644 --- a/test/jalview/ext/htsjdk/VCFReaderTest.java +++ b/test/jalview/ext/htsjdk/VCFReaderTest.java @@ -29,7 +29,8 @@ public class VCFReaderTest // "https://storage.cloud.google.com/gnomad-public/release/2.0.1/vcf/exomes/gnomad.exomes.r2.0.1.sites.vcf.gz"; /** - * A test to exercise some basic functionality of the htsjdk VCF reader + * A test to exercise some basic functionality of the htsjdk VCF reader, + * reading from a non-index VCF file * * @throws IOException */ @@ -156,4 +157,44 @@ public class VCFReaderTest System.out.println(next.toString()); return next.getStart(); } + + // "https://storage.cloud.google.com/gnomad-public/release/2.0.1/vcf/exomes/gnomad.exomes.r2.0.1.sites.vcf.gz"; + + /** + * Test the query method that wraps a non-indexed VCF file + * + * @throws IOException + */ + @Test(groups = "Functional") + public void testQuery_plain() throws IOException + { + File f = writeVcfFile(); + VCFReader reader = new VCFReader(f.getAbsolutePath()); + + /* + * query for overlap of 5-8 - should find variant at 7 + */ + CloseableIterator variants = reader.query("20", 5, 8); + + /* + * INDEL G/GA variant + */ + VariantContext vc = variants.next(); + assertTrue(vc.isIndel()); + assertEquals(vc.getStart(), 7); + assertEquals(vc.getEnd(), 7); + Allele ref = vc.getReference(); + assertEquals(ref.getBaseString(), "G"); + List alleles = vc.getAlleles(); + assertEquals(alleles.size(), 2); + assertTrue(alleles.get(0).isReference()); + assertEquals(alleles.get(0).getBaseString(), "G"); + assertFalse(alleles.get(1).isReference()); + assertEquals(alleles.get(1).getBaseString(), "GA"); + + assertFalse(variants.hasNext()); + + variants.close(); + reader.close(); + } } diff --git a/test/jalview/io/vcf/VCFLoaderTest.java b/test/jalview/io/vcf/VCFLoaderTest.java new file mode 100644 index 0000000..7319123 --- /dev/null +++ b/test/jalview/io/vcf/VCFLoaderTest.java @@ -0,0 +1,164 @@ +package jalview.io.vcf; + +import static org.testng.Assert.assertEquals; + +import jalview.datamodel.AlignmentI; +import jalview.datamodel.DBRefEntry; +import jalview.datamodel.GeneLoci; +import jalview.datamodel.Mapping; +import jalview.datamodel.Sequence; +import jalview.datamodel.SequenceFeature; +import jalview.datamodel.SequenceI; +import jalview.gui.AlignFrame; +import jalview.io.DataSourceType; +import jalview.io.FileLoader; +import jalview.io.gff.Gff3Helper; +import jalview.io.gff.SequenceOntologyI; +import jalview.util.MapList; + +import java.io.File; +import java.io.IOException; +import java.io.PrintWriter; +import java.util.Arrays; +import java.util.List; + +import org.testng.annotations.Test; + +public class VCFLoaderTest +{ + // columns 9717- of gene P30419 from Ensembl (modified) + private static final String FASTA = ">ENSG00000136448/1-25 chromosome:GRCh38:17:45051610:45051634:1\n" + + "CAAGCTGGCGGACGAGAGTGTGACA\n" + // and a 'made up' mini-transcript with two exons + + ">ENST00000592782/1-18\n--AGCTGGCG----AGAGTGTGAC-\n"; + + private static final String[] VCF = { "##fileformat=VCFv4.2", + "##INFO=", + "##reference=GRCh38", + "#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO", + // SNP A/T in position 2 of gene sequence (precedes transcript) + "17\t45051611\t.\tA\tT\t1666.64\tRF\tAC=15;AF=5.08130e-03", + // SNP G/C in position 4 of gene sequence, position 2 of transcript + // this is a mixed variant, the insertion G/GA is not transferred + "17\t45051613\t.\tG\tGA,C\t1666.64\tRF\tAC=15;AF=3.08130e-03" }; + + @Test(groups = "Functional") + public void testLoadVCF() throws IOException + { + AlignmentI al = buildAlignment(); + VCFLoader loader = new VCFLoader(al); + + File f = makeVcf(); + + loader.loadVCF(f.getPath(), null); + + /* + * verify variant feature(s) added to gene + */ + List geneFeatures = al.getSequenceAt(0).findFeatures( + 2, 2); + assertEquals(geneFeatures.size(), 1); + SequenceFeature sf = geneFeatures.get(0); + assertEquals(sf.getFeatureGroup(), "VCF"); + assertEquals(sf.getType(), SequenceOntologyI.SEQUENCE_VARIANT); + assertEquals(sf.getScore(), 5.08130e-03, 0.000001f); + assertEquals("A,T", sf.getValue(Gff3Helper.ALLELES)); + + /* + * verify variant feature(s) added to transcript + */ + List transcriptFeatures = al.getSequenceAt(1) + .findFeatures(4, 4); + assertEquals(transcriptFeatures.size(), 1); + sf = transcriptFeatures.get(0); + assertEquals(sf.getFeatureGroup(), "VCF"); + assertEquals(sf.getType(), SequenceOntologyI.SEQUENCE_VARIANT); + assertEquals(sf.getScore(), 3.08130e-03, 0.000001f); + assertEquals("G,C", sf.getValue(Gff3Helper.ALLELES)); + + /* + * verify variant feature(s) computed and added to protein + * first codon AGC varies to ACC giving S/T + */ + SequenceI peptide = al.getSequenceAt(1) + .getDBRefs()[0].getMap().getTo(); + List proteinFeatures = peptide.findFeatures(1, 6); + assertEquals(proteinFeatures.size(), 1); + sf = proteinFeatures.get(0); + assertEquals(sf.getFeatureGroup(), "VCF"); + assertEquals(sf.getBegin(), 2); + assertEquals(sf.getEnd(), 2); + assertEquals(sf.getType(), SequenceOntologyI.SEQUENCE_VARIANT); + assertEquals(sf.getDescription(), "p.Ser1Thr"); + } + + private File makeVcf() throws IOException + { + File f = File.createTempFile("Test", ".vcf"); + f.deleteOnExit(); + PrintWriter pw = new PrintWriter(f); + for (String vcfLine : VCF) + { + pw.println(vcfLine); + } + pw.close(); + return f; + } + + /** + * Make a simple alignment with one 'gene' and one 'transcript' + * + * @return + */ + private AlignmentI buildAlignment() + { + AlignFrame af = new FileLoader().LoadFileWaitTillLoaded(FASTA, + DataSourceType.PASTE); + + /* + * map gene sequence to chromosome (normally done when the sequence is fetched + * from Ensembl and transcripts computed) + */ + AlignmentI alignment = af.getViewport().getAlignment(); + int[][] to = new int[][] { new int[] { 45051610, 45051634 } }; + List toRanges = Arrays.asList(to); + SequenceI gene = alignment.getSequenceAt(0); + List fromRanges = Arrays.asList(new int[][] { new int[] { + gene.getStart(), gene.getEnd() } }); + ((Sequence) gene).setGeneLoci(new GeneLoci("human", "GRCh38", "17", + new MapList(fromRanges, toRanges, 1, 1))); + + /* + * map 'transcript' to chromosome via 'gene' + * transcript/1-18 is gene/3-10,15-24 + * which is chromosome 45051612-45051619,45051624-45051633 + */ + to = new int[][] { new int[] { 45051612, 45051619 }, + new int[] { 45051624, 45051633 } }; + toRanges = Arrays.asList(to); + SequenceI transcript = alignment.getSequenceAt(1); + fromRanges = Arrays.asList(new int[][] { new int[] { + transcript.getStart(), transcript.getEnd() } }); + ((Sequence) transcript).setGeneLoci(new GeneLoci("human", "GRCh38", + "17", new MapList(fromRanges, toRanges, 1, 1))); + + /* + * add a protein product as a DBRef on the transcript + */ + SequenceI peptide = new Sequence("ENSP001", "SWRECD"); + MapList mapList = new MapList(new int[] { 1, 18 }, new int[] { 1, 6 }, + 3, 1); + Mapping map = new Mapping(peptide, mapList); + DBRefEntry product = new DBRefEntry("", "", "ENSP001", map); + transcript.addDBRef(product); + + return alignment; + } + + @Test(groups = "Functional") + public void testLoadVCF_reverseStrand() throws IOException + { + // TODO a test with reverse strand mapping of + // gene and transcript to chromosome + } +} -- 1.7.10.2