Functional Splicing Network Reveals Extensive Regulatory Potential of the Core Spliceosomal Machinery

Pre-mRNA splicing relies on the poorly understood dynamic interplay between>150 protein components of the spliceosome. The steps at which splicing can be regulated remain largely unknown. We systematically analyzed the effect of knocking down the components of the splicing machinery on alternative splicing events relevant for cell proliferation and apoptosis and used this information to reconstruct a network of functional interactions. The network accurately captures known physical and functional associations and identifies new ones, revealing remarkable regulatory potential of core spliceosomal components, related to the order and duration of their recruitment during spliceosome assembly. In contrast with standard models of regulation at early steps of splice site recognition, factors involved in catalytic activation of the spliceosome display regulatory properties. The network also sheds light on the antagonism between hnRNP C and U2AF, and on targets of antitumor drugs, and can be widely used to identify mechanisms of splicing regulation.


INTRODUCTION
Pre-mRNA splicing is carried out by the spliceosome, one of the most complex molecular machineries of the cell, composed of five small nuclear ribonucleoprotein particles (U1, U2, U4/5/6 snRNP) and about 150 additional polypeptides (reviewed by Wahl et al., 2009). Detailed biochemical studies using a small number of model introns have delineated a sequential pathway for the assembly of spliceosomal subcomplexes. For example, U1 snRNP recognizes the 5' splice site, and U2AF-a heterodimer of 35 and 65 kDa subunits-recognizes sequences at the 3' end of introns. U2AF binding helps to recruit U2 snRNP to the upstream branchpoint sequence, forming complex A. U2 snRNP binding involves interactions of pre-mRNA sequences with U2 snRNA as well as with U2 proteins (e.g., SF3B1). Subse-quent binding of preassembled U4/5/6 tri-snRNP forms complex B, which after a series of conformational changes forms complexes Bact and C, concomitant with the activation of the two catalytic steps that generate splicing intermediates and products. Transition between spliceosomal subcomplexes involves profound dynamic changes in protein composition as well as extensive rearrangements of base-pairing interactions between snRNAs and between snRNAs and splice site sequences (Wahl et al., 2009). RNA structures contributed by base-pairing interactions between U2 and U6 snRNAs serve to coordinate metal ions critical for splicing catalysis (Fica et al., 2013), implying that the spliceosome is an RNA enzyme whose catalytic center is only established upon assembly of its individual components.
Differential selection of alternative splice sites in a pre-mRNA (alternative splicing, AS) is a prevalent mode of gene regulation in multicellular organisms, often subject to developmental regulation (Nilsen and Graveley, 2010) and frequently altered in disease (Cooper et al., 2009;Bonnal et al., 2012). Substantial efforts made to dissect mechanisms of AS regulation on a relatively small number of pre-mRNAs have provided a consensus picture in which protein factors recognizing cognate auxiliary sequences in the pre-mRNA promote or inhibit early events in spliceosome assembly (Fu and Ares, 2014). These regulatory factors include members of the hnRNP and SR protein families, which often display cooperative or antagonistic functions depending on the position of their binding sites relative to the regulated splice sites. Despite important progress (Barash et al., 2010;Zhang et al., 2010), the combinatorial nature of these contributions complicates the formulation of integrative models for AS regulation.
To systematically capture the contribution of splicing regulatory factors, including components of the core spliceosome (Clark et al., 2002;Park et al., 2004;Pleiss et al., 2007;Saltzman et al., 2011) to AS regulation, we set up a high-throughput screen to evaluate the effects of knocking down each individual splicing component or regulator on 36 functionally important AS events (ASEs) and developed a framework for data analysis and network modeling. Network-based approaches can provide rich representations for real-world systems that capture their organization more adequately than more traditional approaches or mere cataloguing of pairwise relationships. Numerous studies have utilized them as a platform for deriving comprehensive transcriptional, metabolic, and physical interaction maps in several organisms and contexts (e.g., Chatr-Aryamontri et Franceschini et al., 2013;Wong et al., 2012). Their application has also proven invaluable for deciphering regulatory relationships in transcriptional circuits, for studying their properties upon perturbation, and for connecting physiological responses and diseases to their molecular underpinnings (e.g., Kim et al., 2012;Watson et al., 2013;Yang et al., 2014).
Key to our approach is the premise that the distinct profiles of splicing changes caused by perturbation of regulatory factors depend on the functional bearings of those factors and can, therefore, be used as proxies for inferring functional relationships. We derive a network that recapitulates the topology of known splicing complexes and provides an extensive map of >500 functional associations among ~200 splicing factors (SFs). We use this network to identify general and particular mechanisms of AS regulation and to identify key SFs that mediate the effects on AS of perturbations induced by iron (see accompanying manuscript by Tejedor et al., 2014, in this issue of Molecular Cell) or by antitumor drugs targeting components of the splicing machinery. Our study offers a unique compendium of functional interactions among SFs, a discovery tool for coupling cell-perturbing stimuli to splicing regulation, and an expansive view of the splicing regulatory landscape and its organization.

