text: user16454: 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=<ID=1,length=248928329>
  ##contig=<ID=2,length=242157857>
  ##contig=<ID=3,length=198157342>
  ##contig=<ID=4,length=190106194>
  ##contig=<ID=5,length=181288140>
  ##contig=<ID=6,length=170595721>
  ##contig=<ID=7,length=159331484>
  ##contig=<ID=8,length=145075470>
  ##contig=<ID=9,length=138211490>
  ##contig=<ID=10,length=133659511>
  ##contig=<ID=11,length=135075227>
  ##contig=<ID=12,length=133261768>
  ##contig=<ID=13,length=114341522>
  ##contig=<ID=14,length=106873816>
  ##contig=<ID=15,length=101922277>
  ##contig=<ID=16,length=90172890>
  ##contig=<ID=17,length=83203771>
  ##contig=<ID=18,length=80252738>
  ##contig=<ID=19,length=58586567>
  ##contig=<ID=20,length=64323573>
  ##contig=<ID=21,length=46679699>
  ##contig=<ID=22,length=50776369>
  ##INFO=<ID=PR,Number=0,Type=Flag,Description="Provisional reference allele, may not be based on real reference genome">
  ##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
  #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  ./.
  
  
  user298: 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
  
  
original_text: 16454: <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>
   298: <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>
  