Regulatory Network Analysis of Mutated Genes Based on Multi-Omics Data Reveals the Exclusive Features in Tumor Immune Microenvironment Between Left-Sided and Right-Sided Colon Cancer

Left-sided colon cancer (LCC) and right-sided colon cancer (RCC) have distinct characteristics in tumor immune microenvironment (TIME). Although existing studies have shown a strong association between gene mutations and TIME, whether the regulatory mechanisms between gene mutations and TIME are different between RCC and LCC is still unclear. In this study, we showed the fractions of CD8+ T cells were higher while those of regulatory T cells were lower in RCC. Besides, a stronger association between gene mutations and TIME was observed in RCC. Specifically, using multi-omics data, we demonstrated the mutations of most top mutated genes (TMGs) including BRAF, PCLO, MUC16, LRP2, ANK3, KMT2D, RYR2 made great contributions to elevated fraction of immune cells by up-regulating immune-related genes directly or indirectly through miRNA and DNA methylation, whereas the effects of APC, TP53 and KRAS mutations on TIME were reversed in RCC. Remarkably, we found the expression levels of several immune checkpoint molecules such as PD-1 and LAG3 were correlated with corresponding DNA methylation levels, which were associated with the mutations of TMGs in RCC. In contrast, the associations between gene mutations and TIME were less significant in LCC. Besides, survival analyses showed APC mutation had adverse impact on immunotherapy while patients with BRAF mutation were more suitable for immunotherapy in colon cancer. We hope that our results will provide a deeper insight into the sophisticated mechanism underlying the regulation between mutations and TIME, and thus boost the discovery of differential immunotherapeutic strategies for RCC and LCC.


INTRODUCTION
Colon cancer (CC) is the third leading cause of cancer and the second most common cause of cancer-related deaths worldwide (1). Due to the wide application of colonoscopy, the incidence rate of CC has seen a decline in recent years. However, patients with CC are still faced with a high risk of disease recurrence, leading to a major cause of CC mortality (2). Fundamentally, the colon develops from two separate embryonic sections of primitive gut: the midgut, which develops into the cecum, ascending colon, hepatic flexure and proximal two-thirds of transverse colon, is defined as the right (proximal) colon; and the hindgut, which gives rise to the distal third of the transverse colon, splenic flexure, descending colon and sigmoid colon, and is defined as the left (distal) colon (3). Tumor development likewise differs depending on its colon location, therefore, CC is classified separately into left-sided colon cancer (LCC) and rightsided colon cancer (RCC), and are recognized as two different types of CC with distinct clinicopathological characteristics and molecular features (4). As reported, RCC patients tend to be older, mucinous, undifferentiated, and have shorter survival time compared to the LCC patients (5). Additionally, the preferred metastasis sites differ between patients with LCC and RCC too. LCC patients are prone to have liver and lung metastasis, while RCC patients tend to have peritoneal carcinomatosis (6).
In the last decade, biologists and bioinformaticians have been trying to understand the heterogeneity of the etiology underlying the distinct clinical manifestation between LCC and RCC based on their multi-omics data including gene expression, DNA methylation, mutation and so on. As an example, a comprehensive investigation of multi-omics data showed that RCC is more hypermethylated and more aggressive relative to LCC (7). Recently, another multi-omics study revealed crucial crosstalk between calcium homeostasis and immune/GPCR signaling process was found only in LCC (8). However, the interplays among different kinds of biomolecular events are not well studied to date.
Tumor immune microenvironment (TIME), which includes immune cells and immune-related molecules in a tumor microenvironment, is an important factor to consider when trying to determine the choice of treatment strategies and to form reliable prognosis biomarkers for cancer (9,10). Previous studies have shown that the degree of immune infiltration in RCC is higher than that in LCC (11,12), which may be a major contributor to the differences of drug-sensitivity and prognosis. Tumor immunity is often regulated by genetic and epigenetic alternation. For example, the down-regulation of SLC64A in RCC results in the loss of immunosuppressive capacity (8). Meanwhile, the TP53 mutation often plays a negative role in tumor immunity in colon cancer and gastric cancer but reversely in breast and lung cancer (13). Furthermore, CC patients with both high tumor mutation burden (TMB) and BRAF mutation always present a high level of immune infiltration with an abundant expression of immune checkpoint molecules such as PD-L1 (14), and BRAF is a biomarker for immunotherapy in melanoma (15). Additional, immune infiltration-associated lncRNAs can also act as biomarkers for immunotherapy such as LINC01184 (16) and LINC02256 (17) et al. However, whether the regulatory mechanisms between gene mutation and TIME are different between RCC and LCC remains unclear.
In this study, we performed a systematic association analysis on the top mutated genes (TMGs) and TIME in both types of CC and constructed the TMG-associated regulatory networks based on multi-omics data, which includes gene and miRNA expression, DNA mutation and methylation, to identify the potential molecular mechanisms behind these phenotypic differences. Furthermore, we explored the effects of APC, BRAF and TP53 mutation on immunotherapy for CC patients. Our results will help us to better understand the tumorigeneses mechanisms of RCC and LCC, and discover new molecular prognostic and therapeutic targets.

