Functional enrichment of alternative splicing events with NEASE reveals insights into tissue identity and diseases
Genome Biology volume 22, Article number: 327 (2021)
Alternative splicing (AS) is an important aspect of gene regulation. Nevertheless, its role in molecular processes and pathobiology is far from understood. A roadblock is that tools for the functional analysis of AS-set events are lacking. To mitigate this, we developed NEASE, a tool integrating pathways with structural annotations of protein-protein interactions to functionally characterize AS events. We show in four application cases how NEASE can identify pathways contributing to tissue identity and cell type development, and how it highlights splicing-related biomarkers. With a unique view on AS, NEASE generates unique and meaningful biological insights complementary to classical pathways analysis.
Alternative splicing (AS) boosts transcript diversity in human cells  and thus contributes to tissue identity , cell development , and pathology in, e.g., cardiomyopathy , muscular dystrophy , or autoimmune diseases . It is estimated that up to 30% of disease-associated genetic variants affect splicing . RNA sequencing technologies (RNA-seq) allow the quantification of different types of AS events and detect splicing abnormalities in disorders. However, RNA-seq utility is currently limited by our incomplete understanding of the functional role of specific exons or the transcripts they contribute to.
A major challenge in AS analysis is the functional interpretation of a set of events, including isoform switching events and differentially spliced exons. The usual approach is to perform gene set enrichment or overrepresentation analysis [8,9,10]. This approach treats all genes affected by AS equally, neglecting that some AS events may not be functionally relevant at the protein level  or result from noise in the splicing machinery . Furthermore, functional differences between protein isoforms remain uncertain in many cases. A promising strategy to identify relevant AS events is to focus on those that lead to meaningful changes in the protein structure. Recent studies have shown that AS has the potential to rewire protein-protein interactions by affecting the inclusion of domain families  and linear motifs  or by activating nonsense-mediated decay .
This motivated the creation of databases and tools that predict the consequences of individual AS events or isoform switches. IsoformSwitchAnalyzeR , tappAS , DoChaP , and Spada  support transcript-level (as opposed to exon-level) analysis to identify isoform switches and their impact on the translation and the resulting isoforms features, such as domains, motifs, and non-coding sites. Exon Ontology  and DIGGER  support exon-level analysis to identify exon skipping events and their possible impact on the protein structure and function. Spada and DIGGER further consider the impact of AS on protein-protein interactions.
Most existing tools allow investigating AS-driven changes in an explorative fashion but tools for systematic analysis of functional effects of AS are lacking. Exon Ontology performs statistical tests to identify enriched features within a set of skipped exons. One example is domain families affected by AS across proteins more frequently than expected. However, none of the existing tools offer a systems biology view to specifically highlight functional consequences of AS events.
To tackle these limitations, we developed the first tool for functional enrichment of AS events. NEASE (Network-based Enrichment method for AS Events) first detects protein domains affected by AS and then uses protein-protein interactions (PPI) integrated with domain-domain interactions (DDI) , residue-level, and domain-motif interactions (DMI)  to identify interaction partners likely affected by AS. Next, it employs an edge-level hypergeometric test for gene set overrepresentation analysis. This approach is new in the way genes are selected for the enrichment test. Rather than considering only differentially spliced or expressed genes, which is currently the most common strategy, NEASE uses network information to select genes that are likely affected in the interactome. This is also superior to a simple network enrichment analysis, as we consider only those edges for which an AS contribution seems relevant and for which false positive results are less likely. We evaluated NEASE using multiple datasets from both healthy and disease cohorts. We show that the NEASE approach complements gene-level enrichment, and even outperforms it in scenarios where gene-level enrichment fails to find relevant pathways. Moreover, NEASE generates unique and meaningful biological insights on the exact impact of AS. Furthermore, since the statistical approach is network-based, NEASE can prioritize (differentially) spliced genes and find new disease biomarkers candidates in case of aberrant splicing. The NEASE Python package, freely available at https://github.com/louadi/NEASE, provides multiple functions for a deeper analysis and visualization of affected protein domains, edges, and pathways (individually or as a set).
Overview of NEASE
NEASE uses a hybrid approach that combines biological pathways with PPIs and DDIs to perform functional enrichment of AS. First, we use the structural annotation of known isoforms by mapping protein domains from the Pfam database  to the corresponding exons (Fig. 1A). Second, we construct a structural joint graph as previously reported  by enriching the BioGRID PPI  with DDIs (from DOMINE  and 3did ), DMIs from the Eukaryotic Linear Motif resource (ELM) , and interface residues from the Protein Data Bank (PDB)  (see Methods). In the joint graph, protein features such as domains, motifs, and residues are mapped to their mediated interactions. Thus, NEASE provides an exon-centric view of the interactome and addresses the limited exon-level annotation. Exons are represented by the features they encode, and interactions between features are represented by edges. In this way, the impact of AS can be seen as an edgetic change in the network. Analyzed AS events are viewed as a set of affected edges that represent gained or lost PPIs.
We then perform statistical tests to find enriched pathways and most likely responsible genes (Fig. 1B). Following, (differential) splicing analysis, a one-sided hypergeometric test is used to test for enrichment of a given pathway or gene set by considering all edges affected by AS in an experiment. A similar test is applied for each spliced gene to prioritize the most relevant events/genes that are affecting a pathway. We further introduce a weighted score (NEASE score) that penalizes hub nodes that are more likely to be connected to the pathway of interest by chance. Notably, this approach also considers genes that are not part of the existing pathway definition but show a significant number of interactions with the pathway, highlighting new putative biomarkers (see Methods and Additional file 1: Figure S6, for details).
The Python package provides an interactive analysis. Using a list of exons or events, users can run a general enrichment on 12 different pathway databases (collected from the ConsensusPathDB resource ), followed up by a specific analysis and visualization for a single affected pathway or module of interest (Fig. 1C). To provide analysis for individual isoforms and events, we linked NEASE to our previously developed database DIGGER, which provides an isoform- and exon-centric view of the interactome .
To check if the structurally annotated PPI is more biased to hubs than the standard PPI network, we computed the node degree distribution of the network before and after filtering for the structure evidence. As shown in Additional file 1: Figure S1, the two histograms show similar trends with an overall smaller number of edges in the structurally annotated PPI. The latter has a maximum node degree equal to 424, compared to 2887 in the full PPI. This observation shows that the structurally annotated PPI does not increase the bias towards hub genes of the interactome.
NEASE gives insights into the role of the muscle- and neural-specific exons
Recent studies suggest that the regulation of AS occurs in a tissue-specific manner and leads to remodeling of protein-protein interactions . Understanding the functional impact of co-regulated exons is critical in understanding gene regulation. We applied NEASE to tissue-specific exons reported in VastDB, a resource that provides information on multiple types of AS events detected by RNA-seq from different tissue types and developmental stages . We extracted 2831 exon skipping events and Percent Spliced In values (PSI) from 12 different human tissue types (Additional file 7, see Methods). We then performed hierarchical clustering on the z score standardized PSI values (Fig. 2A). The heatmap shows two distinct clusters, where neural-specific and muscle-specific (merged with heart-specific) exons are dominant.
Next, we extracted 56 skipped exons with a high PSI in the muscle tissues and 62 skipped exons with a high PSI in the neural tissues (z score ≥ + 2, see the “Methods” section). We checked how many of these events are overlapping with protein features. As shown in Fig. 2B, 27% of the upregulated exons in muscle tissues (13) and 28% of the upregulated exons in the neural tissues (17) overlap with protein features. NEASE also provides statistics of how many of these domains have known binding partners in the joint graph. In the two sets, around 60% of the affected domains have known interactions in our joint graph: 8 binding domains in the muscle tissues and 10 binding domains in the neural tissues (Additional file 2: Tables S4, S5 and Additional file 3: Tables S8, S9). We further identified one affected motif in the gene ATP2B1 in neural exons. For these groups of events, the exact protein complexes involved can be identified, and NEASE statistical analysis can be performed to determine affected pathways. However, it is important to keep in mind that not all affected domains are necessarily interacting domains but could also be regulating gene expression by binding to DNA or RNA .
First, we ran a gene set overrepresentation analysis (one-sided hypergeometric test), which we refer to as gene-level enrichment, to detect enriched pathways (see the “Methods” section). Next, we applied NEASE to the same genes to detect pathways affected by AS. Unlike the gene-level enrichment, the results obtained from NEASE in both sets better explain the functional role of the regulated exons (Fig. 2C, D). We also compared with the results from the Network Enrichment Analysis method NEA . NEA is a PPI-based approach that considers all edges for statistical tests. In contrast, NEASE considers only AS-affected edges. For a fair comparison, we run NEA with the same PPI network (BioGRID) and same pathways databases (see the “Methods” section). The results of NEA did not improve over the classic gene-level enrichment (Additional file 1: Figure S2), which suggests that our exon-specific approach helps to narrow down the exact complexes/pathways affected by AS and reduces false positives.
To further validate the robustness of the enrichment obtained by NEASE, we further conducted permutation tests. Here, our null hypothesis is that the tissue-specific exons are not different from a random set of exons in terms of the quality of the functional enrichment (measured as the p values of the hypergeometric test). For a more realistic scenario, our background set of exons considers only exon skipping events that can actually be found in these tissues (see the “Methods” section for details). This approach will also help evaluate our methods against known and unknown biases. The empirical p values of the permutation test, which indicate the chance of finding an enrichment p value as low or lower than the one reported by NEASE, are 0.0008 and 0.0001 for neural and muscles upregulated exons, respectively. These results further demonstrate the robustness of our analysis.
The upregulated exons in heart and muscle tissues were enriched in “Muscle Contraction” pathways (Fig. 2C and Additional file 2: Table S7 ), while, in the gene-level enrichment, the pathways were related to very common subcellular functions such as the Golgi apparatus, which also is an organelle for collecting, modifying or destroying protein products (Fig. 2 C and Additional file 2: Table S6). NEASE provides detailed information about the affected domains and their interaction partners (Additional file 1: Table S1). The domain Tropomyosin (Pfam id: PF00261), which is part of the gene TPM1, e.g., is involved in the regulation of muscle contraction via actin and myosin. GAS2 (Pfam id: PF02187) is a domain of DST, a dystonin encoding gene, which plays a role in maintaining the integrity of the cytoskeleton. AS affects its binding with the gene CALM1 that encodes a calcium-binding protein involved in various calcium-dependent pathways like muscle contraction .
The exons upregulated in neural tissues showed enrichment in the synaptic vesicle cycle pathway responsible for the communication between neurons (Fig. 2D). Gene-level enrichment performed on par with NEASE, resulting in the same pathway but with a lower rank and significance (adjusted p values: 1.494631e−16 using NEASE and 0.0039 using gene-level, Additional file 3: Tables S10 and S11). Notably, NEASE also detected an enrichment in “oxidative phosphorylation”, which is the initiator for powering all major mechanisms mediating brain information processing . The neuron’s energy demands are remarkable both in their intensity and in their dynamic range and quick changes [35,36,37,38]. Therefore, AS could modify oxidative phosphorylation to serve tissue-specific needs. Experimental studies have also found that several key enzymes in “oxidative phosphorylation” are spliced, e.g., pyruvate kinase (PKM) that shifts from the PKM2 to the PKM1 isoform [39, 40]. NEASE also provides a detailed view on the affected mechanisms, such as an exon skipping event in the gene ATP6V0A1 overlapping with the V_ATPase_I domain (PFAM id: PF01496) and affecting the binding with seven other proteins from the complex vacuolar ATPase (V-ATPase) (p value: 5.853289e−17, Fig. 3, Additional file 1: Table S2 ). V-ATPase is required for synaptic vesicle exocytosis  The a1-subunit of the V0 domain in ATP6V0A1 was recently shown to be highly expressed in neurons and to be essential for human brain development [42, 43]. In another example, NEASE identified two co-regulated events of the genes CLTA and CLTB (Fig. 3). CLTA and CLTB genes are involved in Clathrin-dependent endocytosis which forms clathrin-coated vesicles. Both genes play a major role in forming the protein complex of the coated vesicle. Both events affect the same domain Clathrin light chain (Pfam id: PF01086). The Clathrin light chain domain binds to CLTC and CLTCL1 which are the Clathrin heavy chain genes (p value: 6.943483e−05). These results suggest that the formation of this complex is co-regulated by AS. A similar finding about the role of the Clathrin light chain in neurons was also described in . NEASE highlights these co-regulated events at the network level (Fig. 3). As a sanity check, we manually checked the PSI values of these critical events identified by NEASE in the Genotype-Tissue Expression data set (GTEx), a comprehensive resource for tissue-specific gene expression and regulation . VastDB includes the quantification of PSI values from 8378 samples (49 tissues and 543 individuals) from GTEx version 6 on their website (https://vastdb.crg.eu/). As shown in the examples in Additional file 1 Figures S4 and S5, the exons are confirmed to be highly upregulated in their respective tissues. The analysis generated from VastDB using NEASE agrees with the latest studies at transcriptomics and proteomics levels that emphasize the crucial role of AS in the function and development of brain and heart tissues [46,47,48].
NEASE reveals splicing-related differences of reticulated and mature platelets
AS does not only drive tissue-specific regulation but also plays a major role in cell differentiation and maturation. To illustrate an example of the utility of NEASE in such studies, we used the RNA-seq data set from  which compares the transcriptome profiles of reticulated platelets and mature platelets from healthy donors. Reticulated platelets are younger , larger in size, and contain more RNA . Moreover, they have a prothrombotic potential and are known to be more abundant in patients with diabetes, acute or chronic coronary syndrome, and in smokers [51,52,53]. Additionally, elevated levels of reticulated platelets in peripheral blood are predictors of insufficient response to antiplatelet therapies (e.g., aspirin and P2Y12 inhibitors) and are promising novel biomarkers for the prediction of adverse cardiovascular events in different pathological settings [52, 54]. A strong enrichment of pro-thrombotic signaling in reticulated platelets was observed in healthy donors . Comparative transcriptomic analysis revealed a differential expression of several pathways in addition to an enrichment of prothrombotic pathways and transcripts of transmembrane proteins as the collagen receptor GPVI, the thromboxane receptor A2 and the thrombin receptors PAR1 and PAR4. Gene set enrichment analysis indicated an upregulation of entire prothrombotic activation pathways as the thrombin PAR1 and integrin GPIIb/IIIa signaling pathway in reticulated platelets.
Since AS has been described to occur in platelets , we wanted to investigate the splicing patterns between the previously defined reticulated and mature platelet subgroups. Using MAJIQ  (see the “Methods” section), we found 169 differentially spliced genes. From 25 affected protein domains, 17 have known interactions (68% of affected domains, Fig. 4A, Additional file 4: Tables S12 and S13). Other affected protein features include 6 residues involved in PPIs and one linear motif in the gene PAWR.
We observed that the enrichment at the gene-level using the Reactome  database ranks general cellular pathways higher, including “Membrane Trafficking” and “Vesicle-mediated transport,” and “Golgi-to-ER retrograde transport.” An exception is the “Circadian Clock” pathway, which is hypothesized to be related to platelet activation  (Fig. 4B). The pathway “Platelet activation, signaling and aggregation” was less significant in gene-level enrichment (adjusted p value: 0.061, Additional file 4: Table S14) compared to NEASE enrichment (adjusted p value: 0.004, Additional file 4: Table S15). Using NEASE, we obtained more meaningful results and unique pathways. As shown in Fig. 4C, the most significant pathways in reticulated platelets are G Protein-Coupled Receptor-related. G proteins are essential in the second phase of platelet-dependent thrombus formation . Furthermore, GPCR isoforms are known to have distinct signaling properties . Other relevant pathways associated with platelet activation are “Hemostasis,” “Thromboxane signaling through tp receptor,” and “Platelet homeostasis.” The full tables for enrichment at the gene level and using NEASE are available in the Additional file 4: Tables S14 and S15. The upregulation of these pathways in reticulated platelets emphasizes their previously described prothrombotic phenotype and their involvement in several downstream signaling processes.
We also looked at the individual AS events driving this enrichment. For each affected feature, NEASE tests if it significantly interacts with the GPCR downstream signaling pathway (Additional file 4: Table S16, see the “Methods” section). Figure 4C illustrates affected genes and their p value ranking. The top gene is GNAQ (G-protein subunit alpha q), which is known to be involved in signal transduction in platelets leading to platelet activation . The regulation of the G-protein alpha subunit can be an indication that compared to mature platelets, reticulated platelets are more involved in various signal transduction pathways related to, e.g., pro-thrombotic processes . PRKCA, which also showed different splicing patterns between the two platelet subgroups, plays a major role in the platelet formation process by modulating platelet function , megakaryocyte function, and development  and negatively regulates pro-platelet formation . Moreover, the regulation of PRKCA binding in reticulated platelets might refer to the young nature of reticulated platelets, which have undergone the pro-platelet formation process more recently than mature platelets [50, 65].
NEASE characterizes complex disorders such as Multiple Sclerosis
Multiple sclerosis (MS) is a chronic inflammatory demyelinating disease of the central nervous system. Early in the disease course, MS is characterized by focal lesions in the brain induced by an influx of systemic inflammatory cells. These active lesions infiltrated by immune cells and activated microglia are characterized by inflammatory demyelination and axonal loss . The surrounding white matter tissue is termed normal-appearing white matter due to diffuse pathology without focal lesion activity and dense immune activity . The etiology of MS remains unknown. Recently, a systematic literature review found 27 genes that were alternatively spliced in MS patients .
We used RNA-Seq of macrodissected areas from postmortem white matter tissue of patients with progressive MS . We compared normal-appearing white matter and active lesions regions from postmortem white matter brains of MS patients. We found 109 differentially spliced genes with 19 affected domains and one linear motif with known interactions, in addition to 6 known interacting residues. In total, NEASE identified 156 affected interactions (Additional file 5: Tables S17 and S18).
Gene-level enrichment ranks high pathways likely irrelevant that are involved in muscle contraction, cardiac conduction, and membrane trafficking, with the exception of Ca2+ ion flow across membranes (Additional file 5: Table S19). Ca2+ is an essential signal molecule for all cell activity. Although deregulation of calcium signaling is related to the pathogenesis of multiple diseases , including neurological disorders , it is not specific to neuronal tissues. In line with the neurodegenerative and immune-mediated features of MS, NEASE found unique enriched pathways related to brain network signaling and neuronal pathways “Neurotransmitter receptors and postsynaptic signal transmission,” “Transmission across Chemical Synapses,” “Activation of NMDA receptor and postsynaptic events,” “MAPK family signaling cascades,” “Neuronal System”), as well as pathways related to immune responses (“interleukin-17 signaling,” “Toll-Like Receptor 10 (TLF10) Cascade”) (Table 1 and Additional file 5: Table S20). Two other pathways were related to the uptake of anthrax or bacterial toxins. This could be a result of clean-up from toxic inflammatory processes or increased presence of invaders due to the leaky brain-blood-barrier in MS [72,73,74]. Additionally, it also supports the theory of infections as the trigger of lesion damage in MS .
As shown in Table 1, the pathway “Uptake and function of anthrax toxins” has the best overall adjusted p value, calculated only based on the total number of edges affecting the pathway. When we also included the number of significant genes and calculated NEASE scores (see the “Methods” section), NEASE ranks the pathway “Neurotransmitter receptors and postsynaptic signal transmission” first, and moves pathways such as “Transmission across Chemical Synapses” and “Neuronal System” higher in the rank. These observations illustrate the usefulness of the NEASE score as a complement to the global edge-based enrichment.
Two of the most significant genes in the “Neurotransmitter receptors” pathway were GRIN1 and GRIA1 (Additional file 5: Table S21). GRIN1 encodes GluN1, which is one of the two obligatory subunits for the NMDAR1 receptor, whereas GRIA1 encodes the AMPAR1 subunit. Their ligand is glutamate, and they are both ionotropic receptors and have been associated with MS disease severity [76,77,78]. Interestingly, AS of MAP2K4 appeared in both brain-related and immune-related pathways, significantly enriched in active lesions vs normal-appearing white matter (Table 1). MAP2K4 is a mitogen-activated protein kinase (MAPK) orchestrating multiple biological functions [79, 80]. AS of MAP2K4 has been found in rheumatoid arthritis , as well as in pathways of patients with other autoimmune diseases . MS also precedes autoimmune attack, and therefore AS of MAP2K4 in active lesions detected with NEASE may represent dysregulated immune responses originating from the infiltrating immune cells or inflammatory-activated brain cells. This is supported by previous studies that found (i) overactivity of MAPK pathways in microglia (the resident immune cell of the brain) during neurodegeneration [83, 84], and (ii) increased phosphorylation of MAPK kinases in the systemic immune cells of MS patients [85, 86]. A recent study also characterized activated MS-specific pathways in immune cells from blood using phosphoproteomics. Here, MAP2K4 and its interaction partners (e.g., TAK1) were present in MS-specific signaling activity . Future functional studies on the AS of MAP2K4 may help explain if AS could be the reason for increased phosphorylation and overactivity detected in MS. AS of MAP2K4 could result in switching protein conformation, increasing susceptibility to phosphorylation, or changing the downstream protein cascade.
With NEASE, we were able to specifically detect AS of genes and related pathways already known to be dysregulated within MS from excitotoxicity to inflammation. The detected AS genes in active lesions vs normal-appearing white matter demonstrate how major components in signaling activities may be fine-tuned/changed from regulation of a homeostatic state to an inflammatory state. Combining NEASE with functional experiments to understand the biological impact of AS could fuel new therapeutic opportunities for complex neurological diseases as MS. Novel developments in genome-editing tools and gene-specific strategies have made it possible to use antisense oligonucleotides or small modulators for splice modification. This is already used in the rare neuromuscular disease, spinal muscular atrophy, where an antisense oligonucleotide binds to a site near splicing to ensure the inclusion of an exon during the splicing event .
NEASE finds new biomarker candidates for dilated cardiomyopathy
AS might play a role in driving dilated cardiomyopathy (DCM) . DCM is a common heart muscle disease that is often diagnosed with structural abnormalities resulting in impaired contraction. Previous studies have shown a large number of differentially used exons in DCM patients [4, 10]. In this analysis, we used a list of 1212 differentially used exons between DCM patients and controls as reported by Heinig et al. . 29% of these exons overlap with protein features, including 230 domains and 15 linear motifs. (Additional file 6: Tables S22 and S23). In this exon set, both the gene level enrichment and NEASE show very similar results (Additional file 6: Tables S24 and S25). In both methods, we found that the list of exons was enriched in the dilated cardiomyopathy (DCM) pathway from KEGG, as well as, “Adrenergic signaling in cardiomyocytes” and “Regulation of actin cytoskeleton”.
In contrast to gene-level enrichment analysis, NEASE is able to score the contribution of alternatively spliced genes that are interacting with but are not part of the DCM pathway, allowing us to highlight putative biomarkers (Table 2, Additional file 6: Table S26, Additional file 1: Figure S3). The Myosin head domain from the gene MYO19 interacts with 6 other genes associated with DCM: (1) MYL2, which triggers contraction after Ca+ activation ; (2–5) TPM1/TPM2/TPM3/TPM4, which encode the TPM protein—the main regulator of muscle contraction ; and (6) ACTG, which encodes actin. Interestingly, MYO19 has not been investigated for its role in DCM, while its interacting genes are associated with DCM [91,92,93,94]. Additionally, the gene OBSCN has one affected interaction with the TTN gene . The TTN gene itself is also differentially spliced and associated with DCM . OBSCN was recently reported as a new DCM candidate [96, 97]. Another interesting example is CACNA1C (Calcium Voltage-Gated Channel Subunit Alpha1 C), an already known DCM candidate . The differentially spliced exon overlaps with the domain Ion_trans (Pfam id: PF00520) which is essential for myocyte contraction . The affected interaction identified is with the ryanodine receptor 2 (RYR2). In striated muscles, the excitation-contraction coupling is mediated by this complex . Both CACNA1C and RYR2 are part of the KEGG DCM pathway . Alterations in ryanodine receptors were repeatedly reported to be related to heart failure [102,103,104].
In spite of its importance for biomarker and therapeutic target discovery, differential AS is still not a routine part of transcriptome analysis. A key reason for this could be the lack of suitable methods and software tools for AS-specific functional analysis. Our method NEASE closes this gap and provides a unique view on the impact of AS complementary to functional insights gained from traditional gene-level enrichment analysis. We applied NEASE to four diverse data sets and show that its results generate novel disease-relevant insights and provide valuable context to prior findings on altered RNA- and protein-expression levels consistent with recent literature.
In many cases, NEASE improves over gene-level enrichment analysis focusing on differentially spliced genes. One potential reason for this could be that not all AS events are necessarily functional [11, 12]. NEASE mitigates this by focusing on AS events that affect protein domains. However, it is important to keep in mind that this is not the only way to define functional AS events. AS also affects interacting disordered regions  or facilitates nonsense-mediated decay .
AS events could also lead to completely different functions or interactions , e.g., two isoforms can have different interaction partners depending on the inclusion or loss of a single domain . Such changes in the interactome can not be captured with gene-level enrichment which has a strict focus on nodes rather than edges. With NEASE, we could show that integrating structural information at the exon level and PPI networks helps to identify the functional impact of differentially spliced and co-regulated exons. In practice, we consider both approaches as complementary and recommend running gene-level and edge-level enrichments together (both supported by the NEASE package). Note that while our analysis focuses on exon skipping events as the most studied event type, our method is generally agnostic to the event type.
NEASE relies on structurally annotated interactions and existing pathway annotations from databases such as KEGG  and Reactome . Leveraging reliable structural information and established pathways likely removes many false positive PPI from considerations. While the structural annotations are generally of high quality, it should be noted that their coverage is still limited and, thus, the number of exons considered in our method is comparably low. For instance, the percentage of considered exons, in our example datasets, ranges between 15 and 30%, which is still far from being a global analysis of AS. Expanding the annotations at isoform-level and more widespread availability of structural information will greatly raise the usefulness of NEASE in the future. We also emphasize that while all events can potentially affect protein interactions on the domain level, not all AS events yield functional isoforms and other processes such as nonsense-mediated decay need to be considered as well. In the future, further progress is urgently needed to link transcriptomics and proteomics for better characterization and understanding of the exact impact of AS events. With our current approach, a large fraction of the PPI network remains unexplored, suggesting that adapting de novo network enrichment methods such as KeyPathwayMiner  towards AS could be a promising research direction to uncover previously unknown disease mechanisms. NEASE currently considers the immediate neighborhood of a pathway in the PPI network. When carefully considering the expected increase in false positives, one could also increase the size of the pathway neighborhood using, e.g., a fixed radius for shortest paths. While these are attractive approaches, the biases of the PPI towards hubs, as well as the high number of false (or missing) edges of PPI, in its current form, make such approaches hard to control and statistically challenging. Even though NEASE is relatively conservative, we demonstrated that it is simple, robust, and generates meaningful and interpretable results. Thus, it provides an unprecedented opportunity to understand the functional impact of tissue-, developmental- and disease-specific AS in a system biology manner.
While a plethora of gene set enrichment methods have been proposed in recent years, AS is typically not addressed specifically. Thus, NEASE closes an important gap in functional enrichment analysis of transcriptomics data. The analyses described here, confirm the widespread impact of AS in multiple biological processes and disorders. In the future, we plan to extend NEASE with further model organisms and to add structural annotations covering more types of AS events. Finally, we plan to integrate NEASE with the DIGGER web tool  for seamless downstream analysis of AS in the web browser with the vision of establishing functional AS event analysis as a routine step in the transcriptomic analysis.
NEASE data sources
We construct a human structurally annotated PPI as described previously . Briefly, we integrate DDI and PPI information into a joint network where DDIs were obtained from 3did (v2019_01 ) and DOMINE (v2.0  including high- and mid-confidence interactions) and PPIs were obtained from BioGRID 3.5 . In summary, out of 410,961 interactions from the human interactome 52,467 have at least one domain interaction. The linear motif instances and their interactions were downloaded from the ELM website and mapped to the respective exons. We found 3926 PPIs that are confirmed by at least one source of DMI. Position-specific PPI based on experimentally resolved structure from the PDB was obtained from . In total, 16,161 PPIs were enriched by at least one residue-level interaction. From the combination of all these resources the final structurally annotated graph contained in total 60,235 interactions. Each one of these interactions is annotated with one or multiple levels of evidence (DDIs, DMI, residues). The mapping of exons to their protein features was performed using the Biomart mapping table, Pfam, and ELM annotations [22, 23, 108]. We obtain the biological pathways with their gene list from KEGG  and Reactome  integrated into the ConsensusPathDB database .
Statistical tests and pathway scores
Gene-level enrichment is performed using a hypergeometric test from the package GSEAPY (a Python wrapper for Enrichr ) by considering all genes with (differential) AS events. Network enrichment analysis at the gene level was performed using the EviNet web server (www.evinet.org), which is an implementation of the randomization algorithm NEA . To achieve a fair comparison with NEASE, we run NEA using BioGRID as a PPI database and Reactome and KEGG as pathways references to match the exact conditions of the NEASE analysis.
For NEASE enrichment, we filtered the PPI graph G=(V, E), where V is the set of genes and E is the set of edges, to a subgraph G′=(V′, E′) containing only structurally annotated interactions E′ and their nodes V′. An interaction is considered structurally annotated if it is supported by at least one of these resources: domain-domain interactions, motif-domain interactions, or residue-level interactions. For a submitted query list of exons, NEASE first identifies affected domains, linear motifs, and residues that overlap with the exons and their interactions. Let N be two times the number of edges in G′ (the degree of the network) and n be the number of affected edges from the query. These edges are then considered using a test modified from . For every pathway P with degree K, let k be the number of affected edges that are connected to P. We model X whose outcome is k as a random variable following a hypergeometric distribution:
where k is considered as the number of observed successes out of n draws, from a population of size N containing K success. Subsequently, NEASE tests if the number k is significant using a one-sided hypergeometric test (over-representation). In contrast to the test proposed in , our test only includes structurally annotated edges and the ones likely to be impacted by AS in order to improve the signal-to-noise ratio. For illustration purposes, in the example of Fig. 1B, the overall number of affected edges by AS is n=7, and K=11 is the total degree of pathway A (11 possible success), the number of affected edges that are linked to the pathway is k=4. The enrichment p value of pathway A corresponds then to the significance of this last number. After testing for multiple pathways, the obtained p values for the edge-level enrichment are corrected, using the Benjamini-Hochberg method . The detailed pseudocode of this algorithm is explained in Additional file 1: Figure S6, Algorithm 1.
For a pathway of interest, a similar test can be applied to determine if a splicing event significantly affects interactions of a specific gene with this pathway (Additional file 1: Figure S6, Algorithm 2 and Fig. 1B, C). Here, n is the number of all affected interactions (edges) of a spliced gene and k is the number of affected interactions (edges) across genes that are linked to the pathway of interest. In the example of Fig. 1B, C, the gene G2 can be connected to pathway A just by chance due to its high number of affected interactions. For this reason, it is ranked lower than the genes G3 and G4.
As a result, for every pathway, NEASE provides an overall p value, as well as the most significant genes. Since the p value only depends on the overall number of affected edges but not on the number of genes, the p value can be heavily influenced by hub genes. To reduce this influence, an optional score (NEASE Score, Eq.: 1) can be computed by NEASE to scale the natural logarithm of the p value with the total number of significant genes using a cutoff from the user (for instance p value ≤ 0.05):
where g is the total number of significantly connected genes obtained after testing individual spliced genes. Thus, the NEASE Score prioritizes pathways that are affected by a larger number of spliced genes rather than pathways that have a larger number of affected interactions (edges). The user can choose to rank enrichment based on the adjusted p value or by the NEASE score.
For permutation tests, the set of highly confident 2831 exon skipping events obtained after the initial quality filtering is considered as the background set of exons. We compared the enrichment obtained in the pathways “Muscle Contraction” and “Synaptic vesicle cycle” from the set of exons upregulated in muscles/heart and neural tissues respectively, with 10,000 random sets of exons of the same size and then derived distribution of p values. The empirical p values were then obtained by asking how likely it is to obtain a p value as low or more extreme than the one reported by NEASE in the original set of neural and muscles upregulated exons.
VastDB events processing
PSI values of the exon skipping events from VastDB were quantified by the developers using vast tools [30, 112]. In our analysis, we extracted the PSI values for 32 experiments belonging to 12 main tissues: muscles/heart, neural (whole brain, cortex, and peripheral retina), placental, epithelial, digestive (colon and stomach), liver, kidney, adipose, testis, immune-hematopoietic, and ovary. We then filtered out the events with low read coverage (VLOW) and performed hierarchical clustering of standardized values (z scores). For every exon, we calculated the mean of PSI values from the samples of the same tissues. To extract muscles/heart and neural-specific exons and to ensure that we only consider functional events, we applied two filters: namely that the exon PSI value in the relevant tissue is higher than 20 and that the z score is higher than 2.
Raw RNA-Seq reads for two types of platelets and multiple sclerosis patients were downloaded from the GEO repository (access numbers: GSE126448 and GSE138614 ). The number of samples and sequencing depth are reported in Additional file 1: Table S3. RNA-Seq reads were aligned to the reference human genome (hg38) using STAR 2.7  in a 2-pass mode and filtered for uniquely mapped reads. Differential AS analysis was performed by MAJIQ  with default parameters, and with a threshold of P(dPSI > 20%) > 0.95.
NEASE: The Python package
NEASE’s Python package relies on NumPy , pandas , NetworkX , SciPy , and Statsmodels . The gene-level enrichment is also supported in the NEASE package using the Python implementation of Enrichr . To speed up the edge hypergeometric test, the total degree of every pathway in the structural PPI, as well as the overall degree of the network were pre-computed. For visualization, we use the complete PPI (not the structural PPI) and extract connected subnetworks from each pathway as well as spliced genes and their interactions with the extracted modules. The position of nodes is computed using the Fruchterman-Reingold force-directed algorithm implemented in NetworkX . The interactive visualization for individual genes and events is implemented with information from the DIGGER database and the Plotly package.
The package provides the option to automatically filter exons that are likely to disturb the open reading frame of the transcript based on the prediction in . In the case of multiple AS events affecting the same genes, we consider every event individually and identify all protein features. The standard input of the package is a DataFrame object with the exon coordinates and Ensembl IDs of the genes. The package also supports the output of multiple AS differential detection tools such as rMATs , Whippet , and also tools that are event-based such as MAJIQ  where NEASE only considers annotated exons. NEASE is released as open-source under the GPLv3 license and is available at (https://github.com/louadi/NEASE). Step-by-step tutorials for running NEASE are available at (https://github.com/louadi/NEASE-tutorials).
Availability of data and materials
RNA sequencing data for reticulated platelets was provided by the authors  and it is freely available at GEO (access number: GSE126448). Multiple sclerosis Raw sequence was provided by the authors  and freely available at GEO (access number: GSE138614). Dilated Cardiomyopathy raw data is available in the European Genome-phenome Archive (Dataset ID: EGAS00001002454), in our analysis, we used pre-processed data from the manuscript . VastDB dataset for humans (hg19) was downloaded from https://vastdb.crg.eu/wiki/Downloads. The linear motifs instances and interactions were downloaded from (http://elm.eu.org/). The generated joint graphs and the exon mapping databases are available on the DIGGER database website https://exbio.wzw.tum.de/digger/download. NEASE is released as open-source under the GPLv3 license and is available at GitHub  and deposited to Zenodo . All processed datasets, as well as step-by-step tutorials for using NEASE to reproduce the results presented in this paper, are available at https://github.com/louadi/NEASE-tutorials and deposited to Zenodo .
Stamm S, Ben-Ari S, Rafalska I, Tang Y, Zhang Z, Toiber D, et al. Function of alternative splicing. Gene. 2005;344:1–20. https://doi.org/10.1016/j.gene.2004.10.022.
Yeo G, Holste D, Kreiman G, Burge CB. Variation in alternative splicing across human tissues. Genome Biol. 2004;5(10):R74. https://doi.org/10.1186/gb-2004-5-10-r74.
Baralle F, Giudice J. Alternative splicing as a regulator of development and tissue identity. Nat Rev Mol Cell Biol. 2017;18(7) Available from:. https://doi.org/10.1038/nrm.2017.27.
Beqqali A. Alternative splicing in cardiomyopathy. Biophys Rev. 2018;10(4):1061–71.
Douglas AGL, Wood MJA. Splicing therapy for neuromuscular disease. Mol Cell Neurosci. 2013;56:169–85. https://doi.org/10.1016/j.mcn.2013.04.005.
Evsyukova I, Somarelli JA, Gregory SG, Garcia-Blanco MA. Alternative splicing in multiple sclerosis and other 673 autoimmune diseases. RNA Biology. 2010;7(4)462–73. https://doi.org/10.4161/rna.7.4.12301.
López-Bigas N, Audit B, Ouzounis C, Parra G, Guigó R. Are splicing mutations the most frequent cause of hereditary disease? FEBS Lett. 2005;579(9):1900–3. https://doi.org/10.1016/j.febslet.2005.02.047.
Karlebach G, Veiga DFT, Mays AD, Chatzipantsiou C, Barja PP, Chatzou M, et al. The impact of biological sex on alternative splicing. bioRxiv. 2020:490904. https://doi.org/10.1101/490904.
Tollervey JR, Wang Z, Hortobágyi T, Witten JT, Zarnack K, Kayikci M, et al. Analysis of alternative splicing associated with aging and neurodegeneration in the human brain. Genome Res. 2011;21(10):1572–82. https://doi.org/10.1101/gr.122226.111.
Heinig M, Adriaens ME, Schafer S, van Deutekom HWM, Lodder EM, Ware JS, et al. Natural genetic variation of the cardiac transcriptome in non-diseased donors and patients with dilated cardiomyopathy. Genome Biol. 2017;18(1):170. https://doi.org/10.1186/s13059-017-1286-z.
Tress ML, Abascal F, Valencia A. Alternative Splicing May Not Be the Key to Proteome Complexity. Trends Biochem Sci. Elsevier Ltd. 2017;42(2):98–110. Available from: https://pubmed.ncbi.nlm.nih.gov/27712956/. https://doi.org/10.1016/j.tibs.2016.08.008.
Melamud E, Moult J. Stochastic noise in splicing machinery. Nucleic Acids Res. 2009;37(14):4873–86. https://doi.org/10.1093/nar/gkp471.
Yang X, Coulombe-Huntington J, Kang S, Sheynkman GM, Hao T, Richardson A, et al. Widespread Expansion of Protein Interaction Capabilities by Alternative Splicing. Cell. 2016;164(4):805–17. https://doi.org/10.1016/j.cell.2016.01.029.
Buljan M, Chalancon G, Eustermann S, Wagner GP, Fuxreiter M, Bateman A, et al. Tissue-specific splicing of disordered segments that embed binding motifs rewires protein interaction networks. Mol Cell. 2012;46(6):871–83.
da Costa PJ, Menezes J, Romão L. The role of alternative splicing coupled to nonsense-mediated mRNA decay in human disease. Int J Biochem Cell Biol. 2017;91(Pt B):168–75.
Kristoffer V-S. Sandelin A. Genomics: The Landscape of Isoform Switches in Human Cancers; 2017; Available from:. https://doi.org/10.1158/1541-7786.MCR-16-0459.
delafuente L, Arzalluz-luque Á, Tardáguila M, Delrisco H, Martí C, Tarazona S, et al. tappAS: a comprehensive computational framework for the analysis of the functional impact of differential splicing. Genome Biol. 2020;21(1) Available from:. https://doi.org/10.1186/s13059-020-02028-w.
Gal-Oz ST, Haiat N, Eliyahu D, Shani G, Shay T. DoChaP: the domain change presenter. Nucleic Acids Res. 2021;49(W1):W162–8. https://doi.org/10.1093/nar/gkab357.
Ctor Climente-Gonzá Lez H, Porta-Pardo E, Godzik A, Correspondence EE, Eyras E. The Functional Impact of Alternative Splicing in Cancer. Cell Rep. 2017;20:2215–26.
Tranchevent L-C, Aubé F, Dulaurier L, Benoit-Pilven C, Rey A, Poret A, et al. Identification of protein features encoded by alternative exons using Exon Ontology. Genome Res. 2017;27(6):1087–97.
Louadi Z, Yuan K, Gress A, Tsoy O, Kalinina OV, Baumbach J, et al. DIGGER: exploring the functional role of alternative 709 splicing in protein interactions. Nucleic Acids Res. 2020;49(D1):D309-D318. https://doi.org/10.1093/nar/gkaa768.
Kumar M, Gouw M, Michael S. Amano-S ´ Anchez HS´, Pancsa R, Glavina J, et al. ELM-the eukaryotic linear motif resource in 2020. Nucleic Acids Res. 2020;48. Available from: https://academic.oup.com/nar/article/48/D1/D296/5611669. https://doi.org/10.1093/nar/gkz1030.
El-Gebali S, Mistry J, Bateman A, Eddy SR, Luciani A, Potter SC, et al. The Pfam protein families database in 2019. Nucleic Acids Res. 2019;47(D1):D427–32. https://doi.org/10.1093/nar/gky995.
Oughtred R, Stark C, Breitkreutz B-J, Rust J, Boucher L, Chang C, et al. The BioGRID interaction database: 2019 update. Nucleic Acids Res. 2019;47(D1):D529–41. https://doi.org/10.1093/nar/gky1079.
Yellaboina S, Tasneem A, Zaykin DV, Raghavachari B, Jothi R. DOMINE: a comprehensive collection of known and predicted domain-domain interactions. Nucleic Acids Res. 2011;39(Database issue):D730–5. https://doi.org/10.1093/nar/gkq1229.
Mosca R, Céol A, Stein A, Olivella R, Aloy P. 3did: a catalog of domain-based interactions of known three-dimensional structure. Nucleic Acids Res. 2014;42(Database issue):D374–9. https://doi.org/10.1093/nar/gkt887.
Berman HM, Westbrook J, Feng Z, Gilliland G, Bhat TN, Weissig H, et al. The Protein Data Bank. Nucleic Acids Res. 2000;28(1):235–42. https://doi.org/10.1093/nar/28.1.235.
Kamburov A, Stelzl U, Lehrach H, Herwig R. The ConsensusPathDB interaction database: 2013 update. Nucleic Acids Res. 2013;41(Database issue):D793–800.
Ellis JD, Barrios-Rodiles M, Çolak R, Irimia M, Kim T, Calarco JA, et al. Tissue-Specific Alternative Splicing Remodels Protein-Protein Interaction Networks. Mol Cell. 2012;46(6):884–92.
Tapial J, Ha KCH, Sterne-Weiler T, Gohr A, Braunschweig U, Hermoso-Pulido A, et al. An atlas of alternative splicing profiles and functional associations reveals new regulatory programs and genes that simultaneously express multiple major isoforms. Genome Res. 2017;27(10):1759–68.
Seo PJ, Kim MJ, Ryu J-Y, Jeong E-Y, Park C-M. Two splice variants of the IDD14 transcription factor competitively form nonfunctional heterodimers which may regulate starch metabolism. Nat Commun. 2011;2(1):303. https://doi.org/10.1038/ncomms1303.
Alexeyenko A, Lee W, Pernemalm M, Guegan J, Dessen P, Lazar V, et al. Network enrichment analysis: extension of gene-set enrichment analysis to gene networks. BMC Bioinformatics. 2012;13(1):226. https://doi.org/10.1186/1471-2105-13-226.
Tansey MG, Luby-Phelps K, Kamm KE, Stull JT. Ca(2+)-dependent phosphorylation of myosin light chain kinase decreases the Ca2+ sensitivity of light chain phosphorylation within smooth muscle cells. J Biol Chem. 1994;269(13):9912–20.
Hall CN, Klein-Flügge MC, Howarth C, Attwell D. Oxidative phosphorylation, not glycolysis, powers presynaptic and postsynaptic mechanisms underlying brain information processing. J Neurosci. 2012;32(26):8940–51. https://doi.org/10.1523/JNEUROSCI.0026-12.2012.
Sanganahalli BG, Herman P, Blumenfeld H, Hyder F. Oxidative neuroenergetics in event-related paradigms. J Neurosci. 2009;29(6):1707–18. https://doi.org/10.1523/JNEUROSCI.5549-08.2009.
Vergara RC, Jaramillo-Riveri S, Luarte A, Moënne-Loccoz C, Fuentes R, Couve A, et al. The Energy Homeostasis Principle: Neuronal Energy Regulation Drives Local Network Dynamics Generating Behavior. Front Comput Neurosci. 2019;13:49. https://doi.org/10.3389/fncom.2019.00049.
Du F, Zhu X-H, Zhang Y, Friedman M, Zhang N, Ugurbil K, et al. Tightly coupled brain activity and cerebral ATP metabolic rate. Proc Natl Acad Sci U S A. 2008;105(17):6409–14. https://doi.org/10.1073/pnas.0710766105.
Howarth C, Gleeson P, Attwell D. Updated energy budgets for neural computation in the neocortex and cerebellum. J Cereb Blood Flow Metab. 2012;32(7):1222–32. https://doi.org/10.1038/jcbfm.2012.35.
Magistretti PJ, Allaman I. A cellular perspective on brain energy metabolism and functional imaging. Neuron. 2015;86(4):883–901. https://doi.org/10.1016/j.neuron.2015.03.035.
Zheng X, Boyer L, Jin M, Mertens J, Kim Y, Ma L, et al. Metabolic reprogramming during neuronal differentiation from aerobic glycolysis to neuronal oxidative phosphorylation. Elife. 2016;10:5. Available from:. https://doi.org/10.7554/eLife.13374.
Hiesinger PR, Fayyazuddin A, Mehta SQ, Rosenmund T, Schulze KL, Zhai RG, et al. The v-ATPase V0 subunit a1 is required for a late step in synaptic vesicle exocytosis in Drosophila. Cell. 2005;121(4):607–20. https://doi.org/10.1016/j.cell.2005.03.012.
Aoto K, Kato M, Akita T, Nakashima M, Mutoh H, Akasaka N, et al. ATP6V0A1 encoding the a1-subunit of the V0 domain of vacuolar H+-ATPases is essential for brain development in humans and mice. Nat Commun. 2021;12(1):2107. https://doi.org/10.1038/s41467-021-22389-5.
Poëa-Guyon S, Amar M, Fossier P, Morel N. Alternative splicing controls neuronal expression of v-ATPase subunit a1 and sorting to nerve terminals. J Biol Chem. 2006;281(25):17164–72. https://doi.org/10.1074/jbc.M600927200.
Redlingshöfer L, McLeod F, Chen Y, Camus MD, Burden JJ, Palomer E, et al. Clathrin light chain diversity regulates membrane deformation in vitro and synaptic vesicle formation in vivo. Proc Natl Acad Sci U S A. 2020;117(38):23527–38.
Lonsdale J, Thomas J, Salvatore M, Phillips R, Lo E, Shad S, et al. The Genotype-Tissue Expression (GTEx) project. Nat Genet. 2013;45(6):580–5. https://doi.org/10.1038/ng.2653.
Rodriguez JM, Pozo F, di Domenico T, Vazquez J, Tress ML. An analysis of tissue-specific alternative splicing at the protein level. Orengo CA, editor. PLoS Comput Biol. 2020;16(10):e1008287.
Raj B, Blencowe BJ. Alternative Splicing in the Mammalian Nervous System: Recent Insights into Mechanisms and Functional Roles. Neuron. 2015;87(1):14–27. https://doi.org/10.1016/j.neuron.2015.05.004.
Su C-H, Dhananjaya D, Tarn W-Y. Alternative Splicing in Neurogenesis and Brain Development. Front Mol Biosci. 2018;5:12.
Bongiovanni D, Santamaria G, Klug M, Santovito D, Felicetta A, Hristov M, et al. Transcriptome Analysis of Reticulated Platelets Reveals a Prothrombotic Profile. Thromb Haemost. 2019;119(11):1795–806. https://doi.org/10.1055/s-0039-1695009.
Ault KA, Knowles C. In vivo biotinylation demonstrates that reticulated platelets are the youngest platelets in circulation. Exp Hematol. 1995;23(9):996–1001.
Karpatkin S. Heterogeneity of human platelets. II. Functional evidence suggestive of young and old platelets. J Clin Invest. 1969;48(6):1083–7.
Cesari F, Marcucci R, Gori AM, Caporale R, Fanelli A, Casola G, et al. Reticulated platelets predict cardiovascular death in acute coronary syndrome patients. Thromb Haemost. 2013;109(05):846–53. https://doi.org/10.1160/TH12-09-0709.
Guthikonda S, Alviar CL, Vaduganathan M, Arikan M, Tellez A, DeLao T, et al. Role of reticulated platelets and platelet size heterogeneity on platelet activity after dual antiplatelet therapy with aspirin and clopidogrel in patients with stable coronary artery disease. J Am Coll Cardiol. 2008;52(9):743–9. https://doi.org/10.1016/j.jacc.2008.05.031.
Muronoi T, Koyama K, Nunomiya S, Lefor AK, Wada M, Koinuma T, et al. Immature platelet fraction predicts coagulopathy-related platelet consumption and mortality in patients with sepsis. Thromb Res. 2016;144:169–75. https://doi.org/10.1016/j.thromres.2016.06.002.
Nassa G, Giurato G, Cimmino G, Rizzo F, Ravo M, Salvati A, et al. Splicing of platelet resident pre-mRNAs upon activation by physiological stimuli results in functionally relevant proteome modifications. Sci Rep. 2018;8(1):498. https://doi.org/10.1038/s41598-017-18985-5.
Vaquero-Garcia J, Barrera A, Gazzara MR, González-Vallinas J, Lahens NF, Hogenesch JB, et al. A new view of transcriptome complexity and regulation through the lens of local splicing variations. Elife. 2016;5:e11752.
Fabregat A, Jupe S, Matthews L, Sidiropoulos K, Gillespie M, Garapati P, et al. The Reactome Pathway Knowledgebase. Nucleic Acids Res. 2018;46(D1):D649–55. https://doi.org/10.1093/nar/gkx1132.
Scheer FAJL, Michelson AD, Frelinger AL 3rd, Evoniuk H, Kelly EE, McCarthy M, et al. The human endogenous circadian system causes greatest platelet activation during the biological morning independent of behaviors. PLoS One. 2011;6(9):e24549.
Offermanns S. Activation of Platelet Function Through G Protein–Coupled Receptors. Circ Res. 2006;99(12):1293–304. https://doi.org/10.1161/01.RES.0000251742.71301.16.
Marti-Solano M, Crilly SE, Malinverni D, Munk C, Harris M, Pearce A, et al. Combinatorial expression of GPCR isoforms affects signalling and drug responses. Nature. 2020;587(7835):650–6. https://doi.org/10.1038/s41586-020-2888-2.
Jalagadugula G, Dhanasekaran DN, Kim S, Kunapuli SP, Rao AK. Early growth response transcription factor EGR-1 regulates Galphaq gene in megakaryocytic cells. J Thromb Haemost. 2006;4(12):2678–86.
Moore SF, van den Bosch MTJ, Hunter RW, Sakamoto K, Poole AW, Hers I. Dual regulation of glycogen synthase kinase 3 (GSK3)α/β by protein kinase C (PKC)α and Akt promotes thrombin-mediated integrin αIIbβ3 activation and granule secretion in platelets. J Biol Chem. 2013;288(6):3918–28. https://doi.org/10.1074/jbc.M112.429936.
Harper MT, Poole AW. Diverse functions of protein kinase C isoforms in platelet activation and thrombus formation. J Thromb Haemost. 2010;8(3):454–62. https://doi.org/10.1111/j.1538-7836.2009.03722.x.
Williams CM, Harper MT, Poole AW. PKCα negatively regulates in vitro proplatelet formation and in vivo platelet production in mice. Platelets. 2014;25(1):62–8.
Ault KA, Rinder HM, Mitchell J, Carmody MB, Vary CP, Hillman RS. The significance of platelets with increased RNA content (reticulated platelets). A measure of the rate of thrombopoiesis. Am J Clin Pathol. 1992;98(6):637–46. https://doi.org/10.1093/ajcp/98.6.637.
Bö L, Dawson TM, Wesselingh S, Mörk S, Choi S, Kong PA, et al. Induction of nitric oxide synthase in demyelinating regions of multiple sclerosis brains. Ann Neurol. 1994;36(5):778–86. https://doi.org/10.1002/ana.410360515.
Ludwin SK. The pathogenesis of multiple sclerosis: relating human pathology to experimental studies. J Neuropathol Exp Neurol. 2006;65(4):305–18. https://doi.org/10.1097/01.jnen.0000225024.12074.80.
Hecker M, Rüge A, Putscher E, Boxberger N, Rommer PS, Fitzner B, et al. Aberrant expression of alternative splicing variants in multiple sclerosis - A systematic review. Autoimmun Rev. 2019;18(7):721–32. https://doi.org/10.1016/j.autrev.2019.05.010.
Elkjaer ML, Frisch T, Reynolds R, Kacprowski T, Burton M, Kruse TA, et al. Molecular signature of different lesion types in the brain white matter of patients with progressive multiple sclerosis. Acta Neuropathol Commun. 2019;7(1):205. https://doi.org/10.1186/s40478-019-0855-7.
Gissel H. Ca2+ accumulation and cell damage in skeletal muscle during low frequency stimulation. Eur J Appl Physiol. 2000;83(2-3):175–80. https://doi.org/10.1007/s004210000276.
Maléth J, Hegyi P. Ca2+ toxicity and mitochondrial damage in acute pancreatitis: translational overview. Philos Trans R Soc Lond B Biol Sci. 2016;5(1700):371(1700). Available from:. https://doi.org/10.1098/rstb.2015.0425.
Minagar A, Alexander JS. Blood-brain barrier disruption in multiple sclerosis. Mult Scler. 2003;9(6):540–9. https://doi.org/10.1191/1352458503ms965oa.
Claudio L, Raine CS, Brosnan CF. Evidence of persistent blood-brain barrier abnormalities in chronic-progressive multiple sclerosis. Acta Neuropathol. 1995;90(3):228–38.
Ortiz GG, Pacheco-Moisés FP, Macías-Islas MÁ, Flores-Alvarado LJ, Mireles-Ramírez MA, González-Renovato ED, et al. Role of the blood-brain barrier in multiple sclerosis. Arch Med Res. 2014;45(8):687–97.
Ascherio A, Munger KL. Environmental risk factors for multiple sclerosis. Part I: the role of infection. Ann Neurol. 2007;61(4):288–99.
Baranzini SE, Srinivasan R, Khankhanian P, Okuda DT, Nelson SJ, Matthews PM, et al. Genetic variation influences glutamate concentrations in brains of patients with multiple sclerosis. Brain. 2010;133(9):2603–11. https://doi.org/10.1093/brain/awq192.
Strijbis EMM, Inkster B, Vounou M, Naegelin Y, Kappos L, Radue E-W, et al. Glutamate gene polymorphisms predict brain volumes in multiple sclerosis. Mult Scler. 2013;19(3):281–8. https://doi.org/10.1177/1352458512454345.
Wang JH, Pappas D, De Jager PL, Pelletier D, de Bakker PI, Kappos L, et al. Modeling the cumulative genetic risk for multiple sclerosis from genome-wide association data. Genome Med. 2011;3(1):3. https://doi.org/10.1186/gm217.
Keshet Y, Seger R. The MAP kinase signaling cascades: a system of hundreds of components regulates a diverse array of physiological functions. Methods Mol Biol. 2010;661:3–38. https://doi.org/10.1007/978-1-60761-795-2_1.
Chang L, Karin M. Mammalian MAP kinase signalling cascades. Nature. 2001;410(6824):37–40. https://doi.org/10.1038/35065000.
Shchetynsky K, Protsyuk D, Ronninger M, Diaz-Gallo L-M, Klareskog L, Padyukov L. Gene-gene interaction and RNA splicing profiles of MAP2K4 gene in rheumatoid arthritis. Clin Immunol. 2015;158(1):19–28.
Tuller T, Atar S, Ruppin E, Gurevich M, Achiron A. Common and specific signatures of gene expression and protein-protein interactions in autoimmune diseases. Genes Immun. 2013;14(2):67–82.
GJA t B, Bolk J, `t Hart BA, Laman JD. Multiple sclerosis is linked to MAPKERK overactivity in microglia. J Mol Med. 2021; Available from:. https://doi.org/10.1007/s00109-021-02080-4.
Mass E, Jacome-Galarza CE, Blank T, Lazarov T, Durham BH, Ozkaya N, et al. A somatic mutation in erythro-myeloid progenitors causes neurodegenerative disease. Nature. 2017;549(7672):389–93. https://doi.org/10.1038/nature23672.
Kotelnikova E, Kiani NA, Messinis D, Pertsovskaya I, Pliaka V, Bernardo-Faura M, et al. MAPK pathway and B cells overactivation in multiple sclerosis revealed by phosphoproteomics and genomic analysis. Proc Natl Acad Sci U S A. 2019;116(19):9671–6. https://doi.org/10.1073/pnas.1818347116.
Krementsov DN, Thornton TM, Teuscher C, Rincon M. The emerging role of p38 mitogen-activated protein kinase in multiple sclerosis and its models. Mol Cell Biol. 2013;33(19):3728–34.
Bernardo-Faura M, Rinas M, Wirbel J, Pertsovskaya I, Pliaka V, Messinis DE, Vila G, Sakellaropoulos T, Faigle W, Stridh P, Behrens JR. Prediction of combination therapies based on topological modeling of the immune signaling network in Multiple Sclerosis. Genome Medicine. 2021;13(1):1-6.
Maatz H, Jens M, Liss M, Schafer S, Heinig M, Kirchner M, et al. RNA-binding protein RBM20 represses splicing to orchestrate cardiac pre-mRNA processing. J Clin Invest. 2014;124(8):3419–30. https://doi.org/10.1172/JCI74523.
Sheikh F, Lyon RC, Chen J. Functions of myosin light chain-2 (MYL2) in cardiac muscle and disease. Gene. 2015;569(1):14–20.
Matyushenko AM, Levitsky DI. Molecular Mechanisms of Pathologies of Skeletal and Cardiac Muscles Caused by Point Mutations in the Tropomyosin Genes. Biochemistry. 2020;85(Suppl 1):S20–33.
Caleshu C, Sakhuja R, Nussbaum RL, Schiller NB, Ursell PC, Eng C, et al. Furthering the link between the sarcomere and primary cardiomyopathies: restrictive cardiomyopathy associated with multiple mutations in genes previously associated with hypertrophic or dilated cardiomyopathy. Am J Med Genet A. 2011;155A(9):2229–35. https://doi.org/10.1002/ajmg.a.34097.
Brody MJ, Hacker TA, Patel JR, Feng L, Sadoshima J, Tevosian SG, et al. Ablation of the cardiac-specific gene leucine-rich repeat containing 10 (Lrrc10) results in dilated cardiomyopathy. PLoS One. 2012;7(12):e51621.
Gupte TM, Haque F, Gangadharan B, Sunitha MS, Mukherjee S, Anandhan S, et al. Mechanistic Heterogeneity in Contractile Properties of α-Tropomyosin (TPM1) Mutants Associated with Inherited Cardiomyopathies*. J Biol Chem. 2015;290(11):7003–15.
Huang W, Liang J, Yuan C-C, Kazmierczak K, Zhou Z, Morales A, et al. Novel familial dilated cardiomyopathy mutation in MYL2 affects the structure and function of myosin regulatory light chain. FEBS J. 2015;282(12):2379–93.
Herman DS, Lam L, Taylor MRG, Wang L, Teekakirikul P, Christodoulou D, et al. Truncations of titin causing dilated cardiomyopathy. N Engl J Med. 2012;366(7):619–28. https://doi.org/10.1056/NEJMoa1110186.
Marston S, Montgiraud C, Munster AB, Copeland O, Choi O, Dos Remedios C, et al. OBSCN Mutations Associated with Dilated Cardiomyopathy and Haploinsufficiency. PLoS One. 2015;10(9):e0138568.
Marston S. Obscurin variants and inherited cardiomyopathies. Biophys Rev. 2017;9(3):239–43. https://doi.org/10.1007/s12551-017-0264-8.
McNally EM, Mestroni L. Dilated Cardiomyopathy: Genetic Determinants and Mechanisms. Circ Res. 2017;121(7):731–48. https://doi.org/10.1161/CIRCRESAHA.116.309396.
Boczek NJ, Ye D, Jin F, Tester DJ, Huseby A, Bos JM, et al. Identification and Functional Characterization of a Novel CACNA1C-Mediated Cardiac Disorder Characterized by Prolonged QT Intervals With Hypertrophic Cardiomyopathy, Congenital Heart Defects, and Sudden Cardiac Death. Circ Arrhythm Electrophysiol. 2015;8(5):1122–32.
Mouton J, Ronjat M, Jona I, Villaz M, Feltz A, Maulet Y. Skeletal and cardiac ryanodine receptors bind to the Ca(2+)-sensor region of dihydropyridine receptor alpha(1C) subunit. FEBS Lett. 2001;505(3):441–4. https://doi.org/10.1016/S0014-5793(01)02866-6.
Kanehisa M, Goto S. KEGG: kyoto encyclopedia of genes and genomes. Nucleic Acids Res. 2000;28(1):27–30.
Ather S, Respress JL, Li N, Wehrens XHT. Alterations in ryanodine receptors and related proteins in heart failure. Biochim Biophys Acta. 2013;1832(12):2425–31. https://doi.org/10.1016/j.bbadis.2013.06.008.
Yano M, Yamamoto T, Kobayashi S, Matsuzaki M. Role of ryanodine receptor as a Ca2+ regulatory center in normal and failing hearts. J Cardiol. 2009;53(1):1–7. https://doi.org/10.1016/j.jjcc.2008.10.008.
Moccia F, Lodola F, Stadiotti I, Pilato CA, Bellin M, Carugo S, et al. Calcium as a Key Player in Arrhythmogenic Cardiomyopathy: Adhesion Disorder or Intracellular Alteration? Int J Mol Sci. 2019;16(16):20(16). Available from:. https://doi.org/10.3390/ijms20163986.
Jaffrey SR, Wilkinson MF. Nonsense-mediated RNA decay in the brain: emerging modulator of neural development and disease. Nat Rev Neurosci. 2018;19(12):715–28.
Schwerk C, Schulze-Osthoff K. Regulation of apoptosis by alternative pre-mRNA splicing. Mol Cell. 2005;19(1):1–13.
List M, Alcaraz N, Dissing-Hansen M, Ditzel HJ, Mollenhauer J, Baumbach J. KeyPathwayMinerWeb: online multi-omics network enrichment. Nucleic Acids Res. 2016;44(W1):W98–104.
Kinsella RJ, Kähäri A, Haider S, Zamora J, Proctor G, Spudich G, et al. Ensembl BioMarts: a hub for data retrieval across taxonomic space. Database. 2011;2011:bar030.
Kuleshov MV, Jones MR, Rouillard AD, Fernandez NF, Duan Q, Wang Z, et al. Enrichr: a comprehensive gene set enrichment analysis web server 2016 update. Nucleic Acids Res. 2016;44(W1):W90–7. https://doi.org/10.1093/nar/gkw377.
Signorelli M, Vinciotti V, Wit EC. NEAT: an efficient network enrichment analysis test. BMC Bioinformatics. 2016;17(1):352. https://doi.org/10.1186/s12859-016-1203-6.
Benjamini Y, Hochberg Y. Controlling the False Discovery Rate: A Practical and Powerful Approach to Multiple Testing. J Royal Stat Soc: Ser B (Methodological). 1995;57:289–300. Available from:. https://doi.org/10.1111/j.2517-6161.1995.tb02031.x.
Irimia M, Weatheritt RJ, Ellis JD, Parikshak NN, Gonatopoulos-Pournatzis T, Babor M, et al. A highly conserved program of neuronal microexons is misregulated in autistic brains. Cell. 2014;159(7):1511–23. https://doi.org/10.1016/j.cell.2014.11.035.
Dobin A, Davis CA, Schlesinger F, Drenkow J, Zaleski C, Jha S, et al. STAR: ultrafast universal RNA-seq aligner. Bioinformatics. 2013;29(1):15–21.
Van Der Walt S, Chris Colbert S, Varoquaux G. The NumPy array: a structure for efficient numerical computation. arXiv [cs.MS]. 2011; Available from: http://arxiv.org/abs/1102.1523.
McKinney W, Others. Data structures for statistical computing in python. Proceedings of the 9th Python in Science Conference. 2010;445:51–6. https://doi.org/10.25080/Majora-92bf1922-00a.
Hagberg A, Swart P, S Chult D. Exploring network structure, dynamics, and function using networkx. Proceedings of the 7th Python in Science Conference (SciPy2008). 2008;11-15.
Virtanen P, Gommers R, Oliphant TE, Haberland M, Reddy T, Cournapeau D, et al. SciPy 1.0: fundamental algorithms for scientific computing in Python. Nat Methods. 2020;17(3):261–72. https://doi.org/10.1038/s41592-019-0686-2.
Seabold S, Perktold J. Statsmodels: Econometric and statistical modeling with python. Proceedings of the 9th Python in Science Conference. https://doi.org/10.25080/Majora-92bf1922-011.
Fruchterman TMJ, Reingold EM. Graph drawing by force-directed placement. Softw Pract Exp. 1991;21:1129–64. Available from:. https://doi.org/10.1002/spe.4380211102.
Shen S, Park JW, Lu Z-X, Lin L, Henry MD, Wu YN, et al. rMATS: robust and flexible detection of differential alternative splicing from replicate RNA-Seq data. Proc Natl Acad Sci U S A. 2014;111(51):E5593–601.
Sterne-Weiler T, Weatheritt RJ, Best A, Ha KCH. Whippet: an efficient method for the detection and quantification of 944 alternative splicing reveals extensive transcriptomic complexity. bioRxiv. 2017. https://doi.org/10.1101/158519
Louadi Z. NEASE: A network-based approach for the enrichment of alternative splicing events. Github. 2021. Available from: https://github.com/louadi/NEASE. Accessed 22 Nov 2021.
Louadi Z. NEASE: v.1.1.6. Zenodo; 2021. Available from:. https://doi.org/10.5281/zenodo.5653490.
Louadi Z. NEASE-tutorials: v1.2. Zenodo; 2021. Available from: https://doi.org/10.5281/ZENODO.5562626
Peer review information
Barbara Cheifet was the primary editor of this article and managed its editorial process and peer review in collaboration with the rest of the editorial team.
The review history is available as Additional file 8.
This work was supported by the German Federal Ministry of Education and Research (BMBF) within the framework of the e:Med research and funding concept [grant 01ZX1908A (Sys_CARE)]. JB was partially funded by his VILLUM Young Investigator Grant (nr.13154). MLE is grateful for financial support from Lundbeckfonden (no. R347-2020-2454). ZI is grateful for financial support from Scleroseforeningen (no. A29926, A 31829, A33600). Open Access funding enabled and organized by Projekt DEAL.
Ethics approval and consent to participate
Consent for publication
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Tables S1-S3. Table S1. Enrichment of the pathway “Muscle contraction” from Reactome for the exons upregulated in the muscles generated by the NEASE package. Table S2. Enrichment of the pathway “Synaptic vesicle cycle” from KEGG for the exons upregulated in the neural tissues generated by the NEASE package. Table S3. RNA-Seq samples used in the study. Figs. S1-S6. Fig. S1. Node degree distribution of the classic PPI and structurally annotated PPI, the latter contains only interactions with evidence from DDIs and DMIs or residue-level evidence from the co-resolved structure. Fig S2. Network Enrichment Analysis using EviNet webtool for exons upregulated in muscles and neural tissues. Fig. S3. NEASE visualization highlights the interactions of differentially spliced genes with the DCM pathway. Fig. S4. The PSI values of two exon skipping events in the genes TPM1 and DST from the GTEx dataset confirm that both the exons are upregulated in muscles and heart tissues. Fig. S5. The PSI values of two exon skipping events in the genes CLTA and CLTB from the GTEx dataset confirm that both the exons are upregulated in neural tissues. Fig. S6. Pseudocode of NEASE algorithm.
Tables S4-S7. The analysis of the upregulated exons in muscles. Table S4: List of spliced domains. Table S5: List of affected edges. Table S6: Gene-level enrichment. Table S7: NEASE enrichment.
Tables S8-S11. The analysis of the upregulated exons in neural. Table S8. List of spliced domains. Table S9. List of affected edges. Table S10. Gene-level enrichment. Table S11. NEASE enrichment.
Tables S12-S16. Differential splicing analysis between reticulated platelets and mature platelets. Table S12. List of spliced domains. Table S13. List of affected edges. Table S14. Gene-level enrichment. Table S15. NEASE enrichment. Table S16. Enrichment of the pathway “GPCR downstream signal”.
Tables S17-S21. Differential splicing analysis between normal-appearing white matter and acute lesion from multiple sclerosis patients. Table S17. List of spliced domains. Table S18. List of affected edges. Table S19. Gene-level enrichment. Table S20. NEASE enrichment. Table S21. Enrichment of the pathway “Neurotransmitter receptors and postsynaptic signal transmission”.
Tables S22-S24. Differential splicing analysis between Dilated Cardiomyopathy patients and controls. Table S22. List of spliced domains. Table S23. List of affected edges. Table S24. Gene-level enrichment. Table S25. NEASE enrichment. Table S26. Enrichment of the pathway “Dilated cardiomyopathy”.
Table S27. The PSI values for VastDB exons in different tissues.
About this article
Cite this article
Louadi, Z., Elkjaer, M.L., Klug, M. et al. Functional enrichment of alternative splicing events with NEASE reveals insights into tissue identity and diseases. Genome Biol 22, 327 (2021). https://doi.org/10.1186/s13059-021-02538-1
- Alternative splicing
- Differential splicing
- Functional enrichment
- Systems biology
- Protein-protein interactions
- Disease pathways
- Platelet activation
- Multiple sclerosis
- Dilated cardiomyopathy