Abstract Prostate cancer (PCa) is a malignancy that affects men and is characterized by metastasis and high rates of morbidity. The objective of this study was to explore novel PCa biomarker with potential diagnostic and therapeutic value and relationships between it and tumor immunity and development. A total of 32 key genes were screened out via LASSO based upon 188 intersection genes obtained from WGCNA and DEGs analysis in [31]GSE32571, and PPP1R14B was further identified by COX regression based on the TCGA database and validated by qRT-PCR. Although it has been reported that PPP1R14B may have a certain correlation with the prognosis of uterine corpus endometrial carcinoma, breast cancer and gastrointestinal cancer, there are none of studies about correlation between PPP1R14B and PCa. Predictive ability analysis showed that PPP1R14B had greatly predictive values in occurrence and prognosis of PCa. Immune analysis revealed that overexpression of PPP1R14B was related to the increase of ALKBH2, UCK2, RAC3 and RAB17 and the decrease of CD40, DKK3, COL17A1 and PGRMC1, which would result in downregulation of plasma cells, upregulation of T regulatory cells and disorder of macrophage proportion to suppress adaptive immune directed against PCa. GSEA analysis showed that PPP1R14B, as an inhibitor of PP1, its overexpression was mainly involved in regulating pathways associated with MYC, E2F, PFN1 and so on, which was participated in the regulation of immune factors such as CD40, RAC3, COL17A, DKK3, as well as biological processes such as proliferation and migration. Patients with higher PPP1R14B expression responded more sensitively to drugs selumetinib and vorinostat, zebularine, azacitidine and VER155008. In summary, PPP1R14B was a potential diagnostic and prognostic biomarker of PCa and its high expression had closely association with tumor immune inhibition, proliferation and migration, providing a new target for drug therapy and immunotherapy in PCa. Keywords: prostate cancer, bioinformatics, biomarkers, tumor immunity, anti-PCa drugs Introduction The prostate is a part of the male reproductive system, located at the bottom of the bladder and surrounding the urethra. The prostate is one of the important male organs and prostate fluid provides nutrition and protection for sperm, helping to maintain their survival and motility. Prostate cancer is one of the most common cancers among men worldwide. Recent data showed that PCa accounts for 26% of new male malignant tumor cases, ranking first in incidence rate, and it ranks second in the mortality rate of male malignant tumors [32]^1^, [33]^2. Early-stage PCa can be controlled through treatments such as radiation therapy and surgery. However, as PCa progresses and invades, early-stage PCa can develop into more aggressive metastatic prostate cancer (mPC) with poor prognosis [34]^3. The main treatment method for metastatic hormone sensitive prostate cancer (mHSPC), which is one of the subtypes of mPC, is androgen deprivation therapy (ADT). However, relevant studies have shown that ADT has a risk of inducing cardiovascular disease [35]^4^, [36]^5. For another subtype of mPC, metastatic castration resistant prostate cancer (mCRPC), abiraterone and docetaxel are commonly used as their main treatment drugs [37]^6^, [38]^7. However, resistance to these drugs is inevitable during the treatment process, and research reports have shown that cross resistance to different drugs may occur during the treatment of mCRPC [39]^8^, [40]^9. Although there have been some advancements in the treatment of PCa, the current options for treating mPC are still inadequate, which made it become a significant medical issue for men. Therefore, it is very important to further study the molecular mechanism of PCa occurrence and provide new diagnostic and therapeutic biomarkers for PCa. In this study, based on gene expression data and clinical sample information obtained from Gene Expression Omnibus (GEO) and the Cancer Genome Atlas (TCGA) databases, we used weighted Gene Co-expression Network Analysis (WGCNA), least absolute shrinkage and selection operator (LASSO) logistic regression and COX proportional hazards regression analysis to screen for biomarker which has strong correlation with PCa. qRT-PCR was employed to validate for the expression difference of biomarker from prostate cells with/without PCa. Subsequently, COX regression analysis, nomograms, Receiver Operating Characteristic Curve (ROC) and Decision Curve Analysis (DCA) were used to validate biomarker to obtain the correlation between biomarker and the occurrence, survival, and prognosis of PCa. In addition, immune infiltration analysis is also used to explore the relationship between biomarker and immune cells, immune checkpoints, and immune chemokines, aiming to gain a deeper understanding of the immune mechanisms involved in the development of PCa. Furthermore, GSEA analysis was performed on biomarker to understand their functions. What's more, we conducted drug sensitivity analysis to identify effective drugs that may have a positive impact on PCa. Finally, the competitive endogenous RNA (ceRNA) network was constructed to describe the regulatory mechanism of biomarkers in PCa. Materials and Methods Data acquisition and analysis The research process of this study was shown in Fig. [41]1. The Series Matrix File of [42]GSE32571 was downloaded from the GEO database ([43]https://www.ncbi.nlm.nih.gov/geo/) using the keyword "prostate cancer". [44]GSE32571 is a dataset based on gene signatures in PCa linked to the Gleason score, containing a total of 98 samples (59 tumor samples and 39 normal samples). Inclusion criteria for tumor samples considered: (i) all patients were male and diagnosed with pathologically confirmed PCa; (ii) Gleason score is 6-10; (iii) Primary tumor tissue undergoing radical prostatectomy. Exclusion criteria include (i) excluding samples from patients with a prior history of malignancy; (ii) Exclude samples from patients with multiple malignant tumors present simultaneously; (iii) Exclude samples from patients with metastatic tumors. The clinical information table of the samples can be found in [45]Supplementary Table 1. Figure 1. [46]Figure 1 [47]Open in a new tab Overview of the study design. After analyzing and processing the data based on a p-value less than 0.05 and an absolute value of log fold change (|logFC|) greater than 1, we screened differentially expressed genes (DEGs) from 98 samples. In addition, volcano plots and heatmaps were created to visualize the differential expression of DEGs using the aforementioned data. Furthermore, TCGA database ([48]https://portal.gdc.cancer.gov/) was accessed to obtain the clinical information of 361 PCa samples and 52 normal samples. The inclusion and exclusion criteria for these tumor samples are consistent with [49]GSE32571. The clinical information table of the samples can be found in [50]Supplementary Table 2. WGCNA analysis Genes in GSE 32571 were analyzed using WGCNA to explore critical module genes. The Pearson correlation coefficients between genes were calculated based on [51]GSE32571. The scale-free network was constructed and the appropriate threshold was selected for node network construction. Using two-step construction, after the adjacency matrix was transformed into a topological overlap matrix, the clustering tree was generated through hierarchical clustering. Furthermore, clustering was combined through a dynamic cut to obtain co-expression modules. The correlation between modules and traits was measured, and modules with a correlation greater than 0.6 were considered critical modules. Based on DEGs and critical module genes, a Venn diagram was established on bioinformatics to screen intersection genes. Enrichment analysis Intersection genes were used for Gene Ontology (GO) enrichment analyses (including cellular components (CCs), molecular functions (MFs), and biological pathways (BPs)) and Kyoto Encyclopedia of Genes and Genomes (KEGG). Statistical analysis was performed using the DAVID ([52]https://david.ncifcrf.gov/) and graphical charting was utilized by the clusterProfiler R package. We subsequently performed a functional enrichment analysis by the Metascape ([53]http://metascape.org). Acquisition of hub gene Machine learning is a breakthrough technology that is flexibly applied in the collection and processing of large amounts of data. It involves learning, identifying, and grasping associations within the data through the analysis of vast amounts of data. Utilizing these learned insights, machine learning enables prediction, classification, optimization, and decision-making [54]^10. In this study, LASSO was employed to further select key genes related to PCa among the intersection genes. The "glmnet" package in R was utilized for LASSO logistic regression analysis. The COX regression model is based on the concept of hazard ratio, which can estimate the risk ratio of the occurrence of an event. It is usually used for survival analysis and prediction of individual survival time. In this study, the "survival" package of R language was used to perform univariate COX regression analysis on selected data from the TCGA database to screen out the hub gene. Validation of hub gene by qRT-PCR qRT-PCR was performed to confirm the expression of hub gene in the PCa and normal control (NC) groups. Human PCa epithelial cell lines (PC-3) and normal prostate epithelial cell lines (RWPE-1) were purchased from the Institute of Biochemistry and Cell Biology of the Chinese Academy of Sciences (Shanghai, China) ([55]Supplementary Fig. 1). PC-3 cells are originated from a 62-year-old white male with grade IV primary prostate cancer (Gleason score > 6). Media for cell culture and fetal bovine serum (FBS) were obtained from GIBCO. The PC-3 and RWPE-1 cell lines were cultured in Ham's F12 nutrient medium (Invitrogen). RPMI 1640 medium (Invitrogen) were supplemented with 10% FBS (Sigma) and keratinocyte serum-free medium (Invitrogen), respectively. All the cells were cultured at 37 ℃ in an atmosphere of 5% CO[2]. All dissected samples were immediately stored in liquid nitrogen until they were prepared for total RNA extraction. Briefly, total RNA was extracted from the PCa and NC groups using TRIZOL reagent (Life Technologies, Gaithersburg, MD, USA) according to the manufacturer's instructions. To detect the purity of the RNA, the OD values were measured using a NanoDrop (NanoDrop One, Thermo Fisher Scientific, Waltham, MA, USA). Reverse transcription was performed on each RNA sample (500 μg) using the PrimeScript RT reagent kit with a gDNA eraser (Takara Bio, Dalian, China) according to the manufacturer's instructions. For qRT-PCR analysis, the reaction mixture (including cDNA, DEPC water, forward primers, and reverse primers) was run in an iQ5 real-time PCR machine (Bio-Rad, Hercules, CA, USA) according to the instructions for SYBR Premix Ex Taq II (Takara Bio, Dalian, China). All the expression levels were normalized to those of the internal control (β-actin). The qRT-PCR cycles were as follows: step1, preparative denaturation (30 s at 95 ℃); step 2, 40 cycles of denaturation (5 s at 95 ℃) and annealing (30 s at 60 ℃); and step 3, dissociation, following the manufacturer's protocol. Forward primer of hub gene was as followed: GCTGGGTTCCTCCGAGGTA; reverse primer of hub gene was as followed: GAGTTGAAAGAGGGTGTGAGAG. The 2^-ΔΔct method was used to calculate relative gene expression levels. Correlation between hub gene and PCa To further validate significance of hub gene in PCa, correlations between hub gene and PCa were investigated with TCGA and GEO datasets. Differentially expression of hub gene was analyzed and visualized using the DESeq2 R package from TCGA database. Based on GSE 32571, the ROC curves of hub gene and known PCa-related genes GDF15 and AMACR and area under the curve (AUC) were plotted by the pROC R package. Nomogram was designed using the rms R package according to GEO datasets to construct the risk profile model of PCa occurrence based on hub gene. Calibration curve and DCA were used to evaluate the predictive accuracy and reliability of the nomogram. Further, according to COX regression analysis, the risk score of the hub gene is obtained, and its survival curve is then plotted using the "survminer" R package. Correlation between PCa and hub gene integrating clinical features In order to explore the correlation between PCa and hub gene binding clinical features, we conducted univariate and multivariate COX regression analyses on hub gene and 5 clinical features using the "survival" package using TCGA datasets, respectively. Further, the prognostic nomogram was constructed by integrating hub gene expression and clinical factors. Moreover, ROC curves, calibration curves and DCA in both training and validation cohorts were performed to assess the accuracy and discrimination of the nomogram. Immune infiltration To analysis immune infiltration of PCa and control groups, we used the GEO datasets to calculate the infiltration proportion of 22 kinds of immune cells with the CIBERSORT algorithm and analyzed differences with a significance threshold of p < 0.05. What's more, correlations between hub gene and immmune cells, checkpoints and chemokines were studied by the “dplyr” and “ggplot2” packages. GSEA analysis GSEA ([56]http://software.broadinstitute.org/gsea/index.jsp) is a powerful tool that can be used to reveal the biological significance of gene expression data, and it can also help researchers identify key biological functions and pathways under different conditions. Based on the hub gene expression data obtained from GEO, GSEA was used to conduct pathway enrichment analysis for hub gene high-expression and low-expression samples respectively. |NES|> 1, p < 0.05 and FDR < 0.25 were considered significant enrichment. Drug Sensitivity Analysis The Cancer Therapeutics Response Portal (CTRP) data and GEO datasets were utilized to calculate the sensitivity of common anticancer drugs targeting PCa and hub gene according to the half maximal inhibitory concentration (IC50) using the “oncoPredic” package. Construction of ceRNA network To better understand the mechanism of how lncRNA regulates gene expression through miRNA, the ceRNA network was constructed based on the miRWalk ([57]http://mirwalk.umm.uni-heidelberg.de), TargetScan ([58]http://www.targetscan.org/) and miRNet ([59]https://www.mirnet.ca/miRNet/home.xhtml) databases. We predicted the miRNAs of hub gene using miRWalk and TargetScan databases, and obtained intersection results from the two databases. Next, the miRNet database was used to predict potential lncRNAs of overlapping miRNAs. Subsequently, a ceRNA regulatory network was constructed by integrating the interaction between miRNA and lncRNA. Statistical analysis All statistical analyses were conducted by R language (version 4.3.1). The Kaplan Meier method was used for survival analysis. The prognostic potential of hub genes was studied using univariate and multivariate COX models. A significance level of P<0.05 is considered statistically significant. Results Attainment of intersection genes WGCNA was applied to analyze genes in [60]GSE32571. A sample clustering tree was obtained in Fig. [61]2a. Then we chose the soft-threshold 8 according to Fig. [62]2b. Finally, 38 modules were identified (Fig. [63]2c) and we ascertained four critical modules including the blue module (cor= 0.82, p = 1e^-24), the orange4 module (cor = 0.73, p = 3e^-17), the violet module (cor = 0.68, p = 5e^-14) and the salmon module (cor = 0.64, p = 4e^-12) (Fig. [64]2d). The orange4 module had a positive correlation with PCa while the other three modules showed negative relationship with this disease. The scatter plots of critical modules were shown in Fig. [65]2e-h. Figure 2. [66]Figure 2 [67]Open in a new tab WGCNA analysis results. a The samples clustering tree, b network topology analysis, c clustering tree, d module heatmap, e scatter plot of blue module, f scatter plot of orange4 module, g scatter plot of violet module, h scatter plot of salmon module. 222 DEGs (including 58 up-regulated genes and 164 down-regulated genes) were acquired from the dataset [68]GSE32571, which was obtained from the GEO database. The volcano plot (Fig. [69]3a) and heat map (Fig. [70]3b) were used to visualize the DEGs. The intersection of critical module genes and 222 DEGs yielded 188 intersection genes (Fig. [71]3c). Figure 3. [72]Figure 3 [73]Open in a new tab Acquisition of intersection genes. a Volcano map of DEGs, b heatmap of DEGs, c Venn diagram of critical module genes and DEGs. Enrichment analysis GO, KEGG and Metascape analyses were performed to investigate the functional mechanisms of 188 intersection genes. A total of 188 intersection genes were enriched in (i) CCs: extracellular space, extracellular exosome and focal adhesion and so on; (ii) MFs: extracellular matrix structural constituent, actin binding and protein homodimerization activity and so on; and (iii) BPs: glutathione derivative biosynthetic process and so on (Fig. [74]4a). The results of KEGG pathway enrichment exhibited that enriched pathways mainly contained focal adhesion, glutathione metabolism and drug metabolism - cytochrome P450 and so on (Fig. [75]4b). The results of Metascape enrichment analysis showed that these intersection genes were mainly enriched in functions and pathways such as supramolecular fiber organization, glutathione metabolism and protein homodimerization activity and so on (Fig. [76]4c-e). Figure 4. [77]Figure 4 [78]Open in a new tab Enrichment analysis results. a GO enrichment for three terms, b KEGG enrichment, c-e Metascape enrichment. Acquisition and validation of hub gene Then, based on the 188 intersection genes, a LASSO regression model was constructed with the optimal λ value of 0.00911853 to further screen key genes. We applied LASSO regression analysis to select 32 predictive genes as key genes (Fig. [79]5a-b). Figure 5. [80]Figure 5 [81]Open in a new tab Acquisition and validation of hub gene. a-b LASSO regression analysis results, c COX regression analysis results, d qRT-PCR results in cells with/without PCa. Furthermore, based on the clinical information obtained from the TCGA database, we conducted a univariate COX regression analysis on 32 key genes. The forest plot showed that the p-value of PPP1R14B was less than 0.05, which was statistically significant, and it had the highest HR value of 3.47 among the 32 genes (Fig. [82]5c). qRT-PCR showed that the expression of PPP1R14B (Fig. [83]5d) in the PCa group ([84]Supplementary Fig. 1a) was significantly higher than that in the NC group ([85]Supplementary Fig. 1b). Therefore, PPP1R14B is considered the hub gene of this study. Correlation between hub gene and PCa The boxplot showed that hub gene PPPP1R14B was more greatly expressed in tumor tissues than in normal tissues based on TCGA datasets (Fig. [86]6a). ROC analysis indicated that PPP1R14B and two known PCa-related genes GDF15 and AMACR had AUC values > 0.8, with PPP1R14B having the largest AUC value (AUC, 0.972), AMACR the middle (AUC, 0.922) and GDF15 the smallest (AUC, 0.889) (Fig. [87]6b). Nomogram prediction models, calibration curve and DCA showed that PPP1R1B was a high-risk factor in the occurrence of PCa (Fig. [88]6c-e). The survival curve of PPPP1R14B drawn according to the risk score showed that the survival probability of the low-risk group was higher when the survival time was the same (Fig. [89]6f). Figure 6. [90]Figure 6 [91]Open in a new tab Correlation between Hub Gene and PCa. a Boxplot of hub gene expression in different groups in TCGA datasets, b ROC analysis results of hub gene, GDF15 and AMACR, c-e nomogram, calibration curve and DCA results of hub gene for PCa occurrence, f survival curve of hub gene. Correlation between PCa and hub gene integrating clinical features Next, we further explored the impact of PPP1R14B combined with clinical factors according to TCGA datasets. The screening results of univariate COX regression analysis showed gleason grade (p<0.05, HR: 2.95, 95% CI: 1.3-6.5) (Fig. [92]7a). Subsequently, multivariate COX analysis was performed by combining gleason grade with PPP1R14B, and the results revealed that PPP1R14B (p<0.05, HR: 3.31, 95% CI: 1.2-9.4) and gleason grade (p<0.05, HR: 2.81, 95% CI: 12.2-6.3) (Fig. [93]7b) had independent prognostic value for overall survival of PCa. Figure 7. [94]Figure 7 [95]Open in a new tab Correlation between PCa and hub gene integrating clinical features based on the training cohort. a Univariate COX regression results, b multivariate COX analysis results, c prognostic nomogram to predict 1-, 3-, and 5-year overall survival, d ROC for 1-, 3- and 5-year overall survival, e-g calibration curves of prognostic nomogram for 1-, 3- and 5-year overall survival, h-j DCA of prognostic nomogram, hub gene and gleason-grade for 1-, 3- and 5-year overall survival. The nomogram was used to study the predictive prognostic value of hub gene PPP1R14B expression combining with clinical features and gleason - grade, showing that high PPP1R14B expression and gleason-grade had shorter overall survival (Fig. [96]7c). The predicted accuracy of the nomogram in both the training and validation cohorts was evaluated by the ROC curve, shown the different AUC values (training cohort: 1-year = 0.9931, 3-year = 0.6987, 5-year = 0.8195; validation cohort: 1-year = 0.993, 3-year = 0.7558, 5-year = 0.7439) (Fig. [97]7d, Fig. [98]8g). The calibration curves showed that the risks predicted by the nomogram were highly consistent with the observed overall survival for 1, 3, and 5 years in both cohorts (Fig. [99]7e-g, Fig. [100]8a-c). Moreover, DCA reflecting the performance of the prediction model confirmed that PPP1R14B expression could greatly predict overall survival and PPP1R14B combining with gleason grade had even better prediction ability in two cohorts (Fig. [101]7h-j, Fig. [102]8d-f). Figure 8. [103]Figure 8 [104]Open in a new tab Correlation between PCa and hub gene integrating clinical features based on the validation cohort. a-c Calibration curves of prognostic nomogram for 1-, 3- and 5-year overall survival, d-f DCA of prognostic nomogram, hub gene and gleason-grade for 1-, 3- and 5-year overall survival, g ROC for 1-, 3- and 5-year overall survival. Immune analysis We first focused on the abundance of 22 types of immune cells using the CIBERSORT algorithm (Fig. [105]9a&b). The results showed that PCa owned a higher abundance of M0 macrophage, M1 macrophage, M2 macrophage and T regulatory cells (Tregs) but had a lower abundance of plasma cells (Fig. [106]9c). In addition, the correlation analysis indicated that hub gene PPP1R14B had weak positive correlation (r < 0.4) with M0 macrophage, T cells CD4 memory resting and neutrophils but negatively correlated with plasma cells (|r| > 0.4) and was lowly negatively related (|r| < 0.4) to T cells follicular helper and mast cells resting (Fig. [107]9d). It was further confirmed that PPP1R14B was strongly positively correlated (r > 0.7) with immune checkpoints in ADSL, RAC3, MRPS2, IPO4, UCK2 and so on, and greatly negatively correlated (r < -0.7) with them in ACTG2, COL17A1, DKK3, PGRMC1, MYH11, CD40, SLC22A17, FERMT2, LPP and so on (Fig. [108]9e). What's more, the association between this hub gene and immune chemokine showed that PPP1R14B had stiff positive relevance (r > 0.7) to RPL29, CNOT10, APEX1, SHMT2, MARCKSL1, SCARB1, HPN, SND1, HOXC6, CCT3 and so on while the hub gene had hard negative correlation (r < -0.7) with ACTG2, COL17A1, COX7A1, DKK3, PGRMC1, MYH11, CD40, CRYAB, LPP and so on (Fig. [109]9f). Figure 9. [110]Figure 9 [111]Open in a new tab Immune analysis results. a-b Immune infiltrating cells estimated by CIBERSORT algorithm, c boxplot displaying differences in scores for infiltrating immune cells, d-f lollipop plot showing the correlation between hub gene and immune cells, checkpoints and chemokines. GSEA analysis GSEA analysis results showed that samples with high expression of PPP1R14B could significantly enrich pathways such as Myc targets V2, Myc targets V1, transcription of E2F targets under negative control by dream complex, E2F targets, G2M checkpoint, uronic acid metabolic process, GCM- PFN1 and so on (Fig. [112]10a-g). Figure 10. [113]Figure 10 [114]Open in a new tab GSEA analysis results. a Myc targets V2, b Myc targets V1, c transcription of E2F targets under negative control by dream complex, d E2F targets, e G2M checkpoint, f uronic acid metabolic process, g GCM- PFN1. Drug sensitivity analysis To explore the possibility of applying hub gene PPP1R14B to the personalized and accurate treatment of PCa, we analyzed sensitivity of different anticancer drugs according to the IC50 value. The results showed that 4 commonly used drugs selumetinib vorinostat, zebularine, azacitidine and VER155008 targeting PPP1R14B showed satisfactory differences in sensitivity, that is, patients in the high-risk group were more likely to be sensitive to them (Fig. [115]11a-h). Figure 11. [116]Figure 11 [117]Open in a new tab Drug sensitivity results targeting hub gene. a-d IC50 of selumetinib and vorinostat, zebularine, azacitidine and VER155008 in low and high expression groups of hub gene, e-h scatter plots between expression of hub gene and IC50 of selumetinib and vorinostat, zebularine, azaacitidine and VER155008. Construction of ceRNA Network The ceRNA regulatory network was constructed to investigate the regulatory mechanism of PPP1R14B. We used miRWalk and TargetScan databases to predict miRNAs associated with PPP1R14B and obtained 21 miRNAs from these two databases. Then, based on the obtained miRNAs, we searched for relevant lncRNAs in the miRNet database. Among the 21 miRNAs, 4 miRNAs (hsa-miR-1914-3p, hsa-miR-4640-5p, hsa-miR-4726-5p, and hsa-miR-5194) were found to be associated with 94 lncRNAs, and a total of 204 miRNA-lncRNA pairs were obtained (Fig. [118]12). Figure 12. [119]Figure 12 [120]Open in a new tab CeRNA network built with miRWalk, TargetScan and miRNet. Every box stands for a miRNA, and each circle represents a lncRNA. Discussion PCa is the second most diagnosed solid-organ cancer, after lung cancer, in men [121]^11. PCa has the complicated pathogenesis, and is easy to metastasize, which is challenging to treat. The available prostate-specific antigen (PSA)-based diagnostic and target-based therapeutic strategies approaches have come up with various false positives and off-target side effects in medical diagnosis and therapeutics of PCa [122]^12. As a result, exploring the pathogenesis and identifying the biomarker of PCa is highly critical to therapy. In this study, GEO gene expression datasets were employed by carrying out DEG, WGCNA and LASSO analyses, and furtherly hub gene was screened by COX regression based on TCGA datasets and validated through qRT-PCR in PCa and normal cells. Nomograms, calibration curves, DCA, ROC analysis and COX regression were used to explore relationship between PCa and hub gene, as well as correlation between PCa and hub gene binding clinical features. We delved potential mechanism of hub gene in PCa using GSEA, immune analysis and ceRNA networks. Finally, through the drug sensitivity analysis, we screened therapeutic drugs targeting hub gene related to PCa. Totally, 188 intersection genes were screened out according to GSE 32571 via differentially expression and WGCNA analyses, and significantly enriched in extracellular matrix structural constituent, glutathione derivative biosynthetic process, protein homodimerization activity, focal adhesion and so on. Extracellular matrix structural constituent and immune cell are important parts of the distinct tumor microenvironment (TME) of PCa, which take part in tumor proliferation and progression of PCa [123]^13. Glutathione is an antioxidant that acts as a free radical scavenger and a detoxifying agent in cells, and useful in a multitude of biological processes including cellular division, differentiation and proliferation and so on, and preserves sufficient levels of cysteine and detoxifies xenobiotics, which can promote the survival and proliferation of tumor cells and also confer resistance to cancer cells [124]^14. Protein homodimerization activity can increase the proliferation ability of cancer cells, reduce cell apoptosis, and promote invasion and metastasis [125]^15. Focal adhesion can regulate cell proliferation, survival and migration, and promote tumor growth [126]^16. In summary, the enrichment results revealed that 188 intersection genes had close relationship with PCa. Based on intersection genes, 32 key genes were further screened through LASSO analysis. Furthermore, among them, we found that PPP1R14B had a remarkable relationship with PCa. Both qRT-PCR and TCGA datasets indicated that hub gene PPP1R14B highly expressed in PCa groups. GEO datasets showed that PPP1R14B had area under the curve (AUC) value > 0.95 that were more than known PCa-related genes AMACR and GDF15, displaying ideal predictive performance for PCa occurrence, as assessed by the ROC curves, and the nomogram models based on hub gene expression level further confirmed this point. Subsequently, COX univariate regression analysis was performed based on clinical sample information obtained from the TCGA database. The hazard ratio (HR: 3.47, 95% CI: 1.3-9, p<0.05) and survival curve of PPP1R14B both reveal that PPP1R14B is a high-risk factor affecting the occurrence and prognosis of PCa. Furthermore, we continued to explore the impact of PPP1R14B combined with clinical factors on PCa. Through univariate COX regression analysis, gleason grade was identified as a clinical factor significantly affecting the prognosis of PCa. Subsequently, multivariate COX regression analysis was conducted based on PPP1R14B combined with gleason grade, and the results showed that both PPP1R14B and gleason grade have a strong correlation with the prognosis of PCa and play an important role in the occurrence and development of PCa. Subsequently, we established the prognostic nomogram of PPP1R14B and gleason grade based on TCGA datasets and verified its excellent prognostic predictive performance through calibration curves and DCA in training and validation cohorts. The results additionally showed that PPP1R14B could function as a vital prognostic factor in PCa. PCa is a multifactorial and complex disease involving several immunological, environmental, physiologic, and genetic factors [127]^9, so we explored the relationships between PCa and immune infiltration. In our study, PCa owned a higher proportion of M0 macrophage, M1 macrophage, M2 macrophage and Tregs but relatively lower ones of plasma cells. Mononuclear cells are recruited to the TME by stromal cells and tumor-secreted chemokines and growth factors secreted by the tumor to differentiate into M0 macrophages [128]^17. M0 macrophages display two phenotypes, namely M1 and M2. M0 polarization is driven toward the M1 phenotype by miR-203, where M1 macrophages actively present antigens to effector cells in the immune system and secreted multiple pro-inflammatory factors, thereby influencing inflammatory responses and anti-tumor immunity during tumor progression [129]^18. M2 macrophages in tumor tissue release cytokines epidermal growth factor and transforming growth factor β1 (TGF-β1), which facilitate tumor cell division, proliferation, and growth [130]^18. Thus, an increase in the abundance of M1 macrophage could promote inflammatory response and anti-PCa effects, thereby promoting PCa cell apoptosis, but an increase in the proportion of M2 macrophage would inhibit these processes. In previous vivo and vitro studies have shown that Tregs are another type of immune cell that suppress natural killer (NK) cells and CD8^+ T lymphocytes by secreting inhibitory cytokines (such as transforming growth factor β, TGF β) [131]^19, which is directly related to a greater likelihood of death from PCa [132]^20. Plasma cell is a kind immune cell that can generate and secret antibody [133]^21. When the proportion of plasma cells decreased, the production of antibodies might be limited, which would lead to impaired immune response against prostate cancer. In summary, our results revealed that inhibition of anti-tumor immunity in PCa was closely associated with the increase of macrophage and Tregs and the decrease of plasma cells, which might lead to improvement of PCa cells anti-apoptosis, activity inhibition of NK cells and CD8+ T cells, and strength frustration of humoral immunity. Furthermore, we explored the correlation between hub gene PPP1R14B and the TME. We found that PPP1R14B had highly negative correlation with immunologic factors CD40, COL17A1, DKK3 and PGRMC1, but had strongly positive relevance to ALKBH2, UCK2, RAC3 and RAB17. CD40, a member of the tumor necrosis factor receptor family, is expressed by B cells, professional antigen-presenting cells, as well as non-immune cells and tumors [134]^22. CD40 signaling of B cells promotes germinal center formation, immunoglobulin (Ig) isotype switching, somatic hypermutation of the Ig to enhance affinity for antigen, and finally the formation of long-lived plasma cells and memory B cells [135]^23. Tests on mice show that in B cells, suppression of ALKBH2 activates NF-κB [136]^24 that is one of the transcription factors of CD40 and promotes the transcription and expression of CD40 [137]^25, which can enhance the degree of differentiation of B cells into plasma cells and promote antibody synthesis. It was demonstrated that downregulation of UCK2 also activated the TNF α/NF-κB signaling pathway [138]^26. Our results revealed that the high expression of PPP1R14B might lead to a decrease in CD40 through direct regulation or possible ALKBH2/UCK2-NF-κB-CD40 signaling, thereby hindering differentiation of B cells and formation of plasma cells, which resulted in insufficient expression of antibodies and strength frustration of humoral immunity directed against PCa. Additionally, research suggested that DKK3 could upregulate the antitumor immune function of CD56^bright NK cells by inducing its differentiation and improving its cytotoxicity [139]^27. The reduced expression of DKK3 can stimulate TGF-β signaling [140]^28, which can help naïve T cells differentiate to Tregs [141]^20, thereby preventing NK and CD8+ T cells from exerting anti-tumor effects [142]^19. UCK2, uridine-cytidine kinase 2, is a key regulator of pyrimidine metabolism, and the increased expression of UCK2 partly stimulated TGF β1, which also inhibited NK cell infiltration and NK-cell-mediated killing [143]^26. When the expression level of PPP1R14B increased, it might cause the downregulation of DKK3 and the upregulation of UCK2, which would facilitate TGF β signaling and further Tregs generation, thereby reducing the cytotoxicity of NK cells and promoting CD8+ T cells exhaustion to inhibit cellular immunity in PCa. What's more, RAC3, a member of the Rac family of small guanosine triphosphatases (GTPase), regulates a variety of cell functions, including the organization of the cytoskeleton, cell migration, and invasion [144]^29, and macrophage morphological plasticity and migration is Rac GTPase signalling dependent [145]^30. RAB17 is a subtype of the Rab GTPase family that can promote M2 polarization through up-regulating PPARγ induced glutamine metabolism [146]^31. The high expression of PPP1R14B might lead to the upregulation of RAC3 and RAB17, further causing to increase the proportion of M2 macrophage to impede the apoptosis of PCa cells. In addition, some immunologic factors are related to regulation of neutrophils. COL17A1 is type XVII collagen that regulates the expression level of the pro-inflammatory chemokine Interleukin-8 (IL-8) that has neutrophil chemotactic activity to inhibit the growth and spread of tumor cells [147]^32^, [148]^33. Recent studies in hepatocellular carcinoma showed that PGRMC1, a heme-binding protein, belonged to the membrane-associated progesterone receptor (MAPR) family of cytochrome b5-related proteins [149]^34 and can enhance IL-8 production in multiple immune cells and promotes the recruitment of neutrophils [150]^35. Thus, increased expression of PPP1R14B could cause a reduction in COL17A1 and PGRMC1, which would hinder immune functions of neutrophils. In summary, overexpression of PPP1R14B could result in the upregulation of ALKBH2, UCK2, RAC3 and RAB17 and downregulation of CD40, DKK3, COL17A1 and PGRMC1, which would stifle humoral immunity against PCa through declining the level of plasma cells, block cellular immunity based on NK and CD8+ T cells by upregulating Tregs, destroy cellular immunity according to macrophages via increasing the proportion of M2 macrophages and repress cellular immunity based on neutrophils by means of blocking neutrophil chemotactic activity conclusively inhibiting adaptive immune directed against PCa. Subsequently, pathway enrichment indicated that the high expression samples of PPP1R14B were mainly related to MYC targets V2, MYC targets V1, transcription of E2F targets under negative control by dream complex, E2F targets, GCM-PFN1, G2M checkpoint, uronic acid metabolic process and so on. PPP1R14B, Protein Phosphatase 1 Regulatory Inhibitor Subunit 14B, is protein phosphatase 1 (PP1) inhibitor that can inhibit the monomer catalytic subunit of PP1, further participate in various phosphorylation dependent signaling pathways, and regulate various cellular processes [151]^36^, [152]^37. The enrichment results showed that the overexpression of PPP1R14B is closely related to the body's immunity. The results revealed that the high expression of PPP1R14B is involved in the regulation of MYC and E2F, which are two important transcription factors. Research has shown that PPP1R14B can inhibit PP1, which leads to excessive phosphorylation of MYC on multiple serine and threonine residues, thereby reducing MYC protein levels [153]^36^, [154]^38. There is a mutual interaction between MYC and NF-κB [155]^39. In tumor cells, MYC can serve as a positive regulatory factor of NF- κB [156]^40, which is one of the transcription factors of CD40 [157]^25. Further, Studies have shown that MYC can reverse regulate the expression of GEFs [158]^41, which is a type of GTPases activator and can activate RAC3 [159]^29^, [160]^42. Therefore, PPP1R14B can inhibit CD40 and promote RAC3 through PP1/MYC/NF-κB and PP1/MYC/GEFs signaling, respectively, thus suppressing humoral immunity and cellular immunity based on macrophages. Moreover, PP1 can connect PP1 regulatory subunit NIPP1 and PKA to restrain the expression of E2F1 gene [161]^43, which is one of the most important E2F subtypes [162]^44. Research has shown that through differentiation specific marker type 1 transglutaminase (TG-1), E2F can mediate the inhibition of downstream target Sp1, which is a transcription factor for COL17A [163]^45^, [164]^46. It can be concluded that PPP1R14B can restrain COL17A through PP1/E2F/SP1 signaling, therefore playing an inhibitory regulatory role in neutrophil mediated cellular immunity. Furthermore, inhibition of PP1 by PPP1R14B can promote the phosphorylation of GSK3β by disrupting the phosphorylation-dephosphorylation balance regulated by AKT2 and PP1, hence activating GSK3β [165]^47, which can inhibit the activation of β-catenin signal [166]^44. β-catenin is a key component of the Wnt signaling pathway, which regulates the transcription of genes such as DKK3 by interacting with transcription factors of the TCF/LEF family [167]^48. Hence, PPP1R14B can inhibit DKK3 through PP1/GSK3β/β-catenin signaling, thus playing a regulatory role on activating Treg cells and further inhibiting CD8+T and NK cell mediated immunity. In summary, PPP1R14B regulates the activity of multiple factors including MYC, E2F and GSK3β by inhibiting the dephosphorylation of PP1, which will adjust the expression of downstream factors such as CD40, RAC3, COL17A, DKK3 and so on. Therefore, overexpression of PPP1R14B can lead to dysregulation of the above signaling, consequently inhibiting adaptive immunity against PCa. In addition, the enrichment results also showed that the pathways related to the overexpression of PPP1R14B were involved in the proliferation, apoptosis and migration of tumor cells. The dysregulated E2F stimulates the transcription of pro apoptotic and tumor suppressor genes and induces cell apoptosis, thereby protecting cells from tumor invasion [168]^49. E2F1 promotes the invasion and migration of PCa cells by regulating CD147. Overexpression of E2F1 can predict the adverse prognosis of human PCa [169]^50. G2M is an important stage in the cell cycle associated with mitosis. Inhibition of PP1 leads to phosphorylation of Cyclin Dependent Kinases, furtherly inactivating pocket proteins through high phosphorylation and hence promoting cell cycle progression [170]^51. which will stimulate the survival and proliferation of tumor cells [171]^52. Besides, Profilin 1 (PFN1) can adjust the migration of tumor cells, including PCa [172]^53. In conclusion, overexpression of PPP1R14B can promote the proliferation and migration, as well as inhibit apoptosis of PCa cells based on E2F, G2M, PFN1 etc. related pathways. Recent efforts in cancer therapeutics have focused on the development of targeted therapies directed toward specific molecules related to cancer cell biology to achieve a more specific and effective treatment. Therefore, we evaluated the correlation between hub gene PPP1R14B expression and commonly used drugs sensitivity and observed that patients with higher PPP1R14B expression responded more sensitively to these drugs selumetinib and vorinostat, zebularine, azacitidine and VER155008. The combination of the mitogen-activated protein kinase/extracellular signal-regulated kinase 1/2 inhibitor, selumetinib, and the histone deacetylase inhibitor, vorinostat exerts potent antitumor effects on in vitro models of colorectal cancer, including synergistic inhibition of cell proliferation and migration, G1 cell-cycle arrest, induction of apoptosis, and spheroid formation in 3D culture [173]^54. Research suggested that zebularine led to increased levels of interferon response, and promoted infiltration of CD8+T and NK cells into tumor and therefore suppressed tumor growth [174]^55. Effects of azacitidine include direct cytotoxicity from inhibition of protein synthesis and DNA damage and re-expression of aberrantly silenced tumor suppressor genes due to DNA hypomethylation [175]^56^, [176]^57. The heat shock protein 70 inhibitor VER155008 suppresses the expression of HSP27, HOP and HSP90β and the androgen receptor, induces apoptosis, and attenuates prostate cancer cell growth [177]^58. This suggests that PPP1R14B may provide a new target for choosing appropriate drugs for PCa patients. More and more evidence suggests that ceRNA plays a crucial role in cancer. However, there is still a lack of ceRNA networks that can regulate PCa. In this study, we identified 21 miRNAs associated with PPP1R14B, which play important roles in various cancers. Research has found that hsa-miR-6090 is upregulated in PCa patients, which may be related to the occurrence of PCa [178]^59. What's more, has-circ-0044520 and has-circ-0044529 can exist stably in cells for a long time. This stability can block has-miR-4726-5p and has-miR-4640-5p, which can potentially activate downstream oncogenes and initiate the development of laryngeal squamous cell carcinoma [179]^60. It is reported that hsa-miR-134-3p plays a key role in the development of non muscle invasive bladder cancer [180]^61. However, the role of these miRNAs in PCa is still unclear. Furthermore, we also reverse predicted 94 lncRNAs associated with miRNAs, some of which are closely related to PCa. Some studies have shown that LINC00963 promotes the transition of PCa from androgen dependent to androgen independent mode by participating in the trans activation of EGFR [181]^62. More and more evidence suggests that TERC plays an important role in telomere maintenance and other functions in human cancer. Knocking down the TERC by siRNA can reduce proliferation of PCa cell lines. In addition, studies have found a correlation between TERC levels and MYC levels in clinical samples, with forced overexpression of MYC leading to an increase in TERC levels; After MYC silencing, the activity of the human TERC promoter decreases [182]^63. DANCR targeting miR-185-5p via FAK/PI3K/AKT/GSK3 β/snail pathway upregulates LIM and SH3 protein 1, thereby affecting cell migration and proliferation, promoting prostate cancer [183]^64. In addition, it has been reported that HCG11 is abnormally expressed in various tumors and can serve as a key regulatory factor for various cancers such as PCa [184]^65, non-small cell lung cancer [185]^66, and hepatocellular carcinoma [186]^67. The expression level of HCG11 in PCa is significantly reduced, which is directly proportional to the patient's survival rate [187]^68. The ceRNA network indicates the potential regulatory mechanisms of ceRNA on PPP1R14B and PCa. Although some research has been conducted on the regulatory mechanisms in PCa, further exploration is still needed in the future, especially in the ceRNA network. Conclusion Although PPP1R14B may have a certain correlation with the prognosis of uterine corpus endometrial carcinoma, breast cancer and gastrointestinal cancer [188]^69^-[189]^71, its relationship with PCa has not been studied yet. In our study, hub gene PPP1R14B was screened out through WGCNA, differentially expression analysis, LASSO and COX regression, and had significantly predictive values of diagnosis and prognosis in PCa. Immune and pathway enrichment analysis revealed that overexpression of PPP1R14B could inhibit the dephosphorylation of PP1 to result in the regulation of multiple signaling MYC, E2F, GSK3β and so on in a mess, which accelerated the disorder of several immune- and tumor-related factors, ultimately leading to inhibition of adaptive immunity directed against PCa and promotion of PCa cell proliferation and migration. In summary, our study first proposed that PPP1R14B was a biomarker for evaluating immune characteristics, proliferation and migration in PCa, and had a close relationship with development and prognosis of PCa, which would provide useful insights and a theoretical basis for creating personalized and accurate treatment strategies and drug decisions for patients with PCa. Supplementary Material Supplementary figure. [190]jcav15p6545s1.pdf^ (155.3KB, pdf) Acknowledgments