RESULTS
Results from a genome-wide siRNA screen to identify regulators of Fas/CD95 AS revealed that knockdown of a significant fraction of core SFs (defined as proteins identified in highly purified spliceosomal complexes assembled on model, single-intron pre-mRNAs; Wahl et al., 2009) caused changes in Fas/CD95 exon 6 inclusion (Tejedor et al., 2014). The number of core factors involved and the extent of their regulatory effects exceeded those of classical splicing regulatory factors like SR proteins or hnRNPs. To systematically evaluate the contribution of different classes of splicing regulators to AS regulation, we set up a screen in which we assessed the effects of knockdown of each individual factor on 38 ASEs relevant for cell proliferation and/or apoptosis ( Figures 1A and 1B). The ASEs were selected due to their clear impact on protein function and their documented biological relevance (see Table S1 available online). A library of 270 siRNA pools (Table S2) was used, corresponding to genes encoding core spliceosomal components as well as auxiliary regulatory SFs and factors involved in other RNA-processing steps, including RNA stability, export, or polyadenylation. In addition, 40 genes involved in chromatin structure mod-ulation were included, to probe for possible functional links between chromatin and splicing regulation (Luco et al., 2011).
The siRNA pools were individually transfected in biological triplicates in HeLa cells using a robotized procedure (Supplemental Experimental Procedures). Seventy-two hours posttransfection, RNA was isolated, cDNA libraries generated using oligo-dT and random primers, and the patterns of splicing probed by PCR using primers flanking each of the alternatively spliced regions ( Figure 1A). PCR products were resolved by highthroughput capillary electrophoresis (HTCE) and the ratio between isoforms calculated as the percent spliced in (PSI) index (Katz et al., 2010), which represents the percentage of exon inclusion/alternative splice site usage. Knockdowns of factors affecting one ASE (Fas/CD95) were also validated with siRNA pools from a second library and by a second independent technique (Illumina HiSeq sequencing; Tejedor et al., 2014). Conversely, the effects of knocking down core SFs (P14, SNRPG, and SF3B1) on the ASEs were robust when different siRNAs and RNA purification methods were used ( Figure S1A). Several pieces of evidence indicate that the changes in AS are not due to the induction of cell death upon depletion of essential SFs. First, cell viability was not generally compromised after 72 hr of knockdown of 13 individual core SFs ( Figure S1C). Second, cells detached from the plate were washed away before RNA isolation. Third, splicing changes upon induction of apoptosis with staurosporine were distinct from those observed upon knockdown of SFs ( Figure S1D).
Overall the screen generated a total of 33,288 PCR data points. These measurements provide highly robust and sensitive estimates of the relative use of competing splice sites. Figure 1D shows the spread of the effects of knockdown for all the factors in the screen for the ASEs analyzed. Three events (VEGFA, SYK, and CCND1) were excluded from further analysis, because they were not affected in the majority of the knockdowns (median absolute DPSI <1, Figure 1D). From these data we generated, for every knockdown condition, a perturbation profile that reflects its impact across the 35 ASEs in terms of the magnitude and direction of the change. Figures 1E and 1F show examples of such profiles and the relationships between them: in one case, knockdown of CDC5L or PLRG1 generates very similar perturbation profiles (upper diagrams in Figures 1E and 1F), consistent with their well-known physical interactions within the PRP19 complex (Makarova et al., 2004). In contrast, the profiles associated with the knockdown of SNRPG and DDX52 are to a large extent opposite to each other (lower diagrams in Figures 1E and 1F), suggesting antagonistic functions. Perturbation profiles and relationships between factors were not significantly changed   . Information about ASEs, (a) type (ce, cassette exon; 5ss, alternative 5' splice site; me, mutually exclusive exon; complex, multiexonic rearrangement; 3ss, alternative 3' splice site; ir, intron retention) and (b) involvement in apoptosis and/or cell proliferation regulation, as well as information on the category of genes knocked down, is color coded. (B) Boxplot representation of the mean absolute Z score of AS changes induced by the knockdown of particular classes of factors, including core and noncore SFs, other RNA processing factors, and chromatin remodeling factors. Median and spread (interquartile range) of AS changes of mock siRNA conditions are represented as a green line and box, respectively. (C) As in (B) for spliceosomal factors that assemble early and stay as detectable components through the spliceosome cycle; factors that are present only transiently at early, mid, or late stages of assembly; and components of the exon junction complex (EJC). See also Figure S3, Table S2, Table S3, and Table S4. when knockdown of a subset of factors was carried out for 48 hr ( Figure S2A), arguing against AS changes being the consequence of secondary effects of SF knockdown in these cases (see also below the discussion on IK/SMU1). Pervasive and Distinct Effects of Core Spliceosome Components on Alternative Splicing Regulation The output of the screening process is summarized in the heatmap of Figure 2A. Three main conclusions can be derived from these results. First, a large fraction of the SF knockdowns have noticeable but distinct effects on AS, indicating that knockdown of many components of the spliceosome causes switches in splice site selection, rather than a generalized inhibition of splicing of every intron. Second, a substantial fraction of AS changes correspond to higher levels of alternative exon inclusion, again arguing against simple effects of decreasing splicing activity upon knockdown of general SFs, which would typically favor skipping of alternative exons that often harbor weaker splice sites. These observations suggest extensive versatility of the effects of modulating the levels of SFs on splice site 5 selection. Third, despite the diversity of the effects of SF knockdowns on multiple ASEs, similarities can also be drawn (Figure 2A). For example, a number of ASEs relevant for the control of programmed cell death (e.g., Fas, APAF1, CASP9, or MCL1) cluster together and are similarly regulated by a common set of core spliceosome components. This observation suggests coordinated regulation of apoptosis by core factors, although the functional effects of these AS changes appear to be complex ( Figure S3A). Another example is two distinct ASEs in the fibronectin gene (cassette exons EDA and EDB, of paramount importance to control aspects of development and cancer progression; Muro et al., 2003) which cluster close together, suggesting common mechanisms of regulation of the ASEs in this locus.
To explore the effects of different classes of SFs on AS regulation, we first grouped them in four categories: genes coding for core spliceosome components, noncore SFs/regulators, RNA-processing factors not directly related to splicing, and chromatin-related factors. Core SFs display the greater spread and average magnitude of effects, followed by noncore SFs (Figure 2B). Factors related with chromatin structure and remodeling showed a narrower range of milder effects ( Figure 2B), although their values were clearly above the range of changes observed in mock knockdown samples. Of interest, knockdown of core factors leads to stronger exon skipping effects, while noncore SFs and other RNA processing and chromatin-related factors show more balanced effects ( Figures S3B and S3C).
To further dissect the effects of core components, we subdivided them into four categories depending on the timing and duration of their recruitment to splicing complexes during spliceosome assembly (Wahl et al., 2009): (1) persistent components (e.g., 17S U2 snRNPs, Sm proteins) that enter the assembly process before or during complex A formation and remain until completion of the reaction, (2) transient early components that join the reaction prior to B complex activation but are scarce at later stages of the process, (3) transient middle components joining during B-act complex formation, and (4) transient late components that are only present during or after C complex formation. Our results indicate that knockdowns of factors that persist during spliceosome assembly cause AS changes of higher magnitude than those caused by knockdowns of transient factors involved in early spliceosomal complexes, which in turn are stronger and more dispersed than those caused by factors involved in middle complexes and in complex C formation or catalytic activation ( Figure 2C). As observed for core factors in general, knockdown of persistent factors tends to favor skipping, while the effects of knockdown of more transient factors are more split between inclusion and skipping (Figures S3B and  S3C).
Interestingly, knockdowns of components of the exon junction complex (EJC) display a range of effects that resembles those of transient early or mid splicing complexes ( Figure 2C). These results are in line with recent reports documenting effects of EJC components in AS, in addition to their standard function in postsplicing processes (Ashton-Beaucage et al., 2010).
To test whether the knockdown of core factors causes a general inhibition of splicing, we measured the levels of introns (both in constitutive and alternatively spliced regions) relative to exons in four genes, upon the knockdown of four SFs (BRR2, SF3B1, SNRPG, and SLU7). The results shown in Figure 3 reveal a complex picture in which retention of certain introns is clearly increased and retention of other introns is very mild, and there is even one case (the two introns flanking the fibronectin EDI alternative exon) in which the low levels of intron detected under control conditions are reduced upon knockdown of core SFs (notably, this ASE typically displays a distinct pattern of AS upon knockdown of core SFs, Figure 2A). Therefore, rather than a general uniform inhibition of splicing, reduction in the levels of core SFs induces differential effects in different introns, a concept reminiscent of the differential effects observed for antitumor drugs targeting core components like SF3B1 (Bonnal et al., 2012).

