<p>I don't know about tools, but i've used the following python code to calculate the ratio of reads that overlap the 5' or 3' ends of introns or that are spliced. We sum these across all introns in a gene set (we actaully use this for iCLIP analysis to see if RNA binding proteins bind pre-mRNA or spliced RNA).</p>

<pre class="lang-py prettyprint-override"><code>import pysam
from collections import Counter
from CGAT import GTF, IOTools

def calculateSplicingIndex(bamfile, gtffile, outfile):

    bamfile = pysam.AlignmentFile(bamfile)

    counts = Counter()

    for transcript in GTF.transcript_iterator(
            GTF.iterator(IOTools.openFile(gtffile))):

        introns = GTF.toIntronIntervals(transcript)

        for intron in introns:
            reads = bamfile.fetch(
                reference=transcript[0].contig,
                start=intron[0], end=intron[1])

            for read in reads:
                if 'N' in read.cigarstring:
                    blocks = read.get_blocks()
                    starts, ends = zip(*blocks)
                    if intron[0] in ends and intron[1] in starts:
                        counts["Exon_Exon"] += 1
                    else:
                        counts["spliced_uncounted"] += 1
                elif (read.reference_start &lt;= intron[0] - 3
                      and read.reference_end &gt;= intron[0] + 3):
                    if transcript[0].strand == "+":
                        counts["Exon_Intron"] += 1
                    else:
                        counts["Intron_Exon"] += 1
                elif (read.reference_start &lt;= intron[1] - 3
                      and read.reference_end &gt;= intron[1] + 3):
                    if transcript[0].strand == "+":
                        counts["Intron_Exon"] += 1
                    else:
                        counts["Exon_Intron"] += 1
                else:
                    counts["unspliced_uncounted"] += 1

    header = ["Exon_Exon",
              "Exon_Intron",
              "Intron_Exon",
              "spliced_uncounted",
              "unspliced_uncounted"]

    with IOTools.openFile(outfile, "w") as outf:

        outf.write("\t".join(header)+"\n")
        outf.write("\t".join(map(str, [counts[col] for col in header]))
                   + "\n")
</code></pre>

<p>Unfortunately this uses a bunch of libraries you may or may not have, including CGAT (for the GTF parser and IOTools package) and pysam.</p>

<p>Once you've got these statistics you can calculate the "splicing index" as 
the log2 ratio of 2 times the number of spliced reads divided by the number of reads overlapping the 3' and 5' ends of introns. </p>

CC BY-SA 3.0
2017-06-13T19:14:44.307Z
2017-06-13T19:14:44.307Z
Ian Sudbery
<p>Take this first comment with a grain of salt, since this isn't an area I've worked in much, but: is binary classification possible? If a gene has 3 introns, and 2 are spliced out but 1 is retained, is this "spliced" or "unspliced". My first impression is that an analysis would be a bit more nuanced than a binary classification.</p>

<p>That said, I'm aware of a tool called Keep Me Around for intron retention analysis (<a href="https://arxiv.org/abs/1510.00696" rel="nofollow noreferrer">preprint</a>, <a href="https://github.com/pachterlab/kma" rel="nofollow noreferrer">code</a>). I've never used this software, but it looks like it's actively maintained and was created by a research group with a lot of experience and influence in this area, for whatever that's worth. :-)</p>

<p><strong>UPDATE</strong>: It looks like the original post was modified since I submitted my answer, but I'll leave the text about binary classification around for posterity. :-)</p>

CC BY-SA 3.0
2017-06-09T16:48:33.833Z
2017-06-09T16:48:33.833Z
Daniel Standage
<p>Is there a computational tool for measuring what percentage of RNA is spliced in an RNAseq experiment?</p>

<p>I'm not particularly interested in complicated analyses that give ratios for all possible alternative splicing variations. I'd rather have a binary classification. For example: (RNA1 = 35% spliced, 65% unspliced)</p>

<p><strong>Edit: (Additional information)</strong></p>

<p>We're very early in the exploratory stage, so a specific hypothesis is still forming. But we're interested in the sequence content of lncRNAs, many of which are poorly spliced. So we want to get a sense for how prevalent unspliced version of these transcript are. </p>

<p>Because this is currently exploratory, estimations and noise are fairly acceptable. A very basic approach might be to calculate the ratio of reads that map to exon-to-exon junctions, verses those that map from exon-to-intro junctions. </p>

<p>A more advanced approach I'm also considering is <a href="https://combine-lab.github.io/salmon/" rel="nofollow noreferrer">Salmon</a>.</p>

CC BY-SA 3.0
2018-03-09T00:26:22.560Z
Did you ever settle on something? I have the same question as you.
kmace
CC BY-SA 3.0
2018-03-09T03:14:59.720Z
We ended up using aligning with STAR then quantifying with salmon. We gave salmon all known splice variants, including the unspliced annotation. Then we summed all spliced isoforms to get "spliced vs. unspliced".

