Depth Of Coverage Reported By Igv And Samtools Do Not Agree
4
8
Entering edit mode
11.1 years ago
Anjan ▴ 830

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 • 14k views
ADD COMMENT
1
Entering edit mode

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

ADD REPLY
15
Entering edit mode
11.1 years ago

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.

ADD COMMENT
3
Entering edit mode

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.

ADD REPLY
0
Entering edit mode

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

ADD REPLY
1
Entering edit mode

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

ADD REPLY
4
Entering edit mode
9.1 years ago
-_- ★ 1.1k

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

ADD COMMENT
3
Entering edit mode
10.3 years ago
User6891 ▴ 330

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

ADD COMMENT
0
Entering edit mode
10.6 years ago

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

ADD COMMENT

Login before adding your answer.

Traffic: 2023 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6