Upstream analysis of alternative splicing: a review of computational approaches to predict context-dependent splicing factors

Abstract

Alternative splicing (AS) has shown to play a pivotal role in the development of diseases, including cancer. Specifically, all the hallmarks of cancer (angiogenesis, cell immortality, avoiding immune system response, etc.) are found to have a counterpart in aberrant splicing of key genes. Identifying the context-specific regulators of splicing provides valuable information to find new biomarkers, as well as to define alternative therapeutic strategies. The computational models to identify these regulators are not trivial and require three conceptual steps: the detection of AS events, the identification of splicing factors that potentially regulate these events and the contextualization of these pieces of information for a specific experiment. In this work, we review the different algorithmic methodologies developed for each of these tasks. Main weaknesses and strengths of the different steps of the pipeline are discussed. Finally, a case study is detailed to help the reader be aware of the potential and limitations of this computational approach.

splicing factors, alternative splicing, RNA-binding proteins, bioinformatics

Introduction

Alternative splicing (AS) is the mechanism by which a single pre-mRNA molecule can lead to different mature mRNA molecules, called isoforms or transcripts. In this process, exons can be either included or excluded, shortened or lengthened and skipped or retained. The transcriptome is the complete set of mRNA isoforms in an organism. The phenomenon of AS was first described by Berget et al. [1], where it was shown that one adenovirus produced several transcripts during its infectious cycle.

The number of discovered isoforms increases as the study of an organism improves. In humans, around 95% of multi-exonic genes present AS events in diverse conditions [2, 3]. The paradigm ‘one gene-one protein’ has switched to the present situation in which most genes encode several proteins because of AS [4].

The functions altered by AS can be different: biomass generation, induction of angiogenesis, loss of genomic stability or deterioration of the immune system among others [5]. Besides, the influence of AS on neoplasms and other diseases is well known [6, 7]. In fact, studies suggest that approximately one-third of all disease-causing mutations modify splicing [8]. The regulation of AS has become a therapeutic strategy, and it is also revealing new therapeutic targets [6]. It has also been shown that all the cancer paradigms [9] have their counterpart in aberrant splicing [5].

The mechanism of splicing involves a complex biological machinery with several elements, such as RNA-binding proteins (RBPs), other trans-acting factors or specific sequences in the mRNA, which are the target signals of RBPs. Even epigenetics has been shown to play a role in AS [10–12]. The detection of AS itself is not trivial and requires specific algorithms and software to identify and label AS events.

Multiple approaches have been developed to understand the link between AS events and splicing regulatory elements (SREs) in different diseases. Several works analyzed brain-specific splicing factors (SFs) such as NOVA1 and NOVA2 [13–15]. Besides, one of the first global analysis for the identification of cancer-associated AS events and regulators was performed by Danan-Gotthold et al. in 2015 [16]. As these influential works, numerous strategies have emerged to decipher the splicing mechanism associated with diverse pathologies.

The scope of this review is the description of the computational approaches to detect splicing events and to predict their regulatory elements, i.e. upstream analysis of AS. We will not discuss the functional effects of AS but only its detection and potential regulation. Finally, a case study is detailed to help the reader be aware of the potential and limitations of these computational approaches.

Overview of splicing

Splicing is a posttranscriptional process in which nucleotide sequences, called introns, are removed from the pre-mRNA. The resulting product is a mature mRNA molecule that includes 5′ and 3′ untranslated regions (5′/3′ UTRs) and coding regions (exons) joined together to form a single mRNA strand. The splicing process is orchestrated by the spliceosome, a complex machinery made up of different subunits known as small nuclear ribonucleoproteins (snRNPs) and other protein complexes. snRNPs are noncoding and non-polyadenylated RNA–protein complexes that carry out their functions in the nucleoplasm [17]. A deeper explanation of the splicing mechanism has been included in the Supplementary Material (Section 1: splicing mechanism).

Splicing events in eukaryote cells can be classified into two main groups: constitutive splicing events, which always occur and give rise to the same isoforms independently of the tissue or pathological situation; and AS events, which lead to different isoforms. AS events have been divided into several canonical classes, as shown in Supplementary Figure S1.

Regulatory elements of AS

The mechanisms that control and regulate AS are still subject to active research [3, 11, 18–20]. Here, we focus on two key elements in the regulation of AS: cis-acting RNA elements and trans-acting factors. A scheme of the elements that take part in the AS process is shown in Figure 1A.

Figure 1

(A) Overview of AS process. An example of a cassette exon with its regulating elements is shown. The branch point (BP), and the polypyrimidine tract (PY) are also represented. (B) General pipeline to detect SFs: (1) identification of AS events in a specific condition; (2) identification of binding motifs of RBPs; and (3) prediction of splicing regulatory factors. This task is, in turn, divided into: (3.1) scanning SF’s motifs in splicing regions and (3.2) motif enrichment analysis using PSI, expression levels of RBPs and other sources of information (driver mutations, CNVs, coexpression networks, etc.). Each of the boxes corresponds to a section in the main text.

Cis-acting RNA elements

Cis-elements are RNA sequences (or motifs) in the pre-mRNA that allow the recognition of specific exonic/intronic regions by the spliceosome. Mainly, they comprise the canonical splicing signals and the SREs.

Splicing signals are essential sequences (the 5′ splice site, the 3′ splice site and the adenine branch point) for recognition by the spliceosome. SREs are divided into exonic splicing enhancers, exonic splicing silencers (ESSs), intronic splicing enhancers and intronic splicing silencers. The activity of SREs depends on the recruitment of molecules, which impact positive or negatively in the splicing reaction steps [21]. One common example is the polypyrimidine tract-binding protein (PTB), which causes exon skipping after binding to an ESS by avoiding the formation of the exon definition complex [22]. The effect of cis-acting RNA elements can be altered because of factors such as decoy splice-sites [17] or the surrounding context [19].

Trans-acting splicing regulators

Trans-acting splicing regulators are analytes—usually proteins—that, by interacting with the mRNA, modulate its AS. Most of these trans-acting factors are RBPs.

RBPs are proteins that bind to single- or double-stranded RNA and play key roles in posttranscriptional gene regulation, such as regulation of AS, mRNA stabilization, mRNA location, polyadenylation or translation [23]. They usually have modular designs and consist of various repeats of just a few basic RNA-binding domains, which have, in turn, different strategies to RNA binding. The ability to selectively recognize and bind mRNAs is crucial for the correct functionality of RBPs [3, 11, 18–20]. Most frequent RNA-binding domains, namely, RNA-recognition motif, heterogeneous nuclear ribonucleoprotein K-homology, double-stranded RNA binding domain and Zinc fingers—are described in the Section 1 of the Supplementary Material. We refer the reader to [19, 20, 24–31] for further details.

These domains are the main players of RBPs-mRNA interactions. In turn, these interactions have been extensively studied and compiled in several databases [32–36], but many of them are still unknown. The ATtRACT database [37], which contains curated and validated data from the main databases of RBPs-mRNA interactions (CisBP-RNA [32], SpliceAid-F [38] and RBPDB [39]), collects the binding information of 370 RBPs, which represent about 30% of the ∼1300 RBPs that have currently been discovered [31, 40].

