Constructing a lactylation-related genes prognostic model to effectively predict the disease-free survival and treatment responsiveness in prostate cancer based on machine learning

doi:10.21203/rs.3.rs-3599127/v1

Download PDF

Research Article

Constructing a lactylation-related genes prognostic model to effectively predict the disease-free survival and treatment responsiveness in prostate cancer based on machine learning

https://doi.org/10.21203/rs.3.rs-3599127/v1

This work is licensed under a CC BY 4.0 License

Version 1

posted

You are reading this latest preprint version

Background: More and more studies have revealed that protein lactylation is an important mechanism for lactate to fulfill its duties and participate in important biological processes, which can regulate gene expressions through histone lactation, thereby promoting tumor spread, metastasis and immunosuppression. However, protein lactylation has been poorly studied in prostate cancer.

Methods: This study aimed to identify potential novel lactylation biomarkers of prostate cancer by biomarker analysis and to explore immune cell infiltration and treatment responsiveness. By downloading mRNA-Seq data of TCGA prostate cancer data for differential analysis, we obtained the differential genes related to Lactylation in prostate cancer. Five machine learning algorithms were used to screen for lactylation-related key genes for prostate cancer. The five overlapping key genes screened by five machine learning algorithms were used to construct a survival prognostic model by lasso cox regression analysis. Further, the relationships between the model and related pathways, tumor mutation and immune cell subpopulations and drug sensitivity were explored. Then, two risk groups were established according to the risk score calculated by the five lactylation-related genes (LRGs). Subsequently, a nomogram scoring system was established to predict disease-free survival (DFS) of patients by combining clinicopathological features and lactylation-related risk scores. Besides, the mRNA expression levels of five genes were verified in prostate cancer cell lines by qPCR.

Results: We identified 5 key LRGs (ALDOA, DDX39A, H2AX, KIF2C, RACGAP1) and constructed the LRGs prognostic model. The AUC values for 1 -, 3 -, and 5-year DFS in the TCGA dataset were 0.762, 0.745, and 0.709. The risk score was a better predictor of disease-free survival than traditional clinicopathological features in prostate cancer. The prostate cancer patients in the high-risk group have a higher proportion of regulatory T cells and M2 macrophage, a higher tumor mutation burden, and a worse prognosis. The high-risk group had a lower IC50 for certain chemotherapeutic drugs, such as Docetaxel, Paclitaxel.

Conclusions: The lactylation-related genes prognostic model can effectively predict the disease-free survival and therapeutic responses in prostate patients.

lactylation

prostate cancer

prognostic

tumor immune environment

machine learning

Prostate cancer is one of the most common malignant tumors of the male genitourinary system. According to the latest cancer clinical data in the United States in 2022, the incidence of prostate cancer in male tumors ranks first, and the mortality rate of prostate cancer ranks second [1]. Although the overall survival of prostate cancer is longer than other cancers, the recurrence rate of localized prostate cancer after radical prostatectomy is still high, and the recurrence progresses and eventually causes death. How to identify patients with high risk of recurrence early, and early intervention treatment may prolong the survival of patients and improve the quality of life of patients. Therefore, developing an effective Prognostic Signature that is required to improve patient prognosis and guide patients to assess cancer risk, conduct, precision medicine treatment (such as individualized chemotherapy and immunotherapy), is an urgent issue.

Lactylation is a novel posttranslational modification first reported in 2019, which modulates histones through the addition of lactyl groups to their lysine residues and facilitates specific gene transcription [2]. Yang and colleagues observed extensive protein lactylation in HCC tissues by immunostaining, and further confirmed that the emulsification mainly occurred on histones and non-histone proteins by mass spectrometry [3]. Wan et al. reported that protein lactylation widely exists in various human normal tissues and cancer tissue in 2022 Jul [4]. Protein lactation is a crucial mechanism through which lactate performs its duties and is required in key biological processes such as glycolysis related to cell function (5), macrophage polarization (6), neurodevelopment (7), and regulation of tumor spread (8). It has been reported in the occurrence and progression of a variety of tumors, such as liver cancer, breast cancer, colorectal cancer, etc., all play an important role. In prostate cancer, histone lactylation affects cancer cell plasticity and leads to transcriptional surges of neuroendocrine genes [9]. Besides, HIF1α lactylation enhances KIAA1199 transcription to promote angiogenesis and vasculogenic mimicry [10]. Based on these inspiring findings, we determined that these lactylation-related genes (LRGs) were closely connected to prostate cancer cell malignancy and spreading. Therefore, systematic evaluation of the relationship between differentially expressed LRGs and the prognosis, immune microenvironment, and treatment response of PCa is still worth further exploration. Here, we collected the lactylation-related data published thus far and systematically analyzed the expression levels of these genes in different databases to develop a novel prognostic signature based on LRGs to systematically explore the relationships between the signature and clinicopathological characteristics and disease progression in patients with prostate cancer (PCa). In addition, we further investigated its correlation with the tumor microenvironment (TME), mutation profiles, and the patient’s response to immunotherapy and chemotherapy in PCa. A lactylation-related gene signature that effectively predicts both prognosis and treatment responsiveness in PCa has not been reported to date.

Data acquisition

Figure 1 displays the overall procedure of the study. The mRNA transcriptome profiles of 502 prostate cancer samples and 52 normal samples were downloaded from The Cancer Genome Atlas (TCGA) database, https://protal.gdc.cancer.gov/ (accessed on 7 may 2023). The corresponding clinical information for TCGA-PRAD download from Prostate Adenocarcinoma (TCGA, Firehose Legacy) and could be directly downloaded from the following: http://www.cbioportal.org/. When a patient would correspond to more than one sample, only samples with a sample sequence number (Vial) of A were retained, and samples with no disease-free survival were removed, resulting in 491tumour samples. External validation data were downloaded from Prostate Adenocarcinoma (MSK, Cancer Cell 2010), and 140 prostate cancer samples remained after prostate Adenocarcinoma samples without disease-free survival were deleted. A scale method-based normalization was performed with the gene expression profiles using the R package “limma” (v4.2.2). The lactylation-related genes were gathered from Gene card (lactylation was entered into the Gene card website to search for 92 genes, and 79 related genes with Relevance score > 0.6 were obtained) and from previously published studies [11, 12, 13]. Wan et al. reported a novel diagnostic ion-based strategy and identified various lactylation sites in human cell lines [4]. We got a total gene set of 401 lactylation-related genes (LRGs) after deleting duplicate genes, found in Supplementary Table S1.

