Using TCGA data for 11 solid tumor types (Supplemental Table S1), we analyzed the differential gene expression between normal and tumor sample pairs of 1348 genes encoding known and predicted RNA binding proteins (RBPs) (Methods; Supplemental Table S2). The majority of these genes (1143, 84.8%) show significant differential expression in at least one tumor type (Supplemental Fig. S1; Supplemental Table S3). Examining in detail 162 RBP genes annotated as known or putative splicing factors (SFs), we find they can be separated into three groups. One group is frequently up-regulated, another one down-regulated, and a third one shows opposite patterns in the three kidney tumor types (KICH, KIRC, KIRP) compared to other tumor types (Fig. 1A). Moreover, 132 (80%) of them are differentially expressed in at least one tumor type and 45 were previously associated with oncogenic or tumor suppressor activities (Fig. 1A, labeled in red; Supplemental Table S4). We also found apparent discrepancies with previous literature. For instance, although SRSF5 was described as oncogenic (Huang et al. 2007) , it is down-regulated in six tumor types; and TRA2B, reported as oncogenic in breast cancer (Watermann et al. 2006) , is up-regulated in LUSC but down-regulated in KICH and thyroid carcinoma (THCA). Also, the oncogenic SRSF2, SRSF3, and SRSF6 (Xiao et al. 2007; Jia et al. 2010; Jensen et al. 2014) are down-regulated in KICH, while the oncogenic SRSF1 and the tumor suppressor RBM4 do not show any significant expression changes. New patterns also emerge, including up-regulation of genes from the RBM family and down-regulation of the genes from the MBNL family (Fig. 1A). Unsupervised clustering of all the 4442 tumor samples using normalized expression values per sample for SFs (Supplemental Fig. S2) or for all RBPs (Supplemental Fig. S3) largely separates samples by tumor type. This pattern is reproduced when using a different gene set of similar size (Supplemental Fig. S3). Closer inspection of the BRCA and COAD samples reveals that the expression patterns of RBPs also reproduce tumor subtypes (Supplemental Figs. 4, 5). Collectively, expression analyses indicate that RBPs are frequently and specifically altered in human tumors, and they represent one of the multiple expression profiles that characterize the intrinsic properties of the tumor types and their tissue of origin.

To further define the extent to which RBP genes are altered in human cancer, the TCGA data were analyzed for protein-affecting mutations and copy number variations (CNVs). Although most of the RBP genes are mutated in tumors (Supplemental Table S4), they are generally mutated in fewer samples compared to candidate cancer drivers (listed in Supplemental Tables S5, S6) and other genes (Fig. 1B). We confirmed 5.4% (25/458) of LUAD tumor samples to have protein-affecting mutations for RBM10, in agreement with previous studies (Imielinski et al. 2012) . Using this case as a reference, we observed 205 (15.2%) of all RBPs (13 or 8% of the 162 SFs) mutated in >5% of samples in a given tumor type (Table 1; Supplemental Table S3). In general, there is a weak association between mutations and expression changes, whereas a number of genes show mutual exclusion of mutations and expression Downregulated

Opposing

Upregulated C BRCA COAD HNSC KICH KIRC KIRP LIHC LUAD LUSC PRAD THCA CNV Gain CNV Loss Gene type

RBP

Driver

Other 1PABPC 3BSHRD 1FXR 1SEPR 2TABR I22FBPG 1LBNM L3FEC 93BRM 38BRM LAYR 6FSSR PRNUN 2TAAR 21PABRN I23FBPG 2LBNM 4TNCO 3APBPC 2L1PANR 1EXNO 82BRM 2SKPR 3SRRM FSU 14L2FAU 2SRRM 24BRM 1FXBRO 45SPR 24BRM 1FSSR 3FSSR 1PKSR LPNNRH I12FBPG APSRN 2FSSR 71BRM 3FXBRO 4LFEC 3APNNR LL2EAV L12PANR 314ZCH 2FXR 2YBX 7BRM 1SERW 3PBPC 13FSA 2PBEC 74BRM 12FAU K H NH NH H NH changes (Table 1; Supplemental Table S3). On the other hand, CNVs are highly recurrent across samples, with gains more frequent than losses (Fig. 1C; complete list of CNVs available in Supplemental Table S3). Up-regulated RBPs had generally more CNV gains than nonregulated RBPs (Mann-Whitney U test P-value < 2.2 × 10−16), and 90% of the up-regulated RBPs had CNV gains in >50% of the samples from the same tumor type. The splicing factors ESRP1, PABPC1, and KHDRBS3, which are near the reported 8q24 amplification (The Cancer Genome Atlas Network 2012b) , show CNVs in BRCA, COAD, HNSC, LIHC, and LUAD (Fig. 1C), with ESRP1 and PABPC1 showing frequent association with upregulation (Table 1). We also observed the frequent amplification of TRA2B, IGF2BP2, and FXR1 in HNSC and LUSC, as reported before (The Cancer Genome Atlas Network 2012a, 2015) . However, the association with up-regulation is only observed in LUSC (Table 1). CNV losses are less frequent than gains and show weaker associations with down-regulation (Table 1; Supplemental. Table S3). Finally, only a low fraction of the CNVs were detected as focal (Table 1; Supplemental Table S3). These results indicate that many of the detected expression alterations of RBPs may be explained by CNVs, although most often in association with large chromosomal events, thereby highlighting the potential relevance of these alterations in shaping the tumor phenotypes.