Collection of Multi-Omics Data of CC Patients
All the multi-omics data of CC patients, including DNA mutation, gene expression, miRNA expression, DNA methylation, and clinical characteristics, were collected from The Cancer Genome Atlas (TCGA, https://cancergenome.nih. gov/) portal. Primary tumors located in the appendix to the transverse colon are defined as RCC, whereas tumors located from the splenic flexure to the sigmoid colon are categorized as LCC. 1134 CC patients with gene mutations were selected from the cBioPortal database, in which 1085 patients with prognosis and location information were used for survival analysis. Of these patients, 80 (37 RCC and 43 LCC) were treated with immunotherapy (https://www.cbioportal.org/).
For DNA mutation, we collected the Mutation Annotation Format (MAF) files that record the information of somatic variants processed from whole-exome sequencing data of 322 CC patients. For gene expression, the expression profiles of 428 CC samples and 39 normal samples were obtained. For miRNA, the expression profiles of 407 CC samples and 8 normal samples were obtained from the small RNA-seq data. For DNA methylation, we collected the beta values of 272 CC samples and 35 normal samples, which represent the ratio of the methylated probe intensity and total intensity detected from the Illumina HumanMethylation450 BeadChip Arrays. The clinicopathologic features such as age, gender, cancer stage, TNM stage, and survival data for 434 CC patients were also collected. The sample sizes for each type of dataset in the LCC and RCC cohort were shown in Supplementary Table 1.

Molecular Subtypes of CC
Microsatellite Instability (MSI) and chromosomal instability (CIN) subtypes information was downloaded from the cBioPortal database (https://www.cbioportal.org/). Among them, the samples were divided into three groups levels of MSI, denoted as microsatellite instability-high (MSI-H), microsatellite instability-low (MSI-L), and microsatellite stability (MSS) respectively. While CIN was classified into CIN and non-CIN. CpG island methylator phenotype (CIMP) was determined by the DNA methylation level of five markers (CACNA1G, IGF2, NEUROG1, RUNX3, and SOCS1) using a similar method as previously reported (18). In brief, the intensities of 307 CpG sites within the five genes were extracted from the DNA methylation profiles and used for unsupervised consensus clustering analysis. To confirm the superiority of CIMP subtype classification based on these five genes, we compared the results with those based on 7 genes of DNA methyl-transferases (DNMT1, DNMT3A, DNMT3B and DNMT3L) and demethylases (TET1, TET2 and TET3), which were well-acknowledged DNA methylation-associated genes. According to the clustering result, the CIMP phenotype can be divided into two categories, high level of CIMP (CIMP-H) and low level of CIMP (CIMP-L).

Evaluation of TIME for CC Patients
To evaluate the infiltration level of immune and stromal cells infiltration, we used ESTIMATE (Estimation of Stromal and Immune cells in Malignant Tumor tissues using Expression data) algorithm (19) to calculate the immune and stromal score for each sample using the expression profile, in which a higher immune/stromal score would represent a higher degree of immune/stromal cells filtration. The fractions for 22 immune cells were obtained from the expression profiles using CIBERSORTx software (20). The samples with p-value <0.05 were considered significant and selected for analyses. To further compare the proportions of different types of T cells between the two cohorts, ImmuCellAI was employed to evaluate the abundance of 18 T cells (21). Wilcoxon test was used to compare the difference of immune cell proportions between LCC and RCC.

Association Analysis of Clinical Features
For the association analysis of the clinical features and location, t-test was used to compare two groups of measurement data, Pearson's chi-squared test was applied for categorical data while Mantel-Haenszel's chi-squared test was used to compare two groups of ordinal data. For association analyses of the clinical features and immune score, Wilcoxon test was used to compare two groups while Kruskal Wallis H test was applied for multiple groups. The accepted level of significance is p<0.05. Kaplan-Meier survival curves were performed for prognostic differences and compared using the Tarone-Ware test.

TMGs Identification, Mutated Pairs Characterization, and Visualization
The analysis of the MAF files was conducted using the 'maftools' R package, which has various analytical and visualization methods that are commonly used in cancer genomic studies (22). Next, we obtained the mutation rate of each gene in patients with LCC and RCC respectively. The genes were sorted in descending order by the mutation rates, and the top 30 frequently mutated genes were obtained and defined as TMGs in each cohort. Furthermore, the mutated pairs of TMGs in a mutually exclusive or co-occurring manner were characterized using the Fisher's exact test on a 2×2 contingency table containing the frequencies of mutated and non-mutated samples. The p-values were adjusted using the false discovery rate (FDR) method. FDR<0.05 was a cutoff to select significant pairs. The TMGs associated networks were visualized using the Cytoscape software.

Tumor Mutation Burden Calculation
For each sample, we counted the number of somatic variants occurring throughout the coding regions of the sequenced gene and defined it as the somatic mutation frequency (SMF). Next, we normalized SMF per 30 megabase pairs (Mbp) and divided it by the total genomic territory sequenced. Therefore, TMB was calculated using the following formula:

Differential Expression Analysis of Genes and miRNAs
For gene expression, we first removed the low expressed genes with an average FPKM less than 0.1 or the samples with more than 20% of genes lacking expression values. Next, a differential expression analysis between cancer samples and corresponding normal samples was performed in RCC and LCC respectively using the 'DESeq2' R package. The genes with FDR<0.05 and |log2(Fold Change (FC)) | >log2(1.5) in each cohort were considered as the differentially expressed genes (DEGs). As for miRNA, we removed the miRNAs with low or absent expression and obtained the differentially expressed miRNAs (DEmiRNAs) in RCC or LCC by using the same method and cutoff as above-mentioned.

Differential Methylation Analysis
The processing of DNA methylation profiles was conducted using the Chip Analysis Methylation Pipeline (ChAMP) R package (23) in RCC and LCC respectively. First, the probes with more than 20% of samples lacking signals and the samples with more than 10% of probes lacking values were removed, in total, 395413 (81.4%) probes and 320 (100%) samples were eligible for further analysis. The normalization of probe signals was performed using peak-based correction (PBC) method (24). The probes with both FDR<0.05 and beta-value difference (△ß)>0.15 were considered as differentially methylated probes (DMPs).

Construction of TMG-Centric Regulatory Network
For each TMG, we divided the CC patients into two groups: mutation carriers and wildtype (non-mutation) carriers. Then, ttest was applied to reveal the TMG-associated DEGs by comparing the TMG-mutated and TMG-wildtype group in each cohort, and P-values were adjusted using the FDR method. The genes with FDR<0.05 and Log2(FC)>Log2(1.5) were defined as the TMGassociated genes (i.e., TMG-DEG pairs). After which, we applied the same method to identify the TMG-associated miRNAs (i.e., TMG-DEmiRNA pairs). To obtain TMG-associated DNA methylation probes, we used the t-test to determine whether the DNA methylation level of each probe was significantly different between the mutation and wildtype group for each TMG. The probes with FDR<0.05 and △ß >0.15 were defined as the TMGassociated probes (i.e., TMG-DMP pairs). Finally, the TMGcentric regulatory network was created using TMG-DEG, TMG-DEmiRNA, and TMG-DMP relationships, which were divided into two parts: 1) the negative relationships, i.e., the levels of the DEGs, DEmiRNAs, or DMPs in mutation carriers are lower than those in non-mutation carriers, which indicated that TMGs regulate these molecules negatively; 2) the positive relationships, where the regulation are in the reverse direction. Network visualization was performed using Cytoscape software.

