Iterating through paired-end reads in samtools / pysam
3
0
Entering edit mode
3.0 years ago
olavur ▴ 130

I want to process both reads in a paired-end read simultaneously, but don't know how to do this efficiently.

I'm currently doing it the following way, using the pysam Python wrapper for for the htslib API (Samtools).

samfile = AlignmentFile(filename, 'rb')  # BAM file reader.
# Check that the read has a pair that is mapped and not a duplicate.
# Get the other read in the pair.


This is really slow however.

From what I understand, the samfile.mate(read) call makes the BAM file reader jump to another place in the BAM. So it sounds like it would be a good idea to have two readers open and using one to iterate through the reads and the other to fetch the mate read, but that doesn't speed things up. I don't know whether paired-end reads are adjacent in a sorted and indexed BAM file, in that I case I would be able to take advantage of that, and do something along the lines of:

reads = samfile.fetch()
for i in range(number_of_pairs):


Anyone know how I can do this in a way that won't take forever to run?

BAM SAM Python • 4.9k views
1
Entering edit mode

Anyone know how I can do this in a way that won't take forever to run?

sort your bam on query name

0
Entering edit mode

I tried this, and from what I understand you can't index a BAM that isn't sorted by position. This makes life quite difficult, as many tools need the BAM to be indexed. Thanks for the suggestion. I'm not sure how to make it work.

0
Entering edit mode

okay, why do you need to fetch the mate ? what is the aim of your project ?

0
Entering edit mode

I have linked-read data where the reads that are close together have the same barcode. I want to try to check that the paired reads have the same barcode, with the assumption that a barcode mismatch means there were spurious alignments. Among other things, I think this might help me get a better estimate of insert size, assuming that the spurious alignments give false predictions of insert size. I think I will put this project on hold as I feel I lack the understanding to do this properly, but nevertheless I think it would be useful to be able to obtain the mate pairs in an efficient and fast manner.

2
Entering edit mode
2.7 years ago
ggydush ▴ 20

My method puts reads into a dictionary and will output reads once a pair is found. Usually pairs are near each other so the dictionary shouldn't get that big.

from collections import defaultdict
import pysam

"""
Generate read pairs in a BAM file or within a region string.
"""
continue
else:
else:
else:

bam = pysam.AlignmentFile(filename, 'rb')
# do stuff

0
Entering edit mode

Thank you, an elegant solution. Have you used this method a lot?

0
Entering edit mode

I use this function anytime I need to iterate in pairs... I've used it on a number of BAMs.

Hope it helps!

0
Entering edit mode

This will not work for a BAM file with multiple alignments for the same read pair, am I correct?

0
Entering edit mode

I modified a bit based on the @gizone1 's reply. If convenient, pls take a trial:  def read_pair_generator(bam, region_string=None): """ Generate read pairs in a BAM file or within a region string. Reads are added to read_dict until a pair is found. """ read_dict = defaultdict(lambda: [None, None]) for align in bam.fetch(until_eof=True, region=region_string): if not align.is_paired: continue qname = align.query_name if qname not in read_dict: if align.is_read1: read_dict[qname][0] = [align] else: read_dict[qname][1] = [align] elif align not in read_dict[qname][0] and align not in read_dict[qname][1]: if align.is_read1: read_dict[qname][0].extend(align) else: read_dict[qname][1].extend(align) else: if align.is_read1: yield read_dict[qname][0], read_dict[qname][1] else: yield read_dict[qname][0], read_dict[qname][1] del read_dict[qname]  I change the name from "read" to "align" since each row is an alignment record instead of a read record. Instead of return one pysam.AlignmentSegment object, I returned a list of pysam.AlignmentSegment objects.

And you can call it like this: for read1_aligns, read2_aligns in read_pair_generator(ex_sam):

In this case, read1_aligns will be a python list including all the pysam.AlignmentSegment objects corresponding to this read1. Pls try this out and see if this works.

1
Entering edit mode
2.8 years ago
Michi ▴ 980

I agree that it should be a bit simpler, given that paired end bam files are so common.

I used the following solution, which is only applicable if the .bam file has been previously sorted by read names samtools sort -n!

samfile = AlignmentFile(filename, 'rb')  # BAM file reader.

continue

else:
continue

print("found a pair!")

0
Entering edit mode

Pierre Lindenbaum did suggest this in a comment above. Did you try it? Is it fast, and does it get all reads (assuming they're a proper pair, non-duplicate, etc.)? It may be worth the overhead of sorting by query name.

0
Entering edit mode
13 months ago
u3005579 • 0

This is based on @gizone1's answer.

def read_pair_generator(bam, region_string=None):
"""
Generate read pairs in a BAM file or within a region string.
"""
for align in bam.fetch(until_eof=True, region=region_string):
if not align.is_paired:
continue
qname = align.query_name
else:
else:
else:
else:


I change the name from "read" to "align" since each row is an alignment record instead of a read record. Instead of return one pysam.AlignmentSegment object, I returned a list of pysam.AlignmentSegment objects.

In this case, read1_aligns will be a python list including all the pysam.AlignmentSegment objects corresponding to this read1. Pls try this out and see if this works.

0
Entering edit mode
def read_pair_generator(bam, region_string=None):
"""
Generate read pairs in a BAM file or within a region string.
"""
bam.reset()
for align in bam.fetch(until_eof=True, region=region_string):
if not align.is_proper_pair:
continue
qname = align.query_name
# print(align, type(align))
align_list = [align]
# print(align_list, type(align))
else:
else:
else:
else:


Sorry, fixed something. This is the working one.