Any fast options to query large VCF bed intervals?

Question

I'm doing some analysis and I need to subset a large VCF file (~8GB gziped) given a bed interval and identify within a subset of rsid.

Unfortunately, both my normal choices to do this analysis (snpSift and bedtools) are taking way to long or failing due to memory issues in my local computer and a remote server.

Do you guys know any other options or suggestions to speed up this process?

Follow the commands I use:

bedtools intersect -a <myvcf>.vcf.gz -b <myinterval>.bed -wa | \
    java -Xmx10g -jar snpSift.jar filter --set <myrsid>.txt "ID in SET[0]"

or

gzcat <myvcf>.vcf.gz | \
    java -Xmx10g -jar snpSift.jar intervals <mybed>.bed | \
    java -Xmx10g -jar snpSift.jar filter --set <myrsid>.txt "ID in SET[0]"

The bedtools command usually fail due to unknown reason and SnpSift runs over 6 hours even given 10GB of ram. My local machine have 8GB of RAM, but the server have 32GB.

I forgot to ask about versions. Maybe a newer version is faster ? What are those "unknown reasons" is there any error or message or just a freezing computer? (Maybe there's a bug). In the second command can't you combine the two snpSift commands in one? How big is the region you want? — llrs, Sep 21 '17 at 14:16
@Llopis I'm using the latest versions of both programs. The bedtools error is usually unrecognized format and the command fails. As far as I know, I can't combine the commands but I never had problem with this scheme. And the region was approximately 32-40Mb — andremrsantos, Sep 21 '17 at 15:56

score 4 · Answer 1 · answered Sep 21 '17 at 15:20

4

this:

bedtools intersect -a <myvcf>.vcf.gz -b <myinterval>.bed -wa | \
    java -Xmx10g -jar snpSift.jar filter --set <myrsid>.txt "ID in SET[0]"

can be replaced with one GATK SelectVariants https://software.broadinstitute.org/gatk/documentation/tooldocs/current/org_broadinstitute_gatk_tools_walkers_variantutils_SelectVariants.php

java -jar gatk.jar -T SelectVariants -R ref.fa \
  -L myinterval.bed --keepIDs myrsid.txt -V myvcf.vcf.gz -o out.vcf

and furthermore, you can also work in parallel by splitting your bed file (e.g. by chromosomes), and merge the results.

answered Sep 21 '17 at 15:20

Pierre

1,536
7
11

thanks for the suggestion, but on my experience GATK usually takes longer to run this kind of process. I've also been able to complete the filtering by inverting the snpSift selection order. – andremrsantos Sep 21 '17 at 15:58
2

@andre if you could post the solution it might help other people – llrs Sep 21 '17 at 17:17
@Llopis sure, please find it bellow – andremrsantos Sep 22 '17 at 11:02

score 4 · Answer 2 · answered Sep 21 '17 at 17:52

Generically, with BEDOPS:

$ vcf2bed < <(gunzip -c snps.vcf) | bedops -e 1 - myRegions.bed > answer.bed

Or:

$ vcf2bed < <(gunzip -c snps.vcf) | bedmap --echo --echo-map-id --delim '\t' myRegions.bed - > answer.bed

Etc.

BEDOPS supports input streams; use them where you can.

That said, if you're going to query the VCF data frequently, convert once and use it in downstream operations.

$ vcf2bed < <(gunzip -c snps.vcf) > snps.bed

Then:

$ bedops -e 1 snps.bed myRegions.bed > answer.bed

Etc.

If you have a computational cluster, you can further parallelize work by splitting SNPs on chromosomes with bedextract:

$ for chr in `bedextract --list-chr snps.bed`; do echo ${chr}; bedextract ${chr} snps.bed > snps.${chr}.bed; done

Then submit a batch job for each chromosome:

$ bedmap --chrom ${chr} --echo --echo-map-id --delim '\t' myRegions.bed snps.${chr}.bed > answer.${chr}.bed

Using bedextract and adding --chrom ${chr} to your bedmap call will yield large speed improvements over other toolkits, as these options focus work on the chromosome of interest and nothing else.

score 3 · Answer 3 · answered Sep 22 '17 at 11:01

I was tinkering with the command and was able to complete the execution inverting the order of rsid and bed intervals filtering. The command is as follows:

gzcat <myvcf>.vcf.gz | \
    java -Xmx10g -jar snpSift.jar filter --set <myrsid>.txt "ID in SET[0]" | \
    java -Xmx10g -jar snpSift.jar intervals <mybed>.bed

Maybe my rsid subset is smaller than my BED intervals and thus the process is quicker and more efficient.

Anyway thanks everyone

score 1 · Answer 4 · edited Mar 01 '19 at 15:34

Here is another way to do it in 2019:

$ bcftools view -T regions.bed -i 'ID=@<myrsid>.txt' input.vcf > output.vcf

Things can be sped up if we compress and index the input.vcf. Then bcftools will have near-immediate random access to the position.

$ bgzip -c input.vcf > input.vcf.gz
$ tabix input.vcf.gz
$ bcftools view -T regions.bed -i 'ID=@<myrsid>.txt' input.vcf.gz > output.vcf

Internally bcftools works with bcf. Whenever the input and output is vcf it converts vcf>bcf>vcf. When dealing with large files, this produces a large overhead. Converting vcf to bcf first can reduce runtime for further downstream analyses dramatically. In my experience 2-3x faster.

$ bcftools view -Ob input.vcf > input.bcf
$ bcftools index input.bcf
$ bcftools view -Ob -T regions.bed -i 'ID=@<myrsid>.txt' input.bcf > output.bcf

Helpful comment about VCF > BCF > VCF conversion! – Daniel Standage Mar 01 '19 at 15:35 — Daniel Standage, Mar 01 '19 at 15:35

Any fast options to query large VCF bed intervals?

4 Answers4