Answers.Body: <p>Generically, with <a href="http://bedops.readthedocs.io/en/latest/" rel="nofollow noreferrer">BEDOPS</a>:</p>
  
  <pre><code>$ vcf2bed &lt; &lt;(gunzip -c snps.vcf) | bedops -e 1 - myRegions.bed &gt; answer.bed
  </code></pre>
  
  <p>Or:</p>
  
  <pre><code>$ vcf2bed &lt; &lt;(gunzip -c snps.vcf) | bedmap --echo --echo-map-id --delim '\t' myRegions.bed - &gt; answer.bed
  </code></pre>
  
  <p>Etc.</p>
  
  <p>BEDOPS supports input streams; use them where you can.</p>
  
  <p>That said, if you're going to query the VCF data frequently, convert once and use it in downstream operations.</p>
  
  <pre><code>$ vcf2bed &lt; &lt;(gunzip -c snps.vcf) &gt; snps.bed
  </code></pre>
  
  <p>Then:</p>
  
  <pre><code>$ bedops -e 1 snps.bed myRegions.bed &gt; answer.bed
  </code></pre>
  
  <p>Etc.</p>
  
  <p>If you have a computational cluster, you can further parallelize work by splitting SNPs on chromosomes with <code>bedextract</code>:</p>
  
  <pre><code>$ for chr in `bedextract --list-chr snps.bed`; do echo ${chr}; bedextract ${chr} snps.bed &gt; snps.${chr}.bed; done
  </code></pre>
  
  <p>Then submit a batch job for each chromosome:</p>
  
  <pre><code>$ bedmap --chrom ${chr} --echo --echo-map-id --delim '\t' myRegions.bed snps.${chr}.bed &gt; answer.${chr}.bed
  </code></pre>
  
  <p>Using <code>bedextract</code> and adding <code>--chrom ${chr}</code> to your <code>bedmap</code> call will yield large speed improvements over other toolkits, as these options focus work on the chromosome of interest and nothing else.</p>
  
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-09-21T17:52:45.343Z
Answers.LastActivityDate: 2017-09-21T17:52:45.343Z
Answers.UserDisplayName: Alex Reynolds
Answers.Body: <p>this:</p>
  
  <pre><code>bedtools intersect -a &lt;myvcf&gt;.vcf.gz -b &lt;myinterval&gt;.bed -wa | \
      java -Xmx10g -jar snpSift.jar filter --set &lt;myrsid&gt;.txt "ID in SET[0]"
  </code></pre>
  
  <p>can be replaced with one GATK SelectVariants <a href="https://software.broadinstitute.org/gatk/documentation/tooldocs/current/org_broadinstitute_gatk_tools_walkers_variantutils_SelectVariants.php" rel="nofollow noreferrer">https://software.broadinstitute.org/gatk/documentation/tooldocs/current/org_broadinstitute_gatk_tools_walkers_variantutils_SelectVariants.php</a></p>
  
  <pre><code>java -jar gatk.jar -T SelectVariants -R ref.fa \
    -L myinterval.bed --keepIDs myrsid.txt -V myvcf.vcf.gz -o out.vcf
  </code></pre>
  
  <p>and furthermore, you can also work in parallel by splitting your bed file (e.g. by chromosomes), and merge the results.</p>
  
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-09-21T15:58:33.687Z
Answers.Text: 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.
Answers.UserDisplayName: andremrsantos
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-09-21T17:17:46.623Z
Answers.Text: @andre if you could post the solution it might help other people
Answers.UserDisplayName: llrs
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-09-22T11:02:27.980Z
Answers.Text: @Llopis sure, please find it bellow
Answers.UserDisplayName: andremrsantos
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-09-21T15:20:23.450Z
Answers.LastActivityDate: 2017-09-21T15:20:23.450Z
Answers.UserDisplayName: Pierre
Answers.Body: <p>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:</p>
  
  <pre><code>gzcat &lt;myvcf&gt;.vcf.gz | \
      java -Xmx10g -jar snpSift.jar filter --set &lt;myrsid&gt;.txt "ID in SET[0]" | \
      java -Xmx10g -jar snpSift.jar intervals &lt;mybed&gt;.bed
  </code></pre>
  
  <p>Maybe my rsid subset is smaller than my BED intervals and thus the process is quicker and more efficient.</p>
  
  <p>Anyway thanks everyone</p>
  
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-09-22T11:01:52.813Z
Answers.LastActivityDate: 2017-09-22T11:01:52.813Z
Answers.UserDisplayName: andremrsantos
Answers.Body: <p>Here is another way to do it in 2019:</p>
  
  <pre><code>$ bcftools view -T regions.bed -i 'ID=@&lt;myrsid&gt;.txt' input.vcf &gt; output.vcf
  </code></pre>
  
  <p>Things can be sped up if we compress and index the <code>input.vcf</code>. Then <code>bcftools</code> will have near-immediate random access to the position.</p>
  
  <pre><code><span class="math-container">$ bgzip -c input.vcf &gt; input.vcf.gz
  $</span> tabix input.vcf.gz
  $ bcftools view -T regions.bed -i 'ID=@&lt;myrsid&gt;.txt' input.vcf.gz &gt; output.vcf
  </code></pre>
  
  <p>Internally <code>bcftools</code> works with <code>bcf</code>. Whenever the input and output is <code>vcf</code> it converts <code>vcf</code>><code>bcf</code>><code>vcf</code>. When dealing with large files, this produces a large overhead. Converting <code>vcf</code> to <code>bcf</code> first can reduce runtime for further downstream analyses dramatically. In my experience 2-3x faster.</p>
  
  <pre><code><span class="math-container">$ bcftools view -Ob input.vcf &gt; input.bcf
  $</span> bcftools index input.bcf
  $ bcftools view -Ob -T regions.bed -i 'ID=@&lt;myrsid&gt;.txt' input.bcf &gt; output.bcf
  </code></pre>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2019-03-01T15:35:15.873Z