Identification of differentially expressed and prognostic genes

The mRNA levels of lactation-modified genes were extracted from tumor samples and adjacent tissues samples in the TCGA prostate cancer data for difference analysis by R package “limma” [14]. The criteria for differentially expressed genes (DEGs) identification were a false discovery rate (FDR) < 0.05 and |log2 FC| ≥ 0.585. By applying these criteria, genes expressed at more than 1.5-fold levels in tumor tissues and adjacent tissues were screened out at a false discovery rate of less than 0.05. Univariate Cox analysis to assess the association between disease-free survival (DFS) status was performed with a threshold of p < 0.05, and 16 prognosis-related genes were screened out. The PCA and t-SNE analysis were applied using the R packages “Rtsne” (v0.16) and “ggplot2” (v3.4.0).

Unsupervised cluster analysis

Unsupervised cluster analysis was performed with expression of 16 prognostic genes in them for patients (TCGA, n = 491) using the “ConsensusClusterPlu” R package [15]. The proportion of items per sample was 0.8, and the proportion of features per sample was 1. Partition Around Medoids was used. The adjacency distance matrix was determined as (Pearson correlation coefficient). Default settings were used for other parameters.

Selection of key genes and construction of the prognostic model

The 491 prostate cancer patients in TCGA cohort were randomly divided into training cohort and test cohort at a ratio of 1:1. To check whether the TCGA training set, test cohort and MSKCC cohort meet the conditions for the Cox proportional hazards model, we performed Schoenfeld residual check on TCGA training cohort, TCGA test cohort and MSKCC cohort to check whether the proportional hazards' assumption was met by R package "survminer" and "survival". Five different machine learning algorithms: the least absolute shrinkage and selection operator (LASSO), random forest, the optimal extreme gradient boosting (XGBoost), Boruta algorithm and Support Vector Machine Recursive Feature Elimination (SVM-RFE) were used for the eigengene screening. LASSO was implemented as a dimensionality reduction method to perform variable screening and complexity adjustment while fitting a generalized linear model. LASSO analysis was implemented with a penalty parameter utilizing a 10-fold cross-verification via the “glmnet” package. The XGBoost uses a gradient boosting framework and makes improvement on objective optimization function which is to optimize the loss function and complexity punishment by the R package “xgboost”.The Boruta algorithm was used to reduce the dimension of molecules according to clusters and to distinguish important features from irrelevant features (R-package Boruta, + 1 to the FPKM matrix value, and then take the log2, finally the function scale was used to standardize the matrix, doTrace = 2, ntree = 500).The SVM-RFE is a supervised machine learning method for support vectors that can find the best variables by halving the features each round when there are many features by the R package “e1071”.LASSO-penalized Cox regression analysis was performed on prognostic lactylation-related genes to find the potential gene set suitable for the prognostic signature, and the optimal penalty parameter λ and gene coefficients of the risk scoring formula were obtained [16]. The risk score model trained from the TCGA data was constructed as follows:

$$\text{R}\text{i}\text{s}\text{k}\text{s}\text{c}\text{o}\text{r}\text{e}={\sum }_{i=1}^{n}\left(\text{e}\text{x}\text{p}\times \text{c}\text{o}\text{e}\text{f}\right)$$

where N is the number of model genes; exp represents the expression value of genes; coef is the coefficient of each gene. 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 groups according to the median risk value. Kaplan-Meier analysis was conducted to evaluate differences in disease free survival (DFS) time between high-risk and low-risk groups. To evaluate the stability and specificity 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) [17]. The TCGA test cohort and MSKCC cohort were used to validate the prognostic model.

Gene set variation analysis (GSVA)

“c2. kegg.v7.4.symbols” and “c5.go.v7.4.symbols” gene sets were applied to performed GSVA to investigate the difference of the biological function between high- and low-risk groups in TCGA cohort by the R package “GSVA”.

Analysis of tumor mutation burden and immune cell infiltration

R package "Maftools" was used to analyze and visualize the somatic variation data (MAF files) of high and low-risk group prostate cancer samples in TCGA. We utilized CIBERSORT, ESTIMATE, MCP counter, TIMER algorithms, Immunophenoscore algorithms (IPS), EPIC algorithm and xCell algorithms to evaluated the immune cells levels of the two groups. Furthermore, we applied ssGSEA algorithm to quantify the subgroups of the infiltrating immune cells between the two groups. The differences in immune response under different algorithms were reveled using a Heatmap.

Drug sensitivity prediction of the risk model

We used the R package “OncoPredict” [18]to predict the half-maximal inhibitory concentration (IC50) of 198 drugs based on the Genomics of Drug Sensitivity in Cancer (GDSC) databases in vivo drug responses between high- and low-risk groups.

Construction and evaluation of the nomogram

With the “rms” and “survival” packages in R, a nomogram for predicting the 1-, 3-, and 5-year disease-free survival of PCa was constructed using the risk model with clinicopathological parameters such as T-stage, risk score, Gleason score. To evaluate the nomogram’s effectiveness, the Calibration curve and time-dependent ROC curves analysis was performed.

Cell Lines and Cell Culture

Human prostate cancer cells (PC3, C4-2, 22Rv1, DU145) were purchased from Cell Bank of Type Culture Collection of Chinese Academy of Sciences, Shanghai Institute of Cell Biology, Chinese Academy of Sciences. All PC cell lines were cultured in RP1640 medium (RP1640, Gibco) supplemented with 10% fetal bovine serum (FBS, Gibco), streptomycin at 37°C in 5% CO2.

RNA Extraction and Quantitative Analysis

