Question: Iterating through paired-end reads in samtools / pysam
0
gravatar for olavur
6 months ago by
olavur60
Tórshavn, Faroe Islands
olavur60 wrote:

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.
# Iterate through reads.
for read in samfile:
    # Check that the read has a pair that is mapped and not a duplicate.
    if read.is_paired and not read.mate_is_unmapped and not read.is_duplicate:
        # Get the other read in the pair.
        read_mate = samfile.mate(read)

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):
    read1 = next(reads)
    read2 = next(reads)

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

sam bam python • 520 views
ADD COMMENTlink modified 6 weeks ago by gizone110 • written 6 months ago by olavur60
1

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

sort your bam on query name

ADD REPLYlink written 6 months ago by Pierre Lindenbaum112k

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.

ADD REPLYlink written 6 months ago by olavur60

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

ADD REPLYlink written 6 months ago by Pierre Lindenbaum112k

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.

ADD REPLYlink written 6 months ago by olavur60
1
gravatar for Michi
3 months ago by
Michi940
Barcelona
Michi940 wrote:

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.
# Iterate through reads.
read1 = None
read2 = None

for read in samfile:

    if not read.is_paired or read.mate_is_unmapped or read.is_duplicate:
          continue

    if read.is_read2:
        read2 = read
    else:
        read1 = read
        read2 = None
        continue

     if not read1 is None and not read2 is None and read1.query_name == readf2.query_name:
          print("found a pair!")
          ## do your stuff
ADD COMMENTlink modified 3 months ago • written 3 months ago by Michi940

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.

ADD REPLYlink written 9 weeks ago by olavur60
1
gravatar for gizone1
6 weeks ago by
gizone110
gizone110 wrote:

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


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 read in bam.fetch(region=region_string):
        if not read.is_proper_pair or read.is_secondary or read.is_supplementary:
            continue
        qname = read.query_name
        if qname not in read_dict:
            if read.is_read1:
                read_dict[qname][0] = read
            else:
                read_dict[qname][1] = read
        else:
            if read.is_read1:
                yield read, read_dict[qname][1]
            else:
                yield read_dict[qname][0], read
            del read_dict[qname]

bam = pysam.AlignmentFile(filename, 'rb')
for read1, read2 in read_pair_generator(bam):
    # do stuff
ADD COMMENTlink modified 4 weeks ago • written 6 weeks ago by gizone110

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

ADD REPLYlink written 15 days ago by olavur60

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

Hope it helps!

ADD REPLYlink written 1 day ago by gizone110
Please log in to add an answer.

Help
Access

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