<p>First of all, bcftools might be a bit of overkill if you want one specific variant. All you need is GNU <code>grep</code> (the default on Linux):</p>
<pre><code>zgrep -P '^(#|1\t630128\t)' file.vcf.gz 
</code></pre>
<p>Or, with any <code>grep</code>:</p>
<pre><code>zgrep -E '^(#|1'$'\t''630128'$'\t'')' file.vcf.gz 
</code></pre>
<p>That will print all lines that start with a <code>#</code> (the headers) or that start with <code>1</code>, followed by a tab, then <code>630128</code> and another tab.</p>
<p>Alternatively, you can use <code>awk</code>:</p>
<pre><code>zcat file.vcf.gz |
    awk -F'\t' '/^#/ || ( $1==1 &amp;&amp; $2==630128)'
</code></pre>
<hr />
<p>All that said, yes, you can also do it with <code>bcftools</code> and in most cases that is the right tool for the job so you might as well use it here too to get used to it. However, you need the right command and that's <code>view</code>, not <code>filter</code>:</p>
<pre><code>bcftools view -r 1:630128 /home/oem/Documents/GitHub/Genomics_prac_guide/trio/files/genome_file.vcf.gz 
</code></pre>
<p>Note that the file needs to be indexed, so if you have not done so already, run</p>
<pre><code>tabix file.vcf.gz
</code></pre>
<hr />
<p>Also, the reason your unquoted command failed as it did is because the <code>&amp;</code> has special meaning to the shell, it means &quot;run this command in the background&quot;. Therefore, when you ran <code>bcftools filter -e CHROM=1&amp;POS=63018 file.vcf.gz</code> that actually ran two separate commands:</p>
<pre><code>bcftools filter -e CHROM=1&amp;
</code></pre>
<p>And then:</p>
<pre><code>POS=63018 file.vcf.gz
</code></pre>
<p>Which means &quot;set the variable <code>POS</code> to the value <code>63018</code> and execute the command <code>file.vcf.gz</code>. Since the file isn't executable, you got:</p>
<pre><code>sh: 1: file.vcf.gz: Permission denied
</code></pre>

CC BY-SA 4.0
2022-11-07T15:50:59.613Z
First of all thanks for the quite detailed answer! This gives me a good introduction into bcftools environment. So when I use view to filter SNP's the reason for  `-filter` is then to apply some standard filter to the file? 
Second question: If I want to remove this single variant from the file is the correct way to do this `bcftools view -r 1:630128 -e file.vcf.gz` ?
mugdi
CC BY-SA 4.0
2022-11-07T16:08:35.750Z
I think `filter` is more to get things like "show me variants where the AB is X and the GT is Y" and things like this, while in this case you are just focusing on a specific region, albeit a region of length 1, so `view` is the command for this. Reversing the match is a completely different question. I don't know if bcftools can do that, but `zcat file.vcf.gz | grep -vP '^1\t630128\t'` will get rid of all variants not matching that specific position.
terdon
CC BY-SA 4.0
2022-11-07T14:56:26.803Z
2022-11-07T14:56:26.803Z
terdon
<p>I want to filter a SNP, specifically <code>CHROM1:POS:630128</code>. I tried to use <code>bcftools</code> for this.</p>
<p>I came up with the following variant of a bcftools call:</p>
<pre><code>bcftools filter -e &quot;CHROM=1&amp;POS=63018&quot; /home/oem/Documents/GitHub/Genomics_prac_guide/trio/files/genome_file.vcf.gz 
</code></pre>
<p>which leads to:</p>
<pre><code>Error: cannot use arithmetic operators to compare strings and numbers
</code></pre>
<p>and the call with unquoted parameters leads to</p>
<pre><code>$ bcftools filter -e CHROM=1&amp;POS=63018 /home/oem/Documents/GitHub/Genomics_prac_guide/trio/files/genome_file.vcf.gz 

