Question: How to use bcftools to calculate AF INFO field from AC and AN in VCF?
5
gravatar for Nicholas Blackburn
5.0 years ago by
Tasmania, Australia
Nicholas Blackburn110 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 • 11k views
ADD COMMENTlink modified 11 months ago by Ks-seq10 • written 5.0 years ago by Nicholas Blackburn110

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 11 months ago by Ram32k • written 15 months ago by zimazamiz0
1

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 15 months ago by Ram32k

(can you delete it here)

ADD REPLYlink modified 11 months ago • written 11 months 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 11 months ago by Ram32k
7
gravatar for Kevin Blighe
2.8 years ago by
Kevin Blighe70k
Republic of Ireland
Kevin Blighe70k 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 22 months ago • written 2.8 years ago by Kevin Blighe70k

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 11 months ago by Kevin Blighe70k • written 11 months ago by Ks-seq10

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

ADD REPLYlink written 11 months ago by Kevin Blighe70k

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 11 months 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 11 months ago by Ram32k

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 11 months ago by Kevin Blighe70k
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 11 months 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: 977 users visited in the last hour
_