This is a test version of Biostars. For the public version, visit https://www.biostars.org.
fast way to get last position in a large indexed VCF?

I want to get first and last position in a VCF for chunking purposes

First position is easy to get:

bcftools query -f '%CHROM\t%POS\n' my_vcf  | head -n 1

I realize I could stream the whole vcf to find last position, but I was hoping there was some faster way to get it.

The only thing I can think of is some kind of binary search like approach, but there has to be a better way

bcftools
• 5,076 views • link
updated by bgulko2 • bgulko2 • written by curious curious

Maybe I didn't understand the question, but if the VCF is sorted, why not just use tail?

• 0 views • link
written by liorglic liorglic

getting to the tail of a 10GB compressed file will take quite a bit and means unpacking the entire file.

As for the OP, perhaps reading the index file with a custom program might be able to tell you the last coordinate then you access that directly.

• 0 views • link
written by Istvan Albert Istvan Albert

Do you any idea how to do this, I looked for this approach even before positing this. I cant cat or zcat the index and get anything readable and I cant find any helpful info with 'bcftools index' or 'tabix' options

• 0 views • link
written by curious curious

when poking around this way I found this, which is already doing my chunking to boot:

from sgkit.io.vcf import partition_into_regions
partition_into_regions("CEUTrio.20.21.gatk3.4.g.vcf.bgz", num_parts=10)

['20:1-10108928', '20:10108929-10207232', '20:10207233-', '21:1-10027008', '21:10027009-10043392', '21:10043393-10108928', '21:10108929-10141696', '21:10141697-10174464', '21:10174465-10190848', '21:10190849-10207232', '21:10207233-']

https://pystatgen.github.io/sgkit/latest/vcf.html#partitioning

• 0 views • link
written by curious curious

I have been playing with this problem a bit in Python, and I concluded that you would need a bgzip aware solution for it to work. It does not seem like sgkit is such a tool though.

When using the builtin gzip stream Python will proceed to uncompress the entire stream thus, even if you figured out the correct offset it won't be able to jump to that location without unzipping everything to reach the end.

Alas, I was unable to locate a bgzip library in Python that also allows seeking a file to an offset.

• 0 views • link
written by Istvan Albert Istvan Albert

That's an interesting question but I'm not sure a tabix index is relevant to what position happens to be at the end of the file. Unless it's position sorted, the vcf could end with chr1:123, for instance, instead of chrM:16569. The tabix index doesn't really care.

• 0 views • link
written by Jeremy Leipzig Jeremy Leipzig

for tabix to work, the file has to be position sorted, and grouped by chromosome (usually we sort by both):

(grep "^#" in.gff; grep -v "^#" in.gff | sort -t"`printf '\t'`" -k1,1 -k4,4n) | bgzip > sorted.gff.gz;

tabix -p gff sorted.gff.gz;

tabix sorted.gff.gz chr1:10,000,000-20,000,000;

as seen in: http://www.htslib.org/doc/tabix.html

• 0 views • link
written by Istvan Albert Istvan Albert

I am also looking for this. GATK can dump an index via DumpTabixIndex, but it is not quite clear how the result this maps to contig positions. Did you find an answer to this that avoids traversing the entire VCF file?

• 0 views • link
written by bgulko2 • bgulko2

why do you want this? if you describe the purpose it might help people to answer

• 0 views • link
written by cmdcolin cmdcolin

Purpose is to get the position of the last variant in a compressed, indexed, VCF file without streaming the entire file (e.g. without using "tail" ).

• 0 views • link
written by bgulko2 • bgulko2

2 answers

If the file is uncompressed (or if you find a bgzip compatible stream.seek for your programming language) you can jump near the end of the file quite efficiently:

import os

fname = "subset.vcf"
stream = open(fname)
size = os.path.getsize(fname)

# Move back a few kb
shift = 10 * 1024
stream.seek(size-shift)

for line in stream:
    pass

print(line.split())

it will be pretty fast:

 time python wind.py
['NC_000001.10', '25202832', 'rs1640276932', 'TGGGGTGGGGCGGGGCGGGGAGGGGC', 'T', '.', '.', 'RS=1640276932;SSR=0;VC=INDEL;GNO;FREQ=dbGaP_PopFreq:1,0']

real    0m0.043s
user    0m0.016s
sys     0m0.016s

versus streaming the whole file:

time cat subset.vcf| tail -1
NC_000001.10    25202832        rs1640276932    TGGGGTGGGGCGGGGCGGGGAGGGGC      T       .       .       RS=1640276932;SSR=0;VC=INDEL;GNO;FREQ=dbGaP_PopFreq:1,0

real    0m1.492s
user    0m0.313s
sys     0m1.578s
• 0 views • link
written by Istvan Albert Istvan Albert

How does running bctools query -f'%CHROM %POS\n' | tail -n1 compare?

• 0 views • link
written by 4galaxy77 4galaxy77 
time bcftools query -f'%CHROM %POS\n' subset.vcf | tail -n1
NC_000001.10 25202832

real    0m29.221s
user    0m6.281s
sys     0m23.922s

for the record, the file has 9 million variants

cat subset.vcf | wc -l
9199371
• 0 views • link
written by Istvan Albert Istvan Albert

There is a utility called bgztail available at https://github.com/circulosmeos/bgztail which will provide the end of a b-gzipped vcf via indexed lookup without decompressing the entire file. It did take me awhile to find. It also does not require an index (csi/tbi/gzi). So, for example:

time bgztail -2 SomeFile.vcf.gz | cut -f1-2
chrX    156029827
chrX    156030794

real    0m0.042s
user    0m0.029s
sys     0m0.016s

vs

time bcftools view  SomeFile.vcf.gz | tail -2 | cut -f1-2
chrX    156029827
chrX    156030794

real    0m4.623s
user    0m4.773s
sys     0m0.675s
• 0 views • link
written by bgulko2 • bgulko2

Log in to answer this question.

More posts like this