Next, we investigated the patterns of differential splicing that may occur as a consequence of the alterations described above. To determine those possibly related to expression alterations in RBPs, we first evaluated the significant splicing changes in tumors compared to normal tissue, considering five major event types: skipping exon (SE), alternative 5′ splice-site (A5), alternative 3′ splice-site (A3), mutually exclusive exon (MX), and retained intron (RI) events (Fig. 2A; Supplemental Tables S7, S8). We examined the splicing patterns of cancer drivers in more detail. From the 937 drivers collected (Supplemental Tables S5, S6), 653 (69.7%) have annotated events and 292 (31.2%) have at least one differentially spliced event, with a number of them occurring in multiple tumor types (Supplemental Fig. S6). Moreover, seven of the 11 tumor types show enrichment of differentially spliced events in drivers (Fig. 2B, in red; Supplemental Table S9). Additionally, various cancer hallmarks (Liberzon et al. 2015) are enriched in differentially spliced events but not in differentially expressed genes (Fig. 2C; Supplemental Fig. S7; Supplemental Table S10). These results suggest that alternative splicing contributes to cancer development independently of mutations or expression alterations.

To determine the splicing changes related to mutations in RBPs, we compared the inclusion levels (percent spliced in, PSI) of the events between samples with or without protein-affecting mutations for each RBP (Fig. 2D,E; Supplemental Fig. S8; Supplemental Table S11). We found less differentially spliced events compared with the tumor-normal comparison, and some of the detected events affect cancer drivers (Table 2). HNRNPL had 16 mutations in COAD, 13 of them indels in an RNA recognition motif, causing frameshifts (Fig. 2F). These mutations had 42 events associated with them (Z-score = 41.35), including the cancer driver CASP8 (Fig. 2G). In COAD, we also found the putative RBP MACF1 (Baltz et al. 2012) with mutations along the entire protein (42 events, Z-score = 45.35) (Supplemental Fig. S9). MACF1 is a component of the WNT-pathway known to affect splicing (Bordonaro 2013) . In LUAD, we found 53 events associated with RBM10 (Z-score = 52.35), 83 events associated with U2AF1 (Z-score = 81.35), and 49 events associated with the transcription factor and putative RBP TCF20 (Castello et al. 2012) (Z-score = 50.35). Most of the TCF20 alterations consist of an insertion in a glycine-rich region at the N terminus (Supplemental Fig. S9). As a comparison, we analyzed SF3B1 and found 64 differentially spliced events (Z-score = 63.35) in BRCA, despite being mutated in only 1.7% of the tested samples (16 out of 956 for which we had mutation and RNA-seq data) (Supplemental Fig. S9). Notably, compared with the proportions of the different event types, there is a significant enrichment of A3 events associated with SF3B1 (Fisher’s exact test P-value = 1.45 × 10−11), of A5 events for TCF20 (P-value = 2.09 × 10−6), and of SE events for RBM10 (Pvalue = 0.003) (Supplemental Fig. S9; Supplemental Table S12). There is also a significant depletion of SE events for TCF20 (P-value = 0.0014), of SE events for SF3B1 (P-value = 2.24 × 10−8), and of A5 events for U2AF1 (P-value = 1.67 × 10−3) and HNRNPL (P-value = 0.005) (Fig. 2H). Furthermore, 17 (20%) of the genes detected to change splicing with SF3B1, 32 (38%) associated with U2AF1, and 21 (30%) related with RBM10 mutations coincide with those obtained in recent analyses (Supplemental Table S13; Brooks et al. 2014; Darman et al. 2015; Maguire et al. 2015) . The limited overlap could be due to the use of different methods for differential splicing analysis and the use of different sample sets. In summary, we validated known RBP mutations affecting alternative splicing and described new ones, whose mutations potentially affect splicing patterns. Our analyses reveal a rich source of new information about alternative splicing events associated with RBP alterations with a potential relevance in cancer.

