Comprehensive Analysis and Identification of an Immune-Related Gene Signature with Prognostic Value for Prostate Cancer

Introduction

Prostate cancer (PCa) is the second most frequent malignancy and the fifth leading cause of cancer death in men, with almost 1.4 million new cases and 375,304 deaths annually worldwide in 2020.¹ With the popularization of prostate-specific antigen screening, the rate of detection and diagnosis of PCa has substantially increased in the past 20 years. According to the clinical tumor stage, serum prostate-specific antigen level, and Gleason score, the risk stratification of PCa is classified as low-, intermediate-, and high-risk categories, which can be used to guide treatment and determine prognosis.^2–4 Current treatments for PCa include conservative treatment, surgery, endocrine therapy, radiotherapy, and chemotherapy.⁵ Despite improvements in clinical treatment strategies, the prognosis for PCa remains very poor. Approximately 25% of the patients with PCa suffer from biochemical recurrence within 5 years after radical prostatectomy.⁶ Moreover, up to one-third of patients with localized tumors eventually fail on local treatment and progress to advanced or metastatic PCa within 10 years.⁷ At present, there are few methods to evaluate and predict the prognosis of PCa. Therefore, the primary focus of the research strategy is on the urgent need to determine a non-invasive and effective method for predicting the prognosis of such patients.

Previous research has shown that T cells are the dominant lymphocytes in the normal prostate gland, whose function and CD4/CD8 location distribution are strictly regulated, and no significant difference is found in the immunosurveillance level of the prostate according to patient race, age, and anatomic zone.⁸ The infiltration of cytotoxic lymphocytes into the tumor microenvironment (TME) plays an important role in anti-tumor immunity. However, genome-wide studies have disclosed that PCa undergoes molecular alterations in cells of the TME throughout disease onset and progression.^9,10 Expression of negative regulatory molecules on the surface of tumor-infiltrating lymphocytes contributes to ineffective immune responses and the immunosuppressed phenotype. Thus, it is still necessary to define molecular regulators of the PCa microenvironment for facilitating the development of personalized targeted therapies.

The TME is a complex mixture that is composed of infiltrating immune cells, cancer-associated fibroblasts (CAFs), endothelial cells, and signaling molecules and extracellular matrix proteins.¹¹ Thereinto, CAFs are a key component of the TME with diverse functions, including crosstalk with infiltrating leukocytes, and matrix deposition and remodelling. Research has shown that CAFs are more genetically stable and less likely to acquire therapy resistance than tumor cells.¹² For this reason, precision strategies have been developed and are in progress to target CAFs.^13–15

To evaluate the infiltration degree of stromal and immune cells, ESTIMATE algorithm has been developed to calculate tumor purity for giving an overall view of TME.¹⁶ Using this algorithm, stromal and immune scores of solid tumors can be generated based on single sample gene set enrichment analysis (GSEA). Several studies have shown that immune and stromal scores can predict the prognosis of a variety of cancers, including lung adenocarcinoma,¹⁷ clear cell renal cell carcinoma,¹⁸ and bladder cancer.¹⁹ However, Zhao et al²⁰ found that the immune and stromal scores were not significantly associated with overall survival of PCa patients.

In this study, we downloaded the immune and stromal scores of PCa samples from The Cancer Genome Atlas (TCGA) database from the online tool ESTIMATE and correlated these scores to clinicopathological characteristics and recurrence-free survival (RFS) of PCa patients. Based on these scores, the TME-related differentially expressed genes (DEGs) were identified and selected to perform functional enrichment analysis. Then Cox regression analyses were performed to identify TME-related prognostic genes and establish a predictive risk model. Moreover, we also evaluated the relationship between risk score and immune pathway.

Materials and Methods

Data Collection and Analysis

