Hello everyone. I am using pysam to fetch the reads overlapping a given region, however, I am finding that pysam returns a list of reads whose positions do not correspond with the region coordinates that I input. The BAM file is indexed with samtools index 'file.bam'
For example, I want to fetch this region: chr3:50617849-50618103
I do so with the following code (inside python command line):
The first coordinate printed is 50335428, which is outside the scope of my input region completely.
Does anybody know why this is happening? The BAM file was created from a STAR alignment of paired-end reads against the hg19 human reference genome in case that matters.
For example, I want to fetch this region: chr3:50617849-50618103
I do so with the following code (inside python command line):
Code:
import pysam
bam = pysam.AlignmentFile('file.bam')
reads = bam.fetch('chr3',50617849,50618103)
for read in reads:
print read.pos
Does anybody know why this is happening? The BAM file was created from a STAR alignment of paired-end reads against the hg19 human reference genome in case that matters.
Comment