The above results suggest that mutations in RBPs are not the single cause of splicing changes in tumors. Therefore, we further characterized the splicing changes between tumor and normal samples to determine their association with the expression of RBPs. We first identified common patterns of splicing changes between pairs of tumor types by selecting events with a strong correlation with a differentially expressed SF in the two tumor types (Methods). PSI changes (ΔPSI) for these events show high correlation between tumor pairs (Fig. 3A) and indicate potential common regulators (Fig. 3B; Supplemental Fig. S10; Supplemental Table S14). BRCA and LUAD share 229 events associated with various factors, including QKI and SRSF5 (Fig. 3B, upper left panel), whereas KIRC and PRAD share 78 anti-correlating events associated with ESRP2 and MBNL1 (Fig. 3B, upper right panel). RBM47, described as a tumor suppressor (Vanharanta et al. 2014) , appears as the main common SF between KIRC and HNSC, with 61 associated events (Fig. 3B, lower left panel). HNSC and LUSC share 141 events, 63 of them associated with RBM28 (Fig. 3B, lower right panel). These results suggest functional connections between RBPs and the splicing changes detected in cancer. Next, we also studied whether there are tumor-specific events (Methods) and found 380 events that largely separate the 4442 tumor samples by type (Fig. 3C; Supplemental Fig. S11; Supplemental Table S15). These splicing changes may be indicative of tumor-type–specific oncogenic mechanisms.

To further understand the link between the observed splicing patterns and the RBPs, we tested the enrichment of binding motifs in differentially spliced events in each tumor type. We assigned

Jac | Mex 0.02 | 0.10 0.05 | 0.03 0.07 | 0.02 0.09 | 0.07 0.00 | 0.10 0.10 | 0.16 0.14 | 0.10 0.06 | 0.13 0.05 | 0.07 0.11 | 0.06 0.04 | 0.05 0.08 | 0.08 0.11 | 0.02 0.04 | 0.07 0.05 | 0.08 3.66 0.73 3.77 5.24 2.38 1.43 0.95 0.61 1.58 2.11 0.53 1.05 6.11 2.40 1.09 0.22 7.30 3.37 3.37 2.81 4.49 0.81 0.10 0.10 0.53 0.44 0.56 For each RBP and each tumor type, we indicate the log2-fold change (logFC) and adjusted P-value (Pval) of the differential expression analysis between tumor and normal samples, the frequency of the alteration (Frequency), and the association of the alteration with the expression alteration (Assoc. Freq.). For mutations, we show those cases with mutation frequency >5% and with a Jaccard (Jac) or mutual exclusion (Mex) score >0.05. For CNV gains, we show those cases that have significant up-regulation, CNV gain, and association frequencies in >5% of the samples, with one or more of the CNVs being focal (Focal. Freq). For CNV losses, we show those with expression frequency and association with down-regulation in >0.5% of the samples, with one or more of the CNVs being focal (Focal Freq). Complete data for all RBPs are given in the Supplemental Table S3. binding motifs from RNAcompete (Ray et al. 2013) to 104 of the analyzed RBPs and tested their enrichment (see Methods; Supplemental Figs. S12, S13). Considering enriched motifs of differentially expressed RBPs in the same tumor type, we observed the CELF, RBFOX, and MBNL motifs among the most frequently enriched across tumor types, as well as in luminal breast tumors (Fig. 4A; Supplemental Figs. S14–S18). Down-regulated RBP genes in inclusion events and up-regulated ones in skipping events show more frequently enriched motifs in upstream and exonic regions, consistent with a positional effect (Fig. 4B). On the other hand, down-regulated RBP genes in skipping events and up-regulated ones in inclusion events show enriched motifs most frequently on exons, suggesting that RBPs more often enhance the inclusion of bound exons.

To define candidate target events for differentially expressed RBPs, we selected differentially spliced events whose PSI correlates with the RBP expression (|R| > 0.5, Spearman) and that contain the corresponding RNA binding motif. We could assign between 20% and 80% of the differentially spliced events to at least one RBP (Fig. 4C; Supplemental Tables S16, S17). We ranked the RBPs in each expression group according to the total number of cancer drivers with differentially spliced events across all tumor types (Fig. 4D; Supplemental Fig. S19). The results of this analysis emphasize the relevance of CELF2, ESRP1, ELAVL1, RBFOX2, PTBP1, QKI,