Reconstruction of a Functional Network of Splicing Factors
To systematically and accurately map the functional relationships of SFs on AS regulation, we quantified the similarity between AS perturbation profiles for every pair of factors in our screen using a robust correlation estimate for the effects of the factors' knockdowns across the ASEs. This measure captures the congruence between the shapes of perturbation profiles, while it does not take into account proportional differences in the magnitude of the fluctuations. Crucially, it discriminates between biological and technical outliers and is resistant to distorting effects of the latter (Supplemental Experimental Procedures).
We next employed glasso, a regularization-based algorithm for graphical model selection Supplemental Experimental Procedures), to reconstruct a network from these correlation estimates ( Figures 4A and S4). Glasso seeks a parsimonious network model for the observed correlations. This is achieved by specifying a regularization parameter, which serves as a penalty that can be tuned in order to control the number of inferred connections (network sparsity) and therefore the false discovery rate (FDR) of the final model. For a given regularization parameter, both the number of inferred connections and the FDR are a decreasing function of the number of ASEs assayed ( Figure 4B). Random sampling with different subsets of the real or reshuffled versions of the data indicates that network reconstruction converges to ~500 connections with a FDR <5% near 35 events ( Figure 4B), implying that analysis of additional ASEs in the screening would have only small effects on network sparsity and accuracy.
The complete network of functional interactions among the screened factors is shown in Figure 4A. It is comprised of 196 nodes representing screened factors and 541 connections (average degree 5.5), of which 518 correspond to positive and 23 to negative functional associations. The topology of the resulting network has several key features. First, two classes of factors are easily distinguishable. The first class encompasses factors that form densely connected clusters comprised of core spliceosomal components ( Figure 4A). Factors physically linked within U2 snRNP and factors functionally related with U2 snRNP activity, and with the transition from complex A to B, form a tight cluster (orange nodes in the inset of Figure 4A). An adjacent but distinct highly linked cluster corresponds mainly to factors physically associated within U5 snRNP or the U4/5/6  tri-snRNP (blue nodes in Figure 4A, inset). Factors in these clusters and their mutual links encompass 23% of the network nodes but 63% of the total inferred functional associations (average degree 14.5). A second category includes factors that form a periphery of lower connectivity that occasionally projects to the core. This category includes many of the classical splicing regulators (e.g., SR proteins, hnRNPs) as well as several chromatin factors. An intermediate category-more densely interconnected but not reaching the density of interactions of the central modules-includes multiple members of the RNAdependent DEAD/X box helicase family. While the number of connections (degree) is significantly higher for core spliceosomal factors in general ( Figure 4C), this is especially clear for factors that represent persistent components of splicing complexes along the assembly pathway ( Figures 4D and 4E). Persistent factors are particularly well connected to each other (50% of the possible functional connections are actually detected, Figure 4E). Factors that belong to consecutive complexes are significantly more linked to each other than factors that belong to early and late complexes, and factors that transiently assemble toward the final stages of spliceosome assembly display the least number of connections ( Figures 4D and 4E).
Collectively, these results indicate that core SFs, particularly those that persist during assembly, display highly related functions in splice site recognition and fluctuations in their levels or activities have related consequences on the modulation of splice site choice. Of the observed connections, 20% can be attributed to known physical interactions and 82% of these are among core components ( Figure 4A). Examples include tight links between the two subunits of U2AF, between the interacting partners PLRG1 and CDC5L or IK and SMU1. We estimate that about 50% of the functional links detected in our network could be predicted by previous knowledge of the composition or function of splicing complexes.
Of relevance, a substantial number of the functional associations closely recapitulate the composition and even-to some extent-the topology of known spliceosomal complexes (Figure S4; see also below and Discussion). These observations suggest that differences in the sensitivity of ASEs to SF perturbations reflect the specific role of these factors in the splicing process. Conversely, these results warrant the use of splicing perturbation profiles as proxies for inferring physical or functional links between factors in the splicing process.