miRNA-Target Relationships Prediction
We obtained the predicted miRNA-target relationships by combining the results from two sources: predicted using miRANDA software with default parameters and data that was directly downloaded from the miRDB database. Due to the inhibitory effect of miRNA on gene expression, we finally retained miRNA-target pairs with opposite differential expression trends for downstream analysis.

Functional Enrichment Analysis
GO functional enrichment analysis was performed by "ClusterProfiler" R package. The terms with FDR<0.05 and the gene number corresponding to this term not less than 3 were considered to be statistically significant.

Identification of Immune-Related Genes, miRNAs and Methylation Probes Globally
First, the CC patients in LCC and RCC cohort were classified into either the high-immunity group (immune score>=median) or the low-immunity group (immune score<median) based on the median value of their immune score in each cohort. The significant immune-related DEGs, DEmiRNAs and DMPs were determined by t-test between high-immunity and low-immunity groups with thresholds of FDR<0.05 and |log2(FC)| > |log2(1.5)|. The significant immune-related DEGs with log2(FC)>0 were considered as immune-promoting molecules, indicating these molecules are up-regulated in the high-immunity group compared to the low-immunity group. While the significant immune-related DEGs with log2(FC)<0 were thought as immunosuppressive molecules, meaning these molecules are down-regulated in the high-immunity group compared to the low-immunity group. Finally, by taking the intersection for the TMG association network (including TMG-DEG, TMG-DEmiRNA, TMG-DMP) obtained from the previous analysis and the significant immune-related molecules (DEGs, DmiRNAs, DMPs), we obtained immune regulatory network. Network visualization was performed using Cytoscape software.

Clinical Features and Tumor Immune Microenvironment of RCC and LCC Patients
A total of 434 CC patients from the TCGA-COAD portal were divided into LCC group (n=176) and RCC group (n=258) according to the location of their primary tumor. Then, we compared the clinical characteristics between the two cohorts. In line with existing research (25), RCC patients tend to be older and have a worse overall prognosis compared to LCC patients ( Table 1 and Figure S1). Next, we analyzed the associations between the location and other kinds of molecular subtype including microsatellite instability (MSI), chromosomal instability (CIN), and CpG island methylator phenotype (CIMP) which are classified based on their distinct genomic and molecular characteristics. Consistent with existing study (26), we also found that RCC group has a higher proportion of MSI-H and a lower proportion of CIN patients compared to LCC group ( Figures S2A, B and Table 1). As for CIMP, we determined the CIMP subtype for each patient using the b value of DNA methylation panel or 7 methylation-related genes respectively. The DNA methylation panel, which consists of five markers, CACNA1G, IGF2, NEUROG1, RUNX3, and SOCS1 (18) showed CIMP distribution of CC patients well. However, the 7 genes of DNA methyl-transferases or demethylases, including DNMT1, DNMT3A, DNMT3B, DNMT3L, TET1, TET2 and TET3, showed poor CIMP clustering result ( Figure S3B). Therefore, we distinguished CIMP subtypes based on the above five markers. The results showed that RCC group had a higher proportion of CIMP-H patients ( Figures S2C, S3A and Table 1). TIME plays a crucial role in tumor progression (27). First, we obtained the immune and stromal scores using the ESTIMATE algorithm (19), which measures the levels of immune and stromal cells infiltration in samples respectively. Next, we compared the two kinds of score in cancerous tissues between LCC and RCC patients. We found the immune scores of cancer tissues in RCC were higher than those in LCC ( Figure 1A). However, the differences in stromal scores were not significant ( Figure 1B), suggesting that immune cells infiltration may play a specific role in tumor growth, prognosis and therapy for RCC. Then, the associations between the immune scores and the clinical factors in both types of CC were investigated. In RCC, the immune scores were significantly higher in patients with MSI and non-CIN subtypes (P<0.001, Figures 1D, E), and female patients had slightly higher immune scores than males' (P=0.01, Figure 1C). However, there were no associations between immune scores and gender, MSI, and CIN in LCC ( Figures  1C-E). Furthermore, the immune scores were slightly different between patients with and without lymphatic invasion in RCC (P=0.041) but not in LCC ( Figure 1F). As the immune score is used to evaluate the overall level of immune cells infiltration, we estimated the fractions of 22 immune cells in detail to study the difference in TIME between LCC and RCC patients. As a result, CIMP-H, CpG island methylator phenotype-high; CIMP-L, CpG island methylator phenotype-low; CIN, Chromosome instability; Non-CIN, Nonchromosomal instability; Tis: In situ, non-invasive (confined to epithelium); T1: Small, minimally invasive within primary organ site; T2: Larger, more invasive within primary organ site; T3: Larger and/or invasive beyond margins of primary organ site; T4: Very large and/or very invasive, spread to adjacent organs; M0: No distant metastases; M1: Distant metastases present; N0: No lymph node involvement; N1: Regional lymph node involvement; N2: Extensive regional lymph node involvement; *P<0.05.  we found that the fractions of CD8+ T cells, activated NK cells, and M1 Macrophages were significantly higher in RCC than LCC, while those of immune repressive cell types such as regulatory T cells, M0 Macrophages were significantly lower ( Figure 1G). As T cells have many subsets with specific function, we further compared 18 T cell types between the two cohorts. As shown in Figure 1H, the proportions of CD4+ and CD8+ naïve, exhausted, Th17, central memory, effector memory, nature kill, and CD8+ T cells were all higher in the RCC group compared to the LCC group, while CD4+ T cells, cytotoxic, and regulatory T cells including Tr1 and iTreg, gamma delta were more abundant in LCC compared to RCC. These results revealed the different associations between TIME and molecular subtypes in LCC and RCC.

