3
$\begingroup$

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  ./.
Improve this question
$\endgroup$
2
  • $\begingroup$ Can you give us an example of your file? $\endgroup$
    – terdon ♦
    Commented Nov 7, 2022 at 14:39
  • $\begingroup$ Done! See edited post. $\endgroup$
    – mugdi
    Commented Nov 7, 2022 at 14:47

1 Answer 1

Reset to default
3
$\begingroup$

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
Improve this answer
$\endgroup$
2
  • $\begingroup$ 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 ? $\endgroup$
    – mugdi
    Commented Nov 7, 2022 at 15:50
  • 1
    $\begingroup$ 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. $\endgroup$
    – terdon ♦
    Commented Nov 7, 2022 at 16:08

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.