Generality of the Functional Interactions in the AS Regulation Network
The above network was derived from variations in 35 ASEs selected because of their relevance for cell proliferation and apoptosis, but they represent a heterogeneous collection of AS types and, possibly, regulatory mechanisms. To identify functional connections persistent across AS types, we undertook a subsampling approach. We sampled subsets of 17 from the original 35 ASEs and iteratively reconstructed the network of functional associations from each subset (see Supplementary methods). Following 10,000 iterations, we asked how many functional connections were recovered in at least 90% of the sampled networks. The retrieved set of connections are shown in Figure 5A and they can be considered to represent a ''basal'' or indispensable splicing circuitry that captures close similarities in function between SFs regardless of the subset of splicing substrates considered. The majority of these connections correspond to links between core components, with modules corresponding to major spliceosomal complexes. For example, the majority of components of U2 snRNP constitute the most prominent of the core modules ( Figure 5A). Strikingly, the module is topologically subdivided in two submodules, the first corresponding to components of the SF3a and SF3b complexes and the second corresponding to proteins in the Sm complex. The connectivity between U2 snRNP components remarkably recapitulates known topological features of U2 snRNP organization, including the assembly of SF3a/b complexes in stem loops I and II and the assembly of the Sm ring in a different region of the U2 snRNA (Behrens et al., 1993) (Figures 5A and S4). That the connectivity between components derived from assessing the effects on AS of depletion of these factors correlates with topological features of the organization of the snRNP further argues that the functions in splice site selection are tightly linked to the structure and function of this particle, highlighting the potential resolution of our approach for identifying meaningful structural and mechanistic links.
A second prominent ''core'' module corresponds to factors known to play roles in conformational changes previous to/ concomitant with catalysis in spliceosomes from yeast to human, including PRP8, PRP31 or U5-200. Once again some of the topological features of the module are compatible with physical associations between these factors (e.g., PRP8 with PRP31 or U5-116K or PRP31 with C20ORF14). That knockdown of these late-acting factors, which coordinate the final steps of the splicing process (Bottner et al., 2005;Hä cker et al., 2008;Wahl et al., 2009), causes coherent changes in splice site selection strongly argues that the complex process of spliceosome assembly can be modulated at late steps, or even at the time of catalysis, to affect splice site choice. Other links in the ''core'' network recapitulate physical interactions, including the aforementioned modules involving the two subunits of the 3' splice site-recognizing factor U2AF, the interacting partners PLRG1 and CDC5L and IK and SMU1.
To further evaluate the generality of our approach, we tested whether the functional concordance derived from our network analysis and from the iterative selection of a ''core'' network could be recapitulated in a different cell context and genomewide. We focused our attention in the strong functional association between IK and SMU1, two factors that have been described as associated with complex B and its transition to Bact (Bessonov et al., 2008) ( Figure 5A). Analysis of the effects of knockdown of these factors in the same sets of ASEs in HEK293 cells revealed a similar set of associations, despite the fact that individual events display different splicing ratios upon IK/SMU1 knockdown in this cell line compared to HeLa ( Figures 5B and 5C). To explore the validity of this functional association in a larger set of splicing events, RNA was isolated in biological triplicates from HeLa cells transfected with siRNAs against these factors, and changes in AS were assessed using genome-wide splicing-sensitive microarrays. The results revealed a striking overlap between the effects of IK and SMU1 knockdown, both on gene expression changes and on AS changes (61% overlap, p < 1E-20) covering a wide range of gene expression levels, fold differences in isoform ratios and AS types ( Figures 5D-5F, S5B, and S5C). The overlap was much more limited with RBM6, a splicing regulator  located elsewhere in the network. These results validate the strong functional similarities between IK and SMU1 detected in our analysis. The molecular basis for this link could be explained by the depletion of each protein (but not their mRNA) upon knockdown of the other factor, consistent with formation of a heterodimer of the two proteins and destabilization of one partner after depletion of the other ( Figure S5A). Notably, the best correlation between the effects of these factors was between IK knockdown for 48 hr and SMU1 knockdown for 72 hr, perhaps reflecting a stronger functional effect of IK on AS regulation ( Figures S2B and S2D).
Gene ontology analysis revealed a clear common enrichment of AS changes in genes involved in cell death and survival (Figure S5D), suggesting that IK/SMU1 could play a role in the control of programmed cell death through AS. To evaluate this possibility, we tested the effects of knockdown of these factors, individually or combined, on cell proliferation and apoptotic assays. The results indicated that their depletion activated cleavage of PARP, substantially reduced cell growth, and increased the fraction of cells undergoing apoptosis ( Figures S5E-S5G). We conclude that functional links revealed by the network can be used to infer common mechanisms of splicing regulation and provide insights into the biological framework of the regulatory circuits involved.

