Pan-cancer repository of validated natural and cryptic mRNA splicing mutations

We thank the reviewer for their valuable comments. Our responses follow:
 
1. The authors write that their IT-based framework to predict slicing variants "has been validated in multiple studies" and they refer to numerous papers. However, in all of them they act as co-authors, showing that their method is mostly used by themselves and their collaborators. This is not necessarily a problem, but it would certainly strengthen the study if the authors would perform a comparative assessment of their performance with other available methods to predict splicing mutations, for example those in dbNSFP. This will provide a less biased interpretation of the final results.
 
Response: We have previously compared our mutation prediction methods with others. Mucaki et al. Hum. Mut. 34:556-65 (2013) showed that IT-based single mutation and exon definition methods performed as well or better than MaxEntScan and Human Splice Finder websites. MaxEntScan computes relative entropy which is similar to IT, except it applies a correction for local base composition (which does not measure free energy, in contrast with IT: J. Theor. Biol. 201:87-92, 1989). Human Splice Finder does not measure changes in binding affinity and its basis is ad hoc. We also reviewed all articles (300+) which have used IT-based tools to predict changes in splicing (Caminsky et al., F1000Research 3:282, 2014). This reference covers the vast majority of studies that used IT-based bioinformatic tools for mutation analysis and, compared results obtained using these tools with other available software. The cited studies included a large proportion of mutations that we, ourselves, did not coauthor, or were analyzed by others, removing an obvious source of bias.
 
Regarding potential bias in our results, the IT-based position weight matrices (iPWMs) of splice recognition sites that we derived and use are based on a comprehensive set of splice sites spanning all known coding genes (see appendix of Rogan et al. Pharmacogenetics & Genomics 13(4):207-18, 2003). Other bioinformatic methods for splice site detection are based on many fewer splice sites for PWMs and are much more likely to be subject to bias based on how those sites were chosen. Also, the determination of information content in natural and mutated splice sites obeys the second law of thermodynamics (Schneider, J. Theor. Biol. 189:427-41, 1997); information contents have been formally proven to be related to binding affinities of splice site to splicesomes and splicing factors. MaxEntScan differs from IT because it applies a correction for local base composition, which is energetically and biochemically irrelevant to binding site affinity, and is therefore, biased.
 
As comparisons of IT with other methods have been covered previously, and the F1000Research Data Note article format is intended to specifically present results using the new resource we describe, in our opinion, reevaluation of other algorithms would not add anything of value to this manuscript.
 
2. Somehow related to the previous point, the authors mention that their results "contrast with another TCGA study that investigated alternative mRNA splicing". In my opinion this point should be further explored: what are the main differences and what is the extent of overlap in concordant predictions? What are the possible reasons for these differences? This is important because the cited paper in Cancer Cell analysed the same dataset of mutations.
 
Response: The paper being referred to in the reviewer’s comment is Kahles et al., Cancer Cell., 34(2):211–224.e62018. This paper reports: “In a joint analysis of cis and trans associations with 50% prior on each type, we identified 32 cis- and seven trans-sQTL (Bonferroni corrected p < 0.05).” Information regarding these cis- sQTLs can be found here:
https://api.gdc.cancer.gov/data/3920a044-6874-4049-8010-55f0922243b7.
 
However, the document that provided does not indicate that these are actual splicing mutations. The document contains SNP coordinates, but not sequence changes. An assessment of known rsIDs at these coordinates only accounted for only 7 of the variants in our database. In fact, none of the substitutions (nor the rsIDs) in the article could be related back to changes in mRNA splice strength. It is not possible to comparison the results presented in this paper to those in Kahles et al. (2018).
 
 
3. The authors notice that the number of variants which activate cryptic splicing exceed the number reported in a recently published study in Cell Reports. Similarly to before: what is the extent of overlap between the two datasets? Stating that a dataset is bigger than another one is not necessarily an indication that it is better.
 
Response:  Jayasinghe et al. Cell Rep.;23(1):270-281.e3 (2018) identified 2056 variants in TCGA patients which activated cryptic splice sites, of which 1964 were confirmed by manual review of RNAseq data using the Integrated Genome Viewer (IGV). The data provided (Supplementary Table S1; “Passed Manual Review” tab) was sufficient to perform a comparison with our results. We scanned their manually reviewed variants using the Shannon Pipeline (SP). Despite reporting 1964 manually reviewed mutations, 50 mutations are shared between multiple patients (there are 1914 unique mutations total). 1510 of these mutations were found to alter at least one mRNA splice site (natural or cryptic), of which 1176 met our SP filtering criteria (either decreased natural site strength or cryptic splice sites strengthened by ≥2 bits and exceeding the strength of the nearest natural site of the same polarity).
 