Identification of TMGs in LCC and RCC
The exome sequencing data of 129 LCC patients and 193 RCC patients were used to identify TMGs. First, comparing with LCC, RCC patients have a higher TMB ( x = 5.54 in LCC and x = 18.54 in RCC, P<0.05 by t-test), which was consistent with a previous report (28). Next, we selected the top 30 frequently mutated genes to be TMGs in each cohort. Among them, 18 are common, including some reported cancer-associated genes such as APC, TP53, TTN, KRAS, and PIK3CA (Figures 2A, B and Supplementary Table 2). As expected, most of these genes have higher mutation frequency in RCC compared to LCC, except for APC and TP53 (84% vs. 71%, 70% vs. 57% in LCC and RCC respectively, Figure S4), suggesting the potential dominant roles of APC and TP53 in LCC as previously shown (7). Previous studies have found that the driver genes in dysregulated pathways are often mutated in a mutually exclusive manner, while co-mutated genes linked by synergism promote cancer development in a synergistic manner (29). Next, we investigated the interactions among TMGs in LCC and RCC respectively. We found that the TP53/DNAH11 pair was the only mutually exclusive pair in LCC ( Figure 2C). In contrast, 31 mutually exclusive pairs were identified in RCC ( Figure 2D), of which APC, TP53, and KRAS were the most major exclusive genes and these three TMGs might play different roles from other TMGs in the initiation and progression of RCC. Remarkably, both BRAF and ANK3 were mutually exclusive in all three genes ( Figure 2D), which indicated that BRAF and ANK3 may contribute greatly to the development of RCC through the involvement of multiple key pathways. Regarding the co-occurring pairs, UNC13C, NEB, and ZFHX4 were the  major co-occurring genes in LCC ( Figure 2C), whereas a majority of TMGs in RCC interacted with each other in a cooccurrence manner ( Figure 2D), further proving the complicated tumorigeneses mechanism of RCC. Unexpectedly, no interactions were found for APC in LCC and AMER1 in RCC, suggesting they might have an independent role in LCC and RCC respectively. In summary, our results demonstrated that somatic interactions were significantly more frequent in RCC compared to LCC, which might be one of the factors underlying the poorer prognosis observed in RCC patients.

Association Analysis of TMGs and TIME
Next, we explored the association between the mutation status of each TMG and TIME in each cohort. In RCC, mutations of 20 TMGs (20/30, 66.7%) were observed to be related to the immune score, and that patients with mutations in most of TMGs have a higher immune score, except for the mutations of TP53, APC, and KRAS ( Figure 3A). As previously stated, we found that TP53, APC, and KRAS were main mutually exclusive TMGs in RCC ( Figure 2D), suggesting the three genes have different roles on TIME in RCC. However, only 2 TMGs (APC and USH2A, 2/ 30, 6.7%) were associated with the immune score in LCC ( Figure  3B), suggesting TMGs in RCC have a greater impact on TIME compared to LCC. Strikingly, a positive trend of immune score, although not significant (P=0.125), was observed in LCC patients with TP53 mutation (Figure 3B), which was contrary to its effect in RCC.
To better understand the association between TMGs and TIME, we studied the relationships between the mutation status of each TMG and the fractions of 22 immune cells or 18 T cells. In RCC, 27 TMGs were found to be associated with at least one type of immune cell ( Figure 3C), and 28 TMGs were associated with at least one type of T cell ( Figure 3E). For example, the infiltration of CD8+ T cell was positively correlated with 11 TMGs such as PCLO, CSMD3, ANK3 and LRP2 et al., while that of resting CD4 memory T cell was negatively correlated with 12 TMGs such as DNAH11, KMT2D, CSMD3 and LRP2 et al. (Supplementary Table 3). Besides, we found that the proportions of cytotoxic, exhausted, nTreg, iTreg, Th1 and gamma delta cells were positively correlated with 11 TMGs,   22 TMGs were found to be associated with at least one type of immune cell ( Figure 3D), and 26 TMGs were associated with at least one type of T cell ( Figure 3F). For example, the fraction of resting CD4 memory T cell was positively correlated with USH2A and ARID1A mutation, while that of CD8+ T cell was negatively correlated with CSMD1 mutation (Supplementary Table 3). In addition, we also observed the relationships between 18 T cell types and TMGs in LCC were significantly poorer than those in RCC. Obviously, the fractions of cytotoxic, exhausted, nTreg, iTreg and Th1 cells were positively correlated with partial TMGs, while those of NKT, CD4+ T cells and CD8+ naïve T cells were negatively correlated with some TMGs (Supplementary Table 3). In conclusion, TMGs may play a crucial role in RCC TIME.

