How do I filter out reads where Read 2 doesn't map to the end of a fasta sequence

Question

[cag104]

For simplicity let us say I have one sequence that is 400bp in length. I used STAR to map reads to this sequence and now I am interested in only keeping the read pairs where read 2 maps to the end of the sequence (position ~440-450 of 450). What is the best way to do this? I've been trying to see if its possible using PySAM but am having some issues.

Thanks!

Answer 1

[swbarnes2]

I'd do this in two steps. 1, as above, use samtools with a bed filter to get the reads where read2 do map to the region you want.

Then pull out the read names of those reads. Then use picard filterSamReads to get the read1 and read2 data belonging to those reads.

Answer 2

[i.sudbery]

I'd definitely do this with pysam. What issues are you having? I'd do something like:

import pysam
inbam = pysam.AlignmentFile("my.bam")
outbam = pysam.AlignmentFile("my.filtered.bam", "w", template=inbam)

for read in inbam.fetch("sequence_name"):

    if read.reference_end >= 440:
        outbam.write(read)

outbam.close()