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.


  • $\begingroup$ How are your reads stored? Are they BAM or are they BED. If they are BED how is the pairing encoded? $\endgroup$ Jan 12 '18 at 14:08
  • $\begingroup$ Hi, My reads are in bam file, I can convert them into bed file. I am trying to follow what you did with your written scripts. Although I am having issues following how you calculated effective length under null model. I am trying to first make a non overlapping set of introns just like DEXseq manual says for exons $\endgroup$ Jan 12 '18 at 17:04

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


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



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

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.