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

Question

How to use bcftools to calculate AF INFO field from AC and AN in VCF?

7

Entering edit mode

8.1 years ago

Nicholas Blackburn ▴ 140

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 • 26k views

ADD COMMENT • link updated 2.0 years ago by ross_whetten ▴ 10 • written 8.1 years ago by Nicholas Blackburn ▴ 140

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 • link updated 4.1 years ago by Ram 43k • written 4.5 years ago by zimazamiz • 0

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 • link 4.5 years ago by Ram 43k

Kevin Blighe · Accepted Answer · 2018-05-07

8

Entering edit mode

6.0 years ago

Kevin Blighe 87k

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 • link 5.0 years ago by Kevin Blighe 87k

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 • link updated 4.1 years ago by Kevin Blighe 87k • written 4.1 years ago by Ks-seq ▴ 10

0

Entering edit mode

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

ADD REPLY • link 4.1 years ago by Kevin Blighe 87k

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 • link 4.1 years ago by Ks-seq ▴ 10

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 • link 4.1 years ago by Ram 43k

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 • link 4.1 years ago by Kevin Blighe 87k

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 • link 4.1 years ago by Ks-seq ▴ 10

0

Entering edit mode

Hi Kevin,

I'm trying this approach to add AF info as I'm running mouse exome sequencing samples through GATK/Mutect. I'm using the mm39 variant information from Ensembl (from here: http://ftp.ensembl.org/pub/release-105/variation/vcf/mus_musculus/). This is an example of what the VCF file reads as:

enter image description here

I have tried ./bcftools plugin fill-tags variant.vcf.gz -- -t AF | tail -5

However the output doesn't seem to contain AF. I've also tried adding All tags but the output remains the same. Any suggestions for what to try?

ADD REPLY • link 2.1 years ago by Vincent • 0

0

Entering edit mode

Please do not paste screenshots of plain text content, it is counterproductive. You can copy paste the content directly here (using the code formatting option shown below), or use a GitHub Gist if the content volume exceeds allowed length here.

code_formatting

ADD REPLY • link 2.1 years ago by Ram 43k

0

Entering edit mode

The syntax for bcftools plugins requires a + symbol rather than the word "plugin" before the plugin name, e.g. bcftools +fill-tags test.vcf -- -t AF

ADD REPLY • link 2.0 years ago by ross_whetten ▴ 10