Depth Of Coverage Reported By Igv And Samtools Do Not Agree

Question

[Anjan]

Hi all, I have several samples that have been sequenced on a MiSeq platform using a paired end protocol. Each read set of each sample was individually aligned using BWA (to the relevant reference sequence). Paired read sets were merged using samtools (option sampe) to produce the .sam files. The .sam files were then compressed, sorted and indexed (again using samtools) to produced the .bam; .sorted.bam and .sorted.bam.bai files. Now, I use the mpileup command to create the pileup file (sample command: samtools mpileup -B -f $reference $sorted_bam_file > $out_pileup). Here is the problem: When I visualize a region of the alignment using IGV, the depth of coverage reported by IGV does not agree with that reported by the mpileup command in any given position. Has anyone else faced this problem? Is there a solution? TIA, Anjan

[Manu Prestat]

Could you try Tablet to see if it gives you a 3rd answer? ;-) http://bioinf.hutton.ac.uk/tablet/

Answer 1

[Sean Davis]

Samtools mpileup will ignore bases with Q<13. You could try using -Q 0 to disable filtering of low-quality bases. Also, note that samtools mpileup will ignore all duplicates. Again, this will affect the reported depth. I do not know of a way to disable this feature (and really wouldn't want to, anyway).

Edit: As pointed out by others, there are filters both on the samtools side and the IGV side. In general, when comparing the outputs of software, one needs to first determine if the outputs should be the same. In this case, there is good reason to believe that they should not be identical.

[Jorge Amigo]

to complement this answer a little bit, this probably has to do with filters in general, on both sides of the road. Sean has pointed out a filter applied in samtools' mpileup by default, but also IGV implements filters by default (like not showing reads with mapping quality < 10). just make sure you are not using any filtering or that you are using the same filters, in case you want to see the same information on both programs.

[Biomonika (Noolean)]

I can only agree with both posts above, I have same experience.

[Anjan]

Sean, Thanks for the feed back. Indeed I am using the default filter in the mpileup command to get rid of low quality bases. However I am looking for documentation for your second assertion that "mpileup" also removes duplicates. The samtools manual does not explictly state this. Thanks again. Anjan

Answer 2

[-_-]

I am just cross-referencing another thread from SEQAnswers, http://seqanswers.com/forums/showthread.php?t=19702

Quote: "I found out that samtools filters reads before including them in the pileup; it reads the flag field in the bam file and discards reads that

1. are not paired
2. not properly mapped
3. mate is not mapped
4. alignment is not primary
5. reads fail quality control of vendor
6. is marked as PCR duplicates.

If the filters (1) and (3) are not desired, you can use the parameter -A. In addition, samtools performs realignment unless the parameter -B is used, and discards low quality reads unless -Q0 is used. Finally, it stops at a certain number of reads unless the -d parameter is invoqued."

Answer 3

[User6891]

I encountered the same problem as stated by Anjan. When adding the -A option to samtools mpileup, the per base coverage reported by samtools was the same as seen in IGV.

Command I used:

samtools mpileup -ABQ0 -d 1000000000 -f ref.fasta -l exome.bed mybam.bam

Answer 4

[fatemeh.dorri]

Anjan, I encounter the same issue. I am looking for a way to set samtools to consider duplicate reads. I was wondering if you have found how to do that? Thansk, Fatemeh