Cell Differentiation Trajectory Predicts Prognosis and Immunotherapeutic Response in Clear Cell Renal Cell Carcinoma

Abstract

Clear cell renal cell carcinoma (ccRCC) is the main type of malignancy in kidney related to glucose metabolism. Primary single cell culture and single cell sequencing are novel research technologies. In this study, we explored the differentiation status of ccRCC cells and its significance in prognosis and immunotherapeutic response through bioinformatics. We characterized distinct differentiation states and differentiation-related genes (DRGs) in ccRCC cells through single cell RNA sequencing (scRNA-seq) analysis. Combined with bulk RNA-seq data, we classified patients into two clusters and found that this classification was closely correlated with patient prognosis and immunotherapeutic responses. Based on machine learning, we identified a prognostic risk model composed of 14 DRGs, including BTG2, CDKN1A, COL6A1, CPM, CYB5D2, FOSB, ID2, ISG15, PLCG2, SECISBP2, SOCS3, TES, ZBTB16, and ZNF704, to predict the survival rate of patients and then constructed a nomogram model integrating clinicopathological characteristics and risk score for clinical practice. In the study of immune checkpoints, we found that patients in the high-risk group had a disposition to get worse prognosis and better effects of immune checkpoint blocking therapies. Finally, we found the expression level of model DRGs was associated with a tumor-immune microenvironment (TIME) pattern and the response of 83 compounds or inhibitors was significantly different in the two risk groups. In a word, our study highlights the potential contribution of cell differentiation in prognosis judgment and immunotherapy response and offers promising therapeutic options for ccRCC patients.

Figure 1

Preprocessing of the scRNA-seq data: (a) violin plots of the RNA information of processed scRNA-seq data, (b) scatter plot of the correlation between the numbers of detected genes and sequencing depth, (c) scatter plot of the batch effect after correction, and (d) scatter plot of 1,000 highly variable genes.

Figure 2

Cell clustering and trajectory analysis based on the scRNA-seq data: (a) scatter plot of 23 clusters processed by the UMAP algorithm, (b) scatter plot of 10 cell types obtained through annotation, (c–d) dot plot of the expression of major marker genes in different clusters and cell types, and (e–g) differentiation trajectory of 3 branches with diverse pseudotime, cell types, and states.

Figure 3

Classification constructed for ccRCC patients according to DRGs: (a–b) construction of classification through NMF algorithm, (c) Kaplan–Meier survival analysis for patients in 2 clusters, and (d) scatter plot of the classification of patients through PCA.

Figure 4

Comparisons of clinicopathological characteristics and differentially analysis of enriched functions between 2 clusters: (a) stacked histograms of the proportion of clinicopathological characteristics and (b) heatmaps of the diverse functional annotations.

Figure 5

Identification of diverse TIME, immune gene pattern, and immunotherapeutic response in 2 clusters: (a) violin plots of 4 indicators of TIME between C1 and C2, including immune score, stromal score, ESTIMATE score, and tumor purity; (b) K–M survival analysis for high and low TIME score ccRCC patients; (c) scatter plots of the correlation between immune, stromal score and tumor purity; (d) bar plots of the results of functional enrichment analysis; (e) box plots of the abundances of different infiltrating immune cells; (f) box plots of the expression levels of immune checkpoints; and (g) violin plots of the immunotherapy scores.

Figure 6

Construction and validation of the prognostic risk model based on DRGs: (a) nature of the network topology constructed with unique power values and the relationship between power values and average connectivity, (b) discrete modules of obtained through DRGs clustering, (c) module diagram of the correlation between clinicopathological characteristics and identified modules, (d) coefficient profile plot of the log (lambda) sequence of the LASSO model, (e–f) heatmaps of the expression levels of prognostic genes, curves of the risk score, and scatter plots of survival status in the training and validation cohort, and (g) box plots of the correlation between clinical variables and risk score.

Figure 7

Depth analysis according to the prognostic risk model: (a) bar plots showing the results of functional enrichment analysis based on DRGs, (b) K–M survival analysis of different risk patients in training and validation cohort, (c) ROC curves of the predictive efficiency of our model in training and validation cohort, and (d–f) ROC curves of the diverse predictive efficiency between our model and other published models in training and validation cohorts.

Figure 8

Development and validation of prognostic nomogram: (a–b) forest plots of the univariate and multivariate cox regression in training and validation cohorts, (b) nomogram composed of age, grade, stage, and risk score with the prediction of overall survival rate, (d–e) ROC and calibration curves of nomogram for the prediction of overall survival rate in training and validation cohorts.

Figure 9

Prediction of immunotherapeutic response and drug sensitivity: (a–b) heatmaps of correlation between the expression levels of 14 DGRs and the abundance of infiltrating immune cells in training and validation cohorts and (c–d) violin plots of immunotherapeutic response in different risk groups based on training and validation cohorts.