TMGs Influence TIME Through Regulation of Gene Expression
To illustrate the distinct regulatory mechanisms of TMGs influencing TIME, we further explored the genes whose expression levels were associated with TMGs in the two types of CC.  Table 4). Among them, 88 TMG-DEG pairs were common in both LCC and RCC. The number of TMG-DEG pairs in RCC is much higher compared to LCC, suggesting a stronger relationship between TMG and gene expression regulation in RCC. Furthermore, even the common TMGs between LCC and RCC have a different number of associated genes ( Figure 4C), further hinting that the regulatory network of TMGs is different between LCC and RCC. We then attempted to find the TMGs associated processes that were altered in LCC and RCC, we performed functional enrichment analysis on the TMG-upregulating and TMG-downregulating genes in both types of CC. The results showed that there were 464 TMGupregulating genes in RCC, and they were mainly associated with the immune system including T cell activation, regulation of immune effector process, and positive regulation of cytokine production ( Figure 4D), while TMG-downregulating genes involved with the genes that were related to various metabolic, blood and heat-related processes such as steroid metabolic process, cellular hormone metabolic process, regulation of blood pressure, heat generation, etc. ( Figure S5A). As for LCC, TMG-downregulating genes showed significant enrichment of terms associated with muscle contraction, chondrocyte differentiation, cartilage development, etc. (Figure S5B), and no enrichment results were found in the TMG-upregulating genes, which was likely because only 19 genes were involved. To further explore the functions of DEGs positively regulated by TMG mutations in LCC, we picked the top 1000 genes with highest FC to perform functional enrichment analysis ( Figure  S5C), finding that these genes were also related to the humoral immune response. The results revealed that mutations of TMGs can up-regulate immune-related genes in both LCC and RCC, but the effect is much wider and stronger in RCC than in LCC. The results from a global association analysis of TMGs mutation and TIME led us to speculate that some TMGs may have a strong association with TIME and this relationship may be different between LCC and RCC. To link TMGs to specific immune processes, we conducted functional enrichment analyses for the up-regulated genes and down-regulated genes of each TMG in LCC and RCC respectively. We found that 12 TMGs in RCC and only TP53 in LCC were associated with at least one immune-related process based on the enrichment results of corresponding up-regulated genes for each TMG ( Figures 4E, F and Supplementary Table 5). However, based on the enrichment results from the down-regulated genes of each TMG, no TMG was found to be associated with the immunerelated process in both LCC and RCC. These results further indicated that more TMGs in RCC were related to an elevated expression level of immune-related genes. For example, BRAF had the most associated genes in RCC. We found that 310 upregulated genes associated with BRAF mutation were enriched in the genes related to T cell activation, response to interferongamma, positive regulation of cytokine production, and other immune-related processes ( Figure S6A). While its 637 downregulating genes were mainly involved in neuron recognition, negative regulation of the Wnt signaling pathway, lipid transport, etc. ( Figure S6B). Some associated genes of BRAF mutation had been reported to have important functions in CC. For instance, lncRNA SATB2-AS1 was found to be up-regulated when BRAF mutation occurs in RCC in our study, and it was recently reported to inhibit tumor metastasis and affect tumor immunity in CC (30), indicating that BRAF plays a crucial role in TIME. Surprisingly, we found that TP53 mutations were positively associated with 11 genes in LCC, which were enriched with the genes relating to T cell activation and leukocyte differentiation ( Figure 4F). Considering the negative role of TP53 to TIME in RCC as shown above ( Figure 3A), we hypothesized that TP53 may play a positive role on TIME only in LCC.

TMGs Influence TIME Through Regulation of miRNAs
The regulatory network in biology is exceedingly complex as any alterations in gene expression could be initiated or influenced by multiple regulators such as transcription factors, miRNAs, etc. MiRNA is a kind of ncRNA that modulates the expression levels of up to 60% of the protein-coding genes at both the transcription and post-transcription level (31,32). So, the associations of TMG-DEG relationships identified above may be indirect and connected through miRNAs. Based on small RNA-seq data, we first investigated the DEmiRNAs in LCC and RCC by comparing their expressions level in cancer tissues to those in the normal control tissues. We found that 209 and 167 miRNAs were up-regulated and down-regulated in RCC, while only 119 up-regulated and 121 down-regulated miRNAs were found in LCC. Next, we identified the TMG-DEmiRNA relationships in LCC and RCC using a similar method which was used to identify TMG-DEG pairs. As expected, the number of significant TMG-DEmiRNA pairs in RCC was much more than that in LCC (109 vs. 5), further suggesting the tumorigenesis mechanism in RCC may be more complex. In RCC, 13 positive and 96 negative TMG-DEmiRNA associations were obtained ( Figure 5A and Supplementary Table 6), of which 4 DEmiRNAs were common. While in LCC, only 1 positive and 4 negative TMG-DEmiRNA pairs were found ( Figure 5B and Supplementary Table 6). Among them, miR-552-5p was associated with the most TMGs in RCC. Although no study has reported on the function of miR-552-5p in CC until now, its strong relationship with cancerous pathways such as cell proliferation and metastasis and its potentiality in prognostic prediction have been widely acknowledged (33,34), suggesting that miR-552-5p likely plays an important role in CC.
To further understand the roles of TMGs regulating DEmiRNAs, we attempted to identify the potential targets of these DEmiRNAs. Considering that most known miRNA-target pairs were experimentally verified in cancer-related processes, there may be some biases if we included all the known targets of DEmiRNAs. Therefore, we only considered the miRNA-target relationships that were predicted by miRANDA (35) or recorded in miRDB in this study (36). To improve prediction accuracy, we also required that the DEmiRNA and its target within the same pair should be differentially expressed in the opposite direction in each cohort (If DEmiRNAs were up-regulated, target genes were down-regulated, vice versa). Finally, 1475 and 4836 miRNA-target relationships were obtained for the up-regulated and down-regulated DEmiRNAs of TMGs mutation in RCC respectively, among them, 1093 and 2433 target genes were involved respectively (Supplementary Table 7). However, only 854 miRNA-target relationships with 741 genes were found in LCC (Supplementary Table 7), of which 148 and 641 are targets of up-regulated and down-regulated DEmiRNAs of TMGs mutation. In RCC, 788 genes were the common targets of both TMG-upregulating and TMG-downregulating DEmiRNAs. Functional enrichment analysis showed that the enriched terms were Wnt signaling pathway and some neuron-related biological processes ( Figure 5D middle), suggesting that these pathways may be important in tumorigenesis and are regulated by multiple factors. Specifically, we found the unique targets of TMG-downregulating DEmiRNAs in RCC were involved in immune response ( Figure 5D left). However, no immunerelated biological process was over-represented in the targets of DEmiRNAs in LCC ( Figure S7), suggesting that some immunerelated genes may be up-regulated through the down-regulation of DEmiRNAs by TMGs mutations in RCC.
To determine which mutation of TMG affects TIME through the regulation of miRNAs, we also performed a functional enrichment analysis on the targets of each DEmiRNA. We found the following, 15 out of 38 TMG-associated DEmiRNAs (39.5%) in RCC and only 2 TMG-associated DEmiRNAs (40%) in LCC were related to at least one immune-related process (Supplementary Table 8 and Figure 5C). In detail, 6 immunerelated DEmiRNAs were up-regulated in the mutation groups of 5 TMGs including APC, TP53, DNAH11, BRAF, and KMT2D, while down-regulation of 11 immune-related DEmiRNAs were associated with mutations of 15 TMGs in RCC ( Figure 5E). The major immune-related processes of these DEmiRNAs were T cell activation, T cell migration, leukocyte activation, etc. ( Figure 5E). As an example, miR-183-5p was annotated with T cell activation, T-helper 17 type immune response, T-helper 17 cell differentiation, etc. Existing reports have also confirmed its role in the regulation of Th17 differentiation (37). In this study, we found the mutation of ZFHX4 promoted down-regulation of miR-183-5p in RCC. Additionally, miR-148a-3p has been reported to regulate PD-L1 in CC (38) and functions as an important immune regulator (39). We also found that the mutation of PCLO was associated with low expression of miR-148a-3p. However, only two associations (APC and miR-92A-3p, USH2A and miR-196a-5p) were found in LCC ( Figure 5E). Our results showed a larger crosslink between genetic alternation and miRNA expression aberrance for TIME regulation in RCC.