Answers.Text: Helpful comment about VCF > BCF > VCF conversion!
Answers.UserDisplayName: Daniel Standage
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2019-01-31T05:04:07.483Z
Answers.LastActivityDate: 2019-03-01T15:34:47.657Z
Answers.LastEditDate: 2019-03-01T15:34:47.657Z
Answers.UserDisplayName: finswimmer
Body: <p>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. </p>
  
  <p>Unfortunately, both my normal choices to do this analysis (<code>snpSift</code> and <code>bedtools</code>) are taking way to long or failing due to memory issues in my local computer and a remote server.</p>
  
  <p>Do you guys know any other options or suggestions to speed up this process?</p>
  
  <p>Follow the commands I use:</p>
  
  <pre><code>bedtools intersect -a &lt;myvcf&gt;.vcf.gz -b &lt;myinterval&gt;.bed -wa | \
      java -Xmx10g -jar snpSift.jar filter --set &lt;myrsid&gt;.txt "ID in SET[0]"
  </code></pre>
  
  <p>or</p>
  
  <pre><code>gzcat &lt;myvcf&gt;.vcf.gz | \
      java -Xmx10g -jar snpSift.jar intervals &lt;mybed&gt;.bed | \
      java -Xmx10g -jar snpSift.jar filter --set &lt;myrsid&gt;.txt "ID in SET[0]"
  </code></pre>
  
  <p>The <code>bedtools</code> command usually fail due to unknown reason and <code>SnpSift</code> runs over 6 hours even given 10GB of ram. My local machine have 8GB of RAM, but the server have 32GB.</p>
  
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-09-21T14:16:00.627Z
Comments.Text: 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?
Comments.UserDisplayName: llrs
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-09-21T15:56:50.970Z
Comments.Text: @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
Comments.UserDisplayName: andremrsantos
Tags: <vcf><bed><filtering><interval>
Title: Any fast options to query large VCF bed intervals?
ThreadText: # Any fast options to query large VCF bed intervals?
  
  (Asked by: username_0 on 2017-09-21)
  
  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.
  
  **Tags:** vcf, bed, filtering, interval
  
  > **Comments:**
  > **username_1**: 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?  
  > **username_0**: @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  
  
  ---
  
  ## 4 answers
  
  ### Answer by username_2 on 2017-09-21. Score: 4
  
  Generically, with [BEDOPS](http://bedops.readthedocs.io/en/latest/):
  
  ```
  $ 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.
  
  
  ---
  
  ### Answer by username_3 on 2017-09-21. Score: 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.
  
  > **Comments:**
  > **username_0**: 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.  
  > **username_1**: @andre if you could post the solution it might help other people  
  > **username_0**: @Llopis sure, please find it bellow  
  
  ---
  
  ### Answer by username_0 on 2017-09-22. Score: 3
  
  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
  
  
  ---
  
  ### Answer by username_4 on 2019-01-31. Score: 1
  
  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
  
  ```
  
  > **Comments:**
  > **username_5**: Helpful comment about VCF > BCF > VCF conversion!  
  
  ---
  
  