Exon inclusion signatures enable accurate estimation of splicing factor activity

Domain: biology
Published: 2026-03-12
Editorial outcome: major_revision

Abstract
Splicing factors control exon inclusion in messenger RNAs, shaping transcriptome and proteome diversity. Their catalytic activity is regulated by multiple layers, making single-omic measurements on their own fall short in identifying which splicing factors underlie a phenotype. Here, we posit that splicing factor activity, defined as a splicing factor’s ability to modulate exon inclusion, can be estimated from changes in exon inclusion signatures. To test this hypothesis, we benchmark methods for constructing splicing factor → exon networks and estimating splicing factor activity. We find that combining RNA-seq perturbation-based networks with VIPER (Virtual Inference of Protein Activity by Enriched Regulon analysis) accurately captures splicing factor activity as modulated by multiple regulatory layers. This approach integrates splicing factor regulation into a single score derived solely from exon inclusion signatures, allowing functional interpretation of heterogeneous conditions. As a proof of concept, we identify recurrent cancer splicing programs, revealing associations with oncogenic- and tumor suppressor-like splicing factors missed by conventional methods. These programs correlate with patient survival and key cancer hallmarks: initiation, proliferation, and immune evasion. Altogether, we show splicing factor activity can be accurately estimated from exon inclusion changes, enabling comprehensive analyses of splicing regulation with minimal data requirements.

Introduction
Alternative splicing is regulated by splicing factors that selectively modulate intron removal from transcribed precursor messenger RNAs (mRNAs), leading to different patterns of exon inclusion. Together with alternative transcription start and termination sites, splicing factors determine each gene’s pool of mature mRNAs, impacting overall mRNA stability, translation, protein activation, and protein-protein interactions1. A splicing factor’s ability to modulate exon inclusion, its splicing factor activity, is controlled by multiple mechanisms, including regulation of its transcription and pre-mRNA splicing, protein synthesis and degradation, post-translational modifications, subcellular localization, and interactions with other proteins to form functional complexes (Fig. 1). Accurate assessment of splicing factor activity is thus challenging as it requires accounting for highly diverse regulatory mechanisms.

Splicing factor activities have been typically inferred from the presence of activating and inactivating somatic mutations, or differential gene expression2–5. However, these represent just two of the many possible mechanisms contributing to splicing factor activity modulation. This partial picture limits identifying aberrant splicing factor activity, especially in complex genetic contexts, such as cancer, where multiple molecular alterations contribute to cancer hallmarks4,6–9. A possible solution could be collecting more omic measurements on splicing factors. However, this would not only require large-scale efforts and resources but also methods to integrate these data and calculate splicing factor activities, which are still elusive.

The challenge of assessing protein activity, generally understood as the extent to which a protein performs its biological function under specific cellular conditions, has been effectively addressed for transcription factors. Rather than estimating transcription factor activity from their own omic profiles, the changes in the expression of their target genes more accurately capture the combined regulatory and signaling mechanisms that determine their activity10. This approach requires addressing two problems: first, assessing the transcriptional targets of every transcription factor, and second, estimating their activity based on the differential expression of these targets. Transcription factor → gene interactions (i.e., regulatory networks) can be inferred computationally from gene expression datasets by reverse engineering algorithms such as ARACNe11,12, as well as empirically, by mining curated databases generated from publicly available experiments13–15. Then, given a high-quality regulatory network that includes the probability (likelihood) and directionality (mode of regulation) of each interaction, transcription factor activity can be accurately assessed from expression signatures through bioinformatic methods like VIPER16 and others17, which statistically integrate regulatory interactions into a unified activity score. Although this approach only requires transcriptomic data, it has been proven highly valuable in disentangling the common functional states of transcriptional regulatory circuits in heterogeneous molecular contexts like cancer6.

Following the same reasoning, we posit that exon inclusion signatures report on the activity of splicing factors. If this is the case, frameworks like VIPER could be extended to estimate splicing factor activity from changes in exon inclusion across their target mRNAs, represented as networks of splicing factor → exon interactions (splicing factor networks, for simplicity) (Fig. 2a). To achieve this, we must determine which network-construction and protein activity-inference approaches are best suited to splicing factors. Once defined, these strategies are expected to reveal new insights from heterogeneous datasets, such as cancer samples, where multiple regulatory layers are dysregulated simultaneously.

Here, to repurpose protein activity analysis for splicing factors, we assess different methods for splicing factor network inference and activity estimation. We find that combining empirically derived splicing factor networks with the VIPER framework yields the most accurate activity estimates across our benchmark datasets. Notably, this approach remains effective in independent experiments that are more complex than those used to construct the networks. To illustrate its broader applicability in the biomedical context, we apply the framework to cancer-related phenotypes. This uncovers two biologically relevant and complementary splicing programs that are not detected by standard splicing factor expression analysis. Integration with transcriptomic, proteomic, and phosphoproteomic data suggests that these cancer-associated splicing programs are coordinately regulated through diverse molecular alterations spanning multiple regulatory layers. Altogether, splicing factor activity analysis is a powerful method for generating testable hypotheses about the splicing factors that drive complex cellular phenotypes.

Results
Evaluating splicing factor activity estimation in single-perturbation experiments
As illustrated in Fig. 2a, estimating splicing factor activity with an adaptation of the VIPER framework requires: (i) a network linking each splicing factor to its target exons; and (ii) an enrichment-based activity estimator that quantifies how strongly each factor contributes to the observed delta PSI pattern. Although this framework performs well for transcription factors, its suitability for splicing factors has not been evaluated. We therefore benchmarked both the network construction strategies and the activity estimation methods.

To assess the accuracy of splicing factor activity estimation methods, we created a set of five benchmark datasets based on splicing factor perturbation experiments. This is complicated because some proteins like CPSF6 (polyadenylation machinery) and DNMT1 (methyltransferase) have been reported as splicing factors, although they perform functions that indirectly can influence splicing. As a result, there is no universally accepted definition of a splicing factor, and published lists show limited overlap (Supplementary Fig. 1a). Thus, we opted for an inclusive approach that resulted in a consensus list of 509 splicing factors from landmark studies (Supplementary Data 1 and 2)8,18–21.

We then searched for RNA-seq studies perturbing any of these factors and identified 474 bulk RNA-seq experiments covering 248 distinct splicing factors. In this study, we focused on cassette exons, which provide a robust and widely used readout of splicing regulation. Other splicing event types were collected but not analyzed further (see Data Availability and Discussion). For each experiment, we used vast-tools22 to quantify exon inclusion (PSI, percentage spliced in) across conditions and computed exon inclusion signatures as delta PSI by subtracting control PSI values from those in the perturbed condition. These uniformly processed datasets enabled us to evaluate the accuracy of different splicing factor activity estimation methods at recovering the perturbed splicing factor solely from exon inclusion changes (Supplementary Note 1; Supplementary Fig. 1b).

For the first VIPER framework requirement, we compared computational and empirical network reconstruction approaches. Computational networks were inferred from large RNA-seq datasets using ARACNe (pairwise associations between splicing factor expression and exon inclusion) and multivariate regression (combined linear associations between the expression of all splicing factors and PSI of each exon). Empirical networks were derived directly from the uniformly processed perturbation datasets by selecting exons with a threshold of 15% absolute delta PSI change upon perturbation for each splicing factor. This threshold was chosen after testing performance across a range of values (Supplementary Fig. 1c; see below definition of performance metrics). To ensure a fair comparison, we excluded any activity estimates produced from the same dataset used to construct a given empirical network (see “Methods”).

VIPER
Then, we evaluated four approaches to estimate splicing factor activity: gene set enrichment analysis (GSEA), Spearman and Pearson correlations, and VIPER. These methods interpret networks differently: GSEA treats targets as an unweighted set; correlations incorporate interaction likelihood and direction; and VIPER combines both using a rank-based enrichment analysis (aREA).

For each combination of network type and activity estimation method (producing a total of 8 configurations), we calculated splicing factor activity in benchmark datasets, producing matrices of splicing factors by perturbation experiments (Fig. 2b). We quantified how well activities identified the perturbed splicing factor in each experiment using ROC-AUCs. In this framework, the perturbed splicing factor in each experiment is considered a true positive, and all others are false positives. We assessed performance in two complementary ways: across splicing factors (how well predicted activity identifies which splicing factor was perturbed in each experiment) and across perturbations (how well predicted activity identifies the experiments that perturbed each splicing factor) (Supplementary Note 1).

Empirical networks combined with VIPER consistently achieved the highest performance, outperforming computational networks and other activity estimation methods (Fig. 2c; and Supplementary Note 2). Additionally, we found that the poor performance of computational networks is largely due to their inaccurate prediction of regulatory direction, which limits their utility for activity estimation (Supplementary Fig. 1d–f; and Supplementary Note 3).

Altogether, these results indicate that interpreting exon inclusion signatures with VIPER using networks derived from perturbation data estimates the most accurate and informative splicing factor activities.

Empirical networks and VIPER enable robust splicing factor activity estimation
To assess the robustness of splicing factor activity estimation methods, we evaluated performance across the diverse studies and experimental contexts represented in our benchmark datasets. Given the poor performance of computational networks, we restricted the following analyses to empirical networks. Here, we provide a high-level overview of these assessments; full details are presented in Supplementary Note 4.

Despite the well-documented context-specificity of certain splicing changes, we found that empirical networks combined with VIPER consistently achieve high accuracy (ROC-AUC ~0.9), indicating strong generalizability (Supplementary Fig. 2a, b). To better understand the basis of this robustness (i.e., the stability, reliability, and consistency of the method across diverse conditions), we examined the structure of empirical networks. As described above, an exon is assigned as a target of a splicing factor if, upon perturbation of that splicing factor, its inclusion changes by more than 15%. Most exons are regulated by only a few splicing factors (median = 2), indicating that splicing factor → exon interactions in empirical networks are generally specific, even though a minority of exons (6 out of 22,049) can be influenced by many factors (up to 90 or more; Fig. 2d). In line with this, reducing the number of interactions had little effect on the prediction performance, showing that only a small subset of interactions is sufficient for reliable activity estimation (Supplementary Fig. 2c).

We also assessed how robustness varies by splicing factor class. Core spliceosomal factors, those directly involved in the splicing reaction, tend to be more highly connected (Fig. 2e). Grouping performances by factor class revealed that activity estimates for core splicing factors achieve higher accuracy than those for other types such as RNA-binding proteins (RBPs) (Fig. 2f). This likely reflects that the functional targets of core factors are more stable across conditions, whereas RBP targets are more context-dependent and may therefore vary across datasets, reducing inference accuracy (median ROC-AUC: core = 0.97; RBP = 0.92).

Of note, while most splicing factors (i.e., those perturbed in more than one dataset) achieved ROC-AUC values > 0.5, in certain cases the perturbed factor showed AUC values below 0.5 (e.g., Supplementary Fig. 2a). This may indicate that, in those conditions, activities estimated from the held-out dataset are inconsistent with the expected effect. Biologically, this could be due to co-regulation: perturbing one splicing factor may elicit stronger effects on upstream regulators or downstream effectors. Technically, evaluating across datasets could introduce batch effects or experimental noise. These single-factor perturbation evaluations, therefore, need to be complemented by validations using orthogonal perturbation approaches, as we describe below.

As we did not find independent replicates for all perturbed splicing factors, we chose to retain all available factors in our evaluations, even when quality metrics in this particular assessment may be suboptimal, to provide a broad, comparative view of method performance, highlight incomplete networks or pathways, and enable hypothesis generation in less-characterized regions of the splicing machinery.

Finally, to further confirm the consistency of our findings, we repeated all evaluation analyses using an alternative metric: the rank percentile of the perturbed factor among all factors. These rank percentiles showed strong agreement with the ROC-AUC-based results (Supplementary Fig. 3), reinforcing the robustness of our conclusions.

Empirical networks capture functional regulation
Splicing factors influence exon inclusion through diverse mechanisms, including RNA binding, modulation of spliceosome composition, and chromatin- or transcription-associated processes8. Empirical splicing factor networks provide a unified functional readout of these mechanisms by assigning exons whose inclusion changes upon perturbation of each factor.

To compare this functional definition with an alternative based on direct RNA binding, we built CLIP-based (crosslinking and immunoprecipitation data) networks for splicing factors23 (Supplementary Fig. 2d). For factors with both data types, most empirical interactions were not detected by CLIP, indicating that CLIP peaks and perturbation-responsive exons capture complementary aspects of splicing regulation (Supplementary Fig. 2e). Overlapping interactions performed well in activity inference (Supplementary Fig. 2f). VIPER also maintained high performance even when using only empirical interactions lacking CLIP support (Supplementary Fig. 2f). Full details on CLIP network construction and analysis are provided in Supplementary Note 5.

These results suggest that empirical networks derived from exon inclusion changes robustly capture functional regulatory relationships, regardless of their mechanistic origin, and provide a reliable basis for estimating splicing factor activity.

Estimated splicing factor activity recapitulates regulatory mechanisms of increasing complexity
Because empirical networks combined with VIPER reliably recover the perturbed factor across diverse evaluations, we henceforth adopt this approach for estimating splicing factor activity.

We hypothesized that this strategy could extend beyond single-factor genetic perturbations and capture more complex, multi-omic regulatory mechanisms that ultimately influence splicing factor activity solely through exon inclusion signature analysis (Fig. 1). To illustrate this, we evaluated splicing factor activity analysis across several experimental validations involving perturbations of diverse regulatory mechanisms.

