1. Identification of Differentially Expressed IMRGs and ccRCC Molecular Subtypes.
The TCGA-KIRC database was utilized to obtain transcriptome and clinical information of 539 ccRCC samples and 72 normal kidney samples. The quantification of gene expression levels pertaining to immune function (as elucidated in Supplementary Table 1) and metabolic processes (as explicated in Supplementary Table 2) has been meticulously extracted. After carrying out Wilcoxon signed-rank test, a comparative analysis was conducted between tumor samples and normal samples. Based on the criteria of |log2FC| > 2 and FDR < 0.05 (Fig. 1A, 1B), following rigorous analysis, a comprehensive set of 519 differentially expressed immune- and metabolism-related genes (IMRGs), including 374 overexpressed and 145 underexpressed genes (as listed in Supplementary Table 3), were successfully identified. Furthermore, employing univariate Cox regression, a subset of 125 IMRGs that were deemed to be prognostically relevant for the studied condition were also detected (detailed in Supplementary Table 4).
Subsequently, the NMF clustering algorithm was applied to the prognosis-related IMRGs, resulting in the identification of two molecular subtypes (Fig. 1C). The indicators of cophenetic, silhouette, and dispersion were mainly used to determine the optimal rank value (Supplementary Figure S1). The heatmap demonstrated the expression of the prognosis-related IMRGs in both molecular subtypes (Fig. 1D). On the basis of the statistical analysis, the survival outcomes in terms of overall survival (OS) and progression-free survival (PFS) were significantly superior for cluster 1 when compared to cluster 2, as illustrated by Fig. 1E and 1F. A further investigation revealed that cluster 1 exhibited substantially lower Immune Score, Stromal Score, and ESTIMATE Score values than those of cluster 2, as illustrated by Fig. 1G. Additionally, results from the cell composition analysis depicted that cluster 1 manifested significantly lower levels of B cells, T cells, and fibroblasts as compared to cluster 2. In contrast, cluster 1 demonstrated higher proportions of endothelial cells and neutrophils than those observed in cluster 2 (as illustrated in Fig. 1H). It is suggested that the higher infiltration of immune cells in cluster 2, but worse prognosis, may be related to the highly immune escape state in these tumors. These findings suggest that the two clusters may have distinct biological characteristics that could influence their prognosis and response to therapy.
2. Construction of a Nine-Gene-Based IMRGs Prognostic Signature
The TCGA-KIRC cohort underwent a random division into two groups with a ratio of 7:3, which were referred to as the training cohort and test cohort. Utilizing IMRGs that were associated with prognosis, TCGA training cohort was utilized for the purpose of developing a prognostic prediction model. The construction of this model involved implementing LASSO Cox regression analysis, which allowed for the identification of key predictors and the development of an accurate prognostic model. Supplementary Figure S2 presents a comprehensive visualization of the coefficients attributed to each independent variable in the LASSO regression, along with the optimal logarithmic value of lambda. It was attempted to develop a prognostic signature comprising of nine genes, including CTSE, KLRC2, PDIA2, HAMP, PGLYRP2, ORM2, CHGA, UCN, which are immune-related genes, and ADCY2, a metabolism-related gene. Through the implementation of multivariate Cox regression analysis, a risk score (RS) was calculated to represent an all-encompassing index of immune and metabolic status for the IMRGs signature. In order to obtain the required score, it was attempted to employ a method that involved multiplication of the expression levels of each gene by their corresponding coefficients (Fig. 2A):
RS = -0.139435 × CTSE + 0.115443 × KLRC2 + 0.125361 × PDIA2 + 0.198961 × HAMP + 0.126838 × PGLYRP2 + 0.183848 × ORM2 + 0.103308 × CHGA + 0.205241 × UCN + -0.110254 × ADCY2
By utilizing the median of RS, patients from the TCGA training cohort were partitioned into two groups based on their level of risk to investigate the correlation between RS and prognosis. Results from survival analysis displayed that patients categorized as high-risk had considerably worse prognosis compared to those in the low-risk group (Fig. 2B). In Fig. 2C, it can be observed that the prognostic model's AUCs were evaluated in the entire TCGA-KIRC cohort, with scores of 0.813, 0.751, and 0.779 for 1-, 3-, and 5-year, respectively. These outcomes suggest that the prognostic model has a commendable prognostic value. Additionally, the heatmap depicted in Fig. 2D demonstrated the differential expression profiles of risk genes in the IMRGs signature, highlighting their relevance to both high-risk and low-risk groups. Furthermore, the scatter plot illustrated the distribution of RS and its association with survival outcomes, providing insights into the potential clinical applications of RS in cancer prognosis (Fig. 2E-F). Supplementary Figure S2 will include the aforementioned results for both training and test cohorts.
3. Prognostic Differences Between High- and Low-Risk Groups Unearthed Through Subgroup Analysis
To determine whether the signature of IMRGs was associated with clinical features, it was attempted to employ independent t-tests for the purpose of evaluating the differences in RS in TCGA-KIRC cohort. As the level of age (Fig. 3A), pathological stage (Fig. 3B), Fuhrman grade (Fig. 3C), tumor size (Fig. 3D), lymph node metastases (Fig. 3E), and distant metastases (Fig. 3F) increased, so did RS. Irrespective of the clinical features such as age (Fig. 3G), pathological stage (Fig. 3H), Fuhrman grade (Fig. 3I), and TNM as described above (Fig. 3J-3L), significant prognosis-based differences between high-risk and low-risk groups were identified through subgroup analysis. However, in patients with lymph node metastasis, the absence of a significant difference between the high- and low-risk groups could be attributed to the limited sample size. These results demonstrate that the prognostic signature of IMRGs has robust prognostic predictive capabilities across different clinical characteristics.
4. An evaluation of the prognostic signature of IMRGs through a comparative analysis with signatures published beforehand
To investigate the predictive performance of the model that associates with immune system and metabolism, we conducted a comparison with ten prognostic signatures that were previously published within the last three years. These included five immune-related prognostic models, namely Immune-related_Cao, Immune-related_Lin, Immune-related_Liu, Immune-related_Wang, and Immune-related_Zhang, as well as five metabolism-related prognostic models, involving Metabolism-related_Guo, Metabolism-related_Huang, Metabolism-related_Tan, Metabolism-related_Wang, and Metabolism-related_Wei [14–23]. In order to ensure comparability of the signatures, we utilized the same method for computing and converting the Risk Scores (RS) across TCGA-KIRC cohort. The ability to classify ccRCC samples into high- and low-risk groups was observed in all ten published signatures, with statistically significant differences noted in outcomes for both groups (as depicted in Figs. 4A). However, after conducting ROC curve analysis, it was found that the AUCs of our model were superior to those of the ten published signatures. Specifically, our model's performance was evaluated based on the calculation of AUCs for survival rates at 1 year, 3 years, and 5 years. The resulting AUC values were 0.813, 0.751, and 0.779 respectively, outperforming the published signatures (refer to Figs. 4B). Additionally, our model had the highest C-index value at 0.741, in contrast to the other ten models that exhibited values ranging between 0.578 and 0.697 (as shown in Fig. 4C). These findings consistently demonstrated the superior performance of the IMRGs prognostic signature.
5. Development of a Nomogram Utilizing the IMRGs Prognostic Signature and Assessment of Its Clinical Significance
To figure out the applicability of the IMRGs prognostic signature for clinical purposes, it was attempted to conduct Cox regression analysis on the complete TCGA-KIRC cohort to assess its independence. The univariate and multivariate regression analyses demonstrated significant correlations between RS and prognosis, underscoring the potential value of RS as a prognostic marker for predicting clinical outcomes, as illustrated in Fig. 5A. It was also observed that age, stage, and grade could serve as independent prognostic factors, along with gender. Therefore, gender was not included in the subsequent construction of the nomogram.
A predictive tool in the form of a nomogram was constructed to estimate the likelihood of survival risk, which is depicted in Fig. 5B. To calculate the total score for each patient, the points assigned to each variable on the point scale were added together. On the basis of total points calculated in the previous step using the bottom scale, the likelihood of survival at 1-, 3-, or 5-year intervals was predicted. To ascertain the agreement between the survival predicted by the nomogram and the actual survival, calibration curves were employed (as presented in Fig. 5C). In addition, the AUCs for the nomogram were found to be higher than separate RS and other clinical variables (refer to Figs. 5D-F). In light of the findings, it can be concluded that the nomogram developed using the IMRGs prognostic signature was effective in predicting disease progression and had a substantial association with prognosis.
6. Assessing the Predictive Value of IMRGs Prognostic Signature in Immunotherapy Response.
In order to investigate the influence of the IMRGs prognostic signature on the effectiveness of immunotherapy, an examination was undertaken to indicate the associations between RS and the infiltration of immune cells within TME. Figure 6A displays a significant decline in the proportions of B lineage, CD8 + T cells, cytotoxic lymphocytes, fibroblasts, monocytic lineage NK cells, and T cells in the low-risk group. Conversely, there was a significant increase in the proportions of endothelial cells and neutrophils in this group. In light of the data presented in Fig. 6B, it can be inferred that there exist notable distinctions in the expression levels of immune checkpoints, specifically PD-1, CTLA4, TIGIT, and LAG3 [24–27], between the high-risk and low-risk groups. The high expression of inhibitory immune checkpoints in the high-risk group suggests a stronger immune escape. Hence, Fig. 6C illustrates our examination of the interdependence between two crucial prognostic indicators for immunotherapy response, namely RS and IPS. When using PD1 or CTLA4 inhibitors alone, after completion of the analysis on high-risk and low-risk groups, no noticeable difference could be identified. However, when both PD1 and CTLA4 were blocked simultaneously, significant differences between the high- and low-risk groups emerged, reflecting the clinical value of using ICI in ccRCC in combination. These findings are in line with earlier studies that have reported positive therapeutic outcomes with combined Nivolumab (PD1 inhibitor) and Ipilimumab (CTLA4 inhibitor) treatment, which has received FDA approval as a first-line treatment for medium to high-risk advanced RCC [28]. Thus, the outcomes demonstrated the potential of the prognostic signature of IMRGs to indicate the status of immune infiltration and anticipate the reaction to immunotherapy.