{'fraction_non_alphanumeric': 0.031049125404184034, 'fraction_numerical': 0.01749191631929565, 'mean_word_length': 5.152924256951103, 'pattern_counts': {'":': 0, '<': 13, '<?xml version=': 0, '>': 3, 'https://': 5, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 7, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Comprehensive Analysis and Identification of an Immune-Related Gene Signature with Prognostic Value for Prostate Cancer

Background The tumor microenvironment (TME) has recently been proven to play a crucial role in the development and prognosis of tumors. However, the current knowledge on the potential of the TME in prostate cancer (PCa) remains scarce. Purpose This study aims to elucidate the value of TME-related genes for PCa prognosis by integrative bioinformatics analysis. Materials and Methods We downloaded the immune and stromal scores of PCa samples via the ESTIMATE and correlated these scores to clinicopathological characteristics and recurrence-free survival (RFS) of patients. Based on these scores, the TME-related differentially expressed genes were identified for functional enrichment analysis. Cox regression analyses were performed to identify prognostic genes and establish a predictive risk model. Moreover, gene set enrichment analysis (GSEA) was performed to evaluate the relationship between risk score and immune pathway. Results The stromal and immune scores were associated with clinicopathological characteristics and RFS in PCa patients. In total, 238 intersecting differentially expressed genes were identified. Functional enrichment analysis further revealed that these genes dramatically participated in the immune-related pathways. The immune-related risk model was built with C-type lectin domain containing 7A (CLEC7A) and collagen type XI alpha 1 chain (COL11A1) using Cox regression analyses. Kaplan–Meier survival analysis showed that the expression levels of CLEC7A and COL11A1 were significantly associated with the RFS. Further, the RFS time in high-risk group was significantly shorter than that in low-risk group. The areas under the curve for the risk model in predicting 3- and 5-year RFS rates were 0.694 and 0.731, respectively. GSEA suggested that immunosuppression existed in high-risk PCa patients. Conclusion CLEC7A and COL11A1 were selected to build a predictive risk model, which may help clinicians to assess the prognosis of PCa patients and select appropriate targets for immunotherapy.


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. 1 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][3][4] Current treatments for PCa include conservative treatment, surgery, endocrine therapy, radiotherapy, and chemotherapy. 5 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. 6 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. 7 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. 8 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. 11 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. 12 For this reason, precision strategies have been developed and are in progress to target CAFs. [13][14][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. 16 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, 17 clear cell renal cell carcinoma, 18 and bladder cancer. 19 However, Zhao et al 20 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.

Data Collection and Analysis
TCGA level 3 gene expression data of PCa samples were downloaded from the Genomic Data Commons (https://por tal.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 21 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. 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 immunerelated 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. 22 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). 23 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.

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).
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 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).
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. 24 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.

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.

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 TMErelated 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, DovePress 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.

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 lowrisk 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).

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 highrisk 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.

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. 25 The immune microenvironment and TMErelated genes have reportedly been involved in tumor progression and prognosis and are being pursued as potential therapeutic targets. [26][27][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. 20 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 DovePress 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. 29 As an important part of TME, CAFs provide cancer cells with a suitable microenvironment to promote tumor growth and metastasis. 30 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. 34 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 TMErelated 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 cytoskeletonrelated 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. 37 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. 38 In contrast, dectin-1-mediated recognition of sialic acid-specific glycan structure on tumor cells can effectively lead to activation of macrophage tumoricidal response. 39 Dectin-1 can also direct dendritic cells to prime antitumor Th9 cells. 40 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, 43 pancreatic cancer, 44 ovarian cancer, 45 and breast cancer. 46 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 highrisk 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 highrisk score was significantly associated with lower RFS. GSEA suggested that immunosuppression existed in highrisk 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.

Author Contributions
All authors made a significant contribution to the work reported, whether that is in the conception, study design, execution, acquisition of data, analysis and interpretation, or in all these areas: took part in drafting, revising or critically reviewing the article; gave final approval of the version to be published; have agreed on the journal to which the article has been submitted; and agreed to be accountable for all aspects of the work.

Funding
There is no funding to report.

Disclosure
The authors declare that they have no competing interests.