Question: Depth Of Coverage Reported By Igv And Samtools Do Not Agree
7
gravatar for Anjan
6.2 years ago by
Anjan770
United States
Anjan770 wrote:

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 samtools igv • 9.2k views
ADD COMMENTlink modified 4.2 years ago by -_-780 • written 6.2 years ago by Anjan770
1

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

ADD REPLYlink written 6.2 years ago by Manu Prestat3.9k
15
gravatar for Sean Davis
6.2 years ago by
Sean Davis25k
National Institutes of Health, Bethesda, MD
Sean Davis25k wrote:

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 COMMENTlink modified 6.2 years ago • written 6.2 years ago by Sean Davis25k
3

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 REPLYlink written 6.2 years ago by Jorge Amigo11k

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

ADD REPLYlink written 6.2 years ago by Biomonika (Noolean)3.0k
1

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 REPLYlink written 6.2 years ago by Anjan770
4
gravatar for -_-
4.2 years ago by
-_-780
Canada
-_-780 wrote:

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 COMMENTlink modified 4.2 years ago • written 4.2 years ago by -_-780
3
gravatar for User6891
5.3 years ago by
User6891250
Europe
User6891250 wrote:

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 COMMENTlink written 5.3 years ago by User6891250
0
gravatar for fatemeh.dorri
5.7 years ago by
fatemeh.dorri0 wrote:

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 COMMENTlink written 5.7 years ago by fatemeh.dorri0
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 2126 users visited in the last hour