Prognosis and diagnosis of prostate cancer based on hypergraph regularization sparse least partial squares regression algorithm

Data source

We downloaded gene expression data of prostate cancer patients (TCGA-PRAD) from the TGCA database (https://portal.gdc.cancer.gov/), which included 499 prostate cancer tissue samples and 52 normal tissue samples. The methylation data and corresponding clinical information of prostate cancer patients in the TCGA-PRAD cohort were downloaded from the UCSC Xena database (https://xenabrowser.net/datapages/). In addition, we downloaded the GSE116918 dataset from the Gene Expression Omnibus (GEO) database for external validation of the prognostic model. A total of 248 patients with prostate cancer who received radical radiotherapy were included in the GSE116918 dataset. We divided the training set and test set by the ratio of about 7:3 on the samples of The Cancer Genome Atlas (TCGA) dataset. Finally, 428 training set samples and 123 test set samples were obtained.

This paper uses transcriptome and methylation data from the same batch of prostate cancer samples for association analysis. The weight vectors of the two data are obtained through the proposed HR-SPLS algorithm. Then the modules are divided according to the set number of modules. Modules with smaller objective function values were selected for various bioinformatic analyses. Finally, the disease samples of the test set and the independent test set GSE116918 data set were used to verify the prognostic-related genes and the model.

Partial least squares regression (PLS)

The partial least squares algorithm can simultaneously model multiple independent and dependent variables, especially when multicollinearity exists in both. The objective function is as follows.

$$\underset{g,d}{\text{max}}\text{ cov }(Xg,Yd)$$

$$\text{s.t. }{g}^{T}g=1,d^{T}d=1$$ (1)

Among them, $X\in {ℝ}^{n\times p}$ represents the expression matrix of the first data. $Y\in {ℝ}^{n\times q}$ represents the expression matrix of the second data. n represents the total number of samples. p and q represent the first and second data characteristics, respectively. $\text{cov }(·)$ represents the difference between the co -party. g and d represent the two right vectors. Further, this article introduces two potential variables: u and v. Among them, u = Xg, v = Yd. Formula (1) The covariance between the two potential variables u and v through maximizing the two potential variables u and v.

Sparse partial least squares regression (SPLS)

PLS does not meet the needs for high-dimensional biological chip data analysis. Therefore, SPLS is proposed to solve the feature selection of high-dimensional data. SPLS adds model punishment items to the model of the suitable vector and, based on PLS, helps the algorithm selection of more representative and essential features. The target function of the SPLS algorithm is shown below.

$$\underset{g,d}{\text{max }}\text{ cov }(Xg,Yd)-{\lambda }_1\parallel g{\parallel }_1-{\lambda }_2\parallel d{\parallel }_1$$

$$\text{s.t. }{g}^{T}g=1,d^{T}d=1$$ (2)

Among them, λ₁ and λ₂ control the constraint strength of the l1 norm of the weight vectors g and d, respectively. It can be used to select variables with better biological interpretability.

Hypergraph learning

The simple graph can represent the pair relationship between the objects. The vertex can be expressed as an object, and the edge represents the relationship between the apex. However, the complex relationship may not be represented in a simple graph, which may cause information loss. Hypergraph can connect to two or more vertices through the hyper edge. As an extension of a simple graph, each side of the hyper edge can be connected to multiple vertices, called a hypergraph. G (V, E, w) represents hypergraph. Among them, $V=\{v_1,v_2,\dots ,v_N\}\in {ℝ}^N$represents the vertex in the hypergraph. $E=\{e_1,e_2,\dots ,e_M\}\in {ℝ}^M$ represents the hyper edge in the hypergraph. $w={(w(e_1),w(e_2),\dots ,w(e_M))}^T\in {ℝ}^M$ is the weight of E. Next, this paper introduces the associated matrix H to characterize the relationship between V and E. The element at row i and column j in H can be expressed as:

$$H_{ij}=\{\begin{array}{ll}1,\hfill & \text{ if }{v}_i\in e\hfill \\ 0,\hfill & \text{ if }{v}_i\notin e\hfill \end{array}$$ (3)

Further, define the degree matrix D_v and D_e of the edge and vertex and the diagonal matrix W, as shown below. In addition, this article gives a simple hypergraph example (Figure 1A, 1B).

$$D_v=\left(\begin{array}{ccc}{\sum }_{e_j\in ℰ}w\left(e_j\right)H_{1j}& \cdots & 0\\ ⋮& \ddots & ⋮\\ 0& \cdots & {\sum }_{e_j\in ℰ}w\left(e_j\right)H_{Nj}\end{array}\right)$$ (4)

$$D_e=\left(\begin{array}{ccc}{\sum }_{v_i\in \mathcal{V}}{H}_{i1}& \cdots & 0\\ ⋮& \ddots & ⋮\\ 0& \cdots & {\sum }_{v_i\in \mathcal{V}}{H}_{iM}\end{array}\right)$$ (5)

$$W=\left(\begin{array}{ccc}w\left(e_1\right)& \cdots & 0\\ ⋮& \ddots & ⋮\\ 0& \cdots & w\left(e_{\text{M}}\right)\end{array}\right)$$ (6)

Then this paper defines the similarity matrix S of the hypergraph G and the Laplacian matrix L_H of the hypergraph.

$$S=HW{D}_e^{-1}{H}^T\in {ℝ}^{N\times N}$$ (7)

Similar to the definition of the symmetrical Laplace matrix of the simple graph, the symmetry of the hypergraph is defined below.

$$L=I-D_v^{-\frac{1}{2}}S{D}_v^{-\frac{1}{2}}$$ (8)

Sparse partial least squares regression algorithm based on hypergraph regularization (HR-SPLS)

Hypergraphs can characterize a high-level relationship between complex objects. In this paper, we innovatively added the hypergraph to the SPLS algorithm as the priority information and proposed the HR-SPLS algorithm. Before defining the target function of this paper, first of all, the hypergraph definition of this article is given.

$$\Omega (g)=g^T{L}_{H1}g$$ (9)

$$\Omega (d)=d^T{L}_{H2}d$$ (10)

Among them, L_H1 and L_H2 are hypergraph Laplace matrix.

$$L_{H1}=I-D_{v1}^{-\frac{1}{2}}{S}_1{D}_{v1}^{-\frac{1}{2}}$$ (11)

$$L_{H2}=I-D_{v2}^{-\frac{1}{2}}{S}_2{D}_{v2}^{-\frac{1}{2}}$$ (12)

D_v1 and D_v2 represent the degree matrix of X and Y, respectively, and S₁ and S₂ represent the similarity matrix of X and Y, respectively. We can get the target function of the HR-SPLS algorithm.

$$\begin{array}{c}\underset{g,d}{\text{max }}\text{ cov\hspace{0.17em}}(Xg,Yd)-{\beta }_1\Omega (g)-{\beta }_2\Omega (d)\\ -{\lambda }_1\parallel g{\parallel }_1-{\lambda }_2\parallel d{\parallel }_1\end{array}$$

$$\text{s.t. }{g}^{T}g=1,d^{T}d=1$$ (13)

Among them, β₁ and β₂ control the strength of the hypergraph regularization respectively. λ₁ and λ₂ control the sparseness of the two weight vectors, respectively. Further, we can rewrite it:

$$\begin{array}{c}\underset{g,d}{min}-\frac{1}{p}{g}^T{X}^{T}Yd+{\beta }_1\\ {\displaystyle \sum _{1\le i<j\le n}{s}_{1_{ij}}{\left(\frac{g_i}{\sqrt{l_{H{1}_i}}}-\frac{g_j}{\sqrt{l_{H{1}_j}}}\right)}^2}\\ +{\beta }_2{\displaystyle \sum _{1\le i<j\le n}{s}_{2_{ij}}{\left(\frac{d_i}{\sqrt{l_{H{2}_i}}}-\frac{d_j}{\sqrt{l_{H{2}_j}}}\right)}^2}\\ +{\lambda }_1\parallel g{\parallel }_1+{\lambda }_2\parallel d{\parallel }_1\\ \text{s.t. }{g}^{T}g=1,d^{T}d=1.\end{array}$$ (14)

Here, $s_{1_{ij}}$ and $s_{2_{ij}}$ represent the i – th row and j – th column of the X and Y similarity matrices S₁ and S₂. $l_{H{1}_i}$ and $l_{H{2}_j}$ represent the i – th row and j – th column of $l_{H{1}_i}$ and $l_{H{2}_j}$, respectively. Similar to literature [13], this paper uses the coordinate descent algorithm to find the local maximum of this problem by alternately updating the variables g and d.

The objective of HR-SPLS is to discover a low-dimensional representation containing the most relevant information from both X and Y. Specifically, the algorithm achieves integration by identifying latent variables, denoted as g and d, between X and Y. These latent variables are obtained by projecting X and Y onto a new coordinate system, aiming to maximize their covariance in this new coordinate space $(-\frac{1}{p}{g}^T{X}^{T}Yd)$. This ensures the preservation of the most relevant information between X and Y in the latent variables. The introduction of regularization terms (l₁ norm and hypergraph regularization) mitigates overfitting caused by an excessive number of features, thereby enhancing the model’s generalization capability.

In addition, this paper uses the solution of PLS as the initial solution of the current algorithm. Specifically, any column of X and Y is first randomly selected as the initial values of u and v. The following formula is then used to iteratively select the objective function value for g, d, u and v.

$$g:=\frac{X^{T}v}{v^{T}v};g:=\frac{g}{\parallel g{\parallel }_2};u:=Xg$$ (15)

$$d:=\frac{Y^{T}u}{u^{T}u};d:=\frac{d}{\parallel d{\parallel }_2};v:=Yd$$ (16)

Next, the weight vectors g and d are updated alternately. First, fix d, and obtain the partial derivative of g to obtain the iterative update rule of g.

$$g_j\leftarrow \frac{\text{sign}\left(z_g\right){\left(\left|z_g\right|-{\lambda }_1\right)}_{+}}{2\left({\beta }_1+{\delta }_1\right)}j=1,2,\dots ,n({\delta }_1>0)$$ (17)

Among them, the $z_g=t_{g_j}+2{\lambda }_1{\sum }_{i=1}^n\frac{s_{ij}{g}_i}{\sqrt{l_{H{1}_i}{l}_{H{1}_j}}}$, $t_g=\frac{1}{p}\left(X^{T}Yd\right)=\frac{1}{p}\left(X^{T}v\right)$. $t_{g_j}$ is the j – th element of t_g. Finally, the iterative update rule of d can be obtained by fixing g and taking the partial derivative for d.

$$d_j\leftarrow \frac{\text{sign}\left(z_d\right){\left(\left|z_d\right|-{\lambda }_2\right)}_{+}}{2\left({\beta }_2+{\delta }_2\right)}j=1,2,\dots ,n({\delta }_2>0)$$ (18)

Among them, $z_d=t_{d_j}+2{\lambda }_2{\sum }_{i=1}^n\frac{s_{ij}{d}_i}{\sqrt{l_{H{2}_i}{l}_{H{2}_j}}}$, $t_d=\frac{1}{p}\left(Y^{T}Xg\right)=\frac{1}{p}\left(Y^{T}u\right)$ is the j – th element of t_d.

Module membership confirmation and correlation analysis

After obtaining the final weight vectors g and d, this paper uses the Z-score method to confirm whether features of X and Y are eligible to enter the module. Specifically, g and d are first normalized using the Z-score method to obtain g^* and d^*. Then, the corresponding features whose g^* and d^* are more significant than the artificially set threshold T are selected for subsequent analysis. In addition, this paper normalizes the l₂ norm of u and v and obtains the normalized u^* and v^*. Then normalize (u^*+v^*), and select the sample when (u^*+v^*)>T. The threshold set in this paper is T = 1. After running the algorithm to get the first module, this paper subtracts the module signal from the input data:

$$X:=X-u{p}^T,p=\frac{X^{T}u}{u^{T}u}$$ (19)

$$Y:=Y-v{q}^T,q=\frac{Y^{T}v}{v^{T}v}$$ (20)

To confirm the correlation of two kinds of data in the same module, this paper defines the Pearson correlation coefficient (PCC) as follows.

$$PCC\left(X,Y\right)=\frac{\text{cov}\left(X,Y\right)}{{\sigma }_X{\sigma }_Y}$$ (21)

Among them, σ_X and σ_Y points represent the standard deviation of X and Y. This paper calculates the Pearson correlation coefficient PCC (u, v) between Xg and Yd within each module as a measure of module selection. In addition, this paper also introduces the module error to measure the module’s performance, which is defined as follows.

$$\begin{array}{l}ModuLe\_Error=\text{p}/n\ast {\sum }_{ij}{\left(X-u{p}^T\right)}_{ij}^2\\ +\text{p}/n\ast {\sum }_{ij}{\left(Y-v{q}^T\right)}_{ij}^2\end{array}$$ (22)

Survival analysis

We obtained the overall survival (OS) and disease-free survival (DFS) time of PCa patients from clinical data. Kaplan-Meier (KM) analysis was used to screen out the methylation sites associated with OS in PCa patients. Methylation sites with a p-value less than 0.05 were used as the input for constructing a diagnostic model for PCa patients. Univariate Cox regression analysis was used to screen out candidate features associated with DFS in PCa patients. Genes with p-values less than 0.03 were reserved for further research. Then, we constructed PCa-related prognostic models using multivariate Cox regression analysis and calculated risk scores for PCa patients. risk score=(β_{mRNA_1} * expression level of mRNA_1)+ (β_{mRNA_2} * expression level of mRNA_2)+…+(β_{mRNA_n} * expression level of mRNA_n). Based on the median risk score, we divided PCa patients in the TCGA-PRAD cohort and GSE116918 into two risk subgroups. KM curves were used to compare the differences in DFS of patients in the two risk subgroups. We used Receiver Operating Characteristic (ROC) curves to assess the accuracy of prognostic models for predicting 1-, 3-, and 5-year survival in PCa patients. ROC curves are drawn by the “timeROC” package.

Gene set enrichment analysis (GSEA)

We downloaded the “c2.all.v7.2.symbols.gmt” gene set from the GSEA database (https://www.gsea-msigdb.org/gsea/msigdb/index.jsp). Based on the “c2.all.v7.2.symbols.gmt” gene set, we performed GSEA analysis on high- and low-risk subgroups to identify pathways enriched between the two risk subgroups. The size of the gene set is set from 10 to 500. Gene sets were considered significant pathways when the absolute value of NES was greater than 1.5, p-value < 0.05, and FDR > 0.25.

Analysis of immune infiltration among risk subgroups

To further explore the immune microenvironment between the two risk subgroups, we performed a Cell-type Identification By Estimating Relative Subsets Of RNA Transcripts (CIBERSORT) analysis using the “e1071” package. CIBERSORT analysis was used to calculate the composition and infiltration levels of 22 immune cells (T cells, B cells, macrophages, dendritic cells, Natural Killer (NK) cells, monocytes, mast cells, eosinophils, and neutrophils) between the two risk subgroups. The Pearson correlation coefficient was used to calculate the correlation between prognosis-related genes and immune cells.

Data availability

The gene expression data of prostate cancer patients (TCGA-PRAD) were downloaded from the TGCA database (https://portal.gdc.cancer.gov/). The methylation data and corresponding clinical information of prostate cancer patients in the TCGA-PRAD cohort were downloaded from the UCSC Xena database (https://xenabrowser.net/datapages/). In addition, we downloaded the GSE116918 dataset from the GEO database (https://www.ncbi.nlm.nih.gov/geo/) for external validation of the prognostic model.

Data preprocessing

The Limma package was used for the differential expression analysis of gene expression data and methylation data on samples from the training set. Genes with absolute logFC values greater than 1 and FDR values less than 0.05 were considered differentially expressed. Methylated sites with FDR values less than 0.05 were considered differentially expressed. We obtained 1350 differentially expressed mRNAs (Figure 2A) and 1469 differentially expressed methylation sites (Figure 2B). Differentially expressed genes and expression profiles of methylation sites were used as input data for the algorithm.

Selection of hyperparameters

This paper selects 20 modules. Each hyperedge’s maximum number of vertices when building the hypergraph is chosen (Figure 3A). Our algorithm incorporates methylated sites and genes with significant associations into the same co-expression module. Consequently, the correlation between two modalities of data within the module can be assessed using Pearson Correlation Coefficient (PCC). A higher PCC indicates a stronger correlation among members within the module, providing additional confirmation of the algorithm’s capabilities in association analysis and feature selection. When the maximum number of vertices is 2, pcc(u, v) is the largest. Further, this paper also uses pcc(u, v) to select the four hyperparameters of λ₁, λ₂, β₁ and β₂ from the range of [0.01 0.05 0.1 0.5 1]. We present the PCC (Pearson Correlation Coefficient) of the algorithm for 625 parameter combinations in Figure 3B.

The maximum module correlation can be obtained when the number of neighbors is (Figure 3A). The largest module correlation can be obtained under the 157th set of parameter combinations (Figure 3B). The parameter value corresponding to the 20th group of parameters is λ₁ = λ₂ = β₁ = β₂ = 0.5.

Module description and selection

According to the parameter selection results in Section 3.1, 20 modules are obtained in this paper. To select the salient modules among them, we draw heatmaps of the weights of u and v corresponding to the 20 modules, respectively (Figure 4A, 4B). In addition, a line graph of module error for the 20 modules was also plotted (Figure 4C).

Among the 20 modules obtained by the proposed algorithm, the weights of u and v corresponding to module 2 are higher (Figure 4A, 4B). Furthermore, module 2 has a minor error (Figure 4C). Therefore, module 2 will be analyzed in detail later.

Comparison with other algorithms

To confirm the performance of the HR-SPLS algorithm, this paper introduces the SPLS algorithm and SNPLS algorithm to compare the performance of the three algorithms. The module error and membership correlation were used to compare the performance of the 20 modules corresponding to the three algorithms (Figure 5A, 5B).

Further, this paper counts the mean value of the module error and the mean value of the module correlation in the 20 modules of the three algorithms (Table 1). To confirm the importance of sparse constraints for feature selection in high-dimensional omics data, we compared the objective functions of HR-SPLS obtained by adding sparse constraints or not. Under the same number of iterations, the objective function values without and with sparse constraints were 4.3075 and 0.5667, respectively. Therefore, sparse constraints can significantly accelerate the convergence speed of the algorithm and make the performance of the algorithm better.

The module correlation of HR-SPLS is better than the other two algorithms (Table 1). The mean value of the module error is between the two different algorithms, which further confirms the correlation performance of the algorithm on the two kinds of data.

Additionally, we explored whether elastic net regularization could further enhance the performance of the algorithm proposed in this paper. Specifically, elastic net regularization is a method that combines L1 (Lasso) and L2 (Ridge) regularization. The objective function of applying this regularization to the algorithm in this paper is presented below.

$$\begin{array}{c}\underset{g,d}{mim}-\frac{1}{p}{g}^T{X}^{T}Yd+{\beta }_1\\ {\displaystyle \sum _{1\le i<j\le n}{s}_{1_{ij}}{\left(\frac{g_i}{\sqrt{l_{H{1}_i}}}-\frac{g_j}{\sqrt{l_{H{1}_j}}}\right)}^2}\\ +{\beta }_2{\displaystyle \sum _{1\le i<j\le n}{s}_{2_{ij}}{\left(\frac{d_i}{\sqrt{l_{H{2}_i}}}-\frac{d_j}{\sqrt{l_{H{2}_j}}}\right)}^2}\\ +{\lambda }_1\left|\right|g|{|}_1+{\lambda }_2|\left|d\right|{|}_1+{\lambda }_1\left|\right|g|{|}_2+{\lambda }_2|\left|d\right|{|}_2\\ \text{s.t.\hspace{0.17em}}{g}^{T}g=1,d^{T}d=1.\end{array}$$ (23)

On the basis of the optimal HR-SPLS algorithm, this study fine-tuned the parameters γ₁and γ₂within the range [0.01, 0.05, 0.1, 0.5, 1]. Under the optimal parameters (λ₁ = λ₂ = β₁ = β₂ = 0.5), the algorithm yielded the minimum module error for module 14. The mean of module errors, mean of module correlations, and objective function value were found to be 2.1395, 0.6810, and 28.9895, respectively.

Diagnostic model construction

First, KM analysis was performed on the 105 methylation sites in module 2 to screen out the methylation sites associated with OS in PCa patients. Next, based on the random forest (RF) algorithm, the feature weights were assigned to 25 methylation sites (Supplementary Figure 2) associated with the prognosis of PCa patients (Figure 6A). Further, this paper uses the logistic regression (LR) algorithm, the RF algorithm, and the K-Nearest Neighbor (KNN) algorithm to construct the diagnosis model of prostate cancer. Specifically, we used different numbers of Top features, put them into three classifiers, and compared the AUC of the classifiers (Figure 6B).

Use the characteristics of 10 top LR algorithm (cg20210585 cg12567282, cg11709110, cg13428921, cg24898914, cg11183227, cg22288195, cg24780796, cg24827036, cg11532 655, cg21769117 cg11254726, classify cg25653336), on the test set can reach the highest AUC for the 0.9378. To further validate the effectiveness of the algorithm, this paper introduces two non-negative matrix factorization (NMF) based algorithms, namely MDJNMF [12] and JDSNMF [14]. The ROC curves of the diagnostic models constructed by these two algorithms are presented in Supplementary Figure 1. The diagnostic model constructed by our algorithm achieved the highest AUC.

Construction of mRNAs-related prognostic model

We extracted expression data for 104 mRNAs in Module 2 and clinical information from PCa patients. First, a univariate Cox regression analysis was performed on the expression data of mRNAs in the TCGA-PRAD cohort. According to the p-value of less than 0.03, we screened and obtained 31 mRNAs related to DFA in PCa patients (Supplementary Table 1 and Supplementary Figure 3). Next, we performed a multivariate Cox regression analysis on 31 mRNAs to construct a prognostic model. Finally, we obtained a prognostic model (Supplementary Table 2) associated with 7-mRNAs (MYH11, ACTG2, DDR2, CDC42EP3, MARCKSL1, LMOD1 and MYLK), the risk score of the prognostic model is equal to expression level of MYH11* (-3.648) + expression level of ACTG2* (-4.820) + expression level of DDR2* (-2.740) + expression level of CDC42EP3* (-3.481) + expression level of MARCKSL1 *(0.845) + expression level of LMOD1* (7.614) + expression level of MYLK* (3.542).

We divided PCa patients in the TCGA-PRAD cohort and GSE116918 into high and low-risk subgroups based on the median risk score. The results of the KM analysis showed that patients in the high-risk group in the TCGA cohort (p=0.003) and the GEO (p=0.041) cohort had significantly shorter DFS (Figure 7A, 7B). We further assessed the prognostic model’s predictive accuracy using the ROC curve’s AUC area. The results showed that the constructed 7-mRNAs model could predict the 1-year (AUC=0.725), 3-year (AUC=0.702), and 5-year survival rates (AUC=0.702) of PCa patients in the TCGA cohort with high accuracy 0.761) (Figure 7C). In the external dataset, the AUCs at 1, 3, and 5 years were 0.927, 0.664, and 0.685 (Figure 7D). These results demonstrate that the risk scoring model, validated in the test set, can be used to predict DFS in PCa patients. We also created heatmaps of risk factors in the TCGA cohort (Figure 7E–7G) and the GEO cohort (Figure 7H–7J). The results showed that our risk score divided PCa patients into two risk subgroups, with high-risk patients having a shorter survival time than low-risk patients. MYH11 and ACTG2 were lowly expressed in the high-risk group, while MARCKSL1 was highly expressed in the high-risk group.

In addition, to explore enriched biological pathways between the two risk subgroups, we performed GSEA analysis on high and low-risk subgroups. The enrichment analysis showed that the high and low-risk groups were mainly enriched in immune and inflammation-related pathways (Figure 8A–8L), such as KEGG_NATURAL_KILLER_CELL_MEDIATED_CYTOTOXICITY, WP_B_CELL_RECEPTOR_SIGNALING_PATHWAY, KEGG_TOLL_LIKE_RECEPTOR_SIGNALING_ PATHWAY, and REACTOME_PI3K_AKT_SIGNALING_IN_CANCER.

Analysis of the immune microenvironment among risk subgroups

We performed an immune infiltration analysis of the TCGA-PRAD cohort using the CIBERSORT method to further explore the immune microenvironment between the two risk groups. First, we calculated the proportion of 22 types of immune cells in high- and low-risk patients (Figure 9A). Next, the correlation of 22 immune cells in the TCGA cohort was explored. Dendritic cells activated and T cells CD4 memory activated had the strongest positive correlation (cor=0.42), and T cells CD4 memory resting had a strong negative correlation with T cells CD8 (cor=-0.38) (Figure 9B). Subsequently, this paper compared the infiltration levels of 22 types of immune cells in high and low-risk groups (Figure 9C). The results showed that a variety of immune cells were different between high and low-risk groups, including T cells CD4 memory resting, T cells follicular helper, T cells regulatory (Tregs), NK cells activated, Monocytes, Macrophages M0, Mast cells activated, Eosinophils, Neutrophil, and Dendritic cells resting. Among them, the infiltration level of T cells CD4 memory resting, Monocytes, Dendritic cells resting, and Mast cells activated in the high-risk group was lower in the low-risk group. In addition, this paper compared the differential expression of 47 immune checkpoints between high and low-risk groups, and a total of 33 immune checkpoints were differentially expressed (Figure 9D). Finally, the correlation between seven prognosis-related mRNAs and immune cells was explored (Figure 10A–10Y and Supplementary Figure 4). ACTG2 was positively correlated with Mast cells resting (R=0.16, p=0.0026) but negatively correlated with Macrophages M1 (R=−0.24, p=7.7e−06). CDC42EP3 was positively correlated with T cells CD4 memory resting (R=0.23, p=1.6e−05). DDR2 was positively correlated with B cells naive (R=0.21, p=0.00012). LMOD1 was negatively correlated with Macrophages M1 (R =−0.22, p = 4.4e−05). MARCKSL1 was positively correlated with Macrophages M0 (R=0.19, p=0.00038). MYH11 was negatively correlated with Macrophages M1 (R=−0.17, p=0.0021). MYLK was positively correlated with T cells CD4 memory resting (R=0.23, p=2.4e−05). These results suggest that these immune cells play a crucial role in tumor progression. In addition, we calculated the correlations of T-stage, N-stage, and Gleason score with relation to immune features (“gleason_score.cor_immune “folder in the Supplementary Material). The outcomes revealed a notable correlation between the Gleason score and Macrophages M2, Plasma cells, and T cell regulation (Tregs). Subsequently, we employed the Wilcoxon method to determine differences in distinct immune cells (retaining those with higher immune abundance) across different T and N stages. The findings indicated significant differences in T cell regulation (Tregs) during the N-stages and notable variances in Plasma cells during the T-stages. These results suggest a pivotal role for T cell regulation (Tregs) and Plasma cells in the immune microenvironment of Pca.