We considered the possibility that splicing mutations in Jayasinghe et al. that were not flagged by SP could have instead altered the strength of splicing regulatory factor binding sites (SRFs). The high-throughput IT-based variant analysis tools needed to address this question were not available at the time the TCGA genomic data were processed. We recently introduced a new version of SP which is capable of analyzing variants impacting both constitutive splice sites and SRFs (including SRSF1, SRSF2, SRSF5, SRSF6, hnRNPA1, ELAVL1, PTB and TIA1). Upon analysis of the variants in Jayasinghe et al. for IT-derived changes in natural, cryptic and SRF binding site strengths, 1746 of the 1914 (91.2%) were found to be significant. It is conceivable that remaining unclassified variants may affect binding by splicing factors for which we have not yet derived iPWMs (i.e. SRSF7).
 
Of the 1176 variants meeting our filtering criteria, 824 were flagged by Veridical (70.1%). Interestingly, 27 Veridical-flagged mutations alter the same genomic coordinate in different tumours and another 64 affect the adjacent genomic coordinate within the same splice site in other individuals. We further investigated the remaining mutations (those evaluated, but not flagged) to better understand the discrepancy. In approximately 10% of cases, Veridical did not find any alternative validating cryptic splicing events in the region that contrasted with the read distribution in the set of control transcriptomes. For example, chr9:35389842G>C in TCGA-CN-5370 was expected to abolish the natural acceptor site of UNC13B exon 24, however the read counts for intron inclusion in the RNAseq data were too sparse to be deemed significant (p=0.33).
 
Jayasinghe et al. also found mutations in TCGA patients that were not evaluated in our study. This could occur for either because the RNASeq BAM file for a particular TCGA patient failed to download, or the “key” file that associates the BAM file to its TCGA name was incomplete in some tumors (i.e. TCGA-CG-4436; TCGA-STAD), due to both the BAM file name and header lacking this information. This would not impact the accuracy of the data that is present in our database or that we report, only the level of concordance between our results and those of Jayasinghe et al.
 
When comparing this dataset with our own, we discovered an instance where discrepant RNAseq data for the same tumour in the same TCGA patient led to different IT results. Originally, Veridical did not find a significant splicing change in POLR1B exon 14 at the +1 position, a G>A mutation g.113331138G>A, present in patient TCGA-Z6-AAPN in the ESCA dataset. When we manually reviewed the RNAseq BAM file used in this analysis, alternate or exon skipped splice forms were not observed, confirming the results reported from Veridical. In fact, TCGA deposited two separate BAM files containing RNAseq data for this same patient: “TCGA-Z6-AAPN-01A-11R-A406-31_rnaseq.bam” and “UNCID_2681450.b17c4505-8a84-4cdd-8782-fbc456deb2a6.sorted_genome_alignments.bam”. The latter BAM file shows evidence of both predicted exon skipping and activation of a cryptic acceptor site 12 nt downstream of exon 14 of POLR1B. Comparison of SNPs between these BAM files in other genes did not show any evidence of sample switching or contamination. Although this issue does not appear to be widespread in the TCGA dataset, such discrepancies exceed the scope of our study, and rightfully should be addressed by TCGA.
 
Nevertheless, this discovery did prompt us to re-analyze the TCGA-ESCA variants through Veridical using the second set of BAM files, and we have now included those results to ValidSpliceMut. In this new set, the previously mentioned POLR1B mutation is deemed significant due to increased exon skipping (86 reads showing exon skipping; p=0.0000).
 
4. The authors validate ~27% of predicted splicing variants using the mRNA data (351k validated of the 1.2M predicted). This is a surprisingly low fraction. Later in the manuscript, the authors briefly discuss about the possible reasons of such a discrepancy. One of them is the possible occurrence of nonsense mediate decay which will not confirm the mutations because no or very few reads will be detected. However, as the authors acknowledge, the absence of supporting reads only in mutated individuals as compared to the presence of reads in WT sample would be a strong indication of the effective role of these mutations on splicing. This can be quantified from the same RNAseq data and in my opinion should be done.
 