TMGs Influence TIME Through Regulation of DNA Methylation
DNA methylation is a major type of epigenetic modification and plays a crucial role in tumorigenesis (40). According to the previous reports (41), DNA methyltransferases (DNMT family) and demethylases (TET family) often mediate DNA methylation abnormalities through co-expression with other genes. To explore the potential mechanism of DNA methylation changes, we further analyzed whether DNMT family (DNMT1, DNMT3A, DNMT3B, DNMT3L) or TET family (TET1, TET2, TET3) expression were affected by TMG mutations. In RCC, expression of DNMT1 and TET2 were activated by mutations of 15 and 5 TMGs respectively, while expression of DNMT3A, DNMT3L and TET1 were inhibited by mutations of 1, 12, 12 TMGs respectively ( Figure 6A). However, the associations between methylation-associated genes and TMGs were weaker in LCC than those in RCC ( Figure 6B). Next, we conducted weighted gene co-expression network analyses (WGCNA) based on the expression datasets of LCC and RCC by "WGCNA" R package. In RCC, a total of 6 modules from 3380 genes were identified in the co-expression network, while 9524 genes were not assigned to any modules (presented as grey color). As shown in Figure 6C, all modules were significantly associated with the expression level of TET2, whereas no module was correlated to DNMT3L. Furthermore, we found that TET family had a stronger correlation with most modules than DNMT family. Functional enrichment analysis showed that immune related BPs were associated with yellow module (Figure 6E), and the module was positively correlated with DNMT1 and TET2, which were also up-regulated by TMGs. Therefore, we hypothesized that expression levels of the two methylation-associated genes may be involved in immunerelated processes. In LCC, a total of 25 modules from 12909 genes were identified in the co-expression network, while 1824 genes were not assigned to any modules (presented as grey color, Figure 6D). Based on functional enrichment analysis, we found that immune related BPs were related to blue module ( Figure  6F), and the module was positively correlated with DNMT3A and TET family. Notably, TET2 was positively correlated with immune-related module in both RCC and LCC, which may play an important role in TIME.
Excepting for affecting methylation related genes, the mutations of specific genes may trigger DNA methylation alternation directly (42). After performing an analysis of the DNA methylation profiles in LCC and RCC from TCGA portal, we initially obtained 177996 and 137239 DMPs in RCC and LCC respectively. Next, 5506 and 52858 TMG-DMP pairs were identified in LCC and RCC with FDR<0.05 and △ß value >0.15 respectively. The same as before, the number of TMG-DMP associations identified in RCC was much higher compared to LCC. Considering that DNA methylation disruption usually causes an aberrant expression of the corresponding gene, we calculated the PCC between the methylation level of DMP and the expression level of its corresponding gene in each cohort to further filter the TMG-DMP associations. 911 and 16895 TMG-DMP pairs in LCC and RCC respectively were remained, with the threshold of an adjusted P-value less than 0.05 (Supplementary Table 9). Interestingly, up to 95.5% in RCC and 71.1% in LCC of TMG-DMP associations were positive ( Figure 7A and Figure S8A).
When taking the PCC of DNA methylation signals and the corresponding gene expression levels into account, we found that the DMPs with positive associations with both TMG mutations and corresponding gene expression, and the DMPs with negative association with both TMG mutations and corresponding gene expression were related to immune-related processes such as T cell activation, positive regulation of leukocyte proliferation in RCC ( Figures 7B, E). However, no immune-related processes were over-represented in other groups of TMG-DMP associations in RCC (Figures 7C, D) and any group in LCC ( Figure S8). This suggested that TMG mutation can also up-regulate immunerelated genes via the disruption of the corresponding DNA methylation level in RCC. Figure 7F shows the sub-network of 'T cell activation' formed by positive TMG-DMP associations with DMP also positively correlated with gene expression in RCC. Several known immune-related genes, especially the immune therapy checkpoint PDCD1 (also called PD-1) and its ligand PDCD1LG2 (also called PD-L2) were involved. Our results showed that 4 probes located in the body region of PD-1 (cg10526431, cg09319815, cg07281781, and cg06291111) were up-regulated by 12 TMGs, while 1 probe (cg14351952) located in the 5'UTR of PD-L2 was up-regulated by 5 TMGs, which in turn led to an elevated expression of PD-1 and PD-L2 in RCC ( Figure 7F). Furthermore, another well-known immune checkpoint molecular LAG3 was also found to be up-regulated  Figure 7F). Figure S9 shows the sub-network of 'positive regulation of leukocyte proliferation' caused by TMG-DMP negative associations with DMP inverse correlation with gene expression in RCC. A total of 7 TMGs and 8 immune-related genes were involved. Interestingly, we found several distinct methylation probes in HLA-DPA1 were both regulated by TMG mutation positively and negatively, where cg22941409 was positively regulated by mutations of ADGRV1, ANK3, DNAH10, DNAH11, LRP1B, LRP2 and KMT2D, and cg01804934 was inversely regulated by mutation of BRAF ( Figure S10A). However, the methylation level of cg22941409 was positively correlated with the expression level of HLA-DPA1 ( Figure S10B) but cg01804934 was reversed ( Figure S10C), both resulted in the elevated expression of HLA-DPA1 in RCC. Furthermore, we also found 2 probes located in HLA-DPB1 that were down-regulated by mutations of BRAF and ANK3 while 6 probes were up-regulated by mutations of 13 TMGs ( Figure S11A), all of which also resulted in a higher expression of HLA-DPB1 in RCC ( Figures S11B, C). These results showed that the crosslink between DNA methylation and gene mutation is important for TIME in RCC.