TCGA level 3 gene expression data of PCa samples were downloaded from the Genomic Data Commons (https://portal.gdc.cancer.gov/). The clinical data was extracted from the UCSC Xena project (https://xena.ucsc.edu). The inclusion criteria were as follows: (i) prostate adenocarcinoma diagnosed by pathologic histology, (ii) available survival information, and (iii) available expression data. The immune and stromal scores of PCa samples were downloaded from the online tool ESTIMATE (https://bioinformatics.mdanderson.org/estimate/). Subsequently, the best cutoff values were generated using X-tile plots²¹ and used to divide the PCa samples into high and low score groups. The Log rank test was used to identify the relationship between immune/stromal scores and RFS of patients. A p-value of less than 0.05 was considered to indicate statistical significance. The survival curves were visualized using GraphPad Prism 8 (GraphPad Software Inc., United States). Estimate the Proportion of Immune and Cancer cells (EPIC) database was employed to assess the difference in immune cell infiltration between high and low score groups.

Identification of Differentially Expressed Genes

Differential expression analysis was performed using the R v 3.6.1 software with “limma” package. Genes with a fold-change greater than 2 and a false discovery rate (FDR) below 0.05 were considered significantly differentially expressed. Volcano plots of DEGs were generated using the R software with “ggplot2” package. Furthermore, a list of immune-related genes was downloaded from the Immunology Database and Analysis Portal (ImmPort, https://www.immport.org/). The online Venn diagram tool (https://bioinfogp.cnb.csic.es/tools/venny/) was used to identify the intersection of genes among immune DEGs, stromal DEGs, and immune-related genes. These intersecting genes were selected for further analysis.

Enrichment Analysis of Intersecting Genes

Enrichment analyses of Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) for intersecting genes were performed using the R software with “clusterProfiler” package. A p-value <0.05 was considered statistically significant.

Establishment of a Predictive Risk Model

Further, we established a predictive risk model and identified the TME-related gene signature for predicting the prognosis of PCa patients. Firstly, using the “survival” R package, we evaluated the associations between the expression of selected genes and RFS of patients by univariate Cox regression analysis. The significant genes with P < 0.05 were identified for further analysis. Least Absolute Shrinkage and Selector Operation (LASSO) analysis can not only reduce the estimation variance but also provide an explicable final model.²² Then, we performed the LASSO analysis to identify the independent prognostic genes. Finally, we established a predictive risk model using multivariate Cox regression analysis. The risk score of each patient was generated using the following formula:

In the formula above, “Exp” represents the expression level of gene and “β” is the regression coefficient resulted from the multivariate Cox regression analysis.

Assessment of the Risk Model and Survival Analysis

Next, a time-dependent receiver operating characteristic (ROC) curve with the area under the curve (AUC) values was visualized to evaluate the sensitivity and specificity of the model using the “timeROC” package. Kaplan–Meier survival analyses were performed to explore the associations of risk score and independent prognostic genes with RFS of patients. Survival curves were compared with the use of the Log rank test with a two-sided significance level of 5%.

Relationship Between Risk Score and Immune Pathway

GSEA was performed to evaluate the relationship between risk score and immune pathway using the GSEA software (version 4.0.3).²³ The number of permutation was set at 1000. The significant level of the gene sets was set at the absolute normalized enrichment score (NES) >1 and nominal p-value <0.05.

Results

Immune and Stromal Scores Correlate with Clinicopathological Characteristics and Prognosis of PCa Patients

According to the included criteria, 485 patients with PCa were selected in this study. The immune and stromal scores of these corresponding samples were downloaded from the online tool ESTIMATE. Thereinto, immune score represents the infiltration of immune cells in tumor tissue, and stromal score captures the presence of stroma in tumor tissue. They are calculated using the ESTIMATE algorithm based on the expression matrix of each sample. Statistical analyses were performed with the use of GraphPad Prism 8 and revealed that high immune scores correlated with later T stages (P = 0.007) and lymphatic metastasis (P = 0.043), while there was no significant relationship between immune scores and age (P = 0.055) and Gleason scores (P = 0.176) (Figure 1A–D). However, high stromal scores correlated with an upper age (P = 0.023), later T stages (P = 4.72e-05), and higher Gleason scores (P = 0.002), while there was no significant relationship between stromal scores and M stage (P = 0.567) (Figure 1E–H).

Figure 1 Associations of immune and stromal scores with prostate cancer clinicopathological characteristics. (A–D) Distributions of immune scores among different ages, tumor stages, lymphatic metastases, and Gleason scores. (E–H) Distributions of stromal scores among different ages, tumor stages, lymphatic metastases, and Gleason scores.

To explore whether the immune and stromal scores can predict the prognosis of PCa patients, we evaluated the relationship of these scores and RFS. Based on the best cutoff value, PCa patients were divided into two groups: high score group and low score group. Specifically, 169 patients with higher immune scores and 316 patients with lower immune scores were separated by the best cutoff value of −578.3. Similarly, 298 patients with higher stromal scores and 187 patients with lower stromal scores were separated by the best cutoff value of −820.1. Kaplan–Meier survival analyses revealed that the high immune and stromal scores were strongly associated with poor RFS of patients (P = 0.003) (Figure 2A and B).

Figure 2 Kaplan–Meier survival curves of PCa patients and relative proportions of tumor-infiltrating cells. (A and B) Kaplan–Meier curves of RFS for patients with low vs high immune and stromal scores. (C and D) Relative proportions of tumor-infiltrating cells in the low vs high immune and stromal score groups. *P < 0.05, ***P < 0.001, ****P < 0.0001. Abbreviations: PCa, prostate cancer; RFS, recurrence-free survival; ns, no significance.

To further elucidate the reasons for the poor RFS of patients with higher scores, we evaluated the proportions of infiltrating immune cells in the TME of each sample using the EPIC database.²⁴ Consequently, the proportions of B cells, CAFs, CD4⁺ T cells, endothelial cells, and macrophages were dramatically higher in the high immune score group than that in the low immune score group (P < 0.001), while the proportions of CD8⁺ T cells were dramatically lower (P < 0.001) (Figure 2C). Likewise, the proportions of B cells, CAFs, CD4⁺ T cells, endothelial cells, and macrophages were dramatically higher in the high stromal score group than those in the low stromal score group (P < 0.05) (Figure 2D). These results indicated that the increase of suppressive immune cells and the decrease of killer immune cells might contribute to the poor RFS of patients with PCa.

Identification of Differentially Expressed Genes

To identify the TME-related DEGs, we downloaded level 3 RNA sequencing data of 485 PCa samples from the TCGA database. Using the R software with “limma” package, we identified 994 immune DEGs and 1470 stromal DEGs, respectively. The volcano plots in Figure 3A showed that 878 genes were upregulated and 116 genes were downregulated in the high immune score group compared with the low score groups (fold change >2, FDR < 0.05). In the comparison based on stromal scores, 1361 genes were upregulated and 109 genes were downregulated in the high score group (fold change >2, FDR < 0.05) (Figure 3B). Subsequently, we identified 238 intersecting genes among immune DEGs, stromal DEGs, and immune-related genes obtained from the ImmPort database (Figure 3C). The 238 intersecting genes showed in the Venn diagrams were selected for further analysis.

Figure 3 Identification of differentially expressed genes. (A) Volcano plot of DEGs based on immune score in PCa samples. (B) Volcano plot of DEGs based on stromal score in PCa samples. (C) Venn diagrams showing the overlapping genes among immune DEGs, stromal DEGs, and immune-related genes. Abbreviations: PCa, prostate cancer; DEGs, differentially expressed genes.

Functional Enrichment Analyses of Intersecting Genes

Next, we performed enrichment analyses to excavate the functional properties of the 238 differential TME-related genes. GO and KEGG terms were determined using the “clusterProfiler” package (P < 0.05). The results showed the top 10 GO terms (Figure 4A) and the top 10 KEGG pathway terms (Figure 4B). Go enrichment analysis showed that these TME-related genes were mainly involved in immune responses, especially immune cell activation and proliferation. Moreover, KEGG analysis revealed that these TME-related genes were obviously involved in immune-related pathways, such as cytokine–cytokine receptor interaction and chemokine signaling pathways.

Figure 4 Enrichment analysis of microenvironment-related differentially expressed genes. (A) The top 10 of biological processes GO terms. (B) The top 10 enriched KEGG pathways. Abbreviations: GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes.

Identification of a Prognostic Gene Signature and Survival Analysis

To further identify a prognostic gene signature associated with TME, we firstly performed the univariate Cox regression analysis. The result revealed that 94 differential TME-related genes were significantly associated with RFS of PCa patients (P < 0.05) (Supplementary Table 1). The selected 94 significant genes were subsequently utilized for the multiple LASSO regression analysis (Figure 5A). As shown in Figure 5B, we obtained two key genes, namely, C-type lectin domain containing 7A (CLEC7A) and collagen type XI alpha 1 chain (COL11A1). Table 1 listed the representative GO and KEGG categories in which CLEC7A and COL11A1 were involved. Moreover, Kaplan–Meier analysis revealed that the high expression levels of CLEC7A (P = 2.64e-05) and COL11A1 (P = 1.53e-05) were negatively correlated with RFS (Figure 5C and D), indicating that PCa patients with strong CLEC7A or COL11A1 expression have a higher risk for disease recurrence.

Table 1 Representative GO and KEGG Categories in Which CLEC7A and COL11A1 Were Involved

Figure 5 Identification of a prognostic gene signature and survival analysis. (A) LASSO coefficient profiles of the 45 significant genes from the result of the univariate regression. (B) Feature selection for prognostic biomarkers using the LASSO method. (C) Kaplan–Meier curves of RFS for patients grouped by expression levels of CLEC7A. (D) Kaplan–Meier curves of RFS for patients grouped by expression levels of COL11A1. Abbreviations: LASSO, Least Absolute Shrinkage and Selector Operation; RFS, recurrence-free survival; CLEC7A, C-type lectin domain containing 7A; COL11A1, collagen type XI alpha 1 chain.

Establishment and Assessment of the Predictive Risk Model

The predictive risk model was built with the two key genes. Based on the regression coefficients, the risk scores were generated according to the following formula: (0.415 * CLEC7A expression level) + (0.173 * COL11A1 expression level) (Supplementary Table 2). The score of 1.3 was used as the best cutoff to divide patients into high-risk group (n = 171) and low-risk group (n = 314) (Figure 6A). Compared with the low-risk group, the number of PCa recurrence was higher in the high-risk group (Figure 6B). Moreover, a high-risk score was significantly associated with lower RFS (P = 5.47e-09) (Figure 6C), indicating that the risk score could serve as a good predictive tool. The ROC curve analysis showed that the AUC of the predictive risk model for RFS was 0.694 and 0.731 for the third and fifth years, respectively, indicating that the risk model had medium performance (Figure 6D).

Figure 6 Establishment and assessment of the predictive risk model. (A) Distributions of risk score. (B) The RFS time of patients in high- and low-risk groups. (C) Kaplan–Meier curves of RFS for patients with low vs high risk scores. (D) ROC curves of the risk model for predicting 3- and 5-year survival rates. Abbreviations: RFS, recurrence-free survival; R, recurrence; RF, recurrence-free; ROC, receiver operating characteristic.

The Risk Score Predicts the Involvement of Immune Pathways

Further, we performed GSEA to elucidate the reasons for the poor RFS of patients with higher risk scores. The results showed that the pathways significantly enriched in the high-risk group involved immune response, immune system process, negative regulation of immune effector process, and negative regulation of immune system process (P < 0.05) (Figure 7A–D). These results confirmed the close relationship of the risk signature with the immune microenvironment.

Figure 7 GSEA of the risk score in PCa. The horizontal axis represents genes of the immune response (A), immune system process (B), negative regulation of immune effector process (C), and negative regulation of immune system process (D) gene sets, ranked by decreasing risk score. The vertical axis represents enrichment score. Abbreviations: GSEA, gene set enrichment analysis; PCa, prostate cancer; ES, enrichment score; NES, normalized enrichment score.

Discussion

TME plays an important role in tumorigenesis, development, and response to therapies, especially immunotherapy. Based on genetic abnormalities, tumor cells can interact with their microenvironment to secrete high levels of immunosuppressive cytokines, which inhibit the proliferation and function of effector T cell.²⁵ The immune microenvironment and TME-related genes have reportedly been involved in tumor progression and prognosis and are being pursued as potential therapeutic targets.^26–28 However, research on the impact of tumor immune microenvironment on the prognosis of PCa is still scarce. Therefore, it is necessary to study further.

In this study, we identified a novel TME-related gene signature for predicting the prognosis of PCa via a systematic bioinformatics analysis. Firstly, we obtained the immune and stromal scores of PCa samples from the ESTIMATE tool and analyzed correlations between these scores and clinicopathological characteristics and RFS of these patients. Our results revealed that the stromal and immune scores were markedly associated with age, tumor stage, lymphatic metastasis, and Gleason score. Both immune and stromal scores of later tumor stage group were significantly higher compared to earlier stage group, which is consistent with the result from a previous study.²⁰ In addition, Kaplan–Meier analysis revealed that higher immune and stromal scores were strongly associated with poorer RFS of PCa patients. These findings indicated that the TME composition was a critical factor which affects the clinical outcomes of PCa patients.

Tumor-infiltrating cells have been proved to have a crucial function in cancers and be associated with the development and prognosis of PCa. Compared to the normal prostate tissue, higher B cell infiltration was found in the PCa regions, suggesting that B cells can facilitate the development of PCa and serve as a therapeutic target.²⁹ As an important part of TME, CAFs provide cancer cells with a suitable microenvironment to promote tumor growth and metastasis.³⁰ Moreover, the alteration of gene expression in CAFs can be regarded as prognostic factors.^30,31 Similarly, increased infiltrations of regulatory T cells and M2 macrophages in PCa tissue are also associated with a poor prognosis.^32,33 In this study, we analyzed the proportions of infiltrating cells in the TME of high and low score groups to further elucidate the reasons for the poor RFS of patients with higher scores. Our results showed that the relative abundances of B cells, CAFs, CD4⁺ T cells, macrophages, and endothelial cells were significantly higher in both high immune and stromal score groups than that in the low score groups (P < 0.05), whereas the fraction of CD8⁺ T cells was significantly lower in the high immune score group (P < 0.0001). These results indicated that infiltrating cells in TME played an important role in the recurrence of PCa.

Next, we identified TME-related DEGs and performed enrichment analysis. A total of 238 intersecting genes were obtained for further analysis. The GO enrichment analysis showed that these TME-related genes were mainly involved in immune responses, especially immune cell activation and proliferation. Previous research has shown that cytokines are a family of low molecular weight proteins involved in a variety of immunobiological processes, including host defense, immune response, and inflammation response. Cytokines elicit their effects by binding to receptors on target cells.³⁴ Cytokine and chemokines affect tumorigenesis by directly regulating the biological function of tumor cells, stromal cells, and immune cells, as well as inducing neovascularization in the cancer microenvironment.^35,36 In this study, KEGG analysis revealed that 238 TME-related genes were obviously involved in immune-related pathways, such as cytokine–cytokine receptor interaction and chemokine signaling pathways. These results indicated that the TME-related genes were involved in the proliferation and progression of PCa by affecting the immune microenvironment.

Given the preciseness and accessibility of LASSO regression, we performed the LASSO regression analysis to screen covariates. Based on univariate analysis and regression coefficient, we identified two prognostic TME-related genes, namely, CLEC7A and COL11A1. Further, survival analysis indicated that high levels of CLEC7A and COL11A1 expression were strongly associated with reductions in RFS of PCa patients. Functional enrichment analysis revealed that CLEC7A and COL11A1 dramatically participated in the immune-related and cytoskeleton-related pathways.

CLEC7A, also known as dectin-1, is one of the best characterized C-type lectin receptors involved in various pathophysiological processes, including infection, allergy, regulation of inflammation, cancer, and other diseases.³⁷ Dectin-1 exerts dual functions in cancer, acting both as a tumor promoter and as a tumor suppressor. For example, dectin-1 activation on macrophages by ligating the lectin galectin-9 facilitates pancreatic carcinoma progression and peritumoral immune tolerance.³⁸ In contrast, dectin-1-mediated recognition of sialic acid-specific glycan structure on tumor cells can effectively lead to activation of macrophage tumoricidal response.³⁹ Dectin-1 can also direct dendritic cells to prime antitumor Th9 cells.⁴⁰ However, because the role of dectin-1/CLEC7A in PCa has not been reported in previous studies, little is known about their contribution to the development of PCa. So, more in-depth researches are still needed to explore the role of dectin-1/CLEC7A in PCa.

COL11A1, a minor fibrillary subtype of the collagen family, plays a vital role in growth, metastasis, and chemotherapy resistance of tumors.^41,42 It has been demonstrated that over-expressed COL11A1 was associated with poor prognosis in various cancers, including esophageal squamous cell carcinoma,⁴³ pancreatic cancer,⁴⁴ ovarian cancer,⁴⁵ and breast cancer.⁴⁶ Nonetheless, like CLEC7A, the role of COL11A1 in PCa has not been reported, so the role of COL11A1 in prostate cancer is unknown. Our findings indicated that a high level of COL11A1 expression was strongly associated with a reduction in RFS of PCa patients, and COL11A1 might serve as a potential prognostic marker and therapeutic target for PCa.

Subsequently, we constructed a predictive risk model with the two key genes, which revealed a great capacity for predicting RFS. Patients with PCa were stratified to two distinct groups based on the best cutoff risk score. Compared with the low-risk group (n = 314), the number of PCa recurrence was higher in the high-risk group (n = 171). Specifically, recurrence occurred in 26.9% (46/171) of patients with high-risk scores, while only 8.6% (27/314) of patients experienced recurrence in the low-risk group. In addition, a high-risk score was significantly associated with lower RFS. Therefore, we suggested that PCa patients in the high-risk group should receive more attention and more frequent follow-up after therapy.

Finally, GSEA was performed and confirmed the close relationship between the risk scores and immune pathways. As it turns out, the pathways significantly enriched in the high-risk group involved immune response, immune system process, negative regulation of immune effector process, and negative regulation of immune system process, suggesting that immunosuppression exists in high-risk PCa patients, and is associated with poor outcome.

Conclusion

Based on immune and stromal scores that were significantly correlated with outcomes, a list of TME-related DEGs were identified and assessed via comprehensive bioinformatics analysis. Then two such genes (CLEC7A and COL11A1) were selected to build a predictive risk model using Cox regression analyses. Moreover, a high-risk score was significantly associated with lower RFS. GSEA suggested that immunosuppression existed in high-risk PCa patients, confirming the close relationship of the risk signature with the immune microenvironment. Therefore, the immune-related gene signature could help clinicians to assess the prognosis of PCa patients and select appropriate targets for immunotherapy.