TRA2B, and RBM47 in relation to alterations of splicing in cancer and highlight novel prominent factors, including HNRNPC, SYNCRIP (also known as HNRNPQ), and genes from the CPEB and MBNL families. We confirmed the role of MBNL1, QKI, RBFOX2, PTBP1, RBM47, and ESRP1 in some of these tumors by comparing the ΔPSI values of the events with those obtained by single knockdown or overexpression experiments in cell lines (Fig. 4E, upper panels; Supplemental Figs. S20–S22; Supplemental Table S18).

Additionally, comparing the tumor events with those differentially spliced between human embryonic stem cells (hESC) and differentiated cells or tissues (Han et al. 2013) yielded a high positive correlation of ΔPSI values in BRCA, PRAD, LIHC, and breast luminal tumors, to a lesser extent in COAD and LUAD, and an anti-correlation in KIRC. This is in agreement with the MBNL1 and MBNL2 expression patterns in these tumors and with the majority of the correlating events containing the MBNL motif (Fig. 4E, lower panels; Supplemental Fig. S23; Supplemental Table S19). These results highlight the prominent role of RBP expression changes in the alterations of splicing in cancer and suggest new mechanisms of regulation.

To identify splicing regulatory networks relevant in cancer, we built clusters with the 162 SFs using the correlation between trolled by RBM47, PTBP1, and RBM28 (Fig. 5E). This analysis reveals new roles of RBPs and splicing in cancer-relevant processes. For each tumor type and for each RBP with more than 10 associated differentially spliced events, we indicate the events in cancer drivers predicted to have a significant splicing change. The table includes the PSI change between mutated and nonmutated samples (ΔPSI) and the Pvalue of the comparison after correcting for multiple testing (Pval). SF3B1 is included for comparison. Coordinates for the events are given in Supplemental Table S7. gene expression and event PSI values. We linked these clusters to differentially spliced genes in enriched cancer hallmarks (Methods). This analysis revealed modules of splicing regulation, with one or two genes as the main regulator of each hallmark across different tumor types. Myogenesis and EMT networks are enriched in almost all tumor types and controlled by similar factors (Fig. 5A,B; Supplemental Figs. S24, S25). Other relevant regulatory modules are the G2 checkpoint (G2M), which includes NUMA1, a gene involved in spindle formation (Zheng et al. 2010) , and the WNT/beta-catenin pathway, which includes NUMB, whose alternative splicing is linked to cell proliferation (Bechara et al. 2013) . We predicted the splicing of these genes to be controlled by MBNL1, among other factors (Fig. 5C,D). Interestingly, angiogenesis, which is an enriched hallmark in COAD for splicing but not for gene expression, includes an event in SERPINA5, an inhibitor of serine proteases involved in homeostasis and thrombosis (Suzuki 2008) , which we predict to be con

MBNL1 emerges as a potentially key regulator of splicing for multiple cancer drivers, especially in luminal breast tumors (Supplemental Fig. S26). In particular, MBNL1 is predicted to control an exon skipping event in NUMA1, whose PSI value correlates with MBNL1 expression (Spearman R = 0.65/0.66 in luminal A/B) and contains the MBNL motif (Fig. 6A; Supplemental Fig. S27). The same event has increased PSI value in KIRC (ΔPSI = 0.11, corrected P-value = 3.11 × 10−6), where MBNL1 is significantly up-regulated compared to normal tissues, providing further support for the dependence of NUMA1 splicing on MBNL1. We detected MBNL1 protein in the breast epithelial cell line MCF10A, but not in the luminal-like MCF7 (Supplemental Fig. S28). MBNL1 expression depletion in MCF10A cells with two different siRNAs induced skipping of exon 16 in NUMA1 (Supplemental Fig. S28), recapitulating the pattern observed in the tumor samples (Fig. 6B, upper panel) and in MCF7 (Supplemental Fig. S28). This result was validated in the same samples by real-time quantitative RT-PCR (Supplemental Fig. S29). As a control for NUMA1 gene expression level, we measured the expression of constitutive exon 5 and 7 using semiquantitative (Supplemental Fig. S28) and real-time quantitative (Supplemental Fig. S29) RT-PCR and found no change upon MBNL1 depletion. This is in agreement with the lack of differential expression observed for NUMA1 in all breast tumor samples (corrected P-value = 0.68) and in luminal tumors (corrected P-value = 0.50), compared to the matched normal samples. We also tested the alternative splicing of NUMB exon 9, which we predicted to be dependent on MBNL1 in BRCA luminal tumors. The depletion of MBNL1 recapitulates the NUMB splicing pattern in luminal samples (Fig. 6B, middle panel).

