A novel fatty acid metabolism-related gene prognostic signature and candidate drugs for patients with hepatocellular carcinoma

Data acquisition

The raw RNA sequencing data of HCC were acquired from The Cancer Genome Atlas (TCGA) database on March 14, 2022. Meanwhile, the clinical information was also acquired. The GSE14520 (GPL3921) dataset, including data about gene levels and clinical features of patients with HCC, was extracted from Gene Expression Omnibus (GEO) database. Notably, samples with incomplete clinical data were excluded. Finally, 370 samples with HCC from the TCGA database were applied as the training set and 221 samples with HCC from the GEO database were used as the validation set.

Selecting FAM-related differentially expressed genes

A total of 309 FAM-related genes were extracted in FAM pathways, and then 137 common genes in TCGA and GEO cohorts were identified (Table S1). Moreover, the differentially expressed FAM-related genes (DEFAMGs) between patients with HCC and control individuals were obtained on the basis of the following criteria: False Discovery Rate <0.05 and log FC = 0.585. Then, Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses were carried out on DEFAMGs using the “cluster Profiler” R package.

Constructing and validating the FAM-related prognostic model

In the training set, the genes significantly related to patient overall survival (OS) were selected from DEFAMGs using the Univariate Cox analysis. P value <0.01 was regarded as significant. To avoid overfitting, least absolute shrinkage and selection operator (LASSO) regression analysis was adopted to construct a prognostic model (Liang et al., 2020). With the “glmnet” R package, the LASSO algorithm was utilized for variable selection and shrinkage. The independent variable in the regression was the normalized expression matrix of candidate prognostic DEFAMGs, and the response variables in the training set were OS and patients’ status. According to each gene’s normalized expression and related regression coefficients, the risk score of each patient was computed. The formulae for these calculations have been previously described (Zhao et al., 2021). Then, Kaplan–Meier (K-M) analysis and generation of receiver operating characteristic (ROC) curve were carried out. Next, the independent prognostic performance of the model was assessed by regression analysis. Besides, the principal-component analysis on the basis of FAM-related genes and genes used to build prognostic signature model was performed using the “limma” and “ggplot2” R packages. Finally, the FAM-related prognostic signature was further tested using the K-M analysis and ROC curve in the validation set (GSE14520).

Establishing a predictive nomogram

The nomogram model was created to forecast the OS for patients with HCC via the “rms”, “regplot”, and “survival” packages. Then, the prognostic performance of the nomogram was assessed from various aspects.

Exploration of the model in immunity therapy

