1
$\begingroup$

I'm encountering a discrepancy in the depth of coverage assessment for specific coordinates when comparing results obtained from pileup files and samtools depth.

I've processed pileup files and extracted coverage information for a given coordinate, resulting in a coverage depth of 8028. However, when using samtools depth on the same coordinate, the coverage depth is reported as 88565.

My analysis involves the following steps:

Alignment of paired-end reads using bwa mem. Generation of pileup files and extraction of coverage information. Utilizing samtools depth to calculate coverage depth for the same coordinates. Visual inspection in IGV for validation, where I observe that the number of reads aligning to the region matches the count obtained from samtools depth. I've noticed that in this region, there are no structural variations. The versions of samtools I'm using are samtools 1.19.2.

Could you please help me understand the potential reasons for such big discrepancies? What parameters should I review, and could the differences be attributed to any specific factors?

Thank you for your assistance!

Improve this question
$\endgroup$
2
  • $\begingroup$ Please clarify your specific problem or provide additional details to highlight exactly what you need. As it's currently written, it's hard to tell exactly what you're asking. $\endgroup$
    – Community Bot
    Commented Feb 3 at 21:58
  • $\begingroup$ How are you generating the pileup files? Can you please show the command line steps you use in your analysis to generate these pileup files & depth calculations? $\endgroup$
    – gringer ♦
    Commented Feb 4 at 18:55

1 Answer 1

Reset to default
1
$\begingroup$

This is purely a guess based on what you write given that you are not showing the exact command for samtools etc.

I suspect that you are allowing sequences with MAPQ==0 (supplementary alignments) into your depth calculation. A brief google search tells me that this issue has been previously observed with samtools depth.

In principle, you should be able to see this directly in IGV, by controlling how supplementary alignments are displayed.

This is something that you can directly control with the -Q flag, to restrict it to uniquely mapping reads (where you set the flag to e.g. -Q 1 or higher).

I'd suggest trying to play with samtools depth parameters and seeing if you can bring it more into alignment.

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.