RBPs that participate in the AS regulation are called SFs. These RBPs mainly include serine- or arginine-rich (SR) proteins and heterogeneous nuclear ribonucleoproteins (hnRNPs). The binding profile and functions of a number of SFs have been previously studied [41–44] and reviewed [19]. In general, SR proteins—such as SRSF1 or SRSF2—are considered positive splicing regulators, as they promote exon inclusion [45–47]. In contrast, hnRNPs—such as hnRNP A1, hnRNP A2 or hnRNP B1—seem to have the opposite effect, as they avoid the formation of the splicing machinery.

In addition to the spliceosome, there are other processes that play an important role in the regulation of the AS, such as DNA methylation [12], chromatin status [10], histone modifications [11], phosphorylation of the corresponding RBPs [48] or the secondary structure of the pre-mRNA [49]. In many cases, changes in either of these processes impact on AS.

In the next sections, we discuss the computational approaches to identify AS events and predict their context-dependent regulators. Although some individual parts of these tasks have already been covered by other reviews (Table 1), this work tries to provide a broader view of algorithms developed to unveil the complex regulation of AS.

Table 1

Description of previous reviews of experimental and computational methods related to the upstream analysis of the splicing process

Review	Description	Algorithms reviewed	Reference
A survey of software for genome-wide discovery of differential splicing in RNA-seq data	A review of the software available for analysis of RNA-seq data for differential splicing	Identification of AS sites and events Cuffdiff 2, MISO, DEXSeq, DSGseq, MATS, DiffSplice, Splicing compass, AltAnalyze	Hooper [50]
Advances in the characterization of RNA-binding proteins	Experimental and computational methods for detection of protein–RNA interactions	Experimental methods for detection of protein–RNA interactions	Marchese et al. [34]
		(Protein-centric) RIP, HiTS-CLIP, PAR-CLIP, iCLIP, eCLIP, RNA-compete, SEQRS, RBNS, RNA-MaP, HiTS-RAP, MITOMI.
		(RNA-centric) TRAP/RAT, RaPID, RiboTrap, RNA-assisted chromatography, protein microarray, MS2-BioTRAP, ChIRP, CHART, RAP-MS, Interactome capture
		Motif discovery
		MEME, RBPmap, SeAMotE, RNAcontext
High-throughput characterization of protein–RNA interactions	Review of (a) experimental characterization of RBP-RNA interactions, (b) algorithms to predict RNA secondary structure and (c) motif finding tools	Motif discovery	Cook et al. [24]
		MEME, PhyloGibbs, cERMIT, DRIMUST, StructuRED, TEISER, RNAcontext, GraphProt, CMfinder, RNApromo, MEMERIS
Evaluating tools for TFBS prediction	Review and performance comparison of (a) de novo motif discovery tools and (b) transcription factors binding sites prediction tools	Motif discovery	Jayaram et al. [51]
		rGADEM, HOMER, ChIP- Munk, MEME-ChIP
		Scanning motifs
		Baycis, Cister, MCast, Comet, ClusterBuster, Matrix-Scan, Clover, FIMO, Patser, PossumSearch
Finding the target sitesof RNA-binding proteins	Comprehensive review of resources and methods to detect protein–RNA interactions. It focuses on the importance of the secondary structure of RNA	Resources for RBP binding sites	Li et al. [52]
		ARESITE, CisBP-RNA, CLIPz (no longer available), doRiNA, RBPDB, Rfam, UTRSite
		Motif discovery
		MatrixREDUCE, MEME, MEMERIS, REFINE, AMADEUS, Aptamotif, CMfinder, cERMIT, COVE, FIRE, RNAalifold, RNAcontext, RNApromo

Computational approaches to identify splicing and its regulatory elements

The algorithms will be presented in a conceptual sequential order: first, we discuss the algorithms to detect AS from high-throughput transcriptomic data [RNA sequencing (RNA-seq) and microarrays]. Next, we show different methods to discover the binding motifs of RBPs and, finally, we describe how to combine this information with experiment-specific data (expression of RBPs and relative usage of the exons in an event among others) to predict context-dependent regulators of AS. A summary of the pipeline followed can be found in Figure 1B.

Identification of as sites and events

An AS event is a local alteration of the splicing pattern on a gene, that in turn originates different isoforms. Some of these alterations occur more frequently and are called ‘canonical’ events. These canonical events are described in Supplementary Figure S1. In the case of canonical events, each event has two alternative configurations—the exon is either included or excluded, the 3′ extension can be included or excluded and so on. There are other events that can involve more complex patterns of AS for the same locus in the mRNA—e.g. three exons that can be skipped and are mutually exclusive.

Detecting the AS events is necessarily the first step in the identification of potential regulators: once the events with differential usage in different conditions are identified, using other computational methods it is possible to predict context-dependent regulators.

The task of identifying AS events has already been studied. Hooper [50] reviewed some tools that detect AS with RNA-seq data, but in the past 3 years, there has been a huge development of this family of algorithms. This section includes 33 methods to identify splicing events using either microarrays or RNA-seq data. Table 2 summarizes the reviewed algorithms. In the following paragraphs, we explain the criteria to include a method (rows of the Table 2) and their key characteristics (groups and columns of Table 2).

Table 2

Algorithms for the identification of as events

Algorithm family	General aspects	Algorithm	Operating system	E.	S.	V.	PSI	Information used for quantification	References
RNA-seq novel events	(+) Detect nonannotated events (−) Time-consuming and complexity of the algorithms	AltAnalyze	All	✓	✓	✓	×	Exons and junctions	[53]
		ASpli*	All	✓	✓	×	✓	Only Junctions	[54]
		CASH	All	✓	✓	×	∼	Exons and junctions	[55]
		DEXseq	All	×	✓	×	×	Only Exons	[56–58]
		DiffSplice	Linux	ASM	✓	✓	✓	Exons and junctions	[59]
		EventPointer	All	✓	✓	✓	✓	Exons and junctions	[60]
		Gess	All	CE	✓	×	✓	Only exons	[61]
		JuncBASE	All	✓	✓	×	×	Only junctions	[15]
		Leafcutter	All	×	✓	✓	✓	Only junctions	[62]
		MAJIQ+VOILA	Linux	✓	✓	✓	✓	Junction reads	[63]
		rMATS	Linux	✓	✓	×	✓	Exons and junctions	[64]
		SGSeq	All	✓	×	✓	✓	Exons and junctions	[65]
		SPLADDER	All	✓	✓	×	✓	Exons and junctions	[66, 67]
		SplicePie	Linux	∼	✓	×	✓	Exons and junctions	[68]
		SplicingTypesAnno	All	✓	✓	✓	✓	Exons and junctions	[69]
