Integrative Analysis Reveals Comprehensive Altered Metabolic Genes Linking with Tumor Epigenetics Modification in Pan-Cancer

Background Cancer cells undergo various rewiring of metabolism and dysfunction of epigenetic modification to support their biosynthetic needs. Although the major features of metabolic reprogramming have been elucidated, the global metabolic genes linking epigenetics were overlooked in pan-cancer. Objectives Identifying the critical metabolic signatures with differential expressions which contributes to the epigenetic alternations across cancer types is an urgent issue for providing the potential targets for cancer therapy. Method The differential gene expression and DNA methylation were analyzed by using the 5726 samples data from the Cancer Genome Atlas (TCGA). Results Firstly, we analyzed the differential expression of metabolic genes and found that cancer underwent overall metabolism reprogramming, which exhibited a similar expression trend with the data from the Gene Expression Omnibus (GEO) database. Secondly, the regulatory network of histone acetylation and DNA methylation according to altered expression of metabolism genes was summarized in our results. Then, the survival analysis showed that high expression of DNMT3B had a poorer overall survival in 5 cancer types. Integrative altered methylation and expression revealed specific genes influenced by DNMT3B through DNA methylation across cancers. These genes do not overlap across various cancer types and are involved in different function annotations depending on the tissues, which indicated DNMT3B might influence DNA methylation in tissue specificity. Conclusions Our research clarifies some key metabolic genes, ACLY, SLC2A1, KAT2A, and DNMT3B, which are most disordered and indirectly contribute to the dysfunction of histone acetylation and DNA methylation in cancer. We also found some potential genes in different cancer types influenced by DNMT3B. Our study highlights possible epigenetic disorders resulting from the deregulation of metabolic genes in pan-cancer and provides potential therapy in the clinical treatment of human cancer.


Introduction
Unlike normal cells that regulate their cell division reasonably, cancer cells have the ability to sustain growthpromoting signals, which lead them to obtain enough energy through adjusting epigenetics and metabolism mechanisms to meet the sustained cell division. In recent years, there has been a growing interest in cancer metabolism, especially the theory of carbohydrate metabolism disturbance which was proposed by Otto Warburg in the 1920s [1]. Besides glucose metabolism, other metabolism disorders including lipid and amino acid metabolisms have also been discovered in cancer [2][3][4]. In addition to providing energy, the research showed that disordered genes encoding metabolic enzymes may also promote tumorigenesis through other biological functions like epigenetic modification [5]. e most studied epigenetic alterations associated with carcinogenesis were variation in DNA and histone structure through posttranslational modifications and histone variants. ereinto, DNA methylation and histone acetylation were the most common modifications related to the occurrence and development of cancer. Numerous studies have found that activation or overexpression of oncogene and silencing of cancer suppressor gene were due to the apparent epigenetic modification at their corresponding locations. e cancer suppressor gene BRCA1 was found hypermethylated in ovarian cancer patients [6], and CDH13 was methylated in endometrial carcinoma [7]. In the opposite case, hypomethylation may lead to the activation of normally silenced oncogenes. RAS, involved in cell differentiation, proliferation, apoptosis, cellular adhesion, and migration, was hypomethylated in the promoter in hepatocellular carcinoma [8].
Although separate investigation on these energy metabolism reprogramming and epigenetic modification dysfunction would uncover molecular characteristics, the diagnosis and cure effect that targeted them were not as expected [9,10]. e metabolic network was observed to interplay with many altered signaling pathways in cancer, which increased the difficulty of developing targeted metabolic treatment [11]. In addition to oncogenic signaling, the extensive crosstalk was revealed between the metabolic networks and epigenetics disorders in cancer. For example, with the abnormal condition, acid environment or metabolic rewriting, epigenetic landscape would undergo dysfunction which contributed to gene expression disorder and neoplastic processes [12]. Metabolism has a significant influence on epigenetics through adjusting DNA and histone modification enzymes [13]; in turn, the epigenetics would affect the metabolism by altering the gene expression.
ese studies usually focused on the relationship between metabolism and epigenetics just in one type of disorder or in an isolated tumor type [14,15]. In these specific analyses, the epigenetic alterations were a direct result of the altered activity of a particular metabolic enzyme. However, whether metabolic genes might contribute to epigenetic dysregulation more universally across multiple cancer types is not yet fully clear.
Given that anabolism and catabolism of metabolites and transcriptional dysregulation of metabolic genes are dramatically altered in cancer, there is an urgent issue to explore how such metabolic reprogramming alters epigenetic regulation at a global landscape. Researchers demonstrate that the gene expression patterns indeed reflected metabolic activities [16]. e projects such as the Cancer Genome Atlas (TCGA) provide large amounts of transcriptomic profiles across multiple tumor types with well-annotated human cancer samples [17], which is a good resource to explore metabolic dysfunction at the transcriptional level. e sufficient tumor and normal sample size with transcriptomic profiles are necessary for revealing relatively accurate metabolic dysfunction in each cancer. After searching the TCGA database by sample size for each cancer, we have performed a metabolic gene expression signature analysis covering 11 tumor types (normal samples of each cancer >30) including 5726 samples to identify critical metabolic gene features and their related epigenetic dysfunction, especially the DNA methylation. e different expression genes were also confirmed by Gene Expression Omnibus (GEO), NCBI's publicly available genomics database. To our knowledge, this is the first study that presents metabolic reprogramming plays a major role in the regulation of the epigenome across multiple cancer types. Our objective in this paper is to provide a global understanding of metabolic genes as well as their function in epigenetic modification in tumors and reveal potential targeted therapies for cancer.

