A novel prognosis and drug-susceptibility predictor based inflammatory-related genes signature in prostate cancer

2.1 Data acquisition

In this study, we accessed publicly available data, including transcriptome profiling (RNA-seq) data, clinical information, immune subtypes, and stemness scores derived from DNA-methylation (DNAss) and mRNA (RNAss) analyses, encompassing 475 PCa samples and 52 para-cancerous samples. These datasets were retrieved from The Cancer Genome Atlas (TCGA) database, accessible at https://portal.gdc.cancer.gov/. For external validation, we utilized an additional cohort comprising 268 samples from the cBioPortal cohort [21, 22]. To identify relevant IRGs, we curated a set of 200 genes from the “HALLMARK_INFLAMMATORY_RESPO-NSE” gene set, sourced from the GSEA database, available at http://www.gsea-msigdb.org/ [23].

2.2 Screening and visualizing hub IRGs

Univariate Cox hazards regression analysis was performed on the acquired prognostic IRGs and identified potential hub IRGs using the “survival” R package, and the least absolute shrinkage and selection operator (LASSO) regression was used to refine this set of IRGs. To identify the common genes of interest, a Venn diagram was constructed. The visualization of hub IRGs was accomplished through a heatmap generated with the “pheatmap” R package. Furthermore, we generated a forest plot using the “forest plot” R package to present the hazard ratio (HR), 95% confidence interval (CI), and p-value for each variable.

2.3 Construction and validation of a risk assessment model

Subsequently, we established the risk assessment model through multivariate Cox regression analysis. Patient risk scores were calculated using the following formula: risk score = β1 × 1 + β2 × 2 + … + βi × i, where × i represents the expression level of genes, and βi signifies the corresponding coefficient obtained from the multivariate Cox regression analysis. Patients diagnosed with PCa in both the TCGA and cBioPortal cohorts were categorized into either the low-risk or high-risk group based on their risk scores, which were determined using the median value. To assess the prognostic distinction between these two groups, we conducted Kaplan-Meier survival analysis. Furthermore, we evaluated the predictive performance by calculating the area under the time-dependent receiver operating characteristic (ROC) curve (AUC) using the “time ROC” R package.

2.4 Prognostic nomogram and clinical characteristics value evaluation

To explore the relationship between IRGs and PCa patients, we developed a genomic nomogram using the identified IRGs to predict the 1-, 3- and 5-year survival probabilities for each patient with the “rms” R package. Additionally, we performed univariate and multivariate Cox Hazards regression analyses to determine the significance of the risk score, along with other potential prognostic markers such as age and TNM (Tumor node metastasis) classification; T: Tumor; N: Lymph Node; M: Metastasis), in predicting the prognosis of PCa patients.

2.5 Investigations of the tumor microenvironment (TME), function enrichment analysis, DNA methylation pattern and mRNA expression

Gene Set Enrichment Analysis (GSEA) was conducted to examine differences in enrichment between the high-risk and low-risk groups concerning the activity of immune-related pathways. To assess and quantify variations in immune cell activity and immune function between these groups, we utilized single-sample gene set enrichment analysis (ssGSEA) with the “GSEABase” and “GSVA” R packages [24, 25]. Furthermore, we employed the “ESTIMATE” R package to compute the immune score, stromal score and ESTIMATE score [26].

2.6 Drug sensitivity analysis

To investigate the potential correlation between IRGs and the sensitivity of chemotherapeutic drugs, we acquired the NCI60 drug response data using the CellMiner tool at https://discover.nci.nih.gov/cellminer. This dataset included information from 60 distinct cell lines representing 9 different types of malignancies [27]. Then, Pearson correlation analysis was performed to assess the relationship between IRGs and the efficacy of chemotherapy drugs that have received approval from the Food and Drug Administration (FDA) or are currently undergoing clinical trials.

2.7 Statistical analysis

