Extracting unmapped pairs of mapped reads from a particular chromosomal position

[Xespok]

Hi,

I would like to pool the unmapped pairs of mapped reads to a particular chromosomal location.

Many forum posts address the question of extracting read pairs with one or both reads unmapped using samtools view -f4 -F8 and -f8 -F4 or -f12 combinations. For the whole genome these work well.

But if I try to define a location to narrow down to reads in a particular chromosomal location for example chr1:1-1000, I end up with only the mapped reads in the resulting bam file, the unmapped pairs disappear.

One can get the header info of the resulting file and than go back to the bam file, convert it to some more readable format, and using this information get the unmapped pair information with custom scripts, but with a several GB file this is very slow. Since this needs to be performed several thousand times, this would take ages.

Can anyone suggest an easier way to accomplish this task?

[dpryan]

Just script something with python/perl/whatever-language-you-prefer. If that's not fast enough there's a C API for samtools that's pretty convenient. Those would allow a single pass through the file. One word of advice would be to process reads as pairs, rather than relying on the PNEXT field of unmapped reads (since whether the mate really falls into your target region isn't always discernible from that).

[Xespok]

Thanks, that is what I feared. Unfortunately, my scripting abilities are not sufficient for applying any C API at this moment.