Integrative Analysis of Immune-Related TMG-Centric Regulatory Network
According to the results above, we can make a conclusion that the mutations of some TMGs in RCC are associated with immune-related genes directly or indirectly through miRNA and DNA methylation. We found 20, 12, 18, 24 immunerelated TMGs in the TMG-IS correlation analysis, TMG-DEG correlation analysis, TMG-DEmiRNA correlation analysis, and TMG-DMP correlation analysis respectively. Among them, 7 TMGs (BRAF, PCLO, MUC16, LRP2, ANK3, KMT2D, RYR2) were common in the four kinds of analysis ( Figure 8A). While in LCC, only 3 common immune-related TMGs were observed (APC, TP53, USH2A, Figure 8B).
In order to get a greater insight into the regulatory mechanism of TMGs on TIME, we constructed an immunerelated TMG-centric regulatory network in RCC and LCC respectively. In RCC, the immune-related TMG-centric regulatory network consists of 4289 TMG-DEG, 44 TMG-DEmiRNA, and 213 TMG-DMP (70 TMG-corresponding genes) relationships. As shown in Figure 8, most targets in the positive TMG-DEG relationships were associated with high expression levels in the high-immunity group ( Figure 8C), while the genes within negative TMG-DEG relationships were associated with low expression levels in the high-immunity group ( Figure 8D). In contrast, APC, TP53, and KRAS exhibited opposite relationships, further verifying that the three genes play a negative role on TIME in RCC ( Figures 8C, D). However, in LCC, only 327 TMG-DEG relationships and 1 TMG-DMP relationship were identified ( Figures 8E, F). Of note, while TP53 exhibited a negative correlation with immunity in RCC, it seemed to play a positive role on TIME in LCC as the targets in the positive relationships were upregulated in the high-immunity group ( Figure 8E) while those involved in the negative relationships were down-regulated in the high-immunity group ( Figure 8F). The conclusion runs contrary to previous findings that TP53 mutation suppresses immunity in colon cancer (13).

Effects of APC, BRAF and TP53 Mutations on Immunotherapy Efficacy in CC
The analyses of TMG regulatory network based on multi-omics data above suggest some TMGs such as APC have suppressive effects on cancer immunity while some TMGs such as BARF play opposite roles. It led us to hypothesize that TMG mutations may be served as predictive biomarkers for immunotherapy efficacy in CC. To directly address this point and validate our hypothesis, we collected the survival data of CC patients who received nonimmunotherapy or immunotherapy, with mutations of some TMGs known, from cBioPortal database. Then we investigated the impacts of TMG mutations on the efficacy of immunotherapy through Kaplan-Meier analyses.
As expected, we found immunotherapy can effectively improve the overall survival (OS) for non-APC-mutated patients (P=0.009, Figure 9A) but not for APC-mutated patients in CC (P=0.606, Figure 9A). On the contrary, for CC patients with non-immunotherapy treatment, APC mutation is a favorable prognostic factor (P<0.001, Figure 9A). All the findings were also established in RCC cohort ( Figure 9B). A slight improvement of OS for non-APC-mutated patients with immunotherapy was also observed compared to APC-mutated patients who treated with the same immunotherapy in LCC, although not significant due to the small sample size ( Figure 9C), suggesting non-APC mutation can be considered as a biomarker for immunotherapy efficacy in both LCC and RCC.
In contrary to APC, BRAF mutation is associated with high cancer immunity. Not surprisingly, we found BRAF-mutated patients with immunotherapy have higher OS time than those with non-immunotherapy (P=0.049, Figure 9D). While the OS is not better for non BRAF-mutated patients with immunotherapy than those with non-immunotherapy (P=0.225, Figure 9D). Unfortunately, we didn't find significant differences between the OS of the BRAF-mutated patients with immunotherapy and those with non-immunotherapy in RCC or LCC cohort, it is most likely because of small sample size of BRAF-mutated patients with immunotherapy. However, the similar tendencies were also observed ( Figures 9E, F), indicating BRAF mutation is also a biomarker of immunotherapy efficacy for CC patients.
However, mutation of TP53 had no effect on the efficacy of immunotherapy in all CC patients, RCC or LCC patients ( Figures 9G-I). Therefore, TP53 cannot be considered as a biomarker for immunotherapy efficacy in CC.