For a comparison with MBNL1, we evaluated the role of QKI, which we also observed down-regulated in BRCA luminal tumors and whose product was detected in MCF10A but not in MCF7 cells (Supplemental Fig. S28). Upon QKI expression depletion, NUMA1 exon 16 inclusion changed by a small but reproducible amount in the direction opposite to that with MBNL1 depletion (Fig. 6B, upper panel). Although we did not find a QKI motif on the NUMA1 event, this is consistent with the negative correlation found with QKI expression (R = −0.11) in BRCA. Next, we also tested NUMB exon 9, which we predicted to be controlled by QKI in BRCA luminal tumors. QKI depletion induces exon 9 inclusion, recapitulating the pattern in luminal samples (Fig. 6B, middle panels). NUMA1 alternative splicing is likely controlled by other factors (Fig. 5C). For instance, NUMA1 exon 16 contains an RBM42 motif and its PSI anti-correlates with RBM42 expression in breast luminal tumors; hence, RBM42 potentially acts as a repressor of exon inclusion. Consistent with this, depletion of RBM42 in MCF7 cells leads to increased inclusion of NUMA1 exon 16 (Supplemental Fig. S30).

To measure whether the splicing change in NUMA1 had any effect on cancer cell hallmarks, we used antisense oligonucleotides (AONs) targeting specifically the 5′ and 3′ splice-sites of NUMA1 exon 16. These AONs promote exon skipping, recapitulating in the MCF10A cells the splicing pattern observed in BRCA luminal tumors, with the AON against the 5′ splice-site being more efficient (Supplemental Fig. S31), as validated using real-time quantitative RT-PCR (Supplemental Fig. S32). Importantly, the AONs did not affect NUMA1 expression as measured by semiquantitative

0.7 Supplemental Table S20), and when transfecting cells with the AON against the 5′ splice site (Student’s t-test P-values < 0.05) (Fig. 6C). Using only the 3′ splice site or both AONs, there was also an increase, albeit not statistically significant, possibly due to a smaller effect of the 3′AON on NUMA1 splicing and the shared concentration of both AONs. Furthermore, although QKI has a mild effect on NUMA1 exon 16 splicing, QKI down-regulation U p r e g .

O p p .

D o w n r e g .

U p r e g .

O p p .

D o w n r e g .

B 50 40 ed 30 ich fs r i n t foe om20 % 10 0 1.0 0.5 L C − SC0.0 E I S P Δ−0.5 −1.0 1.0 a0.5 L e H −D0.0 K I S P Δ−0.5 −1.0 U

E

D

U

E

D Expression Downreg. Upreg.

R = 0.72

R = 0.84 −1.0 −0.5 0.0 0.5 ΔPSI Tumor − Normal

BRCA 1.0 −1.0

1.0 −0.5 0.0 0.5 ΔPSI Tumor − Normal

PRAD R = 0.6

R = 0.68 BARC ACDO SHCN IKHC IKRC IPKR ILCH LADU LSCU APDR TACH −1.0 −0.5 0.0 0.5 ΔPSI Tumor − Normal 1.0 −1.0 −0.5 0.0 0.5 ΔPSI Tumor − Normal 1.0 Downregulated has a much stronger effect on NUMB alternative splicing (Fig. 6B), which strongly promotes cell proliferation in a variety of cell types (Misquitta-Ali et al. 2011; Bechara et al. 2013) .

We further studied the possible effects of the alternative splicing of NUMA1 exon 16 on centrosome amplification, a hallmark of breast carcinogenesis. Using the AON against the 5′ splice-site, there is a significant increase in the number of cells with centrosome amplification compared with controls in MCF10A cells (Fig. 6D; Supplemental Table S21). In contrast, the siRNA against MBNL1 did not yield a significant difference, possibly due to the superposition of indirect effects. We also observed an inverse correlation between a signature for chromosome instability (Carter et al. 2006) and the NUMA1 event PSI in luminal tumors (Fig. 6E), which was absent in other tumor types (Supplemental Fig. S33), providing further support for a relationship between NUMA1 alternative splicing and the fidelity of centrosome formation. The exon skipping described in NUMA1 is the only coding difference between the tumor and the normal isoforms. While it is unclear whether this 14-amino acid change alone could explain the observed effects, e.g., we could not detect any protein domain or disordered region ( Dosztányi et al. 2005 ), using GPS (Xue et al. 2011) we predicted loss of a high scoring threonine phosphorylation site

