How Can I Get The Field Of Ac Key-Value Pair In The Info Column In Vcf File?
2
0
Entering edit mode
7.2 years ago
cli2 • 0

Dear all,

I am new to DNA sequence analysis, especially variant calling. I had some sequence data and used bwa to perform the alignment. After I had sam file from bwa, I use samtools and bcftools to perform variant calling. However, when I more the final vcf file, I don't see the AC key-value pairs in the vcf file. I cannot calculate the frequency of alleles at each locus. Would anyone have an idea about what I should do?

My purpose is to find the allele frequencies at each locus. If there are other methods to do this, please give me a help. Thank you.

In fact, when I view the sorted bam file, I can see the alignment and the reference. It looks fine.


I tried the following commands in my workflow:

samtools mpileup -g -f ref.fa short_reads_aligned.sorted.bam > short_reads_variants.bcf

bcftools view -c -v short_reads_variants.bcf > short_reads_variants.vcf


I also tried the following

samtools mpileup -uf ref.fa short_reads_aligned.sorted.bam | bcftools view -bvcg - > var.raw.bcf bcftools view var.raw.bcf > var.flt.vcf

Thanks in advance C

variant calling allele frequency vcf • 2.9k views
ADD COMMENT
0
Entering edit mode

Dear pd3,

I don't know how to use this "vcf-annotate" command. Could you please provide more details? Thank you.

ADD REPLY
1
Entering edit mode
7.2 years ago
pd3 ▴ 340

One of many ways to do it is

vcf-annotate --fill-AC-AN

from http://vcftools.sourceforge.net/perl_module.html#vcf-annotate.

Newer versions of samtools+bcftools (develop branch on github) do output AC, AN natively.

ADD COMMENT
0
Entering edit mode
7.1 years ago
cli2 • 0

I finally use GATK and obtain the AC field in my vcf file. I think the problem is due to that the mpileup does not output the allele counts for each locus in the sequence although mpileup used allele counts in the process of generating bcf and vcf files. So the vcf-annotation can add AC field in the vcf but it is not the correct one.

ADD COMMENT

Login before adding your answer.

Traffic: 2133 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