In the revision of this paper, the fraction of validated predicted splicing variants in the ValidSpliceMut database increased from 27% to 31%, as a result of a series of improvements in the software used for processing. SP was significantly upgraded over the course of this project,. The previous iteration would incorrectly report information changes at pre-existing splice sites adjacent to certain mutations. These sites were characterized by genomic coordinates that were altered by insertions or deletions (indel), regardless of whether the site overlapped or included the sequence change. In such instances, some altered natural splice sites could be designated as cryptic sites. SP also reported changes at cryptic acceptor and donor sites in the first and last exons of a gene, respectively, which were not likely to have a meaningful impact on splicing. Therefore, all datasets processed prior to this upgrade were reanalyzed for indel variants. The TCGA ESCA dataset was reprocessed with a second set of RNAseq BAM files (see above response to point 3), which increased the fraction of flagged mutations for that tumor type. Finally, we processed an additional 7 tumor datasets from ICGC and included the validated mutations in our primary beacon database. The statistics of mutations, their distributions and support have been updated in the present version of the manuscript.
 
From our perspective, the proportion of variants validated is not “a surprisingly low fraction". The results reported here are consistent with our previous published studies (Dorman et al., Sci. Rep. 4: 7063, 2014). Aside from NMD, as demonstrated below, some mutations that significantly alter splice site strength may not have been flagged by Veridical as a consequence of low levels of expression of the gene in the tumor (or controls) itself. Furthermore, Veridical cannot make an accurate assessment of the region of interest in control samples if these lack sufficient read abundance levels to determine the probability (p-value) of observing expression in the mutation-containing vs control samples. Also, our analysis did not take into account other impacts of the variant on sequences that influence exon definition, such as binding to splicing regulatory factors or mRNA secondary structure. In our response to point 3, we described a discrepancy in BAM file sources, which could also lead to a lower fraction of confirmed variants. Finally, miscalled variants (despite the stringent quality control criteria applied to select variants) could contribute to the fraction of variants not supported by Veridical analysis. Such technical artifacts have been shown to be quite common in exome sequencing in areas of the genome characterized by low mappability (from Shi et al., Cell Rep. 2018 Nov 6;25(6):1446-1457, 2018): “Examination of the genomic locations of mutations revealed that 41.1% of the artifactual somatic mutations occurred in regions of low mappability compared with only 6.4% for the validated somatic heterogeneous mutations.”
 
As indicated in our response to point 3, the previous version of this paper only evaluated variants for their impact on constitutional mRNA splice sites and cryptic sites, and excluded impacts of mutations at splicing regulatory factor binding sites (SRFs). The scope and time required to assess SRFs precluded the reanalysis of all datasets for such changes. However, to address the issue raised by the reviewer, we evaluated the degree to which ignoring SRFs would affect the overall discovery of splicing-related variants. The updated version of the Shannon Pipeline (SP) with this capability was used to examine 1050 mRNA splicing variants that have been demonstrated to affect exon recognition (Cheung et al., Mol Cell. 2019 Jan 3;73(1):183-194.e8). These splicing variants were experimentally validated using a high throughput, multiplexed splicing minigene reporter assay in that study. SP reported a change in splicing and/or SRF binding strength for 1017 of these 1050 mutations (96.9%; where change in SRF strength was ≥ 3 bits). After accounting for SRF location (e.g. exonic TIA1 sites were eliminated, since these have not been proven to have splicing effects; see Table 1 of Caminsky et al., F1000Res.;3:282, 2014, for a full description of each SRF), the number of flagged variants was reduced to 940 (89.5%). Based on changes in constitutive splice site strength alone, 447 variants were flagged (435 weaken natural sites, 14 strengthen cryptic sites to a level exceeding the nearest natural site). Therefore, 46.3% of the constitutive mutations at natural or cryptic splice sites were also flagged by SP. This suggests that the lower predictive accuracy of SP in our original submission was, in part, due to the limitations in its ability to detect pathogenic mutations in SRF binding sites.  
 
We addressed the reviewer’s suggestion to compare expression of the same gene in tumours with Veridical-validated mutations with other tumors with SP-predicted mutations in the same gene that were not experimentally validated.  To perform the analysis, we obtained pre-processed mRNA expression data from the same RNAseq sources of TCGA patients from cBioPortal (www.cbioportal.com; provisional datasets were used, which contained largest number of patients for each tumor type). We extracted these gene expression values with a software program we wrote that determined transcripts per million (TPM) for each gene containing a SP-flagged variant.  Expression values for gene present multiple times (due to the presence of multiple splice isoforms) were averaged for the particular tissue from which they were derived.
 