For example, the drug Indisulam downregulates the protein level of splicing factor RBM39 by mediating its proteasomal degradation (Fig. 3a, b)24–27. In contrast, the mRNA levels of its encoding gene increase consistently upon treatment, probably as a compensatory mechanism (Supplementary Fig. 4a and b)24,25. In this case, using changes in splicing factor mRNA levels to evaluate splicing factor activity would lead to the incorrect conclusion that Indisulam increases RBM39’s activity. In sharp contrast, splicing factor activity analysis of transcriptional profiles following drug perturbation identifies RBM39 as the splicing factor with the greatest decrease in activity, thus effectively recapitulating Indisulam’s mechanism of action (Fig. 3c). As a result, splicing factor activity analysis recapitulated protein-level effects, despite using only transcriptomic data.

Splicing factors are often modulated in concert, producing similar exon inclusion signatures. To test whether our approach can still identify perturbed regulators in this regime, we applied activity analysis to our studies in which two to four splicing factors were silenced simultaneously28–31 (Fig. 3d). In all cases, the perturbed factors ranked among the most inactivated, indicating that VIPER can detect multiple regulators contributing to a combined signature.

Because the combinatorial perturbation studies targeted functionally related splicing factors, some unperturbed paralogs or interactors to the perturbed splicing factors showed low activity despite not being perturbed (e.g., CLK1 and CLK2 when CLK3 and CLK4 were perturbed) (Fig. 3d). To test whether this reflects true co-perturbation or simply similar exon-inclusion signatures, we examined the corresponding single-splicing factor perturbations (Supplementary Fig. 4c). Although activity could not be computed for all splicing factors due to limited coverage, the available cases showed that, except for the NCBP1-ZCCHC8 pair, splicing factors with strong activity in combinatorial experiments did not show strong inactivation in the matching single perturbations. This suggests that splicing factor activities can identify co-perturbed splicing factors in pairwise perturbations, but in higher-order perturbations, functionally related factors may still appear strongly inactivated even when they are not directly perturbed.

Splicing factors interact to form a large and dynamic protein complex called the spliceosome21. We thus hypothesized that the VIPER-assessed activity of splicing factors interacting in the spliceosome complex via protein-protein interactions (PPIs) should be correlated. At the beginning of the splicing reaction, during splice site recognition, the spliceosome forms the E, A, and B complexes comprising small ribonucleoprotein (snRNP) complexes (U1, U2, U4, U5, and U6)8. Within the U2 snRNP, the SF3b complex (SF3B1, SF3B2, SF3B3, SF3B4, SF3B5, SF3B6, PHF5A) plays a pivotal role in branch point recognition32. Different splicing drugs perturb this step of the reaction. For instance, H3B-8800 and Pladienolide B bind SF3B133,34; E7107 prevents the binding of U2 snRNP to RNA35; and Isoginkgetin and Spliceostatin A prevent the transition from A to B spliceosome complexes in the splicing reaction36,37. It is thus reasonable to expect that perturbing the SF3b complex with these drugs should strongly decrease the activity of the associated splicing factors, proportionally to their distance to the SF3b complex in the PPI network (Fig. 3e, f).

Consistent with this prediction, treating cell lines with the aforementioned drugs induced notable VIPER-assessed inactivation of splicing factors in the SF3b complex, followed by splicing factors in the U2 snRNP and others splicing factors with direct (shortest path length in PPI space = 1) and indirect (shortest path length in PPI space > 1) interactions with the SF3b complex (Fig. 3g).

Taken together, these data further demonstrate that combining VIPER with empirical splicing factor networks accurately recapitulates splicing factor activity modulation mechanisms across a range of diverse perturbations that were not used to generate the splicing factor networks.

Recurrent aberrant splicing factor activity defines two cancer splicing programs
Cancer co-opts the splicing machinery in favor of tumor progression4,9,38. However, few studies have examined whether these splicing programs are unique to specific cancer cohorts or are shared, with most research focusing on single regulatory layers of splicing factors39,40.

Our approach captures various regulatory mechanisms driving splicing factor activity. This suggests that estimated splicing factor activities reflect the combined output of all regulatory inputs affecting their function. To demonstrate its utility in complex cellular contexts, we investigated whether splicing factor dysregulation in cancer follows a general pattern or varies by cohort. We analyzed 33 cancer cohorts to determine if certain splicing factors are consistently dysregulated, suggesting a recurrent cancer splicing program.

Of the 33 cancer cohorts included in The Cancer Genome Atlas (TCGA), 14 include ≥ 10 exon-level transcriptional profiles from primary tumor (PT) and solid tissue normal (STN) biopsies. We leveraged these samples to assess differential splicing factor activity by comparing the two conditions (Fig. 4a, b; see Methods). We then focused on the splicing factors that were significantly differentially active (FDR < 0.05, two-sided Wilcoxon rank sum test) in PTs compared with their corresponding STNs across all cancer cohorts (Fig. 4b). Interestingly, some splicing factors emerged as recurrently activated or inactivated across a majority of cancer cohorts (Fig. 4c; and Supplementary Fig. 5a). This suggests that oncogenic alterations may cause coordinated and recurrent dysregulation of splicing factor activity, despite their well-described heterogeneity. Considering a threshold of at least 5 out of 14 cancer cohorts to be indicative of recurrent dysregulation (see “Methods”), the analysis identified a set of (n=61) putatively oncogenic splicing factors (i.e., aberrantly activated in ≥ 5 cohorts), and a set of (n=61) putatively tumor suppressor splicing factors (i.e., aberrantly inactivated in ≥ 5 cohorts) (Supplementary Data 3). Of note, PTBP1, a splicing factor with confirmed oncogenic roles in four different cancer types41, emerged as the most frequently activated splicing factor.

If these recurrently dysregulated splicing factors play a role in tumorigenesis and progression, then activation of oncogenic and tumor suppressor factors should be predictive of poor and good prognosis, respectively. Indeed, patient outcome analyses confirmed these expectations (Fig. 4d; see “Methods”). For example, analysis of primary tumor samples across 20 TCGA cohorts with sufficient metadata (see “Methods”) revealed that aberrant activity of a recurrent oncogenic splicing factor (CPSF6) correlates with poor prognosis in 18 cancer cohorts (Cox Proportional Hazards (PH) coefficient > 0) while the activity of a recurrent tumor suppressor splicing factor (DNMT1) correlates with favorable outcome in 16 cancer cohorts (Cox PH coefficient < 0) (Fig. 4e).

Extending the analysis to all 122 recurrently dysregulated splicing factors showed the same trends. Indeed, aberrant activation of oncogenic splicing factors correlates with poor outcomes across significantly more cancer cohorts than tumor suppressor ones (p = 1.3·10-12, two-sided Wilcoxon rank sum test). The opposite happens for tumor suppressor-like splicing factors; their activity correlates with good outcome across significantly more cancer cohorts than the oncogenic ones (p = 8.6·10-13, two-sided Wilcoxon rank sum test) (Fig. 4f). These patterns are conserved if we consider different thresholds to define cancer splicing program recurrence (Supplementary Fig. 5b).

Typically, splicing factor activity is extrapolated from expression. We argue that splicing factor expression is only one of the molecular mechanisms regulating splicing factor activity (Fig. 1). Hence, we expect that our approach should capture activity differences that are not necessarily reflected by differences in mRNA levels. To evaluate this, we repeated the same analysis with splicing factor gene expression. Confirming our expectations, gene expression fold-changes in tumor vs matched normal samples did not correlate with their differential activity (Spearman correlation = 0.047, p = 0.042) (Supplementary Fig. 5c). Candidate oncogenic (n=106) and tumor suppressor splicing factors (n=30) identified by gene expression analysis (Supplementary Fig. 5d) were not predictive of poor and good outcome (p > 0.05 in both cases) (Supplementary Fig. 5e and f). Thus, by focusing on changes in splicing factor targets rather than in splicing factor expression, our approach systematically elucidated the potential contribution of their dysregulation to tumorigenesis and outcome.

To further characterize these cancer splicing programs, we examined their putative roles in terms of associated splicing reactions and in RBP families using SpliceosomeDB42,43 and RBPDB42, respectively. Out of the 122 splicing factors involved in cancer splicing programs, 32 and 54 splicing factors were mapped to each database, respectively. Most were part of the C complex in the splicing reaction – involved in the catalytic steps –, and belonged to the RRMx1 and RRMx2 families of RBP domains (Fig. 4g and h).

Linking cancer splicing program dysregulation to cancer hallmarks
Given the functional relevance of the cancer splicing programs identified through our approach, we wondered which biological processes they may potentially regulate (Fig. 5a). We performed an overrepresentation analysis (ORA) of the genes corresponding to the exons targeted by either oncogenic-like or tumor suppressor-like splicing factors.

To investigate the potential involvement of cancer splicing programs in cancer hallmarks, we performed an over-representation analysis (ORA) on genes corresponding to exons targeted by either oncogenic-like or tumor suppressor-like splicing factors. We selected the top 10 enriched terms (FDR < 0.05) with the most target genes for each program, highlighting both shared and unique terms, while ensuring that all significantly enriched terms from each program were included (Fig. 5b; and Supplementary Data 4). Two of the enriched terms were related to antigen presentation, a critically regulated process in tumor immunology, and six others were related to cell cycle.

Immune evasion hampers the effectiveness of immune checkpoint blockade (ICB) therapies in cancer patients, complicating treatment response predictions44. A total of 181 genes belonging to the immune evasion terms “Antigen Processing Ubiquitination Proteasome Degradation” and “Class I MHC-mediated Antigen Processing Presentation” contain exons regulated by the tumor suppressor-like splicing program (Fig. 5b). To prioritize tumor suppressor-associated exons potentially involved in ICB efficacy, we analyzed exon inclusion differences between ICB-responding and non-responding patients based on RECIST criteria45, using data from Riaz et al.46. Out of 608 exons from genes in immune-related terms, 16 were differentially spliced (p < 0.05, two-sided Wilcoxon rank sum tests). For additional prioritization, we used data from an in vivo genome-wide CRISPR screen comparing ICB response tumors in immunocompetent versus immuno-deficient NOD SCID Il2rg−/− control mice46, identifying exons within genes affecting ICB response (Fig. 5c). Specifically, exon HsaEX1036341 in SEC22B – one of the strongest immune evasion inhibitors identified (impairment score = −1.2) – emerged as a promising candidate. Exclusion of this exon disrupts SEC22B’s open reading frame, potentially reducing functional protein and improving ICB response, as evidenced by its lower inclusion in responder (PSI = 95.9, n = 10) than in non-responder (PSI = 99.4, n = 22) patients (p = 0.024, two-sided Wilcoxon rank sum tests). Survival analysis further supported this finding, with low exon inclusion in ICB-treated patients correlating with better prognosis (p = 0.0079, two-sided Log-rank test; best cutpoint = 98.91 PSI) (Fig. 5d). Examining empirical networks, we identified three splicing factors – PRPF8, SRSF1, and SF3B3 – as regulators of HsaEX1036341. An ORA of their targets revealed enrichment (FDR < 0.05) in cell cycle control pathways and “Class I MHC-mediated Antigen Processing Presentation” (Supplementary Fig. 6a), suggesting that inactivating these splicing factors should impair immune responses as well. In vivo CRISPR data revealed that silencing the mouse ortholog of SF3B3 closely mirrors SEC22B silencing (impairment score = −0.89) compared to silencing SRSF1 and PRPF8 (impairment scores = −0.61 and 0.68, respectively). This suggests SF3B3 as the most likely regulator of HsaEX1036341 in SEC22B. We conclude that analyzing the interactions within the identified cancer splicing programs could unveil exons with potentially relevant roles in ICB efficacy.

To validate the involvement of cancer splicing programs in cell proliferation, we estimated splicing factor activities for 1,015 cell lines in the Cancer Cell Line Encyclopedia (CCLE) and correlated the median activities among either oncogenic-like or tumor suppressor-like splicing factors against MKI67 expression, a marker of cell proliferation47. The activity of the oncogenic-like splicing program correlates positively with the proliferation marker (Pearson correlation = 0.2, p < 6.9·10−11) and the activity of the tumor suppressor-like splicing program correlates negatively (Pearson correlation = −0.32, p < 2.2·10−16) (Fig. 5e). Further, the sign of these relationships is reproducible when computing the correlation by grouping CCLE cell lines by cancer type of origin (Supplementary Fig. 6b). This not only confirms the link between cell proliferative state and the activation of cancer splicing programs but also corroborates the cancer driver-like behavior of the splicing programs identified in an independent dataset as high proliferation marker expression is simultaneously associated with activation of the oncogenic program and inactivation of the tumor suppressor program.

Taken together, we confirmed the biological relevance of the cancer splicing programs identified by splicing factor analyses exploring their associations with two hallmarks of cancer: immune evasion and aberrant proliferation.