RNA-seq known events	(+) Better adapted to compare disparate experiments (+) Faster (−) Non-novel events	ASATP*	All	✓	✓	✓	×	Expression of isoforms involved in event	[70]
		ASprofile	Linux	✓	×	×	×	Expression of isoforms involved in event	[71]
		DSGseq	All	×	✓	×	×	Exons	[72]
		IMAS*	All	×	✓	✓	✓	Exons	[73]
		SpliceR	All	✓	×	×	×	Expression of isoforms involved in event	[74]
		SpliceSEQ	All	✓	✓	✓	×	Exons and junctions	[75]
		SpliceTrap	Linux	✓	✓	×	✓	Exons	[76]
		SplicingCompass	All	×	✓	✓	×	Exons and junctions	[77]
		SplicingExpress	Linux	✓	✓	✓	×	Expression of isoforms involved in event	[78]
		SUPPA	All	✓	✓	✓	✓	Expression of isoforms involved in event	[79]
		Vast-Tools	Linux	✓	✓	✓	✓	Exons and junctions	[80]
Arrays	(+) Good performance (−) Non-novel events	AltAnalyze	All	✓	✓	✓	×	Only Exons	[53]
		EventPointer	All	✓	✓	✓	✓	Exons and junctions	[60]
		ExonPointer	All	CE	✓	✓	×	Exons and junctions	[81]
		IGems	All	×	✓	✓	×	Only Exons	[82]
		MADS+	All	✓	✓	✓	×	Exons and junctions	[83]
		RASA	NA	×	✓	×	×	Exons and junctions	[84]
		TAC 4.0	Windows	✓	✓	✓	×	Exons and junctions	[85]

Notes: *There is not a peer-reviewed reference for this algorithm. E: event classification; S: this method provides statistics; V: visualization; PSI: whether the PSI is returned; CE: cassette exon; ASM: alternative splicing module (any type of event without labeling the canonical ones). It is divided into three groups: algorithms that use RNA-seq to discover novel and non-novel events, and microarray-based algorithms. Other characteristics, such as algorithm’s input data and some comments of each algorithm can be found in the Supplementary Table S1.

Criteria to include a method in the table. We focus on tools that detect AS events. We do not consider pipelines that quantify the expression of (novel) transcripts, such as Cufflinks [86, 87], MISO [88], SpliceGrapher [89] or Stringtie [90]. We do not either consider Nanopore or PacBio, as they are not suitable to pinpoint splicing events but the whole sequences of transcripts [91]. We do include methods that, taking as input isoform concentrations and structure, predicted by the previous or other algorithms, detect the presence of splicing events.

The reason why we only examine algorithms that detect AS events is that the transcriptome reconstruction is a problem much more difficult to solve. In fact, it was shown to be an NP-hard problem [92]. Different heuristics have been proposed, but they are far from perfect. Steijger et al. [93] stated that the best-performing methods have precision and recall around 40–50% at the transcript level for simulated data. This means that <50% of the predicted transcripts are correct and that <50% of the transcripts are recovered. These results worsen for complex genes with many transcripts. The same reference states that precision and recall rise to 80–90% if the analysis is performed at the exon level. Once the exons and the junctions that link them are known (i.e. given the splicing graph), the identification of AS events is straightforward. These facts make it more sensitive and reliable to focus on events than on transcripts to identify the potential regulators of AS. On the other hand, as the SFs bind to specific regions of the pre-mRNA, even if the isoforms were used as input to the algorithms, it would be necessary to perform the analysis at the event level (as many algorithms do).

We include some methods that use arrays. Our experience in the detection of events using microarrays and RNA-seq is that top results using both technologies show strong coherence between them [94]. RNA-seq, of course, has an edge on its ability to detect novel events, but the required computing resources using microarrays are much smaller. Therefore, we have decided to include microarray’s methods that can be applied to junction arrays. This filter also helps to remove methods that are no longer maintained.

Key characteristics.Table 2 is split into three groups: methods based on RNA-seq that discover novel events, methods based on RNA-seq that do not discover novel events and methods based on arrays. The boundary between RNA-seq methods that discover novel events and methods that do not is blurred. Most algorithms that detect annotated events use as input the transcript structure (GTF file) and the estimated transcript expression. If this information is generated by an isoform deconvolution software such as Cufflinks, they can also be used to detect novel events, as Cufflinks (and other methods) predict novel transcripts given the RNA-seq data. In this case, these methods would be unraveling an non-deterministic polynomial-time hard combinatorial problem (isoform deconvolution) to solve a much easier one (event detection and quantification).

The main methodologies proposed to quantify AS events are the percent spliced-in (PSI or Ψ) and the splicing index (SI). PSI [95] is an estimate of relative usage of each alternative path (specific configuration of exons and/or junctions) of an AS event. Estimates of PSI can be validated using a third technology such as PCR (either quantitative or standard). On the other hand, the SI states the relative signal/coverage of an exon or a junction compared with the whole gene.

The SI has two drawbacks. AS every exon or junction has its own SI, the coherence of the SI change of the different exons and junctions involved in an event is not taken into account. For example, for a cassette event in which the exon is skipped in a tumoral condition, the SI of the junction that skips the cassette exon will be positive, the SI of the junctions of the cassette exon and the cassette exon itself will be negative. As the SI is not summarized for the whole event, it corresponds to the researcher to state the coherence between these signals. On the other hand, SI is difficult to be validated using PCR because it would require to run a PCR for every exon and junction to measure the average value. In contrast, the PSI value can easily be validated using PCR. Finally, SI may show spurious changes even for constitutive exons. Algorithms that return the PSI are therefore preferred. In Figure 2 it can be seen both PSI and SI calculations for a cassette event.

Figure 2

(A) Toy example of an exon cassette with differential splicing across two conditions (normal and tumor). Coverage of exons and junctions in both conditions are included. (B) Computation of SI and PSI for a toy example. SI is computed using the log ratio of the coverage of each exon/junction with the average coverage of the whole gene-hear simplistically consider as the average of the coverage of the exons and junctions of the gene. On the other hand, PSI considers the ratio of the mean coverages of the exons and junctions that include the cassette exons (J1, J2 and E2) and the sum of the coverages of both isoforms.

Some of RNA-seq methods use only the exon or only the junction reads to quantify the splicing events. Either of them are theoretically inferior to integrating both sources of information. Junction reads tend to be more scarce and more difficult to map than exon reads [96]. Methods based on junction reads are more sensitive to the characteristics of the aligner than methods that integrate them with exon reads [96]. On the other hand, exon reads alone can miss changes in isoforms that correspond to the less expressed path. The coverage of the junction that skips a cassette exon would be especially informative to state the change if its isoform is weakly expressed. Methods that exploit both sources of information are preferred.

Discussion. Mats—predecessor of rMats [64]—and DEXseq [58] were the first algorithms developed with this purpose. Both of them are actively maintained and several improvements, and new functionalities have been included in them. For example, the ability to detect novel events was included in rMats in late 2016, and DEXseq has recently improved its underlying statistical analysis. The statistics related with these two methods are briefly described in the additional material. rMats is based on the PSI, and DEXseq performs a statistical analysis indirectly based on SI.

Some published methods show a comparison with other algorithms. For example, Spladder [67], rMats [64], SpliceGrapher [89] and JuncBase [15] were compared by the developers of Spladder. Using simulated data, the number of detected events using JuncBase or SpliceGrapher is larger than using rMats. On the contrary, rMats shows better true discovery rates than JuncBase when top-ranked events were selected. Spladder outperforms all the others in any respect in this comparison.

