Question: Pysam, how to calculate how many bases on the reference are covered by reads?
gravatar for Dgg32
7 months ago by
Dgg3290 wrote:

I am using Pysam to calculate how many bases of my reference contig (assembled from PacBio) are covered by the reads (Illumina sequenced from the same sample). I have tried both pileup and count_coverage. But their results are somehow different:


import pysam

input_file = mem_pe.bam"

samfile = pysam.AlignmentFile(input_file, "rb")


n = 0

for pileupcolumn in samfile.pileup("tig00000003", min_base_quality=0):

    n += 1

print (n)

Output is 8817


n = 0

### for 4 bases
for i in range (4):

    n += len([x  for x in samfile.count_coverage("tig00000003", quality_threshold = 0)[i] if x != 0])

print (n)

Output is 16620.

So, why the two numbers are different please?

And how can I calculate the amount of bases covered by reads by using samtools please?

pysam samtools bam • 277 views
ADD COMMENTlink written 7 months ago by Dgg3290
Please log in to add an answer.


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