Counting reads within the intron

Question

I have paired-end RNA-seq samples and a list of intron coordinates in bed format. I want to count the intronic reads such that:

1. The reads should overlap the intron by atleast 25 bp 
2. If both the reads of the pair are within the intron, count it as 1.

What tool can I use to compute it.

Thanks

Answer 1

You could use BEDOPS bedmap to map reads to introns, using 1) the --count operator to do counting of reads overlapping by your criteria; and, 2) the --indicator operator to do a true/false operation, where reads are contained entirely within the intron.

For instance, to count reads that overlap introns by at least 25 bases, use --count to count and --bp-over 25 to specify the minimum overlap threshold:

$ bedmap --echo --count --bp-ovr 25 --delim '\t' introns.bed reads.bed > answer.counts.bed

To filter for introns where reads are contained entirely within the intron, use --indicator to get introns that meet the criteria of full containment of reads, via --fraction-map 1:

$ bedmap --indicator --echo --fraction-map 1 --delim '\t' introns.bed reads.bed | awk '($1==1)' | cut -f2- > answer.indicator.bed

Once you have answer.counts.bed and answer.indicator.bed, you can compare those two results with bedmap --indicator --exact to see where there are commonalities (true/1) and differences (false/0), and increment the final count where there are identical introns, e.g.:

$ bedmap --indicator --echo --exact --delim '\t' answer.counts.bed answer.indicator.bed | awk '{ $NF += $1; print $0; }' | cut -f2- > answer.final.bed

If your reads are in a different format, i.e. BAM, you can use bam2bed in a process substitution, replacing reads.bed with <(bam2bed < reads.bam), e.g.:

$ bedmap --echo --count --bp-ovr 25 --delim '\t' introns.bed <(bam2bed < reads.bam) > answer.counts.bed

Etc.

You can also just convert separately, so that you only need to do it once:

$ bam2bed < reads.bam > reads.bed
$ bedmap --echo --count --bp-ovr 25 --delim '\t' introns.bed reads.bed > answer.counts.bed

Etc.

Answer 2

You can almost do what you want using featureCounts from the subread package. The only problem is that featureCounts requires GTF rather than bed.

If you had a GTF of your introns, you could do:

featureCounts -a introns.gtf my_bam.bam -f -p --minReadOverlap=25 -o counts.tsv

---

The problem with this then, is getting the GTF when you currently have a BED. One way to get the GTF would be to use CGAT.

cgat gtf2gtf -I ensembl.gtf.gz --method=merge-exons -S genes.gtf.gz
cgat gtf2gtf -I genes.gtf.gz --method=exons2introns -S introns.gtf

or you could put this together into one pipeline:

cgat gtf2gtf -I ensembl.gtf.gz --method=merge-exons  | cgat gtf2gtf --method=exons2introns -S introns.gtf

This will first merge all the exons of a gene together to make a single, collapsed transcript. It will that invert this to find that parts that are in the introns of this collapsed transcript. What you end up with is parts of the gene that are never exonic in any of the transcripts.

alternatively, you could convert your bed regions into GTF:

cgat bed2gff -I introns.bed -S introns.gtf --as-gtf