All statistical analyses were conducted using R version 4.0.3 (https://www.R-project.org/). To evaluate the association between risk scores and cancer stemness scores, Spearman’s test was employed. Pearson’s test was utilized to assess the correlation between gene expression and drug sensitivity. Additionally, the Wilcoxon rank-sum test, a commonly used nonparametric statistical test for comparing two groups, was applied. A significance level of p < 0.05 was considered statistically significant.

3.1 Identification of prognostic IRGs

Recognizing the significant influence of the inflammatory response on the development and progression of PCa, we compiled a dataset of human inflammatory response genes from the GSEA database, comprising 200 genes referred to as IRGs. To comprehensively understand the relevance of these genes in PCa, we initially conducted a univariate Cox proportional hazards regression analysis, which revealed 44 IRGs significantly associated with overall survival (p < 0.05) (Table 1). Subsequently, we applied LASSO regression analysis to the IRGs, utilizing 10-fold cross-validation to identify optimal tuning parameter values. As shown in Fig. 1A, all coefficients retained non-zero values, and 22 IRGs were stable (Fig. 1B). Then, a total of 17 IRGs were identified as common factors, represented in the Venn diagram (Fig. 1C), and their expression patterns of these 17 IRGs in PCa are shown through a heatmap (Fig. 1D).

3.2 Construction and validation of a risk assessment model in PCa

Next, all the factors derived from the LASSO-Cox regression analysis were included in the multivariate Cox regression analysis (Fig. 1E), based on which a model was established involving 17 IRGs, and a risk score formula was derived as follows: Risk score = Exp (BTG2) × (−0.3188) + Exp (C3AR1) × (0.6855) + Exp (CCL24) × (0.3293) + Exp (CCRL2) × (−1.5837) + Exp (CD14) × (0.3979) + Exp (CLEC5A) × (−0.7633) + Exp (DCBLD2) × (−0.7625) + Exp (GABBR1) × (0.2648) + Exp (KCNJ2) × (0.3135) + Exp (OPRK1) × (0.2719) + Exp (OSM) × (0.5280) + Exp (SCARF1) × (0.8879) + Exp (SCN1B) × (0.5578) + Exp (SEMA4D) × (0.5921) + Exp (SLC7A1) × (−0.3884) + Exp (SPHK1) × (0.5419) + Exp (TNFRSF1B) × (−0.6780) (Table 2). The correlations among these 17 IRGs are shown in Fig. 1F. Subsequently, a nomogram was constructed using these 17 IRGs to predict the 1-, 3- and 5-year survival of PCa patients (Fig. 2A). The cumulative score of each gene was used to calculate the total score, and vertical lines were drawn downwards at the corresponding positions on the total score axis to determine the relative survival rates at 1-, 3- and 5-year intervals. Based on the calculated median risk score, patients were categorized into two groups: the low-risk and high-risk groups. Kaplan-Meier survival curves demonstrated a statistically significant difference between the two cohorts, indicating that those with high-risk scores had a less favorable prognosis compared to those with low-risk scores (Fig. 2B). Subsequently, a time-dependent ROC analysis was conducted to assess the prognostic performance of the IRGs. The AUC for the prognostic model at 1-, 3- and 5-year intervals were determined to be 0.856, 0.764 and 0.828, respectively (Fig. 2C). To validate the prognostic power of the IRGs, an external validation cohort from the cBioPortal dataset was utilized. Consistently, patients classified as high-risk exhibited significantly shorter overall survival compared to those classified as low-risk (Fig. 2D). The AUC values for 1-, 3- and 5-year intervals were 0.756, 0.768 and 0.787, respectively (Fig. 2E). Additionally, TCGA-PCa patients were divided into high-risk and low-risk groups based on the median risk score, which served as the threshold (Fig. 2F). A high-risk score was associated with an increased likelihood of mortality, as evident from the distribution of IRG risk scores and patient survival status (Fig. 2G). These results suggest that the risk score obtained from these 17 IRGs demonstrates a high degree of accuracy and substantial predictive value for assessing the overall survival of PCa patients.

Fig. 1.Prognostic inflammatory-related gene identification process in prostate cancer. (A,B) Key inflammatory-related genes were selected using LASSO regression analysis. (C) A Venn plot highlights the 17 common prognostic genes associated with inflammatory responses. (D) The heatmap displays the expression of these 17 genes in PCa and normal tissues. (E) A forest plot shows hazard ratios for the 17 prognostic genes, forming the basis for a prognostic model. (F) The correlation plot of signature genes, with red indicating a positive connection and blue indicating a negative association. LASSO: Least Absolute Shrinkage and Selector Operation; PCa: prostate cancer.

Fig. 2.Construction and evaluation of the prognostic model using hub IRGs. (A) A nomogram designed for predicting 1-, 3- and 5-year OS of PCa patients. Vertical lines extend upward to determine values for each variable. The total score is obtained by summing individual variable values, and the anticipated likelihood of OS is determined by extending a vertical line downward from the Total Points axis. (B) Kaplan-Meier curves illustrating the outcomes of low- and high-risk groups in TCGA. (C) Kaplan-Meier curves showing the outcomes of low- and high-risk groups in the cBioPortal cohorts. (D) ROC curves assessing the predictive ability of the risk score for 1-, 3- and 5-year OS in TCGA. (E) ROC curves assessing the predictive ability of the risk score for 1-, 3- and 5-year OS in the cBioPortal cohorts. (F) PCa patients categorized into high- and low-risk groups based on the risk score. (G) Distribution of IRG risk scores and patient survival status. OS: overall survival; PCa: prostate cancer; TCGA: The Cancer Genome Atlas; ROC: receiver operating characteristic; IRG: inflammatory-related genes.

3.3 Independent prognostic value of the IRG risk score

We included the risk score and clinical factors in both univariate and multivariate analyses to assess the potential of the IRGs as independent predictors. In the TCGA population, pathological (T) stage, PSA value and risk score were identified as having a significant association with overall survival (OS) based on both univariate and multivariate Cox regression models (Fig. 3A,B). In the univariate analysis (p < 0.001), the hazard ratio (HR) for the risk score was 1.140 with a 95% confidence interval (CI) of 1.107–1.175. In the multivariate analysis (p < 0.001), the HR for the risk score was 1.128 with a 95% CI of 1.094–1.164. Similarly, in the cBioPortal cohort, univariate and multivariate Cox regression analyses revealed a significant correlation between OS and T stage as well as the risk score (Fig. 3C,D). In the univariate analysis (p = 0.002), the HR for the risk score was 1.002 with a 95% CI of 1.001–1.004. In the multivariate analysis (p = 0.005), the HR for the risk score was 1.002 with a 95% CI of 1.001–1.003. Collectively, these findings provide strong evidence that the IRG risk score serves as a valuable independent prognostic factor for PCa.

Fig. 3.Univariate and multivariate Cox regression analyses of the study cohort. (A,B) Univariate and multivariate Cox regression analyses for the TCGA dataset, assessing the impact of various clinicopathological factors and the risk score on prognosis. (C,D) Univariate and multivariate Cox regression analyses for the cBioPortal dataset, evaluating the influence of different clinicopathologic factors and the risk score on prognosis. PSA: prostate-specific antigen; TMB: tumor mutation burden.

3.4 Tumor microenvironment assessment using the IRG risk score

We used single-sample gene set enrichment analysis (ssGSEA) to quantify 16 immune cell subsets and 13 immune-related functions and evaluate the potential of the IRG risk score in assessing immunological features and elucidating the relationship between the risk score and overall immune status. The analysis of immune cell infiltration revealed a significant positive correlation between the risk score of IRGs and the abundance of various immune infiltrating cells, including tumor-infiltrating lymphocytes (TIL), CD8+ (cluster of differentiation 8) T cells, dendritic cells (DCs), macrophages, plasmacytoid dendritic cells (pDCs), and T helper cells (Fig. 4A). Furthermore, the analysis of immune-related pathways indicated that immunological responses were significantly more pronounced in the high-risk group compared to the low-risk group (Fig. 4B). We also examined the immunophenotyping distribution of different tumor sample types in the TCGA database, revealing risk scores for four distinct immune types (Fig. 4C). To gain a deeper understanding of the underlying functionality of IRGs and their associated signal transduction pathways, we conducted an extensive investigation using Gene Set Enrichment Analysis (GSEA). The results demonstrated significant differences in pathway enrichment between the low-risk and high-risk groups (Fig. 4D). Notably, the high-risk group exhibited positive associations with immune pathways such as cytolytic activity, HLA (human leukocyte antigen), MHC (major histocompatibility complex) class I, and others. These findings suggest that PCa patients with high-risk scores may be more responsive to immunotherapy due to their heightened immune regulation. This insight may offer a new perspective to enhance the effectiveness of immunotherapy in the context of PCa treatment.

Considering the TME and the presence of stromal cells and immune cells play significant roles in cancer, we conducted a correlation analysis to explore the relationship between the risk score and the TME, aiming to better understand the impact of TME on PCa patients. The results reveal a significant positive correlation between the risk score and the presence of stromal cells (p < 0.05) as well as immune cell infiltration (p < 0.05) (Fig. 4E,F). Furthermore, we conducted a Spearman correlation analysis to investigate the association between the risk score and cancer stemness scores based on the stem cell score derived from DNA methylation (DNAss) and RNA sequencing (RNAss) (Fig. 4G,H). The findings demonstrate a statistically significant positive relationship between DNAss and the risk score (p < 0.05). Consequently, these results suggest a potential and significant correlation between the risk score of the prognostic model and the activity level of cancer stem cells, providing valuable insights into the role of cancer stemness in PCa.

Fig. 4.Correlation analysis between the risk score and immune infiltration status. (A) Comparison of immune cell subsets in low- and high-risk groups. (B) Comparison of immune function and pathways in the low- and high-risk groups. (C) Differences in immune classification among PCa patients with different risk scores. (D) Functional enrichment in risk groups revealed by GSEA analysis. (E) Scatterplot depicting the correlation between the stromal cell score. (F) Scatterplot depicting the correlation between the immune cell score. (G) Scatterplot depicting the correlation between DNAss. (H) Scatterplot depicting the correlation between RNAss. (*p < 0.05, **p < 0.01, ***p < 0.001). PCa: prostate cancer; GSEA: Gene Set Enrichment Analysis; DNAss: DNA methylation-based stemness scores; RNAss: RNA expression-based stemness scores.

3.5 The relationship between IRGs and drug sensitivity

We utilized the NCI-60 database, which comprises 60 distinct human cancer cell lines, to investigate the relationship between the expression levels of signature genes and the susceptibility of these cell lines to various drugs (Fig. 5). The study findings revealed a positive correlation between elevated levels of C3AR1 (Complement Component 3a Receptor 1), CLEC5A (C-Type Lectin Domain Family 5 Member A), OSM (Oncostatin M), SEMA4D (Semaphorin 4D) and SCARF1 (Scavenger Receptor Class F Member 1) and the sensitivity of cancer cells to several chemotherapeutic drugs, including Denileukin Diftitox Ontak, Artemether, ABT-199, Nelarabine, Methylprednisolone, Fluphenazine, Zalcitabine, Fludarabine and Ribavirin. Notably, the expression of SCARF1 demonstrated a positive association with the sensitivity of cancer cells to Denileukin Diftitox Ontak, Fludarabine, Methylprednisolone, Nelarabine, Zalcitabine, and Ribavirin. Conversely, the upregulation of C3AR1 was linked to increased resistance to Irofulven in cancer cells (Fig. 5). Table 3 shows the correlation analysis results, indicating significant associations with a p-value of < 0.001. These findings provide empirical evidence supporting the efficacy of the risk score in accurately predicting the responsiveness of cancer cells to these drugs, supporting the potential use of our proposed model in clinical settings to provide guidance in enhancing the precision of drug selection and treatment strategies.

Fig. 5.Investigation of gene-drug sensitivity using the cellminer database. In this figure, we explored gene-drug sensitivity utilizing the CellMiner database and identified the top 16 medications that exhibit the strongest connection with gene expression in an inflammatory-related prognostic model. Cor: correlation.