Total RNA was extracted from cell lines with TRIzol Reagent (Invitrogen, USA). Total RNA was reverse-transcribed into cDNA with PrimeScript RT Master Mix (Takara, USA) and then used to perform quantitative real-time PCR (qRT–PCR) with SYBR qPCR Master Mix (Vazyme, China). GAPDH was used as an internal control for gene quantification. The 2 − ΔCT was calculated for every sample and normalized to GAPDH. The primer sequences used are shown in Supplementary Table S2.

Statistical analysis

All statistical analyses were conducted based on R (v4.2.3). Statistical significance was defined as p-value < 0.05 when there is no special description for above methods. Statistical significance is indicated with asterisks (*). A two-sided p value of < 0.05 was considered as statistically significant (* p < 0.05, ** p < 0.01, *** p < 0.001).

Identifying prognostic lactylation-related genes and subtypes in PCa

The clinical and pathological characteristics of PCa patients in the TCGA and MSKK cohorts are listed in Table 1. Using the R package “limma” with an absolute log2-fold change (FC) > 0.585 and an adjusted P value < 0.05 to perform differential expression analysis, 81 lactylation-related genes were differentially expressed between tumor tissues and adjacent nontumorous tissues in TCGA cohort (Fig. 2A). By univariate cox analysis of the relationship between 81 differentially expressed genes and disease-free survival(P < 0.05) of prostate cancer from TCGA, we obtained 16 of the 81 genes were significantly related to the prognosis of PCa (Fig. 2B). Close correlation was observed among the 16 lactylation-related DEGs (Fig. 2C). CNV status analysis showed a frequent alteration in 16 lactylation-related DEGs (Fig. 2D). It was noted that ALDOA only had copy number deletion and most alterations were losses in copy number. Unsupervised consensus clustering analysis of 16 prognostic lactylation-related genes could obtain two clusters with two different lactylation signature (Fig. 2E, F). In this study, consensus clustering based on differentially expression of lactylation-related genes was achieved using the R package “Consensus ClusterPlus”. The optimal clustering value was k = 2. Kaplan–Meier method suggested that, compared with the DFS of prostate cancer patients in C1 and C2 lactylation cluster, PCa patients in C2 lactylation cluster had significantly shorter DFS time (Fig. 2G; P = 0.045, HR = 1.53, 95% CI = 1.01–2.34). PCA analysis and t-Distributed Stochastic Neighbor Embedding suggested significant differences between C1 and C2 lactylation cluster(p = 0.001) (Fig. 2H, I). The heatmap for the association between 16 lactylation-related DEGs and clinicopathological manifestations was also analyzed (Fig. 2J).

Selection of key LRGs by machine learning and construction of a LRGs prognostic model with good performance

Prostate cancer patients with DFS in TCGA were randomly divided into training cohort (n = 245) and test cohort (n = 246) at a ratio of 1:1. Because we need to construct the Cox proportional hazards regression model later, the main premise for constructing the cox proportional hazards regression model is the assumed hazard ratio. We check each collaborators variable (T stage, N stage, metastasis status, PSA, Residual tumor, risk score) is in line with the cox proportional hazards model assumptions with Schoenfeld residuals. Each schoenfeld individual test p was not statistically significant (p > 0.05), the global schoenfelds were not statistically significant (Supplementary Fig. 1). Therefore, we can assume that this cox model meets the proportional hazards' assumption. Then, we used LASSO regression analysis, random forest, XGBoost, Boruta algorithm and SVM-RFE to find key genes associated with DFS in the TCGA training cohort. Five lactylation-related genes, ALDOA, DDX39A, H2AX, KIF2C and RACGAP1, were identified as key genes by the LASSO regression algorithm. The five variables were obtained based on the optimal value of λ = 0.069(Fig. 3A), and the lowest partial likelihood deviance is shown in Fig. 3B. According to the correlation map between the number of decision trees and the model error, we chose ntree was 600, the error rate of the model tended to stabilize and the importance of variables was ranked according to the VIMP method (Fig. 3C). Six key genes, ALDOA, DDX39A, H2AX, KIF2C, RACGAP1 and H2BC4, were selected by the minimum depth variable selection function in the R package “randomForestSRC”. SVM-RFE algorithm identified 8 key genes significantly associated with the DFS (Fig. 3D). The Boruta algorithm identified 6 key genes significantly associated with the DFS (Fig. 3E). XGBoost algorithm screened 10 most important signature genes affecting the disease-free survival of prostate cancer (Fig. 3F).The intersection of the RF, LASSO, SVM-RFE, Boruta and XGBoost results was shown in a Venn diagram in Fig. 3G. We identified 5 overlapping key genes, including ALDOA, DDX39A, H2AX, KIF2C, RACGAP1. The overlapping key genes screened by the five machine learning algorithms were consistent with those screened by the lasso algorithm. We wanted to use the five signature genes to construct a survival prognosis model. Then, the five key genes were used to construct a survival prognostic model by Lasso cox regression analysis in TCGA training set, and validated in TCGA test set and external validation set MSKCC. The risk score model was as follows: Risk score = (0.0880 × expression level of ALDOA) + (0.0108 × expression level of DDX39A) + (0.0199 × expression level of H2AX) + (0.2848 × expression level of KIF2C) + (0.0181× expression level of RACGAP1). Univariate and multivariate Cox regression analyses demonstrated that T stage, Gleason score, and risk score based on the signature of five lactylation-related genes were independent predictors of prognosis in patients with PCa (Fig. 3H、I).According to the median value of the risk scores in TCGA training cohort, the patients were divided into two groups, either high-risk group or low-risk group. KM analysis indicated that significantly poorer DFS in TCGA training cohort was detected among patients in the high-risk group compared to patients in the Low-risk group (log-rank test, p < 0.0001, HR = 6.04,95%CI:2.83–12.9,5-yearsDFS:0.49,95%CI:0.38–0.63) (Fig. 3J). A receiver operating characteristic (ROC) curve was constructed to estimate the model and assess the reliability of the risk score, and the areas under the curve (AUCs) in TCGA training cohort at 1 year, 3 years and 5 years DFS were 0.82, 0.797 and 0.769, suggesting that the risk model can be useful in predicting prognosis (Fig. 3K).The ROC curve analysis shows the highly sensitive and specific predictive performance of the risk score in TCGA training cohort (Fig. 3L).So, the Prognostic model had a better predictive performance in the TCGA training cohort.

