# How to filter intervals (reads or genomic coordinates) that have the exact same 5' or 3' ends?

Say I have reads that overlap some genes that produce small RNAs, but I want only those reads that start at exactly the TSS of the loci. In other words, reads whose 5' end match the 5' end of a genomic feature.

 5'....3' ...+++++++... Gene ...0000000... R1 ...00000000.. R2 ..000000000.. R3 ..00000000... R4 

The output should be R1 and R2.

I believe this is a fairly common bioinformatics operation but somehow it doesn't seem to be an option in the general use tools I looked at (py-bedtools, samtools, bedops).

A solution I thought of would be:

1. reduce the gene coordinates to its 5' end, and
2. using bedtools bedClosest to annotate the distance between reads and genes, and finally
3. select reads that overlap (d=0) and write them to a (bam) file.

This involves quite a bit of wrangling data with pybedtools or awk/bash, and I wonder is there is a more elegant solution to this. Solutions in R/Bioconductor are also welcome.

• This method won't work as overlapping segments will also have a distance of 0. You will have to clip the read coordinates as well as the gene coordinates. From the bedtools documentation: "-d In addition to the closest feature in B, report its distance to A as an extra column. - The reported distance for overlapping features will be 0." Aug 15 '17 at 18:29

Presuming you have deepTools installed:

#/!usr/bin/env python
from deeptoolsintervals.parse import GTF
import pysam

anno = GTF(["some/gtf/or/bed/file", "or/more/than/one/if you like/"])
bam = pysam.AlignmentFile("something.bam")
e = s + 1
overlaps = anno.findOverlaps(read.reference_name s, e, matchType=4)
# For a GTF, parse the results to see if this is the correct type.
# For a BED3/BED6 file there will be no exons, so anything other than
# None indicates an appropriate overlap

do something


deeptoolsintervals in the GTF/BED parsing and custom interval tree library used by deepTools. anno holds the interval tree and the associated information (gene/transcript name or ID, exons, strand, etc.). the matchType=4 is the pivotal part since that specifies that only overlaps having the exact same start position should be included. Assuming you have a stranded library, you might instead do something like switch between matchType=4 (the same start position) and matchType=5 (the same end position) depending on the orientation of the alignment and the the strand= of the gene.

Note that you'll need to make a set() to hold read names to find the mates (that or jump around the file to get them, which is probably slower).

• Hey @devon-ryan, is this documented somewhere? I searched the DeepTools docs but did not find it. I am looking for other, more tuned, ways to overlap reads with features. Apr 20 '18 at 13:17
• @fridaymeetssunday This is a semi-separate module, you can find the documentation for it here. It's shipped with deepTools but I never included it as part of the standard API documentation, since I figured it'd be internal use only. Apr 20 '18 at 13:24
• Perfect. I was looking to see if something like the option GTF_MATCH_WITHIN existed (it does :)), though in the meantime I figure that good ol' bedtools intersect will solve my problem for now. Apr 20 '18 at 13:43
• You can also do that with findIntersect() in R, I believe. Apr 20 '18 at 14:55