The transcriptome data [the level 3 mRNA expression data (FPKM), normalized using log 2 ( F P K M + 1 ) ] of normal tissues (52 cases) and tumour tissues with complete clinical information (379 cases) were extracted from The Cancer Genome Atlas (TCGA) database of prostate adenocarcinoma (PRAD). The mRNA expression profiles contained in the GSE116918 (Jain et al., 2018), GSE29079 (Börno et al., 2012), and GSE6956 (Wallace et al., 2008) datasets, which were normalized by their corresponding providers, were downloaded from Gene Expression Omnibus (GEO) database. A total of 248 PCa cancer samples with clinical information were included in the GSE116918 dataset. The GSE29079 dataset contained 48 normal samples and 47 PCa samples, while the GSE6956 dataset had 18 normal samples and 69 PCa samples. However, neither GSE29079 nor GSE6956 contains clinical information. The BioGRID database offered 253 unique interactors of pinin with experimental pieces of evidence (Oughtred et al., 2021). TSVdb offered PNN splicing variants expression (Sun et al., 2018). For PNN expression in pan-cancer, we downloaded the standardised pan-cancer dataset TCGA TARGET GTEx (PANCAN, N = 19131, G = 60499) from the UCSC (https://xenabrowser.net/) database and further extracted the expression data of PNN gene in each sample. In addition, we filtered out the samples with zero expression levels, and further transformed each expression value with log2 (x + 0.001), finally, we excluded those with less than three samples in a single cancer species.

The Human Protein Atlas (HPA) provides the protein expression of pinin in normal prostate (via https://www.proteinatlas.org/ENSG00000100941-PNN/tissue/prostate) and tumour tissues (via https://www.proteinatlas.org/ENSG00000100941-PNN/pathology/prostate+cancer) (Uhlén et al., 2015). All images of tissues in HPA database are stained by immunohistochemistry. We extracted the immunohistochemistry images directly from the HPA database.

Correlation analysis of PNN expression and clinicopathological characteristics was performed. The expression of PNN between the subgroups was compared based on the following clinicopathological features: age (<60 or ≥60 years old), N stage (N0, N1), M stage (M0, M1), T stage (T2, T3, T4), surgical margin (R0, R1, R2, RX), prostate-specific antigen (PSA) level (<10 or ≥10 years), and Gleason score (6, 7, 8, 9, 10). Univariate and multivariate Cox regression analyses were implemented to identify independent predictors of survival in the TCGA-PRAD and GSE116918 datasets.

We downloaded GSE38241 (Aryee et al., 2013) and GSE25136 (Sun and Goodison, 2009) datasets (the authors processed normalisation) from GEO. For the merging of these datasets, we used the method of COMBAT (Johnson et al., 2007), implemented in the R package inSilicoMerging (Taminau et al., 2012) to obtain the expression matrix. Finally, the PNN expression was compared using the Kruskal-Wallis test.

We calculated the Pearson correlation of all genes (RNA-seq) in the TCGA dataset with PNN using the Linkomics database (http://www.linkedomics.org/) and selected the genes with correlation coefficients > 0.8 and p < 0.05 as PNN co-expressed genes.

The “clusterProfiler” R package was utilised to conduct Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis (Yu et al., 2012). GO enrichment analysis mainly described the biological processes (BP), cellular components (CC), and molecular functions (MF) correlated with genes. The threshold for significant enrichment was set as a p-value < 0.05 or FDR < 0.05, as stated. Single sample gene set enrichment analysis (ssGSEA) enrichment scores were calculated in each sample using the “GSVA” package of R (Hänzelmann et al., 2013).

In this research, potential drug (or molecules) was predicted using the Drug Signatures database (DSigDB) via Enrichr (https://maayanlab.cloud/Enrichr/) based on the PNN gene as well as the positively co-expressed gene with PNN (correlation coefficient > 0.8 and p < 0.05) (Chen et al., 2013; Kuleshov et al., 2016; Xie et al., 2021).

The CpG sites in the promoter of PNN and PNN’s co-expressed genes were obtained from the MEXPRESS database (Koch et al., 2015; Koch et al., 2019). A univariate Cox analysis in R was used to determine the association between methylation levels at each CpG site and progression-free survival (PFS) for each patient, and p < 0.01 was considered statistically significant. Candidate prognostic CpG sites were selected using the Least Absolute Shrinkage and Selection Operator (LASSO) algorithm. Based on the candidate CpG sites generated from the above algorithm, a multivariate Cox regression model was used to construct a prognostic signature. The RiskScore of each recipient was calculated using the following formula: R i s k S c o r e = Σ i = 1 n   β i × M e t h i In which β refers to coefficient, and M e t h refers to the level of methylation.

Patients were divided into the high-risk ( R i s k S c o r e ≥ m e d i a n ) and low-risk groups ( R i s k S c o r e < m e d i a n ) in the TCGA dataset. Then, we performed ROC analysis using the R software package pROC (version 1.17.0.1) to obtain the AUC. The R package “survival” was used to perform the two risk groups’ Kaplan-Meier (KM) survival analysis.

To inspect the different signalling pathways between the PNN low- and high-expression groups in the TCGA-PRAD dataset, Gene Set Enrichment Analysis (GSEA) was conducted by the “clusterProfiler” package in R software (Subramanian et al., 2005). Pathways with a p-value < 0.05 were considered significantly enriched.

The Tumor and Immune System Interaction Database (TISIDB) (http://cis.hku.hk/TISIDB) database was utilised to analyse the associations between PNN and tumour-infiltrating lymphocytes (TIL), immunosuppressors, and chemokines (Ru et al., 2019).

Statistical analysis was performed using the R software package (version 3.6.1). The differential mRNA expression of PNN between tumour tissues and normal controls was compared using Student’s t-test. The expression of PNN among the clinicopathological parameters groups was compared using Student’s t-test and ANOVA. The area under the curve (AUC) of receiver operating characteristic (ROC) was utilised to determine the diagnostic ability of PNN and was calculated using the “pROC” R package (Malone et al., 2015). KM curves of disease-free survival (DFS or PFS) of the patients were performed by setting the median expression of PNN as the cut-off in the ‘survival’ R package. The log-rank test was used to assess statistical differences, and a cut-off p-value < 0.05 was deemed statistically significant.

The expression levels of PNN between PCa and control samples were compared in the TCGA-PRAD, and the PNN expression level was validated with GSE29079 and GSE6956 datasets. As shown in the violin plots, the mRNA expression level of PNN was significantly higher in the PCa group in all datasets (Figures 2A–C). Next, we used the same datasets to evaluate the diagnostic value of the PNN gene. The accuracy of the diagnostic model was evaluated by ROC curve analysis (Figure 2D). As a result, the AUC of the PNN diagnostic model was greater than 0.7 in all three datasets, indicating that the PNN gene can be used to discriminate cancer from normal tissues. Moreover, we also observed that the abundance of pinin protein was higher in PCa tissue than in normal tissue (Figures 2E,F).

To explore the relationship between PNN expression and the clinicopathological characteristics in PCa, we compared the PNN expression levels according to sample clinical information. The high PNN expression was found in the advanced stage of PCa (Figure 3B), and the Gleason scores were strongly correlated with the PNN expression levels in PCa patients in both TCGA-PRAD datasets (p = 6.3 × 10 − 9 ) and GSE116918 dataset (p = 0.001) in Figures 3E,I. Collectively, the Gleason score was highly positively correlated with PNN expression. Different the surgical margins (R0/1/2/X) found different PNN expression (Figure 3D). It has been found that the PNN gene expression level was significantly higher in tumors than that of the primary tissue (Figure 3J, data process in Supplementary Figure S1), suggesting this gene can be used for diagnostics in metastatic patients. Age (Figures 3A,F), T stage (Figures 3C,H), or PSA level (Figure 3G) are not correlated with the PNN expression’s significance.

Univariate and multivariate Cox analyses were conducted to investigate the independent prognostic factors in TCGA-PRAD and validated with GSE116918 datasets. The univariate analysis in the TCGA-PRAD dataset indicated that the surgical margin, T stage, N stage, Gleason score, and PNN expression were associated with the prognosis of PCa patients (Figure 4A). In contrast, multivariate Cox regression analyses in the same dataset demonstrated that only the Gleason score could be used independently to predict the prognosis of patients ( Figure 4B). Similarly, the PSA levels, Gleason score, T stage, and PNN expression were found to be significant risk factors by univariate Cox analysis in the GSE116918 dataset (Figure 4C). In the same dataset, multivariate Cox regression analyses demonstrated that T stage and PNN expression could be used independently to predict the prognosis of patients (Figure 4D). We then validated these findings by analysing the DFS curves of the PNN high- and low-expression groups, which showed that the PNN high-expression group had remarkably worse survival rates than the low-expression group in both the TCGA-PRAD and the GSE116918 datasets (Figures 4E,F). The hazard ratio of PNN was greater than 1 in both datasets. Taken together, it suggested that PNN was a risk factor in the prognosis of PCa. However, the independent prognostic value of PNN needed further investigation and confirmation.

To identify pharmaceutical molecules with DsigDB database and further uncover the biological processes PNN participated, the co-expression pattern of PNN in PCa was explored. All co-expressed genes are listed in Supplementary Table S2.

BioGrid hosted 243 proteins interacting with pinin extracted from published literature. A total of 368 genes were co-expressed with pinin following the criteria of r > 0.6 and p < 0.05, of them, twenty-five genes overlapped with 243 interactive proteins of pinin (25UC for short). Those 25UC genes were enriched in RNA splicing and RNA/mRNA processing based on GO enriched analysis (Figure 5A) and enriched in the spliceosome, mRNA surveillance pathway, and RNA transport based on KEGG enrichment analysis (Figure 5B). These results suggest that PNN is mainly linked to the RNA process and RNA transport in PCa. PNISR, RBM39, DDX39B, SF3B1, SRSF11, CPSF6, CLK2, and SNRPB2 have the function of splicing or process of RNA; ACIN1 and NKTR participate in cell apoptosis and immune response. The protein-protein interaction network can be found in Figure 5C.

To explore the potential therapeutic targets in PCa, we focused on those genes that strongly positively (r > 0.8 and p < 0.05) correlated with upregulated PNN, including FNBP4, TCERG1, RBM39, DDX39B and DMTF1. Ten possible pharmaceutical molecules were identified using the Enrichr package from the DsigDB database, based on their p-value. Table 1 lists the effective drugs from the DsigDB database for PCa.

After excluding missing values, a total of 180 CpG sites in the PNN and its co-expressed genes FNBP4, TCERG1, RBM39, DDX39B and DMTF1 (r > 0.8 and p < 0.05) promoter regions were extracted from TCGA-PRAD methylation data. Univariate Cox regression analysis showed that 25 CpG sites were significantly correlated with PFS. Following the LASSO algorithm, 16 CpG sites were selected (lambda = 0.009914324, Figures 6A,B). A model was then constructed with multivariate Cox regression. We constructed a risk score system Eq. 1 with seven CpG sites. R i s k   S c o r e = − 17.451 × c g 04787786 − 2.182 × c g 09878914 + 3.097 × c g 16316344 − 38.055 × c g 16408528 + 1.652 × c g 17114847 − 46.828 × c g 17439097 + 1.943 × c g 25800328 (1)

Eq. 1

The areas under the ROC curves (AUC) of 1-, 2-, and 3-year PFS were 0.80, 0.74 and 0.71, respectively (Figure 6C), indicating the good performance of the risk score signature. We noticed that this risk score was linked to the PFS status of the PCa patients (Figure 6D), indicating that this risk score could be used to predict the progression of PCa. Multivariate Cox regression confirmed that the risk scores could also be an independent prognostic factor (Figures 6E,F). In addition, the expression level of 47 immune checkpoint genes (ICG) proposed by Danilova et al. (2019) were compared between high and low-risk groups with the Wilcoxon test based on the signature constructed above. As a result, CTLA4, CD276, CD80, NRP1, TNFRSF18, TNFSF18 and TNFRSF14 were found to be significantly higher in the high-risk group, while the expression levels of BTNL2, CD160, CD200, CD244, CD274, CD40, CD44, LAG3, TMIGD2, TNFSF15, VSIR and VTCN1 were reduced significantly in the high-risk group (Figure 6G). Then, a ssGSEA was performed using the KEGG database to explore different molecular mechanisms between the high- and low-risk groups. Among significantly enriched pathways (p < 0.05), the top 10 were compared between high and low-risk groups. Between the two risk groups, the splicing factors genes, such as SF3B1, SRSF11 and SRSF7, are significantly higher in high risk group based on the median (Figure 6H). The expression of SRSF10 (p = 0.05) and SRSF4 (p = 0.07) were marginally higher. Moreover, the splicing isoforms expression was also significantly increased in the highly risky group based on the median. Besides prostate cancer and other cancer pathways, the risk model also found significant different enrichment scores in spliceosome and biogenesis and degradation pathways (Figure 6I). The correlation between the pathways and the PNN expression is illustrated in Figure 6J.

To infer the pathways by which PNN genes were involved in the development of PCa, GSEA enrichment was performed on the PNN high- and low-expression groups. Among the enriched pathways (adjusted p < 0.05) (Supplementary Table S1), we noted that immune-related pathways were enriched, including the IL-17 signalling pathway, the T cell receptor signalling pathway, Th1 and Th2 cell differentiation, Th17 cell differentiation, and the TNF signalling pathway. In addition, cancer-related pathways, such as the cell cycle, choline metabolism in cancer, PD-L1 expression and PD-1 checkpoint pathway in cancer, and proteoglycans in cancer were also enriched in PCa (Figure 7A). We calculated the expression difference between normal and tumour samples in each tumour, and observed significant upregulation in 14 tumours (Figure 7B). Subsequently, the correlation between PNN and immune infiltration was executed to broaden the cognition of the correlation between PNN and TIL, immune inhibitors, and chemokines in PCa. As to TIL, PNN expression was negatively correlated with iDC, monocyte, NK cell, and Tgd in Figure 7C (rho < −0.3 and p < 0.05). Figure 7D showed the correlations between PNN expression and chemokines, of which CCL14 was negatively correlated with PNN (r < −0.3 and p < 0.05).