SUPPA developers compared it against MISO and MATS. In this study, MATS (an event-based method) was shown to outperform MISO (a transcript deconvolution method) and FIMMO to be the algorithm with best performance [79]. SGSeq developers compare it against methods based on the deconvolution of the transcriptome. The conclusion of both references is that PSI estimation is more reliable when using event-based methods than when using methods based on transcriptome deconvolution. These conclusions must be taken with prudence, as there can be confirmation bias.

Some algorithms do not discover novel events. Despite this disadvantage, these algorithms can be advisable in some situations. First, they can be better adapted to compare disparate experiments as long as the same reference transcriptome is used across such experiments. Second, these methods can provide results in much shorter time, using a fast isoform-quantification algorithm such as Kallisto [97] or Salmon [98]. It is up to the user to decide if the discovery of novel events is worth the additional burden of time and storage or the difficulties in performing meta-analyses.

Table 2 can be used as a guide to select the proper AS detector. If the analysis requires the detection of novel events, only the algorithms of the firs group can be used. Among them, those that use information of exons and junctions and provide the PSI and event classifications should be preferred (i.e. rMats, EventPointer-SGSeq, SPLADDER or SplicingTypesAnno).

If novel events are not required, using algorithms in the second group—probably after quantifying the isoforms using Kalisto or Salmon—is preferred. Among them, several algorithms return the PSI and classify the corresponding events [i.e. SpliceTrap (69), SUPPA (29) or Vast-Tools(687)].

Finally, the number of algorithms using arrays is smaller. EventPointer, AltAnalyze and IGEMs have been recently deployed and showed their performance in real data. Only EventPointer returns the PSI.

Discovering binding motifs of SFs

There are two main approaches to identify pairs of proteins–RNAs that have affinity to bind together: using the output data of biological experiments that characterize protein–RNA interactions [such as cross-linking and immunoprecipitation sequencing (CLIP-seq), photoactivatable ribonucleoside-enhanced crosslinking and immunoprecipitation (PAR-CLIP), RNA immunoprecipitation sequencing (RIP-seq)] and analyzing protein and RNA structures to predict potential binding sites. The first method combines experimental data with computational algorithms (Table 3), whereas the second is purely computational.

Table 3

Computational methods aimed at discovering motifs

Algorithm subtype	Main algorithms	Algorithm with best performance (reference)
Based on RNA primary structure	MEME [101], cERMIT [102], phyloGibbs [102], GLAM2 [103], HOMER [104], ChIP- Munk [105], DREME [106], rGADEM [107], MEME-ChIP [108], DRIMUST [109], RBPmap [110], SeAMotE [111]	rGADEM (Jayaram et al. [51])
Based on RNA secondary structure	MEMERIS [112], RNApromo [113], StructRED [114], RNAcontext [115], CMfinder [115], TEISER [116], mCarts [117], GraphProt [118]	RNAcontext* and MatrixREDUCE (Kazan et al. [115])

Notes: *Note that the authors of RNAcontext algorithm are also the authors of its corresponding comparative review. These methods extract binding motifs using CLIP-, RIP- or CHIP-seq experiments.

Only two references lie on the second family of methods [99, 100]. The computational burden of these two methods make them non-suitable to be applied genome-wide but to check the interaction of specific pairs of RBP-RNAs. These methods are not used in any of the computational approaches to find context-dependent SFs (Table 5).

Table 5

Main approaches to find context-specific SFs

Reference	Title	Inputs	Ouputs	Platform	Novel events	Software	Type of AS events (algorithm)	Data access (samples)	RBPs	Description and comments
Danan-Gotthold et al. [16]	Identification of Recurrent Regulated Alternative Splicing Events across Human Solid Tumors	E, Ψ, RxE	1, 2	RNA-seq	Yes	No	Only cassettes	TCGA	RBFOX2, QKI, CELF2, MBNL1, MBNL2 and PTBP1	A large-scale study of AS in human solid tumors
Sebestyán et al. [149]	Large-Scale Analysis of Genome and Transcriptome Alterations in Multiple Tumors Unveils Novel Cancer-Relevant Splicing Networks	E, Ψ, RxE, O	1, 2, 3	RNA-seq	No	No	CE, IR, A3, A5, AF, AL, MX (SUPPA)	TCGA	1348 RBPs (104 with motifs, CISbp)	Analysis of widespread alterations in the expression of RBP genes, novel mutations and copy number variations in association with multiple AS changes in cancer drivers and oncogenic pathways
Aghamirzaie et al. [150]	CoSpliceNet: a framework for co-splicing network inference from transcriptomics data	E, Ψ*, RxE	1, 2, 3	RNA-seq and Arrays	No	Yes (open- source)	Isoform-specific analysis	GEO: GSE74692	Defined by the user	A tool for co-splicing network inference, which can be used to identify SFs and their candidate targets pre-mRNAs
Zhang et al. [151]	MYCN Controls an Alternative RNA Splicing Program in High-Risk Metastatic Neuroblastoma	E, Ψ*, RxE	1, 2	RNA-seq	No	No	Isoform-specific analysis	dbGaP: phs000868	RBFOX1, RBFOX3, CELF2, CELF6, PTBP1 and HNRNPA1	Analysis of transcription factors that regulate SF genes by analyzing splicing patterns
de Miguel et al. [148]	A Large-Scale Analysis of Alternative Splicing Reveals a Key Role of QKI in Lung Cancer	Ψ, RxG	1	Arrays (HJAY)	No	No	Only cassettes (ExonPointer)	Non-Public	Targets of QKI. Manually curated.	Analysis of AS of lung cancer and QKI
Correa et al. [152]	Functional Genomics Analyses of RNA-Binding Proteins Reveal the Splicing Regulator SNRPB as an Oncogenic Candidate in Glioblastoma	E, Ψ, O	1	RNA-seq and Arrays (Human U219)	No	No	CE, MX, IR, A3, A5, AF, AL	TCGA & SRA	1542 RBPs (CISbp)	A systematic study of RBPs in GBM
Sveen et al. [5]	Aberrant RNA Splicing in Cancer; Expression Changes and Driver Mutations of Splicing Factor Genes	E, Ψ, O	1	RNA-seq, Arrays	Yes	No	CE, IR, A3, A5, AF, AL, MX (SpliceSeq)	TCGA, MediSapiens	261 RBPs	Review of aberrant RNA splicing and its regulation in several cancer types
Brooks et al. [15]	Conservation of an RNA regulatory map between Drosophila and mammals	Ψ, RxE	1, 2	RNA-seq	Yes	No	CE, IR, A3, A5, AF, AL, MX (JuncBASE)	Non-Public	NOVA1, NOVA2	Analysis of the conservation of RNA regulatory elements of NOVA1 and NOVA2

