Question: How Can I Get The Field Of Ac Key-Value Pair In The Info Column In Vcf File?
gravatar for cli2
6.5 years ago by
cli20 wrote:

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

ADD COMMENTlink modified 6.4 years ago • written 6.5 years ago by cli20

Dear pd3,

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

ADD REPLYlink written 6.4 years ago by cli20
gravatar for pd3
6.5 years ago by
pd3350 wrote:

One of many ways to do it is

vcf-annotate --fill-AC-AN


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

ADD COMMENTlink modified 10 months ago by RamRS28k • written 6.5 years ago by pd3350
gravatar for cli2
6.4 years ago by
cli20 wrote:

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 COMMENTlink written 6.4 years ago by cli20
Please log in to add an answer.


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