sh: 1: /home/oem/Documents/GitHub/Genomics_prac_guide/trio/file/genome_file.vcf.gz: Permission denied
Failed to read from standard input: unknown file type
</code></pre>
<ul>
<li>Q1: Is bfctools filter the correct tool here to use?</li>
<li>Q2: Regarding <a href="https://github.com/samtools/bcftools/issues/1028#issuecomment-497651048" rel="nofollow noreferrer">this comment</a> the first error might indicate that I have to annotate the CHROM &amp; POS Tags in the header file. Is this correct?</li>
</ul>
<p>An example of the file:</p>
<pre><code>##fileformat=VCFv4.2
##fileDate=20221104
##source=PLINKv1.90
##contig=&lt;ID=1,length=248928329&gt;
##contig=&lt;ID=2,length=242157857&gt;
##contig=&lt;ID=3,length=198157342&gt;
##contig=&lt;ID=4,length=190106194&gt;
##contig=&lt;ID=5,length=181288140&gt;
##contig=&lt;ID=6,length=170595721&gt;
##contig=&lt;ID=7,length=159331484&gt;
##contig=&lt;ID=8,length=145075470&gt;
##contig=&lt;ID=9,length=138211490&gt;
##contig=&lt;ID=10,length=133659511&gt;
##contig=&lt;ID=11,length=135075227&gt;
##contig=&lt;ID=12,length=133261768&gt;
##contig=&lt;ID=13,length=114341522&gt;
##contig=&lt;ID=14,length=106873816&gt;
##contig=&lt;ID=15,length=101922277&gt;
##contig=&lt;ID=16,length=90172890&gt;
##contig=&lt;ID=17,length=83203771&gt;
##contig=&lt;ID=18,length=80252738&gt;
##contig=&lt;ID=19,length=58586567&gt;
##contig=&lt;ID=20,length=64323573&gt;
##contig=&lt;ID=21,length=46679699&gt;
##contig=&lt;ID=22,length=50776369&gt;
##INFO=&lt;ID=PR,Number=0,Type=Flag,Description=&quot;Provisional reference allele, may not be based on real reference genome&quot;&gt;
##FORMAT=&lt;ID=GT,Number=1,Type=String,Description=&quot;Genotype&quot;&gt;
#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  FAM001_ID001
1   69869   rs548049170 T   .   .   .   PR  GT  0/0
1   74792   rs13328684  N   .   .   .   PR  GT  ./.
1   630128  rs9283150   A   .   .   .   PR  GT  0/0
1   791532  i713426 N   .   .   .   PR  GT  ./.
1   792461  rs116587930 G   .   .   .   PR  GT  0/0
1   817341  rs3131972   G   .   .   .   PR  GT  0/0
1   818725  rs12184325  C   .   .   .   PR  GT  0/0
1   820888  rs12567639  A   .   .   .   PR  GT  0/0
1   823656  rs114525117 G   .   .   .   PR  GT  0/0
1   841166  rs12124819  G   A   .   .   PR  GT  0/1
1   858952  rs12127425  G   .   .   .   PR  GT  0/0
1   866156  rs79373928  T   .   .   .   PR  GT  0/0
1   880041  rs72888853  N   .   .   .   PR  GT  ./.
</code></pre>

CC BY-SA 4.0
2022-11-07T14:39:17.783Z
Can you give us an example of your file?
terdon
CC BY-SA 4.0
2022-11-07T14:47:19.367Z
Done! See edited post.
mugdi
<vcf><filtering><bcftools>
Filter with bcftools
# Filter with bcftools

(Asked by: username_0 on 2022-11-07)

I want to filter a SNP, specifically CHROM1:POS:630128. I tried to use bcftools for this.
I came up with the following variant of a bcftools call:
bcftools filter -e "CHROM=1&POS=63018" /home/oem/Documents/GitHub/Genomics_prac_guide/trio/files/genome_file.vcf.gz 

which leads to:
Error: cannot use arithmetic operators to compare strings and numbers

and the call with unquoted parameters leads to
$ bcftools filter -e CHROM=1&POS=63018 /home/oem/Documents/GitHub/Genomics_prac_guide/trio/files/genome_file.vcf.gz 

sh: 1: /home/oem/Documents/GitHub/Genomics_prac_guide/trio/file/genome_file.vcf.gz: Permission denied
Failed to read from standard input: unknown file type


