question about samtools and snp calling

[guang918]

Hi, everyone,

I am trying to call SNPs from a single .bam file by Samtools. From the .vcf output file, I found the allele frequency (AF1) is either 1 or 0.5. I can't figure out whether I did something wrong here. Has any one seen similar situation?

Call SNPs from a single alignment:
$> samtools mpileup –uf myreference.fasta myalignment.bam > tmp_output.bcf
$> bcftools view –bvcg tmp_output.bcf > raw.bcf
$> bcftools view > raw.vcf
$> vcfutils.pl varFilter –D100 raw.vcf > my_samtool_snps.vcf

thanks very much.

best,

shouhong

[vivek_]

Hi, everyone,

I am trying to call SNPs from a single .bam file by Samtools. From the .vcf output file, I found the allele frequency (AF1) is either 1 or 0.5. I can't figure out whether I did something wrong here. Has any one seen similar situation?

Call SNPs from a single alignment:
$> samtools mpileup –uf myreference.fasta myalignment.bam > tmp_output.bcf
$> bcftools view –bvcg tmp_output.bcf > raw.bcf
$> bcftools view > raw.vcf
$> vcfutils.pl varFilter –D100 raw.vcf > my_samtool_snps.vcf

thanks very much.

best,

shouhong

If you are calling SNPs on a single individual, the genotype can either be homozygous or heterozygous, if its homozygous, the alternate allele appears twice out of two possible chances in one sample, so the allele frequency is 1, if it is het, it appears once out of two times and the allele frequency is 0.5.

Allele frequency is not particularly informative unless you have more than one sample being represented in the VCF.

[guang918]

Thanks very much for the reply. I am still confusion a little bit.

When Samtools call a sequence variant a SNP with AF1=0.5. That sequence variant may represent 20% or 40% or even 60% of the reads at that site. Therefore, we can't differentiate sequence error or a real SNP in this case. Is this thought correct?

Thanks.

[swbarnes2]

SAMTools mpileup pretty much assumes that you are looking at a single diploid organism. It will not ever tell you that it thinks your sample is 40% some alternate letter. It will tell you that it thinks it's 50%, but have a lower quality score, because 40% doesn't quite look like 50%.

So you will either have to examine samtools mpileup output very closely to characterize SNPs at percentages other than 50% and 100%, or you will need some other software to call them.

[gringer]

The VCF files also include raw counts of forward, reverse, alternate forward and alternate reverse. If you're looking at heteroploidy on a haploid genome (or similar), then it's better to use the raw counts rather than the estimated allele frequency.