Stratification of risk based on immune signatures and prediction of the efficacy of immune checkpoint inhibitors in prostate cancer

Prostate adenocarcinoma (PRAD) is a major threat to male health worldwide with a high mortality rate. New therapeutic strategies for the treatment of this malignant disease are of tremendous significance. Much attention has been paid to the involvement of immune cells in the prevention and treatment of cancer as well as how their regulatory systems contribute to effective cancer treatment. In this study, we constructed prognostic immune profiles based on The Cancer Genome Atlas (TCGA)-PRAD data sets and tested their predictive power on total and internal data sets. Then, we looked at how the lymphocyte of tumor invasion varied between the high-risk group and the low-risk group. Five immune-related genes made up the immune marker, which was an independent predictive factor in patients with PRAD. Patients in the low-risk score group had a higher rate of overall survival and a stronger infiltration of immune cells in the tumor microenvironment, which was highly related to clinical outcomes but required prospective validation.


Introduction
Prostate adenocarcinoma (PRAD) is a leading contributor to deaths in men globally [1].In order to effectively treat this cancer, new therapeutic approaches are crucial [2].Much focus has been paid to the involvement of immune cells in cancer treatment and prevention, as well as how their regulatory mechanisms provide effective cancer therapy.At present, Immunotherapy is a hot spot in clinical tumor therapy.A large number of clinical studies or basic experiments have proved that immune checkpoint inhibitors (ICIs) can play a wonderful role in tumor therapy [2].Nevertheless, only a small proportion of cancer patients can benefit from ICIs, and the overall response rate is low [1].In colorectal, ovarian and hepatocellular carcinoma, several immune-related genes (IRGs) have been identified as prognostic biomarkers.These genes have been linked to immune cells and relevant molecular mechanisms in the tumor microenvironment, and to some extent, they reflect the status of the immune response in the tumor microenvironment [3].More recently, using transcriptome sequencing data, computer algorithms like Tumor Immune Estimation Resource (TIMER) [4] and CIBERSORT have made it possible to determine the immune spectrum of cancer.The effect of immunotherapy is revealed by the immune landscape of cancer, which is also highly linked to patient prognoses.It is important to develop a comprehensive immune spectrum model to evaluate the immune-related efficacy of therapies.Few studies have conducted immune profiling and investigated IRGs for their predictive significance in PRAD.In this work, we developed an immune marker based on IRGs, and we looked at how it related to clinic-pathological characteristics and clinical outcomes in PRAD patients [5].In addition, we investigated the molecular functions, cellular components and possible pathways of Different expressed immune-related genes (DEIRGs) enrichment in the prediction model by molecular enrichment analysis.

