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.
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
bedClosestto 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
awk/bash, and I wonder is there is a more elegant solution to this. Solutions in
R/Bioconductor are also welcome.