inconsistency between samtools mpileup and igvtools

[ynwh]

Hi all,
I'm checking reads mapping coverage on rice genome using IGV but find something strange:
(1) I generate .tdf file for IGV using the following command:
igvtools count -z 10 -w 1 a.sorted.bam a.tdf rice.genome
where rice.genome is the genome file generated by igv

(2) I generate coverage file using samtools using the following command:
samtools mpileup -BQ0 -d10000000 -f rice.fasta a.sorted.bam >a.mpileup

However, I find regions where a.tdf shows reads mapped but a.mpileup does not. The following is an example:

a.mpileup (mpileup.jpg) suggests no reads mapped between 16097520 and 16103121

a.tdf (igvtools.jpg) suggests several regions between 16097520 and 16103121 are covered by reads.

Can anybody explain and solve this problem?

Many thanks

Hao

[swbarnes2]

IGV and mpileup are two different tools, they are probably filtering things differently.

For starters, I bet mpileup is not showing anomalous pairs, and I bet IGV is.

So start by getting some read names from IGV, and seeing what their sam entry looks like. Maybe then you can figure out why mpileup isn't putting them where IGV does.

[dariober]

I don't if it is relevant to your case... I think mpileup is hard-coded to skip reads with flag 1024 (duplicates).

Dario

[ynwh]

Thanks for your advices. Now what I want is to get accurate coverage on each base, what tools and what parameter settings should I use?

[GenoMax]

Thanks for your advices. Now what I want is to get accurate coverage on each base, what tools and what parameter settings should I use?

http://www.biostars.org/p/5165/

http://seqanswers.com/forums/showthread.php?t=1835