Computational construction of TME-related lncRNAs signature for predicting prognosis and immunotherapy response in clear cell renal cell carcinoma

Abstract

Background: The tumor microenvironment (TME) is closely related to clear cell renal cell carcinoma (ccRCC) prognosis, and immunotherapy response. In current study, comprehensive bio-informative analysis was adopted to construct a TME-related lncRNA signature for immune checkpoint inhibitors (ICIs) and targeted drug responses in ccRCC patients.

Methods: The TME mRNAs were screened following the immune and stromal scores with the data from GSE15641, GSE29609, GSE36895, GSE46699, GSE53757, and The Cancer Genome Atlas (TCGA)-kidney renal clear cell carcinoma (KIRC). And the TME-related lncRNAs were recognized using correlation analysis. The TME-related lncRNAs prognostic model was constructed using the training dataset. Kaplan-Meier analysis, principal-component analysis, and time-dependent receiver operating characteristic were used to evaluate the risk model. The immune cell infiltration in TME was evaluated using the single-sample gene set enrichment analysis (ssGSEA), ESTIMATE, and microenvironment cell populations counter algorithm. The immunophenoscore (IPS) was used to assess the response to immunotherapy with the constructed model.

Results: In the current study, 364 TME-related lncRNAs were selected based on the integrated bioinformatical analysis. Six TME-related lncRNAs (LINC00460, LINC01094, AC008870.2, AC068792.1, and AC007637.1) were identified as the prognostic signature in the training dataset and subsequently verified in the testing and entire datasets. Patients in the high-risk group exhibited poor overall survival and disease-free survival than those in the low-risk group. The 1-, 3-, and 5-year areas under the curves of the prognostic signature in the entire dataset were 0.704, 0.683, and 0.750, respectively. The risk score independently predicted ccRCC survival based on univariate and multivariate Cox regression. GSEA analysis suggested that the high-risk group was concentrated on immune-related pathways. The high-risk group were characterized by high immune cell infiltration, high TMB and somatic mutation counters, high IPS-PD-1 + CTLA4 scores, and immune checkpoints expression upregulation, reflecting the higher ICIs response. The half inhibitory concentrations of sunitinib, temsirolimus, and rapamycin were low in the high-risk group.

Conclusion: The TME-related lncRNAs signature constructed could reliably predict the prognosis and immunotherapy response and targeted ccRCC patients' therapy.

Keywords: clear cell renal cell carcinoma; immunotherapy; lncRNA; prognostic signature; tumor microenvironment.

FIGURE 1

Identification of TME‐related lncRNAs in KIRC. Volcano plot (A) and heatmap (B) of DRGs between high and low immune scores. Volcano plot (C) and heatmap (D) of DRGs between high and low stromal scores. (E) The scale‐free fit index (left) and the mean connectivity (right). (F) Heatmap of the relationships between the state of immune and model eigengenes. (G) Venn diagram of TME‐related mRNAs. (H) The results of gene ontology analysis. (I) The top 15 most significant KEGG pathways

FIGURE 2

Construction of the TME‐related lncRNAs risk model. (A) Volcano plot showing that 203 upregulated genes (red) and 29 downregulated genes (blue) between KIRC and normal kidney specimens. (B) Heatmap showing the expression levels of the top 20 upregulated and downregulated TME‐related lncRNAs. (C) The optimal values of the penalty parameter. Univariate (D) and multivariate (E) Cox regression analysis of the six TME‐related lncRNAs in the risk model

FIGURE 3

Validating the performance of the risk model. (A) PCA analysis for the risk model in the training dataset. Distribution of risk score (B) and survival status (C) in the training dataset. Kaplan–Meier curves of OS (D) and DFS (E) for the training dataset. The AUCs of the time‐dependent ROC curves for the training dataset (F). Kaplan–Meier curves of OS and DFS for the testing dataset (G, H) and entire dataset (I, J). The AUCs of the time‐dependent ROC curves for the testing dataset (K) and entire dataset (L)

FIGURE 4

Assessment of the independent prognostic value of the risk model. Univariate (A) and multivariate (B) Cox regression analyses of the risk score and clinical characteristics. (C) Distribution landscape of clinical characteristics and the expression profiles of six TME‐related lncRNAs between the high‐ and low‐risk groups. Discrepancies in risk scores by grade (D), stage (E), and T (F). (G) The nomogram combining the risk score with the age in predicting 1‐,3‐, and 5‐year OS for KIRC patients given a total risk score. (H) Calibration curves. ***p < 0.001; **p < 0.001; *p < 0.05

FIGURE 5

Gene set enrichment analyses between the high‐ and low‐risk groups. GO enrichment analyses in the high‐risk group (A) and the low‐risk group (B). KEGG pathway analyses in the high‐risk group (C) and the low‐risk group (D)

FIGURE 6

The relationships between immune cells and the risk signature. (A) The heatmap of 29 immune‐related cells and functions, immune score, stromal score, ESTIMATE score, and tumor purity between two risk groups. (B‐E) The differences of immune score, stromal score, ESTIMATE score and tumor purity between two risk groups. (F) The scores of immune cells comparing high‐ and low‐risk groups by ssGSEA. (G) The scores of immune functions comparing high‐ and low‐risk groups by ssGSEAS. The correlation of the immune cells and the risk scores by MCP counter. (I‐N) The differences of immune cells between high‐ and low‐risk groups by MCP counter. ***p < 0.001; **p < 0.01; *p < 0.05; ns: no significance

FIGURE 7

Somatic mutational analyses between high‐ and low‐risk groups. (A) Difference in TMB between the high‐ and low‐risk groups. (B) Correlation between TMB and risk scores. (C) Survival analysis of OS between H‐ and L‐TMB groups. (D) Survival analysis of OS stratified by TMB and risk groups. (E) Difference in somatic mutation counts between the high‐ and low‐risk groups. (F) Correlation between somatic mutation counts and risk scores. (G) Survival analysis of OS between BAP1 mutation and wild groups. (H) Survival analysis of OS stratified by BAP1 mutation status and risk groups. The top 20 frequently mutated genes in the high‐risk group (I) and low‐risk group (J). (K) Proportion of BAP1 mutation status in high‐ and low‐risk groups. (L) Correlation between the risk score and common immune checkpoints. (M) Expression levels of the common immune checkpoints between the high‐ and low‐risk groups. ***p < 0.001; *p < 0.05; ns: no significance

FIGURE 8

Response to immunotherapy and sensitivity to targeted therapy between the high‐ and low‐risk groups. (A‐D) The association between IPS and the risk signature of KIRC patients. (E) Correlation between the risk score and common immune checkpoints. (F) Expression levels of the common immune checkpoints between the high‐ and low‐risk groups. IC50 values between the high‐ and low‐risk groups for axitinib (G), sunitinib (H), sorafenib (I), rapamycin (G), pazopamib (K), and emsirolimus (L). ***p < 0.001; *p < 0.05; ns: no significance