Extract reads that are contained entirely within a give region

Question

[george.wiggins]

This seems like a very simple task to do and I'd assume that samtools view would be able to do this. Basically all I want to do is to return all reads that a entirely contained within a given region.

samtools view bamfile "chr1:100-1000" returns any read the overlaps that region. For example a read may start before 100 or extend beyond 1000. I only what those that map entirely inside that region.

Answer 1

[Alex Reynolds]

Use bedops and convert2bed (bam2bed) with a process substitution that specifies your region of interest:

$ bam2bed < reads.bam | bedops -e 100% - <(echo -e "chr1\t100\t1000") > answer.bed

Using bedops -e 100% asks for reads from read.bam that are contained entirely within chr1:100-1000.

Any reads that are found are written to answer.bed.

For RNA-seq data, you might wish to use the --split operand with bam2bed. See the documentation for more information.

References:

• https://bedops.readthedocs.io/en/latest/content/reference/set-operations/bedops.html#element-of-e-element-of
• https://bedops.readthedocs.io/en/latest/content/reference/file-management/conversion/bam2bed.html