An internal validation set (TCGA test cohort) and an external validation set (MSKCC cohort) were used to verify the prediction performance of the prognostic model

To verify the prediction model had a powerfully predictive performance in the other cohort, the TCGA test cohort and the MSKCC cohort were used to verify the predictive performance. According to the median value of the risk scores in TCGA test cohort, the patients were divided into high-risk group and low-risk group. KM analysis indicated that significantly poorer DFS in TCGA test cohort was detected among patients in the high-risk group compared to patients in the Low-risk group (log-rank test, p = 0.0003, HR = 3.31,95%CI:1.16–6.61,5-years DFS:0.64,95%CI:0.52–0.77) (Fig. 4A). In the TCGA test cohort, the AUC for the 1-year, 3-year, and 5-years DFS were 0.704,0.674 and0.644(Fig. 4B). The ROC curve analysis shows the highly sensitive and specific predictive performance of the risk score in TCGA test cohort (Fig. 4C). To better assess its predictive accuracy, the above results were replicated in the entire TCGA cohort. We found that the DFS of the high-risk group was shorter than that of the low-risk group (P < 0.001) and the area under the ROC curves (AUC) of the prognostic model for DFS was 0.762 at 1 year, 0.745, at 3 years, and 0.709 at 5 years in entire TCGA cohort (Fig. 4D,E). Also, in the ROC curve containing risk score, stage, Gleason score, residual tumor and prostate-specific antigen (PSA), the AUC value of the risk score was higher than other indicators in the entire TCGA cohort (Fig. 4F) Similarly, patients with PCa in MSKCC cohort were classified into a high-risk group and a low-risk groups according to the median risk score. KM analysis indicated that the DFS was shorter in patients with high-risk scores than those with low-risk scores in MSKCC cohort (log-rank test, p = 0.0021, HR = 2.91,95%CI:1.43–5.92,5-years DFS:0.62,95%CI:0.50–0.76) (Fig. 4G). The areas under the curve (AUCs) of the ROC curve at 1 year, 3 years and 5 years DFS were 0.835, 0.701 and 0.661, suggesting good predictive ability (Fig. 4H). The risk score was statistically different between TCGA subgroups C1 and C2(Fig. 4I). The risk score of each patient in the entire TCGA cohort and MSKCC cohort was calculated based on the expression levels and regression coefficients of these five genes. The distribution of risk scores and survival status in the entire TCGA cohort and MSKCC cohort is shown in Fig. 4J and Fig. 4K. The risk curves reflect the relationship between the risk score and disease-free survival status of patients with prostate cancer, and we found that the probability of recurrence or progression was higher in the high-risk patients than the low-risk patients.

Establishment and evaluation of the predictive nomogram

The independent predictors, including T stage, risk score and Gleason score, which affect the DFS of PCa patients, were incorporated into the nomogram model (Fig. 5A). Time-dependent C-index curves of different variables based on TCGA cohorts show the optimum performance of nomogram compared with other single factors (Fig. 5B). The nomogram AUCs in the TCGA-PRAD cohort for the 1-year, 3-year, and 5-year OS probability were 0.796, 0.755, and 0.746, respectively (Fig. 5C). To evaluate the accuracy of this nomogram, the ROC curves were utilized to compare this signature with several available clinical traits and showed that the predictive value of this nomogram model was more optimal compared to several clinical traits. (Fig. 5D). Moreover, the calibration curve of the nomogram model showed the actual DFS was close to the predicted DFS (Fig. 5E). DCA showed that the nomogram provided net benefit over different threshold probability ranges, respectively, suggesting that might have potential applications in different contexts (Fig. 5F).

Analysis of tumor immune microenvironment of high- and low-risk groups

First, we estimated the composition of infiltrating immune cells in high- and low-risk groups of TCGA cohort by CIBERSORT, ESTIMATE, MCP counter, TIMER algorithms, IPS algorithms, EPIC algorithms and xCell algorithms. The results of Immune cell Infiltration in high and low-risk groups were presented as a heatmap (Fig. 6A). From the heatmap, we have known that regulatory T cell and M2 Macrophages were significantly higher in the high-risk group compared to the low-risk group, while Neutrophils and Fibroblasts were significantly lower in the high-risk group. The increase of regulatory T cells and M2 macrophages are conducive to tumor immune escape and tumor recurrence, which suggested that the infiltration of these immune cell subtypes into the tumor microenvironment might have a significant impact on the prognosis of PCa patients. Based on risk grouping, in terms of immune cells, we discovered that the content of natural killer T cell, natural killer cell, immature dendritic and neutrophil was significantly lower in patients with high-risk group compared to the patients with low-risk group (Fig. 6B). In terms of immune functions, the Immune functions of Type II IFN response, CCR and MHC class I level were significantly downregulated in the patients with high-risk group compared to the patients with low-risk group with ssGSVA algorithm analysis (Fig. 6C). Then, we explored the correlation between infiltrating immune cells and risk group and prognostic genes with the CIBERSORT algorithm. We found that the risk score of LRGs was significantly correlated with regulatory T cells, M2 Macrophages, T cells CD4 memory resting and T cells follicular helper (Fig. 7A). Moreover, the risk score had a significantly positive correlation with the regulatory T cells, M2 macrophages and follicular helper T cells, and it had a significantly negative correlation with the T cells CD4 memory resting (Fig. 7B-E). Considering the importance of checkpoint inhibitors in clinical treatment, we further investigated potential changes in immune checkpoint expression between the two groups. Combined with the bar chart showing immune checkpoint expression, where most canonical markers of exhausted T including TIM3/HAVCR2, LAG3 and CTLA4 are highly expressed in high-risk group of patients, we could conclude that the abundance of exhausted T cells infiltration in high-risk group of patients may be higher than that in low-risk group [19]. Besides, we found that the expression of B2M, JAK1, JAK2, CD40, TNFSF4, CD28 in the high-risk group were significantly lower than those in the low-risk group (Fig. 7F).

Functional enrichment analysis and drug sensitivity