Notes: The table headers display the following fields: Inputs (E: expression of RBPs; Ψ: percent spliced-in; Ψ*: isoform relative usage; RxG: RBP–gene relationships; RxE: RBP–event relationships; O: other genetic sources of information such as copy number variations, or mutations). Outputs (1: SFs that regulate a specific condition; 2: SFs that regulate a specific event; 3: interaction network of RBPs). Platform: platform for mRNA measuring. Type of events (CE: cassette event; IR: intron retention; A3: alternative 3′; A5: Alternative 5′; AF: alternative first; AL: alternative last; MX: mutual exclusive; C: complex event). RBPs covered in each work.

Discovering binding motifs using protein-centric experiments

Experimental methods that characterize protein–RNA interactions can be divided into protein-centric methods, which identify binding RNAs for a particular protein, and RNA-centric methods, which discover the proteins that interact with a specific RNA region. RNA-centric methods are not strictly high-throughput, as only a few proteins bind to a specific locus in the genome. Therefore, we focus the discussion on the protein-centric ones (see Section 2 of Supplementary Material for details).

Most protein-centric experimental methods rely on RNA immunoprecipitation, namely, a protein antigen is precipitated using a protein-specific antibody, followed by RNA identification using either microarrays or RNA-seq. RNA immunoprecipitation (RIP) and CLIP are the main methodologies [14, 119]. They differ on the protocol for immunoprecipitation.

In RIP techniques, a protein antigen is precipitated using a protein-specific antibody, and RNAs are identified using either microarray (RIP-chip) or RNA-seq (RIP-seq) [119]. The main limitation of RIP relies on the low resolution and background noise that causes the detection of nonspecific interactions.

Protein-centered methods were improved with ultraviolet cross-linking and denaturing techniques (CLIP) and can be measured by RNA-seq (high-throughput sequencing of RNA isolated by crosslinking immunoprecipitation (HiTS-CLIP) [14]). Figure 3A shows an example of the output data of a CLIP-seq experiment (BED file). More details about these experimental methods can be found in the review of Marchese et al. [34].

Figure 3

(A) Common output data (BED file) of RBP-RNA interaction experiments such as CLIP-seq. It consists of the genomic ranges in which RBP bind to, the sequence name, a binding score and the strand. (B) Input data of most motif discovery algorithms: a FASTA file of sequences of RNA. (C) Output of MEME-ChIP algorithm: logos of PWMs, the algorithm used to find each motif, E-values of the discovered PWM, similar known motifs, centered distribution of the motif in input sequences and other options to perform additional analyses.

Protein-centric techniques usually return a collection of RNA sequences attached or close to the binding sites of RBPs. The corresponding genomic regions for these sequences are long (50–500 nt) compared with the loci where the RBPs bind, typically a few nucleotides. There are several databases that compile information of these experiments, such as DoRiNA [120], CLIPdb [121] or POSTAR [122] (Table 1).

Identifying short recurring motifs makes it possible to predict RBP-mRNA potential binding sites without the need of a RIP or CLIP experiment. Algorithms to find recurring motifs in multiple, unaligned and long sequences are known as motif discovery algorithms. The input of these programs is a collection of sequences (FASTA files) given by an RIP or a CLIP experiment. The output is a set of motifs that appear recurrently in the given sequences. These motifs are usually represented as position weight matrices (PWMs). An example of the input and output of this family of algorithms is shown in Figure 3.

RNAcompete (with the companion database CisBP) performs a different approach: it uses a pool of nucleotide k-mers randomly generated to determine the preferred RNA sequence of an RBP [123]. Once RNAs bind to a tagged RBP, they are pulled down with a fluorescent label and measured by microarrays. This method outputs k-mers of nucleotides whose affinity to the RBPs is especially high. Using this method, the step of discovering motifs can be skipped, as RNAcompete directly provides the binding motif with high affinity to the RBPs. Only a procedure to merge motifs with high-affinity that are similar is required.

Motif discovery algorithms have been reviewed before [24], so we only briefly describe their main features to provide a wide view of the complete pipeline. In the cited review (Table 1) and in the Section 3 of the Supplementary Material, the reader can find a deeper description of these algorithms.

Most of these methods have been borrowed from the detection of transcription factor-binding sites (TFBSs). Only later, they were applied to the detection of splicing factor binding sites (SFBSs). Among all the features that differentiate TFBS from SFBS—such as SFBS specific motifs, motif length or preferable location in RNA—SFBS discovering algorithms consider, almost exclusively, the primary and secondary structure of RNA, as stated in [24].

The simplest approach to detect motifs is to use only nucleotide sequences. Within this group, one of the most widely used algorithms is MEME [101]. This tool uses probabilistic models based on the maximum likelihood estimation to look for recurring and fixed-length motifs from unaligned sequences.

The MEME algorithm belongs to a broad set of motif-based tools called MEME-suite [124], which contains several variants of this software. DREME [106] uses other models for discovering motifs, GLAM2 [103] allows finding gapped motifs with arbitrary insertions or deletions and MEME-ChIP [108] is an algorithm that performs a comprehensive motif analysis. The MEME-ChIP algorithm, additionally, incorporates other useful motif-based functions, such as analyzing the similarity of predicted motifs with known motifs (TomTom [125]), automatically grouping predicted motifs by similarity (CentriMo [126]), predicting preferred spacing between pair of motifs (SpaMo [127]) and creating a GFF file for visualizing the predicted motifs in integrative genomics viewer [128] or any genome browser.

Other algorithms, using similar approaches, enable finding ungapped motifs: phyloGibbs (which incorporates phylogeny) [102], SeAMotE [111] and cERMIT [102]; and gapped motifs: HOMER [104], ChIP- Munk [105], rGADEM [107], MatrixREDUCE [129], DRIMUST [109] and RBPmap [110].

The MEMERIS algorithm [112] is an extension of MEME that combines primary and secondary structure to find motifs. It uses the single-strandedness information of sequences as prior knowledge in the MEME’s expectation maximization model.

Other algorithms are StructRED [114], which uses mRNA expression levels in addition to the FASTA files, RNAcontext [115], which is available on the RBPmotif Web server [115], GraphProt [118], which uses a graph-based encoding, CMfinder [115], mCarts [117], TEISER [116] and RNApromo [113]. These methods are deeply described in the reviews cited in Table 1.

Several authors have compared motif-discovering algorithms [51, 130]. Jayaram et al. [51] evaluated their performance using ChIP-Seq data. In this analysis, they showed that rGADEM was the best-performing tool for discovering motifs.

Prediction of splicing regulatory factors

The final step in the pipeline is the identification of the RBPs that induce differential AS events across the conditions of the study. This section is split into two parts: scanning the SF motifs in the splicing regions and identifying the potential regulators by using some type of enrichment analysis. These tasks are depicted in Figure 4A and B.

Figure 4

Example of the pipeline for predicting context-dependent SFs (SFs). (A) Scanning SFs’ motifs in splicing regions (typically 300–400 nt upstream and downstream the AS events [131]). A set of three PWMs associated to three SFs is shown. PWMs are examined in the splicing regions three different AS events. A statistical analysis is performed to get a table of individual occurrences, which contains the hits of each motif against the events. (B) Performing a motif enrichment analysis using ΔΨ (PSI) of AS events. An example of the statistics is shown. Other tests are also possible. The main output of the pipeline is a ranking of SFs, which are predicted to regulate the splicing pattern under study.

