<p>You could use <a href="https://github.com/bedops/bedops" rel="nofollow noreferrer">BEDOPS</a> <code>bedmap</code> to map reads to introns, using 1) the <code>--count</code> operator to do counting of reads overlapping by your criteria; and, 2) the <code>--indicator</code> operator to do a true/false operation, where reads are contained entirely within the intron.</p>

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

<pre><code>$ bedmap --echo --count --bp-ovr 25 --delim '\t' introns.bed reads.bed &gt; answer.counts.bed
</code></pre>

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

<pre><code>$ bedmap --indicator --echo --fraction-map 1 --delim '\t' introns.bed reads.bed | awk '($1==1)' | cut -f2- &gt; answer.indicator.bed
</code></pre>

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

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

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

<pre><code>$ bedmap --echo --count --bp-ovr 25 --delim '\t' introns.bed &lt;(bam2bed &lt; reads.bam) &gt; answer.counts.bed
</code></pre>

<p>Etc.</p>

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

<pre><code>$ bam2bed &lt; reads.bam &gt; reads.bed
$ bedmap --echo --count --bp-ovr 25 --delim '\t' introns.bed reads.bed &gt; answer.counts.bed
</code></pre>

<p>Etc.</p>

CC BY-SA 3.0
2018-01-11T21:22:03.560Z
Thanks Alex,
Can you explain what is --indicator and --fraction-map doing in generating answer.indicator.bed. I don't get this part?
user3138373
CC BY-SA 3.0
2018-01-11T21:59:22.760Z
I suggest reading the documentation, as it goes into much more detail: http://bedops.readthedocs.io/en/latest/content/reference/statistics/bedmap.html
Alex Reynolds
CC BY-SA 3.0
2018-01-11T19:06:07.157Z
2018-01-11T19:06:07.157Z
Alex Reynolds
<p>You can almost do what you want using <code>featureCounts</code> from the <code>subread</code> package. The only problem is that featureCounts requires GTF rather than bed.</p>

<p>If you had a GTF of your introns, you could do:</p>

<pre><code>featureCounts -a introns.gtf my_bam.bam -f -p --minReadOverlap=25 -o counts.tsv
</code></pre>

<hr>

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

<pre><code>cgat gtf2gtf -I ensembl.gtf.gz --method=merge-exons -S genes.gtf.gz
cgat gtf2gtf -I genes.gtf.gz --method=exons2introns -S introns.gtf
</code></pre>

<p>or you could put this together into one pipeline:</p>

<pre><code>cgat gtf2gtf -I ensembl.gtf.gz --method=merge-exons  | cgat gtf2gtf --method=exons2introns -S introns.gtf
</code></pre>

<p>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. </p>

<p>alternatively, you could convert your bed regions into GTF:</p>

<pre><code>cgat bed2gff -I introns.bed -S introns.gtf --as-gtf
</code></pre>

CC BY-SA 3.0
2018-01-14T12:48:07.493Z
2018-01-14T12:48:07.493Z
Ian Sudbery
<p>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:</p>

<pre><code>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.
</code></pre>

<p>What tool can I use to compute it.</p>

<p>Thanks</p>

CC BY-SA 3.0
2018-01-12T14:08:55.783Z
How are your reads stored? Are they `BAM` or are they `BED`. If they are `BED` how is the pairing encoded?
Ian Sudbery
CC BY-SA 3.0
2018-01-12T17:04:29.377Z
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
user3138373
<rna-seq>
Counting reads within the intron
# Counting reads within the intron

(Asked by: username_0 on 2018-01-11)

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

**Tags:** rna-seq
> **Comments:**
> **username_1**: How are your reads stored? Are they `BAM` or are they `BED`. If they are `BED` how is the pairing encoded?  
> **username_0**: 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  

---

## 2 answers

### Answer by username_2 on 2018-01-11. Score: 3

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.

> **Comments:**
> **username_0**: Thanks Alex, Can you explain what is --indicator and --fraction-map doing in generating answer.indicator.bed. I don't get this part?  
> **username_2**: I suggest reading the documentation, as it goes into much more detail: http://bedops.readthedocs.io/en/latest/content/reference/statistics/bedmap.html  

---

### Answer by username_1 on 2018-01-14. Score: 1

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


---

