Fatty acid metabolism-related genes

Three gene sets related to fatty acid metabolism—Hallmark fatty acid metabolism genes, Kyoto Encyclopedia of Genes and Genomes (KEGG) fatty acid metabolism pathways, and Reactome fatty acid metabolism genes—were extracted from the Molecular Signature Database v7.5.1 (MSigDB; https://www.gsea-msigdb.org/gsea/msigdb), and 309 fatty acid metabolism-related genes were retrieved after removing overlapping genes (Table S1). One of the fatty acid metabolism-related genes (XIST) was unavailable in the TCGA, leaving 308 genes for further analyses.

Identification of differentially expressed genes (DEGs)

DEGs were identified by comparing the mRNA expression of 308 fatty acid related genes between tumor and normal tissue in TCGA datasets using the “limma” package. Genes with a false discovery rate (FDR) <0.05 and |log FC| >0.5 were selected for further analysis. DEGs between the low- and high-risk groups were identified in TCGA after the calculation of the risk score using the “limma” package. Genes with adjusted FDR <0.05 and |log FC| >0.5 were selected for further analysis.

Establishment and validation of the risk predictive model

Univariate Cox regression analysis was used to determine the fatty acid-related genes whose expression levels were significantly correlated with BCR-free survival in patients with PCa (P < 0.05). The least absolute shrinkage and selection operator (LASSO) Cox regression method was applied to these genes. Finally, multiple stepwise Cox regression, with both backward and forward selection, was performed on the candidate fatty acid-related genes, and the prognostic model was determined based on the lowest Akaike information criterion (AIC) value. The fatty acid-related gene score was calculated based on the gene expression level and the corresponding regression coefficient. The formula was established as follows: Risk Score = expression level of gene1*Coef1 + expression level of gene2*Coef2 + … + expression level of genex*Coefx (Jiang et al., 2022). The median score was used as the cut-off value to divide TCGA patients with PCa into high- and low-risk groups. The “survival”, “glmnet” and “survminer” R package were used to establish the risk predictive model. The same formula was applied to the GEO dataset for validation purposes.

Univariate and multivariate Cox proportional hazards regression analyses were performed to test whether the risk score was an independent prognostic factor in the TCGA and GEO databases. The Kaplan–Meier (KM) survival curve was used to compare the BCR-free survival of patients in the high- and low-risk groups. Principal component analysis (PCA) and t-distributed stochastic neighbor embedding (t-SNE) were used to explore the distribution of different populations. Time-dependent receiver operating curve (ROC) curve analysis was used to evaluate the predictive power of the gene signatures and clinicopathological features.

Functional enrichment analysis

The list of DEGs between tumor and normal tissues in the TCGA dataset was subjected to the Metascape online tool (https://metascape.org/) for functional enrichment analysis. The “clusterProfiler” R package was used to perform Gene Ontology (GO) and KEGG analyses based on the DEGs between different risk groups. The infiltration scores of 28 immune cells and the activities of 13 immune-related pathways were calculated using the “GSVA” R package by single-sample gene set enrichment analysis (ssGSEA).

Prediction of responses to immunotherapy and chemotherapy

Tumor immune dysfunction and exclusion (TIDE) score and tumor mutation burden (TMB) were used to assess the clinical response to immune therapy in patients with PCa. A low TIDE score suggested high sensitivity to immunotherapy (Jiang et al., 2018), while a high TMB was suggestive of high sensitivity to immunotherapy (Yarchoan, Hopkins & Jaffee, 2017). The TIDE score was calculated according to algorithm proposed by Jiang et al. (2018) in a previous study. We obtained the TCGA mutation dataset from the TCGA website (https://portal.gdc.cancer.gov/repository), and each patient’s TMB was measured as the total number of non-synonymous mutations per megabase. In addition, chemotherapeutic response to three commonly used drugs (bicalutamide, cisplatin, and docetaxel) for patients with PCa in the TCGA dataset was calculated using the “pRRophetic” R package.

Cell culture

The normal prostate epithelial (RWPE-1) cells (Cellcook, Guangzhou, China) were cultured in Keratinocyte SFM medium (Sciencell, Carlsbad, CA, USA) containing 1% bovine pituitary extract (M&C Gene, Beijing, China) and 1% epidermal growth factor (Sciencell, Carlsbad, CA, USA). PC3 and LNCaP cells (Cellcook, Guangzhou, China) were cultured in RPMI-1640 medium (Gibco, Waltham, MA, USA) containing 10% fetal bovine serum (Bioind, Cromwell, CT, USA). All cells were cultured in a water-saturated atmosphere with 5% CO₂ at 37 °C.

RNA isolation and real-time PCR

Total RNA was extracted using TRIzol reagent according to the manufacturer’s instructions. The concentration and purity of all the RNA samples were determined at an absorbance ratio of 260/280 nm. The PrimeScript RT Reagent Kit with gDNA Eraser (Takara, San Jose, CA, USA) was used for cDNA synthesis. Quantitative real-time PCR (qPCR) using the Fast SYBR Green PCR Master Mix (Applied Biosystems, Waltham, MA, USA) was performed on a StepOne Real-Time PCR System (Applied Biosystems, Waltham, MA, USA). The relative level of mRNA expression of a gene was determined by normalization to GAPDH. The 2^−ΔΔCT method was used for data analysis. The primers used for real-time PCR were as follows:

ALDH1A1-F: 5′-CCGTGGCGTACTATGGATGC-3′;

ALDH1A1-R: 5′-GCAGCAGACGATCTCTTTCGAT-3′;

CA2-F: 5′-AAGGAACCCATCAGCGTCAG-3′;

CA2-R: 5′-TTCTTCAGTGGCTGAGCTGG-3′;

CPT1B-F: 5′-CCTGCTACATGGCAACTGCTA-3′;

CPT1B-R: 5′-AGAGGTGCCCAATGATGGGA′;

CROT-F: 5′-ATTGGCTGGAAGAGTGGTGG-3′;

CROT-R: 5′-GAGTCCCTTCCTTTGGAGGC-3′;

NUDT19-F: 5′-CCCCACAGTTCTACGAAGTGA-3′;

NUDT19-R: 5′-TCTAATGCACGACCCAAACAAA-3′;

Construction of a fatty acid metabolism-related gene signature

Twelve candidate fatty acid-related DEGs associated with BCR-free survival were identified using the univariate Cox regression analysis (Fig. 2A). We then performed LASSO-Cox regression analysis on the identified prognostic DEGs. As shown in Figs. 2B and ​and2C,2C, using the minimum lambda criteria, nine mRNAs were selected for further analysis. The coefficients of these mRNA are shown in Fig. 2D. To identify the fatty acid-related genes with the greatest prognostic value, we conducted a multiple stepwise Cox regression to construct a risk signature for BCR-free survival prediction, and five genes were selected: ALDH1A1, CA2, CPT1B, CROT, and NUDT19 (Fig. 2E). Further, the expression levels of the five genes were validated in PC3 and LNCaP cell lines by qPCR. As shown in Fig. 3B, the expression levels of ALDH1A1 and CA2 were significantly downregulated in PCa cells compared with those in RWPE-1 control cells, while the expression of CPT1B and NUDT19 was higher in PCa cells, consistent with the findings from the TCGA dataset. However, expression of CROT was reversed compared to that in the TCGA database (Fig. 3A); further studies in other PCa cells lines and PCa tissue are needed to confirm these findings.

Open in a separate window

Figure 2

Establishment of the prognostic model in the TCGA cohort.

(A) Results of the univariate cox analysis of the biochemical recurrence-free survival in the TCGA cohort. (B) Ten-time cross-validation for tuning parameter selection in the LASSO Cox-regression model in the TCGA cohort. (C) LASSO coefficient profiles of the nine fatty acid-related genes in the TCGA cohort. (D) Coefficients of the LASSO-Cox model with the minimum criteria in the TCGA cohort. (E) Results of multiple stepwise Cox regression of the biochemical recurrence-free survival in the TCGA cohort.

Open in a separate window

Figure 3

Validation of the mRNA expression levels of the five fatty acid metabolism related-genes in cell lines.

(A) Expression of the five fatty acid metabolism-related genes in normal and tumor samples in the TCGA dataset. (Normal samples: 52; Tumor samples: 381) (B) Expression of the five fatty acid metabolism-related genes in two prostate cancer cell lines (PC3 and LNCAP) and normal prostate cell lines (RWPE-1). (Number of experimental repetitions: 3) *P < 0.05, **P < 0.01, ****P < 0.0001, ns: not significant.

Clinicopathological characteristics associated with the 5-gene risk signature

Risk score was calculated for each patient according to the coefficients of selected genes. Patients in the TCGA training cohort were classified as high- (n = 190) or low-risk (n = 191) based on the median cutoff value of the risk score. As shown in Fig. 4A, patients with high risk scores were more likely to have BCR than those with low risk scores. The results of PCA and t-SNE analyses indicated that patients in the high- and low-risk groups were distributed in discrete directions (Figs. 4B and ​and4C).4C). Additionally, patients in the high-risk group had poorer BCR-free survival than those with low risk scores (Fig. 4D). The area under the ROC curve for 1-year, 3-year and 5-year BCR-free survival were 0.760, 0.806, and 0.793, respectively (Fig. 4E). Furthermore, as indicated in Figs. 4F–4J, tumors with higher risk score were associate with larger tumor size, more likely to have lymph node involvement, and higher Gleason score, albeit no significant associations with age and prostate-specific antigen (PSA) value were observed. Moreover, multivariate Cox regression analyses showed that the risk score remained independently predictive of BCR-free survival after controlling for major confounders including tumor size, PSA at diagnosis, and Gleason score (HR = 1.11, 95% CI [1.05–1.17], P < 0.001, Table 1). In contrast, the Gleason score rendered insignificant after controlling for confounders. These findings suggest that our risk signature might serves as a potential powerful prognostic indicator in guiding clinical practice.

Open in a separate window

Figure 4

Clincal associations of the fatty acid metabolism-related gene score in prostate cancers from TCGA cohort.

(A) The distribution of risk score, biochemical recurrence-free survival status, and the exoression of the five hub genes between the low- and high-risk groups in the TCGA cohort. (B) The PCA analysis of risk score in the TCGA cohort. (C) The t-SNE analysis of risk score in the TCGA cohort. (D) Kaplan–Meier analysis between the low- and high-risk groups in the TCGA cohort. (E) ROC curve and AUC of risk score for predicting the 1/3/5-years survival in the TCGA cohort. (F–J) Risk score by patient age, Gleason score, PSA value, tumor pathological T-stage, and N-stage in the TCGA cohort.

Validation of the 5-gene risk signature

To validate the prognostic significance of the 5-gene risk signature, the corresponding risk score was calculated for each patient in the GEO cohort using the same formula; then, patients were assigned into high- and low-risk group according to the medial value of the score. As shown in Fig. 5A, the expression patterns of the five prognostic genes in GEO cohort were similar to that in the TCGA cohort. Patients with high risk were more likely to have BCR than those with low risk scores. Likewise, the results of PCA and t-SNE analyses indicated that patients in the high- and low-risk groups were distributed in discrete directions (Figs. 5B and ​and5C).5C). High-risk group also had poorer BCR-free survival compared with those in the low-risk group, consistent with the observation in the TCGA dataset (Fig. 5D). The area under the ROC for 1-year, 3-year, and 5-year BCR-free survival was 0.730, 0.734, and 0.716, respectively (Fig. 5E), comparable to that of the TCGA cohort. Moreover, in line with the findings from the TCGA cohort, the multivariate Cox-regression also revealed that the 5-gene risk score was a robust predictor for BCR-free survival independent of tumor size and Gleason score in the GEO cohort. Taken together, our results demonstrated that our fatty acid metabolism-related 5-gene signature is an independent prognostic predictor of BCR in prostate cancer.

Open in a separate window

Figure 5

Clincal associations of the fatty acid metabolism-related gene score in prostate cancers from GEO cohort.

(A) The distribution of risk score, biochemical recurrence-free survival status, and the exoression of the five hub genes between the low- and high-risk groups in the GEO cohort. (B) The PCA analysis of risk score in the GEO cohort. (C) The t-SNE analysis of risk score in the GEO cohort. (D) Kaplan–Meier analysis between the low- and high-risk groups in the GEO cohort. (E) ROC curve and AUC of risk score for predicting the 1/3/5-years survival in the GEO cohort.

Association with sensitivity to immunotherapy/chemotherapy

Immunotherapy has been the focus of attention in cancer treatment; however, only a small proportion of patients respond to immunotherapy, highlighting the unmet need of biomarkers for response prediction. TMB is an emerging biomarker for immunotherapy response and has been approved as a companion diagnostic marker for immunotherapy (Marcus et al., 2021). To determine whether the 5-gene risk score could serve as a potential biomarker for immunotherapy, we evaluated its associations with TMB. As shown in Fig. 6A, a higher TMB score was observed in the high-risk group compared to the low-risk group. Although the TMB score alone was not predictive of BCR-free survival, the combination of TMB and risk score could identify a subgroup of patients with high-risk scores and low-TMB scores who had the worst survival (Figs. 6B and ​and6C).6C). Moreover, significantly lower TIDE scores were observed in the high-risk group than the low-risk group, suggesting a favorable response to immunotherapy associated with a high risk score (Fig. 6D). Next, we determined whether the risk score was also associated with chemotherapy response. The half-maximal inhibitory concentrations (IC₅₀) of three commonly used chemotherapeutic drugs (bicalutamide, cisplatin, and docetaxel) in high- and low-risk tumors were compared. As shown in Figs. 6E–6G, high-risk score was related to high IC₅₀ values for bicalutamide and docetaxel, which indicates a lower sensitivity to chemotherapy. In summary, the above results suggested that patients with high risk score might benefit from immunotherapy rather than chemotherapy.

Open in a separate window

Figure 6

Association with sensitivity to immunotherapy/chemotherapy in the TCGA cohort.

(A) Tumor mutation burdern (TMB) between the high-risk and low-risk groups in the TCGA cohort. (B) Kaplan–Meier curve analysis of biochemical recurrence-free survival is shown for patients with high and low TMB values in the TCGA cohort. (C) Kaplan–Meier curve analysis of biochemical recurrence-free survival is shown for patients classified according to the TMB value and fatty acid-related model in the TCGA cohort. (D) Tumor immune dysfunction and exclusion (TIDE) score between the low- and high-risk groups in the TCGA cohort. (E–G) The drug sensitivity of bicalutamide (E), cisplatin (F), and docetaxel (G) between the low- and high-risk groups in the TCGA cohort.