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
2 answers
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.
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.
Log in to answer this question.
-
Sam To Bam
written by freddy •
-
Removing Phix Data From Ngs Data For Analysis
written by crivenster
-
Why Do We Need To Set Maximal Read Depth In Samtools
written by a.a.kechin •
-
Large Bam, But Small Mpileup (Samtools Mpileup)
written by gen •
-
Reference Genome Path In Samtools
written by win
-
Samtools Call Snps
written by GPR
-
Clarification For Running Bwa In Parallel
written by Mohamed •
-
Predict/Estimate/Find Bwa Parameters From Bam Or Sam File
written by Houkto •
-
Consensus With Mpileup
written by Dunand •
-
How To Interpret -D File In Samtools Variant Calling Pipeline?
written by Ian
Dear pd3,
I don't know how to use this "vcf-annotate" command. Could you please provide more details? Thank you.