Construction of a novel gene-based model for prognosis prediction of clear cell renal cell carcinoma

Abstract

Background: Clear cell renal cell carcinoma (ccRCC) comprises the majority of kidney cancer death worldwide, whose incidence and mortality are not promising. Identifying ideal biomarkers to construct a more accurate prognostic model than conventional clinical parameters is crucial.

Methods: Raw count of RNA-sequencing data and clinicopathological data were acquired from The Cancer Genome Atlas (TCGA). Tumor samples were divided into two sets. Differentially expressed genes (DEGs) were screened in the whole set and prognosis-related genes were identified from the training set. Their common genes were used in LASSO and best subset regression which were performed to identify the best prognostic 5 genes. The gene-based risk score was developed based on the Cox coefficient of the individual gene. Time-dependent receiver operating characteristic (ROC) and Kaplan-Meier (KM) survival analysis were used to assess its prognostic power. GSE29609 dataset from GEO (Gene Expression Omnibus) database was used to validate the signature. Univariate and multivariate Cox regression were performed to screen independent prognostic parameters to construct a nomogram. The predictive power of the nomogram was revealed by time-dependent ROC curves and the calibration plot and verified in the validation set. Finally, Functional enrichment analysis of DEGs and 5 novel genes were performed to suggest the potential biological pathways.

Results: PADI1, ATP6V0D2, DPP6, C9orf135 and PLG were screened to be significantly related to the prognosis of ccRCC patients. The risk score effectively stratified the patients into high-risk group with poor overall survival (OS) based on survival analysis. AJCC-stage, age, recurrence and risk score were regarded as independent prognostic parameters by Cox regression analysis and were used to construct a nomogram. Time-dependent ROC curves showed the nomogram performed best in 1-, 3- and 5-year survival predictions compared with AJCC-stage and risk score in validation sets. The calibration plot showed good agreement of the nomogram between predicted and observed outcomes. Functional enrichment analysis suggested several enriched biological pathways related to cancer.

Conclusions: In our study, we constructed a gene-based model integrating clinical prognostic parameters to predict prognosis of ccRCC well, which might provide a reliable prognosis assessment tool for clinician and aid treatment decision-making in the clinic.

Keywords: Clear cell renal cell carcinoma (ccRCC); Differentially expressed genes (DEGs); GEO; Nomogram; Overall survival (OS); Risk score; TCGA.

Fig. 1

Flowchart of the whole study

Fig. 2

Data processing, screening of the DEGs. a A PCA plot of the data showing no batch effect in the TCGA KIRC dataset. Red nodes represent the normal cluster while blue nodes represent the tumor cluster. b Heatmap of top 20 DEGs in ccRCC. c Volcano plot of differentially expressed genes in ccRCC when compared with normal tissue. Red nodes represent the significantly up-regulated genes with logFC > 4 and p < 0.05. Green nodes represent the significantly down-regulated genes with logFC < -4 and p < 0.05. PCA, principal component analysis; TCGA, The Cancer Genome Atlas; DEGs, differentially expressed genes; ccRCC, clear cell renal cell carcinoma

Fig. 3

Identification of 5 significantly prognostic genes and their expression data in ccRCC. a Venn diagram of 40 OCGs. 9 upregulated genes with HR > 1. 31 downregulated genes with HR < 1. b LASSO coefficients profiles of 19651 protein-coding genes. c LASSO regression with tenfold cross-validation obtained 14 prognostic genes using minimum lambda value. d Multivariate Cox regression analysis of 5 prognostic genes from BSR. e Expression pattern of the five genes between tumor and normal kidney tissue. f IHC of the five genes in tumor and normal kidney tissue. ccRCC clear cell renal cell carcinoma, OCG overlapping candidate genes, HR hazard ratio, LASSO least absolute shrinkage and selection operator, BSR best subset regression, IHC immunohistochemistry. *p < 0.05, **p < 0.01, ***p < 0.001

Fig. 4

Kaplan–Meier survival analysis of PADI1, ATP6V0D2, C9orf135, DPP6, and PLG

Fig. 5

Prognostic analysis of five-gene signature in the training set. The dotted line represented the median risk score and divided the patients into low- and high-risk group. a The curve of risk score. b Survival status of the patients. More dead patients corresponding to the higher risk score. c Heatmap of the expression profiles of the five prognostic genes in low- and high-risk group. d Kaplan–Meier survival analysis of the five-gene signature. e Time-dependent ROC analysis the of the five-gene signature. ROC receiver operating characteristic

Fig. 6

Kaplan–Meier survival analysis of the five-gene risk score level in different subgroups including stage I/II, stage III/IV (a), grade 1/2, grade 3/4 (b), younger than 65 years old, older than 65 years old (c), male, female (d), left tumor, right tumor (e) and with or without recurrence (f)

Fig. 7

Validation of the five-gene signature. GSE29609 was regarded as the external validation set. Kaplan–Meier survival analysis of the five-gene signature in internal validation set (a), the whole set (b) and external validation set (c). Time-dependent ROC analysis of the five-gene signature in internal validation set (d), the whole set (e) and external validation set (f). ROC, receiver operating characteristic

Fig. 8

Identifying the independent prognostic parameters and construction of gene-based prognostic model. a Forrest plot of univariate Cox regression analysis in ccRCC. b Forrest plot of multivariate Cox regression analysis in ccRCC. c Nomogram integrated five gene-based risk score, AJCC-stage, grade and age

Fig. 9

Performance of gene-based nomogram in predicting survival probability and comparison of the predictive power among gene-based nomogram, risk score and AJCC-stage. The calibration plot of the nomogram for agreement test between 1-, 3- and 5-year OS prediction and actual outcome in the training set (a), the internal validation set (e) and the entire set (i). The time-dependent ROC curves of the nomogram, risk score and AJCC-stage in 1-, 3- and 5-year OS prediction in the training set (b–d), the internal validation set (f–h) and the entire set (j–l). OS overall survival, ROC receiver operating characteristic

Fig. 10

Functional enrichment analysis of 399 DEGs. The more genes enriched in the terms, the darker the color. a Top 20 of GO enrichment analysis of the DEGs. b Top 20 of KEGG enrichment analysis of the DEGs. DEGs differentially expressed genes, GO gene ontology, KEGG Kyoto Encyclopedia of Genes and Genomes

Fig. 11

GSEA associated with the five genes expression. The gene set “GLYCOSAMINOGLYCAN_BIOSYNTHESIS_CHONDROTIN_SULFATE” (a), “LYSOSOME” (b), “ADHESION_JUNCTION” (b), and “GLYCOSAMINOGLYCAN_BIOSYNTHESIS_CHONDROTIN_SULFATE” (c) were enriched in ccRCC samples with highly expressed PADI1, DPP6 and ATP6V0D2, respectively. The gene set “P53_SIGNALING_PATHWAY” (d) “PPAR_SIGNALING_PATHWAY” (e) and were enriched in ccRCC samples with lowly expressed PLG and C9orf135, respectively