Exploring the functional annotation between high-risk and low-risk subtypes, we found that the high-risk group was significantly enriched in the lamin filament, U2AF complex, positive regulation of apoptotic DNA fragmentation and positive regulation of DNA catabolic process in the GSVA of gene ontology biological processes (GOBPs) (Fig. 8A). KEGG pathway analyses were performed to explore the biological processes and pathways using the gene set variation analysis (GSVA). We found that the high-risk group was enriched in the base excision repair, DNA replication, homologous recombination, pyrimidine metabolism, mismatch repair, cell cycle (Fig. 8B). The risk score was closely related to the DNA synthesis, repairing and degradation and cell cycle progression of PCa according to the above functional enrichment analysis, and we further analyzed the response of PCa patients in the TCGA cohort to DNA synthesis and repairing and cell cycle-related chemotherapy drugs (5-Fluorouracil, Cisplatin, Docetaxel, Paclitaxel, Gemcitabine, Epirubicin,Talazoparib). The results revealed that these drugs had lower half maximal (50%) inhibitory concentration (IC50) in patients of the high-risk group, implying that these patients may be more sensitive to these drugs (Fig. 8C-I). Based on the immune checkpoint analysis, we know that JAK1 and JAK2 are lower in the high-risk group, which is consistent with our results from the drug prediction analysis that we know that the high-risk group has a worse response to JAK inhibitors (JAK_8517, Ruxolitinib) (Fig. 8J-K).

Relationship between the LRGs signature and clinicopathological traits and mutation landscape

First, we evaluated the association between clinicopathological parameters and the risk score to further investigate the clinical relevance of the prognostic model. The risk score increased significantly with increasing Gleason score, T-stage and N-stage, but there was no statistical difference between the no metastasis (M0) and metastasis (M1) (Fig. 9A). In assessing the relationship between expression profiles of the signature genes in PCa tissues and related clinicopathological characteristics based on TCGA cohort, we further found the 5 key genes were differentially expressed between tumor tissues and normal tissues, and were highly expressed in tumor tissues (Fig. 9B). There were differences in expression of 5 key genes in different T stages and N stages, and the expression in advanced T stages (T3 + T4) and N1 was higher than that in early T stages(T1 + T2) and N0(Fig. 9C-D). Among the five key genes, except ALDOA, the expression was different in different Gleason scores and increased with the increase in Gleason score (Fig. 9E). We further investigated the mutation profiles separately in the high- and low-risk groups of PCa patients from TCGA cohort. Results indicated that Tumor Mutational Burden (TMB) was more widely distributed in the high-risk group, and miscellaneous mutation was the most common variant classification in PCa. The top five most frequent mutation genes in the high-risk group were TP53(20%), SPOP (16%), TITIN (12%), MUC16(10%), FOXA1(8%), while genes such as TITIN (13%), SPOP (7%), MUC16 (6%), MUC4(6%), SYNE1 (6%) had the top five mutation frequencies in low-risk group (Figure. 9F-G). Previous studies have shown that mutations in these genes are often associated with a poor prognosis for PCa [20–23].

The mRNA expression levels of lactylation-related genes in PCa cells

To further evaluate the mRNA expression levels of these five genes, we utilized their expression in the cell lines using RT-qPCR. The results indicated that most signature genes were expressed at a relatively high level in castration-resistant prostate cancer cell lines (C4-2, PC3,22Rv1, DU145) compared with non-castration-resistant prostate cancer cell line (LNCAP) (Fig. 10A-E). Castration-resistant prostate cancer cell lines are more malignant than non-castration-resistant prostate cancer cell line. In sum, the high expression of these five key genes is positively correlated with the malignant degree of prostate cancer, which validated the previous accuracy of our previous bioinformatics analysis.

Prostate cancer is an extremely common malignancy in man. The rates of PCa recurrence and metastasis remain high after radical prostatectomy and were associated with a poor prognosis. Despite the widely used classification system by TNM staging and Gleason score, the heterogeneous nature of PCa poses difficulty onto prognosis and therapeutic decision. Therefore, exploring new methods of PCa subtyping has become the most urgent task in bringing precision medicine into reality. Recently, new biomarkers have emerged for the prognosis of cancer patients, but this is the first time to discuss the establishment of a prognostic model related to Protein lactation in prostate cancer patients.

Lactylation modification was widely occurred on histones and non-histone proteins. Protein lactylation has been found to be related to cell metabolism and cancer immune regulation [24]. More and more evidence suggests that protein lactylation plays an indispensable role in the development and progression of cancer. In this study, the expression profiles of lactylation-related genes were mined from TCGA PRAD RNA-seq data and were initially analyzed followed by filtering out via differential analysis. Then, 81 LRDEs were subsequent to univariate cox regression analysis, and a total of 16 prognostic LRGs were identified. The five key LRGs were selected from the 16 prognostic genes by five machine learning algorithms (LASSO regression analysis, random forest, XGBoost, Boruta algorithm and SVM-RFE). The five key genes (ALDOA, DDX39A, H2AX, KIF2C, RACGAP1) were used to construct a survival prognosis model. Further, the five key LRGs were used to construct a survival prognostic model by Lasso cox regression analysis. The predictive performance of the model was verified using the KM curve and ROC curve and the results indicated that the model had good predictive performance. Meanwhile, the nomogram, which was constructed by combining risk score, T stage and Gleason score, had a good accuracy in predicting 1,3,5 years disease-free survival and can be used as a practical and reliable prognostic tool for PCa patients. Among the five key genes, all were significantly over expressed in prostate cancer tissues than tumor-adjacent tissue, and their expression was positively correlated with TNM stage. Many studies reported that these genes played a significant role in the pathogenesis of cancer. ALDOA, a glycolytic enzyme to catalyze the reversible conversion of fructose-1,6-bisphosphate to glyceraldehyde 3-phosphate and dihydroxyacetone phosphate, was reported to promote the glycolysis of cancer cells and may be significantly associated with the development, metastatic potential, and poor prognosis of various kinds of tumors [25–28]. Nakata et al. reported that DDX39A regulates androgen receptor splice variant AR-V7 generation [29]. H2AX coordinates DNA repair following genotoxic stress and double strand breaks [30]. The loss of H2AFX increased chromosomal instability seen in acute myeloid leukaemia, acute lymphoid leukaemia and Head and neck squamous cell may contribute to tumor development, progression and resistance to therapy in this cancer subtype [31–32]. Besides, H2AX promotes metastatic progression in breast cancer cells by preserving glycolysis via hexokinase-2[33]. Many studies reported that KIF2C enhanced hepatocellular cancer through the Ras/MAPK and PI3K/Akt signaling pathways and accelerated the growth of cervical cancer by blocking the stimulation of the p53 signaling pathway and activated mTORC1 pathway and promotes tumor cell motility and invasion [34–36]. A previous study indicated that RACGAP1 regulated the downstream factors of the PI3K/AKT signaling pathway and that the compensatory activation of the PI3K/AKT signaling pathway was closely associated with ADT drug resistance and neuroendocrine differentiation in PCa [37]. Besides, the results of qPCR verified that the most signature genes were transcribed at higher levels in more malignant castration-resistant prostate cancer cell lines.

