variant calling on a subsample of a BAM file..

[batman]

Hi,

Can someone explain the following phenomenon. I ran bcftools (bcftools call -v -c )on a pileup created by samtools. This resulted in calling of a SNP. Next I extracted reads spanning the SNP, created the pileup and ran bcftools as before. I was expecting that bcftools would give me back the same result but was surprised when no SNP was predicted. Why does this happen?

[dpryan]

Samtools does some filtering of the reads and bases that you probably didn't do manually. Also make sure that you handled overlapping ends of pairs.

[batman]

Hi @dpryan,

What kind of filtering does samtools do: I used something like samtools view file.name chrY:XXXX-XXXX where XXXX is the SNP location. Would this not give all the reads used in calling the variant? does bcftools care about bases neighboring the SNP?

[dpryan]

Aside from adjusting the phred score on bases from overlapping pairs, it discards orphans, any bases with scores < 13 and any alignments marked as duplicate, secondary, or QC failed. So, samtools view is providing a superset of what mpileup would produce.

[batman]

Hi @dpryan,

Thanks for the very detailed answer.
But this still leaves me wondering, since I have the identical superset of reads spanning the SNPs, shouldnt mpileup and bcftools predict the same set of SNP for the whole BAM file and to the extracted reads? Only difference I can think of is that some reads do not have the paired read extracted.

[dpryan]

At least in the most recent versions of samtools, mpileup produced a gvcf, which includes basically everything. bcftools call then uses that as a source for calculating the more believable variants. Presumably it's this step that you're seeing.