{'fraction_non_alphanumeric': 0.03637041872610277, 'fraction_numerical': 0.016133544716963537, 'mean_word_length': 5.196070878274268, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 5, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A comprehensive characterization of cis-acting splicing-associated variants in human cancer

Although many driver mutations are thought to promote carcinogenesis via abnormal splicing, the landscape of splicing-associated variants (SAVs) remains unknown due to the complexity of splicing abnormalities. Here, we developed a statistical framework to systematically identify SAVs disrupting or newly creating splice site motifs and applied it to matched whole-exome and transcriptome sequencing data from 8976 samples across 31 cancer types, generating a catalog of 14,438 SAVs. Such a large collection of SAVs enabled us to characterize their genomic features, underlying mutational processes, and influence on cancer driver genes. In fact, ∼50% of SAVs identified were those disrupting noncanonical splice sites (non-GT-AG dinucleotides), including the third and fifth intronic bases of donor sites, or newly creating splice sites. Mutation signature analysis revealed that tobacco smoking is more strongly associated with SAVs, whereas ultraviolet exposure has less impact. SAVs showed remarkable enrichment of cancer-related genes, and as many as 14.7% of samples harbored at least one SAVs affecting them, particularly in tumor suppressors. In addition to intron retention, whose association with tumor suppressor inactivation has been previously reported, exon skipping and alternative splice site usage caused by SAVs frequently affected tumor suppressors. Finally, we described high-resolution distributions of SAVs along the gene and their splicing outcomes in commonly disrupted genes, including TP53, PIK3R1, GATA3, and CDKN2A, which offers genetic clues for understanding their functional properties. Collectively, our findings delineate a comprehensive portrait of SAVs, novel insights into transcriptional de-regulation in cancer.


Introduction
Comprehensive genomic characterization of multiple cancer types in large-scale genetic studies has increasingly broadened the catalog of somatic alterations that dictate cancer evolution, including single nucleotide variants (SNVs), small indels (insertions and deletions), and copy number alterations [1][2][3] . Moreover, it has also revealed disturbances in transcriptional regulation, such as expression changes and splicing defects that underlie cancer pathogenesis 1,2 . However, there has been only a little progress in the understanding of how somatic alterations in cancer genomes exert direct transcriptional consequences.
In cancer transcriptomes, splicing defects play important roles in many aspects of cancer development and progression [4][5][6][7] . Discovery of somatic variants affecting RNA splicing factors, such as SF3B1 and U2AF1, which induce extensive alterations in RNA splicing (trans-acting regulation) in several kinds of cancers, highlights the relevance of RNA mis-splicing in cancer pathogenesis [7][8][9] . Another mechanism, which is the focus of this paper, is cis-acting regulation, in which somatic variants directly cause abnormal splicing of the affected gene. For example, somatic variants in canonical splice sites (highly conserved GT-AG dinucleotides at exon-intron boundaries) have long been reported to cause dysregulation of cancer-related genes 4,5 . These variants can induce different forms of abnormal splicing, such as exon skipping, intron retention, and activation of cryptic splice sites (SSs). Recent pan-cancer studies showed that SNVs causing aberrant intron retention in exon-intron boundaries are enriched in tumor suppressor genes (TSGs), especially TP53 10,11 . However, the complexity of splicing systems and the perplexing relationship between somatic variants and splicing alterations have limited the opportunities for systematic analyses of the extent and consequences of splicing-associated variants (SAVs), particularly those located other than at canonical sites (non-canonical sites).
To overcome these limitations, we have developed a novel algorithm, SAVNet

Overview of SAVNet framework
The overview of the proposed framework (SAVNet) is summarized in Fig. 1a. First, we collected evidence of abnormal splicing from tumor-derived RNA-seq data. Exon skipping and alternative 5′SS and 3′SS usage were extracted by capturing abnormal splicing junctions demarcated by split-aligned sequencing reads, whereas intron retention was identified by detecting sequencing reads spanning exon-intron boundaries (Fig. 1b). To obtain reliable and interpretable results, we focused exclusively on either (1) a somatic variant located at or close to an authentic exon-intron boundary (registered in the Refseq database), in which normal splicing is disrupted (SS disruption) or (2) a somatic variant located within a newly created SS inferred by an alternative SS usage event (SS creation). As demonstrated later, somatic variants and splicing alterations may have complicated relationships: a somatic variant occasionally generates different abnormal splicing events, whereas several different somatic variants sometimes cause the same splicing event. To represent these intricate relationships, we constructed a bipartite graph representing all potential associations between somatic variants and abnormal splicing events for each gene. Next, based on a probabilistic model for the number of abnormal splicing-supporting reads and the presence of a somatic variant, we deduced significant causal relationships through the evaluation of a Bayes factor incorporating a Bayesian model averaging framework 12,13 . A simulation study investigating the effect of the number of variant-splicing associations validated that the proposed framework can utilize the information from multiple associations for the sensitive identification of SAVs ( Supplementary Fig. 1a, b).
In the TCGA cohort, we compiled a total of 4,825,046 SNVs and 523,236 indels from 8,976 samples across 31 cancer types that underwent both whole-exome sequencing (WES) and RNA-seq using our in-house pipeline (Supplementary Tables 1 and 2 and Supplementary Methods). Initially, to determine the relevant positions within authentic SSs, we applied SAVNet to these sequencing data and assessed the accuracy of SAVNet for each position by calculating position-wise false discovery rates (FDRs) using a permutation of combinations of WES and RNA-seq data. Within authentic SSs, SS-disrupting variants at positions −3 through +6 of donor sites and −1 through +6 (except for position +4) of acceptor sites had low FDR values (below 20%), whereas much higher FDRs were observed at other positions (Fig. 1c). This observation prompted us to focus on somatic variants at these positions in the subsequent analysis. In addition, to control the overall FDR at these positions below 5%, we employed a threshold of !.! or greater for the Bayes factor, depending on cancer type ( Supplementary Fig. 1c). To evaluate the sensitivity of SAVNet under these settings, we compared our framework with a previous study consisting of six cancer types from the TCGA 10 . In the overlapping patient population (n = 929), SAVNet detected a markedly higher number of SAVs, including two-thirds of those found in the previous study ( Supplementary Fig. 1d, e). These results demonstrate the excellent detectability and satisfactory accuracy of SAVNet.

Landscape of SAVs in human cancers
With this optimized setting, we identified 14,438 somatic variants (13,414 SNVs and 1,024 indels) responsible for 18,036 splicing alterations in the TCGA samples ( Fig. 2a and Supplementary Table 3). A total of 11,153 SNVs and 875 indels disrupted splicing donor (n = 6,799) or acceptor (n = 5,229) motifs, of which 4,406 SNVs and 359 indels were not located within GT-AG canonical sites. In addition, 2,261 SNVs and 149 indels were detected to create novel splicing donor (n = 1,566) and acceptor (n = 844) sites. Thus, 7,175 (49.7%) somatic variants would not be expected to be identified by conventional methods that concentrate on SAVs involving canonical sites. Although the number of SAVs per sample was generally low (median of 1), there were quite a few samples with more instances of SAVs, particularly in cancer types with high somatic variant rates, such as lung and skin cancers (Fig. 2b).
Although the vast majority of SAVs caused a single splicing alteration, 2,778 (19.2%) variants induced multiple splicing alteration events ( Fig. 2c and Supplementary Fig. 2). The transcriptional consequences substantially differed according to the somatic variant pattern (donor vs. acceptor and disruption vs. creation). Exon skipping and intron retention were caused by variants disrupting both donor and acceptor sites (Fig. 2a). As expected, donor disruptions tended to generate an alternative 5′SS, whereas acceptor disruptions more frequently gave rise to an alternative 3′SS. Interestingly, exon skipping was the most frequent consequence of donor disruptions, whereas alternative 3′SS accounted for more than one-half of acceptor disruptions. Many new splice donor and acceptor sites were created by variants outside authentic SSs. Aberrant splicing events associated with variants in trans-acting splicing factors 7 showed no overlap with those detected by SAVNet (Supplementary Tables 4 and 5). canonical GT-AG sites (at positions +1 and +2) had the highest ratios of splicing aberrations (18.3%-24.0%). In donor SSs, non-canonical sites showed a comparable total number of SAVs with canonical sites, whereas most SAVs in acceptor SSs were present at canonical sites. Together with the last exonic bases (−1) of donor sites, whose relevance was pointed out in the earlier study 10 , the fifth intronic bases (+5) also had a relatively high ratio of abnormal splicing, followed by the third intronic bases (+3). Besides GT dinucleotides, these bases are well conserved and relevant to the interaction with U1 and U5 small nuclear RNAs 14,15 . In fact, using minigene splicing assays 16 , we experimentally demonstrated that not only canonical but also non-canonical site variants cause abnormal splicing (Fig. 3b). The transcripts harboring variants at positions +5 as well as −1 showed abnormal splicing, such as exon skipping or intron retention, with comparable efficiency to canonical site variants (+1), while the wild-type transcripts were largely normally spliced.

Mutational signatures associated with SAV generation
There were occasional discrepancies between the efficiency of somatic variants to cause abnormal splicing and the actual number of SAVs. For instance, position +2 of acceptor sites showed only a moderate number of SAVs, albeit the highest SAV ratio (Fig. 3a). These discrepancies may be attributed to the overall number of somatic variants (including those not associated with splicing alterations) and their substitution patterns at each position, which reflect both the unique base composition at SSs and the signatures of mutational processes that have been operative 17,18 . In fact, positions at −1, +1, and +5 of donor sites as well as +1 of acceptor sites, which were dominated by G bases, showed frequent G > A and G > T substitutions, suggestive of age-and smoking-related mutational processes, respectively (Fig. 3c,upper). In contrast, position +2 of donor and acceptor sites, which predominantly consist of A/T bases, showed a relatively low frequency of somatic variants. Among them, variants at canonical GT-AG sites caused splicing alterations, regardless of their base substitution pattern, whereas almost all SAVs at positions −1 and +5 of donor sites occurred at G bases, indicating almost no effect of substitutions from other bases on splicing (Fig. 3c, lower). In addition, positions having a smaller fraction of abnormal splicing were more strongly affected by the base substitution pattern. For example, G > A substitutions were the most common at position +3 of donor sites, but did not result in splicing aberrations. Moreover, despite their low frequency of overall variants, C > G substitutions (compatible with the APOBEC cytidine deaminase mutational pattern as shown below) accounted for a considerable proportion of SAVs at position +3 of acceptor sites. These findings are consistent with relatively limited conservation of splicing motifs at these positions.
To evaluate the underlying mutational process for each SAV, we analyzed mutational signatures found in the current sample set using pmsignature algorithm 17 . Among the five major mutational signatures (processes generating a large number of somatic variants) ( Supplementary Fig. 3b), smoking signature showed the largest contribution to SAV generation, followed by APOBEC and aging signatures. Signatures related to ultraviolet exposure and altered activity of the error-prone polymerase POLE had less impact ( Fig. 3d and Supplementary Fig. 3c). These differences can partly be explained by the predominance of G bases at highly affected positions (−1, +1, and +5 of donor sites and +1 of acceptor sites) and the transcriptional strand bias of several mutational signatures, i.e., the smoking signature preferentially affects G bases, whereas the ultraviolet signature frequently alters C bases. Reflecting these differences among mutational processes, lung squamous cell carcinomas (LUSC) had more SAVs than expected from the overall somatic variant rate, whereas cancers frequently affected by POLE alterations, such as uterine corpus endometrioid carcinomas (UCEC) and colon adenocarcinomas (COAD), as well as ultraviolet-associated skin cutaneous melanomas (SKCM) showed a relatively lower number of SAVs ( Fig. 3e and Supplementary Fig. 3c).

Characteristics of SAVs creating alternative SSs
Our analysis also revealed the positional distribution of SAVs creating alternative donor and acceptor sites. Newly created donor sites were widely distributed in both exons and introns, whereas abnormal acceptor sites were created predominantly within the polypyrimidine tract ( Fig. 3f, upper), likely reflecting the involvement of additional conserved elements in introns, such as branchpoint sequences and polypyrimidine tracts 14,15 . Apparently similar distributions were also seen for cryptic SSs activated by variants disrupting the authentic SSs (Fig. 3f, lower). However, unlike newly created acceptor sites, a biased localization of cryptic acceptor sites toward exons was observed, which can be plausibly explained by a depletion of AG dinucleotides in the polypyrimidine tract.
We also evaluated the substitution pattern of somatic variants creating new splicing sites based on their relative position within the newly created SSs. Most newly created donor sites resulted from GT canonical site generation through C > T substitutions at position +2, whereas variants associated with acceptor creation tended to form a new YAG (Y, pyrimidine) motif at positions +3 through +1 ( Fig. 3g and Supplementary Fig. 3d). These results suggest that, showing a strong bias toward particular base substitutions, these SAVs generate additional consensus donor or acceptor motifs that are more efficient for splicing than those in authentic SSs, as implicated by stronger splicing strength (assessed by MaxEnt 19 or H-bond 20 scores) (Supplementary Fig. 3e-g).

Features of genomic sequences associated with SAVs
Splicing outcomes mediated by SAVs appear to be context-dependent: Somatic variants within authentic SSs can cause different forms of splicing aberration, while the same substitutions at the same relative position frequently do not alter splicing. To elucidate the factors determining the potential of somatic variants within authentic SSs to alter splicing, we compared the genomic features of SSs between normal (those not identified as SAVs) and abnormal splicing groups (SAVs) (Supplementary Fig. 4a). Generally, SS-disrupting SAVs attenuated the splicing strength more than variants that induced no abnormal splicing, regardless of the substituted position and consequent splicing alteration type (

Enrichment of SAVs in TSGs
To evaluate the role of SAVs during cancer development, we investigated which genes are frequently altered by SAVs. Strikingly, the majority of frequently affected genes (present in ≥ 10 samples across the entire cohort) were well-established TSGs (Fig. 5a, b and Supplementary Fig. 6a, b). In consistent with this study 10 , in which intron retention was argued to be a major mechanism of SAV-induced TSG inactivation, SAVs that caused intron retention showed the strongest enrichment of TSGs, regardless of the cancer gene sets 1,25,26 . However, SAVs associated with exon skipping and alternative SS usage also had a greater proportion of TSGs, even when compared with nonsense variants, accounting for 88% of SAVs affecting TSGs (Fig. 5c). These findings suggest that, together with intron retention, exon skipping and alternative SS usage play crucial roles in TSG inactivation. In contrast, oncogenes were less frequently affected by SAVs, comparable to missense variants. In total, 1,684 SAVs in candidate cancer-related genes 25 were identified in 14.7% of the TCGA samples (1,315 of 8,976). Particularly, as many as 914 SAVs found in 9.3% of samples targeted well-known TSGs 1 , of which 341 were not located at canonical sites. Moreover, SAVs accounted for 9.5% of loss-of-function variants in these genes. Therefore, SAVs represent an important but previously underestimated mechanism for TSG inactivation, irrespective of splicing outcome.
Like nonsense variants, splicing alterations are thought to trigger nonsense-mediated decay (NMD), a surveillance mechanism that selectively degrades abnormal transcripts containing a premature termination codon 6,10 . To clarify the effects of SAVs on gene expression through NMD in TSGs, we investigated the whole transcript level across different types of abnormal splicing. In line with the previous report 10 , transcripts with intron retention showed a substantially lower expression level than normal transcripts, which was comparable to those with nonsense variants (Fig. 5d). The expression of transcripts with exon skipping or alternative SS usage was also reduced when their splicing alterations caused frameshift changes.

Genes frequently altered by SAVs
Among genes frequently targeted by SAVs, TP53 was the most frequently altered gene, affecting 233 samples in 22 cancer types (Fig. 5a). Although the last bases of exons 4, 6, and 9 were reported to be frequently mutated 10,11 , we identified a number of recurrent variants at splice donor and acceptor sites of introns 3 through 9, with prominent base-level and/or SS-level hotspots at donor and acceptor sites of intron 4 (Supplementary Tables 6 and 7).
Approximately one-half of recurrent SAVs simultaneously produced different types of abnormal splicing, while identical abnormal splicing events were generated by different SAVs, such as retention of introns 7, 8, and 9 caused by donor and acceptor SAVs of each intron ( Fig. 6, upper left). Most of these SAVs induced frameshift splicing alterations, likely leading to mRNA degradation through NMD, whereas other SAVs generated in-frame exon skipping or alternative SS usage, such as exon 5 acceptor variants activating cryptic 3′SS, followed by 15-amino acid deletion. PIK3R1, encoding the p85 regulatory subunit of phosphatidylinositol 3-kinase, ranked second (39 samples), approximately one-half of which were found in UCEC (Fig. 5a). The majority of these SAVs caused in-frame exon skipping (mainly involving exon 4) resulting in a deletion of the iSH2 domain, which is also affected by the small deletions typically observed in this gene 27 (Fig. 6a, upper right).
In many genes, particularly NF1 and RB1, most SAVs and consequent splicing alteration events were diverse and widely distributed throughout the entire gene ( Supplementary Fig.   7a, upper and middle), whereas several genes displayed prominent hotspots of SAVs (Supplementary Tables 6 and 7). Among the latter genes, SAVs affecting the same SSs tended to generate identical splicing consequences ( Fig. 6b and Supplementary Fig. 7b). A typical example was CDKN2A, a well-known TSG that encodes the p16 INK4A and p14 ARF proteins, which was recurrently affected by SAVs targeting exon 2 common to both proteins.
Other instances included GATA3 SAVs found in breast invasive carcinomas (BRCA), in which most of them were the identical CA dinucleotide deletion at the acceptor site of exon 5, thus activating a cryptic 3′SS (7 nucleotides downstream), (Fig. 6a, lower left). Utilization of this cryptic splice acceptor caused a reading frameshift, resulting in loss of the second zinc finger (ZnF2) domain 28 . As was the case with GATA3, several genes showed tissue specificity, such as FUBP1 and ATRX in lower grade gliomas (LGG), although most of the frequently altered genes were relevant across many cancer types (Fig. 5a).
Together with these well-established TSGs, SAVNet identified several recurrently altered genes (found in ≥10 samples) which had not been included in the cancer-related gene list 25 but reported or predicted to function in a tumor suppressive manner, including KANSL1 29 , NIPBL 30 , CUL3 31 , MYH9 32 , SMCHD1 33 , and HUWE1 34 (Fig. 5a). Thus, SAVNet may have potential to identify putative TSGs that are more prone to be affected by splicing aberrations.
Conversely, MET, which encodes a hepatocyte growth factor receptor, was the only frequently affected oncogene, whose variants in the exon 14 donor site caused in-frame exon skipping known to activate c-Met 35 ( Fig. 5a and Supplementary Fig. 7a, lower left).
Additionally, SAV hotspot analysis also identified recurrent SAVs occurring at the donor site of exon 3 of MIEN1, a putative oncogene located on the HER2 amplicon 36 (Fig. 6b and Supplementary Fig. 7a, lower right). Although the underlying mechanisms need to be clarified, SAVs may contribute to the activation of several oncogenes.

Discussion
The development and application of SAVNet have led to the systematic detection of a substantial number of SAVs that had been overlooked by earlier studies 10,11 , although we focused only on those disrupting or creating splicing donor or acceptor motifs. Following previous studies 10,11,37 , our comprehensive and thorough analysis revealed the landscape of cis-acting somatic variants affecting splicing, characterized their positional differences and genomic features as well as underlying mutational processes in details, and argued that they include driver mutations, especially those inactivating TSGs. In particular, we demonstrated that exon skipping and alternative SS usage were more frequently involved in SAV-mediated TSG inactivation than intron retention. In addition, we clarified the relevance of SAVs at non-canonical sites, including the previously unrecognized position +3 and +5 of donor sites.
The proposed framework with FDR control, which can dissect complex variant-splicing associations based on the Bayesian approach, is applicable to identify additional classes of somatic variants that disrupt splicing regulatory elements, including exonic and/or intronic splicing enhancers and silencers, although further elaboration of association rules will be required. Based on our findings, not only exonic but also intronic SNVs near exon-intron boundaries should be carefully evaluated as pathogenic variants, irrespective of the presence of amino acid changes. In the era of precision medicine, our framework will be critical to capture all driver variants, including previously overlooked SAVs, in cancer patients.

Detection of somatic SNVs and short indels
Our approach for detecting somatic SNVs and short indels consists of the following five steps: (1) Identification of candidate somatic SNVs and short indels using the approach based on Fisher's exact test (as previously described 8 ), which is currently implemented in GenomonFisher (https://github.com/Genomon-Project/GenomonFisher).

Quantification of expression values for each gene from RNA-seq data
To quantify gene expression, we used our in-house software GenomonExpression

Collection of evidences of different types of abnormal splicing
We consider four types of abnormal splicing: exon skipping, alternative 5′ splice site (SS), alternative 3′SS usage and intron retention. The first three types (exon skipping, alternative 5′SS, alternative 3′SS) are extracted using splicing junctions (defined as pairs of start and end positions demarcated by spliced-aligned reads). We first extract abnormal splicing junctions (not registered in RefSeq genes (http://hgdownload.cse.ucsc.edu/goldenPath/hg19/database/refGene.txt.gz)) with ≥ 2 supporting reads (number of uniquely mapped reads crossing the junction) in at least one sample in the cohort by processing SJ.out.tab files generated as by-products of the STAR alignment step. Then, using our in-house program (junc_utils, https://github.com/friend1ws/junc_utils), we classify each splicing junction into exon skipping, alternative 5′SS or alternative 3′SS by the following criteria.
• Exon skipping: Two ends of the splicing junction correspond to annotated intron start (splicing donor) and end (splicing acceptor) sites of a gene, respectively.
• Alternative 5′SS: One end of the splicing junction corresponds to an annotated intron end (splicing acceptor) site of a gene, whereas the other end is located within the gene, but not at an annotated intron start (splicing donor) site of the gene. ii. The average number of supporting reads for the -th splicing alteration in samples with the -th somatic variant is at least three times larger than those in samples without any somatic variants of the gene in consideration: iii.
The median number of supporting reads for the -th splicing alteration in samples without any somatic variants of the gene in consideration is zero.
We For each configuration , we assume that the supporting reads ! are generated by Poisson distributions whose parameters are dependent on the activity status of somatic variants and multiplied by sample weights. The parameter of the Poisson distribution for the n-th sample is set to ! ! when it has only inactive variants on the -th splicing alteration Additionally, we assume that ! and ! are generated by Gamma distribution with shape and rate parameters ( ! , ! ), ( ! , ! ), respectively. In this study, we set ( ! , ! ) = (1,1) and ( ! , ! ) = (1,0.01). Therefore, the likelihood of ! given is and the likelihood of whole the data is Pr ! ( | ) = For each variant , we perform Bayesian model comparison to determine whether the somatic variant has any causal relationships with any splicing alterations (∃ , !,! = 1) or not (∀ , !,! = 0). Typically, to perform model comparison, we evaluate Bayes factor between the two distinct models. Here, as there are often many distinct null and non-null models, we aggregate these models through Bayesian Model Averaging and evaluate the Bayes factor between the aggregated null and non-null models: where Pr ( ) is set to be the uniform distribution. The variant is identified as a SAV if its logarithm of Bayes factor is above the threshold (the default value is set to 3). Also, the splicing alterations caused by the variant are identified by selecting the best model by

Evaluation of FDR by permutation
To evaluate FDR, we permute the pairs of genomic and transcriptome data so that somatic variants and splicing alterations from different patients are coupled, and perform the same procedures (step 1 to 3). Assuming that !"#$%! and ! !"#$ ( = 1, ⋯ , ) are the numbers of SAVs identified in the original step (correct combinations) and in the -th permutation procedure, respectively, then FDR is estimated as In this paper, we performed 100 permutation trials ( = 100).

Post-processing and rescuing SAVs
SAVs causing alternative intronic 5′ or 3′SSs are generally accompanied with intron retention at the original authentic SSs. Therefore, in these cases, we removed intron retention and retained only alternative 5′SS or 3′SSs in this paper. To sensitively detect recurrent SAVs, we performed additional screening and adopted variants satisfying the following criteria: (A) the combination of the same somatic variants (the same substitution at the same position) and the same splicing alterations was identified in other samples by the SAVNet procedure

Evaluation of influences of spliceosome variants on abnormal splicing
First, in the TCGA cohort, we searched for previously known somatic variants of splicing factors, including missense variants at K700, K666, H662, R625, E622, G740, G742, N626, and E902 of SF3B1, S34 and Q157 of U2AF1, and P95 of SRSF2, as well as missense, nonsense and frameshift variants of ZRSR2 7 . First, for each cancer type, we extracted splicing alterations with ≥ 2 supporting reads in at least one sample. Then, we identified splicing factors affected in ≥1% samples within each cancer type, and compared the number of RNA-seq reads supporting each splicing alteration between samples with and without the splice factor variants to derive P-value using t-test. Finally, we calculated Q-value for each splicing alteration using qvalue R package and splicing alterations with Q-value < 0.05 were considered to be associated with splice factor variants.

Estimation of mutational signatures and membership of SAVs
We used pmsignature for estimating the signatures of mutational processes operative in the entire cohort and each cancer cohort as described in the previous paper 17 . Then, the extracted mutation signatures were classified into any of COSMIC signatures (http://cancer.sanger.ac.uk/cosmic/signatures) using minimum centered cosine similarity.
Mutation signatures with centered cosine similarities to all the COSMIC signatures < 0.75 were classified to "other." The estimates of membership (conditional probabilities attributed to each mutation signature) for each variant are provided by the following equation: where the notation for each variable is described in the previous paper 17 . Finally, variant-level membership estimates were aggregated according to the COSMIC signatures and the presence of association with splicing alterations, so that the total numbers of variants and SAVs caused by each mutation signature (e.g., tobacco, ultraviolet) were estimated.

Cell line
HEK293T cells were obtained from RIKEN Cell Bank. Cell lines were authenticated by the provider and routinely tested for mycoplasma infection.

Minigene splicing assay
For each region on interest, exonic and ~300 bp flanking fragments containing either wild-type or variant sequences were synthesized (GeneArt, ThermoFisher Scientific) and cloned into the NheI and BamHI sites of the plasmid H492 (a kind gift from Prof. Masafumi Matsuo, Kobe University) using In-Fusion HD cloning kit (TaKaRa). Each construct was transiently transfected into HEK293T cells using X-tremeGENE 9 DNA Transfection Reagent (Roche) in six-well tissue culture plates. Forty-eight hours after transfection, total RNA was isolated with RNeasy Mini kit (QIAGEN), and used to synthesize cDNA with ReverTra Ace qPCR RT Kit (TOYOBO). Each cDNA was amplified with KOD FX Neo DNA polymerase (TOYOBO) using the primers (forward, 5′-ATTACTCGCTCAGAAGCTGTGTTGC-3′, and reverse, 5′-AAGTCTCTCACTTAGCAACTGGCAG-3′), which correspond to sequences of exons of H492. PCR products were separated by electrophoresis on 2% agarose gel and visualized with a UV transilluminator (UVP). To confirm the sequence of each band, the PCR products were gel purified and analyzed by Sanger sequencing.        SSs (f) among the splicing groups. The significant differences were quantified by the negative logarithm of the P-value of one-sided Wilcoxon rank-sum tests whose alternative hypotheses are that either the splicing pattern 1 is less (blue) or greater (red) than splicing pattern 2. The smaller P-value between the two alternative hypotheses is shown. Schematics depicting methods to quantify the fraction of the most frequent relative to total associated splicing outcomes for each SS-level SAV hotspot. The number at each line represents the count of supporting reads.