Recently, immune-based treatment has emerged for patients with PCa, which has revolutionized cancer therapy and improved the survival of patients with many types of solid tumors. However, prostate cancer is recognized as a poorly immunogenic tissue with immunological ignorance showing low levels of antigen-presenting process and cytotoxic T-cell activation, high levels of immune checkpoint molecules and immunosuppressive cytokines/chemokines, and recruitment of immunosuppressive cells [38]. Due to the immunosuppressive microenvironment and heterogeneous nature of PCa, it is important to identify the molecular subtypes and characterize the TME that can predict response to immunotherapy and identify high-risk patients for early intervention. Therefore, we investigated the status of TME between the high-risk and the low-risk subtypes. With 7 immune algorithms, we found that the proportion of regulatory T cells and M2 macrophages were significantly increased in high-risk than low-risk group, while Myeloid dendritic cells was significantly decreased in the high-risk group. Previous studies have shown that lactylation modification can regulates inflammatory to reparatory macrophage transition and promotes immunosuppression of tumor-infiltrating myeloid cells [6, 39]. Sadasivan et al. found that that patients with high levels of infiltrating M2 macrophages had an almost 5-fold increased risk of PCa recurrence [40]. Our results obtained by immune algorithm were consistent with previous reports that regulatory T cells and M2 macrophages were high in prostate cancer in the high-risk group, so the prognosis of the high-risk group was worse, which may be related to its immunosuppression. Besides, patients in the high-risk group had significantly lower levels of type II IFN response, CCR and MHC class I compared to patients in the low-risk group by ssGSVA algorithm analysis. These findings strongly implied the potential roles of lactylation -related genes in reshaping the TME in PCa. The results of the GSEA analysis showed that the high-risk score group enriched in positive regulation of DNA catabolic process in gene ontology biological processes (GOBPs), and in base excision repair and DNA replication in KEGG. Then, we explored how patients responded to diverse treatments and tried to offer different recommendations based on LRGs risk groups. Using the R package “OncoPredict” to predict the IC50 of 198 drugs, we found that high-risk group patients were sensitive to DNA synthesis and repairing and cell cycle-related chemotherapy drugs, such as Docetaxel, Cisplatin, 5-Fluorouracil.

Our study systematically evaluated the expression of RLGs and their potential prognostic value in PCa. Moreover, a risk model of five RLGs was established in TCGA train cohort and validated using the TCGA test cohort and MSKCC cohort. We found that the high-risk group was correlated with elevated TMB and chemotherapy response, shorter DFS. Undeniably, there are still some limitations of our study. First, since research in the protein lactylation field is still in its preliminary stages, it is not yet clear how lactylation modification affects protein function in tumor cells has yet to be elucidated. Second, due to the lack of appropriate lactylation antibodies, the experimental validation of signature genes is currently hampered by practical problems. Third, the majority of the data are obtained from the public data collection.

We developed a novel lactylation-related prognostic model to predict the prognosis of prostate cancer patients, which had good efficacy in predicting the DFS of PCa patients and provided new insights into personalized therapies for PCa patients. Moreover, further research of these hub genes may contribute to molecular targeted therapy of prostate cancer.

DEGs differentially expressed genes;

DFS disease-free survival;

LRGs lactylation-related genes

AUC area under curve

ROC receiver operating characteristic

TME tumor microenvironment

PCa prostate cancer

ssGSVA single sample gene set variation analysis

CCR chemotactic cytokines receptors

RF random forest

LASSO the least absolute shrinkage and selection operator

XGBoost the optimal extreme gradient boosting

SVM-RFE Support Vector Machine Recursive Feature Elimination

Acknowledgment

The authors acknowledge the TCGA project, cBioPortal database and other groups for providing invaluable datasets for statistical analyses.

Ethics approval and consent to participate

Not applicable to this study.

Consent for publication

Not applicable to this research

Competing interests

We declare that there is no competing interests regarding the publication of this paper.

Funding

The authors declare no funding provided for this study.

Author Contributions

JP, JZ and JL conceived and designed the study. JP and JZ were responsible for data collection, collation and statistical analysis with bioinformatics methods. YC carried out data interpretation and chart drawing. JP and JL wrote the manuscript, which was further polished and confirmed by ZZ. All authors contributed to the article and approved the submitted version.

Availability of data and materials

The gene expression profile data of TCGA-PRAD were obtained from The Cancer Genome Atlas (TCGA, GDC (cancer.gov)). The corresponding clinical information for TCGA-PRAD downloaded from Prostate Adenocarcinoma (TCGA, Firehose Legacy) and could be directly downloaded from the following: http://www.cbioportal.org/. External validation data were downloaded from Prostate Adenocarcinoma (MSK, Cancer Cell 2010) and could be directly downloaded from the following: http://www.cbioportal.org/. All data in the research is available in publicly available databases. This study complies with its data use and publication rules.

Footnotes

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Jinyou Pan, Jianpeng Zhang and Jingwei Lin contributed equally to this work and share the first authorship.

Contributor information

Jinyou Pan, Email: [email protected]

Jianpeng Zhang, Email: [email protected]

