Hope Townsend
unread,Oct 25, 2022, 4:50:18 PM10/25/22Sign in to reply to author
Sign in to forward
You do not have permission to delete messages in this group
Either email addresses are anonymous for this group or you need the view member email addresses permission to view the original message
to Pysam User group
Hello,
I am using pysam & pysamstats to get the coverage count on each strand, but it seems that for the negative strand I am facing some problems where pysam is stating that there are reads when IGV is not according to the bedgraphs originating from the same bams.
When looking at the data, the bedgraph makes more sense based on the genome annotation and since pysam seems to be saying that the positive and negative strands are having extremely similar coverage.
For example, I did both a direct pysam.pileup and pysamstats_coverage_strand for a region that IGV showed read counts of 2 on the negative strand and read counts ~ 100 on the positive strand. But pysam is continuously showing that the regions have counts of about 45 split between the two strands.
I've done this for more than one region, all seeming to show the same trend.
Example usage: for rec in pysamstats.stat_coverage_strand(sam_file, chrom='NC_000913.3',
start=450, end=500,
pad=True, truncate=True):
The sam_file is produced by pysam.AlignmentFile(bam_file_name, "rb")
When printing rec, the reads_rev_pp & reads_fwd_pp are what I look at (no difference between pp & not pp in my case)
Any guidance or insight as to why this might be happening would be appreciated, thank you!!!