Ancillary Functional Network Interactions Reveal Alternative Mechanisms of AS Regulation
The iterative generation of networks from subsets of ASEs explained above could, in addition to identifying general functional interactions, be used to capture functional links that are specific to particular classes of events characterized by distinct regulatory mechanisms. To explore this possibility, we sought the set of functional connections that are robustly recovered in a significant fraction of ASEs subsets but are absent in others (see Supplemental Experimental Procedures). Figure 6A shows the set of such identified interactions, which therefore represents specific functional links that only emerge when particular subsets of ASEs are considered. The inset represents graphically the results of principal component analysis (PCA) showing the variability in the strength of these interactions depending on the presence or absence of different ASEs across the subsamples. The most discriminatory set of ASEs for separating these interactions is also shown (see also Figure S6). While this analysis may be constrained by the limited number of events analyzed and higher FDR due to multiple testing, it can provide a valuable tool for discovering regulatory mechanisms underlying specific classes of ASEs as illustrated by the example below.
A synergistic link was identified between hnRNP C and U2AF1 (the 35 kDa subunit of U2AF, which recognizes the AG dinucleotide at 3' splice sites) for a subset of ASEs ( Figures 6A and S6A). We explored a possible relationship between the effects of knocking down these factors and the presence of consensus hnRNP C binding sites (characterized by stretches of uridine residues) or 3' splice site-like motifs in the vicinity of the regulated exons. The results revealed a significant correlation between the two factors (0.795), associated with the presence of composite elements comprised of putative hnRNP C binding sites within sequences conforming to a 3' splice site motif (see Experimental Procedures) upstream and/or downstream of the regulated exons ( Figures 6B and S6B). No such correlation was observed for ASEs lacking these sequence features ( Figures 6C and S6C). These observations are in line with results from a recent report revealing that hnRNP C prevents exonization of Alu elements by competing the binding of U2AF to Alu sequences (Zarnack et al., 2013). Our results extend these findings by capturing functional relationships between hnRNP C and U2AF in ASEs that harbor binding sites for these factors in the vicinity of the regulated splice sites, both inside of Alu elements or independent of them ( Figures 6B and 6D). Our data are compatible with a model in which uridine-rich and 3' splice site-like sequences sequester U2AF away from the regulated splice sites, causing alternative exon skipping, while hnRNP C displaces U2AF away from these decoy sites, facilitating exon inclusion. In this model, depletion of U2AF and hnRNP C are expected to display the same effects in these ASEs, as observed in our data ( Figures  6B and 6D). In contrast, in ASEs in which hnRNP C binding sites are located within the polypyrimidine tract of the 3' splice site of the regulated exon, the two factors display antagonistic effects ( Figures 6C and 6D), as expected from direct competition between them at a functional splice site. This example illustrates the potential of the network approach to identify molecular mechanisms of regulation on the basis of specific sequence features and the interplay between cognate factors. Positive or negative functional correlations are represented by green or red edges, respectively. Edge thickness signifies the strength of the functional interaction, while node size is proportional to the overall impact (median Z score) of a given knockdown in the regulation of AS. Node coloring depicts the network's natural separation in coherent modules (see Supplemental Experimental Procedures). Known physical interactions as reported in the STRING database (Franceschini et al., 2013) are represented in black dotted lines between factors. The inset is an expanded view of factors in the U2 and U4/5/6 snRNP complexes.
(B) Number of recovered network edges using actual or randomized data sets (red and green points respectively) as a function of the number of ASEs used. Lines represent the best fit curves for the points.  Figure S4E and Table S5.   Mapping Drug Targets within the Spliceosome Another application of our network analysis is to identify possible targets within the splicing machinery of physiological or external (e.g., pharmacological) perturbations that produce changes in AS that can be measured using the same experimental setting.
The similarity between profiles of AS changes induced by such perturbations and the knockdown of particular factors can help to uncover SFs that mediate their effects. To evaluate the potential of this approach, cells were treated with drugs known to affect AS decisions and the patterns of changes in the 35 ASEs induced by these treatments were assessed by the same robotized procedure for RNA isolation and RT-PCR analysis used to locate their position within the network.
The results indicate that structurally similar drugs like Spliceostatin A (SSA) and Meayamycin cause AS changes that closely resemble the effects of knocking down components of U2 snRNP ( Figure 7A), including close links between these drugs and SF3B1, a known physical target of SSA (Kaida et al., 2007;Hasegawa et al., 2011) previously implicated in mediating the effects of the drug through alterations in the AS of cell cycle genes . Of interest, changes in AS of MCL1 (leading to the production of the proapoptotic mRNA isoform) appear as prominent effects of both drugs ( Figure 7B), confirming and extending recent observations obtained using Meayamycin (Gao et al., 2013) and suggesting that this ASE can play a key general role in mediating the antiproliferative effects of drugs that target core splicing components like SF3B1 (Bonnal et al., 2012). A strong link is also captured between each of the drugs and PPIH ( Figure 7A), a peptidyl-prolyl cis-trans isomerase which has been found to be associated with U4/5/6 tri-snRNP (Horowitz et al., 1997;Teigelkamp et al., 1998) but had not been implicated directly in 3' splice site regulation. Indeed, the network indicates a close relationship between PPIH and multiple U2 snRNP components, particularly in the SF3a and SF3b complexes, further arguing that PPIH plays a role in 3' splice site definition closely linked to branchpoint recognition by U2 snRNP. Treatment of the cells with the Clk (Cdc2-like) kinase family inhibitor TG003, which is also known to modulate AS (Muraki et al., 2004) and causes changes in ASEs analyzed in our network ( Figure 7B), led to links with a completely different set of SFs ( Figure 7A), attesting to the specificity of the results. These data confirm the potential of the network analysis to identify bona fide SF targets of physiological or pharmacological perturbations and further our understanding of the the underlying molecular mechanisms of regulation.

DISCUSSION
The combined experimental and computational approach described in this manuscript provides a rich source of information and a powerful toolset for studying the splicing process and its regulation. Similar approaches could be useful to systematically capture information about other aspects of transcrip-tional or posttranscriptional gene regulation. The data provide comprehensive information about SFs (and some additional chromatin factors) relevant to understand the regulation of 35 ASEs important for cell proliferation and apoptosis. Second, the network analysis reveals functional links between SFs, which in some cases are based upon physical interactions. The reconstruction of the known topology of some complexes (e.g., U2 snRNP) by the network suggests that the approach captures important aspects of the operation of these particles and therefore has the potential to provide novel insights into the composition and intricate workings of multiple Spliceosomal subcomplexes. Third, the network can serve as a resource for exploring mechanisms of AS regulation induced by perturbations of the system, including genetic manipulations, internal or external stimuli, or exposure to drugs. As an example, the network was used in the accompanying manuscript by Tejedor et al. (2014) to identify a link between AS changes induced by variations in intracellular iron and the function of SRSF7, which can be explained in terms of iron-induced changes in the RNA binding and splicing regulatory activity of this zinc-knuckle-containing SR protein. Other examples include (1) the identification of the extensive regulatory overlap between the interacting proteins IK/RED and SMU1, relevant for the control of programmed cell death; (2) evidence for a general role of uridine-rich hnRNP C binding sequences (some associated with transposable elements; Zarnack et al., 2013) in the regulation of nearby alternative exons via antagonism with U2AF; and (3) delineation of U2 snRNP components and other factors, including PPIH, that are functionally linked to the effects of antitumor drugs on AS.
Our data reveal changes of AS regulation mediated either by core splicing components or by classical regulatory factors like proteins of the SR and hnRNP families. While the central core of functional links is likely to operate globally, alternative configurations (particularly in the periphery of the network) are likely to emerge if similar approaches are applied to other cell types, genes, or biological contexts. We report considerable versatility in the effects of core spliceosome factors on alternative spicing. Precedents exist for a role of core SFs on the regulation of splicing efficiency or splice site selection in yeast, Drosophila, and mammalian cells (Clark et al., 2002;Park et al., 2004;Pleiss et al., 2007;Saltzman et al., 2011). Despite their assumed general roles in the splicing process, depletion or mutation of particular core factors led to differential splicing effects in these studies, and at least in some cases the differential effects could be attributed to particular features of the regulated pre-mRNAs. For example, Clark et al. (2002) found that the core components Prp17p and Prp18p are dispensable for splicing of yeast introns with short branchpoint to 3' splice site distances. Pleiss et al. (2007) found that introns of yeast ribosome protein genes are particularly sensitive to mutation of various core factors, including two DEAD/X box family of RNA helicases (PRP2 and PRP5) involved in spliceosome conformational transitions and PRP8, a highly conserved protein involved in catalytic activation 13 (Collins and Guthrie, 1999;Mozaffari-Jovin et al., 2012). Distinct effects of the mutations on different ribosomal genes were also observed, arguing again for specific effects on splicing of certain genes. DEAD/H box proteins as well as components of U1, U2, and U4/6 snRNP were found to modulate particular ASEs in an RNAi screen of RNA binding proteins in Drosophila cells (Park et al., 2004). In contrast with the standard function of classical splicing regulators acting through cognate binding sites specific to certain target RNAs, core factors could both carry out general, essential functions for intron removal and in addition display regulatory potential if their levels become limiting for the function of complexes, leading to differences in the efficiency or kinetics of assembly on alternative splice sites. This could be the case for the knockdown of SmB/B 0 , a component of the Sm complex present in most spliceosomal snRNPs, which leads to autoregulation of its own pre-mRNA as well as to effects on hundreds of other ASEs, particularly in genes encoding RNA processing factors (Saltzman et al., 2011).
The results of our genome-wide siRNA screen for regulators of Fas AS (Tejedor et al., 2014) provide unbiased evidence that an extensive number of core SFs have the potential to contribute to the regulation of splice site selection. This is, in fact, the most populated category of screen hits. Similarly, the results of our network highlight coherent effects on alternative splice site choice of multiple core components, revealing also that the extent of their effects on alternative splice site selection can be largely attributed to the duration and order of their recruitment in the splicing reaction.
Remarkably, effects on splice site selection are associated with depletion not only of complexes involved in early splice site recognition but also of factors involved in complex B formation or even in catalytic activation of fully assembled spliceosomes. While particular examples of AS regulation at the time of the transition from complex A to B and even at later steps had been reported (House and Lynch, 2006;Bonnal et al., 2008;Lallena et al., 2002), our results indicate that the implication of late factors is quite general. Such links, e.g., those involving PRP8 and other interacting factors closely positioned near the reacting chemical groups at the time of catalysis, like PRP31 or U5 200KD, are actually part of the persistent functional interactions that emerge from the analysis of any subsample of ASEs analyzed, highlighting their general regulatory potential in splice site selection. How can factors involved in late steps of the splicing process influence splice site choices? It is conceivable that limiting amounts of late spliceosome components would favor splicing of alternative splice site pairs harboring early complexes that can more efficiently recruit late factors, or influence the kinetics of conformational changes required for catalysis. In this context, the regulatory plasticity revealed by our network analysis may be contributed, at least in part, by the emerging realization that a substantial number of SFs contain disordered regions when analyzed in isolation (Korneta and Bujnicki, 2012;Chen and Moore, 2014). Such regions may be flexible to adopt different conformations in the presence of other spliceosomal components, allowing alternative routes for spliceosome assembly on different introns, with some pathways being more sensitive than others to the depletion of a general factor.
Kinetic effects on the assembly and/or engagement of splice sites to undergo catalysis would be particularly effective if spliceosome assembly is not an irreversible process. Results of single molecule analysis are indeed compatible with this concept (Tseng and Cheng, 2008;Abelson et al., 2010;Hoskins et al., 2011;Hoskins and Moore, 2012;Shcherbakova et al., 2013), thus opening the extraordinary complexity of conformational transitions and dynamic compositional changes of the spliceosome as possible targets for regulation. Interestingly, while depletion of early factors tends to favor exon skipping, depletion of late factors causes a similar number of exon inclusion and skipping effects, suggesting high plasticity in splice site choice at this stage of the process.
Given this potential, it is conceivable that modulation of the relative concentration of core components can function as a physiological mechanism for splicing regulation, for example during development and cell differentiation. Indeed, variations in the relative levels of core spliceosomal components have been reported (Wong et al., 2013). Furthermore, recent reports revealed the high incidence of mutations in core SFs in cancer. For example, the gene encoding SF3B1, a component of U2 snRNP involved in branchpoint recognition, has been described as among the most highly mutated in myelodysplastic syndromes (Yoshida et al., 2011), chronic lymphocytic leukemia (Quesada et al., 2012), and other cancers (reviewed in Bonnal et al., 2012), correlating with different disease outcomes. Remarkably, SF3B1 is the physical target of antitumor drugs like Spliceostatin A and-likely-Meayamycin, as captured by our study. These observations are again consistent with the idea that depletion, mutation, or drug-mediated inactivation of core SFs may not simply cause a collapse of the splicing process-at least under conditions of limited physical or functional depletion-but rather lead to changes in splice site selection that can contribute to physiological, pathological, or therapeutic outcomes. SF3B1 participates in the stabilization and proofreading of U2 snRNP binding to the branchpoint region , and it is therefore possible that variations in the activity of this protein result in switches between alternative 3' splice sites harboring branch sites with different strengths and/or flanking decoy binding sites. Similar concepts may apply to explain the effects of mutations in other SFs linked to disease, including, for example, mutations in PRP8 leading to retinitis pigmentosa (Pena et al., 2007).
In summary, the data and methodological approaches presented in this study provide a rich resource for understanding the function of the spliceosome and the mechanisms of AS regulation, including the identification of targets of physiological, pathological, or pharmacological perturbations within the complex splicing machinery.   Table S2). Endogenous mRNAs were purified 72 hr posttransfection by using oligo dT-coated 96-well plates (mRNA catcher PLUS, Life Technologies) following the manufacturer's instructions.
RNA samples corresponding to treatments with splicing arresting drugs were isolated with the Maxwell 16 LEV simplyRNA kit (Promega).
RT-PCR and High-Throughput Capillary Electrophoresis Cellular mRNAs were reverse transcribed using Superscript III Retro Transcriptase (Invitrogen, Life Technologies) following the manufacturer's recom-mendations. PCR reactions for every individual splicing event analyzed were carried out using forward and reverse primers in exonic sequences flanking the alternatively spliced region of interest and further reagents provided in the GOTaq DNA polymerase kit (GoTaq, Promega). Primers used in this study are listed in Table S6. HTCE measurements for the different splicing isoforms were performed in 96-well format in a Labchip GX Caliper workstation (Caliper, Perkin Elmer) us-ing a HT DNA High Sensitivity LabChip chip (Perkin Elmer). Data values were obtained using the Labchip GX software analysis tool (version 3.0).