Scanning SFs’ motifs in splicing regions

Once the SFs’ motifs are known, they are scanned across the splicing regions. This approach can potentially save costs, as binding sites can be predicted without having to use protein-centric experiments. On the other hand, the predicted binding sites can be used to make sound hypothesis on the potential regulators to be validated by an ulterior RIP or CLIP experiment.

Algorithms to scan motifs in nucleotide sequences have been deeply studied and reviewed, as they are a key element for unveiling TFBS [51, 132–134]. As it occurs with motif discovery algorithms, these methods were adapted from algorithms developed to scan TFBS. These methods can be divided into methods to find individual occurrences and methods to discover clusters of binding sites (Table 4, Supplementary Figure S2).

Table 4

Computational methods aimed at scanning motifs against DNA/RNA regions

Algorithm subtype	Main algorithms	Algorithm with best performance (reference)
Individual occurrences	Clover [142], PossumSearch [143], FIMO [135], Matrix-scan and Patser [144–146]	FIMO (Jayaram et al. [51])
Cluster of binding sites	Cister [138], Comet [147], MCAST [136], Cluster-Buster [139], CisModule [140], EMCModule [141], BayCis [137]	MCAST (Jayaram et al. [51])

FIMO [135] is a software of the MEME-Suite, which allows finding individual occurrences of motifs in DNA, RNA or protein sequences. It computes a log-likelihood ratio for each motif in each position in the given sequences and calculates the associated q-values assuming a model in which sequences are randomly generated. This method was found to outperform others when detecting TFBS [51].

Cis-regulatory modules (CRMs) are sets of RBP motifs locally enriched in the given sequences. CRM discovering algorithms return a single score for each CRM that combines the matches of its RBP motifs. A set of RBP motifs with a global significant score could provide evidence that they are acting together.

The MCAST algorithm [136] yields a list of predicted CRMs ranked by E-value. Each CRM represents a group of PWMs that frequently appear together in query sequences. MCAST was found to outperform any other algorithm to find TFBS clusters in [51].

Other algorithms used to discover CRMs such as BayCis [137], Cister [138], Cluster-Buster [139], CisModule [140] or EMCModule [141] were also reviewed in [51].

Motif enrichment analysis and refinement of results

Once the putative binding sites of RBPs in the transcriptome are known, RBPs can be associated to a gene (RxG) or to a splicing event (RxE). Methodologies that predict SFs can be divided into two broad groups depending on these relationships: methods focused on genes (RxG) and methods focused on individual AS events (RxE).

The first group consists of finding in the literature (or using the mapping) the relationships between RBPs and genes (RxG) and comparing the recurrence of RBPs in genes with spliced and non-spliced events.

Using the RxG strategy, de Miguel et al. [148] discovered the key role of protein quaking (QKI) in the regulation of splicing in non-small cell lung cancer (NSCLC). They identified the events (only cassette exons) using ExonPointer [81] and performed an enrichment analysis of genes with differentially spliced exons in different gene sets of putative regulators. QKI was found to be the most significantly enriched gene set. Experimental work showed the functional implications of the depletion of QKI in NSCLC cell lines.

A straightforward refinement is the study of individual splicing events and their potential RBP regulators—inferred by a motif scanning algorithm (RxE). RBP motifs hit regions where AS events occur. RBPs whose hits are significantly enriched in differentially spliced events are potential regulators of AS (Figure 4B).

Following this methodology, Danan-Gotthold et al. [16] analyzed splicing events with a potential role in solid tumors and predicted putative regulators of the splicing patterns for each tumor type. They developed their own algorithms to perform the motif scanning and the estimation of PSI for exon events. To assess statistical significance, they compared the frequency of occurrences between spliced and non-spliced events using a Fisher’s exact test.

In a similar approach, Sebestyén et al. [149] carried out a study of the alterations of RBPs in cancer and associated splicing changes. They performed a comprehensive analysis of 1300 RBPs in multiple tumors of TCGA. They analyzed mutation, copy number and gene expression patterns combined with AS changes and the binding motif enrichment analysis of spliced events. AS events were identified from a transcript quantification based on a known annotation using SUPPA [79]. They used FIMO [135] to scan motifs in splicing regions considering a hit if the P-value was <0.001. They evaluated the differential number of hits between spliced regions and non-spliced regions of the same size controlling for the G + C content. Finally, they measured the possible influence of RBPs by relating their expression with the splicing pattern of each event. With this methodology, they discovered that MBNL1—an SF associated with cell differentiation—controls the AS of several genes involved in the cell.

Aghamirzaie et al. [150] developed a different method called CoSpliceNet, which is based on coexpression networks of transcripts and SFs. They found RBPs that are strongly correlated with transcripts. Then, they used MEME to find conserved motifs in intron and exon sequences adjacent to events (i.e. in a cassette exon they differentiated four regions: Intron-3′, Exon-3′, Intron-5′ and Exon-5′) and found motifs in each of them separately. Finally, they identified significantly enriched motifs and constructed a co-splicing network.

We depict the application of these and other similar approaches in Table 5. This table includes the required inputs, the output of the algorithms as well as other characteristics.

Case study

To illustrate the whole pipeline and the difficulties that appear in each step, we include a worked case study. References in Table 5 perform its own analysis making a succession of decisions, such as selecting the AS detecting algorithm, downloading an RBP motif database, choosing a motif scanning algorithm or performing an enrichment test.

This case study is performed on a previous experiment with some collaborators (GSE 76902) [60]. In this experiment, the SF SRSF1 is knocked down using small interfering RNA (siRNA) on the A549 lung adenocarcinoma cell line. The experiment includes three conditions: cells treated only with the vehicle of the transfection (Lipofectamine 2000, Invitrogen), cells treated with scramble siRNA (i.e. a sequence that will not lead to the specific degradation of any cellular mRNA) and cells transfected with an siRNA that targets SRSF1. These three groups are referred to as Control, SCR and KO-SRSF1, respectively. Each condition has three biological replicates that, in turn, are hybridized three times (nine hybridizations on HTAv2 microarrays).

As a preliminary step, we compared the expression changes (Figure 5A and B) of 1243 genes that code RBPs between conditions SCR and KO-SRSF1 to confirm the knock down effect of SRSF1 and to evaluate the expression changes of other RBPs. Aroma.affymetrix pipeline was performed to summarize the expression values for each gene. Differential expression was performed using LIMMA [153]. P-values were adjusted for multiple testing using the Benjamini–Hochberg procedure. Interestingly, not only SRSF1 significantly changed its expression but also other RBP genes.

Figure 5

(A) Volcano plot of RBP genes corresponding to a LIMMA analysis that compares KO-SRSF1 versus SCR. Top five genes are highlighted. (B) Heatmap of log2 expression of the 20 most enriched RBP genes among three conditions (Control, SCR and KO-SRSF1).