The gene sets, which included numerous different types of human immune cell subtypes and immune-related activities, were gathered for evaluation of immune-related characteristics as previously described (Ding et al., 2021). Single-sample gene-set enrichment analysis (ssGSEA) was performed using the “GSEABase” and “GSVA” R packages to observe the immune-related infiltration in patients with HCC. To explore the differences of biological processes between low- and high-risk score groups, GSVA was performed on the gene profile using the “GSVA” R package. The reference gene sets were “c2.cp.kegg.v7.4.symbols” from the Molecular Signatures Database (https://www.gsea-msigdb.org/gsea/msigdb).In addition, Tumor Immune Dysfunction and Exclusion (TIDE) was applied to forecast the feasibility of immunotherapy (Xu et al., 2020). Moreover, the Wilcoxon test was used to test intergroup differences in the levels of potential immune checkpoint molecules, such as PD-1 and CTLA-4. Statistical significance was set at FDR >0.05.

Exploration of chemotherapy response based on risk model

In order to predict the therapeutic effect of commonly used chemotherapeutic drugs, including sorafenib, gemcitabine, 5-fluorouracil, paclitaxel, and lapatinib, the half-maximal inhibitory concentration (IC50) of chemotherapeutic agents was detected using the “pRRophetic” R package. The relationship between the risk score and IC50 of chemotherapeutic drugs was also assessed.

Cell culture and real-time polymerase chain reaction

The human hepatocyte cell line MIHA and HCC cell lines, MHCC-97H and Huh-7, were obtained from the BeNa Culture Collection (Suzhou, China), and cultured in Dulbecco’s modified Eagle medium with 10% fetal bovine serum in a 37 °C humidified incubator containing 5% CO2. Total RNA was isolated from cells using the TRI-ZOL reagent (Takara, Otsu, Japan). Then, the cDNA synthesis was performed using the PrimeScript™ RT Master Mix Kit (Takara). Real-time polymerase chain reaction (RT-PCR) was carried out using the SYBR-Green PCR master mix (Takara) in the Roche Light Cycler R480 (Roche, Basel, Switzerland). The Glyceraldehyde-3-phosphate dehydrogenase (GAPDH) was used as the internal reference. The relative expression of mRNA levels was calculated by 2^−ΔΔCt and the primer sequences are listed in Table S2.

Western blot analysis

Proteins were extracted from cells with RIPA lysis buffer (Beyotime, Shanghai, China) containing phenylmethyl sulfonyl fluoride. The concentration of protein was measured by bicinchoninic acid protein assay kit (Beyotime). Then, a suitable quality of protein samples was subjected to sodium dodecyl sulfate polyacrylamide gel electrophoresis and transferred to PVDF membranes (Millipore, Billerica, MA, USA). The 8% skimmed milk was used to block the PVDF membranes at room temperature for 1 h, and the PVDF membranes were incubated with primary antibodies at 4 °C overnight. The corresponding secondary antibodies was used at room temperature for 1 h and the protein signals were detected with FluorChem FC3 system (ProteinSimple, Santa Clara, CA, USA) using an enhanced chemilusystem reagent (Thermo Fisher). Image J was used to measure the gray values of the target and reference bands. Then, the target/ reference ratio was regarded as the relative expression level of target. The antibodies used in this study were as follows: PON1 (1:1000, A3441, Abclonal, Wuhan, China), CYP2C9 (1:4000, GB114267, Servicebio, Wuhan, China), ACACA (1:1000, A19627, Abclonal), ACADS (1:1000, A0945, Abclonal), ME1 (1:1000, A3956, Abclonal), ACAT1 (1:1000, A20943, Abclonal), ELOVL1 (1:1000, A07612, Boster, Wuhan, China), SMS (1:1000, A9380, Abclonal), ADSL (1:1000, A6278, Abclonal), UGDH (1:1000, A1210, Abclonal), HSP90AA1 (1:1000, A5006, Abclonal), S100A10 (1:1000, A13614, Abclonal), GAPDH (Abclonal), horse radish peroxidase (HRP)-linked goat anti rabbit IgG, and HRP-linked goat anti mouse IgG (1:4000, Antgene, Wuhan, China).

Statistical analysis

The statistical analyses were carried out via R version 4.1.3 (R Core Team, 2021) and SPSS 25.0 (SPSS, Inc., Chicago, IL, USA). The ssGSEA scores of immune cells or pathways between the high- and low-risk groups were compared using the Mann–Whitney test with P values corrected using the BH technique. K-M analysis was used to evaluate the OS differences between different groups. Univariate and multivariate Cox regression analyses were performed to identify independent predictors of OS. The Chi-squared test was used to compare proportional differences. Gene and protein levels in different cell lines were analyzed by independent sample t-test or correction t-test depending on the homogeneity of variance analysis. Statistical significance was set at P < 0.05.

Enrichment analysis of DEFAMGs

Differential expression analysis using data from the TCGA database exhibited that there were 170 DEFAMGs, with 110 genes and 60 genes being upregulated and downregulated, respectively in HCC samples (Figs. 1A–1B; Table S3). Then, GO analysis showed that FA metabolic process, long-chain FA metabolic process, and FA biosynthetic process were highly enriched (Figs. S1A–S1B). Meanwhile, the KEGG analysis also exhibited that FAM, degradation, and elongation process were highly concentrated (Figs. S1C–S1D). These findings revealed that FAM played an important part in HCC.

Building a FAM-related risk signature

There were 42 genes associated with OS (P < 0.01; Fig. 2A; Table S4). After the LASSO regression (Figs. 2B–2C), the following 12 FAM-related genes, PON1, CYP2C9, ACACA, ACADS, ME1, ACAT1, ELOVL1, SMS, UGDH, ADSL, HSP90AA1, and S100A10, were selected to build the prediction model. The formula used was as follows: risk score = (−0.05089 × PON1_expression) + (−0.0225 × CYP2C9_expression) + (0.08121 × ACACA_expression) + (−0.0002 × ACADS_expression) + (0.06192 × ME1_expression) + (−0.0144 ×ACAT1_expression) + (0.08595 × ELOVL1_expression) + (0.25743 × SMS_expression) + (0.0033 × UGDH_expression) + (0.05056 × ADSL_expression) + (0.08092 × HSP90AA1_expression) + (0.03604 × S100A10_expression), which is shown in Table S5.

Validating the FAM-related risk signature

The risk scores of HCC were computed using the above formula. Using the median risk score, patient samples from the TCGA database were classified into low-risk (n = 185) and high-risk (n = 185) sets (Table S6). Correspondingly, patients from the GEO database (GSE14520; validation set) were divided into low-risk (n = 106) and high-risk (n = 115) groups (Table S7). Patients with higher risk had a worse prognosis both in the training (P < 0.001; Fig. 3A) and validation (GSE14520) sets (P = 0.006; Fig. 3E). In the training set, the AUCs for 1-, 3-, and 5-year survival were 0.790, 0.685, and 0.698, respectively (Fig. 3B), while in the validation set, they were 0.644, 0.608, and 0.601, respectively (GSE14520) (Fig. 3F). The risk score was also verified to be an effective prediction model both in the univariate (P < 0.001; Fig. 3C) and multivariate Cox regression (P < 0.001; Fig. 3D). Besides, the principal-component analysis displayed that the model had a better performance than FAM-related genes in distinction of the different-risk patients (Figs. 3G–3H).

Stratification analysis of the prediction model

Patients in the TCGA dataset were separated into several subgroups by different clinical parameters and the OS difference was then explored. We exhibited high-risk patients 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) (Figs. 4A–4J).

