Overexpressed somatic alleles are enriched in functional elements in Breast Cancer

Asymmetric allele content in the transcriptome can be indicative of functional and selective features of the underlying genetic variants. Yet, imbalanced alleles, especially from diploid genome regions, are poorly explored in cancer. Here we systematically quantify and integrate the variant allele fraction from corresponding RNA and DNA sequence data from patients with breast cancer acquired through The Cancer Genome Atlas (TCGA). We test for correlation between allele prevalence and functionality in known cancer-implicated genes from the Cancer Gene Census (CGC). We document significant allele-preferential expression of functional variants in CGC genes and across the entire dataset. Notably, we find frequent allele-specific overexpression of variants in tumor-suppressor genes. We also report a list of over-expressed variants from non-CGC genes. Overall, our analysis presents an integrated set of features of somatic allele expression and points to the vast information content of the asymmetric alleles in the cancer transcriptome.

Expression imbalance of point mutations is particularly informative for regions with no CNAs, where potential effects on the transcription can be directly linked to the underlying nucleotide change 14 . Therefore, quantitative integration of allele signals between same-source DNA and RNA is instrumental for tracking chromosome-of-origin effects. The latter, in turn, can be used to search for new genes whose allele behavior follows the pattern of known cancer drivers and is thus indicative for potential carcinogenicity. Therefore, the few studies that quantitatively integrate allele abundance from matching DNA and RNA sequencing sources are very informative 10 .
Herein, we apply a software that we recently developed -RNA2DNAlign 9 -to systematically quantify the allele expression of somatic variants in breast cancer samples from The Cancer Genome Atlas (TCGA). RNA2DNAlign counts variant and reference sequencing reads derived from compatible RNA and DNA datasets, and tests for allelic imbalance; it also calls positions with extreme allele distributions, including somatic over-Expression (SOM-E) or loss (SOM-L). We compute and compare the somatic variant allele fraction (VAF) of mutations in genes from the Cancer Gene Census (CGC) 17 to those in the rest of the genes in our samples. We also report a list of non-CGC genes with over-expressed somatic variants. Overall, we present an integrated set of somatic allele-specific expression features, in the context of their potential underlying functionality.