Quantification of AS Changes from HTCE Measurements
Robust estimates of isoform ratios upon siRNA or pharmacological treatments were obtained using the median PSI indexes of the biological triplicates for each knockdown-AS pair. From this value xi the effect of each treatment was summarized as a robust Z score (Birmingham et al., 2009):

Robust Sample Correlation Estimation
Robust correlation estimation is based on an iterative weighting algorithm that discriminates between technical outliers and reliable measure-ments with high leverage. The weighting for mea-surements relies on calculation of a reliability index that takes into account their cumulative influence on correlation estimates of the complete data set. Algorithmic details are provided in Supplemental Experimental Procedures.

Network Reconstruction
The process of network reconstruction relies on the regularization-based graphical lasso (glasso) algorithm for graphical model selection . The glasso process for network reconstruction was implemented using the glasso R package by Jerome Friedman, Trevor Hastie, and Rob Tibshirani (http://statweb.stanford.edu/~tibs/glasso/). A detailed exposition of the algorithmic steps for network reconstruction and analysis, glasso, and its use for graphical model selection is provided in Supplemental Experimental Procedures.
Principal Component Analysis of Ancillary Connections PCA for ancillary edges was performed using R-mode PCA (R function princomp). The input data set is the scaled 95 3 35 matrix containing the absolute average correlation for each of the 95 edges over the subsets of the 10,000 samples that included each of the 35 events. The coordinates for projecting the 95 edges are directly derived from their scores on the first two principal components. The arrows representing the top ten discriminatory events (based on the vector norm of the first two PC loadings) were drawn from the origin to the coordinate defined by the first two PC loadings scaled by a con-stant factor for display purposes.

Splice Site Scoring, Identification of Probable HNRPC Binding Sites and Alu Elements
For the identification and scoring of 3' ss we used custom-built position weight matrices (PWMs) of length 21-8 intronic plus three exonic positions. The matrices were built using a set of human splice sites from constitutive, internal exons compiled from the hg19 UCSC annotation database (Karolchik et al., 2014). Background nucleotide frequencies were estimated from a set of strictly intronic regions. Threshold was set based on a FDR of 0.5/kbp of random in-tronic sequences generated using a second-order Markov chain of actual in-tronic regions.
We considered as probable HNRPC binding sites consecutive stretches of U{4} as reported in Zarnack et al. (2013). Annotation of ALU elements was based on the RepeatMasker track developed by Arian Smit (http://www. repeatmasker.org) and the Repbase library (Jurka et al., 2005).
More detailed experimental procedures are provided in Supplemental Experimental Procedures.

ACKNOWLEDGMENTS
We thank many CRG colleagues, EURASNET members, and Quaid Morris for their advice, encouragement, and critical reading of the manuscript, and Drs. Koide and Yoshida for reagents. We acknowledge the excellent technical support of the CRG Robotics and Genomics facilities. J.R.T. was supported by a PhD fellowship from Fondo de Investigaciones Sanitarias. Work in our lab is supported by Fundació n Botín, Consolider RNAREG, Ministerio de Economía e Innovació n, and AGAUR.

ATTACHED BELOW
RT--PCR assays (lower two panels) were used to detect changes in AS of Fas/CD95 and MCL1 endogenous transcripts.   Table S1. Information on AS events analyzed in this study. The information includes gene symbol, id, exon sequence and genomic coordinates for the AS exons analyzed. Table S2. Information about the custom siRNA library utilized to knockdown splicing and chromatin factors. Information includes siRNA pools catalog numbers, gene symbol, id and accession numbers for the targeted genes. It also includes information about the categories in which their protein products were classified according to two publications (Wahl et al., 2009 andZhou et al., 2002), the classifications used in our network analysis and also information about their functions and known interactions.   The files include a schematic representation for every AS event; information regarding the functional role of the gene, as retrieved from uniprot / genecards database; the protein domain affected by the AS event; the functional effect of each of the isoforms towards cell proliferation or apoptosis; and 102 references consulted for the selection.

Cell lines
HeLa CCL--2 cells and HEK293 cells were purchased from the American Type Culture Collection (ATCC). Cells were cultured in Glutamax Dulbecco's modified Eagle's medium (Gibco, Life Technologies) supplemented with 10% fetal bovine serum (Gibco, Life Technologies) and penicillin/streptomycin antibiotics (penicillin 500 u/ml; streptomycin 0.5 mg/ml, Life Technologies). Cell culture was performed in cell culture dishes in a humidified incubator at 37 o C under 5% CO2.

Drug Treatments
Pharmacological treatments with splicing arresting drugs were performed as described in . HeLa cells were treated for 3 hours with Spliceostatin A 260 nM or 8 hours with Meayamycin 20 nM or TG003 10 µM. RNA was isolated using the Maxwell 16 LEV simplyRNA kit (Promega) and RT PCRs for the 35 splicing events were performed as described above.

Cell proliferation assays
HeLa cells were forward--transfected with siRNAs targeting IK, SMU1 or both genes by plating 2500 cells over a mixture containing siRNA--RNAiMAX lipofectamine complexes, as previously described. Indirect estimations of cell proliferation were

Cell apoptosis assays
Hela cells were transfected with siRNAs against IK, SMU1 or both genes, as described above. Cells were trypsinized and collected for analysis 72 hours post transfection.
Pellets were washed once with complete DMEM medium (10% serum plus antibiotics), and three times with PBS 1X. Cells were then resuspended in 200 µl Binding buffer provided by the Annexin V--FITC Apoptosis Detection Kit (eBiosciences), resulting in a final density of 4x10 5 cells/ml. 5 µl of Annexin V--FITC was added to each sample followed by 10 min incubation at room temperature in dark conditions. Cells were washed once with 500 µl binding buffer and incubated with 10 µl propidium iodide (20 µg/ml) for 5 minutes. FACS analysis was performed using a FACSCalibur flow cytometry system (BD Biosciences).

Removal of sparse variables and imputation of missing values
Estimation of the sample covariance matrix requires complete data. However the p x n matrix that contains the measurements of the p variables (genes) in n conditions (ASEs) includes missing (NA) values. In order to overcome this problem, first sparse variables (#NA values > n/2) are removed. The remaining missing values are subsequently filled--in by applying k--nearest--neighbor based imputation (Hastie et al., 1999, R function impute available as part of the bioconductor project at http://www.bioconductor.org/). For the purposes of this work we set k=ceiling(0.25 sqrt(p)).

Removal of uninformative ASEs.
ASEs not affected to a significant degree (median absolute (ΔPSI)<1) by the KDs offer little information on the genes' covariance and can potentially introduce noise and are thus discarded from subsequent analyses.

Data Scaling.
Prior to sample covariance estimation the data matrix is standardized along the gene dimensions. Standardization of the events ensures that only the shape and not the magnitude of the fluctuations is taken into account for covariance estimation and it is equivalent to using correlation in lieu of covariance in all downstream analyses.

Robust sample correlation estimation.
When outliers are present they dominate the correlation estimates. In order to overcome the sensitivity of sample correlation to outliers we derive robust correlation estimates based on a weighting algorithm that attempts to discriminate between technical outliers and reliable influential measurements.
In particular, starting with an M = p x n dataset of p standardized variables and n samples we wish to derive a robust correlation p x p matrix. We calculate the influence of a sample u in the estimation of the correlation of two variables X a , X b using the deleted residual distance: Where p(X a , X b ) and p --u (X a ,X b ) are the Pearson's correlation estimates before or after removing the observations M au , M bu coming from sample u.
A measure of the reliability !,! ! of these measurements for the estimation of p(X a , X b ) given all the observed data is then based on their cumulative influence on high correlates of X a , and X b : where , is the indicator function: T being a minimum correlation threshold for considering a pair of variables.
The final robust correlation estimates are finally calculated as the weighted Pearson's correlation p(X a , X b ;w) where w is a vector of weights of length n containing the values !,! ! for every observation u.
Intuitively, if an observation M au distorts correlation estimates among all highly correlated partners of X a it is considered an outlier and is down--weighted. On the other hand if including this observation yields similar correlation estimates for most of these pairs it is considered reliable even if its influence on the individual estimate of p(X a , X b ) is high. The advantage of this procedure is that it derives weights for individual measurements (as opposed to a single weight for every sample) while taking into account the complete dataset.
The algorithm can be used iteratively, however in our experience the estimates converge within 1 to 2 iterations. R code for implementing this algorithm is available upon request. Examples of robust correlation estimates and regression based on this metric and taken from our dataset are shown below: Graphical lasso (gLasso) provides an attractive solution to the problem of covariance estimation for undirected models, when graph sparsity is a goal . The gLasso algorithm has the advantage of being consistent, i.e. in the presence of infinite samples its parameter estimates will be arbitrarily close to the true estimates with probability 1.

Graphical model selection for Network Reconstruction using Graphical Lasso
L 1 regularization (lasso) is a smooth form of subset selection for achieving sparsity. In the case of gLasso the model is constructed by optimizing the log-likelihood function: Where Θ is an estimate for the inverse covariance matrix Σ --1 , S is the empirical covariance matrix of the data, ||Θ|| 1 is the L1 norm i.e the sum of the absolute values of all elements in Θ --1 , and r is a regularization parameter (in this case selected based on estimates of the FDR, see below). The solution to the glasso optimization is convex and can be obtained using coordinate descent.

Partition of Network into modules.
Network modules or communities can be defined loosely as sets of nodes with a more dense connection pattern among their members than between their members and the remainder of the. Module detection in real--world graphs is a problem of considerable practical interest as they often capture meaningful functional groupings.
Typically community structure detection methods try to maximize modularity, a measure of the overall quality of a certain network partition in terms of the identified modules (Newman, 2006). In particular, for a network partition p, the network modularity M(p) is defined as: where m is the number of modules in p, l k is the number of connections within module k, L is the total number of network connections and d k is the sum of the degrees of the nodes in module k.
Here, we identify modules of genes that exhibit similar perturbation profiles among the assayed events by maximizing the network's modularity using the greedy community detection algorithm (Clauset et al., 2004) implemented in the