Construction and assessment of a predictive nomogram

As shown in Fig. 5A, the nomogram was created via clinical information and risk score. Besides, calibration plots exhibited that the observed OS was approximately consistent with the predicted OS (Fig. 5B). Moreover, the nomogram (AUC = 0.759) was illustrated to have 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) (Fig. 5C). Then, Figs. 5D and 5E displayed that the nomogram was effective for the prognostic evaluation (P < 0.001 and P = 0.001, respectively).

Exploring the immune-related characteristic

Figures 6A and 6B showed that the fractions of CD8+ T cells were higher in the low-risk set, while in the high-risk set, higher fractions were observed in macrophages M0. In addition, as for the immune-related pathways, the scores of APC co-stimulation, CCR, Check-point, HLA, MHC-class-I, para-inflammation, T-cell-co-stimulation in the high-risk scores set was higher than those of the low-risk set. However, the score of Type-II-IFN-Response represented higher activity in the low-risk set.

Role of the FAM model in predicting the response to immunotherapy

The heatmap analysis showed that a majority of metabolic pathways, such as FAM pathway, were intensive in the low-risk set, while the high-risk was associated with the cancer-related and cell signaling-related pathways (Fig. 7A). The TIDE showed higher risk patients may be more suitable for immunotherapy (Fig. 7B). Furthermore, the levels of PD-1 (P < 0.001, Fig. 7C) and CTLA-4 (P < 0.001, Fig. 7D) were higher in high- than in low-risk patients.

Response to chemotherapy

We identified the candidate drugs by the FAM-related risk signature (Fig. 8). Patients from the high-risk set were more sensitive to sorafenib (P < 0.001), gemcitabine (P < 0.001), 5-fluorouracil (P < 0.001), paclitaxel (P < 0.001), and the sensitivity to these chemotherapeutic drugs was negatively connected with the risk score. However, lapatinib (P < 0.001) was more suitable for the low-risk patients and there was a positive relationship between the sensitivity of chemotherapeutic drugs and risk scores.

Validating the genes in cells by RT-PCR and western blot

The mRNA levels of 12 genes between hepatocyte cell line and HCC cell lines were evaluated (Fig. 9). In comparison to MIHA cell line, the expressions of PON1, CYP2C9, ACADS, and ACAT1 were significantly lower in the MHCC97H and Huh-7 cell lines. In addition, the levels of ME1, HSP90AA1, S100A10, ACACA, ELOVL1, SMS, and ADSL were higher in the MHCC97H cell line, and the levels of ME1, HSP90AA1, S100A10, SMS, and UGDH were higher in the Huh-7 cell line compared with MIHA cell line.

In Fig. 10, we assessed the protein levels of the risk-score genes in these cell lines by western blot. The levels of ACACA, ME1, ELOVL1, SMS, UGDH, ADSL, HSP90AA1, and S100A10 were significantly higher in MHCC97H cells than in MIHA cells, and the levels of ME1, ELOVL1, SMS, UGDH, HSP90AA1, and S100A10 were higher in Huh-7 cells than in MIHA cells. In addition, the levels of PON1, CYP2C9, ACADS, and ACAT1 were lower in MHCC97H cells than in MIHA cells, and the levels of PON1, CYP2C9, and ACADS were significantly lower in Huh-7 cells than in MIHA cells.