text: user108: 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:
  
  reduce the gene coordinates to its 5' end, and 
  using bedtools bedClosest to annotate the distance between reads and genes, and finally 
  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.
  
  user77: 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).
  
original_text: 108: <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>
   77: <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>
  