DOI: https://doi.org/10.21203/rs.3.rs-1552212/v1
Background: Hepatocellular Carcinoma (HCC) is one of the most common malignant tumors with high mortality. Fatty acid metabolism is associated with the development and treatment of HCC. This study aimed to build a fatty acid metabolism-related gene model to evaluate the prognosis of HCC patients and provide guidance for individualized treatment.
Methods: RNA-sequencing data of HCC from The Cancer Genome Atlas (TCGA) and Gene Expression Omnibus (GEO) database were extracted as the training set and validation set. Fatty acid metabolism-related genes predictive signature was built using univariate and least absolute shrinkage and selection operator (LASSO) Cox regression analyses. The predictive value of prognostic model was evaluated by Kaplan-Meier (K-M) curves, ROC curves, univariate and multivariate Cox regression. In addition, the immune infiltration and drug sensitivity were also explored.
Results: A twelve-gene fatty acid metabolism-related risk signature was constructed. Patients were divided into low-risk and high-risk groups according to risk scores. High-risk patients with the higher risk scores had a poor prognosis. Then, univariate and multivariate Cox regression showed that risk score was an independent risk factor for the overall survival of HCC. K-M and ROC curves further verified the predictive value of the signature. Besides, nomograms were drawn and exhibited the good performance in the prognostic prediction. Finally, the immune cell infiltration and sensitivity to chemotherapy drugs exhibited a correlation with different risk levels.
Conclusion: This novel fatty acid metabolism-related gene predictive signature may be useful for improving prognosis evaluation and individualized treatment for HCC patients.
Liver cancer is reported to be the third leading cause of cancer-related death worldwide, only behind lung cancer and colorectal cancer[1], and even more than 1 million patients will die of live cancer in 2030 from the prediction of World Health Organization[2]. As the most common type of liver cancer, hepatocellular carcinoma (HCC) accounts for 90% of primary liver cancer[3]. Presently, treatment for HCC mainly concludes surgical therapy, chemotherapy, immunotherapy, and radiotherapy, but the actual long-term survival rate of HCC patients remains unsatisfactory[4]. Studies showed the genomic and transcriptional heterogeneity among HCC patients may lead to relapse, drug resistance and the poor prognosis[5]. Therefore, the valid biomarkers are needed for individualized treatment and prognostic prediction for patients with HCC.
Metabolism rewiring can disrupt the cell homoeostasis, cause the excessive cell growth and proliferation, and further promote the cancer. Changes in lipid metabolism acclimatize cancer cells to the various microenvironment by resisting oxidative stresses, adjusting intercellular communication, regulating the immune responses, and sustaining key oncogenic functions[6]. Therefore, lipid profiles are emerging as the biomarkers for the prognostic prediction in patients with cancer. Zhu et al developed a signature based on lipid metabolism-related genes and verified its good performance for prognostic prediction in patients with bladder cancer[7]. Recently, as the important element of lipid, abnormality of fatty acid metabolism in cancer attracts increasing attention[8], which affects cancer cell biology in various aspects. Fatty acid not only regulates the synthesis of biological membranes, but also takes part in the modulation of oncogenic signaling. In addition, fatty acid also serves as substrates for mitochondrial energy synthesis and the form of energy storage[9]. Lu et al. revealed the contribution of fatty acid metabolism in maintaining HCC cells survival to resist energy deprivation[10], so targeting the fatty acid metabolism may be a promising approach for the treatment of HCC. However, due to the diversity of genotype and oncobiology, it is difficult to stratify in HCC patients and thus limits the drug discovery based on the fatty acid metabolism. Notably, the development of transcriptomics and metabolomics reveal the novel approach in HCC grouping, and may provide clinical guidance targeting the fatty acid metabolism. Ding et al. founded a prognostic model in colorectal cancer with fatty acid metabolism-related genes and showed that low-risk patients were more sensitive to chemotherapeutics[11]. Consequently, targeting the fatty acid metabolism may contribute to novel signature for the prognostic prediction and clinical medication for HCC patients.
Here, we aimed to develop a fatty acid metabolism-related gene signature to measure the risk scores, improve the accuracy of prognosis prediction, explore the drug sensitivity and immune cell infiltration, and further provide the guidance for clinical treatment in patients with HCC.
The raw RNA sequencing (RNA-seq) data of HCC patients were downloaded from The Cancer Genome Atlas (TCGA) database on 14 Mar 2022. Meanwhile, the clinical information such as gender, age, histological grade, TNM stage, survival time and survival status were also acquired. Besides, the GSE14520 (GPL3921) dataset, including the gene expression and clinical information of HCC patients, was downloaded from the Gene Expression Omnibus (GEO) database. Notably, patients with incomplete clinical data were excluded. Finally, 370 HCC samples from TCGA database were selected as the training set and 221 HCC samples from GEO database were used as the validation set.
Kyoto Encyclopedia of Genes and Genomes (KEGG) gene sets were downloaded from the Molecular Signature Database (MSigDB) and fatty acid metabolic genes were extracted in fatty acid metabolism pathways. A total of 309 fatty acid metabolism-related genes were extracted, and then 137 common genes in TCGA and GEO cohorts were identified (Table S1). Moreover, the differentially expressed fatty acid metabolism-related genes (DEFAMGs) between normal liver tissues and HCC patients were obtained by the “limma” R package on the basis of following criteria: False Discovery Rate (FDR) < 0.05 and |log FC| ≥ 0.585.
To further illustrate the major biological features and cell functional pathways, this study made the Gene Ontology (GO) and KEGG pathway enrichment analyses on DEFAMGs using “cluster Profiler” R package. In addition, “enrichplot” and “ggplot2” R packages were employed to visualize the results. P value < 0.05 was considered statistically significant difference.
The univariate Cox regression analysis was performed to distinguish the genes significantly related to patient overall survival (OS) from DEFAMGs in the training set. P value < 0.01 was considered significant. Then, least absolute shrinkage and selection operator (LASSO) regression analysis was carried out by “glmnet” R package to avoid overfitting. The finally selected genes were used to build the prognostic signature and the formulae used to calculate the risk scores was as follows:
Risk score = (Coefficient mRNA1 × expression of mRNA1) + (Coefficient mRNA2 × expression of mRNA2) +··· + (Coefficient mRNAn × expression of mRNAn)[12]. HCC patients were classified into low-risk and high-risk groups according to the median value of risk scores. Kaplan-Meier (K-M) analysis was conducted using the “survival” and “survminer” R packages to compare the OS difference between low-risk and high-risk groups. Moreover, to assess the predictive performance of risk score model, a time-dependent receiver operating characteristic (ROC) curve was drawn with the “timeROC” package. In addition, the univariate and multivariate Cox regression analyses were carried out to evaluate the independent prognostic value of the signature. Besides, the principal-component analysis (PCA) based on all fatty acid metabolism-related genes and genes used to build prognostic signature model respectively for the low-risk and high-risk groups was performed by “limma” and “ggplot2” R packages. Finally, the fatty acid metabolism-related prognostic signature was further tested by K-M analysis and ROC curve in the validation set.
The nomogram model was created to predict the OS for the HCC patients on the basis of the risk score and clinical characteristics (age, gender, histologic grade, and pathological stage) using the “rms” “regplot” and “survival” packages. The calibration curves, univariate and multivariate Cox regression analysis, and ROC curve were used to evaluate the prognostic performance of the nomogram from various aspects.
Single-sample gene-set enrichment analysis (ssGSEA) was carried out to evaluate the immune-related infiltration of HCC patients by using “GSEABase” and “GSVA” R packages. Then, we calculated the enrichment scores to assess the relative level of immune-related characteristic expression using the ssGSEA algorithm and compared the difference of enrichment scores between the low-risk and high-risk score groups. In addition, Tumor Immune Dysfunction and Exclusion (TIDE) was applied to predict the feasibility of immunotherapeutic treatment response[13]. P value < 0.05 was indicated statistically significant.
In order to predict the therapeutic effect of commonly used chemotherapeutic drugs, including sorafenib, gemcitabine, 5-Fluorouracil, paclitaxel, Lapatinib, the half-maximal inhibitory concentration (IC50) of chemotherapeutic drugs was calculated by using the “pRRophetic” R package. In addition, the correlation between the risk score and IC50 of chemotherapeutic drugs was also performed.
All statistical analyses were carried out by using R software (Version 4.1.3). P value < 0.05 was considered statistically significant.
Differential expression analysis in TCGA database exhibited that there were 170 DEFAMGs, with 110 genes up-regulated and 60 down-regulated in the HCC samples (Fig. 1A-B; TableS2). Then, GO analysis showed that fatty acid metabolic process, long-chain fatty acid metabolic process, and fatty acid biosynthetic process were highly enriched GO terms (Fig. 2A-B). Meanwhile, the KEGG analysis also exhibited that fatty acid metabolism, degradation, and elongation process were highly enriched KEGG terms (Fig. 2C-D). These findings revealed that fatty acid metabolism played an important role in HCC development.
The univariate Cox proportional hazard regression analysis showed 42 genes related to OS (p < 0.01; Fig. 3A; Table S3). After the LASSO regression (Fig. 3B-C), 12 fatty acid metabolism-related genes, including PON1, CYP2C9, ACACA, ACADS, ME1, ACAT1, ELOVL1, SMS, UGDH, ADSL, HSP90AA1, and S100A10, were finally selected to build the prognostic risk score model. The risk score formula was as following: risk score = (-0.05089 × expression of PON1) + (-0.0225× expression of CYP2C9) + (0.08121 × expression of ACACA) + (-0.0002 × expression of ACADS) + (0.06192 × expression of ME1) + (-0.0144 × expression of ACAT1) + (0.08595 × expression of ELOVL1) + (0.25743 × expression of SMS) + (0.0033 × expression of UGDH) + (0.05056 × expression of ADSL) + (0.08092 × expression of HSP90AA1) + (0.03604 × expression of S100A10), which was shown in Table S4.
The risk scores of HCC patients were calculated according to the above formula. Then, the patients from TCGA database were classified into low-risk (n = 185) and high-risk (n = 185) groups by the median risk score in the training set (Table S5). Correspondingly, patients from GEO database were divided into low-risk (n = 106) and high-risk (n = 115) groups by the median risk score in the validation set (Table S6). The K-M curves showed that patients of high-risk group had a worse prognosis than that of low-risk group both in the training set (P ＜ 0.001; Fig. 4A) and validation set (P = 0.006; Fig. 4E). From the time-dependent ROC, we found the AUCs for 1-, 3-, and 5- year survival was 0.790, 0.685, and 0.698 in the training set (Fig. 4B). Correspondingly, the AUCs for 1-, 3-, and 5- year survival was 0.644, 0.608, and 0.601 in the validation (Fig. 4F). The risk score was also verified as an effective prognostic model both in the univariate (P ＜ 0.001; Fig. 4C) and multivariate Cox regression (P ＜ 0.001; Fig. 4D). Besides, the PCA analysis demonstrated that the prognostic model has a better performance than fatty acid metabolism-related genes in the distinction between the high-risk and low-risk patients (Fig. 4G-H).
Patients in the TCGA dataset were separated into several subgroups by different clinical parameters and the difference of OS was explored between the low-risk and high-risk patients. From Fig. 5A-J, K-M analysis displayed that patients with high-risk had a worse prognosis under the condition of Stage I–II (P ＜ 0.001), Stage III–IV (P ＜ 0.001), G1 + 2 (P ＜ 0.001), G3 + 4 (P = 0.01), T1 + 2 (P ＜ 0.001), T3 + 4 (P = 0.009), age ＜ 60 (P = 0.039), age ≥ 60 (P ＜ 0.001), male (P ＜ 0.001) except for female (P = 0.059).
As seen from Fig. 6, age, gender, grade, stages and risk score were selected to create the nomogram. Besides, calibration plots exhibited that the observed OS was found to be approximately consistent with the predicted OS. Univariate and multivariate Cox regression analysis showed that nomogram was effective for the prognostic assessment (P ＜ 0.001 and P = 0.001, respectively). Moreover, ROC curves illustrated that the nomogram (AUC = 0.759) had better performance in the prediction of OS than a single marker, such as risk score (AUC = 0.696), age (AUC = 0.587), Gender (AUC = 0.450), Grade (AUC = 0.539), and pathological stage (AUC = 0.663).
From Fig. 7, we found that the fractions of CD8 + T cells in the low-risk group was higher than that of high-risk group, while the fractions of Macrophages M0 was higher in the high-risk group. In addition, as for the immune-related pathways, the scores of APC co-stimulation, CCR, Check-point, HLA, MHC-class-I, parainflammation, T - cell-co-stimulation in the high-risk scores group was higher than that of low-risk group. However, the score of Type-II-IFN-Response represented higher activity in the low-risk group.
From Fig. 8, the heatmap analysis showed that a majority of metabolism pathways, such as fatty acid metabolism pathway was enriched in the low-risk score group, while the high-risk score group was connected with the cancer-related and cell signaling-related pathways (Fig. 8A). As for the TIDE, there was significant difference between the low-risk group and high-risk group, and the patients with high-risk may be more suitable for the immunotherapy. (Fig. 8B). Furthermore, patients with TP53 mutation (P ＜ 0.001) had higher risk scores (Fig. 8C). Besides, low-risk group exhibited a higher degree of progression free survival (P ＜ 0.001; Fig. 8D).
In addition, we identified the candidate drugs by the fatty acid metabolism-related risk signature. Patients from the high-risk scores group were more sensitive to sorafenib (P ＜ 0.001; Fig. 9A), gemcitabine (P ＜ 0.001; Fig. 9C), 5-Fluorouracil (P ＜ 0.001; Fig. 9E), and paclitaxel (P ＜ 0.001; Fig. 9G), and the sensitivity of chemotherapeutic drugs was negatively associated with the risk score. However, Lapatinib (P ＜ 0.001; Fig. 9I) was more suited for the low-risk groups and there was positive relation between the sensitivity of chemotherapeutic drugs and risk scores.
In this study, we analyzed the DEFAMGs between HCC patients and controls in the TCGA database. Then, the prognosis-related genes were screened and a novel prognostic signature for HCC patients was built. Next, we validated the model as an effective biomarker in the GEO database. Furthermore, we demonstrated that the risk score model could predict the OS, identify the response to immunotherapy, select the suited chemotherapy drugs, and provide the guideline for the clinical work.
In this study, a total of 12 genes (PON1, CYP2C9, ACACA, ACADS, ME1, ACAT1, ELOVL1, SMS, UGDH, ADSL, HSP90AA1, S100A10) were used to construct the prognostic signature. PON1 involved in maintaining the function of endothelial cells, inhibiting the adhesion of leukocytes, reducing the chronic inflammation, and suppressing the tumor invasion or metastasis[14]. ACADS was reported to participate in the proliferation and metastasis of HCC[15]. SMS could influence the metabolism process and involved in oncogenesis and drug resistance[16]. Overexpression of ACACA, ME1, ADSL or S100A10 promoted HCC growth, migration or invasion, and was associated with the poor prognosis[17–19]. As for the ACAT1 and CYP2C9, overexpression suppressed the proliferation and migration of tumor cells[20, 21]. As regards the ELOVL1, UGDH and HSP90AA1, studies only showed the correlation with cancer and needed further exploration to elucidate the underlying mechanisms[22–25]. As in our study, the level of PON1, CYP2C9 ACADS, ACAT1 decreased in the HCC and was positively associated with the OS of HCC patients. The expression of SMS, HSP90AA1, ADSL, UGDH, ACACA, ME1, ELOVL1, S100A10 was increased in the tumor tissues and was negatively associated with the OS.
In recent years, concerns have been raised in HCC because of the increasing incidence and mortality. However, the diagnosis and prognosis analysis of HCC patients are main on the basis of the conventional staging, which is not sensitive enough[26]. Therefore, it is necessary to identify reliable prognostic markers to improve the clinical outcomes of HCC. Recently, abnormal metabolism was reported to be closely related to the course of tumorigenesis and metastasis[27], and the altered fatty acid metabolism in cancer drew the renewed interest in particular[8]. Some studies have revealed that orchestrating fatty acid metabolism can regulate the occurrence and development of HCC[28–30]. Notably, many metabolism-related genes have been revealed to be valuable prognostic biomarkers and metabolism-related risk signatures was built to predict the OS of HCC[12, 31]. However, little study concerning fatty acid metabolism-related gene risk model has been done. Only one study selected 6 fatty acid metabolism-related genes to build a prognostic model of HCC without further drug exploration[32]. Therefore, we screened the fatty acid metabolism-related genes again and finally identified 12 genes to construct the prognostic signature more accurately. We evaluated the predictive capacity of the model and revealed that it may be used as an independent prediction factor in HCC patients. In addition, the prognostic signature showed a higher AUC value and indicated a better predictive power than other clinical characteristics.
Immune system is considered to play an essential role in preventing people from cancer. In recent years, tumor immunotherapy is considered to be the promising adjuvant therapy for HCC[33]. Tumor immunotherapy aimed at strengthening or weakening the abnormal immune state to control tumor growth or kill tumor cells[34]. However, only part of HCC patients was reported to be sensitive to the tumor immunotherapy. Wu et al. found that fatty acid metabolism could regulate the phenotype and function of immune cells, and thus affect the effect of immunotherapy[35]. In this study, according to the fatty acid metabolism-related genes model, high-risk patients responded better to immunotherapy than low-risk patients, which can guide the patient selection for immunotherapy. Besides, tumor mutational burden was regarded as a promising novel biomarker to select patients for the tumor immunotherapy in recent years[36]. TP53 mutation is the most common mutation in HCC, and it influences the progression and prognosis of HCC. Long developed a TP53-associated immune prognostic model for HCC and found the model may have important implications for prediction of OS[37]. In this study, TP53 mutation also showed higher risk scores in HCC patients. Collectively, the fatty acid metabolism-related genes model may provide the valuable information to select the HCC patients fit for immunotherapy.
It is reported that different patients respond inconsistently to chemotherapy drugs. Although sorafenib is the important first-line chemotherapeutic therapy, some HCC patients showed the sorafenib resistance and therapeutic efficacy was unsatisfying[38]. Meanwhile, the same phenomenon has been seen in other chemotherapeutic drugs[39, 40], which really deteriorated the prognosis of patients. Therefore, the individualized treatment has received great interest and may be an important approach to improve the clinical effect. Ding et al. constructed the fatty acid metabolism-related risk signatures for predicting the effect of 5-Fluorouracil in colorectal cancer and revealed patients in low-risk score group were more sensitive to 5-FU[11]. In our study, we demonstrated that HCC patients from the high-risk scores group were more suitable for sorafenib, gemcitabine, 5-Fluorouracil, and paclitaxel, while Lapatinib was more sensitive to the low-risk groups. Therefore, according to the risk score, appropriate chemotherapy drug can be chosen for the HCC patients in the future.
The main strength of this study is the construction and validation of a novel prognostic signature based on fatty acid metabolism-related genes, which has good performance in prognostic prediction and provides guidance for clinical medication. However, this study also has limitation without the experimental validation in vivo and vitro. Therefore, further exploration showed be carried out to evaluate the functions of fatty-acid metabolism-related genes in patients of HCC in the future.
In a word, we screened out fatty acid metabolism-related genes by comprehensive bioinformatic analyses in HCC patients. Then, 12 genes involved in prognosis were used to construct and validate the risk signature, which was independently associated with the OS for HCC patients. The findings provided a novel, effective and accurate prognostic predicting model for HCC patients. Moreover, the risk model could provide the guidance for personalized immunotherapy and chemotherapy.
HCC: Hepatocellular carcinoma; TCGA: The Cancer Genome Atlas; GEO: Gene Expression Omnibus; KEGG: Kyoto Encyclopedia of Genes and Genomes; MSigDB: Molecular Signature Database; DEFAMGs: Differentially expressed fatty acid metabolism-related genes; FDR: False Discovery Rate; GO: Gene Ontology; OS: Overall survival; LASSO: Least absolute shrinkage and selection operator; K-M: Kaplan-Meier; ROC: Receiver operating characteristic; PCA: Principal-component analysis; ssGSEA: Single-sample gene-set enrichment analysis; TIDE: Tumor Immune Dysfunction and Exclusion; GSVA: Gene set variation analysis.
Ethics approval and consent to participate
Not applicable.
Consent for publication
Not applicable.
Availability of data and materials
The raw data of this study are downloaded from TCGA database (https://portal.gdc.cancer.gov/) and GEO data portal (https://www.ncbi.nlm.nih.gov/
geo/), which are publicly available databases.
Competing interests
The authors declare that they have no competing interests.
Funding
This study was supported by the National Natural Science Foundation of China
(NO.81770536).
Authors' contributions
J Y and X Y designed the study, analyzed the data and wrote the article. J G analyzed the data and reviewed the manuscript. S L contributed to the concept, designed the study, and reviewed the manuscript.
All authors read and approved the manuscript.
Acknowledgements
Not applicable.