A

Myogenesis

EMT 1BEPC 4LFEC 4LLVAE L3LVEA L1BNM 1PBAPC 3PPABC 1AVNO 3FXBRO 4AASDM 2TPBP IKQ 5PABPC 02BRM 3BSRM 24BRM L2FEC 3PEBC 2BSKRDH 74BRM 1FAC 2BPEC PCRNNH 7FSSR 1TPBP FSPQ 21PABNRNH IPSYNCR FSU 22FAU 31BPG 82BRM SPANR 1APRNNH 2L1PANRNH 1PPCR 07PSRNN 1FBXRO 6LFEC 4BEPC I32FPBG I1TA E Our study reveals that alterations in genes encoding RNA binding proteins are pervasive in cancer, characterize the different tumor types, and can explain many of the alternative splicing changes observed in tumors. Many of the RBP expression changes appear related to large copy number alterations. In contrast, mutations in RBPs are particularly low in frequency compared to mutation rates in other genes. We measured the potential relevance of RBP alterations by looking at the associated splicing changes. Only a few RBPs beyond the previously described cases show mutations associated with genome-wide effects on alternative splicing. One of the novel cases is HNRNPL, for which we predict that frequent indels in an RNA binding domain would affect the splicing of CASP8, a gene involved in programmed cell death (Yang et al. 2008) . As HNRNPL also controls CASP9 alternative splicing (Goehe et al. 2010) , this suggests a relevant role of HNRNPL in apoptosis. On the other hand, the expression changes in RBP genes appear as major contributors of the alternative splicing changes observed in tumors, since the number of events affected by RBP mutations is lower than those changing between tumor and normal samples. The splicing changes detected are predicted to have an impact on many cancer hallmarks, some of them different from those linked to gene expression changes. This agrees with previous reports on the impact of splicing on cancer hallmarks (Oltean and Bates 2013) and highlights the relevance of alternative splicing as a complementary molecular mechanism to explain tumor development. An important implication of our study for prognostic and clinical studies is that the definition of the functional impact of somatic genetic and epigenetic alterations should be expanded to include changes in the alternative splicing of the genes in cancer-relevant pathways.

We identified many RBPs whose expression alteration has some potential relevance in cancer. For instance, TRA2B is frequently amplified and up-regulated in LUSC and its binding motif enriched in differentially spliced events, including an SE event in the DNA damage-response gene CHEK1, reported recently to be controlled by TRA2 proteins (Best et al. 2014) . Despite the similar genetic alterations between squamous tumors (Hoadley et al. 2014) , only LUSC shows overexpression of TRA2B, which could explain some of the differences in splicing patterns between these two tumor types. We also found that genes from the MBNL family are frequently down-regulated in tumors and their binding motif is enriched in differentially spliced events, especially in breast and prostate tumors. Additionally, splicing changes related to MBNL1 recapitulate the splicing patterns of undifferentiated cells, in agreement with recent studies describing MBNL1 and MBNL2 as regulators of a stem cell-related splicing program (Han et al. 2013; Venables et al. 2013) . MBNL1 potentially controls multiple genes that participate in cancer-related pathways, including the mitotic gene NUMA1. We described that NUMA1 alternative splicing leads to higher proliferation and increased centrosome amplification in normal cells. Notably, the alternative splicing change is not associated with an expression change in NUMA1, as measured from the tumor and normal samples and by semiquantitative and quantitative real-time RT-PCR in cells after MBNL1 depletion or modification of the NUMA1 splicing by anti-sense oligonucleotides. Further, NUMA1 does not change expression between MCF7 and MCF10 cells (Allo et al. 2014) . Although the NUMA1 locus was related before to breast cancer risk (Kammerer et al. 2005) , a clear mechanism explaining its relevance in cancer is still lacking. NUMA1 produces a protein component of the nuclear matrix, which is dependent on threonine phosphorylation to regulate the orientation of mitotic spindles and ensure symmetric cell division (Kotak et al. 2013) . NUMA1 alternative splicing removes a potential threonine phosphorylation site; hence, one attractive hypothesis is that this splicing change affects NUMA1 phosphorylation, thereby impairing correct spindle positioning, leading to increased genome instability. Besides MBNL1, our analyses provided evidence for the control of NUMA1 splicing by RBM42, in agreement with its proposed role in cell cycle (Suvorova et al. 2013) and pointing to an involvement of this factor in alternative splicing and cancer worth investigating further.

