The Intron Exon Retention Estimator (IntEREst) facilitates estimation and comparison of splicing efficiency of transcripts across several samples. In particular, it can estimate the intron-retention levels or the exon junction levels in the transcripts. Our method estimates the Intron-retention by counting the number of rna-seq reads that have been mapped to the Intron-exon junctions of the genes, and it can estimate the exon junction levels by counting the reads that have been mapped to the exon-exon junctions. In addition, it is possible to limit the analysis to reads that are mapped to the intron-exon or exon-exon junctions only (See the junctionReadsOnly parameter in the interest() and interest.sequential() functions). However, by default this limitation is not taken into account, i.e. the reads that are fully mapped to the introns or exons are also considered. The method is similar to the Intron retention analysis used by Niemelä et al.(2014). The package accepts standard BAM files as input and produces tab separated text files together with SummarizedExperiment objects as results. To improve the performance and running time, the processing of each single BAM file can be distributed and run on several computing cores. The results can also be plotted and statistically analyzed to check the distribution of the intron retention levels, and compare the retention levels of U12 type introns to the U2 type across the studied samples. Note that although we mainly use this package to compare the splicing efficiency of the transcripts containing U12-type introns, some functions can be used with U2-type introns as well. The functions u12NbIndex(), u12Index(), u12Boxplot(), u12BoxplotNb(), u12DensityPlot()and u12DensityPlotIntron() are specifically used for the splicing efficiency analysis of the transcripts with U12-type introns. A diagram of the running pipeline is shown in figure 1.
Figure 1: Diagram of IntEREst running pipeline
Creating reference
The first step is to build a reference which will be used for the summarization of the sequencing reads, and for downstream intron retention and/or exon-exon junction analysis. This can be carried out by referencePrepare the function. The resulted reference includes coordinates of introns and exons of the genes and they can be extracted from various sources i.e. UCSC, biomaRt or a user defined file (e.g. GFF3/GTF). The exons with overlapping genomic coordinates can be collapsed (if collapseExons parameter is set as TRUE) to avoid assigning reads mapping to any alternatively skipped exons to their overlapping introns. An example of the process is shown in figure 2.
Figure 2: Running reference isoforms a, b, c and d leads to gene A with exons that are results of collapsing of exons of all isoforms.
Here we build a reference data frame from a manually built GFF3 file that includes exonic coordinates of the gene RHBDD3.
## Annotating U12 type introns {#annoU12} It is possible to annotate the U12 type introns in a reference using the annotateU12 function. U12-type introns (also known as minor type introns) are detected and spliced by the U12 splicing machinery as opposed to the majority of the introns (known as major type or U2 type) which are spliced by the U2 splicing machinery. U12-type introns also feature evolutionary conserved splice sites which are distinguished from the splice sites of U2 type introns hence they can be detected by mapping a position weigh matrix (PWM) to their spice sites and measuring their match score based on the PWM. The following scripts re-annotates introns of the genes RHBDD2 and YBX2.
Intron retention and exon-exon junction level estimation
The normalized intron retention levels and exon junction levels can be estimated using any of the two RNA-Seq read summarization functions, interest() and interest.sequential(). The interest() function is more robust since it distributes the reads in the .bam file over several computing cores and analyze the distributed data simultaneously. Note that regions in the genome with repetitive sequence elements may bias the mapping of the read sequences and the retention analysis. If you wish to exclude these regions from the analysis you can use the getRepeatTable() function, however We did not find repetitive DNA elements in particular biasing our results therefore we do not routinely use this function. As for instance, if you wish to exclude the coordinates in the genome housing Alu elements, you can run the reads summarization functions with the repeatsTableToFilter= getRepeatTable(repFamilyFil= "Alu") parameter setting. Also, to only consider the reads that map to intron-exon or exon-exon junctions set junctionReadsOnly= TRUE, however we recommend setting junctionReadsOnly= FLASE when measuring the intron retention levels (i.e. method=IntRet) and setting junctionReadsOnly= TRUE when measuring the exon-exon junction levels.
The interest() and interest.sequential() functions write output text files, additionally they can return a summarizedExperiment object for every sample they analyze. As shown in the following test script we, however usually prevent the individual runs to return any objects (by setting the returnObj=FALSE); instead, after running the analysis for all samples we generate a single summarizedExperiment object that includes results of all analyzed samples. To build such object from the output text files the readInterestResults() function can be used. In the following scripts a bam file from a single MDS sample with mutated ZRSR2 is used which includes all the reads mapped to the gene RHBDD3 only. We run two analysis that results to the number of reads mapping to the introns of the gene RHBDD3 and the number of reads mapping to the exon-exon junctions. eventually a SummarizedExperiment object is built for each of the 2 analysis that includes the read counts together with the coordinates and annotations of the introns and exons. The same analysis can be run on multiple .bam files to obtain SummarizedExperiment objects that include results for all analyzed .bam files .
As a demo we ran the IntEREst pipeline on 16 .bam files that each includes reads mapped to U12 genes located in chromosome 22. These bam files were results of mapping RNA-Seq data from bone-marrow samples published by Madan et al. (2015) to the Human genome (hg19). The studied samples were extracted from 16 individuals; out of which 8 were diagnosed with Myelodysplastic syndrome (MDS) and featured ZRSR2 mutation, 4 were diagnosed with MDS but lacked the mutation (referred to as ZRSR2 wild-type MDS samples) and 4 were healthy individuals.
The data is accessible through GEO with the accession number GSE63816 and the scripts that we ran to map the RNA-seq data, modify the bam files, extract the reads mapped to U12 genes in chr22 and build mdsChr22Obj and mdsChr22ExObj objects have been listed in the readme.txt file in scripts folder of the IntEREst package. You can get its full path using this script in R: system.file("scripts","readme.txt", package="IntEREst"). The mdsChr22Obj object is a summarizedExperiment object that includes retention levels of the introns of the genes located in the Chromosome 22 that feature at least one U12-type intron, across the 16 MDS samples. The mdsChr22ExObj object contain the exon-exon junction levels. They include two assays: counts and scaledRetention. Both can be accessed using functions with the same names: counts() and scaledRetention(). The former (counts) returns a data frame which includes the read counts of each intron/exon in each sample, and the latter (scaledRetention) returns a data frame with similar dimensions that includes the FPKM normalized read counts. The result objects also include intron/exon and sample annotations that can be retrieved using rowData() and colData() functions.
## DataFrame with 6 rows and 3 columns
## resultFiles type test_ctrl
## <character> <factor> <factor>
## SRR1691633 ./SRR1691633_ZRSR2Mut/interestRes.tsv ZRSR2mut test
## SRR1691634 ./SRR1691634_ZRSR2Mut/interestRes.tsv ZRSR2mut test
## SRR1691635 ./SRR1691635_ZRSR2Mut/interestRes.tsv ZRSR2mut test
## SRR1691636 ./SRR1691636_ZRSR2Mut/interestRes.tsv ZRSR2mut test
## SRR1691637 ./SRR1691637_ZRSR2Mut/interestRes.tsv ZRSR2mut test
## SRR1691638 ./SRR1691638_ZRSR2Mut/interestRes.tsv ZRSR2mut test
It is possible to plot() the object to check the distribution of the intron retention levels. The following scripts plot the average retention of all introns across the 3 sample types: ZRSR2 mutated MDS, ZRSR2 wild type MDS and healthy. The lowerPlot=TRUE and upperPlot=TRUE parameter settings ensures that both, the upper and lower triangle of the grid are plotted.
Figure 3: Plotting the distribution of the retention levels (\(log_e\) scaled retention) of introns of genes located on chromosome 22. The values have been averaged across the sample types ZRSR2 mutated, ZRSR2 wild type, and healthy.
The following script plots the average retention of the U12 introns across the 3 sample types: ZRSR2 mutated MDS, ZRSR2 MDS wild type and healthy. By default the upper triangle of the grid is plotted only (lowerPlot=FALSE).
Figure 4: Plotting the distribution of the retention levels (\(log_e\) scaled retention) of introns of genes located on chromosome 22. The values have been averaged across the sample types ZRSR2 mutated, ZRSR2 wild type, and healthy.
Comparing intron retention levels across various samples
IntEREst also provides various tools to compare the retention levels of the introns or exon junction levels across various samples. Initially, we extract the significantly higher and lower retained introns by using exactTestInterest() function which employs the exactTest() function from the edgeR package, i.e. an exact test for differences between two groups of negative-binomial counts. Note that exactTestInterest() makes comparison between a pair of sample types only (e.g. test vs ctrl).
As shown in the previous analysis ~50% of U12-type introns (of genes on Chr22) are significantly more retained (i.e. stabilized) in the ZRSR2 mutated samples comparing to the other samples, whereas same comparison shows that only ~1% of the U2-type introns are significantly more retained. For more complex experiments such as comparing samples based on a user defined design matrix other differential expression analysis functions from edgeR package, e.g. Linear Model (GLM) functions, have also been implemented in IntEREst; glmInterest() performs GLM likelihood ratio test, qlfInterest() runs quasi likelihood F-test, and treatInterest() runs fold-change threshold test on the retention levels of the introns/exons. The following commands can be used to extract the data for introns/exons that their retention levels vary significantly across all sample types: ZRSR2 mutation, ZRSR2 wild type, and healthy.
Next, to better illustrate the differences in the retention levels of different types of introns across the studied samples, we first use the bopxplot() method to illustrate the retention levels of all U12-type and U2-type introns in various sample types, and then we use the u12BoxplotNb() function to compare the retention of the U12 introns to their up- and down-stream U2-type introns.
Figure 5: Boxplot of the retention levels of U12 introns vs U2 introns, summed over samples with similar annotations i.e. ZRSR2 mutation, ZRSR2 wild type, or healthy.
Figure 6: Boxplot of retention levels of U12 introns vs their up- and down-stream U2 introns across all samples.
The boxplot clearly shows the increase retention of U12-type introns comparing to all the U2 introns (figure 5) and in particular comparing to the U2-type introns located on the up- or down-stream of the U12-type introns (figure 6). It is also clear that the elevated level of intron retention with U12-type introns is exacerbated in the ZRSR2 mutated samples comparing to the other studied samples. In order to better illustrate the stabilization of the U12-type introns comparing to the U2-type, we plot the density of the log fold-change of the retention (ZRSR2 mutated v.s. other samples) of U12-type introns and compare it to the log fold-change values for randomly selected U2-type introns, and U2-type introns up- or down-stream the U12-type introns.
## Design matrix not provided. Switch to the classic mode.
Figure 7: Density plot of the log fold change of U12-type introns, random U2-type introns and U2 introns (up / down / up and down)stream of U12-type introns.
As shown in figure 7 (and computed after the plot), when comparing the ZRSR2 mutated samples vs the other samples, for all U2-type introns the most frequent log fold-change (median) is ~0 whereas this value for the U12-type introns is noticeably higher (~1.01). It is also possible to run a statistical test to see if the log fold-changes of U12-type introns (ZRSR2 mutated samples vs other samples) are significantly higher than the log fold-changes of U2-type introns. For this purpose we use the jonckheere.test() function, i.e. Jonckheere-Terpstra ordered alternative hypothesis test, from the Clinfun package.
##
## Jonckheere-Terpstra test
##
## data:
## JT = 11209, p-value = 0.001
## alternative hypothesis: increasing
The result of the Jonckheere-Terpstra test with 1000 permutation runs shows that when comparing the samples that lack the ZRSR2 mutation to the the ZRSR2 mutated samples, the null hypothesis that the log fold-changes of the retentions of U12-type and U2-type introns are equally distributed was rejected with p-value 0.001, while the alternative being that the values in the U12-type introns are higher compared to the U2-type.
Our recommended pipeline for differential intron retention analysis
After building the reference as described in the test scripts above, we recommend running interest() (or interest.sequential()) twice: once in intron retention mode with method=IntRet and junctionReadsOnly= FALSE parameter settings; and subsequentlay in exon-exon junction mode with method=ExEx and junctionReadsOnly= FALSE parameter settings. Two example SummarizedExperiment objects resulted from Intron-retention and exon-exon junction mode of interest() test runs are mdsChr22Obj and mdsChr22ExObj. These two objects are limited to genes on chr22 that feature at least one U12 type intron. As shown in the scripts we recommend merging the intron-retention and exon-exon junction objects and running the deseqIntEREst() as shown:
Madan, V., D. Kanojia, J. Li, R. Okamoto, A. Sato-Otsubo, A. Kohlmann, M. Sanada, et al. 2015. “Aberrant splicing of U12-type introns is the hallmark of ZRSR2 mutant myelodysplastic syndrome.” Nat Commun 6:6042.
Niemela, E. H., A. Oghabian, R. H. Staals, D. Greco, G. J. Pruijn, and M. J. Frilander. 2014. “Global analysis of the nuclear processing of transcripts with unspliced U12-type introns by the exosome.” Nucleic Acids Res. 42 (11):7358–69.
