Samtools/bcftools problem

[abelhj]

I'm using samtools to look for SNPs with ultradeep (>18,000X) coverage. For some reason samtools and bcftools does not call some of the SNPs.

gi|251831106|ref|NC_012920.1|	16125	.	G	.	99	.	DP=18552;AF1=0;CI95=1.5,0;DP4=832,646,0,1;MQ=35;PV4=0.44,0.14,1,1	PL	0
gi|251831106|ref|NC_012920.1|	16126	.	T	C	99	.	DP=18770;AF1=1;CI95=1,1;DP4=28,21,777,596;MQ=35;PV4=1,0,0.00042,0.28	PL	219,255,0

For example, at position 16215, IGV gives me:
Total count: 18522
A: 17030 (92%, 8124+, 8906-)
C: 0
G: 1516 (8%, 861+, 655-)
T: 5 (0%, 4+, 1-)
N:1 (0%, 1+, 0-)

Sam/bcftools should be making a SNP call at position 16215 but is not. Also, the first two counts in DP4 almost agree between bcftools and IGV, but the third and fourth do not. The sequences are quality trimmed and the reads containing the A allele are not of low quality.

Anyone know why? Thanks.

[cboustred]

Did you get any response to your post?

I have similar data with ultra deep coverage where the total depth counts seen in IGV match the bcftools DP but the alternate allele counts are close but do not match those in the DP4 bcftools output.

I assume this is to do with the quality of the calls but I thought I had disabled this with -BQ0 when generating the mpileup

Any ideas?

Thanks

[dariober]

I'm using samtools to look for SNPs with ultradeep (>18,000X) coverage. For some reason samtools and bcftools does not call some of the SNPs.

Hi,

Can you post the commands of samtools/bcftools you are using? If you use vcfutils.pl varFilter there is a -D option that does not call SNPs exceeding D depth (the default is 10000000 so probably this is not the problem if you didn't change it!)

Sorry I can't be of more help...

All the best
Dario

[cboustred]

Hi, thanks for the quick response

I am simply running

samtools mpileup -l site.txt -BQ0 -d1000000 -uf ref.fasta test.sorted.bam > test_mpileup.bcf

bcftools view -Acg test_mpileup.bcf

I am just after the number of reads for each allele at a specific position, like the display you can get in IGV that the OP showed

[swbarnes2]

Maybe try:

mpileup -l site.txt -BQ0 -d1000000 ref.fasta test.sorted.bam > test_mpileup.pileup

To see what the quality of your alternate letters are according to mpileup. Then look at the individual sam entries for those reads. Maybe there's a discrepancy there.

[cboustred]

Thanks for the reply, I've just come across VarScan (http://varscan.sourceforge.net/index.html) which works on the pileup output from samtools.

Using this I am able to get the counts per allele per position in the pileup which also matches the IGV calls

All the best

Chris