The origin of the expression changes in RBPs remains to be described. In particular, those cases that cannot be explained by DNA alterations in the gene locus may originate from alterations in the pathways that control the regulation of RBPs. It was recently postulated that RBPs involved in development and differentiation controlled by common enhancers could act as master regulators (Jangi and Sharp 2014) . These included MBNL1, RBFOX2, RBM24, RBM38, RBM20, RBFOX1, ZNF638, and RBMS3, which we found frequently down-regulated in tumors, with enriched motifs in differentially spliced events and potential targets in multiple cancer drivers. This suggests a general mechanism in tumors by which the de-activation of one or more RBPs through the alteration of a common enhancer would lead to the reversal of multiple RNA splicing patterns to trigger or sustain an undifferentiated phenotype. Taken together, our results provide a rich resource of information about novel networks of RBPs that trigger common and specific alternative splicing changes in several solid tumors that may be relevant to understand the molecular basis of, and potentially reverse, the oncogenic properties of tumor cells.

Data sets were downloaded from the TCGA data portal (https:// tcga-data.nci.nih.gov/tcga/) (Supplemental Table S1). Details are provided as Supplemental Material. The 1348 genes coding for RNA binding proteins analyzed include those with high confidence for RNA binding (Baltz et al. 2012; Castello et al. 2012; Kwon et al. 2013) and those annotated as RNA-binding in Ensembl (Cunningham et al. 2014; Supplemental Table S2) . From this set, a subset of 162 known and potential auxiliary splicing factors was selected.

Quantile normalization and voom (Law et al. 2014) transformation were performed on gene read counts. Differential expression was analyzed using the limma package (Smyth 2005) , and P-values were corrected for multiple testing using the Benjamini-Hochberg method. Genes were considered differentially expressed if the log2-fold change > 0.5 in absolute value and the corrected P-value < 0.05. An expression Z-score per gene and per tumor sample was calculated using the quantile-normalized and voom-transformed read-counts.

The frequency of somatic mutations across all samples with available data was calculated per gene and tumor type. The association between expression regulation and mutations was measured using a Jaccard index. Mutual exclusion was measured using the number of samples having an RBP mutation and no expression change and the number of samples having expression change but no RBP mutation. For samples with CNV data available, the overlap of each RBP with the annotated CNVs was calculated, requiring a CNV score > log2(3) or < log2(1) for gain or loss, respectively. We required that the full locus of the RBP fall within a copy number region and defined focal CNVs to be those smaller than 5 Mb. Using a Jaccard index, an association was calculated between the expression up- or down-regulation and CNV gain or loss, respectively. More details are provided as Supplemental Material.

A total of 30,820 alternative splicing events were calculated from the gene annotation using SUPPA (Alamancos et al. 2015) : 16,232 exon skipping events, 4978 alternative 5′ splice-site events, 6336 alternative 3′ splice-site events, 1478 mutually exclusive exon events, and 1787 retained intron events. Differentially spliced events were obtained by comparing the event percent spliced-in value distributions between normal and tumor samples using a Wilcoxon signed rank test, removing samples with missing PSI values, using at least 10 paired samples, and correcting for multiple testing with the Benjamini-Hochberg method. Events were considered differentially spliced if the difference between the tumor and normal median PSIs > 0.1 in absolute value and a corrected P-value < 0.05. Events associated with mutations in RBPs were calculated in a similar way separating samples according to whether they had or did not have mutations in an RBP. Tests were performed using protein-affecting mutations. Only RBPs with mutations in at least 10 samples were tested. An enrichment Zscore was calculated per RBP and tumor type by comparing the number of events changing significantly, with the median value obtained using all RBPs tested. Tumor-type–specific alternative splicing events were calculated by comparing their PSI values for pairs of tumor types using an information theoretic measure. More details of the analyses are given as Supplemental Material.

Annotations for 50 cancer hallmarks were obtained from the Molecular Signatures Database v4.0 (Liberzon et al. 2015) , and a Fisher’s exact test was performed per hallmark using genes with annotated events and genes with differentially spliced events in each tumor type. The lists of 82 genes whose alternative splicing was linked before to cancer (Supplemental Table S5) and of 889 cancer drivers based on mutations and CNVs (34 in common with the previous set) (Supplemental Table S6) were collected from the literature (more details in the Supplemental Material).

Details on how motifs from RNAcompete (Ray et al. 2013) were assigned to the RBP genes are given in the Supplemental Material. The tool FIMO (Bailey and Elkan 1994) was used to scan the motifs in the event regions using a P-value < 0.001. Motif enrichment analysis was performed by comparing the frequency of regions in differentially spliced events with a motif with 100 random subsamples of the same size from equivalent regions in nondifferentially spliced events controlling for similar G + C content. An enrichment Z-score per motif, region, and direction of change (ΔPSI > 0.1, ΔPSI < −0.1) was calculated from the observed frequency and the 100 random control sets. A differentially spliced event was considered to be a potential target of a differentially expressed RBP if the correlation between the event PSI value and the gene expression robust Z-score was |R| > 0.5 (Spearman) and the event contained the RBP binding motif. To assess significance, the same number of differentially spliced events in a tumor type was randomly selected from all events 100 times and events associated with the RBPs calculated each time as described previously. A Zscore was calculated from the mean and standard deviation. Cases with Z-score > 1.96 were considered significant.

Networks of RBPs and events were built using the correlations between RBPs through events. For each pair of RBPs, a correlation was calculated using the Spearman correlation values with all differentially spliced events in the same tumor type. RBP clusters were built by calculating an inverse covariance matrix of these correlations using the glasso algorithm (Friedman et al. 2008) and then searching for dense, highly connected subgraphs with a greedy algorithm (Clauset et al. 2004) . Events were associated with a network if they had |R| > 0.8 (Spearman) or |R| > 0.5 plus motif for any of the RBPs in an RBP cluster.

Details on the cell cultures, siRNA transfections, Western blot analyses, and semiquantitative RT-PCR experiments are described in the Supplemental Material.

2′-O-Methyl RNA oligonucleotides were designed with full phosphorothioate linkage, antisense to the 5′ or 3′ splice sites of NUMA1 alternative exon 16 (hg19 coordinates Chr11: 71723447–71723488) optimizing GC content to 45%–60%. Custom modified and HPLC-purified RNA oligos were ordered in a 0.2 µM scale from Sigma-Aldrich:

NUMA1_ex16_5′ss: 5′-ggcauuacCUGCUUAGUUUGC-3′ NUMA1_ex16_3′ss: 5′-CCUCUAGCUGCUCCACcugu-3′ RANDOM 2′-O-Methyl RNA oligo: 5′-GCAAUGGCGU

CAAGUGUGUCG-3′

Antisense RNA oligos were transfected in triplicate at 20-nM final concentration using 2 µL Lipofectamine RNAiMax (Life Technologies, 13778150) per 1 mL of total volume of transfection in OPTIMEM (Life Technologies, 13778150). After 5 h of treatment, media was replaced by DMEM-F12 containing 10% FBS and Pen-Strep.

2500 MCF10A cells/well were seeded the night before treatment in 96-well plates (NUNC, 167008) in 100 µL complete DMEM-F12 medium. Wells with none, half, or double the amount of cells were also seeded for fluorescence calibration. Cells were transfected with siRNA or AON oligos as described. Resazurin (Sigma, R7017) treatment was performed 72, 96, and 120 h after transfection, in seven replicates and incubated for 4 h in a 37°C incubator. Fluorescence was measured after 4 h of incubation, using a TECAN infinite m200 device with 530-nm excitation wavelength, 590-nm emission wavelength, 30-nm emission bandwidth, and set to optimal gain. The medium was replaced by complete DMEM-F12 after measurements.

The number of centrosomes was determined by immunofluorescence assays using an anti-gamma-tubulin (TUBG1) antibody (clone GTU-88, Sigma-Aldrich; dilution 1:1000). The expected immunostaining pattern of this centrosomal marker in normal cells is one or two foci proximal to the nucleus. The cells were fixed in cold methanol for 10 min and washed in phosphate-buffered saline. The secondary antibody was Alexa Fluor 488 (Molecular Probes, Life Technologies), and the cells were mounted using VECTASHIELD with DAPI. The results correspond to at least five independent fields and >200 cells analyzed. The significance of the results was assessed using the one-sided Mann-Whitney U test (Supplemental Table S21). The chromosome instability signature CIN25 (Carter et al. 2006) was used by calculating the mean value of the normalized expression robust Z-score values for the 25 genes from the signature in each sample.