In our study, we explored the m5C regulatory gene expression profiles between PCa and normal samples based on the TCGA database and GEO database. Then, the univariable Cox regression and LASSO Cox regression analyses were applied to identify prognostic biomarkers and develop a prognostic model. Functional enrichment analysis, clinical characteristics analysis, and immune infiltration analysis were conducted. Furthermore, we identified m5C regulatory gene clusters and m5C immune subtypes using a consensus-unsupervised clustering analysis and investigated the correlation between the subtypes and tumor immunity.

The gene expression data and clinical information of PCa patients were obtained from The Cancer Genome Atlas (TCGA) database (https://portal.gdc.cancer.gov/) and Gene Expression Omnibus database (https://www.ncbi.nlm.nih.gov/geo/) (Supplementary Table 1) (25, 26). A total of 499 PCa cases and 52 normal cases with the gene expression profile (HTSeq-FPKM), copy number variation (CNV), single-nucleotide polymorphism (SNP), and relevant clinicopathological information of prostate adenocarcinoma projects were collected from the TCGA database. Level 3 HTSeq-FPKM data were transformed into transcripts per million (TPM) reads for the subsequent analyses. Three GEO cohorts including GSE3325, GSE55945, and GSE155056 datasets were acquired from the GEO database (27–29). The gene expression data of GSE3325 (tumor=13, normal=6) and GSE55945 (tumor=11, normal=8) were based on the platform of GPL570 (Affymetrix Human Genome U13s3 Plus 2.0 Array) and GSE155056 (tumor=3, normal=3) was based on the platform of GPL28784 (085499_SBC human ceRNA microarray version 1.1). The samples in 3 GEO cohorts were all collected from human prostate benign and malignant tissues.

m5C regulatory genes (TET1, TET3, DNMT3B, YBX1, NSUN2, NSUN6, NOP2) were extracted from the prior studies (18, 30). First, three GEO datasets were integrated and the “Combat” algorithm from “SVA” package of R software was employed to eliminate the batch effects caused by non-biotechnological bias (31). Then, we used the “Limma” package to identify differentially expressed genes (DEGs) between the tumor and the normal samples. |Log2(Foldchange)| >1 and an adjusted P-value <0.05 was regarded as the threshold to indicate DEGs. Based on all m5C regulatory gene expression levels, principal component analysis (PCA) was performed using the “prcomp” function. Results were visualized using the R package “ggplot2”.

The Genomic Identification of Significant Targets in Cancer (GISTIC) algorithm is widely applied to identify both broad and focal (potentially overlapping) recurring events, which are associated with trigger cancer pathogenesis (32). GISTIC 2.0 software was employed to detect genes showing significant deletion and amplification in thousands of samples. An amplification or deletion length >0.1 and P<0.05 were considered as the parameter threshold. We utilized MutSig2.0 software to analyze the Mutation Annotation Format (MAF) file of TCGA mutation data to identify significantly mutated genes. A P-value <0.05 was considered significant.

We applied m5C regulatory genes to construct a prognostic risk score signature. Firstly, the univariate Cox regression analysis was conducted to extract m5C regulatory genes that were associated with the overall survival significantly. The result was visualized by the R package “forestplot”. Then, the least absolute shrinkage and selection operator (LASSO) Cox regression analysis was performed to reduce the dimension of high-latitude data and construct the prognostic model using the R package “glmnet” (33, 34). Ten-fold cross-validation was employed to avoid the overfitting problem and select the penalty parameter (λ) according to the minimum criteria. The prognostic scoring system for PCa patients was developed based on a linear combination of regression coefficients derived from the LASSO Cox regression analysis coefficients multiplied by the expression levels of genes, and then patients were divided into high-risk and low-risk subgroups accordingly. We compared the overall survival probability of high-risk and low-risk subgroups using the R package “survival” and “jskm”. A landmark time of 6 years was set. To evaluate the stability of the prognostic model, the R package “survivalROC” was employed to perform receiver operating characteristic (ROC) analysis and calculate the value of the area under the curve (AUC). The expression levels of m5C regulatory genes were analyzed between different risk subtypes classified by the prognostic model.

Gene set variation analysis (GSVA) is a non-parametric, unsupervised method for estimating the variation of gene set enrichments in the gene expression data, which is commonly used for exploring the variation in the pathway and biological process activity in the samples of an expression dataset (35). We performed GSVA to investigate the difference of the biological function between high- and low-risk subgroups. “c2.kegg.v7.4.symbols” and “c5.go.v7.4.symbols” gene sets were applied to performed GSVA using the R package “GSVA”. The R package “pheatmap” was applied to visualize the results.

Gene set enrichment analysis (GSEA) was conducted by the R package “clusterProfiler” to determine whether prior-defined functional or pathway sets of genes differ significantly between high- and low-risk subgroups in the enrichment of MSigDB Collection (“c2.kegg.v7.4.symbols” and “c5.go.v7.4.symbols” gene sets) (http://www.gsea-msigdb.org/gsea/index.jsp) (36). The enrichments of gene sets with a P-value <0.05 were regarded to be significant.

The R package “ESTIMATE” can calculate the TME scores including the stromal score, immune score, and estimate score using gene expression datasets and evaluate the relative contents of stromal cells or immunocytes in the TME. Higher stromal scores or immune scores suggest higher relative contents of stromal cells or immunocytes. Estimate scores represent the aggregation of stromal scores and immune scores (37). We investigate the correlation between TME scores and high- or low-risk subgroups.

To individualize the predicted PCa patients’ survival probability for 1, 3, and 5 years, a nomogram was developed using the R package “regplot.” The nomogram contained the risk score of the m5C prognostic model and clinical characteristics including the age, M stage, and T stage. Then, we conducted calibration analysis and calculated the optimism-corrected concordance index (C-Index) by 1,000 bootstrap resamples to assess the discrimination of the predictive nomogram.

Based on the expression of m5C regulatory genes, we performed consensus unsupervised clustering analysis to classify PCa patients into distinct gene clusters employing R package “ConsensusClusterPlus”, with the parameters of reps = 1,000 and pItem = 0.8 (38). To verify the stability of classification, PCA was conducted based on the expression of all m5C regulatory genes and the R package “ggplot2” was used to visualize the results. Furthermore, we investigated the relationship of the m5C prognostic model, m5C regulatory gene clusters, and clinical pathological features.

Cell-type Identification by Estimating Relative Subsets of RNA Transcripts (CIBERSORT) is a de-convolution algorithm that uses a set of reference gene expression matrixes to evaluate 22 immune cell-type proportions from bulk tumor sample expression data based on support vector regression. We conducted CIBERSORT analysis using the R package “CIBERSORT” to investigate TME cell infiltration between different groups. The relative contents of immune cells were calculated in distinct risk score subtypes and m5C regulatory gene clusters by the CIBERSORT algorithm.

Based on the relative contents of immune cells, which were significantly different between the high- and low-risk subgroups classified by the m5C prognostic model, as well as between different m5C regulatory gene clusters, we performed consensus unsupervised clustering analysis to classify PCa patients into distinct immune subtypes employing the R package “ConsensusClusterPlus”, with the parameters of reps = 50 and pItem = 0.8 (38). Meanwhile, we explore the expression levels of m5C regulatory genes between the different m5C immune subtypes.

The cytotoxic T lymphocyte–associated protein 4 (CTLA4) gene has been demonstrated to be a key immunotherapy related gene (39). To determine whether CTLA4 plays a role in immunotherapy for PCa patients through TME cell infiltration, we investigated the correlation of the expression of CTLA4 and m5C immune subtypes and infiltration levels for different immune cells.

All statistical analyses were conducted with R software (v4.0.2). To compare the variables between the 2 groups, we employed the independent sample t tests for normally distributed continuous variables and the Wilcoxon rank sum test (Mann–Whitney U tests) for nonnormally distributed continuous variables. All tests were 2 sided, and P < 0.05 was considered statistically significant.

Three GEO datasets were enrolled into one integrated dataset, and the batch effects were eliminated (Figures 1A, B). Figures 1C, D showed the result of PCA conducted in the integrated GEO dataset and TCGA dataset, respectively, which revealed significant differences in the m5C regulatory genes transcription profiles between PCa and normal samples. Seven m5C regulatory genes (TET1, TET3, DNMT3B, YBX1, NSUN2, NSUN6, NOP2) were extracted for subsequent analyses. Differential expression analysis, SNP analysis, and co-expression analysis were performed. CNV frequencies were computed based on the data from TCGA to identify genes with significant amplifications or deletions. Among m5C regulatory genes, NOP2 exhibited the highest amplification frequency and lowest deletion frequency (Figure 1E). Figure 1F exhibited the landscape of SNP and mutations in PCa samples. Figure 1G showed the m5C regulatory gene alterations closely associated with PCa progression, including frameshift insertion, frameshift deletion, missense mutation, synonymous mutation, nonsense mutation, and other non-synonymous or cleavage sites. NOP2, TET3, NSUN6, and DNMT3B were identified as differentially expressed genes (DEGs) (P < 0.05) (Figures 1H, I). Co-expression analysis showed that varying degrees of correlation existed among m5C regulatory genes (Figures 1J, K).

We conducted univariate Cox regression analysis to investigate the prognostic value of m5C regulatory genes and screened NSUN2 (HR=5.428, 95%CI=1.645-17.905, P=0.005), TET3 (HR=2.956, 95%CI=1.017-8.593, P=0.047), and YBX1 (HR=3.428, 95%CI=1.403-8.347, P=0.007) that were associated with the prognosis of PCa (Figure 2A; Supplementary Table 2).

M5C regulatory genes were fit in the LASSO Cox regression analysis, and three genes including NSUN2, TET3, and YBX1 were selected to develop an m5C prognostic model based on the optimal value of λ (Figures 2B, C). PCa patients were classified into high-risk and low-risk subtypes according to the risk scores calculated by the model. Supplementary Figure 1A showed the distribution of the risk score and cut-off value in the TCGA cohort. As the survival curves crossed, we used landmark survival analysis to compare the difference between the different risk subtypes (Supplementary Figure 1B). The result of landmark analysis survival showed a longer OS in PCa patients in the low-risk subtype within 6 years (P = 0.014). Nevertheless, no significant difference was found in the survival probability beyond 6 years (P = 0.578). Subsequently, the ROC curve was plotted to assess the accuracy of our model’s predictions. As shown in Figure 2D, the AUC of the model was 0.797, which suggested a good efficacy in prognostic prediction. The expression levels of m5C regulatory genes between high-risk and low-risk subtypes were investigated. DNMT3B, NOP2, NSUN2, NSUN6, TET3, and YBX1 were differentially expressed (P<0.001) (Figures 2E–J).

Then, to further explore the functional annotation between high-risk and low-risk subtypes, GSVA was performed (Supplementary Table 3). The results of the GSVA of gene ontology biological processes (GOBPs) showed that the high-risk subtype was significantly enriched in the regulation of protein localization to the chromosome telomeric region, positive regulation of telomerase RNA localization to the cajal body, telomerase RNA localization, IMP biosynthetic process, and protein localization to nucleoplasm (Figure 3A). As for KEGG terms in GSVA, the high-risk subtype was enriched in the nucleobase biosynthetic process, aminoacyl tRNA biosynthesis, spliceosome, mismatch repair, glyoxylate and dicarboxylate metabolism, and RNA degradation (Figure 3B). Gene set enrichment analysis (GSEA) was employed to identify the biological processes and signaling pathways involved in PCa between high- and low-risk subtypes. The most significantly enriched biological processes and signaling pathways are shown in Figures 3C–E and Supplementary Table 4. The risk score based on the m5C prognostic model was enriched in GOBP terms including the G protein coupled purinergic nucleotide receptor signaling pathway, G protein–coupled purinergic nucleotide receptor activity, sialic acid binding, T-cell receptor complex, and negative regulation of interleukin 8 production and KEGG terms including the B-cell receptor signaling pathway, calcium signaling pathway, cell adhesion molecules cams, chemokine signaling pathway, and cytokine–cytokine receptor interaction.

Furthermore, we evaluated the TME scores of high-risk and low-risk subtypes using the Estimation of STromal and Immune cells in MAlignant Tumor tissues using Expression data (ESTIMATE) algorithm. The results demonstrated that the patients in the high-risk subtype have a lower tumor immune infiltration level than those in the low-risk subtype (Figure 4A). Correlation analysis was performed to analyze the relationship between clinical features and risk scores. We observed that the risk score of PCa patients increased with the T stage (Figure 4B). We constructed a nomogram containing the risk score and clinical characteristics including the age, T stage, and M stage (Figure 4C). The risk score, T stage, and M stage were significantly associated with the prognosis of patients.

To further explore the potential biological characteristics of m5C regulatory genes in PCa patients, a consensus clustering algorithm was employed to classify patients into two distinct modification patterns based on the expression of 7 m5C regulatory genes, including 458 cases in modification pattern 1 and 93 patients in the modification pattern 2, which were termed as m5C regulatory gene cluster 1 and 2 (Figures 5A–C). The PCA plots both demonstrated an obvious different distribution between two clusters (Figure 5D). Then, we used a Sankey diagram to visualize the relationship between gene clusters, risk score subtypes, and the survival status (Figure 5E). We observed that the patients in m5C gene cluster 1 corresponded with the high-risk subtype, which dominated most of the PCa patients with a death status. Compared to gene cluster 2, gene cluster 1 had a significantly higher risk score (P<2.2e-16, Figure 5F). In addition, the association between the expression levels of m5C regulatory genes and m5C regulatory gene clusters was investigated. Figures 5G–K showed that DNMT3B, NSUN2, NSUN6, TET1, and TET3 were differentially expressed with a statistical significance between gene cluster 1 and gene cluster 2 (P <0.001).

To analyze the difference of tumor immune infiltration between high-risk and low-risk subtypes, we used the CIBERSORT algorithm to estimate the proportions of 22 distinct immune cell phenotypes between different subtypes. The analysis result of 22 kinds of cells showed B-cell naïve (P<0.05), plasma cells (P<0.05), macrophages M0 (P<0.001), macrophages M1 (P<0.001), and macrophages M2 (P<0.05), mast cells resting (P<0.05) and neutrophils (P<0.01) were statistically different between high-risk and low-risk subtypes (Figures 6A–H).

Moreover, between m5C regulatory gene cluster 1 and m5C regulatory gene cluster 2, we identified 9 immune cell phenotypes that were statistically different, including naive B cells (P<0.01), memory B cells (P<0.05), plasma cells (P<0.01), CD8 T cells (P<0.001), resting CD4 T cells (P<0.001), Tregs (P<0.001), M1 macrophages (P<0.001), M2 macrophages (P<0.01), and resting dendritic cells (P<0.001) (Figures 7A–J). We extracted 4 overlapping immune cell phenotypes (naive B cells, CD8 T cells, M1 macrophages, M2 macrophages) for subsequent analyses.

In order to further explore the correlation between key immune cell phenotypes and m5C in PCa, the R package “ConsensusClusterPlus” was used once again to classify patients based on 4 overlapping immune cell phenotypes (naive B cells, CD8 T cells, M1 macrophages, M2 macrophages) that were observed to be significantly different. PCa patients were divided into two m5C immune subtypes (subtype A, n=385; subtype B, n=166) (Figures 8A–C). Then, we conducted differential expression analysis for the immune checkpoint gene CTLA4 and m5C regulatory genes; we found that the expression of CTLA4 (P<0.01), NSUN6 (P<0.01), TET1 (P<0.05), and TET3 (P<0.01) differed significantly between the m5C immune subtype A and B (Figures 8D–G).

To reveal a potential correlation between the infiltration of immune cell phenotypes and the efficacy of immunotherapy, we performed the CIBERSORT algorithm to assess the association between the immune checkpoint gene CTLA4 and the abundance of immune cell phenotypes. As shown in the scatter diagrams (Figure 9), CTLA4 was positively correlated with many types of innate and acquired immune cell types including memory B cells (R=0.16, P=000049), activated dendritic cells (R=0.15, P=000082), resting dendritic cells (R=0.14, P=000017), M1 macrophages (R=0.18, P=4.2e-05), CD4 T cells (R=0.24, P=6.8e-08), CD8 T cells (R=0.2, P=8.6e-06), T follicular helper cells (R=0.097, P=0.03), gamma T cells (R=0.2, P=8.4e-06), and Tregs (R=0.24, P=8.5e-08). CTLA4 was negatively correlated with M2 macrophages (R=-0.13, P=0.0051), resting mast cells (R=-0.45, P<2.2e-16), monocytes (R=-0.13, P=0.003), and plasma cells (R=-0.43, P<2.2e-16).