This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Depth Of Coverage Reported By Igv And Samtools Do Not Agree

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

next-gen igv samtools
• 17,041 views • link
updated by -_- -_- • written by Anjan Anjan

4 answers

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.

• 318 views • link
written by Sean Davis Sean Davis

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.

• 0 views • link
written by Jorge Amigo Jorge Amigo

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

• 0 views • link
written by Biomonika (Noolean) Biomonika (Noolean)

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

• 0 views • link
written by Anjan Anjan

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."

• 0 views • link
written by -_- -_-

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

• 0 views • link
written by User6891 User6891

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

• 0 views • link
written by fatemeh.dorri • fatemeh.dorri

Log in to answer this question.

More posts like this