How to use bcftools to calculate AF INFO field from AC and AN in VCF?
1
6
Entering edit mode
5.5 years ago

Hi,

I need to add into my VCFs the AF INFO field, which is not output when variant calling with samtools/bcftools.

I think this can be done with the bcftools annotate command and the AF tag is mentioned in the Expressions part of the bcftools man page. (http://www.htslib.org/doc/bcftools.html#expressions)

But I can't determine how to apply this to my vcf.

For example i've tried:

bcftools annotate -a MyVCF.vcf.gz -c FORMAT,INFO,AF -o test.vcf -O v MyVCF.vcf.gz -h annots.hdr

where annots.hdr is a text file with:

##INFO=<ID=AF,Number=1,Type=Integer,Descripion="Allele frequency">

The header is added in but the calculation on AC/AN is not conducted.

Any suggestions?

Thanks.

samtools bcftools VCF AF • 14k views
ADD COMMENT
0
Entering edit mode

Shortcut just to filter on AF:

bcftools view -q 0.01:minor test.vcf

from here

Works as long as you have the fields INFO/AC and INFO/AN already

ADD REPLY
1
Entering edit mode

This is not a relevant answer to the top level question. This command filters on AF when AF information is already available in the VCF, whereas the question is on calculating AF from AC and AN. I'm moving this post to a comment. Please be careful and read the question carefully before adding answers - especially on older questions.

ADD REPLY
0
Entering edit mode

(can you delete it here)

ADD REPLY
1
Entering edit mode

Please do not add an answer unless you're answering the top level question. If you have a question for Kevin, reply to him by adding a comment on his answer, not by adding an answer yourself.

ADD REPLY
7
Entering edit mode
3.4 years ago

The solution is posted on GitHub: https://github.com/samtools/bcftools/issues/345

1, get the latest htslib and BCFtools, and set-up the BCFtools plugins:

git clone git://github.com/samtools/htslib.git
git clone git://github.com/samtools/bcftools.git
cd bcftools; make
sudo make install

Note the directory where the plugins were installed: /usr/local/libexec/bcftools

export BCFTOOLS_PLUGINS=/usr/local/libexec/bcftools

--------------------------------------

NB - if you don't have root / super-user privileges, then you can avoid the make install command and set the BCFTOOLS_PLUGINS environmental variable to the 'plugins' directory, which will be located under bcftools/ where you downloaded the program. For example export BCFTOOLS_PLUGINS=/Programs/bcftools/plugins/

---------------------------------

2, View the VCF without AF

bcftools view test.vcf | tail -5
[W::bcf_hdr_check_sanity] PL should be declared as Number=G
5   135337248   .   CT  C   .   PASS    END=135337249;HOMLEN=3;HOMSEQ=TTT;SVLEN=-1;SVTYPE=DEL;AC=1;AN=2 GT:AD   0/1:205,118
5   135337259   .   AG  A   .   PASS    END=135337260;HOMLEN=4;HOMSEQ=GGGG;SVLEN=-1;SVTYPE=DEL;AC=1;AN=2    GT:AD   0/1:190,220
5   135337259   .   A   AG  .   PASS    END=135337259;HOMLEN=5;HOMSEQ=GGGGG;SVLEN=1;SVTYPE=INS;AC=1;AN=2    GT:AD   0/1:192,71
5   135337264   .   GA  G   .   PASS    END=135337265;HOMLEN=1;HOMSEQ=A;SVLEN=-1;SVTYPE=DEL;AC=1;AN=2   GT:AD   0/1:184,83
5   135337274   .   A   ATTATTGCATCAACTCCTCCGACATCTCTTCCCCTGCAAGAGTTCAGGCCCACAGGTTCTGGTGTGGGCTTGCTCAGCTGGAGGTAGCCTGAGGTGAGCTGGAG    .PASS   END=135337274;HOMLEN=23;HOMSEQ=TTATTGCATCAACTCCTCCGACA;SVLEN=103;SVTYPE=INS;AC=1;AN=2   GT:AD   0/1:130,17

3, Now add the AF

bcftools +fill-tags test.vcf  -- -t AF | tail -5
[W::bcf_hdr_check_sanity] PL should be declared as Number=G
5   135337248   .   CT  C   .   PASS    END=135337249;HOMLEN=3;HOMSEQ=TTT;SVLEN=-1;SVTYPE=DEL;AC=1;AN=2;AF=0.5  GT:AD   0/1:205,118
5   135337259   .   AG  A   .   PASS    END=135337260;HOMLEN=4;HOMSEQ=GGGG;SVLEN=-1;SVTYPE=DEL;AC=1;AN=2;AF=0.5 GT:AD   0/1:190,220
5   135337259   .   A   AG  .   PASS    END=135337259;HOMLEN=5;HOMSEQ=GGGGG;SVLEN=1;SVTYPE=INS;AC=1;AN=2;AF=0.5 GT:AD   0/1:192,71
5   135337264   .   GA  G   .   PASS    END=135337265;HOMLEN=1;HOMSEQ=A;SVLEN=-1;SVTYPE=DEL;AC=1;AN=2;AF=0.5    GT:AD   0/1:184,83
5   135337274   .   A   ATTATTGCATCAACTCCTCCGACATCTCTTCCCCTGCAAGAGTTCAGGCCCACAGGTTCTGGTGTGGGCTTGCTCAGCTGGAGGTAGCCTGAGGTGAGCTGGAG    .PASS   END=135337274;HOMLEN=23;HOMSEQ=TTATTGCATCAACTCCTCCGACA;SVLEN=103;SVTYPE=INS;AC=1;AN=2;AF=0.5    GT:AD   0/1:130,17

4, add all available tags in plugin

bcftools +fill-tags test.vcf | tail -5
[W::bcf_hdr_check_sanity] PL should be declared as Number=G
5   135337248   .   CT  C   .   PASS    END=135337249;HOMLEN=3;HOMSEQ=TTT;SVLEN=-1;SVTYPE=DEL;AC=1;AN=2;NS=1;AF=0.5;MAF=0.5;AC_Het=1;AC_Hom=0;AC_Hemi=0;HWE=1;ExcHet=1  GT:AD   0/1:205,118
5   135337259   .   AG  A   .   PASS    END=135337260;HOMLEN=4;HOMSEQ=GGGG;SVLEN=-1;SVTYPE=DEL;AC=1;AN=2;NS=1;AF=0.5;MAF=0.5;AC_Het=1;AC_Hom=0;AC_Hemi=0;HWE=1;ExcHet=1 GT:AD   0/1:190,220
5   135337259   .   A   AG  .   PASS    END=135337259;HOMLEN=5;HOMSEQ=GGGGG;SVLEN=1;SVTYPE=INS;AC=1;AN=2;NS=1;AF=0.5;MAF=0.5;AC_Het=1;AC_Hom=0;AC_Hemi=0;HWE=1;ExcHet=1 GT:AD   0/1:192,71
5   135337264   .   GA  G   .   PASS    END=135337265;HOMLEN=1;HOMSEQ=A;SVLEN=-1;SVTYPE=DEL;AC=1;AN=2;NS=1;AF=0.5;MAF=0.5;AC_Het=1;AC_Hom=0;AC_Hemi=0;HWE=1;ExcHet=1    GT:AD   0/1:184,83
5   135337274   .   A   ATTATTGCATCAACTCCTCCGACATCTCTTCCCCTGCAAGAGTTCAGGCCCACAGGTTCTGGTGTGGGCTTGCTCAGCTGGAGGTAGCCTGAGGTGAGCTGGAG    .PASS   END=135337274;HOMLEN=23;HOMSEQ=TTATTGCATCAACTCCTCCGACA;SVLEN=103;SVTYPE=INS;AC=1;AN=2;NS=1;AF=0.5;MAF=0.5;AC_Het=1;AC_Hom=0;AC_Hemi=0;HWE=1;ExcHet=1    GT:AD   0/1:130,17

Kevin

ADD COMMENT
0
Entering edit mode

Hi Kevin

This bcf +fill-tags cant fill details in my vcf file. I tried:

 bcftools +fill-tags my.vcf > out.vcf
 bcftools +fill-tags my.vcf -Ov --output fill-all.vcf -- -t AN,AC
 bcftools +fill-tags my.vcf -Ov --output fill-all.vcf -- -t all

It only fills NS: output is below:

DP=1;ADF=0;ADR=0;AD=0;SCR=0;MQ0F=0;AN=0;DP4=0,0,0,0;MQ=.;NS=0

my.vcf is output of RNAseq :

bcftools mpileup my.bam -f $REF --annotate FORMAT/AD,FORMAT/ADF,FORMAT/ADR,FORMAT/DP,FORMAT/SP,FORMAT/SCR,INFO/AD,INFO/ADF,INFO/ADR,INFO/SCR | bcftools call -m -O v --output my.vcf

Any idea? Thanks in advance.

ADD REPLY
0
Entering edit mode

Can you paste a few full line records from the starting VCF?

ADD REPLY
0
Entering edit mode

Here is the file: https://drive.google.com/file/d/1nDwIRlT9Qwfec7yZp8Q7Gpn_LnJGBTuH/view?usp=sharing

over 200 lines might look like a mess here :)

