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.

  • 1
    $\begingroup$ 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." $\endgroup$
    – Bioathlete
    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")
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).

  • $\begingroup$ 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. $\endgroup$ Apr 20 '18 at 13:17
  • 1
    $\begingroup$ @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. $\endgroup$
    – Devon Ryan
    Apr 20 '18 at 13:24
  • $\begingroup$ 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. $\endgroup$ Apr 20 '18 at 13:43
  • $\begingroup$ You can also do that with findIntersect() in R, I believe. $\endgroup$
    – Devon Ryan
    Apr 20 '18 at 14:55

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.