This pipeline was easy, accurate, and we were really satisfied with the result.
Jessime Kirk
<rna-seq><rna-splicing>
Spliced vs. unspliced ratios for transcripts in RNA-seq data
# Spliced vs. unspliced ratios for transcripts in RNA-seq data

(Asked by: username_0 on 2017-06-09)

Is there a computational tool for measuring what percentage of RNA is spliced in an RNAseq experiment?

I'm not particularly interested in complicated analyses that give ratios for all possible alternative splicing variations. I'd rather have a binary classification. For example: (RNA1 = 35% spliced, 65% unspliced)

**Edit: (Additional information)**

We're very early in the exploratory stage, so a specific hypothesis is still forming. But we're interested in the sequence content of lncRNAs, many of which are poorly spliced. So we want to get a sense for how prevalent unspliced version of these transcript are.

Because this is currently exploratory, estimations and noise are fairly acceptable. A very basic approach might be to calculate the ratio of reads that map to exon-to-exon junctions, verses those that map from exon-to-intro junctions.

A more advanced approach I'm also considering is [Salmon](https://combine-lab.github.io/salmon/).

**Tags:** rna-seq, rna-splicing
> **Comments:**
> **username_1**: Did you ever settle on something? I have the same question as you.  
> **username_0**: We ended up using aligning with STAR then quantifying with salmon. We gave salmon all known splice variants, including the unspliced annotation. Then we summed all spliced isoforms to get "spliced vs. unspliced".  This pipeline was easy, accurate, and we were really satisfied with the result.  

---

## 2 answers

### Answer by username_2 on 2017-06-13. Score: 3

I don't know about tools, but i've used the following python code to calculate the ratio of reads that overlap the 5' or 3' ends of introns or that are spliced. We sum these across all introns in a gene set (we actaully use this for iCLIP analysis to see if RNA binding proteins bind pre-mRNA or spliced RNA).
import pysam
from collections import Counter
from CGAT import GTF, IOTools

def calculateSplicingIndex(bamfile, gtffile, outfile):

    bamfile = pysam.AlignmentFile(bamfile)

    counts = Counter()

    for transcript in GTF.transcript_iterator(
            GTF.iterator(IOTools.openFile(gtffile))):

        introns = GTF.toIntronIntervals(transcript)

        for intron in introns:
            reads = bamfile.fetch(
                reference=transcript[0].contig,
                start=intron[0], end=intron[1])

            for read in reads:
                if 'N' in read.cigarstring:
                    blocks = read.get_blocks()
                    starts, ends = zip(*blocks)
                    if intron[0] in ends and intron[1] in starts:
                        counts["Exon_Exon"] += 1
                    else:
                        counts["spliced_uncounted"] += 1
                elif (read.reference_start <= intron[0] - 3
                      and read.reference_end >= intron[0] + 3):
                    if transcript[0].strand == "+":
                        counts["Exon_Intron"] += 1
                    else:
                        counts["Intron_Exon"] += 1
                elif (read.reference_start <= intron[1] - 3
                      and read.reference_end >= intron[1] + 3):
                    if transcript[0].strand == "+":
                        counts["Intron_Exon"] += 1
                    else:
                        counts["Exon_Intron"] += 1
                else:
                    counts["unspliced_uncounted"] += 1

    header = ["Exon_Exon",
              "Exon_Intron",
              "Intron_Exon",
              "spliced_uncounted",
              "unspliced_uncounted"]

    with IOTools.openFile(outfile, "w") as outf:

        outf.write("\t".join(header)+"\n")
        outf.write("\t".join(map(str, [counts[col] for col in header]))
                   + "\n")

Unfortunately this uses a bunch of libraries you may or may not have, including CGAT (for the GTF parser and IOTools package) and pysam.
Once you've got these statistics you can calculate the "splicing index" as 
the log2 ratio of 2 times the number of spliced reads divided by the number of reads overlapping the 3' and 5' ends of introns.


---

### Answer by username_3 on 2017-06-09. Score: 2

Take this first comment with a grain of salt, since this isn't an area I've worked in much, but: is binary classification possible? If a gene has 3 introns, and 2 are spliced out but 1 is retained, is this "spliced" or "unspliced". My first impression is that an analysis would be a bit more nuanced than a binary classification.

That said, I'm aware of a tool called Keep Me Around for intron retention analysis ([preprint](https://arxiv.org/abs/1510.00696), [code](https://github.com/pachterlab/kma)). I've never used this software, but it looks like it's actively maintained and was created by a research group with a lot of experience and influence in this area, for whatever that's worth. :-)

**UPDATE**: It looks like the original post was modified since I submitted my answer, but I'll leave the text about binary classification around for posterity. :-)


---

