Question: How to use bcftools to calculate AF INFO field from AC and AN in VCF?
4
gravatar for Nicholas Blackburn
4.1 years ago by
Tasmania, Australia
Nicholas Blackburn100 wrote:

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.

af bcftools samtools vcf • 6.9k views
ADD COMMENTlink modified 4 weeks ago by Ks-seq10 • written 4.1 years ago by Nicholas Blackburn100

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 REPLYlink modified 28 days ago by RamRS26k • written 5 months ago by zimazamiz0

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 REPLYlink written 5 months ago by RamRS26k

(can you delete it here)

ADD REPLYlink modified 28 days ago • written 4 weeks ago by Ks-seq10
1

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 REPLYlink written 4 weeks ago by RamRS26k
7
gravatar for Kevin Blighe
23 months ago by
Kevin Blighe56k
Kevin Blighe56k wrote:

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 COMMENTlink modified 11 months ago • written 23 months ago by Kevin Blighe56k

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 REPLYlink modified 28 days ago by Kevin Blighe56k • written 29 days ago by Ks-seq10

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

ADD REPLYlink written 28 days ago by Kevin Blighe56k

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 REPLYlink written 28 days ago by Ks-seq10
1

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 REPLYlink written 28 days ago by RamRS26k

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 REPLYlink written 28 days ago by Kevin Blighe56k
1

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 REPLYlink written 25 days ago by Ks-seq10
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 1693 users visited in the last hour