Data Resources.
Samples for different cancer types were downloaded from the TCGA data portal. To eliminate heterogeneity from patients and reveal universality more accurately, the minimum threshold of normal samples was 30 for each cancer type. After the screening, 5137 cancer and 589 normal samples were analyzed. 11 cancer types were selected as follows: breast invasive carcinoma (BRCA), colon adenocarcinoma (COAD), colorectal adenocarcinoma (COADREAD), head and neck squamous cell carcinoma (HNSC), kidney renal clear cell carcinoma (KIRC), kidney renal papillary cell carcinoma (KIRP), liver hepatocellular carcinoma (LIHC), lung adenocarcinoma (LUAD), lung squamous cell carcinoma (LUSC), stomach adenocarcinoma (STAD), and thyroid carcinoma (THCA). For each cancer type, we classified samples into 2 groups: normal samples and tumor samples (Supplementary Tables 1 and 2).
For the validation material, we extracted the microarray expression profile (GEO group) from 2010 to 2019 in GEO database. Datasets were included if the following criteria were fulfilled: containing both tumor samples and adjacent normal tissues and normal sample size of more than 30. e five datasets (GSE13507, GSE87630, GSE37182, GSE76427, and GSE32665) containing adjacent normal tissues as the control were downloaded from the GEO repository. Details of each microarray study, including sample descriptions, are provided in Supplementary Table 1. For this study, we used 2071 human metabolism-related genes assigned to metabolic pathways in the Cancer Cell Metabolism Gene DataBase (ccmGDB), official gene symbols, and gene IDs are included in Supplementary Table 3.

Differential Gene Expression Analysis.
e mRNA expression profiles data were classified into 2 groups: normal and tumor. Gene expression was represented with TPM (transcripts per million) calculated by the counts of mapping reads in each gene. e gene with the mean TPM of less than 1 was excluded. e Cyber-T bayesreg.R, based on a Bayesregularized unpaired t test, was used to analyze differences between tumor samples and normal samples. Differentially expressed genes (DEGs) were selected by FDR (false discovery rate) less than 0.05 and together with the absolute value of fold change value not less than 1.5. According to the metabolic gene list, we filtered the differential expression for each metabolic gene in each cancer type. For the GEO data, limma R package was used to screen the DEGs between tumor and normal samples in each included dataset. We performed gene differential analysis (|LogFC| > 1, adjusted P value (FDR) < 0.05) as the cutoff criteria. Hierarchical clustering analysis was performed for the DEGs using the R packages ("pheatmap").

Functional Annotation.
e metabolic differential expression genes (MDEGs) that were found to differentially express at least 8 cancers were selected. For functional annotation of those metabolic genes, we used GO and KEGG pathway annotations within David bioinformatics database to perform the gene ontology analysis [18,19]. Significant pathways with a threshold of P < 0.05 were selected. e pathway enrichment bubble plot was drawn using R packages ("ggplot2").