DISCUSSION
Tremendous variations in the microbiome, clinical factors, microenvironment, genomic characteristics and molecular events have been reported between LCC and RCC in the last decade (43). Among them, molecular event, which includes multiple kinds of alternations such as DNA variation, gene expression, and DNA methylation, is the most important mechanism underlying between LCC and RCC. However, the relationships among different types of molecular events, and their impact on TIME in LCC and RCC has not been systematically studied. In this study, we revealed the distinct differences in the immune cell infiltration levels especially CD8+ T cells between RCC and LCC. Then, we classified the 30 most frequently mutated genes as TMGs in each cohort, and analyzed the relationship between TMG mutation and the alternations in the gene expression, miRNA regulation, and DNA methylation. Our results showed that the mutations of TMGs in RCC can regulate more genes and affect more pathways compared to LCC, which was consistent with a previous study that demonstrating RCC had a more diverse mechanism in tumorigenesis (7). Additionally, we also showed that TMG mutation was strongly associated with TIME in both cohorts, and especially in RCC. Based on the above relationships, we further validated the conclusion that non-APC mutation and BRAF mutation can serve as biomarkers for immunotherapy but not TP53. Until now, a growing body of evidence suggests that gene mutation is strongly associated with TIME (44). In a previous study on lung adenocarcinoma, several commonly mutated genes were predicted to be neoantigens of the major histocompatibility complex (MHC) Class I and Class II molecules (45), and some of these genes are also immunerelated TMGs in RCC, including MUC16, ZFHX4, FAT3, TTN, and RYR2. Specifically, RYR2 and RYR3, the members of ryanodine receptors (RyR), have been reported to be associated with T cell activation (46). In this study, we found that mutations of most TMGs were associated with a higher expression of immune-related genes, among which 7 immune-related TMGs consisted of BRAF, PCLO, MUC16, LRP2, ANK3, KMT2D, RYR2 were observed in multiple methods. BRAF was a well-known regulator in CC and its mutation has been considered to be a reliable biomarker for prognosis and therapy (47). Regarding ANK3, although it is currently not well studied, existing studies have suggested that it is a potential therapeutic target oncogene in CC (48). However, the genes such as APC, TP53, and KRAS exhibit a negative effect on TIME in RCC. Interestingly, our study showed that TP53 may promote tumor immunity in LCC. Mutation of TP53 was associated with a higher expression of 11 genes, among which, 8 genes were up-regulated in the highimmunity group of LCC and were enriched in the gene set of T cell activation and leukocyte differentiation. Meanwhile, the unique down-regulated gene in the group of TP53 mutation in LCC, FOXA1, which has been reported to inhibit T cells proliferation in lung cancer (49), was also observed to be downregulated in the high-immunity group of LCC ( Figure 8F), indicating that TP53 mutation may exhibit a suppressive function against immunity inhibition in LCC.
Acting as a key regulator at the post-transcriptional level, miRNAs participate in many biological processes including the immune system, and play a major role in the cancer-immune response (50). We found 15 immune-related miRNAs that were associated with TMGs mutations in RCC. Among them, TP53 and APC could only up-regulate immune-related miRNAs, suggesting that they may be able to down-regulate immunerelated genes through miRNAs ( Figure 5E). However, the genes such as KMT2D, DNAH11, and BRAF, which have a positive correlation with tumor immunity level, also up-regulated several immune-related miRNAs ( Figure 5E). They may function together with other miRNAs by forming negative feedback loops to regulate immune-related genes.
DNA methylation is one of the universal epigenetic types which can regulate gene expression and genomic stability. CIMP, one of the CC molecular subtypes, is found to be correlated with BRAF mutation and MSI phenotype in the early years (51), suggesting that there is a relationship between gene mutation and the disruption of DNA methylation in CC. Remarkably, we found that the DNA methylation levels of the immune molecular checkpoints such as PD-1 and LAG3 were regulated by mutations of some TMGs in RCC, which in turn influenced the expression of their corresponding genes and resulted in higher immunogenicity. Therefore, mutations of these TMGs may act as biomarkers of immunotherapy and immune cell infiltration.
The identified studies indicated that POLE and POLD1 mutations are biomarkers for immunotherapy across multiple cancer types (52), and BRAF mutation can improve antitumor activity of immunotherapy in melanoma (15). In our study, we found that APC mutation exerted a suppressive effect on immunotherapy while the role of BRAF mutation on immunotherapy are inversely. Taken together, we propose immunotherapy benefits for CC patients with non-APC mutation or BRAF mutation. However, TP53 mutation had no effect on the efficacy of immunotherapy.

CONCLUSIONS
In conclusion, we systematically compared the associations between the mutational profiles and tumor immunity based on the multi-omics data in LCC and RCC respectively and found several distinct regulatory mechanisms underlying the TIME in both types of CC. Our results will provide valuable clues that would be helpful for the identification of novel targets for immunotherapy.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/ Supplementary Material.

ETHICS STATEMENT
The studies involving human participants were reviewed and approved by School of Medicine, Ningbo University. The patients/participants provided their written informed consent to participate in this study.

AUTHOR CONTRIBUTIONS
QL and TY wrote the manuscript. YZ, YL, SH, and HR did the bioinformatics analysis. QL, XD, and GZ designed the idea of this study. DN improved English writing in manuscript. JL, XF, YYX, YX, JM, LeC, CC, SN provided some suggestions for the study. MY and LvC contributed to the revised manuscript and provided the corresponding data analysis. All authors contributed to the article and approved the submitted version.