We found 134 RBPs (of 1243) with differential changes of expression (adjusted P-value < 0.05 and |log₂-fold change| > 0.5). As expected, SRSF1 had the best P-value (adjusted P-value = 9.1e-29) with a log₂FC of −1.5. The differential expression of RBPs occurred in both directions but not with the same proportion (we found 72 and 28% of RBPs downregulated and upregulated, respectively).

For identifying AS events, we compared the splicing pattern of KO-SRSF1 against SCR cells. EventPointer was used to discover the AS events, as it is the only algorithm that returns the PSI value using arrays. If the experiment had been performed with RNAseq, rMats or Spladder would be the methods of choice as described in the section ‘Identification of as sites and events’. We set a filter based on the expression of genes—if the gene is not expressed, there is no point in discussing splicing. All genes whose expression was under quantile 0.25 in all the samples were discarded. Of the theoretical 97 482 events interrogated by the array, 35 963 pass the expression threshold and 3686 showed a P-value < 0.001 according to the LogFC test (approximately 4% of the events). The application of this expression filter when identifying the AS events is crucial to ignore irrelevant events, as it will be shown later.

The sequence of the neighborhood of the events (400 nt upstream or downstream, equivalent to Figure 4A: splicing regions) was extracted taking into account their corresponding strand for every splicing event interrogated by the array HTAv2. The size of the neighborhood is somehow arbitrary, but it is in a range according to [131].

To assess the validity of the 400 nt selection, we used two different CLIP data sets that target SRSF1: CLIP-seq data of HEK 293 human cell line from Sanford et al. [47] and CLIP-seq data of mouse embryo fibroblasts from Pandit et al. [46]. We mapped these data sets against the human genome and found that most CLIP hits (∼70%) were located within the selected 400 bp window (Supplementary Figure S3). The mapping between the mouse and the human genomes was performed using the liftover tool of UCSC [154, 155].

Motif enrichment analysis was performed using the PWMs from the ATtRACT database. This database contains the largest number of PWMs collected from different resources.

Many RBPs in ATtRACT include several nearly identical annotated PWMs collected from different studies. We grouped similar PWMs into a single motif using the Kullback Leibler (KL) divergence [156] (Figure 6C). If two motifs of a certain RBP are similar (KL < 0.5; Figure 6B), they are merged into a single one (for convenience, we selected the longest one). Following this criterion, we got 487 PWMs (24% of PWMs were lost).

Figure 6

(A) Proportion of Homo Sapiens’ motifs of ATtRACT database according to their length. (B)Two similar PWMs of SRSF1 which were joined together. (C) KL divergence between a set of SRSF1’s motifs. Every pair of motifs with KL divergence <0.5 were merged. (D) Each dot represents a PWM. The influence of information content (IC) and motif length with the number of significant hits (nHits; FIMO’s P-value < 1e-4) and the significance of hits (−log10 of the best hit using FIMO) is shown. PWMs with no significance hits (FIMO’s P-value < 0.001) were discarded and not shown in the figure.

We used the FIMO algorithm to scan these PWMs against the neighboring regions of AS events (as recommended in [51]). We built the background as a one-order Markov model and set the default threshold of P-values to consider a significant hit (P-value < 1e-4). We constrained FIMO to search only in the strand of the corresponding gene.

Consensus binding motifs or RBPs are often too short to get statistically significant matches. In our analysis, motifs with six nucleotides or less gave no significant hits (FIMO’s P-value < 1e-4). In the ATtRACT database, almost 44% of motifs have six nucleotides or less (Figure 6A). Interestingly, the motif length is related to the information content of motifs and to the significance of hits (Figure 6D).

There are motifs whose entropy is high, i.e. they lack well-defined binding sites. For this reason, 10% of motifs with 7 nt or more had no significant hits when scanned against the transcriptome. We finally got significant matches for 445 PWMs that correspond to 125 RBPs.

Once we have the significant hits against the event regions, we studied the significance of RBPs in differentially spliced events by using a Fisher’s exact test. In total, 14 of 125 RBPs were significantly enriched (Fisher P-value < 1e-3; Table 6). SRSF1 was one of them (Fisher P-value = 8.32E-04). However, 12 of the 14 RBPs were even more significant than SRSF1. These findings could be considered false positives, as the only direct interaction was precisely SRSF1. Interestingly, 9 of the 14 RBPs are differentially expressed (Table 6) and 13 of the 14 RBPs have strong relationships—direct or indirect—with SRSF1 according to the STRING database [157] (Figure 7) and [158]. Somehow, these false positives are showing the relationships of these RBPs in the experiment and the tight coupling among the SFs, as the depletion of SRSF1 provokes significant changes in the expression of other SFs.

Table 6

RBPs Predicted to be splicing regulators (fisher P-value < 1e-3)

RBP ranking	RBP/CLIP- seq	PWM	Fisher P-value (PSI)	LIMMA adjusted P-value (differential expression)
	Clip-seq_ SRSF1_union	–	1.25E-23	9.18E-29
	Clip-seq_ SRSF1_Pandit	–	8.30E-19	9.18E-29
1	ELAVL2	ELAVL2.1093	7.60E-15	1.35E-04
2	TIAL1	TIAL1.1287	1.75E-13	2.55E-01
3	TIA1	TIA1.1284	3.78E-11	NA in HTAv2
4	ELAVL4	ELAVL4.1095	9.17E-11	2.04E-01
	Clip-seq_ SRSF1_Sanford	–	2.28E-09	9.18E-29
5	ELAVL1	ELAVL1.161	1.39E-08	2.57E-15
6	AKAP1	AKAP1.97	1.15E-06	5.62E-07
7	ELAVL3	ELAVL3.119	2.82E-05	4.32E-01
8	SSB	SSB.58	2.64E-04	1.80E-01
9	HNRNPH2	HNRNPH2.925	3.03E-04	4.60E-02
10	TRA2A	TRA2A.s77	5.52E-04	2.40E-02
11	SF1	SF1.120	5.66E-04	5.83E-12
12	SRSF2	SRSF2.1311	6.83E-04	4.05E-16
13	SRSF1	SRSF1.1223	8.32E-04	9.18E-29
14	PTBP1	PTBP1.1012	9.13E-04	5.10E-09

Notes: CLIP-seq data of Pandit and Sanford are also included. The ranking of RBPs, the RBPs’ name, the best PWMs, the Fisher P-values of PSI and the adjusted P-values of the enrichment analysis are shown (LIMMA adjusted P-values < 1e-3 in bold).

Figure 7

String’s interactions network of the 15 significantly enriched RBP genes (Fisher P-value < 1e-3; number of nodes: 14; number of edges: 31). TIAL1 does not appear in the STRING database, but it is also related to SRSF1 [158].

It is important to note that, before applying the expression filter described above, the enrichment P-values were inaccurate (for example, the P-value of SRSF1 was nonsignificant).

We evaluated whether SRSF1 promotes exon inclusion or exclusion by comparing PSI values between KO-SRSF1 and SCR samples (selecting just the cassette exon events). SRSF1 was found to be positive splicing regulator (P-value = 5.82e-06), which is in accordance with the bibliography [45, 47].