Strategy.
Our strategy was to first systematically quantify the variant allele fraction of the tumor RNA (VAF{tRNA}), and then to assess for correlation between RNA allele asymmetry and functional features (Fig. 1). Somatic variants (SOM) with a bi-allelic signal in the tumor DNA and a mono-allelic signal in the tumor RNA were classified as SOM-L (VAF{tRNA} ~ 0) or SOM-E (VAF{tRNA} ~ 1; Fig. 2). We assess both absolute VAF{tRNA}, and relative to VAF{tDNA}, for which we introduce the expression V R:D = VAF{tRNA}:VAF{tDNA}. We note that through accounting for the VAF{tDNA}, V R:D reflects the overall genome composition of the sample, including the contribution from large rearrangements, and admixture with non-tumor genomes (i.e. the sample purity). First, we analyzed the allele distribution for mutations in known oncogenes and tumor suppressors from CGC. We evaluated VAF{tRNA} and V R:D for correlation with functional features including conservation, predicted pathogenicity, and location in critical sequence motifs. Next, we assessed these features, in the context of their allelic expression, in the non-CGC dataset, and highlighted variants whose somatic allele patterns follow functionality-associated allele behavior of known cancer drivers.
Overall dataset characteristics. A total of 1238 (1139 unique) mutations in 921 genes, from which 68 were listed in CGC, satisfied the requirements for our analysis (Supplementary Table 1 and Supplementary Figure 1). Between 7 and 51 somatic point mutations in expressed coding regions were assessed per individual sample. Most of the mutations (94%) were singletons (present in only one sample), whereas 44 mutations were seen in 2, 12 in 3, 4 in 4, 2 in 5, and one mutation each was found in 6 and 7 different samples. Notably, all non-singleton mutations shared similar allele expression status across the different samples. A total of 437 somatic mutations (38.3%) were not expressed at all in the transcriptome (SOM-L), and 73 mutations (4.9%) were over-expressed (SOM-E). The analysis of the variant allele fraction showed an overall positive correlation between VAF{tDNA} and VAF{tRNA} (Spearman correlation r = 0.38, Fig. 3A-C). The functional distribution of the predicted consequences on the protein, and the intersection with their allele-expression status is presented on Fig. 3D. The missense, non-coding and stop-codon variants showed clearly different patterns of V R:D with a higher V R:D in the missense mutations, as compared to the non-coding and stop-codon variants (p = 0.0004, Kruskal-Wallis test 18 , Fig. 3E). Notably, we observed distribution towards higher V R:D of the variants predicted to be pathogenic through FATHMM (Functional Analysis Through Hidden Markov Models), Fig. 3F 19, 20 . CGC genes somatic allele expression: overall features. The 68 known cancer driver genes collectively contained 103 (88 unique) somatic mutations qualifying for the analysis (Supplementary Table 2) 17 . Mutations in PIK3CA, MITF, ACVR2A, CLIP1, and TCEA1 were called in more than one sample. In this gene-set, we called 10 SOM-E variants: seven missense substitutions, two synonymous variants, and, notably, the stop-codon R63X Figure 1. Major steps of the analysis of allele distribution for somatic variants in our dataset. V R:D was analyzed for correlation with different functional mutations groups in oncogenes, tumor suppressors, and the rest of the genes. SOM-E and SOM-L variants were compared with the rest of the somatic mutations for predicted pathogenicity and location in functional motifs such as transcription and splice factor binding sites, and highly preserved sequences.
in CDH1. Of note, four of the SOM-E missense substitutions were called in TP53 (See Supplementary Table 2). A higher number -25 -SOM-L variants were completely absent from the transcriptome in the CGC dataset.
Several noticeable observations were made in the CGC subset. First, different V R:D distribution was observed in the CGC variants as compared to the rest of the dataset (p = 0.02, Kruskal-Wallis test 18 , Fig. 4A); the difference due to larger proportion of CGC variants with higher allele expression. Second, the CGC missense mutations showed higher allele expression as compared to the missense mutations in the entire dataset (p = 0.03, Kruskal-Wallis test 18 , Fig. 4B). Notably, a tendency for higher V R:D was also seen for the stop-codon mutations, albeit not reaching statistical significance (Fig. 4C). In contrast, the non-coding variants did not show significant differences between the CGC and non-CGC genes (Fig. 4D). Third, we documented positive correlation between V R:D and predicted pathogenicity assessed by the CADD score (Combined Annotation Dependent Depletion) 21 , (Spearman r = 0.25), FATHMM score (Functional Analysis Through Hidden Markov Models) 19,20 (Spearman r = 0.17), and conservation of the position of the somatic mutation as assessed through GERP (Genomic Evolutionary Rate Profiling, Spearman r = 0.29) [22][23][24][25][26] . Of note, 21% of the variants in the CGC dataset modeled through FATHMM as pathogenic have been reported in cancer-based studies 17 . Collectively, all the above analyses supported preferential expression of functional alleles in the CGC dataset. We then assessed CGC SOM-E and SOM-L mutations in the context of their harboring gene's function and mechanism of action. The first noticeable observation was a tendency for over-representation of genes acting in recessive molecular mode among the SOM-E variants, as opposed to more-frequent dominant mode of action in the genes bearing SOM-L variants (p = 0.15). Recessive mode is traditionally more often associated with tumor-suppressive function, while dominant action is reported frequently for oncogenes 27 . In our study, SOM-E status appears not to result from a genomic DNA loss, as evident by the tumor DNA's biallelic signal (0 < VAF{tDNA} < 1). Both the inhibition of the reference and the enhancement of the mutant allele's expression could result in mutant RNA dominance, and these effects could be independent or related to the functionality of the particular mutation. In the case of the mutations acknowledged as pathogenic in suppressor genes, the observed overexpression is consistent with mutation-driven allele inactivation, possibly favored by positive selection forces. Such interpretation is in line also with the over-expressed stop-codon R63X in CDH1 28 .
For the SOM-L mutations, whether their expressional loss is linked to potential oncogenic action of the host gene, is to be determined on per-gene basis. It is important to recognize that many somatic variants are randomly lost in the tumor transcriptome, and the number of transcribed ones can depend on factors such as Estrogen Receptor (ER) expression levels 1 . While it is possible for a SOM-L variant to reside on a lost allele by coincidence, this is unlikely to explain all SOM-L patterns for variants with known pathogenicity.

