Calling SNP frequencies in BWA alignments

[Nepaliaashish]

Namaste all.

I am using BWA to align my reads. I need to not only call the variants but also their frequencies, for example, if I have 30X coverage at a chromosomal position, the reference has A and of the 30 reads I have 20 As and 10 Ts, how do I do that? any tools available for this?

Any help would be highly appreciated.

[Yilong Li]

Have a look at samtools mpileup and GATK UnifiedGenotyper, I think they can do what you want.

[Nepaliaashish]

Thanks YiLong. Both the tools you mentioned can only call allele frequencies if multiple aligned .bam files are present. They do not call allele frequencies from different reads of the same sample. What I need is: if base (say at chr2 position 135487) is heterozygous A/T, say reference is A, and I have 30 reads at that position. I want to know how many reads are As and how many are Ts. I want to this for every position in the genome for every sample. I am not sure whether mpileup or GATK do such a thing.

[Yilong Li]

In mpileup in the INFO field, there's the AF1 subfield that'll tell you the frequency.

In GATK UnifiedGenotyper the same info can be found in the last field in the AD subfield.

At least this is how I have understood it though I haven't read this from anywhere so use with caution... Also you might have to play around with the parameters to get these fields into the output - I don't remember what parameters I used for generating these .vcf files.