text: user2040: 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
  
  user235: 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
  
  
  user776: 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.
  
original_text: 2040: <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>
   235: <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>
   776: <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>
  