A practical framework RNMF for exploring the association between mutational signatures and genes using gene cumulative contribution abundance

Abstract Background Mutational signatures are somatic mutation patterns enriching operational mutational processes, which can provide abundant information about the mechanism of cancer. However, understanding of the pathogenic biological processes is still limited, such as the association between mutational signatures and genes. Methods We developed a simple and practical R package called RNMF (https://github.com/zhenzhang‐li/RNMF) for mutational signature analysis, including a key model of cumulative contribution abundance (CCA), which was designed to highlight the association between mutational signatures and genes and then applying it to a meta‐analysis of 1073 individuals with esophageal squamous cell carcinoma (ESCC). Results We revealed a number of known and previously undescribed SBS or ID signatures, and we found that APOBEC signatures (SBS2* and SBS13*) were closely associated with PIK3CA mutation, especially the E545k mutation. Furthermore, we found that age signature is closely related to the frequent mutation of TP53, of which R342* is highlighted due to strongly linked to age signature. In addition, the CCA matrix image data of genes in the signatures New, SBS3*, and SBS17b* were helpful for the preliminary evaluation of shortened survival outcome. These results can be extended to estimate the distribution of mutations or features, and study the potential impact of clinical factors. Conclusions In a word, RNMF can successfully achieve the correlation analysis of mutational signatures and genes, proving a strong theoretical basis for the study of mutational processes during tumor development.


| INTRODUCTION
The cancer genomes contain many mutations, which are derived from exogenous and endogenous mutational processes that operate during the cell lineage. 1 These mutational processes are cumulative effects of DNA damage and repair processes, indicating unique patterns in tumorigenesis, namely mutational signatures. [2][3][4][5] Each mutational process from a tumor may involve some particular signatures that their biological combination processes could induce a large number of mutations. 6,7 Via studying the completeness of these mutations and identifying the digital genomic footprints that contribute to the mutation characteristics of tumors, we can not only reveal the potential mutational process information, understand the carcinogenic mechanism of tumor occurrence and development but also provide biomarkers for early diagnosis, accurate cancer stratification and clinical response prediction, and realize individual treatment strategies. [8][9][10][11][12][13] Analysis of mutational signatures may reveal previously unknown mutation mechanisms and mysterious environmental exposure, such as herbal supplements containing aristolochic acid. 13 However, understanding of the pathogenic biological processes is still limited. Therefore, to systematically describe the mutational process leading to cancer, it is necessary to decipher the mutational signatures from the somatic mutation catalog by using mathematical statistical methods, [14][15][16][17][18][19][20][21][22][23] the number of mutations that each feature in a single sample can be attributed to each feature, which annotate the probability of each mutation class in each tumor and the possibility of each feature producing.
Currently, the final referenced mutational signatures are archived in the catalog of cancer of somatic mutations in cancer (COSMIC, https://cancer.sanger.ac.uk/cosmi c/ signa tures). Most of them are common in various tumors, of which are specific to a certain type of tumor, of which are part of normal cell biology, and of which are related to abnormal exposure or tumor progression. [24][25][26][27][28][29][30][31] They may be attributed to known environmental exposure and mutation processes, such as tobacco smoke, ultraviolet radiation, the activity of the APOBEC series of cyclobutylaminases, and DNA mismatch repair defects or mutations in POLE. Besides, as known to us, the association between genes and mutational signatures was confirmed, and much focus were paid to study the role of hotspot mutations in the formation of mutational signatures, which provides a good research idea for the mechanism of tumorigenesis and development. However, at present, there are few tools systematically and in-depth mining the relationship between genes and mutational signatures, 32,33 which undoubtedly does not bring much convenience to the association between mutational signatures and genes.
Here, we first used the R language to design a simple and convenient package RNMF, which can directly start from the mutation data set to realize the correlation analysis of mutational signatures. Then, we pooled 1073 samples from Asian ESCC population, and then used RNMF to verify the practicability of this method framework. During the analysis, we highlighted interesting correlations through association analysis with driving mutations. Finally, deep learning method is used to explore the cumulative contribution abundance (CCA) matrix image data of gene under fixed features, and hierarchical learning of prognosis is done.

| Genomic data collection
All somatic mutations were initially collected from the supplementary data of six previous studies (Table S1) comprising 1073 esophageal squamous cell carcinoma (ESCC) cases, including 508 genome-wide data and 565 exon sequencing data.

| R package link and parsing description
In the analysis process of this software, somatic variants can be imported from a Variant Call Format (VCF) file or a Mutation Annotation Format (MAF) file. Then it relies on the Bioconductor library, such as BS.genome. Hsapiens.UCSC.hg19 or BS.genome.Hsapiens.UCSC. hg38, to acquire an information matrix of mutation types (single-base substitutions [SBS], double-base substitutions [DBS], and insertion and deletion mutation [ID]). Subsequently, the program extracts mutational signatures according to the generated data and finally obtains the result files. The R package is publicly available at https:// github.com/zhenz hang-li/RNMF. The detailed document file also provides some examples of commands usage. In addition, scripts for running the package will also be provided in the R package. Alexandrov,15 an NMF software, was developed for decomposing the characteristic map based on MATLAB language long before. Subsequently, there were some related programs developed based on R language, such as SomaticSignatures. 16 However, these two software only targeted SBS Signatures' analysis. Then SigProfilerExtractor, which was developed based on Python language, appeared, aiming to analyze the three characteristic maps such as SBS, DBS, and ID signatures. Quite fittingly, our designed RNMF software provides the ability to obtain the count matrix of SBS, DBS, and ID directly from mutation data sets by a function interface. To connect the result format obtained by this function more conveniently and select a more optimized initial value, we decided to build a new R interface based on NMF principle and generate the data result format required in downstream analysis. Therefore, based on the model definition in previous reports, 34,35 we can obtain a classical formula V = PS + E, which is used to extract the mutational signatures in human cancers. In this equation, V refers to the observation matrix with size M × N, of which M represents the observed characteristics, and N is the number of samples. Supposing the number of mutational signatures is K, then we can estimate a non-negative mutational signature matrix P with size M × K and a non-negative abundance fractions matrix S with size K × N. Simultaneously, error matrix E that refers to nonsystematic errors and sampling noise is calculated during the processing.

| Optimal mutational signature extraction framework
For an observed mutational catalogs V, these K mutational signatures could be extracted by denovoNMF ( Figure 1) as following: Step 1 (Capture random matrices): We randomly generate matrices P (P ≥ 0) and S (S ≥ 0). Ideally, for a mutational signature, its components are basically fixed, so the sum of its standardized components is equal to 1; for one sample, the total of normalized abundance fractions for each mutational signature should be infinitely close to 1. Hence, here we require that Step 2 (Optimize the initial value): We apply resampling to obtain a new matrix Ṽ Ṽ ≈ V by using the Dirichlet distribution, and straighten matrices P and S by columns and then merge them into a solution vector x. An optimized objective function F I G U R E 1 Overview of RNMF workflows. The automatic process can start from inputting mutation data files in MAF or VCF format to generate the necessary mutation type matrix, including single-base substitution mutation types with 96 different contexts, 78 strand agnostic double-base substitution mutation types, and a compilation of 83 different types of small insertion and deletion (ID) mutation types. Next, the process can use two program interfaces (denovoNMF and InverseNMF) for signature analysis to get the mutational profiles. Then, according to the clinical grouping, non-silent mutation data and the list of genes to be studied, combined with P matrix and S matrix, the cumulative contribution abundance matrix of genes was calculated. At the same time, the corresponding image data are generated according to the specified order, and finally the training and analysis are carried out by machine learning or deep learning methods, so as to obtain the expected results min P,S E(P, S, , which previously reported 34 is used to find the best solution, where α and β are the parameters to control the accelerated convergence. After smoothing x, the updated matrices P and S are restored according to the straightening rules. Finally, we generate a new initial mutational signature matrix P and a new initial abundance fraction matrix S through 20,000 iterations via referring to the previous implementation process. 15 Step 3 (Re-updated the initial value): Perform Steps 1 and 2 for I(I ≥ 5) iterations. Their errors generated by iterations are calculated by formula E = ‖V − PS‖ 2 , and the results of the 5 items with the smallest error are selected. Then we apply the k-means 36 algorithm to the set of matrices P and S to cluster the data into K clusters, respectively. Subsequently, class-center P and S are obtained. Similar to Step 2, we straighten them to calculate the optimal solution space and finally gain the excellent initial value.
Step 4 (Rerun NMF): In this step, we still use the multiplicative update rules to generate the final matrix. The iterative model is as follows: where γ is a parameter to control the accelerated convergence. Iterate until P and S convergence or until the maximum number of 100,000 iterations is reached.
Step 5 (Build a solution space): Repeat the process of steps 1 to 5 with Î Î ≥ 20 times to generate a solution space for all the value of K K ∈ N + . Then suitable selection in this solution space using the silhouette coefficient measure and error gradient. The criteria we chose are as follows: 1) stable, without sudden decline or relatively large gradient of descent and large width of confidence interval and 2) small standard error term and the gradient of standard error between adjacent classes tends to be gentle.
After a series of analysis process, the final mutational profiles of each K clusters can be acquired. Follow the previous experience, 20 we use the cosine similarity to determine the similarity between two mutational signature A and B.

| Deriving the contribution of defined mutational signatures
At the same time, we also developed a reversible method named InverseNMF, which can identify mutational signatures within a small data set or a single tumor sample. Previous reports 20,33 have confirmed the importance of such applications and provided another analytical strategy for tumor characteristic map studies.
In this step, the P matrix is a user-defined matrix (such as the signature matrix provided by COSMIC), and then a feasible sample contribution matrix S can be estimated iteratively by observing matrix V. The iterative model is as follows: where γ is a parameter to control the accelerated convergence. Iterate until S convergence or until the maximum number of 100,000,000 iterations is reached.

| Mutational signature operative in ESCC
We applied our framework to extract mutational signatures from 508 WGS samples of 1073 ESCC tumors. At the same time, the previously reported tool "SigProfilerExtractor" 23 was used to analyze the data. Finally, the similarity between the two results was analyzed ( Figure S1). Our framework optimizes the non negative matrix factorization (NMF) algorithm and takes the change of silhouette coefficient and error gradient as the evaluation index of feature number selection ( Figure S2). Here, the threshold is set to 0.9.
To analyze and explore the potential features of exons of 1073 ESCC samples, We used COSMIC Mutational Signatures (v3.2 -March 2021 and v2.0 -March 2015) and mutational signals extracted from 508 WGS data as the background to obtain the number of mutations in each mutational signatures. Based on the results of three single base mutation data sets (508-WGS cohort renames as WGS508, exon region of 508-WGS cohort renames as EXON508, and exon region of 1073 samples renames as EXON1073), the eigenvalues of mutational signatures for subsequent analysis of ESCC were determined.

| Cumulative contribution abundance of genes
As described earlier, 7 because each mutational process takes into account the mutation category and the generation of mutations in tumor is attributed to each corresponding process, we define that the influence of mutations of category m in tumor n during the mutational process can be expressed as: where the unit P mk S kn represents the influence of mutations of category m attributed to signature k in tumor n. Then the probability effects in mutations of category m due to signature k in tumor n can be expressed as: Let us now consider that all genes with the number of mutations of each mutation type in tumor n can be calculated in matrix that designed as following: where the G represents the number of gene lists in the data set. Here, we define n mg as the impact factor of mutations of category m in sample n due to a gene g, and it can be calculated by the following formula: Then the contribution of mutational signature s to mutations of category m in a gene g in tumor n can be estimated as: Therefore, we consider that the cumulative abundance of a gene g that is attributed to mutational signature s in tumor n can be estimated as follows: On the other hand, considering the influence of gene length, we also define the relative abundance of a gene g that is attributed to mutational signature s in tumor n, which can be calculated as follows: where l g represents the length of the gene g. And (sg, mn) is on behalf of the relative abundance of mutational signature s to mutations of category m in a gene g in tumor n can be estimated as: where the ′n mg can be calculated by the following formula: Finally, based on the research results of the previous workers, 7,14,21,37 we integrated and optimized some methods to evaluate the association between gene mutation and mutational signatures. We hope that through these methods, we can provide a convenient program interface for researchers, and further provide new ideas for the study of tumor mechanism.

| Prognostic analysis
Kaplan-Meier survival analysis and Cox proportional hazards model were used to analyze an association between cancer-related genes and prognosis. Kaplan-Meier survival and Cox regression analyses were carried out with the R survival package (2.40-1). p-value less than 0.05 was considered to be statistically significant.

| Image analysis
To obtain some favorable statistical information, and even get some prognosis evaluation or medication guidance, which is helpful for clinical treatment, we designed a way to transform CCA of gene in each mutational signature into intuitive image information. The detailed operations are as follows: (1) we can sequence these genes through some potential relationships, such as similarity, pathway, or clinical association; (2) according to these CCA matrix results, some standard visual impact images are generated; and (3) we can apply these images to deep learning model, combined with clinical information for analysis and mining.
In this work, we use "hclust" clustering to get the gene sequence based on the CCA matrix of genes. From our mutation data set, we will screen out the list of all non-silent mutation genes in the existing mutation set and do intersection processing with 717 genes. The intersecting gene set is the final data information transformed into image. We organize the data set into a N × N matrix (where N = √ Number of intersecting genes). In this paper, we arrange the genes from top to bottom with the aim of obtaining the effect picture of gene mutation. To study the effect of CCA changes in each mutational signature on the prognosis, the deep learning framework PyTorch 38 was performed to train some publicly available models, such as ResNet 39 (ResNet50, ResNet101), DenseNet 40,41 (DenseNet-121, DenseNet-161), Inception-V4 42 and MoblieNet-V2, 43 V3. 44 First of all, we uniformly limit the pixels of the picture to 500*500.

| Framework analysis of RNMF
Optimizing and improving the derivation of mutational signatures can not only rediscover some known features, but also produce new discoveries that were previously masked by technical and biological confusions. Here, we designed a simple and convenient process framework by R language (Figure 1) and developed a novel R package called RNMF (https://github.com/zhenz hang-li/RNMF). The package RNMF can directly analyze the mutation data set which exist in the form of MAF or VCF format files, rapidly obtain the number of specific mutation types in each sample and then resolve the mutational signatures. This software can achieve 7 kinds of outstanding functions: (1) the mutation rate per M shown in LEGO graph; (2) rapidly extraction of mutational signatures; (3) accurately evaluating the contribution of each sample to the known mutation profile; (4) deeply studying on the similarity of mutational signatures and displaying the heat-map; (5) calculating the CCA matrix of genes; (6) studying the relationship between mutational signatures and genes; and (7) studying on the relationship between hotspot mutation of driving factors and mutational signatures.
One outstanding function of RNMF is defined as SigsInput, which is suitable for large or most small data sets, providing analysis interfaces for different types of genomic DNA changes, such as SBS, DBS, small ID. To prove high performance of our designed RNMF, the mutation data sets of 1073 ESCC samples from Asia (including 508 WGS data) are analyzed thoroughly. These mutation data were extracted from the appendix of published articles (Table S1), and were annotated by Oncotator. 45 As shown in Figure 1, by extracting the mutational signatures with the designed denovoNMF function, the evaluation results exhibits the change of silhouette coefficient and error gradient.
In the Figure 1, an application interface "InverseNMF" can calculate the fraction of one sample which is contributed to each provided mutational signatures. By our designed RNMF, we analyzed the contribution scores of these coding region mutations to each mutational signature which are gotten from the COSMIC Mutational Signatures (v2.0 -March 2015). The results are highly consistent with the analysis ones of the known "MutationalPatterns" 20 and "deconstructSigs" software, 37 as shown in Figure S3A. It is worth noting that operating speed of "InverseNMF" rivals to that of "MutationalPatterns" software, however, is 12 times faster than that of "deconstructSigs" software ( Figure S3B). To go insight into the inner link between mutational signatures and gene mutations, a series of powerful functions are designed according to the the CCA calculation model of gene, such as "cumulativeCA," "ge-nePerMutSigs," "samFisherSigs," "samPerMutSigs," and "eachMutationCA." The detailed analysis description will be provided in the example of our designed RNMF.
According to the previous reports, 7,24,33 we define that if the CCA of a gene in a sample to a signature is less than 6%, the gene in the sample has little effect on the signature. Next, we combined non-silent mutations, sorted gene list and clinical grouping information to obtain gene mutation image sets of samples under different mutational signatures based on gene CCA matrix. Then, deep learning method was used to process these image information for hierarchical learning, which could obtain some biological cognition, so that we could deeply understand some useful information for adjuvant therapy, such as prognosis evaluation, efficacy evaluation, medication guidance and so on (Figure 1).
In a word, as a versatile R package, RNMF can realize parallel operation that helps to study and evaluate the mutational processes during tumor development. Thus, molecular analysis can be performed based on extracted mutational signatures, further revealing the molecular mechanisms and optimizing the diagnosis and treatment decisions.

| Identifying mutational signatures via RNMF
The incidence and mortality of esophageal cancer have always been relatively high, as in China, for example, among which ESCC accounts for 90% of esophageal cancer. 46 A previous study has provided a large genome-wide sequencing cohort of Chinese ESCC population. 47 To demonstrate the role of RNMF in the extraction of mutational singatures, our designed RNMF is used to systematically analyze the mutation data set (single base mutations and INDEL mutations) of this cohort data. In the overall single base mutation pattern of ESCC, it is mainly C>T and C>G mutation, followed by C>A mutation, accounting for 34.9%, 18.81% and 15.94%, respectively ( Figure S4). Besides the insertion of a zero-length 1-bp T base homologous sequence, most other types of deficient incongruities are characterized by long (≥5) thymine mononucleotide repeats.
To verify the accuracy of the above results, the known SigProfilerExtractor software 23 is also used to analyze SBS signatures, and extract 13 features ( Figures S2B, S5). By further making a similarity comparison, we found that whether using V2 or V3 version of COSMIC signature, it shows strong consistency ( Figure 2, Figure S1A,B).
Moveover, the results were basically consistent between RNMF and SigProfilerExtractor ( Figure S1C). Thus, these results enough show the validity of the RNMF. In addition, we also extract ID signature. Nine ID signatures were prominent ( Figure S2C, Figure 3, Table S2C-D), and three of them have been previously reported, including two with known mutational processes. 23 Mutational Signatures description are represented in Table 1.

| Estimation of sample contribution under known mutational signatures by RNMF
To verify the practicability of the functions in our framework, the contribution of samples was analyzed in the known mutational signatures. We collected a total of 115,130 somatic mutations from exon region of 1073 ESCC samples, including 16 of them were hyper-mutated with mutation count more than 500 ( Figure 4A, Table S3). Comparing with 508-WGS cohort, we found that the overall proportion of C>T mutations increased in the exon region, among which *[C>T]G context changed greatly ( Figure S6A). However, the single base mutations in the exon region were mainly C>T (48.19%) and C>G (18.2%), followed by C>A (13.8%), indicating that the mutation pattern of exon region was similar to that of whole genome.
To better explore and analyze the potential features of exon regions, we used Single Base Substitution (SBS) Signatures (v3.2 -March 2021 and v2.0 -March 2015) from COSMIC database as the provided mutational signatures to obtain the number of contribution mutations of samples in each signature. We analyzed three single base mutation data sets: 508-WGS cohort (WGS508), exon region of 508-WGS cohort (EXON508) and exon region of 1072 samples (EXON1073). For the 12 feature maps extracting from WGS data, the trends of WGS508, EXON508 and EXON1073 are basically the same under the same background ( Figure 4B, Figure S6B). However, the results of WGS508 revealed that the signature, which accounted for a large proportion, may not exist in the list of 12 decomposed signatures, and these signatures either have high similarity or low similarity related to those 12 signatures. What's more, the results of WGS508 showed that compared with the results of different versions of COSMIC Mutational Signatures, we found that with the increase of features background analysis, high similar or low similarity signatures with high proportion will be produced, indicating that adding background features may introduce some unimportant or similar features to share the load weight ( Figure 4B, Figure S6B). Meanwhile, under the same background, the results of EXON508 and EXON1073 showed that the number of samples play little role in the proportion of signatures ( Figure 4B, Figure S6B). In addition, previous reports 24,47 showed that each kind of cancer has its own important characteristics. Of which the number or type of signatures is usually different. These different main signatures play a leading role in the occurrence and development of different tumor types.
In this paper, by combining the results of our analysis and the reported results ( Figure S5), we strongly believe that 12 stable SBS signatures, which is extracted from the 508-WGS cohort, should be the outstanding features of ESCC, which are considered as the leading mutation patterns in the occurrence and development of this tumor. Therefore, these 12 features were selected as background walls for the analysis of ESCC tumors (Table S4). We found that the signatures SBS1*, SBS2*, SBS5*, SBS13* and SBS15* were effective in most samples containing more somatic mutations, which are considered as ubiquitous signatures Note: "-" represents there is no match feature or the similarity is generally very low.

| Cumulative contribution abundance analysis reveals potential prognostic features and mechanisms through RNMF
In this work, a total of 717 cancer-related genes are screened from the COSMIC census (https://cancer. sanger.ac.uk/census, Table S5) and their CCAs were calculated via RNMF software with a given model function named "cumulativeCA". Through this program algorithm interface, CCA of a gene on a signature and CCA of a gene on a signature in a sample can be obtained. We found that these cancer-related genes are more enriched in ubiquitous signatures, especially signatures such as SBS1*, SBS2*, and SBS13* ( Figure 4E). There are differences in CCA levels of genes under different mutational signatures, which indicates that genes have their preference for mutational signatures. At the same time, survival analysis found that CCA of several genes was associated with prognosis ( Figure S7A, Table S6). In addition, a multivariate Cox model was confirmed that some of them were still significant ( Figure S7B), such as ARHGAP5, SETD2, RNF213, CDKN2A, NOTCH1, NFE2L2, and so on. Analysis of mutation characteristics showed that TP53 mutation significantly increased exposure to SBS1*, but conversely, the contribution to SBS1* was significantly reduced in samples carrying RNF213 mutation ( Figure 4F, Figure S8A). Furthermore, among the TP53 mutant samples, the contribution of the RNF213 mutant samples to SBS1* was significantly reduced( Figure S8A). Similarly, we found a significant increase in the contribution of APOBEC signatures (SBS2*, SBS13*) in ESCC samples with PIK3CA mutations, and the other three genes (ARHGAP5, SETD2, and UBR5) were also associated with APOBEC signature (SBS13*) ( Figure 4F, Figure S8A). It is noted that NFE2L2 mutation was related to SBS16*. These above results indicate that there is a potential mechanism between gene mutations and mutational signatures in the process of tumorigenesis and development. Thus, we studied 40 representative genes from 717 cancer-related genes, which contained about 3.6% of the total number of nonsilent SNV mutations. Most of these mutations preferred the characteristic SBS1* (15%), SBS2* (11%), SBS5* (12%), , each color represents a mutational signature, the length of each column represents the contribution ratio of mutation to the signature, the red mark represents the signature most similar to the 12 mutational signature, and the green arrow and green font indicate that this signature is very similar to the 12 mutational signature. Yellow means that the similarity between the signature and 12 mutational signatures is very low. (C) Percentage of ESCC tumors in each signature was displayed (top) and mutation rate for each signature in the relevant samples (bottom). If the contribution of a sample assigned to one signature is not less than 20%, we would consider that this signature is present in the sample. (D) The association between mutational activity of single-base substitution (SBS) signatures (SBS5* and SBS16*) and alcohol consumption. The violin compared the difference between drinking and non drinking groups from the mutation count and the contribution of samples to the signatures, and the significant p-value was statistically analyzed by Student's t-test with two-sided.
(E) The heatmap shows the distribution of cumulative contribution abundance of cancer-related genes from the COSMIC census, and the depth of color represents the degree of correlation. (F) Box plot showing that the SBS signatures were associated with cancer-related genes mutations (including SNV and indel), where n represents the number of samples. Statistical significance was tested by rank sum test with two sided. (G) The contribution of non-silent mutations in the coding regions of 40 cancer-related genes is statistically analyzed. Each color represents a class of mutation feature map, and the pie chart shows the proportion of each feature. In the figure, the top figure shows the total CCA of each gene in 1073 samples, while the bottom figure shows the proportion of each mutational signature, and a column represents a gene. Gene selection rules: the number of non-silent mutations is more than 30 and belongs to cancer-related genes from the COSMIC census SBS13* (10%), and SBS18* (12%), but were less distributed in SBS17b* (3%) ( Figure 4G). We found that 66.5% of 1073 ESCC samples had non-silent SNV of TP53, resulting in a highest level of CCA of TP53 gene, and the proportion of three SBS signatures (SBS1*, SBS5*, and SBS18*) was higher, followed by SBS15*. Obviously, different genes have different ratios for different mutated traits. It attracts our attention that PIK3CA gene is obviously in favor of APOBEC signatures (SBS2*), which accounts for more than 50% ( Figure 4G). Hence, all findings above suggest that preference of cancer-related genes for mutational signatures can be defined by CCA, which can further expose some potential prognostic features or mechanisms.

| Cumulative contribution abundance analysis exposing APOBEC are closely associated with PIK3CA helical mutation E545K
According to the reports, [48][49][50] PIK3CA is a typical protooncogene that typically harbors some hotspot mutations in tumors and is enriched in APOBEC characteristics in a variety of cancer types, especially these two most-common and well-studied hotspots: E542K (c.1624G>A) and E545K (c.1633G>A) in the helical domain. In ESCC, the results of previous studies implicitly implicated APOBEC activity as a key driver of PIK3CA mutagenesis. 50,51 In this cohort, we are committed to further study the potential mechanism between PIK3CA biological mutations and APOBEC signatures via CCA model. By the CCA enrichment analysis, it is found that tumors with non-silent mutations in PIK3CA had increased activity of the signature SBS2* (ESCC1073-EXON: 89 tumors with non-silent PIK3CA mutations and a median increase CCA of 0.74 per sample; q < 0.0001, p < 0.0001; regular tumors of ESCC1073-EXON: 88 tumors with nonsilent PIK3CA mutations and a median increase CCA of 0.74 per sample; q < 0.0001, p < 0.0001; WGS508-EXON: 38 tumors with non-silent PIK3CA mutations and a median increase CCA of 0.741 per sample; q < 0.0001, p < 0.0001; Figure 5A, Figure S8B). Analogously, those tumors with non-silent mutations in PIK3CA also had increased activity of the signature SBS13* (ESCC1073-EXON: 89 tumors with non-silent PIK3CA mutations and a median increase CCA of 0.136 per sample; q < 0.0001, p < 0.0001; regular tumors of ESCC1073-EXON: 88 tumors with non-silent PIK3CA mutations and a median increase CCA of 0.1359 per sample; q < 0.0001, p < 0.0001; WGS508-EXON: 38 tumors with non-silent PIK3CA mutations and a median increase CCA of 0.133 per sample; q = 0.0008, p < 0.0001; Figure 5b, Figure S8C). To further prove these connections, we performed mutational signature enrichment analyses and gained the same results ( Figure S8B,C), providing the strongest statistical evidence for the relationship between PIK3CA mutation and APOBEC signatures activity in ESCC. Together, these results further strongly suggest that, although APOBEC signatures activity are present in all tumors, somatic PIK3CA mutations are associated with a significant increase in APOBEC signatures activity.
We also dissect the mutation spectrum of PIK3CA non-silent mutations, revealing the underlying mechanism of during mutation processing. In our series, 65.22% of PIK3CA non-silent mutations were C>T substitution, and those mutations were frequent presenters that mostly contributed to APOBEC signatures (SBS2* and SBS13*) with highest percentage ( Figure 5C). Simultaneously, we investigated the PIK3CA helical (E545K: 36.96%; E542K: 14.13%) and kinase (H1047R: 11.96%) hotspot mutations, and found that only the helical mutations had a high cumulative proportion for APOBEC signatures (SBS2* and SBS13*) ( Figure 5C). Then, We related the PIK3CA helical mutations to each APOBEC signatures, and observed a significant increase for mutational exposure of APOBEC signatures in samples harboring helical domain mutations ( Figure 5D-E). Significantly, tumors carrying a hotspot mutation E545K significantly hold a high degree of contribution fraction of SBS2*, yet hotspot mutation E542K can not bring significant benefits to SBS2*, as well as the other mutations ( Figure 5D), implying that only mutation E545K can affect the benefit of the overall mutation data of the sample on SBS2* compared with other mutations. Similarly, we found that although the E542K mutation significantly increased the benefits of SBS13*, the significant intensity of the increase was not as high as that of E545K (Median: 0.146 vs. 0.172) ( Figure 5E), indicating that the E545K mutation in PIK3CA can accelerate the increase of SBS13* activity. Furthermore, from the perspective of gene itself, the CCA of PIK3CA genetic hotspot mutations for APOBEC signatures was significantly higher than that of other mutations ( Figure 5F-G). However, compared with E542K mutation, the effect of E545K mutation is more significant, indicating E545K among PIK3CA mutation is more closely associated with APOBEC signatures. In a word, PIK3CA helical mutation E545K contributes more significantly to APOBEC signatures, suggesting that they are strongly association with APOBEC.

| Cumulative contribution abundance analysis displaying the relationship between age signature and TP53 typical hotspot mutations
In previous study, the results of TP53 mutations on mutational signatures indicates that driver mutations of TP53 mutations are associated with specific mutation processes in human cancers, such as colon, skin, bladder, lung, and liver cancers. 51,52 They only mentioned that the most frequent TP53 mutations were associated with the most commonly observed age signature which featured by C>T transitions at CpG dinucleotides. It is worth considering that there is no detailed report on the association between TP53 typical hotspot mutations and age signature in the related studies of ESCC, including the previous large cohort analysis of ESCC. 47,53 Here, we analyzed the association between TP53 mutations and age signature confirming that age signature was associated with TP53 mutations via CCA enrichment analysis (Regular tumors of ESCC1073-EXON: 715 tumors with non-silent TP53 mutations and a median increase CCA of 0.0823 per sample; q < 0.0001, p < 0.0001; ESCC1073-EXON: 728 tumors with non-silent TP53 mutations and a median increase CCA of 0.085 per sample; q < 0.0001, p < 0.0001; Figure S8D). Moreover, mutational signature enrichment analyses also revealed the strong relationship between TP53 mutation and age signature activity ( Figure S8D). It is worth noting that TP53 was mainly enriched with C>T substitutions (44.51%) with a large proportion of them were preferentially contributed to age signatures (SBS1*) ( Figure 6A). We screened the six kinds of hotspot mutations with the highest risk rate (R342*: 4.07%; R213*: 3.37%; R282W: 2.95%; R175H: 2.81%; R273H: 2.53%; R248Q: 2.11%) and analyzed their association with mutational signatures. We found that except R175H, the other five hotspots preferred the age signatures (SBS1*) ( Figure 6A), which indicated that there was a potential mechanism between these hotspots and age signatures (SBS1*). Compared with other F I G U R E 5 Cumulative contribution abundance (CCA) enrichment analysis identifies an association between somatic PIK3CA mutations and activity of SBS2* or SNS13* in esophageal squamous cell carcinoma (ESCC). (A-B) Here, we use two data sets: exon regions of 1073 and 508 ESCC cases. First, the median CCA of each gene in the current signature is calculated, and then the contribution importance of each gene is calculated by PERMUTATION test to study the association between gene and signature. The regular patterns in the figure represent the samples with non hypermutated. For genes mutated in >5% of samples, the CCA of genes attributed to SBS2* or SBS13* was compared in tumors with wild-type versus mutated copies of the gene. Genes with FDR q < 0.1 are highlighted in red. (C) Mutation trend and hotspot analysis of PIK3CA non-silent mutations: the pie chart on the top left shows the proportion of six mutation types, and the pie chart on the top right shows the proportion of coding protein, with the name of the protein that accounts for a large proportion. Each color in the figure below represents a mutational signature, and the statistical proportion of the contribution of six mutation types to each mutation signature is on the left, a column represents a mutation type; The figure on the right shows the contribution abundance of classical hotspot mutant protein to each signature, and a column represents a classical mutant protein. (D-E) The violin diagram shows the difference of contribution to SBS2* and SBS13* between the samples with PIK3CA hotspot mutation and other types of samples, whereas (f-g) shows the difference of contribution abundance between the samples with PIK3CA non-silent mutations (E545K and E542K) and other samples. The Wilcoxon rank sum test with two-sided is used here, * represents (0.01 ≤ p < 0.05), * * (0.001 ≤ p < 0.01), ***represents (p < 0.001), and ns represents (p ≥ 0.05) mutations, tumors harboring at least one of these hotspots will significantly increase its contribution to age signature (SBS1*) ( Figure 6B). Actually, although R282W can improve the contribution of the sample to age signature (SBS1*) (R282W vs. other mutation vs. wild-type: median increase of 0.367 vs. 0.306 vs. 0.276), only hotspot mutation R342* can significantly affect the benefit of the whole mutation data of the sample to age signature (SBS1*) (R342* vs. other mutation vs. wild-type: median increase of 0.388 vs. 0.305 vs. 0.276; Figure 6B), indicating that R342* mutation is the primary factor to increase the activity of age signature (SBS1*). However, from the perspective of gene mutation itself, the CCA of TP53 typical hotspot mutations for age signatures was significantly higher than that of other mutations ( Figure 6C), which indicated that except for R342* and R282W mutations, the other four hot spot mutations could not significantly improve the overall mutation contribution level of the sample, but they are still involved in the specific mutational process. In addition, under age signature (SBS1*), tumors with CCA of TP53 no less than 0.06 were significantly associated with poor prognoses (p = 0.44, Figure S7A). Notably, CCA of TP53 in tumors carrying TP53 mutation R248Q were more than 0.06, and those cases were associated with deceased survival outcomes ( Figure S8E). Multivariate cox model shows that TP53 hotspot mutation R248Q are independent prognosticators for poor survival in ESCC ( Figure S8F). Finally, we also found that TP53 small INDEL mutations were related to New3 (ID signature; Figure S8G).

| Image analysis based on cumulative contribution abundance matrix of genes
Considering the limitation of relatively short follow-up time, we set three-year survival as a threshold and divided the sample into five groups (Table S7): G1) Survival status, follow-up time is less than 3 years; G2) survival status, follow-up time is not less than 3 years; G3) death status, follow-up time is not less than 3 years; G4) death status, follow-up time is less than 3 years; G5) others (rejection). Because G1 individuals are living and their follow-up time is less than 3 years, we will use G2-3-4 groups for the next analysis and mining, as shown in Table S7A. Limited by clinical information, and to better evaluate the prognosis, we especially compared G2 and G4 groups. The above grouped samples were respectively converted into images by CCA matrix ( Figure S9A) for deep learning model analysis. For choosing an ideal model, we randomly select a feature as a template for training, and take the high average accuracy as the judgment basis for model selection. Here we choose one SBS signature (New) as an example, as shown in Table S7B. Four of them have higher accuracy, and they are Resnet50 (68.571%), DenseNet121 F I G U R E 6 Cumulative contribution abundance enrichment analysis identifies an association between somatic TP53 mutations and activity of SBS1* in esophageal squamous cell carcinoma (ESCC). (A) Mutation trend and hotspot analysis of TP53 non-silent mutations: the pie chart on the top left shows the proportion of six mutation types, and the pie chart on the top right shows the proportion of coding protein, with the name of the protein that accounts for a large proportion. Each color in the figure below represents a mutation feature, and the statistical proportion of the contribution of six mutation types to each mutation signature is on the left, a column represents a mutation type; The figure on the right shows the contribution abundance of classical hotspot mutant protein to each signature, and a column represents a classical mutant protein. (B) The violin diagram shows the difference of contribution to SBS1* between the samples with TP53 hotspot mutation and other types of samples, while (C) shows the difference of contribution abundance between samples with TP53 hotspot mutations and other samples. The Wilcoxon rank sum test with two-sided is used here, * represents (0.01 ≤ p < 0.05), ** (0.001 ≤ p < 0.01), *** represents (p < 0.001), and ns represents (p ≥ 0.05) (71.429%), MoblieNetV2 (68.571%), and InceptionV4 (68.571%). To test the stability of our model, the above four models were trained for 10 times. In the training process, the model parameter random seed was fixed, other parameters are the same. Finally, through the analysis, we found that DenseNet-121 model has a higher average accuracy (Table S7B; Figure S9B). Consequently, we chose DenseNet-121 as a model to analyze all signatures, the full schematic representation as shown in Figure S9C. In the process, a stochastic gradient descent method 54 was used with an initial learning rate of 0.01, weight decay of 10 -4 and momentum of 0.7 in the process of training. Next, dropout, data augmentation and L2-regularization were applied to prevent overfitting. The above parameter sets were properly tuned for DenseNet-121 model. Then, for testing the stability of data and finding the global optimal solution, the model DenseNet-121 was trained 10 times for each sub-feature data of G2-G4. Random seed was set free in the training process. That some less accuracy than others may be a local optimal solution because stochastic gradient descent method was used as the optimizer in the training process. We found that the accuracy was comparatively stable to each sub-feature and relatively higher in mutational SBS3*, New and SBS17b* ( Figure S9D), suggesting that the beneficial feasibility of this conversion method of the CCA matrix image data of gene. Simultaneously, the results are given the best accuracy of G2-G4 is the sub-feature SBS17b* (77.500%), followed by SBS3* (72.500%) and New (71.429%). Furthermore, the probability distribution over the above three subfeatures of G3 group in the G2-G4 was tested. Interestingly, the distribution of G3 group is more likely to fall on G2 group ( Figure S9E). So further the model DenseNet-121 was also trained 10 times for each subfeature data of G2G3-G4 ( Figure S9D). The results show that the best accuracy of G2G3-G4 is the sub-feature SBS3* (76.190%), followed by New (72.973%), SBS16* (71.111%), and SBS17b* (70.732%), as shown in Table S7C. This illustrate that the G3 group addition has slight effect on the classification results. Additionally, the mutational signatures such as SBS3*, New and SBS17b* still have a high degree of explanation for G2-G4 and G2G3-G4 (Figure 7). Finally, we found that the survival group or the samples with a follow-up time of no less than 3 years had a higher contribution to SBS3* and SBS17b*, respectively ( Figure S9F).