Q1: Is bfctools filter the correct tool here to use?
Q2: Regarding this comment the first error might indicate that I have to annotate the CHROM & POS Tags in the header file. Is this correct?

An example of the file:
##fileformat=VCFv4.2
##fileDate=20221104
##source=PLINKv1.90
##contig=
##contig=
##contig=
##contig=
##contig=
##contig=
##contig=
##contig=
##contig=
##contig=
##contig=
##contig=
##contig=
##contig=
##contig=
##contig=
##contig=
##contig=
##contig=
##contig=
##contig=
##contig=
##INFO=
##FORMAT=
#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  FAM001_ID001
1   69869   rs548049170 T   .   .   .   PR  GT  0/0
1   74792   rs13328684  N   .   .   .   PR  GT  ./.
1   630128  rs9283150   A   .   .   .   PR  GT  0/0
1   791532  i713426 N   .   .   .   PR  GT  ./.
1   792461  rs116587930 G   .   .   .   PR  GT  0/0
1   817341  rs3131972   G   .   .   .   PR  GT  0/0
1   818725  rs12184325  C   .   .   .   PR  GT  0/0
1   820888  rs12567639  A   .   .   .   PR  GT  0/0
1   823656  rs114525117 G   .   .   .   PR  GT  0/0
1   841166  rs12124819  G   A   .   .   PR  GT  0/1
1   858952  rs12127425  G   .   .   .   PR  GT  0/0
1   866156  rs79373928  T   .   .   .   PR  GT  0/0
1   880041  rs72888853  N   .   .   .   PR  GT  ./.

**Tags:** vcf, filtering, bcftools
> **Comments:**
> **username_1**: Can you give us an example of your file?  
> **username_0**: Done! See edited post.  

---

## 1 answers

### Answer by username_1 on 2022-11-07. Score: 3 | Accepted answer

First of all, bcftools might be a bit of overkill if you want one specific variant. All you need is GNU grep (the default on Linux):
zgrep -P '^(#|1\t630128\t)' file.vcf.gz 

Or, with any grep:
zgrep -E '^(#|1'$'\t''630128'$'\t'')' file.vcf.gz 

That will print all lines that start with a # (the headers) or that start with 1, followed by a tab, then 630128 and another tab.
Alternatively, you can use awk:
zcat file.vcf.gz |
    awk -F'\t' '/^#/ || ( $1==1 && $2==630128)'


All that said, yes, you can also do it with bcftools and in most cases that is the right tool for the job so you might as well use it here too to get used to it. However, you need the right command and that's view, not filter:
bcftools view -r 1:630128 /home/oem/Documents/GitHub/Genomics_prac_guide/trio/files/genome_file.vcf.gz 

Note that the file needs to be indexed, so if you have not done so already, run
tabix file.vcf.gz


Also, the reason your unquoted command failed as it did is because the & has special meaning to the shell, it means "run this command in the background". Therefore, when you ran bcftools filter -e CHROM=1&POS=63018 file.vcf.gz that actually ran two separate commands:
bcftools filter -e CHROM=1&

And then:
POS=63018 file.vcf.gz

Which means "set the variable POS to the value 63018 and execute the command file.vcf.gz. Since the file isn't executable, you got:
sh: 1: file.vcf.gz: Permission denied

> **Comments:**
> **username_0**: First of all thanks for the quite detailed answer! This gives me a good introduction into bcftools environment. So when I use view to filter SNP's the reason for  `-filter` is then to apply some standard filter to the file?  Second question: If I want to remove this single variant from the file is the correct way to do this `bcftools view -r 1:630128 -e file.vcf.gz` ?  
> **username_1**: I think `filter` is more to get things like "show me variants where the AB is X and the GT is Y" and things like this, while in this case you are just focusing on a specific region, albeit a region of length 1, so `view` is the command for this. Reversing the match is a completely different question. I don't know if bcftools can do that, but `zcat file.vcf.gz | grep -vP '^1\t630128\t'` will get rid of all variants not matching that specific position.  

---

