Question about best way to extraction mates of reads

Question

[user191006]

When making read extraction with regions in BAM file, also mates are considered. What is the most efficient way also to collect them?

[zengxi.hada]

Hi,

May I ask have you solve the problem? I also encounter the same problem. I want to extract the mate reads of specific reads even if the mate is unmapped. Actually, I have posted a question at How to extract unmapped mate for specific reads using pysam or other module in python?.

Thank you!

Best regards,

Xi Zeng

Answer 1

[Devon Ryan]

Just query-sort the files (samtools sort -n). You don't need to randomly access anything then.

[user191006]

This is not what I want, sorting is not the problem for me. I focused on mates, my aim is to collect all of the mates of reads, thanks for response.

[WouterDeCoster]

By sorting on read names mates will be following each other in the sam file.

[GenoMax]

If your definition of mates is not what we are expecting then you should spell your definition out. mates, to most, are R1/R2 reads from a fragment.

[user191006]

I am going over regions in BED file. I just focus on reads in the region at a time, I try to collect their mates even if they are not in the region.

[Devon Ryan]

If the mates aren't terribly close then you'll end up doing a lot of linear searches even if you get both mates into the same bgzf block. I suspect it'll be faster to query sort and search for overlaps with the BED file.

Answer 2

[Samuel Brady]

I would suggest performing your alignment step using BWA within the Speed-Seq suite. This will automatically extract discordant mate pairs (that map to different chromosomes, for example) in a separate bam file for you, which will enable you to look at them without having to go through a huge bam file of "normal" reads. You can load these into IGV or another browser and look at your discordant reads without having all of the "normal" reads bogging you down.

Speed-Seq can be downloaded here: https://github.com/hall-lab/speedseq.

Here is a simple example of how to run Speed-Seq alignment that will result in the files I just described: https://github.com/samuelwb/tumor-evolution/blob/master/Alignment/SpeedSeqAlign.sh.