/**
* 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.
+ * <p>
+ * Client code should call close() on the iterator when finished with it.
*
* @param chrom
* the chromosome to query
public CloseableIterator<VariantContext> 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<VariantContext> queryUnindexed(
+ final String chrom, final int start, final int end)
+ {
+ final CloseableIterator<VariantContext> it = reader.iterator();
+
+ return new CloseableIterator<VariantContext>()
+ {
+ 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();
+ }
+ };
}
/**
{
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;
+ }
}
// "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
*/
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<VariantContext> 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<Allele> 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();
+ }
}
--- /dev/null
+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=<ID=AF,Number=A,Type=Float,Description=\"Allele Frequency, for each ALT allele, in the same order as listed\">",
+ "##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<SequenceFeature> 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<SequenceFeature> 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<SequenceFeature> 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<int[]> toRanges = Arrays.asList(to);
+ SequenceI gene = alignment.getSequenceAt(0);
+ List<int[]> 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
+ }
+}
git://source.jalview.org / jalview.git / commitdiff