<p>Presuming you have deepTools installed:</p>

<pre class="lang-py prettyprint-override"><code>#/!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")
for read in bam:
    s = read.reference_start
    if read.is_reverse:
        s = read.reference_end
    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
</code></pre>

<p><code>deeptoolsintervals</code> in the GTF/BED parsing and custom interval tree library used by deepTools. <code>anno</code> holds the interval tree and the associated information (gene/transcript name or ID, exons, strand, etc.). the <code>matchType=4</code> 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 <code>matchType=4</code> (the same start position) and <code>matchType=5</code> (the same end position) depending on the orientation of the alignment and the the <code>strand=</code> of the gene.</p>

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

CC BY-SA 3.0
2018-04-20T13:17:51.233Z
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.
fridaymeetssunday
CC BY-SA 3.0
2018-04-20T13:24:04.710Z
@fridaymeetssunday This is a semi-separate module, you can find [the documentation for it here](https://github.com/deeptools/deeptools_intervals). 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.
Devon Ryan
CC BY-SA 3.0
2018-04-20T13:43:07.543Z
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.
fridaymeetssunday
CC BY-SA 3.0
2018-04-20T14:55:08.823Z
You can also do that with `findIntersect()` in R, I believe.
Devon Ryan
CC BY-SA 3.0
2017-08-15T08:33:27.510Z
2017-08-15T08:33:27.510Z
Devon Ryan
<p>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. </p>

<p><code>
   5'....3'
...+++++++... Gene
...0000000... R1
...00000000.. R2
..000000000.. R3
..00000000... R4
</code></p>

<p>The output should be R1 and R2.</p>

<p>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). </p>

<p>A solution I thought of would be:</p>

<ol>
<li>reduce the gene coordinates to its 5' end, and </li>
<li>using bedtools <code>bedClosest</code> to annotate the distance between reads and genes, and finally </li>
<li>select reads that overlap (d=0) and write them to a (bam) file.</li>
</ol>

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

CC BY-SA 3.0
2017-08-15T18:29:52.350Z
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."
Bioathlete
<bedtools>
How to filter intervals (reads or genomic coordinates) that have the exact same 5' or 3' ends?
# How to filter intervals (reads or genomic coordinates) that have the exact same 5' or 3' ends?

(Asked by: username_0 on 2017-08-15)

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.

**Tags:** bedtools
> **Comments:**
> **username_1**: 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."  

---

## 1 answers

### Answer by username_2 on 2017-08-15. Score: 1 | Accepted answer

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")
for read in bam:
    s = read.reference_start
    if read.is_reverse:
        s = read.reference_end
    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).

> **Comments:**
> **username_0**: 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.  
> **username_2**: @fridaymeetssunday This is a semi-separate module, you can find [the documentation for it here](https://github.com/deeptools/deeptools_intervals). 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.  
> **username_0**: 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.  
> **username_2**: You can also do that with `findIntersect()` in R, I believe.  

---

