Research Paper Volume 13, Issue 5 pp 7361—7381

Artificial intelligence prediction model for overall survival of clear cell renal cell carcinoma based on a 21-gene molecular prognostic score system

Qiliang Peng1,2, *, , Yi Shen3, *, , Kai Fu4, *, , Zheng Dai4, , Lu Jin4, , Dongrong Yang4, , Jin Zhu4, ,

  • 1 Department of Radiotherapy and Oncology, The Second Affiliated Hospital of Soochow University, Suzhou, China
  • 2 Institute of Radiotherapy and Oncology, Soochow University, Suzhou, China
  • 3 Department of Radiation Oncology, The Affiliated Suzhou Science and Technology Town Hospital of Nanjing Medical University, Suzhou, China
  • 4 Department of Urology, The Second Affiliated Hospital of Soochow University, Suzhou, China
* Equal contribution

Received: November 2, 2020       Accepted: January 14, 2021       Published: March 3, 2021      

https://doi.org/10.18632/aging.202594
How to Cite

Copyright: © 2021 Peng et al. This is an open access article distributed under the terms of the Creative Commons Attribution License (CC BY 3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

Abstract

We developed and validated a new prognostic model for predicting the overall survival in clear cell renal cell carcinoma (ccRCC) patients. In this study, artificial intelligence (AI) algorithms including random forest and neural network were trained to build a molecular prognostic score (mPS) system. Afterwards, we investigated the potential mechanisms underlying mPS by assessing gene set enrichment analysis, mutations, copy number variations (CNVs) and immune cell infiltration. A total of 275 prognosis-related genes were identified, which were also differentially expressed between ccRCC patients and healthy controls. We then constructed a universal mPS system that depends on the expression status of only 21 of these genes by applying AI-based algorithms. Then, the mPS were validated by another independent cohort and demonstrated to be applicable to ccRCC subsets. Furthermore, a nomogram comprising the mPS score and several independent variables was established and proved to effectively predict ccRCC patient prognosis. Finally, significant differences were identified regarding the pathways, mutated genes, CNVs and tumor-infiltrating immune cells among the subgroups of ccRCC stratified by the mPS system. The AI-based mPS system can provide critical prognostic prediction for ccRCC patients and may be useful to inform treatment and surveillance decisions before initial intervention.

Introduction

Renal cell carcinoma (RCC) is one of the most lethal cancer types in the urinary system, which accounts for 85%-95% of renal malignancies and 2%-3% of all human tumors [1]. Clear cell renal cell carcinoma (ccRCC) is the most common subtype of RCC, accounting for approximately 70%-85% of the pathological types of cases [2]. In the past three decades, the incidence of ccRCC has increased gradually and approximately 30% patients already have metastasis when they were first diagnosed with ccRCC. Currently, treatments for localized ccRCC mainly depend on radio-frequency ablation and partial or radical nephrectomy; nevertheless, when ccRCC progresses to distant metastasis, the curative effect of present targeted drug therapies is unsatisfactory [3]. To this end, it is necessary to identify accurate predictors in ccRCC that are helpful for detecting it early, monitoring tumor progression, revealing survival outcome, promoting personalized therapy and for the purpose of individualized follow-up strategies.

Currently, the American Joint Committee on Cancer (AJCC) tumor node metastasis (TNM) staging system and Fuhrman grade have been the most widely accepted clinical classification systems for prognostic prediction of ccRCC. In addition, the University of California Integrated Staging System (UCISS) and stage, size, grade, and necrosis (SSIGN) score are also applied [4]. However, ccRCC patients with the same clinical stage can still have different survival outcomes, suggesting that these methods cannot accurately evaluate and manage patients with ccRCC [5]. Recent rapid advances in research techniques have promoted the development of various molecular prognostic indicators [6]. A series of genetic or epigenetic alterations have been demonstrated to play a vital role in the initiation and progression of ccRCC and may serve as potential biomarkers for clinical application [7]. Nevertheless, none of the tests developed to date are sufficiently qualified to predict overall survival. For this reason, a more effective and comprehensive tool for understanding of the molecular classification characteristics of ccRCC is urgently required.

Artificial intelligence (AI) has been known as a field of science and engineering involved in the computational understanding of intelligent behavior such as perceiving one’s environment, acting in complex environments, learning and understanding from experience, using reasoning to solve problems and to discover hidden knowledge, and applying knowledge successfully in new situations [8]. Machine learning, a branch of AI technology, has been recognized as a collection of data-analytical techniques aimed at building predictive models from multi-dimensional datasets [9]. With the utilization of AI-based machine learning method, new data could be predicted based on the identified patterns through detecting difficult-to-recognize patterns from complex combinations of multiple biomarkers [10]. Based on these advantages, AI-based machine learning has gained great attention and has been widely applied in medical fields for image detection, diagnosis, and outcome prediction [11].

In the present study, we created a system applying AI-based machine learning approaches and applied it to establish a universal molecular prognostic score (mPS) for overall survival prediction of ccRCC that relies on the expression status of 21 genes. During the construction of mPS, we applied two prediction models, a random forest prediction model and an artificial neural network model, for predicting the survival outcome of ccRCC. We then validated the mPS system by using different patient cohorts. Afterwards, a comprehensive nomogram comprising the mPS system and several independent variables were set up to predict ccRCC patient survival. In addition, we compared the differences among three subgroups of ccRCC samples, which were divided according to mPS scores, regarding mutations and copy number variations (CNVs), and immune cell infiltrations.

Results

Differentially expressed genes related to the prognosis of ccRCC

The study design was plotted in Figure 1. To study the specific genes associated with the prognosis in patients with ccRCC, we downloaded gene expression data from 531 ccRCC samples with clinical survival information from The Cancer Genome Atlas (TCGA) database. The detailed characteristics are presented in Table 1. We divided the 531 ccRCC samples into a high expression group and a low expression group according to the median expression of genes, and then observed the relationships between gene expression level and overall survival (OS) in these two groups. A total of 2317 OS-related genes in the TCGA discovery cohort were identified including Interleukin 20 Receptor Subunit Beta (IL20RB) and KL (Klotho) (Figure 2A, 2B). Then the gene differential expression analysis was conducted on 531 ccRCC samples and 72 healthy control samples in TCGA. Based on the Limma R package, a total of 2275 genes were identified to differentially express between ccRCC patients and healthy controls. After intersection with the prognosis-related genes and differential expressed genes, a total of 275 differentially expressed genes which were associated with the prognosis of ccRCC were retained (Figure 2C). In particular, IL20RB and KL were the most promising prognosis-related genes, with the lowest and highest hazard ratios (HRs), respectively, in the TCGA discovery cohort. The integrated HR for top 40 prognosis related genes in the TCGA data sets were plotted at Figure 2D.

Study pipeline. In this study we first used the TCGA clear cell renal cell carcinoma (ccRCC) cohort to identify a list of 275 genes, which were associated with the prognosis of ccRCC patients according to survival analysis and were also differentially expressed between ccRCC patients and healthy controls. Then, artificial intelligence (AI) methods including random forest and neural network were applied to establish the mPS system based on 21 prognosis-related genes. Then, we validated the mPS system by using the ICGC cohort. Next, we found that the mPS system could be applied to ccRCC subsets. Moreover, we evaluated the mPS system by conducting univariate and multivariate Cox regression analysis of the TCGA dataset and built a nomogram comprising the mPS score and several independent variables to predict ccRCC patient prognosis. Finally, we explored the potential mechanisms underlying the mPS system by performing gene set enrichment analysis (GSEA), mutations, copy number variations (CNVs) and immune cell infiltration analysis.

Figure 1. Study pipeline. In this study we first used the TCGA clear cell renal cell carcinoma (ccRCC) cohort to identify a list of 275 genes, which were associated with the prognosis of ccRCC patients according to survival analysis and were also differentially expressed between ccRCC patients and healthy controls. Then, artificial intelligence (AI) methods including random forest and neural network were applied to establish the mPS system based on 21 prognosis-related genes. Then, we validated the mPS system by using the ICGC cohort. Next, we found that the mPS system could be applied to ccRCC subsets. Moreover, we evaluated the mPS system by conducting univariate and multivariate Cox regression analysis of the TCGA dataset and built a nomogram comprising the mPS score and several independent variables to predict ccRCC patient prognosis. Finally, we explored the potential mechanisms underlying the mPS system by performing gene set enrichment analysis (GSEA), mutations, copy number variations (CNVs) and immune cell infiltration analysis.

Table 1. TCGA cohort ccRCC patient characteristics.

Clinical characteristicVariableTotalPercentages (%)
Age (years)>6026650.09
≤6026549.91
GenderFemale18635.03
Male34564.97
LateralityLeft25047.08
Right28052.73
Others10.19
Histological gradeG1142.64
G222742.75
G320738.98
G47514.12
GX50.94
Unknown30.56
TT127251.22
T26912.99
T317933.71
T4112.07
NN023945.01
N1163.01
NX27651.98
MM042079.10
M17814.69
MX315.84
Unknown20.38
Pathological stageStage I26650.09
Stage II5710.73
Stage III12323.16
Stage IV8215.44
Unknown30.56
Survival statusAlive35967.61
Dead17232.39
Identification of prognosis-related genes in the TCGA ccRCC cohort. (A, B) Kaplan-Meier curves of OS for the TCGA cohort based on IL20RB and KL expression levels, respectively; (C) Venn diagram of prognostic-related genes and differentially expressed genes in TCGA cohort; (D) The integrated HR for Top 40 prognosis-related genes in the TCGA data sets.

Figure 2. Identification of prognosis-related genes in the TCGA ccRCC cohort. (A, B) Kaplan-Meier curves of OS for the TCGA cohort based on IL20RB and KL expression levels, respectively; (C) Venn diagram of prognostic-related genes and differentially expressed genes in TCGA cohort; (D) The integrated HR for Top 40 prognosis-related genes in the TCGA data sets.

AI-based development of the mPS system

In order to verify whether these 275 genes are enough to predict the survival rate of ccRCC patients at 3 years, we used the TCGA ccRCC cohort to construct a molecular prognosis score system as illustrated in detail in the Materials and methods. We first used the random forest algorithm for feature gene screening. When the cutoff of feature importance was set to 0.5, a total of 21 genes were identified and then used for the following artificial neural network construction. The survival curves of the 21 genes were illustrated at Figure 3.

The survival curves of the 21 prognosis-related genes used for constructing the mPS system.

Figure 3. The survival curves of the 21 prognosis-related genes used for constructing the mPS system.

The neural network model consisted of an input layer, two hidden layers and an output layer. The 21 genes were used as the input layer, and the middle 2 hidden layers included 4 neurons and 2 neurons, respectively. The output layer used the softmax activation function and consisted of two dimensions, each of which represented an outcome: alive or deceased. Based on the constructed neural network model, we obtained the weight of each gene, and then built a mPS system, which was measured by summation of “Gene-Score” × “Gene-Weight” for all 21 genes. For the TCGA cohort, the mPS value ranges from 0 to 71.57 (Table 2). Under each mPS value, the samples were stratified into three groups: low-mPS (mPS < 22), median-mPS (22 ≤ mPS < 30), and high-mPS group (mPS ≥ 30). Then, survival analysis was conducted with the Kaplan-Meier (K-M) method and log-rank statistical test. The results indicated that there was a significant difference among the three groups (Figure 4A).

Table 2. The 21 genes necessary and sufficient for calculation of mPS.

SymbolGene IDFull nameScore
(high)
Score
(low)
Weight
DDAH123576Dimethylarginine dimethylaminohydrolase 1103.580
CRABP21382Cellular retinoic acid binding protein 2013.009
TGFA7039Transforming growth factor alpha103.425
SEMA3G56920Semaphorin 3G103.042
SPATA18132671Spermatogenesis associated 18104.169
PTTG19232Pituitary tumor-transforming gene 1013.344
SCGN10590Secretagogin104.200
CYP39A151302Cytochrome P450 family 39 subfamily A member 1103.595
CLDN41364Claudin 4103.541
ZNF39555893Zinc finger protein 395103.977
IL15RA3601Interleukin 15 receptor subunit alpha013.859
APLNR187Apelin receptor103.072
APOLD181575Apolipoprotein L domain containing 1103.207
NTN459277Netrin 4103.723
PABPC1L80336Poly(A) binding protein cytoplasmic 1 like013.670
UBE2C11065Ubiquitin conjugating enzyme E2 C014.217
CDH41002Cadherin 4103.183
GNG72788G protein subunit gamma 7103.515
CEACAM1634CEA cell adhesion molecule 1103.506
PLAUR5329Plasminogen activator, urokinase receptor013.418
SIM26493SIM bHLH transcription factor 2013.364
mPS system can precisely stratify prognosis of ccRCC patients and is applicable to T1 and T2 ccRCC subsets. (A) Kaplan-Meier curves of OS according to mPS for the TCGA cohort; (B) Kaplan-Meier curves of OS according to mPS for the ICGC cohort; (C) Kaplan-Meier curves according to mPS for OS of patients at clinical T1 and T2 stage in the TCGA cohort; (D) Kaplan-Meier curves according to mPS for OS of patients at clinical T1 and T2 stage in the ICGC cohort.

Figure 4. mPS system can precisely stratify prognosis of ccRCC patients and is applicable to T1 and T2 ccRCC subsets. (A) Kaplan-Meier curves of OS according to mPS for the TCGA cohort; (B) Kaplan-Meier curves of OS according to mPS for the ICGC cohort; (C) Kaplan-Meier curves according to mPS for OS of patients at clinical T1 and T2 stage in the TCGA cohort; (D) Kaplan-Meier curves according to mPS for OS of patients at clinical T1 and T2 stage in the ICGC cohort.

Validation of the mPS system in independent cohorts

To assess whether mPS is feasible for the prediction of prognosis in other independent ccRCC cohorts, we investigated another ccRCC data set from International Cancer Genome Consortium (ICGC, https://icgc.org/). The detailed characteristics of ICGC cohort are presented in Table 3. We first calculated the mPS score of each sample in the ICGC data, and then divided the samples into different groups according to the above mentioned methods. Since the sample size of ICGC database is too small, we only divided the ICGC dataset into two groups: low-mPS (mPS < 22) and high-mPS (mPS ≥ 22) groups. The results revealed that there was a significant difference about the prognosis among different mPS groups. These findings demonstrated that mPS can stratify prognosis in different ccRCC cohorts (Figure 4B).

Table 3. ICGC cohort ccRCC patient characteristics.

Clinical characteristicVariableTotalPercentages (%)
Age (years)>604549.45
≤604650.55
GenderFemale3942.86
Male5257.14
StatusAlive6167.03
Deceased3032.97
TT15156.04
T299.89
T32729.67
T411.10
Unknown33.30
NN07986.81
N122.20
NX1010.99
MM08189.01
M199.89
MX11.10
Fuhrman gradeG11314.29
G24852.75
G31516.48
G41415.38
Unknown11.10

Application of the mPS system to ccRCC subsets

We further investigated whether mPS is also applicable to distinguish samples with different prognosis in different ccRCC subsets in TCGA data set. The results indicated that there was a significant difference regarding the prognosis among different mPS groups of patients in low-, median- and high-mPS groups at the clinical T1 and T2 stage in TCGA cohort (Figure 4C). Meanwhile, we also examined the utility of mPS for different TNM tumor stages determined from clinical information in ICGC cohort. The results suggested that the mPS system also stratified OS of patients in low- and high-mPS groups at the clinical T1 and T2 stage in the ICGC cohort (Figure 4D). We thus conclude that the mPS system can further stratify different ccRCC subsets in different ccRCC cohorts.

Univariate and multivariate Cox regression of mPS

To determine whether the prognostic ability of mPS for survival prediction was independent from other clinical parameters as a prognostic factor for ccRCC, we carried out univariate Cox regression analysis and multivariate Cox regression analysis on mPS and other clinical factors (age, gender, grade, stage, metastasis, lymph node, hemoglobin, platelet, and calcium result). The univariate analysis revealed that the variables of age, histological grade, stage, metastasis result, lymph node result, platelet result, hemoglobin and mPS score were significantly related to the prognosis of ccRCC patients (P<0.05). When entered into the significant factors into multivariate Cox regression analysis, the mPS score (HR=0.46, P<0.001), metastasis result (HR=1.83, P=0.002), platelet result (HR=1.74, P=0.016), hemoglobin result (HR=1.60, P=0.007), age (HR=1.58, P=0.003) and stage (HR=2.18, P=0.025) were found to be independently associated with OS in the entire TCGA cohort (Table 4). These results demonstrated that the mPS system could be applied independently to predict the survival outcomes of ccRCC patients.

Table 4. Univariate and multivariate analyses of OS in the TCGA cohort.

VariablesUnivariate analysisMultivariate analysis
HR95%CIPHR95%CIP
mPSHigh vs. Low0.310.22-0.42<0.0010.460.32-0.65<0.001
Grade3&4 vs. 1&22.601.85-3.67<0.0011.410.98-2.020.07
GenderMale vs. Female0.940.69-1.280.700
HemoglobinLow vs. Normal2.371.66-3.39<0.0011.601.14-2.260.007
AgeOlder vs. Younger1.701.26-2.29<0.0011.581.17-2.140.003
MetastasisM1 vs. M04.333.17-5.92<0.0011.831.24-2.700.002
Lymph nodeN1 vs. N03.431.82-6.460.0011.300.67-2.500.44
StageIII&IV vs. I&II3.932.86-5.39<0.0012.181.10-4.320.025
T stageT3&T4 vs. T1&T23.212.39-4.35<0.0010.890.49-1.620.72
CalciumLow vs. Normal0.820.57-1.160.300
PlateletElevated vs. Normal3.762.49-5.68<0.0011.741.11-2.710.016

Establishment of the nomogram model for survival prediction

In order to gain an effective tool to predict the survival probability with ccRCC patients, a nomogram was constructed based on mPS and clinical information using the TCGA cohort. After integrating all significant independent variables, the final model was established and presented as a nomogram to estimate the probability of the 3- and 5-year OS for ccRCC (Figure 5A). C-index was calculated to assess the predictive value of the model. The nomogram had a C-index of 0.79, which revealed relatively high discrimination ability. In addition, the calibration plots exhibited optimal agreement between the model prediction and actual observation for predicting ccRCC survival probability at 3- and 5- year (Figure 5B, 5C).

Nomogram construction results. (A) Nomogram to predict the 3- and 5-year OS for ccRCC patients in the TCGA cohort; (B, C) calibration curves for the nomogram model of the 3- and 5-year OS.

Figure 5. Nomogram construction results. (A) Nomogram to predict the 3- and 5-year OS for ccRCC patients in the TCGA cohort; (B, C) calibration curves for the nomogram model of the 3- and 5-year OS.

Gene set enrichment analysis-based KEGG analysis

Based on the Limma R package, a total of 434 genes were identified to differentially express between low- and high-mPS groups, of which 279 genes were up-regulated and 155 genes were down-regulated in high-mPS group compared with low-mPS group. The principal component analysis result was plotted at Figure 6A while the volcano plot of the differentially expressed genes was illustrated at Figure 6B. In addition, the heat map of the top 50 differentially expressed genes was presented at Figure 6C. The differentially expressed genes were then selected for Gene set enrichment analysis (GSEA) analysis. As a result, a total of 20 prominent Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways including activated and suppressed pathways were selected (Figure 6D). Activated pathways mainly includes ultraviolet (UV) response, transforming growth factor-beta (TGF-beta), heme metabolism, adipogenesis, and bile acid metabolism while suppressed pathways were mainly concentrated on IL6/JAK/STAT3 signaling, DNA repair, allograft rejection, complement, E2F targets, epithelial mesenchymal transition, G2M checkpoint, and inflammatory response. GSEA enrichment plots of representative gene sets on “epithelial mesenchymal transition”, “IL6/JAK/STAT3 signaling”, “complement”, “E2F targets”, and “DNA repair” were also shown in Figure 6E6I.

GSEA analysis results of differentially expressed genes between low- and high-mPS groups. (A) Principal component analysis result; (B) Volcano plot of the differentially expressed genes; (C) The heat map of the top 50 differentially expressed genes; (D) Significantly enriched activated and suppressed KEGG pathways. The vertical items are the names of KEGG terms, and the length of horizontal graph represents the gene ratio. The depth of the color represents the adjusted p-value. The area of circle in the graph means gene counts. (E–I) GSEA-based KEGG-enrichment plots of representative gene sets from activated and suppressed pathways: “IL6/JAK/STAT3 signaling” (E), “DNA repair” (F), “complement” (G), “E2F targets” (H), and “epithelial mesenchymal transition” (I).

Figure 6. GSEA analysis results of differentially expressed genes between low- and high-mPS groups. (A) Principal component analysis result; (B) Volcano plot of the differentially expressed genes; (C) The heat map of the top 50 differentially expressed genes; (D) Significantly enriched activated and suppressed KEGG pathways. The vertical items are the names of KEGG terms, and the length of horizontal graph represents the gene ratio. The depth of the color represents the adjusted p-value. The area of circle in the graph means gene counts. (EI) GSEA-based KEGG-enrichment plots of representative gene sets from activated and suppressed pathways: “IL6/JAK/STAT3 signaling” (E), “DNA repair” (F), “complement” (G), “E2F targets” (H), and “epithelial mesenchymal transition” (I).

Landscape of mutation profiles

The somatic mutation profiles of 531 ccRCC patients were also downloaded from TCGA. According to the above analysis, the TCGA ccRCC samples can be divided into three groups, namely low-mPS group (mPS<22), median-mPS group (22≤mPS<30), and high-mPS group (mPS≥30), and the number of samples in these three groups were 129, 55, 347, respectively. The “maftools” package was applied to visualize the results according to the mutation data with VCF format. The statistically significant mutated genes were identified. Here we mainly focused on the top ten significantly mutated genes for in-depth analyses. A total of 16 significantly mutated genes were then retrieved after taking the combination of the top 10 mutated genes in the three groups. The mutation distributions and annotations of the identified 16 genes in these three groups of samples were illustrated in waterfall plot (Figure 7A7C). The results indicated that there were significant differences regarding the mutations among the three groups stratified by mPS system, which may help us understand the promising utility of mPS to some degree.

Mutations and copy number variations in the three subgroups stratified by mPS system. (A–C) Distribution and phenotype of common gene mutations in the three subgroups stratified by mPS system: (A) low-mPS group (mPSB) median-mPS group (22≤mPSC) high-mPS group (mPS>30). (D) The levels of amplification and deletion of chromosome arms of the three subgroups stratified by mPS system. Objects with pentagrams indicate there is a significant difference in the frequency distribution among the three subgroup (P E–G) The copy number variations located on minimal common regions (MCRs) of the three subgroups stratified by mPS system: (E) low-mPS group (mPSF) median-mPS group (22≤mPSG) high-mPS group (mPS≥30).

Figure 7. Mutations and copy number variations in the three subgroups stratified by mPS system. (AC) Distribution and phenotype of common gene mutations in the three subgroups stratified by mPS system: (A) low-mPS group (mPS<22); (B) median-mPS group (22≤mPS<30); (C) high-mPS group (mPS>30). (D) The levels of amplification and deletion of chromosome arms of the three subgroups stratified by mPS system. Objects with pentagrams indicate there is a significant difference in the frequency distribution among the three subgroup (P < 0.01); (EG) The copy number variations located on minimal common regions (MCRs) of the three subgroups stratified by mPS system: (E) low-mPS group (mPS<22); (F) median-mPS group (22≤mPS<30); (G) high-mPS group (mPS≥30).

Copy number variation analysis

The segmented CNVs in the three subgroups of ccRCC samples were applied to screen the significantly frequent CNVs by the Genomic Identification of Significant Targets in Cancer (GISTIC) method. Specifically, in the low-mPS subgroups (mPS<20), the amplifications of 7q, 7p and deletions of 9p, 9q were identified more significant. In addition, in the median-mPS subgroups (22≤mPS<30), the amplifications of 5q, 5p and deletions of 3p, 4q were considered statistically significant. Moreover, in the high-mPS subgroups (mPS≥30), frequent amplifications were observed in chromosomal arms 7q, and 7p, and deletions were observed in 9p, and 9q. Chi-square test was performed on the frequency distribution of each chromosome arm. The frequency distribution of amplification and deletion of chromosome arms in the three groups of samples obtained by GISTIC analysis were presented at Figure 7D.

Then, the CNVs located on minimal common regions (MCRs) were detected in the three subgroups of ccRCC samples (Figure 7E7G). In the low-mPS subgroups, 26 MCRs of CNVs were identified including 11 amplifications and 15 deletions. The most significant gained regions were located at 5q35.3, while the most significant regions of loss were located at 3p22.1 and 3p25.3. The distribution of these MCRs in the genome was shown in Figure 7. In the median-mPS subgroups, a total of 8 amplifications and 8 deletions of MCRs were detected. In detail, regions 5q23.3 and 5q23.2 with amplifications and 3p22.3, 3p12.3 with deletion were the most frequently altered sites. Moreover, a total of 6 amplifications and 6 deletions of MCRs were detected in the high-mPS subgroups. The most significant amplification regions were 5q35.3, 5q21.3, etc., while the most significant deletion regions were 3p21.31, 3p14.1, etc. Collectively, these findings suggested that the level of chromosomal amplification and deletion was correlated with the mPS score and prognosis, as a worse outcome subgroup (low-mPS subgroup) has more genomic abnormality compared with other groups. In general, a worse outcome may be correlated with more genomic abnormality, revealing that the chromosomal amplification and deletion identified by mPS system can be a patient stratification for immune checkpoint therapy.

Immune cell infiltration in the tumor microenvironment

Immune cell infiltration has been recognized as one of the most important factors that may positively or negatively shape tumor initiation, progression, or survival outcome. In this study, we analyzed the composition of tumor-infiltrating immune cells in ccRCC samples based on the RNA-seq data of ccRCC samples in TCGA using CIBERSORT. When characterizing the abundances of different immune cell types with CIBERSORT, 22 subsets of tumor-infiltrating immune cells in 531 samples were obtained at a threshold of P-value < 0.05, including memory B cells, activated dendritic cells, and M0 macrophage. The distribution results of the inferred fractions of immune cell populations in different samples produced by CIBERSORT were illustrated at Figure 8A. In addition, Figure 8B depicted the differences of the ratio distributions of 22 infiltrating immune cells in the three subgroups of samples including low-mPS, median-mPS, and high-mPS group. The significant immune gene signatures among the three groups included signatures for naive B cells, plasma cells, CD8 T cells, CD4 memory resting T cells, CD4 activated memory T cells, follicular helper T cells, regulatory tregs T cells, resting NK cells, monocytes, M0 macrophages, M1 macrophages, M2 macrophages, resting mast cells, and neutrophils, suggesting that considerable variability existed in the nature of the tumor immune infiltrate across different subgroups of ccRCC stratified by the mPS system, partly accounting for molecular features of the tumor underlying mPS system.

Immune subtypes in patients with ccRCC. (A) Unsupervised clustering of all samples based on immune cell proportions in low-, median- and high-mPS groups; (B) 22 types of adaptive and innate immune cells in low-, median- and high-mPS groups (ns, P>0.05; *,P

Figure 8. Immune subtypes in patients with ccRCC. (A) Unsupervised clustering of all samples based on immune cell proportions in low-, median- and high-mPS groups; (B) 22 types of adaptive and innate immune cells in low-, median- and high-mPS groups (ns, P>0.05; *,P<0.05; **, P<0.01; ***, P<0.001; ***, P<0.0001;by unpaired two-tailed t test).

Discussion

Although the advances of diagnostic techniques and comprehensive treatment have enhanced the survival and life quality of ccRCC patients to some extent, about 30% of patients undergoing curative nephrectomy still developed recurrence and metastasis when the clinical outcomes are poor. In recent years, the utilization of high-throughput sequencing technology has brought up a series of genome-wide biomarkers for ccRCC. However, they are still not enough for accurate prediction and individualized treatment. Therefore, exploring prognostic tools to more accurately manage patients with ccRCC with poor survival is becoming increasingly significant. In this study, we developed a new molecular prognostic signature (mPS system) to predict the overall survival of ccRCC patients.

We first identified 275 differentially expressed genes which were associated with the prognosis of ccRCC. We then used two computational learning models (neural network and random forest) to further screen 21 genes for establishing the mPS system, which combined gene expression level, HR value and gene weight. The survival analysis indicated that there was a significant difference among the three subgroups stratified by the mPS system. Moreover, the prognostic value of the mPS system was validated in independent ccRCC cohorts. Besides, we confirmed that the mPS system is also applicable to distinguish samples with different prognosis in different ccRCC subsets. In our study, we demonstrated that the mPS system may offer valuable information to clinicians regarding variables that are the most useful for patient stratification.

As we all know, ccRCC is a heterogeneous tumor with various confounding factors involved in the initiation and progression of this disease. In this study, univariate and multivariate Cox regression analyses were conducted to assess the correlation between several clinical factors, along with mPS and survival outcome of ccRCC patients. The results showed that mPS score and several clinical factors were independent prognostic factors in patients with ccRCC. Afterwards, a promising clinicopathologic prognostic nomogram was constructed by modeling mPS score and significant factors to predict the 3-, and 5-year OS for patients with ccRCC. Nomograms have frequently served as a vital component of modern medical decision-making under complicated clinical conditions without needing standard guidelines [12]. They can make the results of prediction models clear by presenting visual graphical interfaces, thus providing an opportunity to assess various patient-based variables and predicting a patient’s individualized risk for survival or a specific outcome [13]. Our developed nomogram has indicated moderately high discrimination and sufficient calibration. We suggest using this nomogram to further improve the personalized management of ccRCC patients.

After reviewing the existing literatures, we found that these twenty-one genes are more or less related to tumors. Among them, ten genes were identified directly associated with ccRCC based on PubMed literature reports. For example, Semaphorin 3G (SEMA3G) has been identified as an immune-related gene which is associated with the prognosis of ccRCC patients [14, 15]. Pituitary tumor-transforming gene 1 (PTTG1) is a recently identified oncogene involved in the progression of malignant tumors, and the expression of the PTTG1 oncogene has been confirmed to be related to progression and prognosis of ccRCC patients [16]. There is growing evidence that Claudin 4 (CLDN4) plays a significant role in the occurrence and development of chromophobe renal cell carcinoma and could serve as an independent risk biomarker for the overall survival in patients with chromophobe renal cell carcinoma [17]. Zinc finger protein 395 (ZNF395), a ccRCC master regulator, was identified as a downstream protein of CCDC50-S, and the interaction initiated oncogenic pathways which were highly associated with the pathogenesis of ccRCC [18]. Apelin receptor (APLNR) expression was negatively related to PD-L1 expression by tumor cells in a subset of patients with ccRCC and the expression of APLNR has been recognized as an independent prognostic factor for survival of patients with ccRCC [19]. Recent evidence has identified Ubiquitin conjugating enzyme E2 C (UBE2C) as a pathological stage-relevant gene, which is associated with carcinogenesis and progression in ccRCC, and may provide potential diagnostic, therapeutic and prognostic biomarkers for ccRCC [20]. Previous studies have indicated Cadherin 4 (CDH4) was significantly linked with the survival outcome of ccRCC patients and the transcriptional level of CDH4 may serve as an effective diagnostic and prognostic biomarker for ccRCC patients [21]. Recent findings reveal that G protein subunit gamma 7 (GNG7) is a tumor suppressor in ccRCC progression and has a potential to serve as a new biomarker or therapeutic target in ccRCC [22]. CEACAM1 has been known as a tumor suppressor in various epithelial tumors. Transient expression of CEACAM1 by tumor cells and subsequent homophilic interaction with CEA cell adhesion molecule 1 (CEACAM1) on tumor-infiltrating lymphocytes may provide a novel immune escape mechanism in ccRCC. Moreover, Plasminogen activator, urokinase receptor (PLAUR) was developed and validated as a prognostic immune-associated gene signature in ccRCC along with other six genes [23]. However, eleven genes (CRABP2, DDAH1, TGFA, SPATA18, SCGN, CYP39A1, IL15RA, APOLD1, NTN4, PABPC1L, and SIM2) have not been studied in relation to ccRCC, according to PubMed searches for “GENE and ccRCC”. Further basic and clinical studies are required to validate our observations and the mechanisms underlying the prognostic value of 21 genes on which mPS is based and for the development of novel therapeutic targets to prolong OS of ccRCC patients.

Recently, AI-based algorithms have been successfully applied in medical fields as they are capable of learning patterns from massive, complex datasets and creating useful predictive outputs [24]. Moreover, AI-based technologies are well adapted to rapid, high-volume data processing, promoting tailored and specific management according to each patient’s characteristics [25]. In particular, artificial neural networks are computational models inspired by biological neural networks, which are currently most common practiced models of AI-based algorithms developed for survival prediction and decision-making [26]. Another AI technique is the random forest tool, which is a multivariate prediction model that conducts a computationally extensive and robust data-mining and could handle large sets of proposed variables as inputs for screening variables related to the outcomes of interest [27]. Recently, Nakayama et al. established a universal mPS system that relies on the expression status of 23 genes by using AI models (random forest and neural network), which is applicable to almost all subsets of breast cancer patients [28]. Moreover, Aramini et al. set up a mPS system for stratifying early-stage, resected lung cancer patients, which could be used for predicting risk of distant recurrence and informing treatment and surveillance decisions [29]. In our study, we developed the mPS system based on two computational learning models (AI-based random forest and neural network), which has exhibited promising predictive performance for ccRCC patients.

Based on the mPS system, we preliminarily divided ccRCC patients into low- and high-mPS groups. Then, differentially expressed genes between these two groups were identified and used for GSEA analysis. A series of activated and suppressed pathways were identified including “epithelial mesenchymal transition”, “IL6/JAK/STAT3 signaling”, “complement”, “E2F targets”, and “DNA repair”, which are highly involved in the initiation and progression of ccRCC. For example, the epithelial mesenchymal transition (EMT) has now been recognized as a complicated biological trans-differentiation process, promoting epithelial cells to transiently gain mesenchymal features, containing motility and metastatic potential [30]. EMT has been highly involved in a extensive range of malignant tumor types including ccRCC and can be driven by a conserved set of inducing signals, transcriptional regulators and downstream effectors [31]. IL-6/JAK/STAT pathway is aberrantly hyper-activated in many types of cancer including ccRCC and plays a significant part in the progression of cancer cachexia through regulating the inflammatory response [32]. Studies have convinced the roles of DNA repair deregulation in promoting immune recognition and immune destruction of cancer cells [33]. Emerging evidence has identified that defects in DNA repair machinery play an important part in the pathogenesis and progression of human cancer [34]. Moreover, the information gathered so far indicated that E2F is one of the most important transcription factors that regulates various cellular functions related to cell cycle and apoptosis and involves not only in proliferation and tumorigenesis but also in apoptosis and differentiation [35]. These pathways may provide new insights for the occurrence and development of ccRCC.

In recent years, it is widely acknowledged that genetic or epigenetic alterations play an important role in a variety of transcriptional and non-transcriptional biological processes [36]. Mutations and CNVs have been vital forms of epigenetic modifications, which can cause genomic and molecular phenotype heterogeneity, promoting the initiation and progression of complex diseases, including cancer [37, 38]. In this study, we presented an integrative analysis of multiple types of genomic data, including mutations, and CNVs. We identified 16 significantly mutated genes, the mutation distributions and annotations of which were statistically different in the three groups stratified by mPS system. Moreover, the significantly frequent CNVs in the mPS subgroups of ccRCC samples were identified by GISTIC. In addition, the CNVs located on MCRs were also detected in the three subgroups of ccRCC samples. The results revealed that the CNVs were significantly different among different mPS subgroups. Moreover, a worse outcome is correlated with more genomic abnormality. These findings further supported the promising predictive value of mPS system and provided potential therapeutic targets for ccRCC patients.

Increasing evidence supports the importance of the immune infiltration of tumors on the initiation and progression of various tumors, which could provide potential biomarkers to improve the reliability and precision of diagnosis and prognosis [39]. Recent studies have extensively investigated the dysfunction of immune cell infiltration in promoting the metastasis and invasion of ccRCC through mediating tumor cell budding and disruption of focal basal cell layer [40]. In this study, we applied the newly developed algorithm CIBERSORT to estimate different expressional cell patterns of immune infiltration in different mPS subgroups. The distributions of immune cell fractions in ccRCC tissues were significantly different among different subgroups stratified by the mPS system. Several important immune cells were identified differentially expressed among the mPS subgroups, including naive B cells, CD8 T cells, CD4 memory resting T cells, CD4 activated memory T cells, which may explain the prognosis differences of the low-mPS, median-mPS, and high-mPS group.

The advantage of this study is that, firstly, we developed a simple and cost-effective mPS system by using multivariate models including neural network and random forest, which were beneficial to obtain a stable scoring system. Secondly, the scoring model was verified in the ICGC data set, and this scoring model is also effective for subgroups. Moreover, the sub-population samples were further characterized, adding mutation analysis, CNV analysis, and immune infiltrating cell ratio analysis. The mutations, CNVs and tumor immune infiltrations will provide new insights into carcinogenesis and might promote to identify potential therapeutic targets for ccRCC. We also combined the mPS system with several independent variables into a predictive nomogram, which could provide convincing information for helping clinicians to illustrate the prognosis of patients with ccRCC.

Nevertheless, there are inevitably several limitations of our study that should be acknowledged. To begin with, the number of samples with follow-up data of the validation cohort was limited due to availability, which limits the ability to adjust for confounding factors. Thus, additional large validation datasets are required to further improve prediction performance of the mPS system. Then, we were unable to compare the prognostic value of mPS with that of other representative models. Furthermore, these results from our study were provided through analytic methods and machine learning without experimental validation. Thus, it is of great importance to find the molecular mechanism for further characterization of mPS.

Conclusions

In summary, we developed a new molecular prognostic model with AI-based methods, called the mPS system, based on expression levels of twenty-one prognostic genes. The mPS system obtained from AI-based method can be used as a medical decision support system that may provide critical information for prognostic assessment of patients with ccRCC before initial intervention. In addition, we explored the potential mechanisms underlying the mPS system by performing GSEA analysis and evaluating the differences of mutations, CNVs and immune cell infiltrations among the subgroups stratified by the mPS system, which should be validated to further characterize the molecular background of the mPS system.

Materials and Methods

Data source

The mRNA expression data of TCGA ccRCC discovery cohort (comprising 72 normal specimens and 539 ccRCC specimens) were retrieved from the UCSC Xena database (https://xenabrowser.net/datapages/) [41]. Clinical data pertaining to patients’ age, gender, grade, stage, survival outcome were also acquired from the TCGA data portal by using cgdsr package in R, which contained clinical information of 531 ccRCC patient samples. We downloaded 523 SNP6 copy number samples with corresponding transcriptome data from http://firebrowse.org/. Moreover, the R-based TCGAbiolinks was employed to collect the SNP data of ccRCC, which consists of 380 samples with corresponding transcriptome data. We adopted another ccRCC data set from the ICGC database as the validation cohort, which covers 91 ccRCC samples with mRNA expression and clinical data [42].

Identification of prognosis-associated genes and differential expression analysis

The public data were retrieved from https://xenabrowser.net/datapages/. All human protein-coding genes were assessed the potential utility as a prognostic predictor with the TCGA ccRCC discovery cohort. The K-M method was performed for survival analysis by using R package survival and survminer [43]. According to the median value of gene expression, the ccRCC samples were divided into high or low expression group for a particular gene. The cox regression analysis in the survival package was applied to calculate the HR and 95% confidence interval (CI) of each prognosis-related gene. LIMMA R software package was used to analyze differential expressed genes [44]. The adj.P-value <0.05 and a |log2FC| >2 were regarded as the cutoff values of statistical significance.

Random forest prediction model

Random forest is one of the most widely used nonparametric techniques for data classification and regression analysis [45]. The idea of random forest is to establish multiple trees that use randomly generated samples from existing situations with a bootstrapping technique and to generate an average prediction of the individual trees [46]. The basic unit of random forest is a decision tree. To classify a sample, each tree in the forest is given a random input vector and all vectors in the random forest are independent and identically distributed. Random forest is designed to randomize the column variables and row observations of the dataset, generate multiple classification subgroups, and finally summarize the classification results. In this study, we applied the R package randomForest to build the random forest prediction model by using the ccRCC samples in TCGA. The expression status of 275 selected prognostic-related differentially expressed genes was used as the model input, while the survival status of the sample in the third year (alive=0, deceased=1) was selected as the output, where ntree was set to 500 and maximum depth (max-depth) was set to 10. To assess the generalized performance, we applied a stratified 10-fold cross-validation test (CV = 10). We measured the importance of each feature by referring to the average decrease value of node impurity: the larger the value, the greater importance of the variable. After a 10-fold cross-validation, we identified 21 genes on the basis of feature importance values for further artificial neural network mode (cutoff=0.5).

Artificial neural network mode

We first transformed the expression status of the 21 genes to “Gene-Score” based on the expression level (above or below the median) and integrated HR for each gene with the following step function:

Gene-score matrixGene expression
lowhigh
Integrated HR<110
>101

Subsequently, an artificial neural network system was constructed and trained. At each hidden neuron node, the ReLU (rectified linear unit) was applied as the activation function [47]. In the output layer, two nodes were set including a1 and a2, which represented alive and deceased, respectively. We employed a softmax function for each node, and designated y2 (probability of death; that is, the a2 node) as Y. We applied cross entropy error for the model loss function (E). Moreover, we used the Adam method for the optimization of each weight (learning rate=0.001; epochs=1000). After the neural network model was constructed, the weight of each node in the input layer (Gene-Weight) was used to calculate mPS (the sum of Gene-Score*Gene-Weight for all 21 genes). In our study, the Python-based Keras library was applied for neural network training.

Univariate and multivariate Cox regression analysis

Univariable and multivariable Cox regression analyses were applied to evaluate the associations between overall survival and potential confounders including mPS and other clinical variables [48]. Univariable Cox regression analysis was performed to explore the potential confounders. Afterwards, significant prognostic variables with a P-value < 0.05 in univariate analysis were further entered into a multivariate Cox proportional hazards model along with the corresponding 95% CI for each potential risk factor. P-values less than 0.05 were considered statistically significant.

Construction of prognostic nomogram

A nomogram was established to predict 3- and 5-year OS by including all independent prognostic factors of the multivariate Cox regression analysis by using RMS package in the R Statistical [49]. Afterwards, we illustrated the predicted and observed results in the calibration curves for the visualization of the predictive performance of the nomogram. The discrimination of the prognostic models was estimated and compared by Harrell’s concordance index (C-index). A higher C-index indicated a superior discriminative capacity for prognostic prediction. Generally, a C-index value ≥0.70 suggests a good fit [50].

Gene set enrichment analysis

GSEA was applied to further understand mPS-associated pathways [51]. The ccRCC patients in the TCGA cohort were divided into low- and high-mPS groups based on the median expression value of mPS. Then, the genes differentially expressed between low- and high-mPS groups were identified based on Student’s t-test in the Limma R package. The P-value < 0.05 was regarded as the cutoff value of statistical significance. Afterwards, GSEA was performed based on GSEABase package on the R platform.

Mutational pattern analysis

The mutation data of ccRCC sample in TCGA were processed and visualized by R package maftools. Maftool is an efficient and comprehensive tool, which provides various analysis and visualization modules that are frequently applied in cancer genomic studies, such as driver gene identification, pathway, signature, enrichment, and association analyses [52]. The top 10 frequently mutated genes of each group were identified and compared by performing Chi-square test. Tumor mutation burden (TMB) was defined as the total amount of coding errors of somatic genes, base substitutions, insertions or deletions detected across per million bases [53]. In the present study, the mutation frequency was calculated with number of variants/the length of exons (45 million) for each sample.

Copy number variations analysis

Loci of amplification and deletion were investigated across samples to define MCRs targeted by overlapping events in two or more. In our study, we used the GISTIC method to detect the common CNV regions in all samples based on SNP6 Copy Number segment data, including CNVs at the chromosomal arm level and MCRs between samples [54]. GISTIC is an algorithm designed to distinguish regions of variation that are more possible to trigger cancer pathogenesis. We chose the threshold of q-value<0.1 for the MutSigCV method. In addition, we selected the confidence level of 0.95 when the peak interval was determined. Moreover, the area (greater than 0.98) was applied as the standard when analyzing the variation of the chromosomal arm level. We performed these analyses by using the MutSigCV module in the online GenePattern analysis tool (https://cloud.genepattern.org/gp/pages/index.jsf), which is developed by Broad Research Institute [55].

Evaluation of immune cell infiltration

To quantify the proportions of immune cells among different subtypes of ccRCC, the CIBERSORT algorithm was applied, which is a deconvolution method that uses a series of reference gene expression values for a minimal representation of each cell type and investigates the cell component of complex tissues according to the gene expression data from massive tumor samples with support vector regression [56]. The normalized 531 ccRCC TCGA RNA-seq data were uploaded to the CIBERSORT web portal (http://cibersort.stanford.edu/), running with the 1000 permutations and LM22 signature, which defines 22 immune cell subtypes including B cells, T cells, natural killer cells, dendritic cells, macrophages, and myeloid subsets. At a threshold of P < 0.05, the results of the inferred fractions of immune cell populations identified by CIBERSORT were represented to be statistically significant. For each sample, the final CIBERSORT output estimates were normalized and immune cell type fractions summed up to one.

Ethical statement

All information in this study was retrieved from public datasets; therefore, written informed consent was not necessary. This study meets the publication guidelines provided by the individual public datasets.

Abbreviations

AI: artificial intelligence; mPS: molecular prognostic score; ccRCC: clear cell renal cell carcinoma; AJCC: American Joint Committee on Cancer; TNM: tumor node metastasis; TCGA: The Cancer Genome Atlas; ICGC: International Cancer Genome Consortium; CNVs: copy number variations; UCISS: University of California Integrated Staging System; HR: hazard ratio; CI: confidence interval; GSEA: gene set enrichment analysis; TMB: tumor mutation burden; MCRs: minimal common regions; GISTIC: Genomic Identification of Significant Targets in Cancer.

Author Contributions

QP performed the computational analysis and drafted the manuscript. YS and KF participated in the statistical analysis for constructing the mPS system. ZD and LJ participated in evaluation and exploration of the potential mechanisms of the mPS system. DY and JZ conceived of the study and participated in its design and coordination. All authors read and approved the final manuscript.

Conflicts of Interest

The authors have declared that they have no conflicts of interest.

Funding

This study was supported by the National Natural Science Foundation of China (grant no. 81773221), the Natural Science Foundation of Jiangsu Province (grant no. BK20161222), Suzhou Science and Technology Planed Projects (grant no. SS201857) and the grant for Key Young Talents of Medicine in Jiangsu (grant no. QNRC2016875).

References