Identification of DNA methylation patterns and biomarkers for clear-cell renal cell carcinoma by multi-omics data analysis

Workflow chart and samples preparation

Figure 1 shows the study workflow chart. In this study, DNAm (290 tumor and 160 normal), mRNA (526 tumor and 72 normal) and clinical data (318 cases of ccRCC) were collected from TCGA (https://portal.gdc.cancer.gov/). Of note, 160 cases of ccRCC had both ccRCC tissues and matched adjacent normal tissues in terms of DNAm profiles. In addition, 24 cases of ccRCC, which had both ccRCC tissues and matched adjacent normal tissues in terms of mRNA profiles, had both DNAm and mRNA expression profiles. The human KEGG Markup Language (KGML) files were obtained from the KEGG database (http://www.genome.jp/kegg/), and these files were used to predict key enzyme-coding genes. In addition, the DNAm status for the following five additional cancers were derived from TCGA, including both DNAm and gene expression profiles: BRCA (774 tumor and 82 normal), COAD (292 tumor and 38 normal), LIHC (380 tumor and 50 normal), LUAD (455 tumor and 32 normal) and PRAD (502 tumor, 50 normal). Moreover, the following two cancers were obtained from TCGA, including only DNAm or gene expression profiles: KIRP (274 tumor and 45 normal) and UCEC (425 tumor and 34 normal).

Differential analysis of mRNA expression and DNAm

Based on DNAm profiles, probes with more than 50% missing values were deleted. The “impute.knn” function (R version 3.5.0) was used to impute the remaining missing values. CpG probes in chromosomes X and Y were excluded, and CpG probes identified as non-unique in the genome as well as CpG probes known as single-nucleotide polymorphisms were also excluded (Noushmehr et al., 2010; Selamat et al., 2012). Differentially methylated CpG sites (CpGs) were detected using the non-parametric testing method (Selamat et al., 2012), and p-values were adjusted using the false discovery rate (FDR) method in R (R version 3.5.0; p.adjust (method = “BH”)). CpG probes were filtered as described below, and these filtered CpG probes were considered to be statistically significantly differentially methylated. FDR-values were less than 0.05, and the absolute values of the mean β-value difference between ccRCC and non-ccRCC tissues were greater than 0.2 (Selamat et al., 2012). In addition, the differentially expressed genes were detected based on “limma” function in Bioconductor 3.10. The significantly differentially expressed genes were identified using the following cutoffs: FDR-value < 0.05 and |log2 fold change| > 1.

Finally, we investigated the correlation between DNAm and mRNA expression. Pearson correlation coefficients for each gene or CpG were calculated between mRNA expression and DNAm using the Spearman’s rank correlation coefficient. The correlation was considered statistically significant when the absolute value of Spearman’s Rho was greater than 0.2 and the FDR-values were less than 0.05 (Cancer Genome Atlas Research, 2011; Cheng et al., 2018).

Screening key enzyme-coding genes by met-express

In previous study, the gene co-expression and genome-scale metabolic network were integrated to identify the key enzyme-coding genes in human diseases. Here, we briefly describe the algorithm (Chen et al., 2013). (1) Based on transcriptome data, the Pearson correlation coefficients (PCC) were calculated for all pairs of enzyme-coding genes. The top three associated genes with the highest PCC for each gene were used to construct gene co-expression network. Qcut (Ruan & Zhang, 2008) was used to divide the network into gene co-expression modules, and we only remain some modules including more than 10 genes. (2) We carried out the reconfiguration of metabolic network based on KGML files (http://www.genome.jp/kegg/) containing the information of critical metabolic reactions, key enzyme-coding genes and corresponding metabolites. (3) The importance score was calculated for each gene by integrating with gene cox-expression network and reconstructed metabolic network. Then, the median importance score was used for a threshold to select the candidate key enzyme-coding genes.

Identification of candidate tissue-specific diagnostic biomarkers

For 24 patients who have both mRNA expression and DNAm profiles, first-level feature selection was performed though the following steps. First, DNAm profiles were obtained by Wilcoxon rank-sum test to select differentially methylated CpGs (FDR < 0.05, |Δβ_mean| > 0.2); mRNA expression profiles were obtained by “limma” function (R version 3.5.0) to select differentially expressed genes (FDR < 0.05, |log2 fold change| > 1). Second, the Rho was calculated between DNAm and mRNA expression by Spearman’s rank correlation coefficient to select genes with higher correlation coefficient (FDR < 0.05, |Rho| > 0.2). For genes with multiple CpGs measuring Rho, we selected the CpG with the highest absolute value of Rho. Third, the key enzyme-coding genes in cancers were selected by Met-Express (Chen et al., 2013). Finally, CpGs obtained from the above three steps were intersected to obtain the candidate CpGs, namely the candidate tissue-specific diagnostic biomarkers for first-level feature selection.

To build the effective diagnostic model, it was essential to identify a small set of biomarkers. Thus, some redundant biomarkers were excluded to avoid over-fitting and achieve the best prediction performance (Tang et al., 2018). After the first-level feature selection, the number of those CpGs was still large and redundant. Therefore, we performed a second-level feature selection by the LASSO method (Tibshirani, 1996) to further remove redundant CpGs and identify tissue-specific diagnostic biomarkers.

Evaluation and modeling of candidate tissue-specific diagnostic biomarkers

In this study, the random forest (RF) classifier was built based on the tissue-specific diagnostic biomarkers to evaluate performance of those biomarkers. Apart from 24 patients for the variable selection, 136 patients who have both tumor and adjacent normal tissue were randomly allocated to the training and testing sets. The training set containing 90 patients was used to optimize the parameters of classifiers and build classifiers. The testing set containing 46 patients was used to evaluate performance of the classifiers. Ten-fold cross-validation was used during the training and testing of classifiers, and the performance of classifier was evaluated based on area under the ROC Curve (AUC). Meanwhile, we also introduced the sensitivity and specificity to assist in evaluating the performance of the diagnostic classifier.

Validation based on independent data sets and literatures

In this study, the Illumina HumanMethylation450 BeadChip data was comprehensive queried by using keywords focused on “renal clear cell carcinoma”, “ccRCC” and “methylation”. In the process of data retrieval, the following conditions should be met as much as possible which included the initial experimental research providing a comparison of ccRCC and non-ccRCC tissues, the patients without receive neoadjuvant therapy. We obtained two independent datasets, including GSE70303 (Becket et al., 2016) (46 tumor and non-tumor tissues from GEO database) and E-MTAB-2007 (Sato et al., 2013) (106 tumor and non-tumor tissues from ArrayExpress database), to perform verification of candidate tissue-specific biomarkers.

Construction of prognostic model

Apart from 24 patients for the variable selection, we split the rest of ccRCC patients in a training dataset and a validation dataset and explored to build a predictive model for prognostic and survival analysis. We used 192 cases for training and 98 cases for validation. We applied a sequential model-based variable selection to screen markers for predicting survival outcome (Ding, Chen & Shi, 2019; Xu et al., 2017). Based on the candidate biomarkers, we first fitted a univariate Cox proportional hazards model by using each marker as the covariate. A marker with p-value < 0.05 from the Wald statistic was retained in the dataset. Second, we used a similar strategy in a prognosis marker selecting process based on the LASSO-Cox method to decrease the marker numbers to a reasonable range (less than events). The above analysis generated final markers with non-zero coefficient to construct a prognostic signature using R package “glmnet”. Based on the penalized Cox regression model, we obtained a combined prognostic score (designated as risk score) for each individual. The risk-score for patients is calculated as follows: $$\mathrm{R}\mathrm{i}\mathrm{s}\mathrm{k}\mathrm{S}\mathrm{c}\mathrm{o}\mathrm{r}\mathrm{e}\mathrm{f}\mathrm{o}\mathrm{r}\mathrm{p}\mathrm{a}\mathrm{t}\mathrm{i}\mathrm{e}\mathrm{n}\mathrm{t}\mathrm{s}={\sum }_{\mathrm{i}=1}^{\mathrm{N}}(\mathrm{c}\mathrm{o}\mathrm{f}\mathrm{f}\mathrm{i}\mathrm{c}\mathrm{i}\mathrm{e}\mathrm{n}\mathrm{t}\mathrm{s}\mathrm{o}\mathrm{f}\mathrm{e}\mathrm{a}\mathrm{c}\mathrm{h}\mathrm{C}\mathrm{p}\mathrm{G}\mathrm{s}\times \beta -\mathrm{v}\mathrm{a}\mathrm{l}\mathrm{u}\mathrm{e}\mathrm{s}\mathrm{o}\mathrm{f}\mathrm{e}\mathrm{a}\mathrm{c}\mathrm{h}\mathrm{C}\mathrm{p}\mathrm{G}\mathrm{s})$$

To validate the predictive classifier, we calculated a risk-score for each patient in the validation dataset using the coefficient estimated from the training dataset. By grouping the risk-score based on its median, we formed high and low risk-score groups with roughly equal number of observations. The Kaplan–Meier estimator and log-rank test were performed to test whether the median survival time for two groups was significantly different. All analyses were conducted in R version 3.5.0 with the “glmnet” and “survival” packages.

Functional enrichment analysis based on gene ontology

The CpG and gene annotations were downloaded from TCGA (http://tcga-data.nci.nih.gov/tcga/tcgaPlatformDesign.jsp) and GECODE (https://www.encodeproject.org/files/gencode.v22.annotation/). Genes with abnormal methylation were analyzed for Gene Ontology (GO) analysis or Kyoto Encyclopedia of Genes (KEGG) pathway enrichment analysis by using the web version of DAVID (Da Huang, Sherman & Lempicki, 2009). GO/KEGG-terms have an adjusted p-values which were less than 0.05. The genes were uploaded to the online database GeneMANIA (http://genemania.org; version 3.6) (Warde-Farley et al., 2010) to explore their interactions at the protein level with an interaction score > 0.4 as the cutoff value. Afterward, the network was visualized in software Cytoscape (version 3.6.1).

Statistical analysis

Based on Euclidean distance matrix and complete-linkage method, unsupervised hierarchical clustering was carried out to verify the reliability of the analysis method about CpGs in R version 3.5.0 with the “heatmap” and “dist” functions. To discuss the widespread nature of cancer-specific methylation alterations and the impact of aberrant methylation on gene expression, the Chi-square test and Wilcoxon rank-sum test were used in R version 3.5.0 with the “chisq.test” and “wilcox.text” functions. When p-values were less than 0.05, it was deemed to be statistically significant. The diagnostic models were built in scikit-learn framework (version 0.20.3; Python 3.6..8) with the “sklearn.ensemble.RandomForestClassifier” algorithm. The parameters of models was tuned in scikit-learn framework (version 0.20.3; Python 3.6.8) with the “sklearn.grid_search.GridSearchCV” algorithm. The hyper-parameters included the number of trees (10 to 100), the criterion (“gini” and “entropy”), “oob_score” in the RF. The default values were used when other parameters were not provided in functions or classifier.

Identification and landscape of differentially methylation in ccRCC

Based on the differentially methylated CpGs, the exploratory hierarchical clustering was first performed (Fig. S1). The methylation profiles of ccRCC and non-ccRCC resulted in separate groups, suggesting a potential difference. In our study, we obtained 4367 CpGs (corresponding to 2,490 genes) that were significantly hypermethylated in ccRCC and 5,031 CpGs (corresponding to 3,551 genes) that were significantly hypomethylated (Fig. 2A). Some of the most hypermethylated loci contained HOXA5, involving 31 CpGs in the CGI (29/31) or N_Shore (2/31) and coding for a DNA-binding transcription factor, which regulated gene expressions, cell differentiation and morphogenesis (Shenoy et al., 2015). Some of the most hypomethylated loci contained CXCL8, which promoted angiogenesis and metastasis (Shenoy et al., 2015). To investigate the biological functions of genes with abnormal methylation, the target genes was conducted by DAVID (Da Huang, Sherman & Lempicki, 2009). Most genes were significantly enriched in the important biological processes, including cell adhesion, signal transduction and positive regulation of transcription (Benjamini–Hochberg (BH) adjusted p-value < 0.05) (Table 1), indicating that these genes played an important role in tumorigenesis.

For the significantly differentially methylated CpGs, we also analyzed whether they are related with CGIs and whether they are located in the promoter or gene body (Fig. S2). We identified 9398 CpGs defined as being differentially hyper- or hypomethylated. Moreover, 33.33% and 31.96% of CpGs located in promoter regions, which was defined as 1.5 kb up-or downstream from the transcription start site (TSS), were significantly hyper-or hypomethylated, respectively. We found no statistically significant difference between the degree of methylation and whether the CpGs were located in the gene promoter regions (p-value = 0.17, Fisher’s exact test). However, 31.87% and 3.79% of CpGs located in CGIs were significantly hyper-or hypomethylated (p-value < 2.2 × 10⁻¹⁶, Fisher’s exact test), respectively. These results showed a significant difference between the methylation characteristics of CGIs and non-CGIs for differentially methylated CpGs. Based on CGI regions, the significantly hypermethylated CpGs were highly enriched in CGIs (p-value < 2.2 × 10⁻¹⁶, χ² test), whereas hypomethylated CpGs were enriched in non-CGIs (p-value < 2.2 × 10⁻¹⁶, χ² test) (Fig. 2B). Compared to hypomethylated CpGs, hypermethylated CpGs were located closer to TSSs (p-value = 2.37 × 10⁻⁹, Wilcoxon rank sum test) (Fig. 2C). Compared to CGIs containing hypomethylated CpGs, CGIs containing hypermethylated CpGs were located closer to TSSs (p-value = 2.6 × 10⁻⁶, Wilcoxon rank sum test) (Fig. 2D).

In conclusion, the hypermethylated CpGs were preferentially enriched in CGIs, but CpGs exhibiting abnormal hypomethylation in tumors were frequently occurred in non-CGI regions. Furthermore, the CGI and promoter regions containing hypermethylated CpGs were located closer to TSSs.

Impact of the cancer-associated differentially methylation no gene expression

To identify the tumor-associated methylation alterations with concomitant changes at the transcriptional levels, the correlation analysis of the methylation and transcriptome data were performed. By focusing on CpGs and their nearest neighbor genes (called “cis-interactions”) that fall within 1,500 bp upstream of TSS to the end of gene body, genes with abnormal methylation were studied (Fig. S3). By comparing 1,154 significantly differentially methylated CpGs located in 316 associated genes, 998 (86.48%, 998 of 1,154 CpGs) CpGs with the significant correlations between DNAm and gene expression were identified.

In accordance with DNAm mechanism to silence or enhance local transcription, 68% of CpGs identified as significantly correlated were negative (681of 998 CpGs) (Figs. 3A and 3B). A portion of the CpGs were primarily in gene promoter regions (68.14%, 464 of 681 CpGs). Of these, 336 CpGs (72.41%, 336 of 464 CpGs) were statistically significantly hypermethylated and inhibited gene expression, while 128 CpGs (27.59%, 128 of 464 CpGs) were hypomethylated and caused over-expression of genes, indicating that abnormal DNAm might have functional consequence in approximately 40% (464 of 1,154 CpGs) of CpGs associated with genes. The hypermethylated and under-expressed genes included Secreted Fzd-related protein 1 (SFRP1), which was consistent with previous ccRCC studies (Awakura et al., 2008; Gumz et al., 2007), and PIC3K2G, which was a known cancer susceptibility gene (Brown et al., 2015) (Table S2). The WNT antagonist, SFRP1, was also a tumor suppressor that regulates cell proliferation and blood vessel formation (Atschekzei et al., 2012). CGI methylation located in SFRP1 has also been shown to be indirectly associated with tumor progression and tumor stage (Urakami et al., 2006) as well as to directly lead to the occurrence of poor prognosis in RCC patients (Morris et al., 2010; Saini et al., 2009). PIK3C2G encoded a phosphoinositide three kinase (PI3K) subtype, which was involved in the regulation of signaling pathways, such as cell proliferation, oncogenic transformation, cell survival and cell migration (Brown et al., 2015; Freitag et al., 2015). Conversely, the hypomethylated and over-expressed genes included single transduction-related genes (GUCA2B) as well as functional genes involved in cell differentiation (PPDPFL), PPAR signaling pathway (FABP6), HIF-1-alpha transcription factor network (CA9), and receptor-associated protein activity (SLC6A3) (Table S2). The molecular biological functions uncovered that hypermethylated/under-expressed genes were associated with the immune response, inflammatory response and cell adhesion (Table S3). However, hypomethylated/over-expressed genes were associated with cell or tissue development, intracellular transport and cellular secretion (Table S3).

In our investigation, a certain proportion (31.76%, 317 of 998 CpGs) of positive correlation was observed (Fig. 3B). Of which a lower percentage (10.73%, 34 of 317 CpGs) of the positively correlated CpGs was located in promoter regions, but our analysis indicated that these CpGs might be related to atypical TSS distances compared to negatively correlated CpGs. Specifically, the hypermethylated and over-expressed CpGs tended to lie closer to TSSs (p-value = 1.4 × 10⁻⁵, Wilcoxon rank sum test), whereas the hypomethylated and under-expressed CpGs were located further upstream of TSSs compared to their cognate negatively correlated (hypomethylated and over-expressed) CpGs (p-value = 1.4 × 10⁻², Wilcoxon rank sum test) (Figs. 3C and 3E). The molecular biological functions revealed that the hypermethylated CpGs with positive correlation were mainly enriched in immune processes and cell motility, whereas the hypomethylated CpGs with positive correlation were enriched in cell proliferation and development. Furthermore, CpGs that appeared to be positively correlated were inclined to be located in gene bodies (p-value < 2.2 × 10⁻¹⁶, χ² test) (Fig. 3D). In the gene bodies, the hypermethylated and over-expressed genes included RUNX3, and hypomethylated and under-expressed genes included FBXO2 (Figs. 3F–3K). Intriguingly, several genes showing positive correlation in gene bodies have been previously linked to cancer, including MUC15 (Huang et al., 2009), HEPACAM2 (Klopfleisch et al., 2010), CA10 (Romeo et al., 2009), NRG1 (Huang et al., 2004) and RAB25 (Mitra, Cheng & Mills, 2012). In addition, ~26% of genes that exhibited positive correlation and were located in gene bodies also exhibited negative correlations in promoters (Table S4), including RUNX3 and TMEM30B (Fig. S4). Based on four independent datasets including two methylation datasets (GSE70303 and GSE61441) and two gene expression profiles (GSE40435 and GSE76351), we also found the same events (Table S5). These results suggested that epigenetic alterations had a dual effect on gene expression, that is, genes with hypermethylation in the promoter and hypomethylation in the gene body were suppressed or genes with hypomethylation in the promoter and hypermethylation in the gene body were activated. The occurrence of the “tandem control” mechanisms indicated that methylation alterations located in the promoter might interact with methylation status in the gene body, thereby regulating the expression of cancer-associated genes.

Identification of tissue-specific biomarkers based on multi-omics integration

To obtain cancer tissue-specific CpGs that would distinguish tumor from non-tumor tissues, the correlation analysis was used for the selection of highly tissue-specific biomarkers using 24 tumor samples with both DNAm and matched gene expression profiles (Fig. 1). Based on correlation analysis, 134 CpGs associated with 54 genes were obtained by integrating gene expression, DNAm, enzyme gene networks and gene co-expression networks. For genes with multiple CpGs, the CpGs with the highest correlation between DNAm and gene expression were selected (Table S6). The optimal number of features was determined by the automated searching algorithm of LASSO. Classifiers with higher complexity would contain redundant features with no corresponding improvement in performance, thereby resulting in poor generalization ability. Eventually, the five ccRCC-specific CpGs (Table 2) were obtained using the LASSO algorithm, and they were significantly differentially methylated between the tumor and matched adjacent normal tissues (p-value < 0.001; Wilcoxon rank sum test) (Fig. 4A). Unsupervised hierarchical clustering was further performed in the independent datasets, and the similarity matrix was constructed using Pearson’s correlation coefficient (Fig. 4B). Tumor tissues were well discriminated from the matched adjacent normal tissues in independent datasets, and an accuracy of 100% was obtained. The results showed that five CpGs had better robustness and the potential to distinguish cancer.

Functional enrichment and biological network analysis

To further elucidate the function of the above five CpGs significantly associated with gene expression, the molecular biological function analysis were performed by DAVID (Da Huang, Sherman & Lempicki, 2009) and GeneMANIA 3.6 (Warde-Farley et al., 2010). As shown in Fig. 5A (Table S7), the genes associated with five CpGs were mainly enriched in Morphine addiction, Purine metabolism, Pentose phosphate pathway, Galactose metabolism and Fructose and mannose metabolism. In general, the reprograming of cellular metabolism is a major hallmark of cancer cells, which is a direct and indirect consequence of carcinogenic factors. The common characteristic of cancer cell metabolism is the ability to access adequate nutrients from its surrounding environment and make good use of these nutrients to maintain viability or proliferation (Fouad & Aanei, 2017; Pavlova & Thompson, 2016). Two major nutrients utilized for survival and biosynthesis in cancer cells are glucose and glutamine. Interestingly, PFKP and NAT8 were significantly associated with Glycolysis/Gluconeogenesis, Fructose and mannose metabolism, Galactose metabolism and Glutathione metabolism (p-value < 0.05, Fisher’s exact test). In addition, PDE1A and PDE8B were significantly enriched in Purine metabolism and G protein-coupled receptor signaling pathway (p-value < 0.05, Fisher’s exact test). Apart from playing an indispensable role in synthesizing DNA and RNA, purine metabolites supply some cancer cells with necessary nutrients and cofactors for survival or proliferation (Pedley & Benkovic, 2017; Yin et al., 2018) and they are also being targeted for the treatment of cancers (Fridley et al., 2011). Moreover, the majority of biological responses in both normal and cancer cells are regulated via multiple signaling pathways by G protein-coupled receptors (GPCRs) (Ristic, Bhutia & Ganapathy, 2017). Some research has suggested that GPCRs played an important role in purinosome assembly/disassembly and regulated metabolic flux through de novo purine biosynthesis in cells (Fang et al., 2013; Verrier et al., 2011). Furthermore, the ALOX5 gene was associated with several physiological and pathological processes, including inflammation, oxidative response and carcinogenesis (Tomita et al., 2019). Protein–protein interaction networks from GeneMANIA 3.6 suggested a more complex interaction for PFKP, PDE1A and ALOX5 (Fig. 5B). Most of these genes, including SREBF1 (Audet-Walsh et al., 2018), CALM1 (Fortney et al., 2015), PCBP1 (Zhang et al., 2015), COTL1 (Guo et al., 2017) and BTK (Yue et al., 2017), are tumor-related genes, which have been reported to be associated with tumorigenesis, progression and therapy.

Evaluation and validation of diagnostic accuracy with 10-fold cross-validation

Based on the five CpGs, a RF algorithm was implemented to construct an effectively diagnostic classifier. We obtained a sensitivity of 100% and specificity of 100% for ccRCC in the training dataset of 90 ccRCC and control samples (Fig. 6A) as well as a sensitivity of 97.83% and a specificity of 100% in the testing dataset of 46 ccRCC and control samples (Fig. 6B). Furthermore, the diagnostic classifier was demonstrated to effectively discriminate ccRCC from normal samples both in the training dataset (AUC = 99.99%) and the testing dataset (AUC = 99.95%) (Figs. 6C and 6D). To evaluate the performance of the five-CpG-based classifier, comparison of different independent datasets was required. In our research, two independent datasets, namely GSE70303 (Becket et al., 2016) and E-MTAB-2007 (Sato et al., 2013), were used with AUCs for 100% and 99.68% (Figs. 6E and 6F), respectively. To explore and verify universal applicability of the feature selection method and classifiers, the five most common cancers from TCGA were used to construct diagnostic classifiers, and the independent datasets from GEO were used to evaluate performance (GSE69914, GSE48684, GSE66836, GSE76938 and GSE54503). Remarkably, the AUCs based on the selected cancers were 97.31% (BRCA), 92.68% (COAD), 98.03% (LUAD), 94.90% (PRAD) and 98.32% (LIHC) (Fig. S5). The results showed that screening and identification of biomarkers based on multi-omics analysis could be extended to the application for any other cancers.

To further demonstrate the diagnostic capacity of the five-CpG-based classifier in the early stages, the diagnostic classifier was implemented and evaluated at different cancer stages (stage I–IV). Firstly, the unsupervised hierarchical clustering was implemented in the independent datasets, including four different cancer stages. Strikingly, these patients were divided into two different groups, which separately contained tumor and non-tumor tissues (Fig. 7A). Subsequently, the five-CpG-based classifier was used to distinguish ccRCC from normal tissues at different stages. Remarkably, the five-CpG-based classifier achieved ROC curves with an AUCs of 99.60% (stage I) and 99.95% (stage II–IV) (Fig. 7B). The results indicated that the DNAm levels of the five CpGs were potentially effective biomarkers for distinguishing tumor tissues from normal tissues in the early stages (stage I–II).

Constructing and validating the prognostic model

To establish an effectively prognostic classifier for ccRCC, a series of rigorous criteria were used to filtrate biomarkers (see “Materials and Methods”). Table S8 showed that 21 CpGs were significantly correlated with overall survival based on univariate Cox regression analysis (p-value = 0.035 – 3.8 × 10⁻⁶). Subsequently, the multivariate LASSO Cox regression analysis was performed to screen prognostic biomarkers, and 8 of 21 CpGs were obtained (Figs. 8A–8C). Based on the β-values of the eight CpGs, the weighting coefficients were used to calculate risk scores for patients as follow: risk-score = −1.77 × cg19516340 + 0.43 × cg22884714 + 0.72 × cg09203199 + 1.42 × cg16836311 + 0.14 × cg15357821 −0.18 × cg03415545 + 0.86 × cg15200711 + 0.39 × cg09799983. The risk scores was used to evaluate survival status, patients with higher risk scores had poor survival status comparing with those with lower risk scores in two datasets (Figs. 8D and 8F). The median risk score (−0.13) was used as the cut-off, and patients were significantly grouped into high-risk and low-risk groups in the training datasets. The results showed that the risk-score was significantly correlated with risk of death (HR: 5.84; 95% confidence interval (CI) [3.16–10.8]; log-rank test p-value < 0.001; Fig. 8E). To assess the validity and reproducibility of prognostic classifier, patients in the validation datasets were also divided into high-risk or low-risk groups (HR: 4.63; CI [1.39–15.5]; log-rank test p-value = 2.9 × 10⁻³; Fig. 8G). The genes associated with the eight CpGs were identified, including ADH1C, CES2, CYP1B1, LPCAT1, HOOK3, RRM2, CHEK2 and MAN1C1. The prognostic value of RRM2, CES2, HOOK3, ADH1C and LPCAT1 has been reported and validated in many cancers, such as breast cancer, pancreatic cancer, colorectal cancer, bladder cancer, ovarian cancer, non-small cell lung cancer, head and neck cancer (Liu et al., 2007; Morikawa et al., 2010; Rahman et al., 2013; Xia et al., 2017; Yoshida et al., 2011). MAN1C1 and CHEK2 have been implicated in tumor immune response, cancer cell proliferation and cell-to-cell adhesion (Zhang et al., 2010), and provided potential therapeutic targets for RCC (Li et al., 2018). However, the molecular biological function of CYP1B1 was still unknown. The results showed that those genes with abnormal methylation may paly crucial roles in tumorgenesis, and the DNAm sites might also be potential targets for new therapies.

The clinical prognostic signatures, including age, gender and TNM stage, were used as covariates to adjust prognostic classifier, and the results showed that the risk score calculated by the combination of eight CpGs was an independent prognostic factor with the better performance in the TCGA datasets (Table 3). Based on the subsets of ccRCC patients with different clinical signatures, the stratification analysis of the prognostic classifier was further implemented. When ccRCC patients were grouped according to the clinical signatures (age, gender, TNM stage), the prognostic classifier (Fig. S6) still had the better performance (log-rank test p-value < 0.05). The Fig. S7 shown that the ccRCC patients with the same clinical characteristics were also significantly divided into high-risk and low-risk groups by the prognostic classifier (log rank test p-value < 0.05). To further evaluate the universality of the construction methodology, we used five different cancers with the clinical information as the independent validation datasets (TCGA-BRCA, TCGA-KIRP, TCGA-LIHC, TCGA-LUAD and TCGA-UCEC). Compared to patients in the low-risk group, patients in the high-risk group had significantly shorter overall survival in all five different cancers (log-rank test, all p-value < 0.05; Fig. S8). In summary, these results clearly demonstrated that investigation of the DNAm mechanism with the multi-omics analysis may identify potential pathogenic factors and provide potential biomarkers for prognosis analysis of patients.