Mate Pairs In Pysam

Question

[bioLife]

Case: I am trying to get the mate of the reads that map to eg. chromosome 1 (chr1)

I use the pysam module to extract those reads. While I can extract the position of the mates(rnext), I cannot extract the mate itself (meaning the sequence, qualty etc.)

The below code returns Valuerror: mate not found.

sortedBamfile = pysam.Samfile("sortedmappedReads.bam", "rb")
for alignedread in sortedBamfile.fetch("chr1"):
    print sortedBamfile.mate(alignedread)

However when I print the position of the mate I don't get an error. Which means that the mate exists.

for alignedread in sortedBamfile.fetch("chr1"):
    print sortedBamfile.getrname(alignedread.rnext)

Do you have any ideas/hints? Thanks!!

[Philipp Bayer]

The piece of code that returns ValueError: mate not found, does it print any other mates first or does it immediately crash?

Also, for clarity, I think you mean print sortedBamfile.getrname(alignedread.rnext) instead of print sortedBamfile.getrefname(alignedread.rnext)

[bioLife]

It immediately crashes. It does not print any mate at all. Any ideas what could be wrong? You were right about the typo. I edited the question.

[Philipp Bayer]

Hmm, I'm a bit at a loss here - you could try using Picard's ValidateSamFile to see whether there's anything weird with your bam-file, or try and re-make the index (.bam.bai file) using samtools index.... Which mapper did you use to generate the sam-files? I know there can be problems when you use SOAP's soap2sam.pl

[bioLife]

I used bowtie2. The thing is that I am also trying to get discordant reads. Would that be an issue? (but still I can get the mate's position) .I will validate my sam file.

[John]

rnext is stored in the read you are currently processing. .mate(someread) is a totally different operation that goes and finds the mate somewhere in the BAM.

I'm not familiar with .mate() since i typically read the whole file and match up mates in memory (this is often a lot faster than using the index via .mate), however I would imagine that your error is something due to the fact that pysam doesn't know how to find your mates. Consider indexing the bam, and posting a few lines of the bam (ideally a mate-pair) so we can see why pysam is struggling

EDIT: I just noticed this thread is two years old, lol.

Answer 1

[horticulturesh]

Hello, I am using pysam these days too but I face a problem. I assign a position to get reads aligned to that position and then I use is_read1 and is_read2 to estimate it is read1 or read2. If it is read1, how can I get its corresponding read2's position?

Thanks a lot, FredSamhaak