Methods

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "In this paper, Anglada-Girotto et al. develop a method to infer splicing factor (SF) activity based on their targets. They do this by repurposing VIPER, an algorithm previously used to infer transcription factor activity from their targets’ expression. They use RBP perturbation data from ENCODE, as well as diverse datasets obtained from ENA to build a model of the SF regulatory network, which was used as an input to VIPER to make predictions on SF activity. While network reconstruction with ARCANe and multiple linear regression were ineffective, the authors were able to use a network defined by which exons were differentially spliced following SF perturbation. They validate this algorithm by predicting the effects of an anti-RBM39 drug, the samples with multiple perturbed SFs, and the similarity of interacting spliceosomal proteins. Finally, the authors apply this method to gain insight into SFs' contribution to cancer, identifying SFs with altered activity between cancer and normal tissue in multiple TCGA cohorts. They find these are correlated with prognosis and cell proliferation, and may be relevant for specific cancer pathways including immune evasion in the context of immune checkpoint blockade. Overall, the issue of understanding how SF networks regulate complex splicing programs is an important and interesting topic. This is especially relevant in the context of cancer, where splicing dysregulation is common, yet mechanisms are often poorly understood. We appreciate the authors’ approach to addressing this topic and their use of diverse datasets to do so. However we have major concerns with the methods and results, as presented here. Major concerns: ============ One major issue we found with the paper is lack of basic clarity of the methods description. The paper is simply not self contained. The main text repeats “See Methods”, coupling those with cryptic statements/descriptions. When we go to the Methods section, things do not improve. This leads to a distinct feeling of disregard or neglect which one hopes stayed at the level of writing and did not infiltrate the science itself. Let us give a few examples to make this criticism more concrete: The function of the VIPER algorithm is not clearly described, in the main text, or in the methods, including relevant details such as how it handles missing data. On page 12 the authors describe using ARACNe but adjusted with code from the Zhang lab. That code refers to a pre print - this means they are using unpublished algorithms with no assessment of its accuracy to derive their results, which is quite troubling. They then continue to state that they “ran the standard pipeline with a mutual information of …. Using table….” then state “...estimate the mode of regulation using viper::aracne2regulation”. Simply put: Why on earth would you think this makes sense to anyone who does not already know what is going on? Submitting a manuscript with this lack of clarity on what’s being done is concerning. But there are more severe concerns with how the method is executed and the choices of design and parameters along the way. For example, the authors state they used KNNImputer with K=5 to impute missing values since their method can not handle missing values. Why would that make any sense? What does the method do to the data? If some values are close to 0 and some are close to 1 (as sometimes happens in splicing) would you impute 0.5 (which is rare)? On page 13 the authors state “threshold of 15% leads to highly accurate and reasonably-sized networks” - based on what? What are the metrics? Where are the evaluations? Then they state “...the absolute delta PSI as the likelihood of the interaction” - why is that (a basic assumption in all the analysis) reasonable? I would understand if you are unsure if the change is small that it’s real but for example a change of 0.5 and a change of 0.25 in dPSI are both pretty large/confident. The authors also state “To keep consistency between all types of perturbations, we multiplied the delta PSI in overexpression experiments by -1” - this is another major assumption that is concerning. Clearly if inclusion level is say 0.1 and over expression leads to 0.4 this manipulation will give PSI of -0.3 violating the basic assumptions of PSI in 0-1 range. None of this is discussed, not to mention evaluated in the manuscript. Other examples: What are “extreme values” in page 15? What do you mean “within” and “across” (same paragraph p.15)? Additionally, differential gene expression analysis should be performed using specifically developed tools such as DESeq2. All of the above makes this reviewer uncertain about the validation of the analysis, making it difficult to understand the rationale for using this algorithm, or the value it adds to the analysis. To this point, it would also be helpful to evaluate an alternative approach or baseline method of some type for comparison, to demonstrate the utility of VIPER. For example, how do the results from VIPER compare to just looking at the overlap of differentially spliced exons from ENCODE or other perturbation datasets, or using CLIP-seq/motif based inference of regulation? The methods also mention that some runs were done using metaVIPER, a newer version of the algorithm, but the function of this algorithm is not described either. There is a similar issue for the ARCANE algorithm, as it is difficult for the reader to understand the function of this algorithm, or the rationale for using it, based on the main text and methods of this manuscript. One major issue with the paper is how datasets are integrated. We appreciate the authors integration of diverse RBP perturbation datasets, but batch effects are known to be a significant issue for ENCORE RBP perturbation datasets produced at different times/locations. Further, the authors are integrating many highly heterogeneous datasets obtained from ENA, which are likely to have significant batch effects. Additionally, based on the description in the methods, it is difficult for the reader to understand how the ENA datasets were grouped in constructing the empirical network. In general, the issue of batch effects seems to go completely ignored in this work. Minor comments: ============== Given the authors’ finding of correlation between SF activity and proliferation in CCLE cell lines, it may be worthwhile to analyze DepMap CRISPR screening data to check how these cell lines depend on the relevant SFs for proliferation. Due to the context-dependent effects of many SFs, the authors should address how well their model generalizes between different cell types and conditions. This could provide better justification for the reliability of the model’s predictions on diverse cancer datasets. It would be helpful to show which SF(s) may be responsible for regulating the SEC22B exon which may affect ICB response, and the specific pathways enriched in the targets of each SF. This would be very helpful to demonstrate the method’s ability to understand the specific mechanisms of regulation of cancer-related exons. It would also be helpful to know whether these SFs are themselves implicated in regulating ICB response by the in vivo CRISPR screen. In general, it would be helpful to clarify whether oncogenic/tumor suppressor-like SFs seem to have similar effects or if they have distinct phenotypes. All oncogenic or tumor suppressor-like SFs are lumped together for pathway analysis, when looking at correlation with CCLE proliferation, and when looking at carcinogenesis. How do different oncogenic/tumor suppressor-like SFs affect different cancer phenotypes, or different cancer types? At a high level, we were unsure given all of the above regarding the overall importance/impact. The computational analysis does not seem to offer a new method but rather adaptation of existing ones for the specific task at hand with a variety of parameter/design choices along the way which are unclear. It does not appear generalizable/usable for other researchers but rather a way to perform specific analysis. As such it hinges on the importance and validity of the derived conclusions/results, some of which look promising. In summary, we appreciate this new approach to an important and challenging topic, and the significant work and integration of diverse datasets that went into this paper. However, there are serious issues with the methods and the description of the methods.", "reviewer_2": "Summary In the manuscript “Disentangling the splicing factor programs underlying complex molecular phenotypes,” Anglada-Girotto et al. create a new database of splice-factor regulated exons and combine it with the already existing VIPER tool to estimate splice-factor activity from RNA-seq data. After some validation, the authors apply the data to estimate splice factor activity in several cancer datasets accompanied by selected follow-up analysis. While the database+VIPER approach is ingenious, and the datasets analyzed many, the manuscript in its current form, however, falls short. A substantial concern is that the validation needs to be extended, especially regarding how well it works in cell lines/tissues that are (very) different from the cell line in the database they create - meaning how well the analysis approach generalizes and thereby how trustworthy the cancer results are. In our opinion, the paper would also benefit from being rethought/re-structured. Is it a paper about a new and useful bioinformatic tool or a paper about understanding cancer better? Currently, the paper falls between these two chairs as it does neither sufficiently. To be successful, the tool angle would need more validation – and a more extensive database as only profiling ~1/2 of known splice factors does not feel comprehensive. As a cancer paper, more novel biological insights are needed since it is already very well-established that splicing, and its regulation, are dysregulated in cancer. Major comments • Validation of the splice factor network is too weak: 1. Regarding the splicing factor-> exon networks, please also provide other metrics than recall, such as false discovery rate. Similarly performance on the “completely independent perturbational signatures” should also be shown. 2. The splice factors profiled more than once (e.g. in two different cell lines) should have multiple observations. This also means they can be used to calculate ROC and precision-recall curves (that can be summarized as area under curve) 3. To show generalization, performance also needs to be assessed when the data evaluated originates from cell types/tissues that are (very) different from the ones used to create the SF->exon database. Since the database is (mainly) based on cell-line data, evaluation should also be on complex tissues. Something like evaluating the ENA-based dataset when using the ENCODE data as a database? • More work should be put into doing the computational deconvolution of exon-splice factors as that would increase the usefulness (of the analysis approach tremendously. Profiling ~1/2 of known splice factors does not feel particularly satisfactory or comprehensive for a new tool. What would happen if ARACNe (or similar) was run on e.g. all the TCGA data the authors already quantified? • If this is to be a cancer paper, more novel biological insights are needed since it is already very well-established that splicing, and its regulation, are dysregulated in cancer. Here a low-hanging fruit could be to create an additional database this time with exons affected by transcription factors knock-out and also apply that to the cancer datasets. • The title does not adequately reflect the study done and must be redone. Both title and abstract is to vague. After reading title and abstract neither of us understood what the study was about and what problem was solved. • While the method section is very detailed it needs to be more pedagogical – it is very hard to recreate what was done. • The Github/Code needs extensive updates about installation instructions, dependencies, and environmental requirements. Currently, it is hours of work of trial and error to get it running. Minor Comments • Abstract is sub-par it is a bit to vague in both what the problem is, how this paper helps and what the impact of that is. The impact/perspective is also missing from the discussion • There is a tendency towards using very long sentences reducing redability • Introduction fails to mention how alternative start sites and alternative termination sites are part of isoform generation • Discussion of how splicing-factor→pre-mRNA exon interaction networks are likely cell/tissue/disease specific is needed • Why are there 4 ENA based networks? And how were these 4 defined/sperated • Why is the ARACNe algorithm run on such a small dataset? Would increasing the number and diversity of samples not increase power? Something like running it on GTEx. • How does the approach work with existing databases like Polvèche et al 2023 (doi.org/10.1093/database/baad091) • The introduction should mention how the importance of splicing in cancer is already well established • Survival analysis: 1. How were the factors in Figure 3E chosen? They are not among the genes highlighted in 3C. How many of the “cancer-cohorts” are significant after correcting for multiple testing? I also assuming you tested all factors instead of only the two factors shown here so how do you correct for multiple testing there? 2. To carry any validity surivival analys needs to be corrected for cofactors such as gender age and tumor stage. • The ORA analysis in the “Cancer splicing programs regulate cancer hallmarks” section is highly over-interpreted. A statement like “exons regulated by the oncogenic-like splicing program are not enriched in immune evasion-related terms” does not hold much validity when only the top 10 gene-sets are considered since gene-set ranking is heavily biased by e.g. size of gene-set. In other words gene-set ranking is not the same as biological importance (or completeness) and hence cannot be interpreted this way. Furthermore the difference between the enrichment scores of the “ongogenic” and “tumorsupressor” in figure 4B is very small for most cases where they are both represented? For the ORA what gene-background is used and why? Carefull argumentation is needed as that is essential to the interpretation of the ORA • What gene-annotation (or rather exon annotation) database were quantified? Gencode? Ensemble? Refseq? And how are overlapping but non-identical exons handled (caused by e,g. an alternative donor site?) • Create a detailed guide that walks users through the analysis process. Explain how to execute each workflow, including any prerequisites or configuration steps. Provide examples of command-line instructions for running the Snakemake workflows, as well as expected outputs. • Watch out for causation vs correlation in the text We look forward to read and updated manuscript. Dimitrios S. Kanakoglou PhD Fellow Section for Bioinformatics, Health Technology Technical University of Denmark (DTU) Denmark And Kristoffer Vitting-Seerup, PhD Associate Professor Section for Bioinformatics, Health Technology Technical University of Denmark (DTU) Denmark", "reviewer_3": "In their manuscript ‘Disentangling the splicing factor programs underlying complex molecular phenotypes’ Serrano and colleagues aim to use a bioinformatics approach to characterize alternative splicing programs that are controlled by individual splicing factors. The authors use an extensive set (hundreds) of published RNA-Seq datasets with altered splicing factor expression, e.g. knockdown, knockout or overexpression of specific splicing regulatory proteins, to characterize alternative splicing signatures of splicing factors. In contrast to some previous work they do not rely on alteration of the mRNA of a given splicing factor, which can be misleading due to posttranscriptional/posttranslational regulation of their activity, but derive splicing factor activity from previously published datasets with known splicing factor alterations. They use a training dataset to define splicing patterns caused by alteration of a given splicing factor and then use this to infer misregulation of splicing factors based on global splicing patterns in disease conditions, e.g. different types of cancer. This is an important question and an interesting approach that has been taken for transcription factors before. However, there is a fundamental problem with applying this concept to the regulation caused by a specific splicing factor: There is an extensive amount of cross regulation between splicing factors as has been documented over many years by many labs. That means that altering the activity of one splicing factor (especially by artificial knockdown or overexpression approach as used for the training dataset) will lead to the alteration of the expression and/or activity of (many) other splicing factors. Given this extensive network of splicing factor cross regulation, the approach taken here will not yield meaningful results for individual splicing factors as it is never clear whether an effect is direct or indirect. For such an approach the training dataset would have to rely on other methods, e.g. acute depletion of splicing factors using inducible degrons or knockdown-rescue approaches and the combination with CLIP data to identify high confidence direct targets. With the current datasets and approach the results are too diffuse and potentially nonspecific. As such, the findings do not provide a sufficient advance for publication at this level.", "reviewer_4": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts.", "reviewer_5": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts."}, "author_response": "We thank the reviewers for their time and thoughtful feedback. We appreciate their recognition of the importance of our study. In response to the concerns raised, we have rewritten and clarified the methods section to better explain our approach and its benchmarking. These revisions aim to address the specific points mentioned and improve the overall clarity of the manuscript. Major concerns: ============ One major issue we found with the paper is lack of basic clarity of the methods description. The paper is simply not self contained. The main text repeats “See Methods”, coupling those with cryptic statements/descriptions. When we go to the Methods section, things do not improve. This leads to a distinct feeling of disregard or neglect which one hopes stayed at the level of writing and did not infiltrate the science itself. Let us give a few examples to make this criticism more concrete: The function of the VIPER algorithm is not clearly described, in the main text, or in the methods, including relevant details such as how it handles missing data. We thought that because VIPER was published some time ago and numerous papers have used it, it was not necessary to include a detailed explanation, we apologize for not doing so. In this revision, we have expanded our benchmarking of splicing factor activity estimation by including 1 three additional approaches: normalized enrichment scores (NES) from gene set enrichment analysis (GSEA), and Spearman and Pearson correlation coefficients. In response to the reviewers’ feedback, we have added a dedicated paragraph in the Results section “Empirical splicing factor networks combined with VIPER outperform computational ones in activity estimation” to conceptually introduce these methods. Additionally, we have updated the Methods section to provide detailed explanations of each approach, including how VIPER handles missing data. On page 12 the authors describe using ARACNe but adjusted with code from the Zhang lab. That code refers to a pre print - this means they are using unpublished algorithms with no assessment of its accuracy to derive their results, which is quite troubling. ARACNe was originally designed to compute pairwise mutual information scores between genes in a gene expression matrix. To adapt if for our study and the mentioned preprint, minimal updates were made to the code to compute mutual information scores between splicing factor gene expression and exon inclusion vectors from two different matrices. These updates are documented in the GitHub repository version history (https://github.com/chaolinzhanglab/ARACNe-AP). Without these adjustments, ARACNe would not function in this context. As part of our collaboration with the Zhang lab, we used their updated code rather than modifying it ourselves. We clarified this in the corresponding Methods section. They then continue to state that they “ran the standard pipeline with a mutual information of …. Using table….” then state “...estimate the mode of regulation using viper::aracne2regulation”. Simply put: Why on earth would you think this makes sense to anyone who does not already know what is going on? Submitting a manuscript with this lack of clarity on what’s being done is concerning. But there are more severe concerns with how the method is executed and the choices of design and parameters along the way. For example, the authors state they used KNNImputer with K=5 to impute missing values since their method can not handle missing values. Why would that make any sense? What does the method do to the data? If some values are close to 0 and some are close to 1 (as sometimes happens in splicing) would you impute 0.5 (which is rare)? We apologize for not clarifying our use of KNNImputer in the Methods section. Missing values in the exon inclusion PSI matrices occur when exons have fewer than 10 reads (see Methods section “Quantification of gene expression and alternative splicing with vast-tools”), and 2 computational network inference methods like ARACNe and MLR cannot handle missing values, necessitating imputation for benchmarking. The KNN algorithm imputes missing values by averaging the values in the K nearest neighbors based on Euclidean distance. This approach preserves local data patterns and assumes that similar samples will have similar exon inclusion values. While the reviewer raises a valid concern about potentially imputing intermediate values (e.g., 50%) in cases where the actual values are close to 0 or 100% in neighboring samples, such a scenario would require highly distinct samples to be highly similar (i.e., low Euclidean distance). To address this, we simulated missing data (20% of exons) in the CardosoMoreira2020 PSI table and tested different K values (k=2, k=5, k=10, k=50, k=100). The imputed values closely resembled real PSI values (Pearson correlation = 1, p < 2.2·10-16), with no overestimation around 50%, supporting our choice of K=5 for imputation (new Supplementary Fig. 6). On page 13 the authors state “threshold of 15% leads to highly accurate and reasonably-sized networks” - based on what? What are the metrics? Where are the evaluations? At the end of that sentence, we had referenced old Supplementary Fig. 1c (new Supplementary Fig. 1c). This figure panel presents the evaluation metrics for empirical networks generated using different absolute delta PSI thresholds (from 5% to 45%). In the revised manuscript, we have clarified this reference to make it more explicit: “Using a 15% delta PSI threshold produced networks that, when combined with VIPER, achieved highly accurate activity estimations (median ROC-AUC = 0.91 across both axes) while limiting the number of target exons (n = 22,049) (Supplementary Fig. 1c)”. Then they state “...the absolute delta PSI as the likelihood of the interaction” - why is that (a basic assumption in all the analysis) reasonable? I would understand if you are unsure if the change is small that it’s real but for example a change of 0.5 and a change of 0.25 in dPSI are both pretty large/confident. The likelihood of interaction helps the VIPER algorithm determine which splicing factors best explain observed exon inclusion signatures. These likelihoods are transformed into ranks internally, and ties can lead to random rankings. For example, assigning a likelihood of 1 to all interactions with absolute delta PSI above 25%, would create ties and result in unstable splicing factor activity estimations. Using the raw absolute delta PSI allows the algorithm to better distinguish between splicing factors with similar interaction networks but different likelihoods. The authors also state “To keep consistency between all types of perturbations, we multiplied the delta PSI in overexpression experiments by -1” - this is another major assumption that is concerning. Clearly if inclusion level is say 0.1 and over expression leads to 0.4 this manipulation will give PSI of -0.3 violating the basic assumptions of PSI in 0-1 range. None of this is discussed, not to mention evaluated in the manuscript. 3 We believe there is a misunderstanding. We do not alter the PSI values, but rather the sign of delta PSIs. As noted “...we multiplied the delta PSI in overexpression experiments by -1…”. This means that if the delta PSI is -14%, we change the sign to 14% for overexpression experiments. This does not affect the original PSI values, which remain within the 0-100% range. We have clarified this in the revised manuscript. Other examples: What are “extreme values” in page 15? What do you mean “within” and “across” (same paragraph p.15)? In this context, “extreme values” refers to estimated splicing factor activities with large values, while “within” and “across” were intended to describe evaluation metrics across samples (“across”) or splicing factors (“within”) of a benchmark dataset. However, recognizing that these terms were unclear, we have revised them to “across perturbations” and “across splicing factors” in the main text, and illustrated in new Fig. 2b. These terms are now introduced in a new Results section titled “Benchmark datasets for splicing factor network inference and activity estimation”. Additionally, following reviewer #2’s advice, we have switched to using ROC AUC metrics for a clearer and more informative assessment. We have also rewritten this section of the Methods to improve clarity. Additionally, differential gene expression analysis should be performed using specifically developed tools such as DESeq2. Following the reviewer’s advice, we re-analyzed differential gene expression using DESeq2 rather than Wilcoxon Rank Sum tests. This approach yielded a highly similar number of differentially expressed splicing factors and did not alter our conclusions. We have updated Supplementary Fig. 4 and the Methods to reflect this change. All of the above makes this reviewer uncertain about the validation of the analysis, making it difficult to understand the rationale for using this algorithm, or the value it adds to the analysis. To this point, it would also be helpful to evaluate an alternative approach or baseline method of some type for comparison, to demonstrate the utility of VIPER. For example, how do the results from VIPER compare to just looking at the overlap of differentially spliced exons from ENCODE or other perturbation datasets, or using CLIP-seq/motif based inference of regulation? Following this and other reviewers’ suggestions, we benchmarked the utility of VIPER against three alternative methods: gene set enrichment analysis (GSEA), Pearson correlation, and 4 Spearman correlation (new Fig. 2 and new Supplementary Figs. 1 and 2). VIPER consistently outperformed the other methods, with Pearson correlation achieving reasonably close performance but showing less robustness to variations in splicing factor networks (new Supplementary Fig. 1c and e). GSEA, serving as an overlap-based benchmark, returned metrics close to random in most cases and only outperformed VIPER and correlation methods in analyses using MLR-inferred splicing networks. However, these show median ROC-AUC scores lower than 0.6 in all cases (new Fig. 2c). Regarding using CLIP inference of regulation, our goal was to develop a general method for estimating splicing factor activity, including both RNA-binding and non-RNA-binding splicing factors. Although CLIP datasets focus only on RNA-binding proteins (RBPs), we benchmarked splicing factor activity analysis using splicing factor networks derived from CLIP experiments. Specifically, we constructed networks considering CLIP peaks near exons (500bp window) in the POSTAR3 database. These CLIP-based networks performed poorly, showing near-random evaluation metrics. However, combining CLIP interactions with the likelihood and mode of regulation from empirical networks significantly improved the accuracy of VIPER splicing factor activity estimates (new Supplementary Fig. 2a-c). The methods also mention that some runs were done using metaVIPER, a newer version of the algorithm, but the function of this algorithm is not described either. There is a similar issue for the ARCANE algorithm, as it is difficult for the reader to understand the function of this algorithm, or the rationale for using it, based on the main text and methods of this manuscript. We apologize for the confusion. MetaVIPER1 is an internal enhancement of VIPER that averages estimated activities when multiple networks are available for a given regulator. Since it is part of VIPER itself, distinguishing between them was unnecessary. We have clarified this in the revised Methods section, along with a clearer explanation of the ARACNe algorithm. One major issue with the paper is how datasets are integrated. We appreciate the authors integration of diverse RBP perturbation datasets, but batch effects are known to be a significant issue for ENCORE RBP perturbation datasets produced at different times/locations. Further, the authors are integrating many highly heterogeneous datasets obtained from ENA, which are likely to have significant batch effects. Additionally, based on the description in the methods, it is difficult for the reader to understand how the ENA datasets were grouped in constructing the empirical network. In general, the issue of batch effects seems to go completely ignored in this work. We agree with the reviewers that integrating datasets from different sources may introduce batch 5 effects. To address this, we have added a new Results section, “Stress testing splicing factor activity estimation using empirical networks”. In this section, among other benchmarks, we evaluate metrics by holding out different benchmark datasets from the empirical splicing factor network. VIPER consistently showed high performance across held-out datasets, with an average median ROC-AUC of ~0.9 across both evaluation axes (new Supplementary Fig. 1d). Regarding the description of how ENA datasets were grouped, perturbations targeting the same splicing factor coming from independent studies were kept separate. Following also Reviewer #2’s feedback, we have rewritten the Methods section “Inference of splicing factor networks” to clarify this. Minor comments: ============== Given the authors’ finding of correlation between SF activity and proliferation in CCLE cell lines, it may be worthwhile to analyze DepMap CRISPR screening data to check how these cell lines depend on the relevant SFs for proliferation. We analyzed the distributions of DepMap DEMETER2 gene dependencies considering our classification of splicing factors into cancer splicing programs. We also considered whether the tested cell lines were annotated as cancerous or non-cancerous (fibroblasts (HS729_SOFT_TISSUE) or engineered (HEKTE_KIDNEY, OELE_OVARY, PRECLH_PROSTATE)) and whether the splicing factor was known as a pan-essential gene2. As expected, the distributions for essential splicing factors showed negative gene dependencies across all cases, as shown in new Supplementary Fig. 5d. For non-essential splicing factors, gene dependencies were closer to 0, but no significant differences were observed between tumor suppressor and oncogenic splicing factors in cancerous and engineered cells. Interestingly, in fibroblasts, we observed a significant tendency (p = 0.028) for tumor suppressor splicing factors to exhibit higher gene dependencies than oncogenic ones, aligning with their expected behavior. Due to the context-dependent effects of many SFs, the authors should address how well their model generalizes between different cell types and conditions. This could provide better justification for the reliability of the model’s predictions on diverse cancer datasets. As suggested, we evaluated how well empirical splicing factor networks from different benchmark datasets generalize to held-out ones. This analysis, shown in new Supplementary Fig. 2d, demonstrates consistently high VIPER performance, even when using ENCODE-derived 6 networks (from only two cell types, K562 and HepG2) to evaluate the ENA dataset. Additionally, we assessed whether accuracy differs between different types of splicing factors, namely, core, RBPs, and others (new Fig. 2e). Core splicing factors produced higher evaluation metrics across benchmark datasets, likely because their exon inclusion signatures are more consistent across conditions. However, even for RBPs and other context-specific splicing factors, VIPER’s predictions were still better than random (new Fig. 2f). It would be helpful to show which SF(s) may be responsible for regulating the SEC22B exon which may affect ICB response, and the specific pathways enriched in the targets of each SF. This would be very helpful to demonstrate the method’s ability to understand the specific mechanisms of regulation of cancer-related exons. It would also be helpful to know whether these SFs are themselves implicated in regulating ICB response by the in vivo CRISPR screen. Following Reviewer #1’s comment, we inspected which splicing factors regulate SEC22B exon in our empirical splicing factor networks and performed an over-representation analysis of all targets of these splicing factors. This exon is regulated by three different tumor suppressor-like splicing factors (PRPF8, SRSF1, SF3B3) whose targets are enriched in ReactomeDB terms mostly 7 related to cell cycle control (FDR < 0.05) (new Supplementary Fig. 5a). However, we also observe enrichment in Class I MHC-mediated antigen presentation. This suggests alterations in the activity of these splicing factors could influence tumor ICB response. To further investigate this, we assessed whether these splicing factors themselves impact ICB response using in vivo CRISPR screen data. Among the three, silencing SF3B3 produces the most similar impairment (score = -0.89) in ICB response to SEC22B (score = -1.2), underscoring its likely regulatory role on HsaEX1036341 exon in SEC22B. In general, it would be helpful to clarify whether oncogenic/tumor suppressor-like SFs seem to have similar effects or if they have distinct phenotypes. All oncogenic or tumor suppressor-like SFs are lumped together for pathway analysis, when looking at correlation with CCLE proliferation, and when looking at carcinogenesis. How do different oncogenic/tumor suppressor- like SFs affect different cancer phenotypes, or different cancer types? In old Supplementary Fig. 3a (new Supplementary Fig. 4a), we show the differential activation or inactivation of each splicing factor across all TCGA cancer cohorts. To maintain coherence, splicing factors are sorted in the same order as in old Fig. 3c (new Fig. 3c). While we examined cancer types with similar splicing factor activity patterns, we did not observe a clear structure or clustering of cancer phenotypes based on these differential activities (new Supplementary Fig. 4a). At a high level, we were unsure given all of the above regarding the overall importance/impact. The computational analysis does not seem to offer a new method but rather adaptation of existing ones for the specific task at hand with a variety of parameter/design choices along the way which are unclear. It does not appear generalizable/usable for other researchers but rather a way to perform specific analysis. As such it hinges on the importance and validity of the derived conclusions/results, some of which look promising. We appreciate the reviewer’s feedback. To ensure transparency and reproducibility, we have made the full code and empirical splicing factor networks available in a public GitHub repository, allowing other researchers to apply or adapt our method to similar analyses. It is true that our approach builds on existing algorithms. But, while VIPER has shown to be quite accurate in predicting transcription factor activities, it has never been shown it could be used to predict splicing actor activities. By showing we can repurpose VIPER integrating it with splicing factor network inference, we have created a new method not existing until now for splicing factor activity estimation, particularly in the context of complex regulatory networks. The flexibility of our approach makes it adaptable to different datasets and conditions, providing a generalizable framework for future studies in the field. We have adjusted the tone of the manuscript to reflect its main contribution as a methods paper, emphasizing the broader utility of our approach beyond the specific analysis presented focused in cancer biology. 7 In summary, we appreciate this new approach to an important and challenging topic, and the significant work and integration of diverse datasets that went into this paper. However, there are serious issues with the methods and the description of the methods. We thank the reviewers for their feedback and acknowledge the challenges in understanding our approach. In the revised manuscript, we have provided clearer explanations of the methodology, added further validation, and clarified assumptions. We believe these revisions address the concerns and strengthen the manuscript, offering a more robust framework for interpreting the results. Reviewer #2 (Remarks to the Author): Summary In the manuscript “Disentangling the splicing factor programs underlying complex molecular phenotypes,” Anglada-Girotto et al. create a new database of splice-factor regulated exons and combine it with the already existing VIPER tool to estimate splice-factor activity from RNA-seq data. After some validation, the authors apply the data to estimate splice factor activity in several cancer datasets accompanied by selected follow-up analysis. While the database+VIPER approach is ingenious, and the datasets analyzed many, the manuscript in its current form, however, falls short. A substantial concern is that the validation needs to be extended, especially regarding how well it works in cell lines/tissues that are (very) different from the cell line in the database they create - meaning how well the analysis approach generalizes and thereby how trustworthy the cancer results are. In our opinion, the paper would also benefit from being rethought/re-structured. Is it a paper about a new and useful bioinformatic tool or a paper about understanding cancer better? Currently, the paper falls between these two chairs as it does neither sufficiently. To be successful, the tool angle would need more validation – and a more extensive database as only profiling ~1/2 of known splice factors does not feel comprehensive. As a cancer paper, more novel biological insights are needed since it is already very well-established that splicing, and its regulation, are dysregulated in cancer. We thank reviewers #2, Dimitrios S. Kanakoglou and Kristoffer Vitting-Seerup, PhD, for their careful evaluation of our manuscript. We believe that implementing their constructive feedback has significantly improved the quality and clarity of our work. Please find our detailed point-by- point responses below. Major comments • Validation of the splice factor network is too weak: 8 1. Regarding the splicing factor-> exon networks, please also provide other metrics than recall, such as false discovery rate. Similarly performance on the “completely independent perturbational signatures” should also be shown. We acknowledge the term “recall” was misused in the previous version of the manuscript. In response to this and the following comment, we have recalculated evaluation metrics using areas under receiver operator curves (ROC-AUC) (new Fig. 2 and new Supplementary Fig. 1 and 2). We did not use precision-recall curves, as they could not be stably calculated at the single experiment level, where only one splicing factor was perturbed, leaving only one positive label per experiment. Additionally, we have clarified in the Results section (“Benchmark datasets for splicing factor network inference and activity estimation”) and Methods section (“Evaluation metrics of splicing factor activity estimation”) that all performance metrics were calculated using held-out benchmark datasets, representing completely independent perturbational signatures. 2. The splice factors profiled more than once (e.g. in two different cell lines) should have multiple observations. This also means they can be used to calculate ROC and precision-recall curves (that can be summarized as area under curve) As explained above, the reported evaluation metrics (ROC-AUC) are always calculated using held-out datasets that are independent of those used to generate the splicing factor networks. We have rewritten the corresponding Results and Methods sections to clarify this and added new Fig. 2b, to illustrate how ROC-AUCs are computed across either the perturbations or the splicing factor axis in the benchmark datasets. 3. To show generalization, performance also needs to be assessed when the data evaluated originates from cell types/tissues that are (very) different from the ones used to create the SF- >exon database. Since the database is (mainly) based on cell-line data, evaluation should also be on complex tissues. Something like evaluating the ENA-based dataset when using the ENCODE data as a database? Following this, and other reviewers’ suggestions, we added a new section in Results titled “Stress testing splicing factor activity estimation using empirical networks” that benchmarks how well activity estimation methods perform in the presence of technical and biological artifacts. Specifically, we observed high evaluation metrics for VIPER (median ROC-AUC values of 0.90 and 0.91 across splicing factors and perturbations, respectively) when using ENCODE data to predict ENA perturbation signatures (new Supplementary Fig. 1d). This suggests that combining empirical splicing factor networks with VIPER generalizes well across different data sources and cell types. • More work should be put into doing the computational deconvolution of exon-splice factors as that would increase the usefulness (of the analysis approach tremendously. Profiling ~1/2 of known splice factors does not feel particularly satisfactory or comprehensive for a new tool. What 9 would happen if ARACNe (or similar) was run on e.g. all the TCGA data the authors already quantified? We agree that having networks for all splicing factors would be ideal for a comprehensive analysis. Following the reviewers’ suggestion, we ran ARACNe and MLR to generate computational splicing factor networks using all solid tissue normal (n = 736) and primary tumor (n = 9,854) samples from TCGA. However, these computational networks did not perform as well as empirical networks (new Fig. 2c). In the Results section “Poor mode of regulation prediction limits performance of computational networks”, we show the limited performance stems from the incorrect inference of the splicing factor→exon mode of regulation (new Supplementary Fig. 2d- f). We have highlighted the need for novel algorithms to address this issue in the Discussion section as well. Since reverse-engineering splicing factor networks is a complex challenge and the focus of our article is on splicing factor activity estimation, we proceeded with empirical networks due to their consistent accuracy. While we currently cover ~50% of annotated splicing factors, this represents the largest compendium of splicing factor-caused exon inclusion signatures from perturbation experiments. Additionally, as noted in the Discussion, our approach is easily extendable as new datasets become available. • If this is to be a cancer paper, more novel biological insights are needed since it is already very well-established that splicing, and its regulation, are dysregulated in cancer. Here a low-hanging fruit could be to create an additional database this time with exons affected by transcription factors knock-out and also apply that to the cancer datasets. Our approach builds on previous work that uses gene expression to infer transcription factor activity, with a focus on splicing factors controlling exon inclusion. Transcription factor activity in cancer has already been extensively studied3–6. Our goal is to present a method for accurate splicing factor activity estimation, drawing on established approaches for transcription factor analysis. Cancer is used here as an illustrative example to show how splicing factor activity analysis offers a distinct and relevant perspective on splicing regulation in complex molecular phenotypes. While interesting and potentially addressable in the future, investigating how transcription factors indirectly affect splicing is beyond the scope of our current study. • The title does not adequately reflect the study done and must be redone. Both title and abstract is to vague. After reading title and abstract neither of us understood what the study was about and what problem was solved. We appreciate the feedback regarding the title and abstract. We agree that both need to better reflect the scope and contributions of the study. We have revised the title and abstract to be more specific and clearer about the method developed and the problem addressed. • While the method section is very detailed it needs to be more pedagogical – it is very hard to recreate what was done. 10 We appreciate the reviewer’s suggestion and understand the need for clarity in the Methods section. In the revised manuscript, we have made the Methods section more pedagogical, providing step-by-step explanations to improve reproducibility and ensure that readers can more easily follow and recreate our approach. • The Github/Code needs extensive updates about installation instructions, dependencies, and environmental requirements. Currently, it is hours of work of trial and error to get it running. We have updated the instructions in both GitHub repositories, providing detailed steps on installation, dependencies, and environment requirements. These updates cover running VIPER with empirical splicing factor networks (https://github.com/MiqG/viper_splicing) and reproducing our analyses (https://github.com/MiqG/publication_viper_splicing). Additionally, we encourage reviewers to report any specific issues or questions they encounter in the future as a \"GitHub issue\" so we can promptly address any unforeseen problems. Minor Comments • Abstract is sub-par it is a bit to vague in both what the problem is, how this paper helps and what the impact of that is. The impact/perspective is also missing from the discussion We appreciate the reviewer’s feedback on the abstract and discussion. In the revised manuscript, we have made the abstract more specific, clearly stating the problem addressed, how our method contributes to solving it, and the impact of this work. We have also updated the discussion to better highlight the broader implications and potential impact of our approach, providing a clearer perspective on how this method advances the field of splicing factor activity estimation and its relevance to complex molecular phenotypes like cancer. • There is a tendency towards using very long sentences reducing redability In the revised manuscript, we have shortened and simplified sentences throughout to improve readability and clarity. We hope these changes make the text easier to follow while maintaining the precision needed to convey the findings. • Introduction fails to mention how alternative start sites and alternative termination sites are part of isoform generation In the revised manuscript, we have expanded the introduction to include the role of alternative start sites and alternative termination sites in isoform generation. • Discussion of how splicing-factor→pre-mRNA exon interaction networks are likely cell/tissue/disease specific is needed 11 As discussed in previous responses, our evaluation of performance metrics based on splicing factor type revealed that empirical networks for core splicing factors generalize better than those for RBPs or other splicing factors, likely due to the context-specific regulation of RBPs (new Fig. 2e). • Why are there 4 ENA based networks? And how were these 4 defined/sperated In the ENA dataset, we collected all publicly available splicing factor RNA-seq perturbation signatures. In some cases, we found multiple experiments perturbing the same splicing factor. Hence, each independent signature generates a unique splicing factor→exon network. To maintain consistency with the ENCODE dataset –where networks were separated by cell line (HepG2 or K562) and perturbation type (CRISPR knockout or shRNA knockdown)– we separated the ENA signatures into four groups, ensuring each resulting network reflects a distinct experiment. We apologize for the lack of clarity in the Methods section and have revised it accordingly. • Why is the ARACNe algorithm run on such a small dataset? Would increasing the number and diversity of samples not increase power? Something like running it on GTEx. We addressed this comment in the Major comments section. Briefly, we performed computational splicing factor network inference using TCGA data, but increasing the number and diversity of samples did not substantially improve performance (new Fig. 2c). We demonstrate the main issue lies in accurately inferring the sign of the mode of regulation between splicing factors and target exons, rather than statistical power. • How does the approach work with existing databases like Polvèche et al 2023 (doi.org/10.1093/database/baad091) Povlèche et al.7 generated a database, SplicingLore, of alternatively spliced exons for 75 different splicing factors by downloading and processing publicly available data from perturbation experiments. Following Reviewer #2’s suggestion, we benchmarked the predictive power of using these networks from SplicingLore as input for splicing factor activity estimation. We mapped SplicingLore exons to VastDB annotations to ensure compatibility and computed evaluation metrics against our benchmark datasets. We considered whether perturbation experiments overlapped with benchmark datasets to avoid overfitting. SplicingLore’s networks generated evaluation metrics comparable to those produced by our empirical networks for the same splicing factors, underscoring the robustness of our approach (Response Fig. 1). This also highlights the flexibility of our method, demonstrating that our framework can be adapted to other empirical datasets generated through different quantification pipelines. 12 Response Fig. 1. Benchmark evaluation metrics for empirical splicing factor networks generated from data processed in this study and in SplicingLore7. • The introduction should mention how the importance of splicing in cancer is already well established We have revised the Introduction to emphasize that the importance of splicing in cancer is already well established. • Survival analysis: 1. How were the factors in Figure 3E chosen? They are not among the genes highlighted in 3C. How many of the “cancer-cohorts” are significant after correcting for multiple testing? I also assuming you tested all factors instead of only the two factors shown here so how do you correct for multiple testing there? We apologize for mistakenly using the word “significantly” in the main text describing this analysis which only focused on the sign of the coefficient of association of Cox Proportional Hazards models (Cox PH). We used a Cox PH model to regress the estimated activity of each splicing factor against survival probability across TCGA primary tumor cohorts. We then counted the number of cohorts where the association coefficient was positive (“High risk”) or negative (“Low risk”). Although we performed a statistical test, we did not focus on p-values, but rather on whether the association sign was consistently positive or negative across cohorts, in line with the classification of splicing factors as oncogenic-like or tumor suppressor-like. This is then statistically tested in old Fig 3f (new Fig. 4f). The splicing factors in old Fig. 3e (HNRNPD and SNRNP200) were selected from those in the cancer splicing programs (old Fig. 3c) with the most frequent positive (HNRNPD) and negative (SNRNP200) associations with survival. They were not highlighted in old Fig. 3c as they were not 13 among the top 5 recurrently differentially activated splicing factors across cohorts. It is important to note that differential splicing factor activity (old Fig. 3c) compares primary tumor vs. normal tissue, while Cox PH analysis (old Fig. 3e) correlates splicing factor activity with survival probability in primary tumors, which explains the difference in outcomes. In the updated analysis, DNMT1 replaces SNRNP200 as the splicing factor most recurrently associated with “Low risk” in new Fig. 4e. This change, due to a code error identified while rerunning our pipeline, does not alter the study’s conclusions. We have clarified this in the revised Results section and figure footer. 2. To carry any validity surivival analys needs to be corrected for cofactors such as gender age and tumor stage. We agree that sex, age, and tumor stage may confound associations. Following the reviewers’ suggestion, we re-ran the Cox PH analysis, incorporating these factors as cofactors. However, for 12 cancer cohorts (CESC, DLBC, GBM, LAML, LGG, OV, PCPG, PRAD, SARC, THYM, UCEC, UCS), tumor stage information was unavailable, preventing their inclusion. As shown in Response Fig. 2, the sign distributions of the associations remained consistent, with significant differences observed in the recurrent association patterns of programs defined through splicing factor activities, but not for those defined through their expression, confirming the robustness of our findings. However, we acknowledge that these cofactors may not solely act as confounders but could also directly influence splicing factor activity. Adjusting for them might unintentionally remove relevant biological signals. Thus, we have opted not to include this analysis in the manuscript. Response Fig. 2. Impact of cofactors on the association between splicing factor activity/gene expression and survival probability. Distributions showing the number of cancer cohorts with consistent association signs (“High Risk”>0, and “Low Risk”<0) using Cox PH models for each splicing factor, comparing splicing factor activity (left) and gene expression (right) with survival probability across primary tumor samples. Associations were adjusted for age, sex, and tumor stage. Splicing factors in each program are selected based on recurrent differential activity 14 (left) or up/downregulation (right) in at least 5 cancer cohorts, comparing primary tumors with solid tissue normal (FDR < 0.05). • The ORA analysis in the “Cancer splicing programs regulate cancer hallmarks” section is highly over-interpreted. A statement like “exons regulated by the oncogenic-like splicing program are not enriched in immune evasion-related terms” does not hold much validity when only the top 10 gene-sets are considered since gene-set ranking is heavily biased by e.g. size of gene-set. In other words gene-set ranking is not the same as biological importance (or completeness) and hence cannot be interpreted this way. Furthermore the difference between the enrichment scores of the “ongogenic” and “tumorsupressor” in figure 4B is very small for most cases where they are both represented? For the ORA what gene-background is used and why? Carefull argumentation is needed as that is essential to the interpretation of the ORA We believe there was a misunderstanding regarding how we subsetted the ORA results. To explore the results, we ranked enriched terms (FDR<0.05) by gene ratio to highlight those terms with the largest overlaps for each program, as these summarize the genes targeted by the splicing factors in each cancer splicing program. We then used these 20 terms to subset the ORA results for both programs, ensuring comparability. For example, the term “RAC1 GTPase Cycle” appears in both splicing programs but it is highlighted in old Fig. 4b (new Fig. 5b) because it ranks among the top 10 enriched terms for the oncogenic-like splicing factor program. The statement “exons regulated by the oncogenic-like splicing program are not enriched in immune evasion-related terms” refers to the fact that those immune-related terms are not significantly enriched among the targets of the oncogenic-like program, in contrast with the tumor suppressor program. We have clarified this in the revised Results section “Linking cancer splicing program dysregulation to cancer hallmarks” and included these results in new Supplementary Table 4. Regarding the small differences in gene ratio, we do not claim any specific interpretation of magnitude, but rather point out that some terms with overlapping genes are significantly enriched in either or both programs. We used the default gene background of all genes in ReactomeDB and have now clarified this in the Methods section. • What gene-annotation (or rather exon annotation) database were quantified? Gencode? Ensemble? Refseq? And how are overlapping but non-identical exons handled (caused by e,g. an alternative donor site?) As described in the Methods section, we used vast-tools for exon inclusion quantification, which relies on the VastDB splicing event annotation database. Detailed information about vast-tools and its underlying database can be found in the tool’s documentation (https://github.com/vastgroup/vast-tools) and associated publications8–10, as well as relevant benchmarks11–13. 15 Regarding VastDB’s exon annotation, we consulted Dr. Manuel Irimia, the lead developer, who confirmed that VastDB is based on ENSEMBL and incorporates large amounts of publicly available RNA-seq data to detect previously unannotated exon-exon junctions and isoforms (details in 8). For exons with alternative donor/acceptor sites, vast-tools generates a single PSI value per exon ID by aggregating all exon-exon junctions associated with the exon variants, similar to gene-level quantifications pooling isoform level data. • Create a detailed guide that walks users through the analysis process. Explain how to execute each workflow, including any prerequisites or configuration steps. Provide examples of command-line instructions for running the Snakemake workflows, as well as expected outputs. We thank the reviewers for their thorough consideration of our study. We refactored our analysis pipeline repository to better guide readers through the analysis and facilitate reproducibility with the command line instructions to run the worfkflows and the expected outputs in each case. • Watch out for causation vs correlation in the text We are unsure of the specific instances the reviewers are referring to in this comment. However, in the revised manuscript, we have carefully reviewed the language and ensured that we only use causation where appropriate. We look forward to read and updated manuscript. Dimitrios S. Kanakoglou PhD Fellow Section for Bioinformatics, Health Technology Technical University of Denmark (DTU) Denmark And Kristoffer Vitting-Seerup, PhD Associate Professor Section for Bioinformatics, Health Technology Technical University of Denmark (DTU) Denmark Reviewer #3 (Remarks to the Author): In their manuscript ‘Disentangling the splicing factor programs underlying complex molecular phenotypes’ Serrano and colleagues aim to use a bioinformatics approach to characterize alternative splicing programs that are controlled by individual splicing factors. The authors use an 16 extensive set (hundreds) of published RNA-Seq datasets with altered splicing factor expression, e.g. knockdown, knockout or overexpression of specific splicing regulatory proteins, to characterize alternative splicing signatures of splicing factors. In contrast to some previous work they do not rely on alteration of the mRNA of a given splicing factor, which can be misleading due to posttranscriptional/posttranslational regulation of their activity, but derive splicing factor activity from previously published datasets with known splicing factor alterations. They use a training dataset to define splicing patterns caused by alteration of a given splicing factor and then use this to infer misregulation of splicing factors based on global splicing patterns in disease conditions, e.g. different types of cancer. This is an important question and an interesting approach that has been taken for transcription factors before. However, there is a fundamental problem with applying this concept to the regulation caused by a specific splicing factor: There is an extensive amount of cross regulation between splicing factors as has been documented over many years by many labs. That means that altering the activity of one splicing factor (especially by artificial knockdown or overexpression approach as used for the training dataset) will lead to the alteration of the expression and/or activity of (many) other splicing factors. Given this extensive network of splicing factor cross regulation, the approach taken here will not yield meaningful results for individual splicing factors as it is never clear whether an effect is direct or indirect. For such an approach the training dataset would have to rely on other methods, e.g. acute depletion of splicing factors using inducible degrons or knockdown-rescue approaches and the combination with CLIP data to identify high confidence direct targets. With the current datasets and approach the results are too diffuse and potentially nonspecific. As such, the findings do not provide a sufficient advance for publication at this level. We thank the reviewer for taking the time to provide thoughtful feedback and for recognizing the relevance of our study. We fully acknowledge the extensive cross-regulation between splicing factors, as mentioned by the reviewer and in our manuscript. But we think that the concern about the meaningfulness of our results due to splicing factor cross-regulation, may be due to a misunderstanding. Protein activity estimation approaches like VIPER have proven very accurate at estimating transcription factor activity despite their extensive cross- and post-translational regulation14. Hence, the primary aim of our study was to port the method from transcription factors to splicing factors in order to have an accurate splicing factor activity estimation method than the commonly used approach of relying on splicing factor gene expression as a proxy for its activity. In our empirical splicing factor networks, we include all exons whose inclusion changes upon the perturbation of a specific splicing factor, regardless of whether the target is direct or indirect. This was intentional, as our method is designed to capture the overall splicing program driven by a splicing factor, not limited to its direct targets. Interestingly, we do not observe high cross- regulation between splicing factors in these perturbation-derived networks, as each exon is 17 regulated by a median of only two splicing factors (new Fig. 2d). Despite the indirect nature of these networks, we consistently observe that perturbing specific splicing factors generates characteristic exon inclusion signatures, which enables accurate splicing factor activity estimation across various conditions. Hence, Reviewer #3’s statement that extensive cross-regulation would prevent identifying the splicing factor responsible for a given splicing signature does not hold in our study. Our results demonstrate high accuracy in estimating splicing factor activity, validated not only by single-perturbation benchmarks but also by more complex scenarios. For example, our method successfully recapitulates the mechanism of action of Indisulam, a drug that inactivates RBM39 (new Fig. 3a-c). This would be missed if RBM39 mRNA levels were used as a proxy of its activity, as Indisulam treatment increases RBM39 mRNA levels. Moreover, our approach correctly identifies perturbed splicing factors in combined knockdown experiments (new Fig. 3d), and it accurately predicts the mechanism of action of five drugs targeting core splicing factors by assigning splicing factor activities that correlate with protein-protein interaction networks (new Fig. 3e-g). These validations were conducted using independent studies and datasets, showcasing the reproducibility and robustness of our approach. In response to Reviewer #3’s suggestion, we benchmarked splicing factor networks derived from POSTAR3, a database of CLIP experiments. As CLIP data only covers RNA-binding protein (RBP) splicing factors, the benchmark was restricted to these (new Supplementary Fig. 2a-c). As expected, CLIP data alone shows near random predictive power. But, consistent with the reviewer’s comment, incorporating splicing factor→exon interactions identified by both CLIP and perturbation data improved the accuracy of our method (new Supplementary Fig. 2c). This suggests that, when possible, identifying direct splicing factor→exon interactions enhances the precision of splicing factor activity estimation. However, direct targets can only be identified for a subset of splicing factors due to data limitations. In summary, combining empirical splicing factor networks with VIPER provides reliable splicing factor activity estimation across diverse conditions, even in the presence of extensive splicing factor cross-regulation –as observed for transcription factors. Despite not always being able to identify direct targets, our approach captures the functional output of splicing factors, offering accurate and reproducible estimates of their activity. In addition, the enrichment of oncogenic or tumor-suppressive splicing factors appears mild, and it is unclear whether the observed enrichment in cell cycle related-processes are a cause or consequence of tumorigenesis. In our analysis of cancer samples, the goal was to demonstrate how our splicing factor activity estimation approach can provide new insights into such a complex molecular phenotype as cancer. We defined two cancer splicing programs based on recurrent activation or inactivation using a threshold of 5 cohorts. We are unsure what the reviewer considers “mild enrichment”. As shown 18 in new Fig. 4 and new Supplementary Fig. 4, our classification of splicing factors by splicing program is more functionally relevant (e.g., as shown by recurrent association with patient outcome in new Fig. 4) than an equivalent classification based on differential overexpression or downregulation (new Supplementary Fig. 4). The activity of oncogenic- and tumor suppressor-like splicing factors correlates with high and low risk of death across TCGA cohorts, respectively, matching their expected functional behavior (new Fig. 4f). In contrast, gene expression-based cancer splicing programs defined with the same criteria show no association with patient survival (new Supplementary Fig. 4e). This pattern holds even when increasing the threshold for defining cancer splicing programs (new Supplementary Fig. 4b and f). Additionally, the activity of these cancer splicing programs correlates with cancer cell proliferation, and they coordinate their activation during tumorigenesis in line with their classification (new Fig. 5). For example, oncogenic-like splicing factors correlate positively with MKI67 expression, and tumor suppressor-like splicing factors correlate negatively (new Fig. 5e). Furthermore, oncogenic-like splicing factors are activated, while tumor suppressor-like factors are inactivated, as fibroblasts transition to a tumorigenic state (new Fig. 5f-g). These findings suggest that the cancer splicing programs identified provide a distinct and functionally relevant perspective on splicing factor activity in cancer. While our data cannot definitely determine whether these splicing programs are a cause or consequence of tumorigenesis, addressing this would require further experimentation beyond the scope of our study. Nonetheless, our analyses of splicing factor dependencies (suggested by Reviewers #2) and carcinogenesis samples do not suggest splicing factors as the main drivers of the process. Rather, it indicates that oncogenes rewire the activation of splicing factors, which then implement cancer splicing programs. Our cancer analysis showcases how splicing factor activity estimation can generate novel, testable hypotheses using existing data. While our data cannot definitely determine whether these splicing programs are a cause or consequence of tumorigenesis, addressing this would require further experimentation beyond the scope of our study. Nonetheless, our analyses of splicing factor dependencies (suggested by Reviewers #2) and carcinogenesis samples do not suggest splicing factors as the main drivers of the process. Rather, it indicates that oncogenes rewire the activation of splicing factors, which then implement cancer splicing programs. Our cancer analysis showcases how splicing factor activity estimation can generate novel, testable hypotheses using existing data. The mechanistic cause behind splicing factor alterations remains unclear, as does the connection between these factors—whether they are enriched at specific stages of the spliceosome or involve certain RNA-binding protein families. Uncovering the mechanistic causes behind splicing factor alterations would require additional experimentation, which is beyond the scope of this study. Our primary aim here is to develop an approach for accurate splicing factor activity estimation, with cancer biopsy analysis serving as a case study to illustrate its utility. Following Reviewer #3’s suggestion, we investigated whether splicing factors involved in cancer splicing programs are associated with specific spliceosome stages or RNA-binding protein families. Using SpliceosomeDB15 annotations, we found no enrichment for any particular spliceosome stage, with only 32 of 122 splicing factors in cancer splicing programs associated with core splicing reaction stages (new Fig. 4g). Additionally, we examined RBP family memberships using RBPDB16 annotations and found that 54 splicing factors were part of known RBP families (new Fig. 4h). 19 It is also uncertain if these expression changes are coordinated. We agree with the reviewer that the regulation of splicing factor activity coordination remains unclear, and we believe this is an important avenue for future work. Our findings show that changes in splicing factor activity are indeed coordinated, as evidenced by consistent patterns of recurrent activation/inactivation, associations with patient outcomes, and coordinated shifts during tumorigenesis. However, we found no consistent association between mRNA expression changes and differential splicing factor activity (new Fig. 4 and new Supplementary Fig. 4), indicating that expression alone as happens with transcription factors3 does not account for these activity shifts. To further explore this, we plotted splicing factor expression changes across carcinogenesis (new Supplementary Fig. 5c). Despite widespread mRNA upregulation of splicing factors upon tumorigenesis, these changes do not explain splicing program-specific activity alterations. This reinforces the idea that splicing factor expression alone is insufficient to capture regulatory shifts, underscoring the need for methods, such as ours, that can capture splicing factor activity directly from exon inclusion data. This analysis has been incorporated into the updated Results section, “Cancer splicing programs exhibit coordinated regulation during carcinogenesis”. The carcinogenesis cell culture system could be utilized to investigate the functional significance of splicing programs, for example, through perturbation experiments targeting key oncogenic splicing factors. As noted, the primary aim of our study is to develop accurate splicing factor activity estimation. While conducting functional perturbation experiments targeting key oncogenic splicing factors would be valuable, these experiments fall outside the scope of this study. Although we cannot directly assess the functional significance of the identified cancer splicing programs for carcinogenesis, our method provides a robust framework for identifying and detecting these programs, which can be used to explore such questions in future research. In summary, the findings related to cancer are interesting but lack specificity and as presented do not present a significant advance over current knowledge. As mentioned above, our identification and characterization of cancer splicing programs from the perspective of splicing factor activity is intended to showcase our approach can generate novel, testable hypotheses from existing datasets. We have revised the manuscript to better emphasize this point. Reviewer #4 (Remarks to the Author): 20 I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts. We thank Reviewer #4, for the time and effort dedicated to revising our manuscript. Reviewer #5 (Remarks to the Author): I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts. We thank Reviewer #5, for the time and effort dedicated to revising our manuscript. References 1. Ding, H. et al. Quantitative assessment of protein activity in orphan tissues and single cells using the metaVIPER algorithm. Nature Communications 9, 1–10 (2018). 2. Anglada-Girotto, M. et al. In silico RNA isoform screening to identify potential cancer driver exons with therapeutic applications. Nature Communications 15, 1–19 (2024). 3. Alvarez, M. J. et al. Functional characterization of somatic mutations in cancer using network-based inference of protein activity. Nat. Genet. 48, 838–847 (2016). 4. Jethalia, M., Jani, S. P., Ceccarelli, M. & Mall, R. Pancancer network analysis reveals key master regulators for cancer invasiveness. J. Transl. Med. 21, 1–18 (2023). 5. Paull, E. O. et al. A modular master regulator landscape controls cancer transcriptional identity. Cell 184, 334–351.e20 (2021). 6. Mall, R. et al. Network-based identification of key master regulators associated with an immune-silent cancer phenotype. Brief. Bioinform. 22, bbab168 (2021). 7. Polvèche, H. et al. SplicingLore: a web resource for studying the regulation of cassette exons by human splicing factors. Database 2023, baad091 (2023). 8. Tapial, J. 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. 27, 1759–1768 (2017). 9. Irimia, M. et al. A Highly Conserved Program of Neuronal Microexons Is Misregulated in Autistic Brains. Cell 159, 1511–1523 (2014). 10. Braunschweig, U. et al. Widespread intron retention in mammals functionally tunes transcriptomes. Genome Res. 24, 1774–1786 (2014). 11. Sterne-Weiler, T., Weatheritt, R. J., Best, A. J., Ha, K. C. H. & Blencowe, B. J. Efficient and Accurate Quantitative Profiling of Alternative Splicing Patterns of Any Complexity on a Laptop. Mol. Cell 72, (2018). 12. Jiang, M., Zhang, S., Yin, H., Zhuo, Z. & Meng, G. A comprehensive benchmarking of differential splicing tools for RNA-seq analysis at the event level. Brief. Bioinform. 24, (2023). 13. Parada, G. E. et al. MicroExonator enables systematic discovery and quantification of microexons across mouse embryonic development. Genome Biol. 22, 1–26 (2021). 14. Garcia-Alonso, L., Holland, C. H., Ibrahim, M. M., Turei, D. & Saez-Rodriguez, J. Benchmark and integration of resources for the estimation of human transcription factor activities. Genome research 29, (2019). 21 15. Cvitkovic, I. & Jurica, M. S. Spliceosome Database: a tool for tracking components of the spliceosome. Nucleic Acids Res. 41, D132 (2013). 16. Cook, K. B., Kazan, H., Zuberi, K., Morris, Q. & Hughes, T. R. RBPDB: a database of RNA- binding specificities. Nucleic Acids Res 39, D301–D308 (2010). 22"}, {"round": 2, "reviewer_comments": {"reviewer_2": "The article is now much more readable and understandable. We especially appreciate the improvements in the methods section and their GitHub. It is clear the authors put a lot of work into this revision, and the effort clearly worth it. Major concern (singular since all comments below are about the validation): • Since there is only one observation, AUC/ROC is not adequate. AUROC can be very misleading when data is imbalanced (like it is here). Please use better metrics – with emphasis on the plural. • Seem that a ranking would do nicely as a metric here? If the right one is ranked 4 of 240 SF the ranking would be top 5/240 = top 2.08%. Such a metric would be much more informatice and easy to interpret. But feel free to explore others in suplementary figures. • The ENA data contains K562 and HepG2 data and is therefore not an independent dataset. Please divide the ENA data so that cell lines do not appear in both train/test data. • Specifically “train” on K562 and HepG2 data and apply to other cell lines • For the analysis where ENA (without K562 and HepG2 ) are held out you need to show the data in more detail. Specifically you need to show which cell lines and SF have poor AUCs and also annotate annotate tissue orign of cell types. This is the only way we can is if the method generalizes to other tissues/cell types. • Please speculate why there are <0.5 AUCs both across SF and across datatypes • Should some SF be excluded from all downstream analysis since they are not trustworthy? • Remove the held-out ENCODE data (at least from the main figures) – non of them are indepdent and hence they are misleading. Minor comments (all about survival) • Survival analysis: Change from “high” and “low” to “increased” and “decreased”. A non-significant HR > 0 cannot be called \"High Risk\". • Please also show a similar survival based plot but this time counting significant changes. The defintion of a significant association is that the HR is reliable - now you are overinterpreting point values with very high uncertainty (meaning the sign of the HR is not trustworthy). • We don’t quite understand why the authors do not replace the survival analysis with the new one (the one with co-factors) since we agree the conclusions are unchanged – but now much more reliable? We look forward to read and updated manuscript. Dimitrios S. Kanakoglou PhD Fellow Section for Bioinformatics, Health Technology Technical University of Denmark (DTU) Denmark And Kristoffer Vitting-Seerup, PhD Associate Professor Section for Bioinformatics, Health Technology Technical University of Denmark (DTU) Denmark", "reviewer_3": "We are reviewing the revised version of the manuscript ‘Exon inclusion signatures enable accurate estimation of splicing factor activity’ by Serrano and colleagues. The authors have made a substantial effort to address the initial comments and to improve the manuscript's readability and accessibility. While we acknowledge these efforts, several initial concerns remain and some of the new analysis lead to additional questions. The authors propose that alternative exons respond only to a limited number of splicing factors. However, a recent study (Rogalska et al., Science 2024) estimated this number to be higher and earlier studies came to the same conclusion (e.g. Papasaikas et al., Mol Cell 2015), particularly when considering core components of the spliceosome. It is not clear where this discrepancy comes from. More importantly, even if indirect effects in the present study enable the prediction of altered splicing factor activity, it is not clear what can be learned from this. If the pipeline predicts that the activity of a splicing factor is altered based on (many) indirect targets, any functional consequence, for example in tumorigenesis, may also be indirect and not related to the activity of that particular splicing factor. Our concern of direct vs indirect effects therefore remains. Similarly, we continue to have reservations about the cancer analysis. While we appreciate that the wet lab validations we suggested might fall outside the scope of this study, we believe that such validations are essential to increase the manuscript's broader appeal and accessibility. For example, demonstrating experimentally changes in expression levels or PTMs in the identified splicing factors in a cancer model through Western blotting or proteomics would strongly substantiate the modelling approach. Without this evidence, the data remain too abstract to form a robust basis for new models and concepts and seem too specialized for a broader audience. Finally, the manuscript remains difficult to read, especially for a broad and more general audience and it is still difficult to follow the bioinformatics procedures in detail. Therefore, we suggest targeting a bioinformatics-focused journal, as this would align better with the manuscript's technical nature and present state of validation.", "reviewer_4": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts.", "reviewer_5": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts."}, "author_response": "We thank Reviewer 2 for their thoughtful and constructive feedback. Below, we address each point raised in detail. Major concern (singular since all comments below are about the validation): • Since there is only one observation, AUC/ROC is not adequate. AUROC can be very misleading when data is imbalanced (like it is here). Please use better metrics – with emphasis on the plural. • Seem that a ranking would do nicely as a metric here? If the right one is ranked 4 of 240 SF the ranking would be top 5/240 = top 2.08%. Such a metric would be much more informatice and easy to interpret. But feel free to explore others in suplementary figures. We thank the Reviewer for this insightful comment. Indeed, we initially evaluated performance using ROC AUC and the percentile rank of the perturbed splicing factor and observed that the two metrics were highly correlated (Spearman coefficient = 1; see Fig. 1). Given this redundancy and the familiarity of ROC AUC in the community, we initially opted to report only ROC AUC. Please note when reviewing this analysis that the scale is reversed from the one suggested by the Reviewer, so a better ranking would be 100%. However, in response to the Reviewer’s suggestion, we now include the ranking percentile metric in all relevant analyses (see manuscript Extended Data Fig. 1), which fully supports the same conclusions. Figure 1. Relationship between ROC AUC and Rank Percentile metrics to evaluate different methods for splicing factor activity estimation. Top, Spearman correlation and corresponding p-value. The high R value indicates that ROC-AUC serves as a comparable strategy compared to ranking for the evaluation of different approaches for splicing factor activity estimation. • The ENA data contains K562 and HepG2 data and is therefore not an independent dataset. Please divide the ENA data so that cell lines do not appear in both train/test data. • Specifically “train” on K562 and HepG2 data and apply to other cell lines • For the analysis where ENA (without K562 and HepG2 ) are held out you need to show the data in more detail. Specifically you need to show which cell lines and SF have poor AUCs and also annotate annotate tissue orign of cell types. This is the only way we can is if the method generalizes to other tissues/cell types. We appreciate the Reviewer’s attention to potential data leakage. Regarding this comment, we would like to clarify: - The ENA dataset includes samples from 54 projects unrelated to ENCODE. - The ENA dataset does not contain HepG2 samples, as verified in the “cell_line_name” field of Supplementary Table 2. In all benchmark settings involving empirical splicing factor→exon networks, the dataset used to generate the network is strictly excluded from the evaluation set. This ensures that network construction and activity evaluation occur on non-overlapping and independent datasets. To directly address the Reviewer’s request, we have added a new analysis in which splicing factor activities are estimated using networks derived only from ENCODE (HepG2 and K562), and evaluated non-ENCODE cell lines in ENA grouped by tissue of origin. We observe that VIPER consistently outperforms alternative methods when applied to cell types not used to construct the networks, supporting its cross-cellular context generalizability (median ROC-AUC VIPER by lineage: hematopoietic and lymphoid (n=13) = 1.0; breast (n=13) = 0.71; cervix (n=7) = 0.93; bone (n=4) = 0.74; kidney (n=3) = 0.96; lung (n=2) = 0.4; central nervous system (n=1) = 0.38; large intestine (n=1) = 1.0). These results are reported in the manuscript Supplementary Fig. 2b. • Please speculate why there are <0.5 AUCs both across SF and across datatypes We thank the Reviewer for highlighting this point. AUCs below 0.5 occur when the perturbed splicing factor ranks among the lowest predicted activities, in the opposite direction of the expected effects. There are several possible explanations: - Biological: splicing factors are known to be co-regulated; perturbation of one factor may induce stronger apparent effects on its regulators or downstream partners. - Technical: some networks may be more context-dependent or affected by experimental noise, reducing predictive accuracy in certain cases. While we acknowledge these limitations, we view this benchmark as a useful (though imperfect) approximation. Notably, when we expand the evaluation to real-world applications (Results: “Estimated splicing factor activity recapitulates regulatory mechanisms of increasing complexity”), our approach reveals functionally meaningful patterns missed by expression alone. We included this limitation in the Results section, “Splicing factor activity estimation is robust across studies and contexts”. • Should some SF be excluded from all downstream analysis since they are not trustworthy? As we did not find independent replicates for all perturbed splicing factors, this part of our validation is intended to provide a broad, comparative view of method performance across available perturbation datasets. Performance variability likely reflects a combination of dataset-specific factors and network quality, and we believe it is valuable to retain this heterogeneity in the benchmark. Therefore, we prefer not to exclude any splicing factors for downstream analyses. Even though we know their quality metrics might not be optimal, preserving them in the study can help in highlighting networks/pathways with no complete information to generate testable hypotheses and to extend knowledge on them. We also included this limitation in the Results section, “Splicing factor activity estimation is robust across studies and contexts”. • Remove the held-out ENCODE data (at least from the main figures) – non of them are indepdent and hence they are misleading. We appreciate this suggestion and understand the concern. However, while ENCODE and held-out ENCODE data may share cellular context, they originate from distinct experiments and provide the necessary independence to benchmark how generalizable our approach is across independent experiments in the same context. We have extended manuscript Figure 2c’s footer to clarify that this first evaluation of benchmark datasets offers an overview of the performance in experimentally independent benchmark datasets, and to point the readers to manuscript Supplementary Fig. 2b for the evaluation in independent cellular contexts from the ENA benchmark dataset. Minor comments (all about survival) • Survival analysis: Change from “high” and “low” to “increased” and “decreased”. A non-significant HR > 0 cannot be called \"High Risk\". We have updated manuscript Fig. 4 (panels e and d) and manuscript Supplementary Fig. 4 (panels b, e and f) to reflect this more accurate terminology. • Please also show a similar survival based plot but this time counting significant changes. The defintion of a significant association is that the HR is reliable - now you are overinterpreting point values with very high uncertainty (meaning the sign of the HR is not trustworthy). Thank you for this suggestion. We agree it is important not to over-interpret nonsignificant Cox Proportional Hazards coefficients. However, the purpose of disregarding whether coefficients are significantly different than zero in each particular cohort was to assess whether associations share the same trend across cohorts. To address this, we now provide binomial test results quantifying the consistency of direction effects across cohorts (p-value < 0.05 in all cases) as well as how many of these coefficients are statistically significant (FDR < 0.05) (see manuscript Fig. 4e). • We don’t quite understand why the authors do not replace the survival analysis with the new one (the one with co-factors) since we agree the conclusions are unchanged – but now much more reliable? We appreciate this suggestion and agree that the cofactor-adjusted model is more robust. Although the conclusions were unchanged, we have updated all survival analyses to use the improved model, as suggested (see manuscript Fig. 4e and f; Supplementary Fig. 4b, e and f). We look forward to read and updated manuscript. We thank the Reviewer once again for their insightful and constructive feedback, which helped us strengthen the manuscript further. We hope that the revisions fully address the concerns raised. Dimitrios S. Kanakoglou PhD Fellow Section for Bioinformatics, Health Technology Technical University of Denmark (DTU) Denmark And Kristoffer Vitting-Seerup, PhD Associate Professor Section for Bioinformatics, Health Technology Technical University of Denmark (DTU) Denmark"}, {"round": 3, "reviewer_comments": {"reviewer_6": "This manuscript describes the development of a new computational tool, building on previous work from the Serrano lab, to measure splicing factor activity from RNA-Seq data using delta PSI values. The authors demonstrate reliable prediction efficiency of the tool across various perturbed datasets and subsequently apply it to TCGA datasets, identifying numerous factors with potential oncogenic or tumor-suppressor activity. Although the authors report a correlation between predicted factor activity and patient prognosis, direct experimental evidence supporting the role of these newly identified factors in cancer phenotypes is lacking. The analysis is then extended to additional models, including a multistage carcinogenesis model and CCLE cell lines, correlating protein activity with proliferation markers. This is an interesting study, but I have two major concerns. First, the introduction and early results sections are written in a way that makes the work difficult to follow for a broad readership. In particular, the bioinformatic procedures are not clearly delivered in the results section. Second, as indicated above, the tool is not used to uncover novel mechanistic insights for proteins with high predicted splicing activity in the tested systems. This limits the scope of the work, nonetheless, it provides a valuable framework and has the potential to be useful to the field. Below I provide more specific comments: 1. The first section of the manuscript should clearly present all essential information so the readership can understand the procedures and tool development. For example, the authors should specify how delta PSI values were calculated (e.g., vast-tools), and whether only cassette exons were included or all alternative splicing types (retained introns, alternative splice sites, etc.). It should also be clarified whether 509 or 248 splicing factors were used and why they were selected, especially given that many of them are not splicing factors but chromatin modifiers or regulators of other aspects of RNA metabolism. Figures are cited in the introduction (e.g., Fig. 1 and Fig. 2a), but each should also be discussed in the results. The procedure for developing SF-Viper could be repeated briefly at the start of the results to aid accessibility. The frequent use of “benchmarking” is unclear, as it is not always stated relative to what. ROC curves are central to tool evaluation, yet definitions of true and false positives are missing, which would help non-experts follow. In general, the manuscript would benefit from simplifying overly long sentences (e.g., one paragraph consists of a single 5-line sentence) and from input by non-expert colleagues to improve clarity. While the methods are well written, the essentials of the methodology should be understandable without constant cross-referencing to methods, and notes. 2. The authors state that “each exon is regulated by a small number of splicing factors (median = 2).” This statement requires clarification. How is regulation defined? Are the exons never affected excluded from this analysis? If only exons regulated by at least one SF are used for the analysis (to remove always included exons and focus on “regulated” exons), does the median still remain 2? 3. Figure 2e lacks labeled axes. 4. Although this point may stem from addressing previous reviewer comments, I would emphasize that direct RNA binding is not the only way a splicing factor can (directly) influence exon inclusion. Many bona fide splicing factors do not necessarily bind RNA directly but, when perturbed, can alter spliceosome stoichiometry and thereby affect specific splicing programs. This has been demonstrated in several studies, some of which are already cited in the rebuttal. Similarly, chromatin-related factors can influence splicing indirectly by regulating the expression of splicing regulators, or potentially directly through chromatin modifications that impact transcription kinetics and the recruitment of splicing regulators. Therefore, I think that the CLIP section is somewhat confusing and the term “indirect” shouldn’t be used just due the lack of a CLIP signal. 5. The indisulam experiment is an interesting validation. The authors should consider citing the original Science manuscript (PMID: 28302793) and the Nature Communications study (PMID: 37666821), which describe the mechanism by which RBM39 autoregulates its expression. 6. The simultaneous knock-down experiments are intriguing but not sufficiently controlled. Does VIPER identify multiple perturbed splicing factors because they are truly perturbed, or because they produce highly similar splicing changes? Figure 3d shows additional gray dots near the CLK factors, presumably unperturbed CLKs with similar splicing effects. Similarly, SCAF4 and SCAF8 are paralogs, and NCBP1 and ZCCHC8 are known to physically associate. The authors should include single knock-down datasets and indicate the position of the unperturbed counterparts in the plots. 7. The authors refer to CPSF6 and DNMT1 as splicing factors. CPSF6 is a component of the cleavage and polyadenylation machinery and DNMT1 is a DNA methyltransferase. Their role in splicing is likely indirect (although direct effects cannot be excluded). Unless direct evidence is provided by previously published literature, these should not be referred to as “splicing factors” and their possible indirect effects should be clearly stated.", "reviewer_7": "The authors have addressed all concerns raised by reviewers 2 and 3, or provided comprehensive and detailed explanations. The previously introduced changes were valuable and improved the benchmarking strategy."}, "author_response": "We thank Reviewer #6 for the thorough revision and helpful comments. We have revised the Introduction and early Results sections to improve clarity and accessibility, with clearer descriptions of the bioinformatic workflow. Please find below our point-by-point responses to the comments. Below I provide more specific comments: 1. The first section of the manuscript should clearly present all essential information so the readership can understand the procedures and tool development. For example, the authors should specify how delta PSI values were calculated (e.g., vast-tools), and whether only cassette exons were included or all alternative splicing types (retained introns, alternative splice sites, etc.). It should also be clarified whether 509 or 248 splicing factors were used and why they were selected, especially given that many of them are not splicing factors but chromatin modifiers or regulators of other aspects of RNA metabolism. Figures are cited in the introduction (e.g., Fig. 1 and Fig. 2a), but each should also be discussed in the results. The procedure for developing SF-Viper could be repeated briefly at the start of the results to aid accessibility. The frequent use of “benchmarking” is unclear, as it is not always stated relative to what. ROC curves are central to tool evaluation, yet definitions of true and false positives are missing, which would help non-experts follow. In general, the manuscript would benefit from simplifying overly long sentences (e.g., one paragraph consists of a single 5-line sentence) and from input by non-expert colleagues to improve clarity. While the methods are well written, the essentials of the methodology should be understandable without constant cross-referencing to methods, and notes. We thank the Reviewer for the thorough and constructive assessment. We agree that some sections have become unclear after multiple revisions. Following the Reviewer’s suggestions, we have revised the section “Empirical networks outperform computational ones for activity estimation” to improve clarity and accessibility. Specifically: - We now specify that Delta PSI values were calculated from PSI values computed with vast-tools. - We clarify that only cassette exons were included, as they provide a robust and generalizable signature of splicing factor activity. Including additional types of splicing events is beyond the scope of this study; however, we note in the Discussion that this also represents an interesting direction for future work. This study serves as an initial exploration of splicing factor activity, which we hope will be further developed by the field. - We explain that we curated a list of 509 “splicing factors” due to the lack of a standardized definition. Of these, 248 had available RNA-seq perturbation datasets, which were used to evaluate the method. This point is now stated in Discussion section as well. - Figure 2a is now also cited and discussed in this sections. In particular, Figure 1 illustrates the diverse regulatory mechanisms that can influence splicing factor activity, and it is referenced again in the Results section “Estimated splicing factor activity recapitulates regulatory mechanisms of increasing complexity”. - We have briefly reiterated the VIPER procedure at the start of the Results to aid accessibility and contextualize its adaptation for splicing factors. - We have tried to clarify the use of the term “benchmarking”; we now specify what is being benchmarked every time. - We have explained that true and false positives used to compute ROC-AUCs are defined as the experiment where a splicing factor has been perturbed (true) or not (false) across the compendium of uniformly-processed experiments. 2. The authors state that “each exon is regulated by a small number of splicing factors (median = 2).” This statement requires clarification. How is regulation defined? Are the exons never affected excluded from this analysis? If only exons regulated by at least one SF are used for the analysis (to remove always included exons and focus on “regulated” exons), does the median still remain 2? In this analysis, related to Fig. 2d, we consider the regulatory interactions defined in empirical splicing factor networks. In these networks, we define as a target exon all those exons showing an absolute Delta PSI (|Delta PSI|) of at least 15% upon perturbation of a splicing factor. Exons that do not exhibit such changes upon perturbation were excluded from the analysis. Accordingly, when considering only exons regulated by at least one splicing factor in our empirical networks, the median number of regulating factors remains 2. 3. Figure 2e lacks labeled axes. The Y-axis in Figure 2e represented the same metric as in the left panel of Figure 2d. We have now added the corresponding label to Figure 2e for clarity. 4. Although this point may stem from addressing previous reviewer comments, I would emphasize that direct RNA binding is not the only way a splicing factor can (directly) influence exon inclusion. Many bona fide splicing factors do not necessarily bind RNA directly but, when perturbed, can alter spliceosome stoichiometry and thereby affect specific splicing programs. This has been demonstrated in several studies, some of which are already cited in the rebuttal. Similarly, chromatin-related factors can influence splicing indirectly by regulating the expression of splicing regulators, or potentially directly through chromatin modifications that impact transcription kinetics and the recruitment of splicing regulators. Therefore, I think that the CLIP section is somewhat confusing and the term “indirect” shouldn’t be used just due the lack of a CLIP signal. We thank the reviewer for this important clarification. Indeed, a previous Reviewer suggested distinguishing between direct and indirect effects using CLIP-seq data, and we acknowledge that our wording was misleading. We have now revised the text to clarify that direct regulation can occur through multiple mechanisms, not only via RNA binding. CLIP-seq data allow us to assess one such mechanism (direct RNA binding) versus exon inclusion signatures. We have rewritten the Results section “Empirical networks capture functional regulation” and Supplementary Note 5 correspondingly. 5. The indisulam experiment is an interesting validation. The authors should consider citing the original Science manuscript (PMID: 28302793) and the Nature Communications study (PMID: 37666821), which describe the mechanism by which RBM39 autoregulates its expression. We thank the Reviewer for this suggestion and apologize for the earlier omission. We have now cited both studies in the relevant section discussing RBM39 autoregulation. 6. The simultaneous knock-down experiments are intriguing but not sufficiently controlled. Does VIPER identify multiple perturbed splicing factors because they are truly perturbed, or because they produce highly similar splicing changes? Figure 3d shows additional gray dots near the CLK factors, presumably unperturbed CLKs with similar splicing effects. Similarly, SCAF4 and SCAF8 are paralogs, and NCBP1 and ZCCHC8 are known to physically associate. The authors should include single knock-down datasets and indicate the position of the unperturbed counterparts in the plots. We thank the Reviewer for raising this point. Our activity definition indeed implies that functionally related splicing factors can show similar exon inclusion signatures, which could be consistent with the behavior of paralogs (SCAF4/SCAF8) or physically associated factors, including NCBP1 and ZCCHC8, or the CLK family. The combinatorial perturbation datasets we analyzed were originally designed to study such interactions, and therefore are not ideal (compared with perturbations of unrelated pairs) for assessing how well activity estimation distinguished true co-perturbation from similarity in splicing effects. Although these datasets include single knockdowns, we cannot use them directly to validate our approach because they were also used to construct the empirical splicing factor→exon networks; evaluating the networks with the same signatures used to infer them would be circular. Following the Reviewer’s advice, we assessed, within each single-factor perturbation experiment, the inferred activities of the other splicing factors later perturbed in combination. These results are now shown in Supplementary Fig. 3c, with the relevant splicing factors labeled. Some splicing factors cannot be assessed when their target exons lack sufficient coverage. For example, in PRJNA587741 (prostate-derived cell line), ESRP1 and KHDRBS3 cannot be evaluated in each other’s single perturbations for this reason. Nevertheless, among the splicing factors for which activities can be estimated, except for NCBP1 upon perturbing its interactor ZCCHC8, the splicing factors showing the highest activity in the combinatorial experiments (Fig. 3d) do not rank highly in the corresponding single-perturbation experiments (Supplementary Fig. 3c). Despite the limitations of these experimental designs, this suggests that VIPER identifies truly perturbed splicing factors rather than simply those with similar targets. As requested, we now also indicate the positions of the corresponding unperturbed genes in Figure 3d to visualize their relationship to the perturbed splicing factors. This clarification has been added to the Results section: “Estimated splicing factor activity recapitulates regulatory mechanisms of increasing complexity”. 7. The authors refer to CPSF6 and DNMT1 as splicing factors. CPSF6 is a component of the cleavage and polyadenylation machinery and DNMT1 is a DNA methyltransferase. Their role in splicing is likely indirect (although direct effects cannot be excluded). Unless direct evidence is provided by previously published literature, these should not be referred to as “splicing factors” and their possible indirect effects should be clearly stated. We thank the Reviewer for this important clarification and agree that the list of splicing factors may not be perfect. We have now added a note in the Discussion addressing this limitation. To investigate these specific cases, we reviewed the sources used to compile our list. CPSF6 is annotated as a splicing factor in 3 of the 5 reference studies1–5, as well as in SpliceosomeDB and the Gene Ontology (GO) RNA-binding protein term, though not under the GO splicing factor term. DNMT1 appears in Papasaikas et al.1, and in SpliceosomeDB. In Papasaikas et al., the authors analyzed “a fraction” of the proteins defined by Wahl et al.6, as proteins in highly purified spliceosomal complexes assembled on model single-intron pre-mRNAs. We now explicitly acknowledge in the Discussion that some factors in our list, such as CPSF6 and DNMT1, may influence splicing indirectly, given the available evidence. At the same time, the consistently strong performance of splicing factor activity estimation across our evaluations suggests that exon inclusion signatures could also serve more broadly as a functional readout to identify genes contributing to a molecular phenotype. We now highlight this in the revised Discussion as well."}, {"round": 4, "reviewer_comments": {"reviewer_6": "The authors have improved the clarity of the initial section and have satisfactorily addressed all of my previous comments. I am pleased to recommend the manuscript for publication in Nature Communications."}, "author_response": ""}]

Decision letter