Can you read a bam file by mate pairs?

[kgulukota]

I need to process a bam file by reading mate pairs together i.e. each read of the bam file must return TWO alignments to analyze together.

I have been using the Bam API function bam_read_1 (http://samtools.sourceforge.net/samt...func/bam_read1) which returns the "next" alignment in a bam. And I know that the alignment object has some basic info about the mate, but I am looking for all of the mate's alignment.

I realize that a bam file is usually sorted by genome coordinates and therefore the mate pairs are not together and that "next pair" call to a function like bam_read_1 is ambiguous since there are two positions in the bam file that are being read.

Nevertheless, I wonder, is there a way to get a pointer from each alignment to that of its mate? (I know we then have to worry about whether the "next" alignment has already been read as some other read's mate, but still...)

This seems like it should be a common programming problem in reading bam files. Curious if there is any solution out there.

Thanks,
Gulu

[swbarnes2]

If you sort your file by name, instead of coordinate, the pairs will be next to each other.

[dpryan]

Just name-sort the file and then you can (normally) just take two reads at a time. I mention normally since some programs, such as tophat, will leave orphaned reads and leave them marked as paired.

Edit: I guess I should have updated before replying, since swbarnes2 beat me to it. I'll add that I do this in the methylation extractor in bison, which also uses the samtools C API.

[gt1]

I need to process a bam file by reading mate pairs together i.e. each read of the bam file must return TWO alignments to analyze together.

I have been using the Bam API function bam_read_1 (http://samtools.sourceforge.net/samt...func/bam_read1) which returns the "next" alignment in a bam. And I know that the alignment object has some basic info about the mate, but I am looking for all of the mate's alignment.

I realize that a bam file is usually sorted by genome coordinates and therefore the mate pairs are not together and that "next pair" call to a function like bam_read_1 is ambiguous since there are two positions in the bam file that are being read.

Nevertheless, I wonder, is there a way to get a pointer from each alignment to that of its mate? (I know we then have to worry about whether the "next" alignment has already been read as some other read's mate, but still...)

This seems like it should be a common programming problem in reading bam files. Curious if there is any solution out there.

Thanks,
Gulu

As the previous replies already state, one way to collate all the pairs is to perform a full sorting by query name. For most BAM files this is however a waste of computational resources, as even if the file is not sorted by name, the pairs will be close together. The biobambam package (https://github.com/gt1/biobambam, Linux binaries at ftp://ftp.sanger.ac.uk/pub/users/gt1/biobambam/) contains a tool named bamcollate which reads a BAM file from standard input and outputs a BAM file on standard output in which the ends of the pairs appear consecutively. It tries to find pairs which are close together in the file without sorting and only performs sorting for those which are far apart. There is an API also (see http://arxiv.org/abs/1306.0836).

[kgulukota]

Thank you @dpryan and @swbarnes2. But sorting by name won't work for me. My project involves examining reads mapped to a particular region - so I have to sort by coordinate in order to load index and use bam_fetch.

So, the real question is: when bam_fetch retrieves an alignment in the region of interest, is there a way to quickly skip to the mate's alignment?

Thank you @gt1. I will check it out.

[dpryan]

There's no direct function for that. The best you can do is bam_fetch the appropriate area and then iterate through the results until you find the mate. I presume that this is what the package that gt1 posted does.