Answers.Body: <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>
  
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2018-04-20T13:17:51.233Z
Answers.Text: 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.
Answers.UserDisplayName: fridaymeetssunday
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2018-04-20T13:24:04.710Z
Answers.Text: @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.
Answers.UserDisplayName: Devon Ryan
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2018-04-20T13:43:07.543Z
Answers.Text: 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.
Answers.UserDisplayName: fridaymeetssunday
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2018-04-20T14:55:08.823Z
Answers.Text: You can also do that with `findIntersect()` in R, I believe.
Answers.UserDisplayName: Devon Ryan
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-08-15T08:33:27.510Z
Answers.LastActivityDate: 2017-08-15T08:33:27.510Z
Answers.UserDisplayName: Devon Ryan
Body: <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>
  
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-08-15T18:29:52.350Z
Comments.Text: 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."
Comments.UserDisplayName: Bioathlete
Tags: <bedtools>
Title: How to filter intervals (reads or genomic coordinates) that have the exact same 5' or 3' ends?
ThreadText: # 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:
  
  ```py
  #/!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 @username_2-username_2, 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**: @username_0 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.  
  
  ---
  
  