Differentially expressed IRGs
All original data were obtained from the TCGA database.We examined the differentially expressed genes between PRAD and corresponding normal tissues to screen out IRGs involved in oncogenesis.The abnormally expressed genes were iden-tified using the software edgeR (version 4.0.0,John Chambers and colleagues, Auckland, New Zealand) (https://www.rproject.org/),as described previously (adjusted p value < 0.05 and |log2 (fold change)| > 1).The intersection between differentially expressed genes and the IRGs was termed differentially expressed IRGs.

Pathway and functional enrichment analyses
We carried out Pathway and functional enrichment analyses to investigate the biological significance of these differentially expressed IRGs.Kyoto Encyclopedia of Genes and Genomes (KEGG) [6] and Gene Ontology (GO) [7] enrichment analyses were carried out.Terms and pathways were regarded as significantly enriched objects if their false discovery rate was less than 0.05.The package ggplot2 was used to carry out the visualization.Pathway illustrations were generated using the "pathview" R package.

Constructing an IRGs-related immune signature for PRAD
We identified immune genes associated with prognosis.To verify their robustness and predictive ability, we developed a prognostic risk model.We first screened for immune-related genes with prognostic significance by univariate Cox proportional hazards regression analysis.Genes with a p-value of less than 0.05 were examined by the glmnet software using the minimum absolute contraction and selection operator (Lasso) to prevent over-fitting [8].After screening by the Lasso model, a multivariate Cox proportional hazards model was employed to develop the immune-related risk model as follows: risk score = Gene A level × coefficient a + gene B level × coefficient B + gene C level × coefficient C + ... + gene N level × coefficient N. The prognosis for PRAD patients is represented by the risk score in the model, where the lower the risk score, the better the prognosis.The patients were classified into high-risk and low-risk groups using the median risk score as the cut-off value.With log-rank p-value < 0.05 considered statistically significant (for packet survival and investigator), predictive ability was calculated using Kaplan-Meier survival curves.Using a time-dependent receiver operating characteristic curve (ROC), the prediction ability of this immune signature was evaluated [2,9,10].

Calculating the ratio of tumor-infiltrating immune cells
By using a deconvolution technique, CIBERSORT could determine the ratios of infiltrated immune cells from tissue transcriptional profiles.We computed the ratios of 22 types of tumor-infiltrating immune cells using the TCGA-PRAD transcriptional profiles and the CIBERSORT R script [11,12].

Statistical analysis
Chi-square and Student's t-tests were conducted to examine differences among variables.The effects of multiple clinicpathological characteristics along with the immune signature on patients' survival were evaluated using univariate and mul-tivariate cox regression analyses [13].Using R software programming language (version 0.0.0).Statistical analysis was carried out.The package "pheatmap" was used to create the heat maps.Statistical significance was defined as p-value < 0.05.

Quantitative reverse transcription-polymerase chain reaction (qRT-PCR)
The specific experimental steps are as described in the previous paper [14].In brief, total RNA was extracted from normal tissue and tumor tissue from PRAD patients using TRIzol (15596026, Thermo Fisher Scientific, Waltham, MA, USA).Quantitative reverse transcription-polymerase chain reaction (qRT-PCR) was conducted on the obtained RNA from each sample (2 µg) with FastStart Universal SYBR® Green Master (4913914001, Roche, Basel, Switzerland., USA) on a Q5 PCR System (Roche, USA).The primers were showed in Table 1.

Identification of prostate cancer prognosis-related differentially expressed genes
We first classified the prostate cancer samples in the TCGA database as "Alive" and "Dead" according to patients' overall survival profiles.After adjusting the cut-off value adjusted p value < 0.05 and |log2 (fold change)| > 1).Then, 1309 Differentially expressed genes were screened.At the same time, 2489 immune-related genes with an "Immune" correlation score greater than 2 were screened from the GeneCards database.After crossing, 18 immunologically related prognostic differentially expressed genes were obtained.

The prognosis of PRAD patients can be predicted by the five-gene immune signature
To evaluate the predictive significance of the chosen 18 differentially expressed IRGs and to avoid over-fitting, we further performed a LASSO analysis in which five of the 18 genes were predictors of patient outcome (Fig. 1A,B).We ranked patients' risk scores and outlined their distribution (Fig. 1C Top).A dot plot (Fig. 1C Middle) was used to display the risk scores, survival status and overall survival of the patients.The heat map displays the five IRGs' expression patterns in patients with varying risk scores.In the lowrisk group, Insulin Like Growth Factor 1 (IGF1) and Opioid Receptor Kappa 1 (OPRK1) were expressed at high levels, whereas in the high-risk group, Anti-Mullerian Hormone (AMH), Defensin Alpha 5 (DEFA5) and Serum Amyloid A1 (SAA1) were expressed at high levels.The predictive ability of this immune signature was further verified.Prostate cancer patients were classified into high-risk and low-risk groups (categories).Those with low-risk scores had longer Overall survival (OS) duration, according to the survival analysis (p = 0.004) (Fig. 1D).Moreover, this immune signature's area under the curve (AUC) value was 0.875 (Confidence interval (CI): 0.796-0.953)(Fig. 1E).Time-dependent ROC curves at 1-year, 3-year and 5-year durations (AUC = 1.000, 0.891 and 0.912, respectively) are also displayed (Fig. 1F).

The tumor immune microenvironment and the 5-IRG signature
To generate histograms displaying the proportion of 22 immune cells infiltrating in tissue samples from the PRAD patients, we computed the proportion of 22 immune cell types in each PRAD sample using the CIBERSORT method (Fig. 2A,B).The proportion of immune cells that differed between the groups with high and low-risk scores was compared.The percentage of B cells, Eosinophils, Macrophages, mast cells, Neutrophils, NK (Natural killing cell) CD56 (Neural Cell Adhesion Molecule 1) bright cells, NK cells Tgd, PDC (plasmacytoid dendritic cell), Th1 cells and T helper cells differed remarkably (Fig. 2C) [15].We compared the expression differences of prognostic DEirgs in TCGA-PRAC, with only AMH being significantly overexpressed in the tumor (Fig. 3A) while having the same expression differences in paired samples of tumor tissue (Fig. 3B).The expression of the five genes was significantly different and were demonstrated using a between-group differential gene volcano map (Fig. 3C) and a sequence map (Fig. 3D) based on patient survival status.Subsequently, both univariate and multivariate Cox regression analyses were carried out to investigate the prognostic value of model DEIRGs (Fig. 3E,F & Table 2).

Relationship between survival prognosis and 5-gene expression difference in prostate cancer patients
We explored the association of differences in 5-gene expression and disease-specific survival, overall survival and progression-free interval in prostate cancer patients.Although high AMH expression was not significantly linked to low OS (overall survival) rate.Prostate cancer patients with low expression of SAA1 had poor Overall Survival (p = 0.041) (Fig. 4A), patients with elevated AMH expression levels had worse overall survival (p = 0.127) (Fig. 4B).And patients with OPRK1-low expression had poor disease-specific survival and overall survival (p = 0.033 and 0.022, respectively) (Fig. 4C,D); Progression-free interval and disease-specific survival (p = 0.004, 0.029, respectively) were obtained in patients with prostate cancer with low expression of IGF1 (Fig. 4E,F) [16].

The efficacy of prognostic identification was verified
To better understand how a model gene's expression affects clinical prognosis in prostate cancer, ROC curve analysis was used to verify its efficacy in differentiating different clinical outcomes.We first validated the diagnostic efficacy of a five-gene prognostic model for prostate cancer disease status as DEFA5 (AUC = 0.464), IGF1 (AUC = 0.636), SAA1 (AUC = 0.553), OPRK1 (AUC = 0.493), and AMH (AUC = 0.809), respectively (Fig. 5A).And then, we validated the efficacy of 5-gene prognostic models for the identification of Clinical variable outcomes of Alcohol history, Smoker, Lymphnode neck dissection, Radiation therapy, Age, Clinical stage, Histologic grade, Clinical T stage, pathologic stage, pathologic N stage, pathologic T stage, respectively (Fig. 5B-L).These results suggest that, for the assessment and anticipation of clinicopathological factors in prostate cancer, the 5gene prognostic model may be a potential biomarker.

The prognostic model is associated with prostate cancer
To further explore the association of 5-gene prognostic models with prostate cancer, we explored their association with the expression of well-known molecular markers known to be associated with prostate cancer.We selected some reliable prostate cancer biomarkers for further analysis to explore whether the 5-gene prognostic model is associated with the molecular mechanisms of prostate cancer initiation and progression.It has been reported that transmembrane serine endopeptidase 2 (TMPRSS2) and ETS-related genes (ERG) frequently fuse during Principal component analysis (PCA) [17].A mitochondrial and peroxisome enzyme, Alphamethylacyl-coa racemase (AMACR), participates in fatty acid oxidation and intermediate products of bile acids [18].The over-expression of AMACR on prostatic epithelial cells occurs in about 80% of prostate cancer cases, and the detection of AMACR by immunohistochemistry (IHC) helps to distinguish benign prostate cancer from prostate cancer [19].At the same time, studies have reported that AMACR overexpression in prostate cells predicts poor prognosis (associated with a high Gleason score and high initial PSA levels) and an increased risk of bone metastasis.In addition, there have been reports linking several Solute Carrier Family 45 Member (SLC45A) family genes to prostate cancer [20][21][22].As SLC45A3 (prostate cancer-associated protein 6/P501S/protein) is mainly expressed in the Gleason apparatus of prostate epithelial cells and is highly specific for prostate gland cells, it has been used to distinguish between other tumor types and metastatic prostate cancer [22,23].SLC45A3 may be expressed in PSA-negative prostate tumors, so the combined use of these markers may increase the sensitivity to identify prostate cancer metastases.SLC45A3 is weakly expressed in some aggressive tumors and is associated with increased Gleason score and risk of tumor recurrence.Overall findings revealed that in the 5-gene prognostic model, OPRK1 expression was positively linked to BRCA1 DNA Repair Associated (BRCA1), BRCA2 DNA Repair Associated (BRCA2), and AMACR (Fig. 6A,B,E).The positive correlation between OPRK1 expression and TMPRSS2 was significant, while the negative correlation between SAA1, AMH and TMPRSS2 was significant (Fig. 6C).SAA1 was positively correlated with Hexokinase 3 (HK3), SLC45A1, SLC45A3, SLC45A4 (Fig. 6D,F-H).

Functional enrichment analysis
In addition, we summarised the differential genes between the low-risk and high-risk categories, followed by GO and KEGG enrichment analyses.Based on the results of the GO enrichment analysis (Fig. 7A-D & Table 3), the BP with significance were negative regulation of the reproductive process, ovulation cycle, oogenesis, interleukin-1 production and regulation of Interleukin-1 production, respectively.Cell component (CC) was enriched in transport vesicles, exocytic vesicles, cytoplasmic vesicle lumen, vesicle lumen and secretory granule lumen.Hormone activity, growth factor activity, insulin receptor binding, insulin-like growth factor receptor binding,  and transforming growth factor beta receptor binding were all enriched in MF.Furthermore, KEGG analysis revealed some immune-related pathways.
Subsequently, the 5 genes of the prognostic model were mapped in red letters on the KEGG pathway map.Pathway enrichment analysis suggested that the 5-gene prognostic model was enriched in aldosterone regulation, steroid production, long-term depression, longevity regulatory pathways, and Staphylococcus aureus infection, and the well-known Epidermal Growth Factor Receptor (EGFR), Hypoxia Inducible Factor 1 (HIF-1), The AMP-activated protein kinase (AMPK), and other signaling pathways, suggesting that a variety of cancer-related (Fig. 7E,F & Table 4).Fig. 8 shows the role of 5 genes involved in Prostate cancer, Long-term depression, Endocrine resistance, Transforming Growth Factor Beta (TGFβ) signaling pathway, and EGFR tyrosine kinase inhibitor resistance.

Gene set enrichment analysis based on a 5-gene prognostic model to predict differences between low and high-risk groups for prostate cancer
The potential mechanisms by which prognostic models regulate OS timing in PRAD patients need to be further explored.To find the possible signaling pathways involved in the prognostic model in the two risk groups, we performed GSEA (Fig. 9, Supplementary Table 1).The KEGG high-risk group was primarily enriched in Ribosome, Cell Cycle, Spliceosome, Ascorbate and Aldarate Metabolism, and Base Excision Repair, whereas the low-risk group was primarily enriched in Cardiac Muscle Contraction, Hypertrophic Cardiomyopathy HCM, Dilated Cardiomyopathy, Proximal Tubule Bicarbonate Reclamation, and Primary Immunodeficiency (Fig. 9A,B).In WikiPathways, high-risk groups were enriched in Striated Muscle Contraction Pathway, and Extrafollicular B Cell Acti-  vation by Sarscov2, Folate Metabolism, Selenium Micronutrient Network, and Complement System (Fig. 9C,D).The highrisk groups in the Reactome database were mostly enriched for striated muscle contraction, CD22-mediated B-cell receptor (BCR) regulation, clearance of heme in plasma, and creation of initial triggers for the complement system, such as C4 and C2 activators; The low-risk group was enriched in Striated Muscle Contraction, CD22 Mediated BCR Regulation, Scavenging of Heme from Plasma, spindle checkpoint and sister chromatid, of C 4 and C 2 creation Activators, Initial Triggering of Complement (Fig. 9E,F).In Persistent Identifier (PID), the high-risk group was enriched in Polo Like Kinase 1 (PLK1), ATR Serine/Threonine Kinase (ATR), Aurora B, Fanconi and Aurora a Pathway, while the low-risk group was enriched in NCADHERIN, INTEGRIN3, INTEGRIN2, Interleukin (IL42) pathways (Fig. 9G,H).Finally, the expression of hub gene was confirmed by PCR (Fig. 10).

Discussion
Overall, we analyzed DEIRGs associated with PRAD prognosis by LASSO regression, with predicted prognostic models including DEFA5, IGF1, SAA1, OPRK1 and AMH.During cancer development, t-cell infiltration is impaired, antigen regulation is disrupted, and expression levels of immune checkpoints and their ligands are elevated [24].Hence, patients may not benefit greatly from ICI treatment when immune checkpoints are not the only rate-limiting step [25].Patients in the current study who had low risk scores had a higher level of immune checkpoint molecules.Indirectly, the rise in immune checkpoint levels, including PD-1 and Cytotoxic T-Lymphocyte Associated Protein 4 (CTLA-4), in the low-risk score group suggests pre-existing t-cell activation.Patients with low-risk scores may therefore be more sensitive to ICI treatment [26].
Although we elucidated the role of the 5-gene prognostic model in predicting prognosis and verified the expression levels of DEFA5, IGF1, SAA1, OPRK1 and AMH in tumor  tissues by proteomics, there are still some limitations, to be further explored.To explore the characteristic function of prostate cancer, AMH and OPRK1 were found to be the key genes affecting Tumor microenvironment (TMB) and potential targets for immunotherapy.It was suggested that higher expression of AMH was beneficial to PCA progression, and patients' OS was worse when the OPRK1 level was lower, along with a higher recurrence rate in the early stage of PCA than those with high expression of OPRK1 [23,27].The onset and progression of prostate cancer may be influenced by decreased circulation levels of the testis-derived TGF-family peptide hormone, anti-mullerian hormone (AMH) [28].In a study of 1000 patients, the overall risk of prostate cancer has been demonstrated to be correlated with prediagnostic serum AMH levels [29].Physiological concentrations of endogenous AMH increased the viability of cancer cells [30].Furthermore, serum amyloid A1 (SAA1) is an inflammatory high-density lipoprotein.It is also regarded as a prognostic indicator and a predictor of cancer risk [31].It is crucial for fatty acid oxidation (FAO), Association of Tennis Professiona (ATP) and cell growth in prostate cancer tissues, which are induced by the loss of the cancer biomarkers SUN2 [32,33].The inhibitory effect of the IGF1/Signal Transducer and Activator of Transcription 3 (STAT3) pathway on the invasion of prostate cancer is well known [34].Correlations between IGF1 and Wnt/β-catenin signaling have also been found in human prostate cancer samples [35].Single-cell transcriptomic analysis has shown that AR activation enhances the mecasermin signaling pathway (IGF1) and initiates oncogenic transformation [36].Elevated IGF1 signaling further accumulates the Wnt/β-catenin pathway in transformed cells to promote prostate tumor development [37].The detection rate of DEFA5 in cancer was significantly higher in patients.In addition to the role of defensins in the innate immune system, which is associated with antimicrobial and immune signaling activities, the significance of α-defensin 5 (DEFA5) in the onset and progression of gastric cancer is being discussed in an increasing number of studies [3,38].A study showed mechanistically how DEFA5 binds to Body Mass Index (BMI) directly, then decreases its binding at the Cyclin Dependent Kinase Inhibitor 2A (CDKN2A) locus, and upregulates the expression of the cyclin-dependent kinase inhibitors p16 and p19 to inhibit stomach cancer [1,2].
To further improve the reliability and persuasiveness of guiding clinical strategic decisions, we incorporated expression profiles and survival data from TCGA prostate cancer during the analysis.In addition, it was shown that several specific DEIRGs with remarkable variations in variable risk factors were linked to the advancement of prostate cancer and patient prognosis.More importantly, these DEIRGs, as molecular biomarkers, have an important role in predicting and evaluating survival outcomes in prostate cancer.
Despite some positive results, there are still some problems.First, this immune signature is built on a common data set.Predictive power needs to be further tested in Randomized controlled trials.First, we did not combine metabolomics  with immunogenomics testing.In addition, the practical value of our prognostic gene model for prostate cancer needs to be further validated.To sum up, we comprehensively analyzed and validated 5-gene prognostic models to predict the clinical prognosis of prostate cancer.The findings will help to establish a reliable and referential risk assessment model and provide new insights into immune-related research and treatment strategies.

Conclusions
In this study, we constructed prognostic immunoprofiles based on the TCGA-PRAD dataset and tested their predictive power on both total and internal datasets.There was a significant difference in tumor infiltrating lymphocytes between high-risk and low-risk groups, and an immune marker composed of five immune-related genes was an independent predictor of PRAD patients.Patients in the low-risk category had better overall survival and greater infiltration of immune cells in the tumor microenvironment, which is highly correlated with clinical outcomes but requires prospective validation.

AVA ILABILITY OF DATA AND MATERIALS
The data are contained within this article and supplementary material.

A UTHOR CONTRIBUTIONS
YDS and XYY-designed the research study.YY and ZHLperformed the research.FMK, ZHW and SQR-analyzed the data.YDS, ZHW and ZCL-wrote the manuscript.All authors read and approved the final manuscript.

F
I G U R E 1. Lasso analysis and forest plot show multivariate Cox model results for 5 immune-related genes.(A) Cross validation of optimal parameter selection in LASSO model.(B) LASSO coefficient profiles of 18 prognostic immune-related genes; (C) The risk factor graph shows the distribution of risk scores (Top), the relationship between risk scores and survival time (middle), and the expression patterns of five immune-related genes in the high-and low-risk groups (bottom); (D) Kaplan-Meier curves of overall survival of PRAD patients in high-and low-risk score groups; (E) Receiver performance curve based on risk score; (F) Time-dependent receiver operating characteristic curves.AUC: area under the curve; TPR: True positive rate; FPR: False positive rate; OS: Overall survival; AMH: Anti-Mullerian Hormone; DEFA: Defensin Alpha; IGF: Insulin Like Growth Factor; OPRK: Opioid Receptor Kappa; SAA: Serum Amyloid A. F I G U R E 2. Evaluation of immune infiltration.(A,B) Stacked bar graphs showing the proportion of infiltrating immune cells between the high-and low-risk groups.(C) Boxplot of the difference in immune cell infiltration between the high-and low-risk groups (Wilcoxon's test; *: p < 0.05; **: p < 0.01; ***: p < 0.001).Horizontal coordinates indicate the type of infiltrating immune cells, and the vertical axis indicates the proportion of infiltrating immune cells.Immune cell abbreviation: Th: Helper T lymphocytes; TR: regulatory T cell; NK: Natural killing cell; CD56: Neural Cell Adhesion Molecule 1; DC: Dendritic cell; PDC: plasmacytoid dendritic cell; TFH: Follicular helper T cells; Tgd: gamma delta T.

F I G U R E 3 .
Differential expression of prognostic DEIRGs in TCGA-PRAC.(A,B) The same difference was found in tumor tissues and normal samples.(C,D) Volcano plot and sequence plot showed the difference and significance of the five genes in the prognostic model based on the survival status of the patients.(E,F) Forest plot of 5 modifiers with prognostic value in univariate and multivariate Cox regression models.DEIRGs: Different expressed immune-related genes; TCGA: The Cancer Genome Atlas; PRAC: Prostate adenocarcinoma; AMH: Anti-Mullerian Hormone; DEFA5: Defensin Alpha 5; IGF1: Insulin Like Growth Factor 1; OPRK1: Opioid Receptor Kappa 1; SAA1: Serum Amyloid A1; GAPDH: Glyceraldehyde-3-Phosphate Dehydrogenase.