Survival Analysis.
Kaplan-Meier analysis was used to identify the genes showing clinical relevance with tumor patients according to the TCGA datasets. e tumor samples were divided into two groups according to the median expression of each gene: high expression (with TPM values higher median) and low expression (with TPM values lower median). en, the log-rank test was used to analyze the differences between groups, and the significantly prognostic value was selected with a threshold of the P value <0.05. Univariate Cox regression analysis was used to further identify the independent prognosis factors. Statistical analyses were performed using R package ("survival").

DNA Methylation Analysis.
For DNA methylation analysis, the average DNA methylation value (β values) for all CpG sites correlated with a gene in the genomic region between 0 and 1,500 bps ahead of the transcription start site was calculated. e ChAMP is an R/Bioconductor package that was used to calculate differentially the significantly expressed probes that had a P value <0.05 and |Δβ| > 0.2 between two group samples, respectively, in our analysis (http://www. bioconductor.org/packages/release/bioc/html/ChAMP.html). One group was between normal and tumor, and the other group was the different DNA methylated level between top10% and bottom10% samples by DNMT3B expression. en, the overlap probes were obtained between the two groups. Differentially expressed probes with Δβ > 0.2 were defined as hypermethylation of the corresponding genes, and those with Δβ -0.2 were defined as hypomethylation of the corresponding genes. en, these probes were compared with different expression genes for the following analysis.

Spearman's Rank Correlation
Analysis. Spearman's rank correlation was calculated using the function "cor.test" in R between promoter probe methylation level and their corresponding gene expressions. e significant negative correlations were considered if the correlation coefficient r was <0, and significant positive correlations were considered if r was >0. We calculated the above overlapped genes and their corresponding probes. e genes were obtained whose Δβ > 0.2 and r < -0.2 were regarded as hypermethylated and significantly anticorrelated between methylation and expression results.

Association Analysis between Oncogene/Tumor Suppressor Mutation and Metabolic Pattern.
To identify the association between oncogene/tumor suppressor mutation events and metabolic profile, we first choose the top 15 mutation driver genes in pan-cancer according to the previous comprehensive study [20]. Information regarding mutation frequency and sample matrix in each cancer type was downloaded from the TCGA, an open-access database that is publicly available at cBioportal [21] (http://www.cbioportal. org). We choose the top 3 oncogene/tumor suppressor with a high mutation frequency for each cancer type and divided the sample into two parts according to whether each oncogene/ tumor suppressor is mutated or not. e differential gene expression analysis was conducted between the above grouping samples in each cancer type. e associations were analyzed between the focus metabolic gene of this study as well the whole metabolic spectrum and key oncogene/tumor suppressor mutation.

Global Changes in Metabolic Gene Expression.
To understand metabolic gene expression in different cancers, we analyzed 5137 tumor and 589 normal samples spanning 11 different tumor types. According to the clinical information, we divided samples into 2 subgroups for each cancer type: normal and tumor (Supplementary Table 1). Figure 1(a) presents the total number of expressed genes (EGs), differentially expressed genes (DEGs), and metabolic differential expression genes (MDEGs) in 11 cancers. e number of upregulation and downregulation of MDEGs in 11 cancers is shown in Figure 1(b). To validate the different analysis results in TCGA, we choose the MDEGs (found differential expression at least 8 cancers) to get the differential distribution in GEO database. As shown in Figure 1(c), the GEO data exhibited the same differential variation trend as in TCGA cancer types. To gain a view of metabolic status in cancer, we carried out the MDEGs (were differentially expressed in at least 8 cancers) sets of metabolic pathways based on the GO and KEGG pathway analysis (Supplementary Tables 3 and 4). e pathways, namely, purine metabolism, carbon metabolism, pyrimidine metabolism, and beta-alanine metabolism, were dysregulated in cancers ( Figure 2, Supplementary Table 4).

Disordered Metabolic Genes Related to Histone Acetylation.
In the transacetylation process ( Figure 3, Table 1), the citrate, acetate, and fatty acids were the main sources for generating acetyl coenzyme A (acetyl-CoA) which would provide the acetyl group to histone. In the citrate metabolism, the genes encoding pyruvate dehydrogenase (PDH) and ATP citrate lyase (ACLY) were included. PDH had no significant changes between cancer and normal samples. ACLY was upregulated in 5 cancer types. In the acetate metabolism, the gene-encoded acyl-CoA synthetase short-chain family member 2 (ACSS2) was downregulated in 5 cancer types. In the fatty acid synthetase part, ACACA and FASN were upregulated in 4 cancer types. In the fatty acid oxidation   process, the acylcarnitine would be transported to the mitochondria to mediate fatty acid metabolism by solute carrier family 25 member 20 (SLC25A20). SLC25A20 had a lower expression in 8 cancer types compared with those of normal samples. Citrate produced acetyl-CoA through the ACLY and acetate produced acetyl-CoA through ACSS2 also in the nucleus which indicated that acetyl-CoA existed in the nucleus and cytoplasm of mammalian cells, influencing both metabolism and the global regulation of the gene expression [5]. Histone acetyltransferases (HATs) could acetylate the acetyl group from acetyl-CoA to the histone. Lysine acetyltransferase 2A (KAT2A) was mainly studied as HATs. e gene KAT2A had an upregulated expression in 9 tumor types.

Disordered Metabolic Genes Related to DNA Methylation.
e disordered metabolic genes in transmethylation pathways were filtered out in our results. Accordingly, we assembled a metabolic map depicting the distribution of these changes in the pathway ( Figure 4, Table 1). ere were three metabolites as methyl donors: methionine, folate, and betaine. All of them provided the methyl group to S-adenosylmethionine (SAM) through their metabolic cycle. In the betaine part, the expression of SLC44A4, CHDH, and BHMT was disordered across cancers and had a general downward trend. In the folate cycle, SLC19A1 and MTHFD1 had upregulated expression in 6, 6 tumors and downregulated expression in 3, 1 tumor, respectively. DHFR and MTHFD2 were not downregulated in cancers, and they were upregulated in 6 and 9 cancers, respectively. Methionine and ATP were catalyzed into SAM by methionine adenosyl transferase (MAT). e gene MATIA was downregulated, and MAT2A was upregulated in one cancer. In the transmethylation reactions, the genes encoding histone methyltransferase (HMT) were not significantly different in cancers. Only DNA methyltransferase (DNMT) genes were presented in results: DNMT1, DNMT3A, and DNMT3B were all highly expressed in 8, 6, and 9 cancer samples compared with normal, respectively.

Survival Analysis of Disordered Metabolic Genes.
To analyze the DNA methylation and identify the related genes which may show clinical relevance with tumor patients, we performed the survival analysis on the disordered DNMTs genes using TCGA datasets. Disordered DNMT3B was screened and was found to have an impact on overall survival months. In BRCA, KIRC, KIRP, LIHC, and LUAD, patients whose tissues have a higher expression of DNMT3B  Figure 3: Detailed network map of transacetylation of gene expression levels in the 11 cancer types. Note: each metabolic reaction is marked with the number of tumors (out of 11 considered in our analysis) in which at least one gene in the corresponding reaction is significantly (FDR < 0.05) upregulated (red), downregulated (green), and neutral (black). If unmarked, no statistically significant change in mRNA expression was detected. Cytoplasmic citrate came from mitochondrial citrate which was produced by glucose for glycolysis in the cytoplasm and TCA cycling in the mitochondria. With the catalysis of ACLY, ACACA, and FASN, part of citrate was used to generate long-chain fatty acid which then would be oxidized in the mitochondria. Nuclear citrate and acetate were separately catalyzed by ACLY and ACSS2 to ac-CoA, which provided the acetyl groups to histones that led to histone acetylation. SLC2A1: solute carrier family 2 member 1; SLC25A20: had significantly shorter overall survival compared to those with lower expression. e higher expression of DNMT3B led to poor prognosis in those cancer types, especially in KIRP with the maximum hazard ratio ( Figure 5). e DNMT1 was significantly associated with patient survival times only in one cancer type (not shown), and it was seen that the DNMT3A expression did not affect patient survival. e red plots present the high expression of each individual while the blue plots present the low expression of each individual. P value of less than 0.05 was considered as statistically significant. HR: hazard ratio.

DNMT3B Contributing to the Altered DNA Methylation and Corresponding Gene Expression.
Based on the analysis the DNA methylation disorders resulted from DNMT3B, different DNA methylation probes were obtained in high and low DNMT3B expression samples and then compared to the difference between tumor and normal. e number of altered methylation probes was varied from 26 to 234 across tumor types, and these changes may be contributed by high expression of DNMT3B (Supplementary Table 6). ese probes correspond to genes that were enriched in the different GO terms depending on the cancer types, such as regulation of cell differentiation in KIRC and signal transduction in KIRP (Supplementary Table 7). Since the gene expression level was contributed by DNA methylation alteration in tumors, we overlapped the differential methylation and the differential expression genes data in cancer versus normal samples. 8 to 87 genes were screened out and were dysregulated with the differential methylation and mRNA across cancer types but few of them overlapped (Figure 6(a)). Functional annotation showed that the most frequently hypermethylated genes (Δβ > 0.2) across cancers were enriched in anterior/posterior pattern specification (Figure 6(b)), and the hypomethylated genes were in signal transduction (Figure 6(c)). Based on the calculated correlations between methylation and expression data, the highly anticorrelated genes were retained. After filtering and merging, thereinto, the hypermethylated-repressed and hypomethylated-activated genes are shown in Table 2 (Supplementary Table 8).

Association between Metabolic Gene Expression and
Oncogene or Tumor Suppressor Mutation. Metabolic reprogramming can be partly influenced by oncogenic driver events. To identify somatic alterations that potentially drive metabolic expression, we identified the top 3 oncogene/tumor suppressor with the highest frequency of mutations and assessed whether their mutation status correlated  with the metabolic profile. THCA and KIRP were screened out because sample size with each 3 gene mutation was less than 20 (Supplementary Table 9). For each cancer type, we performed the differential expression analysis of metabolic genes based on the two group samples with oncogene/tumor suppressor mutated or not (Supplementary Table 10). DNMT3B was upregulated in BRCA, LIHC, and STAD in TP53-mutated conditions when compared against nonmutated samples. Upregulated expression SLC2A1 was observed in LIHC and LUAD in TP53-mutated conditions. e other key metabolic genes related to DNA methylation and histone acetylation described in Table 1 were not associated with oncogenes or tumor suppressors within the scope of our study (Supplementary Table 10).
We also analyzed the association between metabolic pathways and mutation patterns. Compared to other cancer types, most of metabolism pathways in BRCA were disordered with TP53 and MYC mutation, respectively. In BRCA, LIHC, LUAD, HNSC, and STAD, groups with TP53 mutations were associated with most of the genes related to nucleotide, lipid, and amino acid metabolisms. e other metabolic process associated with mutation patterns across cancer types is shown in Supplementary Figure 1

Discussion
Dysregulation of metabolism and dysfunction of epigenetic modification is now the established feature of cancer. However, the mechanism of how metabolic gene rewriting influenced the global epigenetic modification levels in pancancer is not yet fully clear. We performed a pan-cancer analysis of transcriptome profiles covering 11 tumor types in a comprehensive collection of 2071 metabolic genes. We also elaborated on the altered metabolic genes which are involved in the dysfunction of epigenetic modifications including acetylation and methylation across cancers. ACLY, SLC2A1, KAT2A, and DNMT3B are the critical genes contributing to epigenetic dysfunction across cancers. ereinto, DNMT3B is upregulated in 9 cancer types, and its higher expression showed a correlation with overall survival in 5 cancer types and it might contribute to the specific DNA methylation disorder. is result suggests a close link between DNMT3B expression, DNA methylation, and carcinogenesis.   Various metabolism pathways responsible for epigenetic modifications have been identified. One relevant pathway of metabolic gene rewiring and epigenetic landscape is the histone acetylation, which plays a critical role in regulating chromatin structure, thus participating in specific gene regulation [22]. A number of researches show that unstable histone acetylation is contributed by changed acetyl-group donor acetyl-CoA [23], which is a key metabolite in the mitochondria and cytoplasm, associated with breakdown of carbohydrates and fats via altered glycolysis and β-oxidation, respectively [24]. e acetyl-CoA synthesis is by ACSS2 [25] and ACLY [5] according to the different substrates. ere is a similar or opposite trend for key gene expression, ACSS2 is decreased in 5 cancers, and ALCY is found in 5 tumors compared to their normal samples. We consider that the citrate to the acyl-CoA pathway plays an important role in histone acetylation in the carcinogenesis of categories we studied. ACLY plays a critical role in determining the histone acetylation, and ACLY knockdown leads to apoptosis and growth suppression in breast cancer cells [26]. Inhibition of ACLY inhibits glucose-dependent histone acetylation to suppress glioblastoma cell proliferation [27]. ACLY was upregulated in COAD, COADREAD, HNSC, KIRC, and LIHC in the present study which indicates ACLY might promote tumor proliferation in these cancer types. Interestingly, the reduction in histone acetylation and defect of SLC2A1 expression are observed in the ACLY knockdown cells [5]. erefore, upregulated ACLY and SLC2A1 in COAD, COADREAD, HNSC, KIRC, and LIHC might result in histone acetylation leading to selective regulation of genes involved in glucose metabolism. us, we speculate ACLY and SLC2A1 may serve as critical linker of metabolic and histone acetylation in these cancer types. For metabolite analysis, lysine acetyltransferase 2A (KAT2A) is an affirmed histone acetyltransferase (HAT) that binds to acetyl-CoA and transfers the acetyl group to histones [28]. As expected, the results of our data show that the KAT2A gene has high expression in 9 cancer types, and it may be explained that KAT2A is important for controlling the gene expression program for adjusting histone acetylation in these cancer types. KAT2A over expression increases the extent of histone acetylation by interacting with E2F1, cyclin D1, and E1 promoters to promote the proliferation of lung cancer cell lines [29]. KAT2A knockdown leads to a significant reduction of DNA synthesis in cervical cancer cells by decreasing histone H3 acetylation in the E2F1 promoter [30]. Our data show that the KAT2A has higher expression in 9 cancer types which indicates the increased histone acetylation level in specific gene promoters among these cancers to be involved in the development of cancer. Further study is needed to investigate ACLY, SLC2A1, and KAT2A in the histone acetylation function across cancer types.
Another epigenetic modification is the DNA methylation, changes in transmethylation influenced by altered methionine, choline, and folate. Unstable methylation usually relates to altered methyl donor [31]. In transmethylation reactions, methyl donor is finally provided by SAM which is produced from methionine by MAT [32]. Increased adenosylhomocysteinase (AHCY) gene expression will lead to an increase in SAM and aberrant DNA methylation [33]. AHCY knockdown leads to SAM depletion and proliferation rate reduction in hepatocellular carcinoma cells [34]. SAM is the methyl donor, and suppression of SAM leads to decreased DNA methylation and slows the growth of pancreatic cancer cells [35]. Upregulated AHCY gene expression is seen in BRCA, COAD, COADREAD, LIHC, LUAD, LUSC, and STAD that may indicate the transmethylation reactions have a more active state in these cancer types. e altered expression of DNMTs, encoding DNMTs (DNMT1, DNMT3a, and DNMT3B) could result in global changes of methylation in leukemia cells and breast cancer stem-like cells [36,37], while DNMT1, DNMT3A, and DNMT3B are all upregulated in most cancer types (8, 6, and 9, respectively) in this study. Overexpressed DNMT3A contributes to the various methylation patterns and is a consequence of AML progression [37]. DNMT1 promotes hypermethylation and downregulation of tumor suppressor gene ISL1 which increases the tumor stem cell population in breast cancer cells [38]. However, only DNMT3B higher expression is correlated with poorer overall survival in BRCA, KIRP, KIRC, LIHC, and LUAD in the present study. We then analyzed the DNA methylation level in 5 cancer types distinguished by high and low expression of DNMT3B. e relation between DNMT3B expression and methylation sites varies in cancer. Studies have shown that DNMT3B can induce DNA methylation in specific CpG islands in colorectal cancer [39] and it could also induce the distinct methylation level in different regions such as CpG and non-CpG [40]. Both hypermethylation and hypomethylation in gene promoters between high and low DNMT3B expression groups are found across cancers, which is consistent with the previous research [40]. We then choose the overlap genes which harbor differential methylation and expression; meanwhile, many of these genes are involved in different functional annotation that depends on the tissue types, suggesting that DNA methylation at different loci happens through DNMT3B in different cancers. Among these genes, we found numerous genes whose function in tumorigenesis and their expression associated with promoter hypermethylation or hypomethylation were well-documented in the literature. ZNF154 and AQP1, whose hypermethylated patterns are biomarkers for distinguishing tumor from normal samples [41,42], are downregulated and hypermethylated in BRCA and LUAD in our study. In addition, hypermethylation and downregulation of PPP2R2B, whose mRNA expression is in a DNMT-dependent manner [43], were found in BRCA. Loss of DNMT3B in the mouse model delayed the melanoma formation suggesting that other DNMTs do not adequately compensate for DNMT3B loss and suggesting nonredundant roles for DNMTs in melanoma and DNMT3B, in particular, may play specific, nonredundant roles [44,45]. Inhibitions of DNMT3B as a novel therapeutic strategy for inhibiting the proliferation of carcinoma cells are being studied [43]. Together with these findings, it is suggested that DNMT3B with encoding protein DNMT3B and its potential methylation substrate genes would be the effective therapy targets to human cancer.
To understand the factors affecting the expression of metabolic genes at the present study, we have analyzed the association between oncogene/tumor suppressor mutation and metabolic gene expression in each cancer type. Several metabolic pathways, such as nucleotide, lipid, and amino acid metabolisms, are associated with specific mutation patterns across cancer types indicating that metabolic reprogramming may partly result from diverse oncogene/ tumor suppressor alterations in different tumor contexts, which is also shown in the previous study [16]. Except for SLC2A1 and DNMT3B, other genes related to histone acetylation and DNA methylation mentioned above were not associated with oncogene/tumor suppressor mutation at the present study. ese results suggest that there might be other factors affecting the expression of metabolic genes, like TME, rather than just oncogene/tumor suppressor mutation. e expression of ACLY and ACSS2 is influenced by glucose concentration [27] and hypoxic conditions [46] in local microenvironment. TME modulates cancer cell metabolism and epigenetic modification contributing to tumor heterogeneity and therapeutic response through limited nutrient supply, acidic, hypoxia, and other characteristics [47][48][49]. erefore, it is necessary to combine metabolic phenotype, epigenetic modification, genomic landscape, and local TME to reveal the relationship between metabolic genes and epigenetics modification in pan-cancer, which may help overcome therapy resistance in cancer patients and guide the rational design of combinational therapies targeting both tumor cells and microenvironment components.

Conclusion
In summary, our research clarifies some key metabolic genes, ACLY, SLC2A1, KAT2A, and DNMT3B, which are involved in epigenetic regulation indirectly in tumors. e uniqueness of our study is that it is first report that presents metabolic gene reprogramming as a factor influencing epigenetic regulation across human pan-cancers. We found some potential genes' methylation and expression in different cancer types influenced by DNMT3B whose expression is strongly relevant to the patient overall survival. Besides, the connection was identified linking oncogenes and tumor suppressors with deregulated cancer cell metabolism. e metabolic phenotype affected by TME should also be gotten attention. us, targeting these tumor metabolic genes might reverse epigenetic dysregulation and improve the effect of targeted therapy. Moreover, the combinational therapies targeting metabolic phenotype and epigenetic modification with genomic landscape and local TME can more effectively eradicate tumor cells from the patients and improve the quality of life. Expressed genes DEGs:

Abbreviations
Differentially expressed genes MDEGs: Metabolic differential expression genes FDR: False discovery rate

Supplementary Materials
Supplementary Table 1: the number of sample information and differential expression genes for 11 different cancer types and 5 GEO data projects. Supplementary Table 2: the differential expression gene list and its fold change in 11 cancer types compared with normal samples. Supplementary Table 3: the metabolic genes and the differential expression of metabolic genes. Supplementary  Table 7: gene ontology analysis of differential expression and methylation genes. Supplementary Table 8: the correlation between expression and methylation of specific genes. Supplementary Table 9: the oncogene/tumor suppressor mutation frequency and the sample size with mutation across cancer types. Supplementary Table 10: the differential expression MDEGs in mutated samples compared with non-mutated samples across cancer types. Supplementary Table 11: the pathway enrichment analysis in each mutation group. Supplementary  Figure 1: relationships among metabolic pathways, cancer types, and oncogene/suppressor gene mutations. (Supplementary Materials)