Jingwei Lin, Email: [email protected]

Yinxin Cai, Email: cyx15625879613#@163.com

Zhigang Zhao, Email: [email protected]

Siegel RL, Miller KD, Fuchs HE, Jemal A. Cancer statistics, 2022. CA Cancer J Clin. 2022 Jan;72(1):7-33. doi: 10.3322/caac.21708. Epub 2022 Jan 12. PMID: 35020204.
Zhang D, Tang Z, Huang H, Zhou G, Cui C, Weng Y, Liu W, Kim S, Lee S, Perez-Neut M, et al. Metabolic regulation of gene expression by histone lactylation. Nature. 2019; 574:575–580.
Yang Z., Yan C., Ma J., Peng P., Ren X., Cai S., Shen X., Wu Y., Zhang S., Wang X., et al. Lactylome analysis suggests lactylation-dependent mechanisms of metabolic adaptation in hepatocellular carcinoma. Nat. Metab. 2023; 5:61–79.
Wan N, Wang N, Yu S, Zhang H, Tang S, Wang D, Lu W, Li H, Delafield DG, Kong Y, Wang X, Shao C, Lv L, Wang G, Tan R, Wang N, Hao H, Ye H. Cyclic immonium ion of lactyllysine reveals widespread lactylation in the human proteome. Nat Methods. 2022 Jul;19(7):854-864.
Li L, Chen K, Wang T, Wu Y, Xing G, Chen M, et al.. Glis1 facilitates induction of pluripotency via an epigenome-metabolome-epigenome signalling cascade. Nat Metab (2020) 2(9):882–92.
Irizarry-Caro R, McDaniel M, Overcast G, Jain V, Troutman T, Pasare C. TLR signaling adapter BCAP regulates inflammatory to reparatory macrophage transition by promoting histone lactylation. Proc Natl Acad Sci United States America. (2020) 117(48):30628–38.
Hagihara H, Shoji H, Otabi H, Toyoda A, Katoh K, Namihira M, et al.. Protein lactylation induced by neural excitation. Cell Rep (2021) 37(2):109820.
Yu J, Chai P, Xie M, Ge S, Ruan J, Fan X, et al.. Histone lactylation drives oncogenesis by facilitating mA reader protein YTHDF2 expression in ocular melanoma. Genome Biol (2021) 22(1):85.
He Y, Ji Z, Gong Y, Fan L, Xu P, Chen X, Miao J, Zhang K, Zhang W, Ma P, Zhao H, Cheng C, Wang D, Wang J, Jing N, Liu K, Zhang P, Dong B, Zhuang G, Fu Y, Xue W, Gao WQ, Zhu HH. Numb/Parkin-directed mitochondrial fitness governs cancer cell fate via metabolic regulation of histone lactylation. Cell Rep. 2023 Jan 30;42(2):112033.
Luo Y, Yang Z, Yu Y, Zhang P. HIF1α lactylation enhances KIAA1199 transcription to promote angiogenesis and vasculogenic mimicry in prostate cancer. Int J Biol Macromol. 2022 Dec 1;222(Pt B):2225-2243.
Zhang D, Tang Z, Huang H, Zhou G, Cui C, Weng Y, Liu W, Kim S, Lee S, Perez-Neut M, Ding J, Czyz D, Hu R, Ye Z, He M, Zheng YG, Shuman HA, Dai L, Ren B, Roeder RG, Becker L, Zhao Y. Metabolic regulation of gene expression by histone lactylation. Nature. 2019 Oct;574(7779):575-580.
Moreno-Yruela C., Zhang D., Wei W., Bæk M., Liu W., Gao J., Danková D., Nielsen A.L., Bolding J.E., Yang L., et al. Class I histone deacetylases (HDAC1-3) are histone lysine delactylases. Sci. Adv. 2022;8: eabi6696.
Cheng Z, Huang H, Li M, Liang X, Tan Y, Chen Y. Lactylation-Related Gene Signature Effectively Predicts Prognosis and Treatment Responsiveness in Hepatocellular Carcinoma. Pharmaceuticals (Basel). 2023 Apr 25;16(5):644.
M.E. Ritchie, B. Phipson, D. Wu, et al.Limma powers differential expression analyses for RNA-sequencing and microarray studies Nucleic Acids Res, 43 (7) (2015), p. e47
Wilkerson MD, Hayes DN. ConsensusClusterPlus: a class discovery tool with confidence assessments and item tracking. Bioinformatics. 2010 Jun 15;26(12):1572-3.
Friedman J, Hastie T, Tibshirani R. Regularization Paths for Generalized Linear Models via Coordinate Descent. J Stat Softw (2010) 33:1. Available at: /pmc/articles/PMC2929880/ [Accessed January 10, 2022]
Heagerty, P.J., Lumley, T., Pepe, M. S. (2000) Time-dependent ROC Curves for Censored Survival Data and a Diagnostic Marker Biometrics, 56, 337 – 344
Maeser D, Gruener RF, Huang RS. oncoPredict: an R package for predicting in vivo or cancer patient drug response and biomarkers from cell line screening data. Brief Bioinform. 2021 Nov 5;22(6):bbab260.
Wherry EJ, Kurachi M. Molecular and cellular insights into T cell exhaustion. Nat Rev Immunol. 2015;15(8):486-499. doi:10.1038/nri3862.
Nyquist MD, Corella A, Coleman I, De Sarkar N, Kaipainen A, Ha G, et al. Combined Tp53 and Rb1 Loss Promotes Prostate Cancer Resistance to a Spectrum of Therapeutics and Confers Vulnerability to Replication Stress. Cell Rep (2020) 31(8):107669.
Zhang P, An Z, Sun C, Xu Y, Zhang Z. Flg Gene Mutation Up-Regulates the Abnormal Tumor Immune Response and Promotes the Progression of Prostate Cancer. Curr Pharm Biotechnol (2022).
Shi Q, Zhu Y, Ma J, Chang K, Ding D, Bai Y, et al.. Prostate Cancer-Associated Spop Mutations Enhance Cancer Cell Survival and Docetaxel Resistance by Upregulating Caprin1-Dependent Stress Granule Assembly. Mol Cancer (2019) 18(1):170.
Jiang J., Huang D., Jiang Y., Hou J., Tian M., Li J., Sun L., Zhang Y., Zhang T., Li Z., et al. Lactate Modulates Cellular Metabolism Through Histone Lactylation-Mediated Gene Expression in Non-Small Cell Lung Cancer. Front. Oncol. 2021; 11:647559.
Chen L, Huang L, Gu Y, Cang W, Sun P, Xiang Y. Lactate-Lactylation Hands between Metabolic Reprogramming and Immunosuppression. Int J Mol Sci. 2022 Oct 8;23(19):11943.
Shi S, Ji S, Qin Y, et al. Metabolic tumor burden is associated with major oncogenomic alterations and serum tumor markers in patients with resected pancreatic cancer. Cancer Lett. 2015; 360:227–233.
Long F, Cai X, Luo W, Chen L, Li K. Role of aldolase A in osteosarcoma progression and metastasis: in vitro and in vivo evidence. Oncol Rep. 2014; 32:2031–2037.
Du S, Guan Z, Hao L, et al. Fructose-bisphosphate aldolase A is a potential metastasis-associated marker of lung squamous cell carcinoma and promotes lung cell tumorigenesis and migration. PLoS One. 2014;9: e85804
Chang L, Ni J, Beretov J, Wasinger VC, Hao J, Bucci J, Malouf D, Gillatt D, Graham PH, Li Y. Identification of protein biomarkers and signaling pathways associated with prostate cancer radioresistance using label-free LC-MS/MS proteomic approach. Sci Rep. 2017 Feb 22; 7:41834.
Nakata D, Nakao S, Nakayama K, Araki S, Nakayama Y, Aparicio S, Hara T, Nakanishi A. The RNA helicase DDX39B and its paralog DDX39A regulate androgen receptor splice variant AR-V7 generation. Biochem Biophys Res Commun. 2017 Jan 29;483(1):271-276.
Bonner WM, Redon CE, Dickey JS, Nakamura AJ, Sedelnikova OA, Solier S, Pommier Y. GammaH2AX and cancer. Nat Rev Cancer. 2008 Dec;8(12):957-67.
Thirman MJ, et al. Rearrangement of the MLL gene in acute lymphoblastic and acute myeloid leukemias with 11q23 chromosomal translocations. N Engl J Med. 1993; 329:909–914.
Parikh RA, et al. Loss of distal 11q is associated with DNA repair deficiency and reduced sensitivity to ionizing radiation in head and neck squamous cell carcinoma. Genes Chromosomes Cancer. 2007; 46:761–775.
Liu Y, Li H, Wilson CN, Bai HJ, Boufraqech M, Weyemi U. Histone H2AX promotes metastatic progression by preserving glycolysis via hexokinase-2. Sci Rep. 2022 Mar 8;12(1):3758.
Yang J, Wu Z, Yang L, Jeong JH, Zhu Y, Lu J, et al.. Characterization of Kinesin Family Member 2c as a Proto-Oncogene in Cervical Cancer. Front Pharmacol (2021) 12:785981.
An L, Zhang J, Feng D, Zhao Y, Ouyang W, Shi R, et al.. Kif2c Is a Novel Prognostic Biomarker and Correlated With Immune Infiltration in Endometrial Cancer. Stem Cells Int (2021) 2021:1434856.
Moon HH, Kreis NN, Friemel A, Roth S, Schulte D, Solbach C, et al.. Mitotic Centromere-Associated Kinesin (Mcak/Kif2c) Regulates Cell Migration and Invasion by Modulating Microtubule Dynamics and Focal Adhesion Turnover. Cancers (Basel) (2021) 13(22):5673.
Hazar-Rethinam M, de Long LM, Gannon OM, Boros S, V argas AC, Dzienis M, et al. (2015). RacGAP1 Is a Novel Downstream Effector of E2F7-Dependent Resistance to Doxorubicin and Is Prognostic for Overall Survival in Squamous Cell Carcinoma. Mol Cancer Ther, 14:1939-1950.
Sun BL. Immunotherapy in treatment of metastatic prostate cancer: An approach to circumvent immunosuppressive tumor microenvironment. Prostate. 2021 Nov;81(15):1125-1134.
Xiong J, He J, Zhu J, Pan J, Liao W, Ye H, Wang H, Song Y, Du Y, Cui B, Xue M, Zheng W, Kong X, Jiang K, Ding K, Lai L, Wang Q. Lactylation-driven METTL3-mediated RNA m6A modification promotes immunosuppression of tumor-infiltrating myeloid cells. Mol Cell. 2022 May 5;82(9):1660-1677.e10.
Sadasivan SM, Chen Y, Gupta NS, Han X, Bobbitt KR, Chitale DA, et al. The Interplay of Growth Differentiation Factor 15 (GDF15) Expression and M2 Macrophages During Prostate Carcinogenesis. Carcinogenesis (2020) 41(8):1074-82.

No competing interests reported.

Download PDF

Version 1

posted

You are reading this latest preprint version

Constructing a lactylation-related genes prognostic model to effectively predict the disease-free survival and treatment responsiveness in prostate cancer based on machine learning

Status:

Version 1

Abstract

Figures

Introduction

Materials and methods

Data acquisition

Identification of differentially expressed and prognostic genes

Unsupervised cluster analysis

Selection of key genes and construction of the prognostic model

Gene set variation analysis (GSVA)

Analysis of tumor mutation burden and immune cell infiltration

Drug sensitivity prediction of the risk model

Construction and evaluation of the nomogram

Cell Lines and Cell Culture

RNA Extraction and Quantitative Analysis

Statistical analysis

Results

Identifying prognostic lactylation-related genes and subtypes in PCa

Establishment and evaluation of the predictive nomogram

Analysis of tumor immune microenvironment of high- and low-risk groups

Functional enrichment analysis and drug sensitivity

Relationship between the LRGs signature and clinicopathological traits and mutation landscape

The mRNA expression levels of lactylation-related genes in PCa cells

Discussion

Conclusion

Abbreviations

Declarations

References

Additional Declarations

Supplementary Files

Status:

Version 1