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

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.

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.

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.

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).

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).

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.

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).

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 6E–6I.

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. (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).

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 7A–7C). 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.

Figure 7. 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 (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); (E–G) 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 7E–7G). 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.

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).

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:

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 a₁ and a₂, which represented alive and deceased, respectively. We employed a softmax function for each node, and designated y₂ (probability of death; that is, the a₂ 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.