Comprehensive analysis of fatty acid metabolism-related gene signatures for predicting prognosis in patients with prostate cancer

Data collection

Tissue sample sequencing data and corresponding clinical follow-up data were obtained from the public database The Cancer Genome Atlas (TCGA, https://portal.gdc.cancer.gov/), including 499 cases of prostate cancer tissue and 52 cases of normal tissue. After excluding samples with incomplete clinical information, such as lack of survival time and some clinical characteristics, 489 tumour samples were included. Specific clinical information is shown in Table S1. Microarray data from the GSE21032 dataset from the Gene Expression Omnibus (GEO, https://www.ncbi.nlm.nih.gov/geo/) contained the sequencing data of 218 prostate cancer tissues and matching clinical information. After excluding samples with incomplete clinical information, 140 cases were included. Specific clinical information is shown in Table S2. The two datasets were normalized to remove batch effects.

The TIMER database (https://cistrome.shinyapps.io/timer/) is a public database for studying the abundance of immune cell infiltration in tumour tissues (Li et al., 2017a). It uses RNA-Seq expression profile data to detect the infiltration of immune cells in tumour tissues, including B cells, CD4+ T cells, CD8+ T cells, neutrophils, macrophages, and dendritic cells. This database was used to examine the effect of model-related genes on tumour immune infiltration.

Genomics of Drug Sensitivity in Cancer (GDSC, https://www.cancerrxgene.org/) provides free and publicly available genomic data on tumour therapy (Yang et al., 2013). It is dedicated to the discovery of potential tumour therapeutic targets to improve tumour treatment, and it is the largest public database of isotypes worldwide. Variations in the tumour genome will affect clinical treatment efficacy, and different targets have different responses to drugs. Therefore, the present study used this database to obtain the sensitivity of prostate cancer cells to anti-tumour drugs and their potential therapeutic targets.

The Human Protein Atlas database (HPA) provides tissue and cellular distribution information for all 24,000 human proteins. The expression of fatty acid metabolism-related proteins involved in the model building was investigated in normal and tumour tissues/organs using the HPA database (Colwill & Gräslund, 2011).

Identification of differentially expressed fatty acid metabolism-related genes and cluster analyses

Fatty acid metabolism-related genes were obtained according to earlier research results and the Molecular Signatures Database, as shown in Table S3. Differentially expressed fatty acid metabolism-related genes in prostate cancer tissues and normal tissues were analysed using the “limma package” in R software. We used the STRING online tool to construct the PPI network of differentially expressed fatty acid metabolism genes (Szklarczyk et al., 2021). To study the role of fatty acid metabolism-related genes in the pathological process of prostate cancer, we used the “Consensus Cluster Plus package” to divide patients into different subtypes based on fatty acid metabolism-related gene expression and survival data, where reps were 50 and pItem was 0.8. The “cluster Profiler package” and the “org.Hs.eg.db package” were used to annotate the differentially expressed genes between different types. |Log2fold change| > 0.58 and a false discovery rate (FDR) <0.05 were the cut-off criteria for differential analysis.

Prognostic model construction and verification

First, differentially expressed fatty acid metabolism-related genes associated with RFS in prostate cancer patients were determined using univariate analysis, and filtering criteria were P < 0.05. Second, the genes related to fatty acid metabolism associated with RFS were used to construct a prognostic model using LASSO Cox regression analysis (glmnet package) (Tibshirani et al., 2012). The following risk score calculation method was used: risk score = (coefficient mRNA 1 × expression of mRNA 1) + (coefficient mRNA 2 × expression of mRNA 2) +… (coefficient mRNA n × expression of mRNA n). The risk score of each sample in the TCGA-prostate cancer dataset was calculated. All samples were divided into high-risk and low-risk groups according to the median value of the risk score. The Kaplan-Meier method was used to analyse the prognostic differences between the high- and low-risk groups of patients. Receiver operating characteristic (ROC) curves were drawn to assess the accuracy and effectiveness of the risk model in predicting the prognosis of prostate cancer patients. Principal component analysis (PCA) and t-distributed stochastic neighbour embedding (t-SNE) methods were used to detect whether the prognostic model had reliable clustering ability.

Independence of fatty acid metabolism-related gene models for the prognostic diagnosis of patients with prostate cancer

To clarify whether the fatty acid metabolism-related gene prognostic model had independent diagnostic significance for the RFS of prostate cancer patients, we performed univariate and multivariate analyses on the TCGA dataset and GEO dataset. Relevant clinical prognostic factors included age, TNM stage, and metastasis. P < 0.05 was considered statistically significant.

Nomogram

A nomogram shows a predictive model intuitively and effectively, and it may be used for multi-index joint diagnosis or prediction of disease onset or progression. Nomograms were drawn to visualize the prognostic diagnostic model of this study. The nomograms quantified the risk score according to the stages of different clinical features and predicted the time of disease-free survival according to the total score. Nomogram plotting and calibration were performed using the R package “rms”.

Gene set enrichment analysis

Gene set enrichment analysis (GSEA) is an efficient and widely used tool for the biological functional annotation of a set of genes (Subramanian et al., 2005). We divided 489 prostate cancer samples into high- and low-risk groups based on prognostic risk scores and performed GSEA. The gene set database was set to “KEGG cell signalling pathways”. P value < 0.05, q-value < 0.25, and |NES| > 1 were considered significantly enriched.

Immune analysis

To examine differences in immune infiltration and the tumour immune microenvironment between the high- and low-risk groups, we performed a series of analyses. First, the stromal score, immune score, estimate score, and tumour purity were evaluated for each sample using the ESTIMATE algorithm, and the differences in relevant scores between the two groups were compared (Siemers et al., 2017). Second, the CIBERSORT algorithm was used to evaluate differences in immune cell infiltration for each sample in the two groups (Yoshihara et al., 2013). Single-sample gene enrichment analysis (ssGSEA) is an immune infiltration estimation method that uses the GSVA, GSEABase, and limma packages of the R language to estimate the degree of immune infiltration by calculating the expression levels of 29 immune gene sets in a single sample (Hänzelmann, Castelo & Guinney, 2013). In this study, ssGSEA was used to quantify the abundance of immune cells in the two groups of samples. Immune checkpoint correlation analysis was also performed to clarify the relationship between the sample risk score and immune checkpoints between the two groups. Immune checkpoints were accessed from previous literature (Tang et al., 2021).

Somatic mutation analysis

Somatic mutation information of prostate cancer samples was obtained from TCGA. The mutation profiles of patients in the high- and low-risk groups of prostate cancer patients were analysed using the “maftool” package (Mayakonda et al., 2018). Tumour mutational burden (TMB) is the total number of somatic gene mutations per million bases. All samples were divided into a high TMB group and a low TMB group according to the median TMB values. The association of TMB with patient prognosis was evaluated using survival analysis and the relationship of TMB with the model risk score.

Cell culture

Prostate cancer cells (22RV1 and DU145) were obtained from the China National Biomedical Experimental Cell Resource Bank. The following culture conditions were used: 22RV1 cells were cultured with RPMI 1640 medium containing 10% fetal bovine serum. DU145 cells were cultured with MEM containing 10% fetal bovine serum. All cells were cultured in a 37 °C incubator with 5% CO2.

Quantitative reverse transcription PCR

Total RNA was extracted from cells using a Total RNA Kit I (R6834-02, Omega, USA). RNA quality and purity were checked using a NanoDrop (Thermo, USA). cDNA was obtained by reverse transcription using a NovoScript Plus All-in-one 1st Strand cDNA Synthesis SuperMix kit (E047-01B; Novoprotein, Shanghai, China). RT-qPCR was used to detect the expression level of ARPC1B using a Novostart SYBR qPCR SuperMix Plus Kit (E096-01A; Novoprotein, Shanghai, China). The primer sequences required for the experiments are shown in Table S4. GAPDH was used as an internal reference, and the 2^−ΔΔCT calculation method was used.

Statistical analysis

Data analyses and processing were performed using R software (version 4.1.1). The Mann-Whitney U test was used to detect differential genes between normal and tumour tissues. RFS survival analysis was performed using the Kaplan–Meier method. To evaluate whether the predictive model could be used as an independent prognostic factor, univariate and multivariate Cox analyses were performed. The Wilcoxon test was used to analyse differences in immune cell types and immune checkpoints between the two groups. P < 0.05 was considered significantly different.

Fatty acid metabolism related DEGs between normal and tumour tissues

Firstly, the research flow chart of our study is shown in Fig. 1. The data for 489 prostate cancer samples were downloaded from the TCGA. To examine the biological function and role of fatty acid metabolism-related genes (FRGs) in the pathological process of prostate cancer, we performed gene differential expression analysis on tumour tissues and normal tissues obtained from the TCGA database. The differential expression results are shown in Fig. 2A. Blue shows low expression, and orange shows high expression. To examine the potential biological connection of these differentially expressed FRGs, we performed a PPI analysis. The results showed the interaction and relationship between the proteins encoded by these genes (Fig. 2B). These initial studies obtained the FAM-related DEGs and explored their connections.

Clustering of tumours based on fatty acid metabolism-related genes

To examine whether prostate cancer may be divided into different subtypes based on FRGs, consensus clustering analysis was performed on 489 prostate cancer patients from TCGA. The results showed that the subgroup classification was stable and credible when the number of clusters was 2 (Fig. 3A). The 489 prostate cancer samples could be divided into two subtypes based on fatty acid metabolism-related genes. Survival analysis was performed to clarify the difference in the prognosis of the two subtypes of patients, and the results showed that C1 had a better prognosis than C2 (Fig. 3B). Cluster analysis was performed for the two subtypes. As shown in Fig. 3C, the gene expression patterns were different between the two subtypes, and the clinical features were significantly different, which indicated that FRG-based tumour classification was feasible.

Construction of a prognostic model of prostate cancer based on fatty acid metabolism-related genes

These results showed that prostate cancer could be divided into two subtypes based on fatty acid metabolism-related genes. However, whether fatty acid metabolism-related genes were related to the RFS of prostate cancer patients and whether a prognostic model could be established were not known. Therefore, we performed univariate analysis for differentially expressed fatty acid metabolism-related genes, and 21 genes were related to the RFS of prostate cancer patients according to the shear criterion P < 0.01 (Fig. 4A). We performed least absolute shrinkage and selection operator (LASSO) regression analysis based on these genes. The results of the algorithm analysis (Figs. 4B and 4C) identified 10 candidate genes (ADH5, RXRA, PTGES3, D2HGDH, NUDT7, EPHX2, HAO2, CPT1C, LTC4S, and SLC25A17) that were selected for the construction of the prognostic model. The following risk score calculation formula was used: Risk score = (−0.25*ADH5 exp) + 0.25*RXRA exp + 0.22*PTGES3 exp + 0.24*D2HGDH exp + (−0.24*NUDT7 exp) + (−0.07*EPHX2 exp) + 0.16*HAO2 exp + 0.56*CPT1C exp + 0.09*LTC4S exp + (0.17*SLC25A17 exp). The risk score of 489 prostate cancer samples was calculated based on the formula. All samples were divided into high- and low-risk groups according to the median score (Figs. 4D and 4E). The PCA and tSNE analysis results showed that the distribution of patients in the two groups had mutual aggregation within the group, and there was good differentiation between the groups (Figs. 4F and 4G). This result indicated that the sample risk grouping had good reliability and discrimination. To clarify the difference in the RFS of the two groups of patients, we performed a survival analysis. The results showed that patients in the high-risk group had shorter recurrence-free survival and poor prognosis (Fig. 4H). The ROC curve showed that the area under the curve (AUC) values were 0.813, 0.799, and 0.790 at 1 year, 3 years, and 5 years, respectively (Fig. 3I). The clinical correlation analysis showed that the risk score increased with increasing clinical T, pathological T and N stages of prostate cancer, and the PSA and Gleason scores of patients in the high-risk group were higher than the low-risk group (Fig. S1). The expression level of the gene-encoded protein constructed by the model was verified using the HPA database (Fig. S2). These results showed that the prognostic risk model was dependable and had some diagnostic value for the prognosis of prostate cancer patients.

Validation of the prognostic model

To further evaluate the validity and reliability of the prognostic model, the GSE21032 dataset containing 140 prostate cancer samples was used for external verification. First, the risk score of the sample in the GSE21032 dataset was calculated using the risk score calculation formula. The sample was divided into high- and low-risk groups based on the median risk score, and patients in the high-risk group had a larger recurrence ratio (Figs. 5A and 5B). Second, the PCA and tSNE analysis results showed good discrimination between the two groups of patients (Figs. 5C and 5D). This result indicated that the prognostic signature had high reliability. Survival analysis showed that the RFS of patients in the high-risk group was significantly worse (Fig. 5E). The ROC curve showed that the 1-, 3-, and 5-year AUC values were 0.890, 0.754, and 0.701, respectively (Fig. 5F.) This result indicated that the model had a certain diagnostic value for the RFS of prostate cancer patients in the GSE21032 dataset. These results showed that the prognostic risk model established in this study was effective for the external dataset GSE21032.

The prognostic model is an independent prognostic factor in prostate cancer patients

The above results showed that the prognostic risk model established in this study was reliable and had a certain diagnostic value for the RFS of patients, but whether it could be used as an independent prognostic factor for prostate cancer patients was not known. Therefore, we performed univariate and multivariate regression analyses on patients from the TCGA and GEO datasets. The results of the univariate analysis showed that the risk score could be used as an independent prognostic factor in both datasets (Figs. 6A and 6C, TCGA: hazard ratio (HR) = 14.841, 95% confidence interval (CI) [7.389–29.808], P < 0.001; GEO: HR = 9.227, 95% CI [3.841–22.165], P < 0.001). Multivariate analysis also revealed that the risk score could be used as an independent prognostic factor in both datasets (Figs. 6B and 6D, TCGA: HR = 12.345, 95% CI [5.908–25.794], P < 0.001; GEO: HR = 6.247, 95% CI [2.691–14.5], P < 0.001). Cluster analysis was performed on model building genes and clinical features of the TCGA dataset. The results showed that the clinical characteristics of the two groups were significantly different (Fig. 6E). Briefly, these studies suggest that the model constructed in this study was an independent prognostic factor for prostate cancer patients.

Construction of a predictive nomogram

To visualize the model, we constructed a model prediction nomogram (Fig. 7A). The nomogram included risk score, age, clinical T stage, pathological N stage, and pathological T stage. The nomogram was used to assess the prognosis of patients with prostate cancer. We performed graphic calibration to verify the validity of the nomogram. The x-axis in the calibration curve represents the nomogram-predicted recurrence rate of prostate cancer patients, and the y-axis represents the actual recurrence rate. The grey line represents the ideal predictive effect, and the orange line represents the actual predictive effect of the nomogram. As shown in Figs. 7B–7D, the coincidence of the two lines was higher, which indicated that the nomogram had a better prognostic prediction effect on patients’ 1-, 3-, and 5-year outcomes.

Gene set enrichment analysis

To examine the potential differences in biological functions and cell signalling pathways between the high- and low-risk groups, we performed GSEA. According to the shear criteria (P value < 0.05 and q-value < 0.25), a total of 31 cell signalling pathways were enriched (Table S5). The top 5 pathways enriched in the high-risk group were “cell cycle”, “glycerolipid metabolism”, “oocyte meiosis”, “primary immunodeficiency” and “ribosome” (Fig. 8A). The top 5 pathways enriched in the low-risk group were “arrhythmogenic right ventricular cardiomyopathy”, “cardiac muscle contraction”, “dilated cardiomyopathy”, “hypertrophic cardiomyopathy” and “propanoate metabolism” (Fig. 8B). In addition to these partial immune-related pathways, multiple immune activity-related cell signalling pathways were enriched in the high-risk group, such as “ECM-receptor interaction”, “cytokine‒cytokine receptor interaction”, and “natural killer cell-mediated cytotoxicity”.

Immune signature and immune checkpoint analysis

To examine the correlation between prognostic models and immune characteristics, we performed a series of analyses. First, the ESTIMATE algorithm was used to clarify the differences in immune and stromal components in the immune microenvironment between the two groups of patients. The results showed that the ESTIMATE score, immune score, and stromal score of the high-risk group were higher than the low-risk group, and tumour purity exhibited the opposite results (Figs. 9A–9D). Survival analysis showed that patients with low tumour purity had a poor prognosis, and the remaining three immune scores negatively correlated with patient prognosis (Fig. S3). These results indicated that there were differences in the proportion of tumour and immune cells in the tumour microenvironment between high- and low-risk groups, which have implications for patient prognosis. The abundance of various immune cell infiltrates in the samples of the two groups is not known. Therefore, we used CIBERSORT for further analysis of the tumour immune microenvironment. The results are shown in Fig. 9E. The differences in plasma cells, CD8 T cells, regulatory T cells (Tregs), gamma delta T cells, resting NK cells, M0 macrophages, M2 macrophages, resting dendritic cells, and neutrophils were statistically significant. Notably, the immune cell score had a certain correlation with the model risk score. Survival analysis revealed that the expression levels of Tregs, CD8 T cells, plasma cells, activated NK cells, M2 macrophages, and resting dendritic cells negatively correlated with patient RFS (Fig. S4).

To assess the relationship between the risk score and immune status, the relevant immune cell gene set was quantified using ssGSEA. The results showed that the immune scores of APC co-stimulation, checkpoint, cytolytic activity, DCs, HLA, inflammation promotion, pDCs, T-cell co−inhibition, T-cell co−stimulation, T helper cells, Th1_cells, TILs, Tregs and type I IFN response were significantly different between the high- and low-risk groups (Fig. 9F). Survival analysis showed that most immune cells or function scores (Type I IFN response, TIL, T-cell co-inhibition, T helper cells, T-cell co-stimulation, pDCs, HLA, inflammation promotion, DCs, cytolytic activity, APC co-inhibition, and checkpoint) negatively correlated with the patient’s RFS. Treg and Th1 cells positively correlated with patient RFS (Fig. S5).

Immune checkpoint inhibitors are an important breakthrough in the field of oncology. These inhibitors have been clinically used in a variety of solid tumours, and the effects are good. To examine the relationship between risk scores and immune checkpoints, we first performed a correlation analysis. The results showed that the expression levels of immune checkpoints, such as CD274, CTLA4, PDCD1, HAVCR2, PDCD1LG2, and IDO1, positively correlated with risk scores (Figs. 10A–10F). We also analysed the expression levels of other immune checkpoints in the two sets of samples. The results showed that the expression levels of the high-risk group were significantly higher than the low-risk group (Fig. 10G). Taken together, these findings suggest that the fatty acid metabolism risk model is a predictor of immunotherapy response.

Tumour somatic mutation characteristics and the relationship between model-related gene copy number variation and immune infiltration

Somatic mutations in the genome continue to accumulate during the pathological process of tumours to drive the occurrence and development of tumours. Tumour gene mutations are of great significance to tumour treatment, drug resistance monitoring, and prognostic evaluation. Therefore, we performed genetic mutation analysis on samples from the high- and low-risk groups, and the top 20 genes with the highest mutation frequency in the two groups are shown in Figs. 11A and 11B. We also performed an analysis of the relationship of TMB with the model risk score and prognosis. Differences in TMB between the high- and low-risk groups are shown in Fig. 11C. The high-risk group had significantly higher TMB than the low-risk group. The relationship between TMB and risk score was positively correlated (Fig. 11D). The survival analysis results showed that high TMB tended to have a poor prognosis (Fig. 11E and 11F). Taken together, these results indicate that there are differences in TMB between patients in the high- and low-risk groups and indirectly indicate that the two groups of patients respond differently to immunotherapy.

We also used the TIMER analysis tool to clarify the relationship between model-related gene copy number variation (CNV) and immune cell infiltration. The results showed that the CNVs of model-associated genes induced changes in the levels of partial immune cell infiltration. As shown in Fig. 12, deep deletion, arm-level deletion, and high amplification reduced the infiltration of some immune cells. We examined the relationship between the expression of model-related genes and the abundance of immune cell infiltration using TIMER, and the results showed a certain correlation. These results suggested that expression level is a potential factor affecting the immune microenvironment of prostate cancer (Fig. S6).

Drug sensitivity analysis

GDSC is a practical database containing the sensitivity and response of various tumour cells to drugs. To clarify whether the model-related genes were potential drug treatment targets for prostate cancer and which drugs were more sensitive to this target, we performed correlation analysis. The results showed that model-related genes in prostate cancer cells were related to the drug resistance or sensitivity of some tumour-targeted drugs (Fig. 13). The expression of CPT1C, D2HGDH, HAO2, LTC4S NUDT7, PTGES3 and SLC25A17 negatively correlated with the sensitivity to rapamycin and imatinib. CPT1C, LTC4S, NUDT7 and RXRA positively correlated with sensitivity to lapatinib and salubrinal. The expression of CPT1C, LTC4S, NUDT7, RXRA and SLC25A17 positively correlated with sunitinib sensitivity. These results provide new potential ideas for targeted therapy of prostate cancer.

Verification of the expression of the model construction gene in PCa cells

To further validate the effectiveness of the signature, we used the primary human prostate cancer cell line 22RV1 and the human prostate cancer metastasis-derived cell line DU145 as the study subjects and detected the expression levels of the signature in prostate cancer cells of different aggressiveness levels using RT-qPCR. As presented in Fig. 14, CPT1C, NUD7T, PTGES3 and RXRA were elevated in DU145 cells compared to 22RV1 cells. D2HGDH, EPHX2, HAO2 and SLC25A17 were expressed at low levels in DU145 cells. Taken together, these model genes showed significant differences in PCa cell lines, which indirectly confirmed the importance and validity of these model genes in PCa.