Allele expression of somatic mutations in the non-CGC genes. The integrated features of somatic
allele expression in the non-CGC genes is presented in Supplementary Table 3. We documented concurrent to the CGC dataset positive correlation between increased allele expression and predicted pathogenicity and conservation scores (Spearman CADD r = 0.11, FATHMM r = 0.12, and GERP r = 0.17 (Supplementary Table 3).
The non-CGC somatic mutations with strong overexpression of the mutation-bearing allele (VAF{tRNA} = 1) are presented in Table 1. We next assessed the SOM-E variants for location within transcription and splicing factor binding sites, including analysis for generation of a new binding site outside of known protein -recognizable sequences 29 . Indeed, 18 out of the 42 non-CGC SOM-E variants positioned outside an existing TFBS were predicted to generate a new motif recognizable by either transcription or a splicing factor 29,30 . Next, we reviewed, on a per-gene basis, the current knowledge on the SOM-E genes and their possible implications in cancer. Despite not being listed in the CGC, some of these genes -such as MSH3 and NUAK1 and NFE2 -have been repeatedly linked to cancer [31][32][33] . Notably, more of the SOM-E genes linked to tumor suppressor features (as opposed to oncogenic, p = 3.8e-4, Metacore), which we concurrently observed in the CGC dataset 34,35 . Another striking observation is that 6 of the genes with SOM-E variants -MSH3, RAD51, TCOF1, TP53BP1, CCNB2, and TOP3B -are directly implicated in DNA damage response and repair [36][37][38][39] which was also the top-enriched pathway in the SOM-E dataset (p = 0.05, Metacore). In contrast, the most represented pathway in the SOM-L group was the immune response (p = 0.05, Metacore). In regards to GO annotations, two differences were detected between the SOM-E and SOM-L groups (Supplementary Figure 2). First, SOM-E variants were more frequently located in genes encoding receptors and signal transducers, while a higher proportion of the SOM-L variants resided in structure-supportive genes. In regards to biological processes, the SOM-E group was enriched in genes involved in response to stimuli.

Discussion
Ultimately, the accurate assessment of the expressed allele fraction is only possible in the context of the corresponding DNA alleles' content. Herein, we integrate matching RNA and DNA allele fraction from bi-allelic DNA regions to identify transcriptome-favored alleles. We focus more specifically on somatic point mutations in breast cancer, which we assess for tumorigenic functionality that can underlie selective transcriptome preference.
The first striking observation from our study is that transcriptome-preferred alleles are enriched in functional features, which are often predicted to alter the original protein function. This correlation was stronger in the group of genes traditionally acknowledged as tumor suppressors. Tumor suppressors are often lost during progression, and their loss is considered a contribution to tumor growth 40 . In our data we see a strong expression preference towards somatically mutated tumor suppressor transcripts, including such bearing a premature stop-codon. Increased allele expression can be either directly caused by mutation-promoted cis transcription  Continued activation, or/and retention of the mutant allele in the transcriptome via positive selection. Both scenarios infer functionality and growth-supportive potential. Conforming with that, highly expressed somatic variants, including SOM-E, were more frequently located in highly conserved and predicted to be functional genomic sequences. Taken together, these data are consistent with gain-of-function mechanism favored by the tumor transcriptome. An active role of over-expressed variants is also supported by the selection for maintaining the expression of a complete, translation-ready transcripts, suggesting a possible role of the altered/shortened proteins in the tumor progression. Indeed, once recognized as tumor suppressors, many of the genes in our SOM-E set, including TP53 are now acknowledged to play more complex roles that include oncogenic action [41][42][43][44] . Both inactivation and altering the protein function can be crucial for the tumor development. Regardless the mechanism of action, the above observations mark allelic overexpression as a highly informative metric that can be used to outline functionally enriched somatic datasets. The proportion of SOM-L alleles in our data is generally consistent with other reports 1 . Under-expressed alleles, including SOM-L, also correlated with functional annotations and regulatory motifs, though did not reach the significance of SOM-E. In contrast to SOM-E, SOM-L variants confer features that imply intolerance of the transcriptional machinery to the harbored variant. In the absence of CNAs, several mechanisms could potentially lower allele expression levels of mutation bearing transcripts. A well acknowledged scenario is the surveillance-driven targeting of transcripts with deleterious variants, the most prominent example of which is NMD 1,45 . A degradation mechanism can also take place where the mutation results in an unstable RNA structure 46 . Finally, a mutation can destroy a binding site for a transcription or splicing factor, thus directly abolishing the expression of the underlying alleles 14 . Additional factors, such as high ER expression levels, are also reported to correlate with a decreased number of expressed somatic mutations 1 . Besides the above mutation-focused mechanisms, SOM-L may result from random under-expression in the tumor transcriptome, and the general infidelity of cancer transcriptional machinery 47,48 . The later confers higher contribution of randomness towards SOM-L loci, which is likely to dilute functional annotations in this group.
Another striking observation from our analysis is the expression pattern of stop-codon mutations. Several recent studies have published decreased expression of stop-codon bearing variants in cancer, and have linked it to NMD 1,49 . Notably, in our data we see stop-codon bearing alleles over-represented as compared to the reference. Whether these expressed RNAs are translated into shorter proteins is subject of further studies, but this possibility is consistent by the presence of premature stop containing, translation-ready transcripts 1 . While NMD is knowledgeably impaired in cancer, our data suggests gene-selective NMD actions [50][51][52] .
Distinguishing pathogenic mutations from the more prevalent neutral variants constitutes one of the greatest challenges of cancer biology, leading to substantial effort towards developing confident analytic strategies. Modern methods integrate traditional frequency based approaches with expression abundance, functional effects, interaction networks and pathway context 13, 53-60 . Here, we integrate somatic allele fraction with most of the above strategies and the knowledge on tumor driving mechanisms, and evaluate the potential of asymmetric allele expression to predict cancer implicated variants. We document distinct allele signatures of cancer drivers at several levels. First, mutations in known cancer genes from our dataset presented more frequently with extreme allele patterns. An example is TP53, mutation in which were frequently either over-expressed or lost. Second, mutations in known cancer-implicated genes presented with higher allele expression. This was also reflected in the higher percentage of SOM-E variants among the known cancer genes. Third, SOM-E mutation sites were enriched in conservation and functional motifs. Cumulatively, these findings highlight the SOM-E status as a potential indicator for cancer-driving functionality. Based on the above, we list the non-CGC genes whose expression status follows the drivers-enriched SOM-E status (see Table 1); albeit not included in the CGC list, some of these genes have been linked to cancer before and are worth further investigation. In summary, our research illustrates an important correlation between asymmetric alleles and cancer-implicated functionality, and functionality in general, and underscores the vast information content of our strategy to systematically outline asymmetrically expressed alleles. This strategy is applicable to all types of cancer and is now enabled by the growing accessibility of matched DNA and RNA sequence data new tools for their integration and analysis 9, 61, 62 .  [65][66][67][68] . From these, we excluded samples with extensive (more than 3 standard deviations) number of somatic mutations, possibly due to clustered genomic rearrangements 69,70 . The remaining 72 samples (Supplementary Table 4) were retained for further analysis. We reviewed the pathology reports and retrieved the available clinical information; data for 41 (57%) of the studied samples was available (See Supplementary Table 4). The highest proportion of the samples were ductal adenocarcinomas, either ER, or ER/PR positive. We did not observe any significantly distinguishing somatic expression patterns, which is likely due to the small sample size. The purity, as assessed by the above-mentioned algorithms, is shown in Supplementary  Functional and enrichment analyses. Functional annotations, conservation scores and modeled pathogenicity were extracted using the SeattleSeq annotation 147 (http://snp.gs.washington.edu/ SeattleSeqAnnotation147/index.jsp). Pathogenicity was modeled using PolyPhen, CADD and FATHMM models, and Conservation was assessed based on Phast, GREP and Grantham Scores 20-26 . Gene Ontology categories, pathway enrichment and network analysis were assessed using Metacore (Claritive Analytics). Transcription factor binding cites were analyzed using TRANSFAC 7.0 29 and splicing motifs were assessed using SpliceAid2 30 .
Statistics. SOM, SOM-E and SOM-L variants were called based on a binomial test for variant and reference sequencing read distribution, as previously described 9 . The distributions of SOM-E and SOM-L across tumor-suppressors, oncogenes, and the rest of the genes in the datasets, as well as the distribution of functional elements across SOM, SOM-E and SOM-L, were assessed using the Fisher exact test, Pearson chi-square test, Kruskal-Wallis rank sum test, linear regression analysis, and the Spearman rank correlation coefficient 18, 74, 75 .
Yates's correction for continuity was applied for tests with less than 5 measurements in any category 76 . The means of the VAF across different mutation types were compared using Student's t-test 77 . P-values below 0.05 were considered significant. For multiple trials, the significance value was corrected using Benjamini-Hochberg False Discovery Rate (FDR) technique.