ADD REPLY
1
Entering edit mode

You can use GitHub gist to post the file. It would show up as a preview yet not be actually stored on the site.

ADD REPLY
0
Entering edit mode

Thanks. The problem seems to be that your VCF is like en empty 'shell' - all variants are like this:

#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  Y.bam
Y   5325732 .   T   .   29.5864 .   DP=1;ADF=0;ADR=0;AD=0;SCR=0;MQ0F=0;AN=0;DP4=0,0,0,0;MQ=.;NS=0   GT:DP:SP:ADF:ADR:AD:SCR ./.:0:0:0:0:0:0
Y   5325733 .   T   .   29.5864 .   DP=1;ADF=0;ADR=0;AD=0;SCR=0;MQ0F=0;AN=0;DP4=0,0,0,0;MQ=.;NS=0   GT:DP:SP:ADF:ADR:AD:SCR ./.:0:0:0:0:0:0
Y   5325734 .   A   .   29.5864 .   DP=1;ADF=0;ADR=0;AD=0;SCR=0;MQ0F=0;AN=0;DP4=0,0,0,0;MQ=.;NS=0   GT:DP:SP:ADF:ADR:AD:SCR ./.:0:0:0:0:0:0
Y   5325735 .   T   .   29.5864 .   DP=1;ADF=0;ADR=0;AD=0;SCR=0;MQ0F=0;AN=0;DP4=0,0,0,0;MQ=.;NS=0   GT:DP:SP:ADF:ADR:AD:SCR ./.:0:0:0:0:0:0
Y   5325736 .   G   .   29.5864 .   DP=1;ADF=0;ADR=0;AD=0;SCR=0;MQ0F=0;AN=0;DP4=0,0,0,0;MQ=.;NS=0   GT:DP:SP:ADF:ADR:AD:SCR ./.:0:0:0:0:0:0
Y   5325737 .   T   .   29.5864 .   DP=1;ADF=0;ADR=0;AD=0;SCR=0;MQ0F=0;AN=0;DP4=0,0,0,0;MQ=.;NS=0   GT:DP:SP:ADF:ADR:AD:SCR ./.:0:0:0:0:0:0
Y   5325738 .   G   .   29.5864 .   DP=1;ADF=0;ADR=0;AD=0;SCR=0;MQ0F=0;AN=0;DP4=0,0,0,0;MQ=.;NS=0

In fact, the entire file just contains reference calls.

So, the BCFtools plugin cannot function correctly.

How did you produce my.bam ?

ADD REPLY
1
Entering edit mode

You are right, bcftools +fill-tags doesnt add anything if there is only ref calls. I run it on different file and got all the tags added. Thank you for correcting me 🙏😀.

ADD REPLY

Login before adding your answer.

Traffic: 2220 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6