Extract R2F1 reads

Question

[IP]

Hi everyone!

I want to extract all reads from a bam file that have the orientation R2F1, as in the following image:

.

I have try to do this with samtools or pysam and I have not been able. I have seen in other posts some scripts for doing this with C++ or java, but, besides I don't know how to code in this languages, I want this to be a part of a larger script in python/bash, and thus, it will be great to have it in python/bash

[h.mon]

What did you try with samtools or pysam? Or you haven't been able to try at all?

[IP]

Or you haven't been able to try at all?

I would not have asked if I haven't tried :)

I have try both samtools and pysam. With samtools, I have try to extract them using the bitwise flag (145 flag) but the output is not the expected one, it includes reads supporting inversions (<- <-, FF orientation). With pysam,l I have tried the following, with very few reads in the output (not expected):

for read in bam.fetch('chr1',195411284,195423136):
    if read.is_reverse == True and read.is_second == True:
        if read.mate_is_reverse == False:
            #save the read and its mate to file

Also, I have tried using gridss ExtractSVreads, but it includes reads that support inversions (<- <-, FF orientation), which I don't want.

Thanks

Answer 1

[IP]

I am the OP,finally I have been able to do this efficiently with pysam doing the following:

NB: The input is a query name sorted bam, this allows to iterate the bam faster

bam = pysam.AlignmentFile("/path/to/input/input.bam", "rb")
tandem_dups = pysam.AlignmentFile("/path/to/output/output.bam", "wb", template=bam)
read1 = ''

for read in bam:
if type(read1) != type(read):
    if read.is_paired and read.is_read1:
        if read.mate_is_reverse and read.is_reverse == False:
            if read.reference_start > read.next_reference_start:
                read1 = read

else:
    tandem_dups.write(read1)
    tandem_dups.write(read)
    read1 = "None"