output differents value between samtools depth and pysam

Question

[electronsource75]

Hello

I use some pysam function to calculate the depth of sequencing and I check the result with samtools depth. Most of the time the result is identical but sometimes there is a shift in the between the samtools result and the pysam result.

I put my code in a zip file with the example and a readme that contains the command to execute and the result files

https://gist.githubusercontent.com/fab-genty/53a814433f6e3d3794d3f2b04c37de66/raw/d550ddbaa6f5e0b88f948df4f2bd5fd0fcc0864a/fct_depth_biostars.py

[WouterDeCoster]

A better place to share code would be a GitHub Gist. I'm not going to download and extract a zip :-)

My first guess is that you see this difference because samtools has a default maximum depth of 8000. I believe if you set -d 0 you disable the maximum (if you have a reasonable recent version of samtools).

[Pierre Lindenbaum]

please use http://gist.github.com or paste the code here, instead of this temporary link.

[electronsource75]

I made the gist but I was unable to join the bam file that provokes the error

[finswimmer]

Hello,

some more general problems that might came up if one try to compare programs which calculate read depths.

• Are any reads filtered out because there MAPQ is below a certain threshold?
• Will duplicates are counted?
• Any threshold for base quality?
• How are overlapping paired reads handled? (Count as 1 or 2?)
• Make sure you are comparing the same region. Keyword: 0-based vs. 1-based coordinates.

fin swimmer

Answer 1

[Devon Ryan]

You'll want to use bamfile.pileup(chr_name, stepper='nofilter'): rather than the default filtering function to get more similar results.

[electronsource75]

Thank you this is exactly what I wanted. The output is now similar to the samtools function depth

[WouterDeCoster]

If an answer was helpful you should upvote it, if the answer resolved your question you should mark it as accepted.