Alternative Splicing (AS) is a regulation mechanism that contributes to protein diversity and is also associated to many diseases and tumors. Alternative splicing events quantification from RNA-Seq reads is a crucial step in understanding this complex biological mechanism. However, tools for AS events detection and quantification show inconsistent results. This reduces their reliability in fully capturing and explaining alternative splicing. We introduce ESGq, a novel approach for the quantification of AS events across conditions based on read alignment against Event Splicing Graphs. By comparing ESGq to two state-of-the-art tools on real RNA-Seq data, we validate its performance and evaluate the statistical correlation of the results. ESGq is freely available at https://github.com/AlgoLab/ESGq.
curate than the transcript-based competitors [
In the last decade, RNA-Sequencing has become the de- of the classical and more strict classification. However,
facto standard for the analysis of alternative splicing and the detection and quantification of classical AS events
a plethora of tools have been proposed in the literature. is an already hard - and not fully solved - problem that
From a very high level point of view, the approaches does not need further complications. Indeed, tools and
for the analysis of alternative splicing available in the methodologies show several limitations [
In this work, we will focus on the second category. ery tool uses its own format to describe the AS events,
This kind of approaches characterizes AS at the most making any downstream analysis quite complex.
ifne-grained level by giving a detailed and strict descrip- In this context, we focus on the detection and
quantion of what happens at the exon-exon (or splice) junction tification of non-novel AS events across conditions and
level. Although being more rigorous in the description provide an extensive comparison of two state-of-the-art
of AS events, this kind of approaches resulted more ac- tools, rMATS [
Experiments on a very recent real RNA-Seq dataset
sequenced from Drosophila Melanogaster flies at two time
points show that, by pairing simple graphs with accurate
mapping, ESGq is able to achieve comparable results to
state-of-the-art without sacrificing its eficiency.
Moreover, our comparison proves once again the inconsistency
of the results obtained by diferent methodologies for the
detection and quantification of AS events.
ESGq starts its computation by extracting the annotated
alternative splicing events from the input gene
annotation. To this aim, it employs a module of the SUPPA2
tool [
Diferently from state-of-the-art approaches, which rely on spliced read alignment to a reference genome or quasi-mapping transcript quantification, ESGq relies on read alignment against a graph-based structure representing the events that need to be quantified. Instead of using a full splicing graph or a pantranscriptome, that are common structures in the literature [13, 14, 19], ESGq limits its computation to smaller and less complex graphs, the Event Splicing Graphs. An event splicing graph is a splicing graph which encodes only the exons and splice junctions involved in an alternative splicing event. Differently from splicing graphs commonly used in the literature, that represent all known transcripts of a gene, and diferently from pantranscriptomes, where entire gene loci and intergenic regions are represented, an Event Splicing Graph encodes only the portions of the two transcripts involved in an event. By using this simpler representation, ESGq is able to achieve great eficiency without sacrificing its accuracy.
ESGq consists of three steps (also depicted in Figure 1): • to model an exon skipping event (SE), ESGq needs to take into account three exons and this implies that the corresponding event splicing graph is composed of three nodes 1, 2, 3. In detail, 2 represents the exon that is spliced out during the event. The canonical isoform is represented by the path involving all three nodes ( = 1 → 2 → 3) while the alternative isoform consists in the skip of 2 ( = 1 → 3); • to model an alternative acceptor site event (A3),
ESGq needs to take into account three exons and accordingly three nodes 1, 2, 3. In detail 1 represents the upstream exon, 2 the canonical downstream exon, and 3 the downstream exon with the alternative acceptor splice site. The canonical isoform is represented by the path involving the shared upstream exon and the canonical downstream exon ( = 1 → 2) whereas the alternative isoform changes the downstream exon with the alternative one ( = 1 → 3); • to model an alternative donor site event (A5),
ESGq needs to take into account three exons and accordingly three nodes 1, 2, 3. In detail, 1 represents the canonical upstream exon, 2 the canonical downstream exon, and 3 the upstream exon with the alternative donor splicing site. The canonical isoform is represented by the path involving the canonical upstream exon and the shared downstream exon ( = 1 → 2) whereas the alternative isoform changes the up
stream exon with the alternative one ( = 3 → 2); • to model an intron retention event (RI), ESGq needs to take into account three exons. However, this case is harder than the previous ones: the three nodes 1, 2, 3 of the graph do not closely correspond to three exons, but one of them (3) correspond to a portion of it. In detail, 1 and 2 represent the two upstream and downstream exons whereas 3 represents the retained intron (i.e., the internal portion of the exon linking the upstream and downstream canonical exons). The canonical isoform is then represented by the path involving the upstream exon and the downstream exon ( = 1 → 2) whereas the alternative isoform includes the retained intron in the path ( = 1 → 3 → 2).
contributes to an event splicing graph, but, from a more practical point of view, ESGq builds a single graph with multiple connected components, one per event.
In the second step, ESGq indexes the graph constructed in the previous step. Since the goal is to align the input RNA-Seq reads to the graph using the Girafe aligner [ 23], ESGq employs the VG toolkit [24] to build the gBWT (graph Burrows–Wheeler Transform) index [25]. Each input replicate is then independently aligned to the graph using Girafe. We note that, since by default vg breaks each node longer than 32bp into smaller nodes (of length ≤ 32), in order to keep the association between the nodes in the input graph (that is the one built by ESGq) and the nodes in the indexed version, ESGq directly breaks the nodes while building the event splicing graphs and links them accordingly to maintain the same two paths, i.e., isoforms. Due to this, in an event splicing graph we can observe two kind of edges: edges linking the smaller (≤ 32) nodes internal to an exon and edges that represent the real splice junction of interest for the AS event. Although ESGq diferentiates between these two kinds of edge, conceptually only the edges representing a splice junction are used by ESGq. For this reason, we decided to omit these edges from Figure 2.
In the third and last step, ESGq computes the value of the events w.r.t. each replicate and then summarize these values by comparing the two conditions and computing a ∆ value per event. This value represent the diferential expression of each event across the two input conditions.
To do so, ESGq analyses the graph alignments computed in the previous step and assign a weight to each edge that represents a splice junction. Since each read is aligned to a path of the graph, computing this weight is straightforward as increasing a counter per edge. Indeed, a read can be aligned to a single node of the graph, hence without using any edge, or to a sequence of nodes, hence using one or more edges. In such a case, ESGq checks every edge used by the alignment and, if an edge is a junction edge, it increases its weight by 1. In other words, since each junction edge represent a splice junction, its weight represents the number of reads that have been spliced aligned over it.
Finally, ESGq uses these weights to compute the value of each AS event following its classical formulation, i.e., the proportion of reads supporting the standard isoform over the reads supporting both isoforms [26]. A3
G TC TA G TC TA 1 E1 E1 1 E1 E1
W1 I1
W2
I2 2 W3 E2 I3 W1
W2 I1
I2 3 E3 E3 2 E2 3 E3 A5
G TC TA G TC TA 1 E1 1 E1 3 E3
W2 3 I1 E
W1 W2 I2
I1
W3 2 E2 2 E2 E2
Diferently from other approaches, which rely on both spliced and not spliced reads, ESGq computation is based only on spliced reads counts, hence the support of an isoform is approximated using only these values and does not take into account its full coverage. We believe that this is a good approximation of the correct value and this is also confirmed by our experimental evaluation. calculation can be summarized as follows (we also refer the reader to Figure 2):
Python and the code is freely available at
https://github.com/AlgoLab/ESGq. We note that
the most computationally intensive steps of the pipeline
(i.e., graph indexing and read-to-graph alignment) rely
on the VG toolkit [24], that is implemented in C++. We
assessed ESGq eficacy and eficiency on a real dataset
of RNA-Seq reads (SRA BioProject ID: PRJNA718442)
• = 1+21+22+23 considers the mean of the tbheatwtceoemneasgferionmg aanredcednitfersetundtiyal[2g8e]noenetxhperceossriroenlatiinon
weights 1, 2 of the canonical isoform and the Drosophila Melonogaster. More precisely, the study
weight 3 of the alternative isoform (in a similar conducted a genome-wide diferential expression
fashion to [27]); analysis at two diferent time points: day 1 and day 60
• 3 = 1+12 and 5 = 1+12 consider lfyes. The dataset consists of three replicates for two
the weight 1 of the canonical isoform and the conditions (the two time points), for a total of 6 Illumina
weight 2 of the alternative isoform Hiseq samples. All samples are paired-end and consist of
• = 1+12+2 3 considers the weight 1 151bp-long reads (see Table 1 for more details).
of the canonical isoform and the mean of the Diferently from the aforementioned study, where the
weights 2, 3 of the alternative isoform focus was the analysis of diferential gene expression, in
this work, we analyze the same dataset from the
perspecStarting from these values (one per event per replicate), tive of diferential quantification of alternative splicing
which summarize the event quantification for each input events. To do so, we applied ESGq and two other
statereplicate, ESGq computes the diferential quantification of-the-art approaches for the diferential quantification
across the two input conditions (∆ ) as the diference be- of alternative splicing events across multiple conditions:
tween the absolute value of the means in the two condi- rMATS [
(a) ESGq vs rMATS (b) ESGq vs SUPPA2 (c) rMATS vs SUPPA2
Condition Replicate n.Pairs Size (GB) methodologies based on diferent filtering criteria [
SRR14101759 26 658 610 19.2 In our analysis we included all those events correctly Day 1 SRR14101760 25 474 257 18.4 quantified by all three tools and considered statistically SRR14101761 28 339 185 22 significant by both rMATS and SUPPA2, i.e., events with SRR14101762 24 985 317 18 p value ≤ 0.05. A total of 933 events resulted from this Day 60 SRR14101763 25 569 084 18.6 iflter: 374 exon skipping (40%), 190 alternative 3’ (20%), SRR14101764 24 605 265 17.8 154 alternative 5’ (17%), and 215 intron retention (23%).
We note that we also tried to include Whippet [15] in Table 1 our evaluation but, due to a diferent AS event representaReal dataset used in our experimental evaluation (SRA Bio- tion that cannot be easily compared to the representation Project ID: PRJNA718442). Table reports the number of reads given by the other tools, we ended up omitting it from and the size in GigaByte of each paired-end replicate. the analysis.
Figure 3 and Table 3 report the results of our analysis. All tools achieved comparable results. Remarkto compare three methodologies based on completely ably, ESGq and rMATS achieved a very good correlation, diferent frameworks: graph-based alignment, reference- with a Pearson correlation coeficient equal to 0.918 (Figbased alignment, and transcript-based quasi-mapping. ure 3a). Although both approaches are based on read rMATS was run starting from the alignments produced alignment, they show two main methodological diferby STAR aligner [30] (version 2.7.10b) whereas SUPPA2 ences that slightly afect their results. Firstly, rMATS uses was run starting from Salmon transcript quantification read counts coming from read alignment to a reference (version 1.10.1). All tools were run using their default genome whereas ESGq uses (spliced) read counts comparameters and 16 threads (when possible). In our analy- ing from alignment to event splicing graphs, that are ses, we considered the reference, gene annotation, and a reduced representation. Secondly, the quantification transcripts provided by FlyBase [31], release 6.52. step implemented in ESGq is quite simplistic and not
We compared the ∆ values reported by the three elaborate as the probabilistic framework implemented in considered tools. ESGq reported 3 276, rMATS reported rMATS. However, none of these two diferences (that are 3 699, and SUPPA2 reported 1 619 AS events correctly a wanted restriction and a current - undesired - limitaquantified. We note that both ESGq and SUPPA2 start tion) seems to afect the results of ESGq. On the other from a list of events and quantify them: each time an hand, SUPPA2 resulted less correlated to the two other event can not be correctly quantified (due to, for instance, approaches, showing a Pearson correlation coeficient no coverage support), they assign ∆ = to that ranging from 0.766 to 0.813 (Figure 3b and 3c). This event. In that case, the event is not considered as cor- was somewhat expected since it is based on transcript rectly quantified, thus excluded from the analysis. The quantification, and not on spliced read alignment. As huge diference in the number of reported events proves proven in the literature, such a diference is expected the complexity of detecting AS events from RNA-Seq since current transcript annotation models may result reads and highlights the inconsistency between diferent inaccurate [15].
Event type Tool1 Tool2 Pearson methodologies are afected by read length, starting from
ESGq rMATS 0.952 the 151-bp paired-end dataset, we manually trimmed the SE ESGq SUPPA2 0.808 input reads to 51bp and 101bp using seqtk and created rMATS SUPPA2 0.836 two additional datasets. Table 3 reports the results of this ESGq rMATS 0.868 analysis. Even with very short reads (i.e., 51bp reads), the A3 ESGq SUPPA2 0.786 three tools achieved the same correlation (with a very rMATS SUPPA2 0.862 marginal diference of < 0.015).
Finally, to evaluate how much paired-end information
ESGq rMATS 0.920 may improve the accuracy of the tools, we merged the A5 rEMSAGTqS SSUUPPPPAA22 00..676068 two pairs of each replicate into a single sample, in order to simulate a single-end dataset. Surprisingly, there is small
ESGq rMATS 0.859 to none diference between the results on paired-end and RI ESGq SUPPA2 0.677 single-end dataset, highlighting the robustness of the rMATS SUPPA2 0.760 considered approaches. For instance, the Pearson corTable 2 relation coeficient between ESGq and rMATS decreased Pearson correlation coeficients between ESGq, rMATS, and by a very marginal 0.009 whereas correlation between SUPPA2 (over the transcript quantification of Salmon ran with SUPPA2 and the other approaches decreased by 0.041 = 31) broken down by event type on the 151bp paired-end (w.r.t. ESGq) and 0.027 (w.r.t. rMATS). dataset. The experimental evaluation has been implemented as a Snakemake workflow [ 32], thus it is fully reproducible and easily replicable. Scripts and instructions are
Table 2 reports the correlation between the three con- available at https://github.com/AlgoLab/ESGq. All the sidered tools broken down by event type. Surprisingly experiments were performed on a 64bit Linux (Kernel there is no clear trend that can be observed. ESGq and 5.15.0) system equipped with two 16-core AMD EPYC rMATS exhibit the highest correlation on exon skipping 7301 2.2GHz processors and 128GB of RAM. events and the lowest on intron retentions. ESGq and SUPPA2, instead, show higher correlation on exon skippings and lower correlation on alternative donor events. 4. Conclusions Finally, rMATS and SUPPA2 exhibit higher correlation on alternative acceptor site and lower correlation on In this paper we introduced ESGq, a novel graph-based alternative donor events. These results are somewhat approach for the AS event quantification across two conunexpected and require a further investigation. ditions. Diferently from state-of-the-art tools, ESGq is
ESGq also resulted very computationally eficient, com- based on read alignment against local graph structures, pleting the analysis in half an hour requiring 1GB of RAM. introduced here as event splicing graphs, that represent Similarly, SUPPA2 ran in less than 10 minutes and used AS events precisely represented in a given gene annota1.5GB of RAM. On the contrary, rMATS resulted the most tion. An extensive exploratory analysis on real RNA-Seq expensive approach, requiring more than 5 hours and dataset showed that ESGq is able to achieve comparable 8GB of RAM. The most expensive step is read alignment results with respect to other approaches based on the with STAR, that required from half an hour to two hours alignment of reads to the reference genome, while being per sample. By pairing simple and precise graph repre- 10x faster. sentation of well localized loci of the genome (i.e., the Future works will be devoted to improving the statistievent splicing graphs) with fast and accurate read align- cal framework behind the and ∆ computation and to ment, ESGq is able to achieve results comparable to the extending the experimental evaluation by assessing the other alignment-based approach, while being 10x faster. actual accuracy of ESGq (e.g., by computing performance
Since SUPPA2 is based on the -mer based quasi- metrics on simulated and RT-PCR validated events). An mapping of Salmon, we also analyzed how -mer size interesting future direction consists in extending the noafects its results. For this reason, we ran Salmon (and, tion of event splicing graphs to include information on consequently, SUPPA2) two additional times with ∈ known genetic variations (SNPs and indels) in order to {13, 21} (we note that 31 is the default value used in improve the quality of read alignment, as proven in a the previous results). As shown in Table 3, the results recent work on pantrascriptomes [19]. Moreover, the of SUPPA2 seems to be unafected by the choice of the detection and quantification of novel AS events from parameter. Indeed, the correlation between SUPPA2 pantrancriptomes remain another interesting open prob(ran with diferent value) and other tools changes lem. marginally.
Moreover, to evaluate if the results of the considered
ESGq rMATS ESGq rMATS ESGq rMATS ESGq rMATS rMATS SUPPA2 (k13) SUPPA2 (k21) SUPPA2 (k31) SUPPA2 (k13) SUPPA2 (k21) SUPPA2 (k31)
rMATS SUPPA2 (k13) SUPPA2 (k21) SUPPA2 (k31) SUPPA2 (k13) SUPPA2 (k21) SUPPA2 (k31)
rMATS SUPPA2 (k13) SUPPA2 (k21) SUPPA2 (k31) SUPPA2 (k13) SUPPA2 (k21) SUPPA2 (k31)
rMATS SUPPA2 (k13) SUPPA2 (k21) SUPPA2 (k31) SUPPA2 (k13) SUPPA2 (k21) SUPPA2 (k31)
[13] A. Kahles, C. S. Ong, Y. Zhong, G. Rätsch, Spladder: variation in the reference, Nature biotechnology 36 identification, quantification and testing of alterna- (2018) 875–879. tive splicing events from rna-seq data, Bioinformat- [25] J. Sirén, E. Garrison, A. M. Novak, B. Paten, ics 32 (2016) 1840–1847. R. Durbin, Haplotype-aware graph indexes, Bioin[14] L. Denti, R. Rizzi, S. Beretta, G. D. Vedova, M. Pre- formatics 36 (2020) 400–407. vitali, P. Bonizzoni, Asgal: aligning rna-seq data to [26] N. L. Barbosa-Morais, M. Irimia, Q. Pan, H. Y. Xiong, a splicing graph to detect novel alternative splicing S. Gueroussov, L. J. Lee, V. Slobodeniuc, C. Kutter, events, BMC bioinformatics 19 (2018) 1–21. S. Watt, R. Colak, et al., The evolutionary landscape [15] T. Sterne-Weiler, R. J. Weatheritt, A. J. Best, K. C. of alternative splicing in vertebrate species, Science Ha, B. J. Blencowe, Eficient and accurate quantita- 338 (2012) 1587–1593. tive profiling of alternative splicing patterns of any [27] K.-T. Lin, A. R. Krainer, Psi-sigma: a comprehencomplexity on a laptop, Molecular Cell 72 (2018) sive splicing-detection method for short-read and 187–200.e6. long-read rna-seq analysis, Bioinformatics 35 (2019) [16] J. Vaquero-Garcia, A. Barrera, M. R. Gazzara, 5048–5054.
J. Gonzalez-Vallinas, N. F. Lahens, J. B. Hogenesch, [28] M. Bajgiran, A. Azlan, S. Shamsuddin, G. Azzam, K. W. Lynch, Y. Barash, A new view of transcrip- M. A. Halim, Data on rna-seq analysis of drosophila tome complexity and regulation through the lens melanogaster during ageing, Data in brief 38 (2021) of local splicing variations, elife 5 (2016) e11752. 107413. [17] Y. I. Li, D. A. Knowles, J. Humphrey, A. N. Barbeira, [29] R. Patro, G. Duggal, M. I. Love, R. A. Irizarry, S. P. Dickinson, H. K. Im, J. K. Pritchard, Annotation- C. Kingsford, Salmon provides fast and bias-aware free quantification of rna splicing using leafcutter, quantification of transcript expression, Nature Nature genetics 50 (2018) 151–158. methods 14 (2017) 417–419. [18] W.-W. Liao, M. Asri, J. Ebler, D. Doerr, M. Haukness, [30] A. Dobin, C. A. Davis, F. Schlesinger, J. Drenkow, G. Hickey, S. Lu, J. K. Lucas, J. Monlong, H. J. Abel, C. Zaleski, S. Jha, P. Batut, M. Chaisson, T. R. Ginet al., A draft human pangenome reference, Nature geras, Star: ultrafast universal rna-seq aligner, 617 (2023) 312–324. Bioinformatics 29 (2013) 15–21. [19] J. A. Sibbesen, J. M. Eizenga, A. M. Novak, J. Sirén, [31] J. Thurmond, J. L. Goodman, V. B. Strelets, H. AtX. Chang, E. Garrison, B. Paten, Haplotype- trill, L. S. Gramates, S. J. Marygold, B. B. Matthews, aware pantranscriptome analyses using spliced G. Millburn, G. Antonazzo, V. Trovisco, et al., Flypangenome graphs, Nature Methods (2023) 1–9. base 2.0: the next generation, Nucleic acids research [20] S. Beretta, P. Bonizzoni, L. Denti, M. Previtali, 47 (2019) D759–D765.
R. Rizzi, Mapping rna-seq data to a transcript graph [32] J. Köster, S. Rahmann, Snakemake—a scalable via approximate pattern matching to a hypertext, in: bioinformatics workflow engine, Bioinformatics Algorithms for Computational Biology: 4th Inter- 28 (2012) 2520–2522. national Conference, AlCoB 2017, Aveiro, Portugal, June 5-6, 2017, Proceedings 4, Springer, 2017, pp.
49–61. [21] M. F. Rogers, J. Thomas, A. S. Reddy, A. Ben-Hur,
Splicegrapher: detecting patterns of alternative splicing from rna-seq data in the context of gene models and est data, Genome biology 13 (2012) 1–17. [22] S. Beretta, P. Bonizzoni, G. D. Vedova, Y. Pirola,
R. Rizzi, Modeling alternative splicing variants from rna-seq data with isoform graphs, Journal of
Computational Biology 21 (2014) 16–40. [23] J. Sirén, J. Monlong, X. Chang, A. M. Novak, J. M.
Eizenga, C. Markello, J. A. Sibbesen, G. Hickey, P.C. Chang, A. Carroll, et al., Pangenomics enables genotyping of known structural variants in 5202 diverse genomes, Science 374 (2021) abg8871. [24] E. Garrison, J. Sirén, A. M. Novak, G. Hickey,
J. M. Eizenga, E. T. Dawson, W. Jones, S. Garg, C. Markello, M. F. Lin, et al., Variation graph toolkit improves read mapping by representing genetic