We separated mRNA expression values for each gene in TCGA patients into the Veridical-flagged vs. non-validated SP-predicted mutation categories, and performed a Student’s t-test on the two groups. The expression values of 58.2% of genes were statistically distinct with 90% confidence; >2 patients per category per gene). With at least 10 patients per category, the number of statistically different genes increased to 69.3%. Among these genes, patients with Veridical-flagged variants had higher overall gene expression in 99.7% of cases. These inherent differences in expression suggest that the failure to validate predicted mutations may be related to little or no expression of these genes in tumour and/or control samples, rather than to accuracy of IT prediction methods. Non-sense mediated decay could be responsible for the decreased expression of these mutated genes in the tumour genomes that carry them, or the failure to validate could be related to low levels of expression of these genes in the particular tissues from which the tumours were derived. This analysis is now described in the revised manuscript (‘Dataset validation and discussion,’ para. 8).
 
Note: The gene expression data from cBioPortal had some limitations: 4196 genes containing variants flagged by SP are not present in the mRNA expression datasets, though the vast majority occurred in non-coding RNAs (i.e. 145 microRNAs, 194 LINC RNAs) or other uncharacterized RNAs (e.g. 324 ‘LOC’ RNAs). Furthermore, certain TCGA patients that we analyzed were not available in cBioPortal among the available expression datasets (2 TCGA-BRCA, 18 TCGA-COAD, 19 TCGA-GBM, 1 TCGA-HNSC, 1 TCGA-LUAD, 119 TCGA-OV, 4 TCGA-READ, 4 TCGA-STAD, 4 TCGA-THCA, and 7 TCGA-UCEC patients). Nevertheless, sufficient data were obtained for the analysis that we carried out.
 
Furthermore, this analysis has other caveats, especially regarding the accuracy of RNAseq quantification between samples, even upon normalization of the data. The vast majority of TCGA tumor samples do not have a matched normal counterpart, and many TCGA tumor datasets (e.g. TCGA-LAML) do not provide any normal control RNAseq data at all. While we could perform an analysis in which the tumor is compared to a set of normals for the same tissue, variation in expression between different individuals would obfuscate evidence of any apparent NMD caused by the mutation. A 50% reduction of expression (or less) may not be observable under these conditions.
 
 
5. In general, the authors seem to exclude that their prediction method could lead to false positives. Rather they justify the poor overlap with limitations of mRNA detection. If this is the case, this should be quantified and probably a comparison with other prediction methods could help.
 
Response: We have quantified RNA expression in tumors for mutations that were validated by Veridical vs those which were not (see response to point 4). Regarding false positives: An extensive comparison of Information Theory to other bioinformatic programs which evaluate variants for splicing impact (MaxEntScan and Human Splice Finder) has been performed (Mucaki et al. 2013; Caminsky et al. 2014; see response to point 1). False positives are extremely rare because strength of binding sites (in bits) is directly related to their binding affinities; we have demonstrated that unused cryptic splice sites in the vicinity of natural splice sites are significantly weaker. Based on our experience and published analyses of a genome-wide site of binding sites (Rogan et al. 2003), such decoy splice sites are nearly always at least 4 bits (2^4 or 16 fold) weaker than sites that are actually recognized by spliceosomes.
 
The predictive accuracy of the IT methodology for detecting expression-validated mutations was determined to be 87.9% (762 of 867 variants from 122 different publications; changes to SRFs were included in this variant dataset; Caminsky et al. 2014). This value is similar to the predictions made to those of Cheung et al. (2019), where we predicted splice site and/or SRFs changes to 89.5% mutations experimentally validated to cause exon definition events.
 
The performance of IT-based methods for predicting splicing mutations has been well established over the past two and a half decades. Re-evaluation of its accuracy is not necessary, and this issue is, at best, only tangentially relevant, to the purpose of presenting the resource described in a Data Note article.
 
6. Of the >351k mutations with an effect on splicing supported by RNA data, only 35 affect CGC genes. Is this only a subset of mutations affecting driver genes or is it the complete list? In the former case, I would suggest that the authors provide the full list as supplementary data. In the latter case, the authors should discuss the implication of such a low number. Considering that there are >700 CGC genes, does it mean that aberrant splicing is very rarely a driver event? Is the overwhelming majority of splicing variants passenger?
 
Response: There are 25 CGC genes indicated in Table 2, however these were never intended to be interpreted to be a complete list of CGC genes with Veridical-flagged mutations. The table has been renamed to indicate that these are a set of representative mutations. In the previous version of this paper, the number of variants (“n=25”) did not indicate the total number of CGC genes. This has been removed and replaced with the actual number of CGC splicing mutations that were validated:
 
     “In Table 2, we highlight a subset of validated splicing mutations which were identified in known driver genes implicated in the COSMIC (Catalogue Of Somatic Mutations In Cancer) Cancer Gene Census catalog (CGC) ²⁷. In total, 543 “Tier 1” CGC genes have at least one Veridical-flagged variant present in the ValidSpliceMut database.”