| DISCUSSION
Here, we provide an integrated mutational signature analysis framework with a CCA model of genes, achieve a meta-analysis of 1073 ESCC samples, and verify the practicability and application value of our framework.
Via this framework, we obtained known and uncovered previously undescribed signatures (including 12 SBS signatures and 9 ID signatures) from 508 WGS tumors of 1073 ESCC cases. And further identified and highlighted an association between PIK3CA helical mutation E545K and activity of APOBEC signatures. Similarly, we also reported that age signature and the hotspot mutation R342* of TP53, and TP53 (R248Q) is a poor predictor for ESCC. In addition, the CCA matrix image data of genes under mutational signatures New, SBS3*, and SBS17b* were calculated. This is helpful for the preliminary evaluation of short-term prognosis.
In addition to feature extraction and sample contribution analysis, we can also assign graph variation features to each sample by the designed framework, or even each gene, and then CCA of gene under a certain signature can be also obtained. Yet, compared with the previously published software, 23 our framework spends more time on a de novo extraction of signature analysis. The reason is that we design a correction process "Re-updated the initial value" and a solution space process "build a solution space.". Hence, this is a weakness in the framework that need to be optimized in the future. However, our framework provides a new idea for understanding the panorama of tumor mutational processes, and help scientific researchers to study the mechanisms of tumor progression. It is of great application value to study the characteristics and statistical distribution of one gene under a certain signature by assigning the signature. In this study, F I G U R E 7 Application of cumulative contribution abundance matrix of genes in prognosis for esophageal squamous cell carcinoma. The sensitivity-specificity curve of set classifier at the test sample sets for urgent versus non-urgent binary classification based on previous reports 23,24,47 and statistical evaluation (Figure 4b, Figure S6B), we confirmed the background of expected mutational signatures of ESCC. Concurrently, 717 cancer-related genes from the COSMIC census (https://cancer.sanger.ac.uk/census, Table S5) were selected as the base to calculate CCA. Due to the inconsistent sequencing background of the data, we uniformly analyzed the data of coding regions, and calculated the contribution of each sample to these signatures for subsequent analysis. However, this analysis will have some limitations. To avoid the impact of this limitation, we are committed to explore the ubiquitous signatures such as APOBEC signatures and age signatures, and discuss those frequently mutated genes that present in ESCC, such as TP53 and PIK3CA.
Previous reports have revealed that mutations in the helix domain and kinase domain of PIK3CA cause activation through different mechanisms, and the mutation process may be related to driving mutations in a variety of cancers. 55 In ESCC studies, it was also mentioned that PIK3CA mutation was associated with the APOBEC signatures. 50,51 Furthermore, we found that there was a close relationship between APOBEC signatures and PIK3CA mutations in the meta-analysis of 1073 ESCC tumors, especially PIK3CA mutation E545K. In conclusion, APOBEC are strongly linked to driver mutations, especially the known hotspot mutation E545K of PIK3CA, suggests that there is a mechanism between APOBEC and PIK3CA, which needs to be deeply verified by a large number of subsequent experiments.
Next, we analyzed the association of TP53 hotspot mutations and mutational signature. Notably, the most frequent TP53 mutations found in ESCC were associated with the most commonly observed mutational signature, age signature, which reflects the natural degradation of 5-methylcyto sine to thymine. 56 In particular, the mutation R342* of TP53 can affect the mutation process of tumor occurrence and development, resulting in a significant increase contribution of the sample to age signature ( Figure 6B). This led us to put forward the hypothesis that the mutation R342* of TP53 in ESCC, which is the primary factor to increase the activity of age signature. In many tumor types, driver mutations of TP53 appear to be strongly associated with multiple signatures, and their probably arises due to the selection of loss-of-function and dominant-negative alleles, which are generated by specific mutational processes. 51,52 In our analysis, TP53 mutations R213* was also shown to be an independent prognostic factor.
In addition to the CCA matrix image, deep convolutional neural network denseNet-121 was used to analyze and the CCA matrix image data of gene in SBS3* and SBS17b*, which can be preliminarily distinguished the shortened survival outcome (follow-up time no less than 3 years). This finding indicates that the results of this method serves as one of the criteria to evaluate the prognosis of 3-year survival. Combined with AI technology, our designed scheme directly provides a new way to explore from the single gene relationship research to the multi gene association analysis. In general, the number of individuals studied in this paper is relatively small, which is one of the shortcomings of model learning. We hope to further achieve the useful information in the era of big data. More excellent learning model is also one of the improvement ideas to obtain accurate results, which needs further research in the future. At the same time, complete clinical information, including treatment methods, medication information, extended follow-up time and so on, is of great clinical significance for the further exploration of this idea. We are reasonably optimistic that in the future, CCA matrix of genes can be used to evaluate the prognosis, metastasis risk, recurrence risk, and even provide medication guidance and suggestions for individuals.
Overall, it is indispensable to understand and explore the mechanism of tumorigenesis and development by studying the relationship between genes and mutational signatures. The potential application of CCA of genes needs to be further studied and explored, such as giving some specific gene lists, forming image pictures, and perhaps evaluating prognosis and guiding medication through deep learning.

DATA AVAILABILITY STATEMENT
RNMF is publicly available at https://github.com/zhenz hang-li/RNMF and will be regularly maintained and updated. A detailed description of installing/running RNMF in the GitHub repository. The data sets used and/or analyzed during the current study are available from the corresponding author upon reasonable request.