Finally, we evaluated the enrichment of AS events in the aforementioned CLIP experiments using also a Fisher test. The enrichment for the union of both experiments was more significant than for any predicted motif (Fisher P-value = 1.25E-23). It is interesting to note that the Pandit’s CLIP-seq data were strongly significant (more than any other test), despite being data from a different organism.

According to these results, the identification of SRSF1 as driver of this change is difficult to pinpoint, as other SFs are even more significant than it. However, SRSF1’s enrichment P-value was strongly significant (in fact, it is in the position 13 of 125 RBPs). In addition, most of the other significantly enriched RBPs were differentially expressed and related to SRSF1. Even in the case of RBPs that are not differentially expressed (TIA1 and TIAL1), there is strong evidence of their interaction with SRSF1 [158] not reflected in the STRING database. Consequently, these RBPs could be participating in the regulation of AS splicing as well. The differential expression of the SFs across the studied conditions helps to filter out some SFs.

Somehow, the described experiment was optimal to ‘discover’ the SF regulating the differentially spliced events. However, the enrichment analysis alone was not sufficient to infer the key role of SRSF1. The combination of the enrichment analysis with differential expression is a must to uncover the key regulators in the experiment. In fact, most of the methods in Table 5 combine both sources of information.

Discussion and conclusion

We have outlined a conceptual computational pipeline to infer AS regulators. The first step is to detect the AS events, the second is to predict RBP-mRNA-binding sites and the last is, using both pieces of information, to predict the context-dependent SFs that regulate splicing in a specific condition.

Regarding the detection of AS events, we have already discussed qualitatively the different algorithms in the corresponding section. It would be desirable a comparison that states their performance also quantitatively. However, this task is not trivial at all: events are difficult to match across algorithms, different outputs (SI and PSI for example) can hardly be compared, a ground truth simulated experiment able to fairly compare the methods should be designed and, of course, the results should be compared in real samples with a proper validation strategy. Nevertheless, the provided comparison is still useful and can be used to guide the researcher to find the algorithms that better suit his/her needs.

Regarding the prediction of RBPs’ binding sites, the review includes different algorithms that discover motifs based on RIP and CLIP methods. These motifs—along with other ones obtained from RNAcompete techniques—are included in databases. Using these databases could potentially save costs, as there is no need to perform additional biological experiments to predict candidate binding sites of a RBP in a specific sample. It is important to point out that the broadest database collects motifs of for only around 30% of known RBPs.

Once the motifs are selected, different software packages identify the loci in the transcriptome where there are putative hits of these motifs. All the methods reviewed to scan PWMs are borrowed from the detection of TFBSs. Although the algorithms developed to scan motifs in DNA sequences can also be used with RNA sequences, this is a simplification. As previously stated, the secondary structure of RNA and the specific characteristics of SFBS play a key role in the binding process. A potential improvement of specificity and sensitivity would be achieved by this information, as some methods to extract the PWMs do it. Specifically, MEMERIS and GraphProt use the secondary structure of RNA to check single-stranded regions.

The PWMs for RBPs are usually short and repetitive and, consequently, prone to have too many potential binding sites in the transcriptome. This fact, in turn, makes it difficult to find hits that are statistically significant. As we have pointed out in the case study, only motifs >7 nt achieve statistical significance. Short motifs—6 nt or less—were discarded by the motif scanning algorithm.

The splicing machinery is complex. The interaction networks and synergistic effects of RBPs should also be considered. The process of regulation of the splicing is guided by a group of RBPs acting as a whole and not only by their individual activity (i.e. a given genomic sequence could be differently recognized by the same RBP depending on the expression of the other RBPs). This fact makes the elucidation of the regulators a much harder problem. The case study illustrates the tight control of the expression of different SFs. More than 100 SFs showed strong differential expression across the conditions. This differential expression makes it difficult to pinpoint which of the differentially expressed SFs is the driver of the change.

Finally, we summarized different works that applied these methodologies for deciphering context-dependent SFs in a certain experiment. The procedure focused on RBP–gene interactions (RxG) is the simplest pipeline, as there is no need to predict or scan RBP motifs against a transcriptome. However, using this approach, two SFs that regulate different AS events of the same genes cannot be distinguished from each other. This drawback can be resolved by analyzing RBP–event relationships (RxE).

All these methods sensibly combine information of the expression of the RBPs with overrepresentation of their putative targets in the corresponding experiment. Overrepresentation alone does not seem to be sufficient to accurately identify the drivers of the changes.

One potential reason is that the specificity and sensitivity of the methods to map PWMs are far from perfect. As stated above, one of the problems is the PWMs themselves: many of them are too short to predict the binding sites accurately. The weakest part of the pipeline is the identification of the binding sites for the RBPs: the computational prediction of these sites is prone to errors (both false positives and false negatives).

In fact, in the case study, the overrepresentation of SFBS using CLIP-seq data instead of motif-scanning was much more significant. In the long term, once the RIP- and CLIP- based techniques are settled down and results for most SF readily available, it makes more sense to use this information than scanning the putative motifs (even for different cell lines).

Despite the concerns described in this work, it is possible to predict splicing regulators with acceptable sensitivity and precision. In fact, different functional studies showed that the predictions were indeed correct. The described methodologies do not substitute RIP- and CLIP- based experiments but complement them by providing some candidates to be driver regulators in the condition under study. Besides, this approach could help the scientific community to understand the regulation networks of SFs and infer groups of SFs that cooperate in the regulation of AS.

Key Points

• Deciphering the regulation of AS is conceptually divided into three steps: detection of AS events, estimation of their interactions with SFs and contextualization for a specific experiment.

• There are many methods to detect and quantify AS events. Most of them are recent. Several algorithms that use RNA-seq data detect novel unannotated events.

• The motifs of the SFs tend to be small and repetitive making it difficult to have good precision pinpointing the binding sites.

• The algorithms that include the contextualization of the results have helped to discover novel roles in the regulation of splicing that were experimentally validated.

Acknowledgements

The authors are grateful to Francisco J. Planes, Lucía Campuzano and Xabier Cendoya for their comments on the preparation of this manuscript.

Funding

This work was supported by the Basque Government with the grant promoting doctoral theses for young pre-doctoral researcher (grant number PRE_2016_1_0194 to F.C.).

Fernando Carazo is a PhD student at the University of Navarra (Spain). He holds an MSc degree in Industrial Engineering. He is interested in deciphering the splicing rules in different diseases and in analyzing high-throughput data in translational medicine.

Juan P. Romero received his PhD in Bioinformatics at Tecnun, University of Navarra (Spain) in 2017 and now works at CIMA in the Onco-Hematology Department. His research interests include the development of algorithms to identify alternative splicing events and machine learning techniques to study the transcriptome.

Angel Rubio is an Industrial Engineer and holds a PhD by the University of Navarra (1999). His research has focused on alternative splicing analysis, analysis of copy number alterations and integration of disparate biological data sources (miRNA and mRNA expression, number of copies and expression, etc.). He has been the Director of the Biomedical Engineering Department of the University of Navarra. At present, he is a Professor of Biostatistics, Next Generation Sequencing and Genomics and Proteomics at the University of Navarra.