Molecular subtypes and scoring tools related to Foxo signaling pathway for assessing hepatocellular carcinoma prognosis and treatment responsiveness

Background: Transcription factors in Foxo signaling pathway influence hepatocellular carcinoma metastasis through epithelial mesenchymal transition-related pathways. Prognostic factors in the Foxo signaling pathway are feasible for HCC prognosis and therapeutic management.

Methods: Based on the differentially expressed genes and Foxo signaling pathway genes in HCC, the ConsensusClusterPlus package was conducted to identify Foxo signaling pathway-related molecular subtypes in HCC. Based on the DEGs in the FMSs, the optimal prognostic factors in HCC were screened by cox and least absolute shrinkage and selection operator (LASSO) cox analysis to form the Foxo prognosis score (FPS). The prognostic predictive effectiveness of FPS was assessed by Kaplan Meier (K-M) analysis and Receiver Operating Characteristic (ROC) analysis. Additionally, tumor microenvironment (TME) score, tumor mutation burden (TMB) and treatment sensitivity differences in FMSs and FPS groups were also evaluated.

Results: There were low, medium and high Foxo signaling pathway activity molecular subtypes in HCC named FMS 1, FMS 2 and FMS 3, respectively. FMS 1 with lowest Foxo signaling pathway activity presented an excellent survival advantage, while FMS 3 with highest Foxo signaling pathway activity exhibited an inhibitory TME status. According to FPS grouping, low FPS exhibited favorable survival, low TMB and anti-tumor activity. Patients in the low FPS group were mostly in the early stage of cancer. Moreover, we found that patients with high and low FPS exhibited different sensitivity to chemotherapy, and patients with low FPS were more sensitive to immunotherapy.

Conclusion: We revealed a novel molecular subtype and prognostic tool based on Foxo signaling pathway signature, which could potentially provide a direction for accurate and effective assessment of potential personalized treatment options and prognostic management for HCC patients.

Introduction

Hepatocellular carcinoma (HCC) is the most prevalent pathological type of liver cancer, with 90% of liver cancer cases being HCC (Llovet et al., 2021). According to the prognostic model of cancer pathology by the World Health Organization (WHO), more than 1 million HCC deaths were expected (Llovet et al., 2021). In current clinical practices, liver resection, transplantation, and chemotherapy remained the dominant options for the treatment of HCC. Surgical resection demonstrated good five-year survival in the treatment of HCC, with five-year survival rates of approximately 70%–80%, however, there were postoperative liver dysfunction and recurrence difficulties, maintaining liver function in HCC patients after surgery was reported to be the key challenge (European Association for the Study of the Liver, 2018). Liver transplantation was accepted as the superior HCC treatment option to surgical resection in clinical practice, with more than half of HCC liver transplant patients achieving postoperative survival of 10 years or more, however, the waiting time for donor in liver transplantation was an uncertainty (Llovet et al., 2021). In contrast, chemotherapy was predominantly compromised by tumor heterogeneity, with different individuals exhibiting diverse treatment progression (Vogel et al., 2022). Therefore, postoperative prognostic management of HCC was the substantial challenge in clinical practice, and accurate prognostic assessment tools were integral to improve the survival rate of HCC.

The transcription factors in the Foxo signaling pathway, FOXO1 and FOXO3, contributed critical proteins to the development of malignant progression of HCC (Yang et al., 2021). In HCC, Akt signaling normally occurs in an active state, with overactive AKT signaling inhibiting FOXO1 transcriptional processes. Regular FOXO1 transcription inhibited epithelial mesenchymal transition (EMT) and transforming growth factor-β (TGF-β) expression, whereas dysregulation of EMT and TGF-β expression promoted HCC cell invasion and metastasis to other tissues when the FOXO1 transcription process was blocked (Dong et al., 2017). FOXO3 exhibited normal expression levels in normal liver tissues but was abnormally highly expressed in liver tissues of HCC patients, with further studies also confirming that FOXO3 contributed to the suboptimal disease-free survival and prognosis of HCC (Ahn et al., 2018; Song et al., 2020).

In this study, we integrated high-throughput sequencing data of HCC from multiple databases (The Cancer Genome Atlas, GENE EXPRESSION OMNIBUS) to conduct a comprehensive theoretical analysis. The intention of this study was aimed to establish molecular subtypes of Foxo signaling pathway activity in HCC, followed by an attempt to construct the prognostic scoring system for HCC based on the specialized genes in the subtypes. To our knowledge, this is the first time that Foxo signaling pathway was studied in the prognosis of liver cancer.

Materials and methods

Data collection

RNA-seq data for the TCGA-LIHC sequencing project were sourced from the Cancer Genome Atlas (TCGA, https://portal.gdc.cancer.gov/) database. This study was conducted using the HCC samples cohort in TCGA as the training set. To validate the results of this investigation, three additional datasets were accessed as validation sets. RNA-Seq data for HCC sequencing researches were accessed from HCCDB (Lian et al., 2018) (http://lifeome.net/database/hccdb/home.html) (project name: ICGC-LIRI-JP) and GENE EXPRESSION OMNIBUS (GEO) database (https://www.ncbi.nlm.nih.gov/gds) (registration number: GSE14520 (Li et al., 2022), GSE76427 (Grinchuk et al., 2018)), respectively. Additionally, clinical information of samples in TCGA-LIHC, ICGC-LIRI-JP, GSE14520 and GSE76427 were also collected. The genes in the Foxo signaling pathway were sourced from the Kyoto encyclopedia of genes and genomes (KEGG) database (https://www.genome.jp/kegg/) (Supplementary Table S1).

Data pre-processing

Clinical data and RNA-Seq data of samples in TCGA-LIHC, ICGC-LIRI-JP, GSE14520 and GSE76427 were imported into the sangerbox database (Shen et al., 2022) for the following processing. Samples with missing clinical data were excluded; samples with survival time >0 were retained; ensemb information was converted to Gene symbol; gene with multiple probe information were averaged as the expression data of the gene. After processing, 355 HCC samples and 50 control samples were included in TCGA-LIHC; 212 HCC samples were maintained in ICGC-LIRI-JP; 221 and 115 HCC samples were maintained in GSE14520, GSE76427. Detailed information was listed in Table 1.

TABLE 1

TABLE 1. Clinical information of HCC samples in TCGA-LIHC, ICGC-LIRI-JP, GSE14520 and GSE76427.

Identification of Foxo Molecular subtypes

In TCGA-LIHC, the limma package (Ritchie et al., 2015) was performed to screen differentially expressed genes (DEGs) in HCC samples based on the expression matrix of HCC samples and control samples (|log2FC|>1 & FDR<0.05). The expression matrix of genes in the Foxo signaling pathway was extracted and the univariate COX model was performed to identify HCC prognosis-related genes (p < 0.05). Overlapping genes in prognosis-related genes in DEGs and Foxo signaling pathway were extracted, and the consistency clustering analysis was performed based on the expression matrix of overlapping genes. The ConsensusClusterPlus R package (Wilkerson and Hayes, 2010) was performed to execute the consistency clustering analysis. The clustering parameters were set as follows: metric distance: km algorithm and euclidean; number of bootstraps: 500; number of clusters range: k = 2–10. Foxo Molecular subtypes (FMS) in HCC were determined according to the consistency matrix and consistency cumulative distribution function for each k-value range.

Protein-protein interaction network analysis

The limma package was used to perform differential analysis to identify DEGs in each FMS according to the gene expression matrix in the different FMSs. The DEGs were imported into the STRING database (Szklarczyk et al., 2021) to construct the Protein-protein interaction network (PPI network) with parameters set to: confidence score > 0.7. The PPI network was imported into Cytoscape software (version: 3.9.1) (Shannon et al., 2003) for MCODE sub-network clustering and topology analysis. Genes in the sub-networks with the highest scores based on MCODE algorithm were included for subsequent analysis.

Construction of Foxo prognosis score

Based on the gene expression matrix in MCODE 1, the coxph function in the survival package (Therneau and Lumley, 2015) was performed for univariate COX model analysis to initially screen for HCC prognosis-related genes. The Least absolute shrinkage and selection operator COX model was performed to reduce the model fit, and the penalty parameter lambda was selected by 10-fold cross-validation method. The model under the best lambda value was selected for multivariate COX model construction, and the step AIC function in the MASS package was performed for optimal model selection, and the model under the minimum value of Akaike information criterion (AIC) was considered the optimal model. The Foxo prognosis score (FPS) was constructed based on the gene regression coefficients (β) and expression data in the model. The evaluation equation was:

$\mathrm{F}\mathrm{P}\mathrm{S}=\sum {\mathrm{\beta }}_i\ast {expression}_i$

Prognostic guidance value of FPS

The FPS of samples in TCGA-LIHC was calculated according to the FPS evaluation equation. Samples with FPS>0 was defined as high FPS group and samples with FPS<0 was defined as low FPS group. The survival package was performed to perform Kaplan-Meier (K-M) survival analysis and graphed K-M survival curves. The timeROC package (Blanche et al., 2013) was performed for Receiver Operating Characteristic (ROC) analysis and ROC curves were graphed. The prognostic guidance value of FPS was validated in the external validation cohorts, CGC-LIRI-JP, GSE14520 and GSE76427.

Mutational landscape analysis of genes in the Foxo signaling pathway

The somatic mutation data corresponding to HCC samples in TCGA-LIHC were also sourced from the TCGA database. The Maftools package (version: 2.8.05) (Mayakonda et al., 2018) was deployed to perform copy number variation (CNV) and Single nucleotide variant (SNV) analysis. The tmb function in the Maftools package was performed to assess the level of Tumor mutation burden (TMB) in HCC. Based on Thorsson et al. (2018), molecular mutation characteristics of TCGA-LIHC were captured to assess the Aneuploidy Score, Fraction Altered, Number of Segments, and Number of Segments levels in HCC samples.

Scoring or abundance analysis of infiltrating immune cells in the immune microenvironment

For HCC samples in TCGA-LIHC, multiple tumor microenvironment (TME) infiltration algorithms were conducted to assess the immune score or infiltration abundance of infiltrating immune cells in these samples. The CIBERSORT algorithm (Chen et al., 2018) was developed to assess the infiltrative abundance of 22 immune cells in HCC samples; the MCP-Count (Becht et al., 2016), and TIMER algorithms (Li et al., 2020) were developed to assess the immune scores of immune cells in TME. The GSVA package (Hanzelmann et al., 2013) was performed to perform single-sample gene set enrichment analysis based on the 28 immune cell gene sets in the study by Barbie et al. (2009).

Pathway analysis

The h.all.v7.5.1.symbols.gmt gene set from MSigDB database was imported from the GSEA website for single sample gene set enrichment analysis (ssGSEA) via the GSVA package. The Hmisc package (Harrell and Harrell, 2019) was employed for spearman correlation analysis between the pathway ssGSEA score and FPS.

Evaluation of immunotherapy/chemotherapy sensitivity

The present study also sought to reveal the sensitivity of HCC patients in different FMS and FPS groups to immunotherapy as well as chemotherapeutic drug treatment. Immune checkpoint gene expression levels in HCC samples were assessed according to the immune checkpoint gene in the study of Auslander et al. (2018). Next, the Tumor Immune Dysfunction and Exclusion (TIDE) score, interferon gamma (IFNG), exclusion score, dysfunction score, and myeloid-derived suppressor cells (MDSCs) score of HCC samples were accessed in the TIDE (http://tide.dfci.harvard.edu/login/) (Jiang et al., 2018) website. Chemotherapeutic drugs, Erlotinib, Saracatinib, TGX221, Roscovitine, GNF-2, CGP-082996, Pyrimethamine, NSC-87877, treated HCC Sequencing data were sourced from the Genomics of Drug Sensitivity in Cancer (GDSC, https://www.cancerrxgene.org/) database. pRRophetic package (Geeleher et al., 2014) was developed to assess the half maximal inhibitory concentration (IC50) values of these drug-treated HCC data.

Assessment of the clinical value of FPS and construction of nomogram

Age, Gender, Stage, Grade, and FPS information of HCC samples in TCGA-LIHC were extracted and univariate COX and multivariate COX analyses were performed to identify pivotal independent prognostic parameters in HCC for Nomogram construction. Meanwhile, the calibration curves were performed to assess the prognostic predictive value of Nomogram based on the 1-year, 3-year, and 5-year survival rates of HCC predicted by Nomogram versus the actual survival rates of recorded HCC. Additionally, Decision curve (DCA) was graphed for assessing the clinical predictive value of FPS, Nomogram based on its information.

Statistical analysis

All statistical analyses were performed using R software (version 4.0.3). Moreover, Sangerbox (http://sangerbox.com/home.html) also assisted in data processing. p value < 0.05 was treated as statistically significant.

Results

Genomic alterations of prognosis-related Foxo signaling pathway genes in HCC

In TCGA-LIHC, 2751 DEGs were screened out between HCC and normal tissues and prognosis-related Foxo signaling pathway genes were screened, with 19 genes co-existing in them (Figure 1A). Among these genes, 5 genes expressed highly in control tissues and 14 genes expressed highly in tumor tissues (wilcox.test, p < 0.05) (Figure 1B). Subsequently, the mutation status of 19 genes in HCC tissues was analyzed with 23 HCC samples (6.32%) in which these genes were mutated. Missense Mutation and Nonsense_Mutation constituted the most frequent type of mutations (Figure 1C). 19-genes were also estimated for CNV, for which we found genes exhibited lower CNV amplification or deletion (Figure 1D). Finally, according to CNV status, HCC samples were divided into CNV amplification, CNV deletion and CNV diploid groups, and 19-gene expression was evaluated in all three groups. We found that these genes exhibited higher expression levels in the CNV amplification group overall (Figure 1E).

FIGURE 1

FIGURE 1. Prognosis-related genomic alterations of Foxo signaling pathway genes in HCC samples. (A) Wayne diagram showing overlapping genes of prognosis-related Foxo signaling pathway genes and DEGs in HCC. (B) 19 - Gene expression in tumor tissue and normal tissue. (C) Mutational landscape of 19 genes in tumor samples. (D) 19-gene CNV landscape in HCC samples. (E) 19-Gene expression levels in the CNV amplification group, CNV deletion group, and CNV diploid group. *p < 0.05, **p < 0.01, ****p < 0.0001, ns, no difference.

Foxo Molecular subtypes in HCC

Following the 19-genes expression matrix, HCC samples in TCGA-LIHC were performed consistency clustering analysis to uncover the molecular subtypes concerning Foxo signaling pathway signature (FMS). When k = 3, the clustering of HCC patients exhibited excellent consistency with low interference between samples (Supplementary Figure S1A–C). Therefore, we defined 3 molecular subtypes in HCC named FMS 1, FMS 2, and FMS 3. The molecular subtypes were further verified in the samples of GSE14520, and the trends were consistent (Supplementary Figure S1D–F). In TCGA-LIHC and GSE14520, FMS 1 both exhibited an excellent survival benefit (Figures 2A, B). In TCGA-LIHC, from the statistical information of clinical information and 19-gene expression heatmap in FMS 1-3, we could clearly observe that FMS 1 exhibited remarkably high expression of protective factors and remarkably low expression of risk factors. FMS 3 exhibited the contrast tendency. Therefore, FMS 1, FMS 2, and FMS 3 were defined as the low activity cluster, medium activity cluster, and high activity cluster of Foxo signaling pathway, respectively. Moreover, we also noted that patients with HCC in FMS 3 exhibited high clinical stage features (Figure 2C).

FIGURE 2

FIGURE 2. Differences in survival, clinical characteristics among FMSs. (A–B) K-M Ssurvival curve of FMSs in TCGA-LIHC, GSE14520. (C) Heatmap showing 19-gene expression and clinicopathological information in FMSs. *p < 0.05, **p < 0.01, ***p < 0.001, ****p < 0.0001.

TME differences in FMSs

According to CIBERSORT results, B cells naïve, T cells CD4 memory resting, Monocytes, Macrophages M1, Mast cells resting expressed high infiltration abundance in FMS 1; T cells regulatory (Tregs), Macrophages M0 expressed high infiltration abundance in FMS 3 (Figure 3A). Depending on the ssGSEA, MCP-Count, and TIMER results, we observed that higher levels of immune cell infiltration scores were exhibited in FMS 3 (Supplementary Figure S2). Additionally, lower levels of immune checkpoint gene expression were measured in FMS 1, and the highest level of immune checkpoint gene expression was demonstrated in FMS 3 (Figure 3B). We further remarked that TIDE score, IFNG score, Exclusion score, MDSC score appeared to be higher in FMS 3 (Figure 3C).

FIGURE 3

FIGURE 3. TME differences between FMSs. (A) CIBERSORT results in FMSs, demonstrating the abundance of 22 immune cell infiltrates. (B) Immune checkpoint gene expression levels in FMSs. (C) TIDE score, IFNG score, Exclusion score, Dysfunction score, MDSC score. *p < 0.05, **p < 0.01, ***p < 0.001, ****p < 0.0001, ns, no difference.

Protein-protein interaction network in HCC

To explore the prognostic differences and TME activity differences, differential analysis was performed in FMS 1-3 to construct PPI networks to explore the protein regulatory network differences among FMSs. Firstly, differential analysis was performed. 576 DEGs were presented in FMS 1, consisting of 188 upregulated expression DEGs and 388 downregulated expression DEGs; 2 upregulated expression DEGs were presented in FMS 2; 1,048 DEGs were presented in FMS 3, consisting of 658 upregulated expression DEGs and 390 downregulated expression DEGs. 1,108 overlapping DEGs were presented in FMSs, and these genes were included in STRING to construct PPI network, and 717 genes were included in PPI network with confidence score >0.7 (Figure 4A). According to the MCODE algorithm in cytoscape software, the sub-network MCODE 1 with the highest large MCODE degree value in the PPI network was retained for functional enrichment analysis, containing 103 genes. These genes were remarkably enriched in cell cycle, DNA replication, oocyte meiosis, progesterone-mediated oocyte maturation, and p53 signaling pathway (Figure 4B). The GO results further evidenced that these genes were closely interlinked with DNA replication transcriptional processes, such as, regulation of mitotic cell cycle, cell division, cell cycle phase transition, and mitotic cell cycle phase transition (Figure 4C).

FIGURE 4

FIGURE 4. PPI network and functional analysis. (A) PPI network for DEGs in FMSs (The darker the color of a single gene, the more it interacts with other genes, indicating its higher importance). (B) KEGG chord diagram of genes in MCODE 1. (C) GO bubble map of genes in MCODE 1.

FPS for assessing HCC prognosis

In TCGA-LIHC, 92 prognosis-related genes were identified by univariate COX analysis, we found that these genes were all Risk genes. According to the trajectory diagram of the change of independent variable coefficients in LASSO COX analysis, the interference of similar genes in the model was minimized at lambda = 0.0172, and the model fit was minimized at this time, when 16 genes were included in the model (Supplementary Figure S3A, B). Finally, CENPA, CDKN3, KPNA2, ARHGAP11A, KIF18A, ASF1B, HMMR, CDCA8, and CCNB2 were screened out by multivariate COX analysis and were considered to be the optimal model to compose the HCC prognostic scoring system with FPS = 0.875*CENPA+(-0.291*CDKN3)+0.414*KPNA2 +(-1.298*ARHGAP11A)+0.694*KIF18A+(-0.496*ASF1B)+ 0.292*HMMR+0.899*CDCA8+(-0.554*CCNB2) (Figure 5A). The prognostic differences between the high FPS group (FPS>0) and the low FPS group (FPS<0) according to the FPS = 0 grouping and the prognostic accuracy of FPS were assessed by K-M survival curves and ROC curves, respectively. Patients with HCC in the low FPS group presented remarkable survival advantage (Figure 5B). The AUC values of FPS predicting 1-year, 2-year, 3-year, 4-year, and 5-year survival of HCC were 0.82, 0.77, 0.77, 0.79, and 0.8 (Figure 5C). The prognostic value of FPS was further validated in ICGC-LIRI-JP, GSE76427, and GSE14520, we observed that patients with low FPS in all three datasets exhibited remarkable survival advantage with FPS demonstrated higher AUC values (Figures 5D–G, Supplementary Figure C, D).

FIGURE 5

FIGURE 5. Predictive value of FPS. (A) Multivariate COX forest plot of 9-prognostic factors. (B, C) K-M survival curve, ROC curve in TCGA-LIHC. (D, E) K-M survival curve, ROC curve in ICGC-LIRI-JP. (F, G) K-M survival curve, ROC curve in GSE76427. *p < 0.05, **p < 0.01, ***p < 0.001.

Prognostic performance of FPS in pathological subgroups

In TCGA-LUAD, the clinicopathological statistical information, 9-gene expression level, and FPS distribution of HCC samples were demonstrated in Figure 6A. We observed that there were remarkable differences in Stage, T stage, and Status of patients in the high FPS group and low FPS group, and the 9-gene expressed higher levels in the high FPS group. In subgroups including Stage, FPS in Age, T stage, Grade, and FMS Clusters showed elevated trends with advanced pathological stage. Also, there is also a slight trendy of differences in FPS between men and women (Figure 6B). In Stage subgroups (Stage I + II, Stage III + IV), Grade subgroups (G1+G2, G3+G4), Age subgroups (Age<=60, Age>60), and Gender subgroups (Female, Male), the low FPS group all exhibited remarkable survival advantage (Figures 6C–F).

FIGURE 6

FIGURE 6. Clinical characteristics in the FPS groups. (A) Heat map of 9-prognostic factor expression in FPS groups combined with clinical characteristics of patients. (B) FPS differences in clinicopathological groups. (C–F) K-M survival curves of patients in high and low FPS groups in Stage I + II, Stage III + IV, G1+G2, G3+G4, Age<=60, Age>60, FEMALE, and MALE groups. *p < 0.05, ****p < 0.0001.

Characterization of pathway and genomic variants in FPS groups

To observe the connection between FPS and biological functions, the spearman correlation between ssGSEA scores and FPS was evaluated for each pathway in the h.all.v7.5.1.symbols.gmt gene set. We found that FPS exhibited remarkably negative correlations with metabolism-related pathways and remarkably positive correlations with cell cycle-related pathways (Figure 7A). Further, we found differences in gene mutations between high and low FPS groups. The top five mutated genes in the high FPS group were TP53, FRAS1, NBEA, SETD2, and TG; the top five mutated genes in the low FPS group were TP53, MUC4, SPTA1, BAP1, and DYNC2H1. TP53 was more frequently mutated in the high FPS group (Figure 7B). Whereas, the high FPS group exhibited higher TMB (Figure 7C). Patients in the low-TMB group also exhibited remarkable survival advantage (Figure 7D). HCC patients with low TMB and low FPS developed remarkable survival advantage compared with those with high TMB and high FPS scores (Figure 7E). Moreover, we found that the high FPS group exhibited higher Aneuploidy Score, Fraction Altered, Number of Segments, and Number of Segments compared to the low FPS group (Figure 7F).

FIGURE 7

FIGURE 7. Genomic mutation statistics and pathway differences in FPS groups. (A) Biological pathways significantly associated with FPS. (B) Mutation landscape in high and low FPS groups. (C) TMB in high and low groups. (D) K-M survival curves of TMB groups. (E) K-M survival curves in patients with TMB combined with FPS groups. (F) Homologous Recombination Defects, Aneuploidy Score, Fraction Altered, Number of Segments in high and low groups. ***p < 0.001, ****p < 0.0001.

Differences in treatment sensitivity and TME in the FPS groups population

The immune cell infiltration characteristics in HCC samples in FPS grouping were then characterized. According to the CIBERSORT results, T cells regulatory (Tregs), Macrophages M0 exhibited higher infiltration abundance in the high FPS; B cells naïve, T cells CD4 memory resting, Mast cells resting exhibited higher infiltration abundance in the low FPS group (Figure 8A). FPS exhibited remarkable positive correlation with immune checkpoint gene expression and 9-gene expression levels (Figure 8B). According to the immune cell immune score in ssGSEA, MCP-Count, and TIMER algorithms, FPS exhibited remarkable positive correlation with most immune cell immune scores (Figure 8C). Moreover, we found that the TIDE score, IFNG score, Exclusion score, Dysfunction score, and MDSC score were remarkably higher in the high FPS group than in the low FPS group (Figure 8D). In addition, the sensitivity of patients with different FPS to chemotherapeutic agents according to its characteristics was estimated. Notably, patients in the low FPS group were more sensitive to treatment with Erlotinib, Saracatinib, TGX221, and Roscovitine; patients in the high FPS group were more sensitive to treatment with GNF-2, CGP-082996, Pyrimethamine, NSC -87877 were more sensitive (Figure 8E).

FIGURE 8

FIGURE 8. TME differences and treatment sensitivity in FPS groups. (A) CIBERSORT results in FPS groups, demonstrating the abundance of 22 immune cell infiltrates. (B) Correlation between FPS and immune checkpoint expression, 9-prognostic factor expression. (C) Correlation of FPS with ssGSEA, MCP-Count, and TIMER immune scores. (D) TIDE score, FNG score, Exclusion score, Dysfunction score, MDSC score in low and high FPS groups. (E) IC50 of Erlotinib, Saracatinib, TGX221, Roscovitine, GNF−2, CGP−082996, Pyrimethamine, NSC−87877 in low and high groups. *p < 0.05, **p < 0.01, ***p < 0.001, ****p < 0.0001, ns, no difference.

Nomogram plots predicting HCC survival

Univariate COX and multivariate COX analyses were performed in TCGA-LIHC by integrating information on Age, Gender, Stage, Grade, and FPS. We identified FPS and Stage as independent clinical prognostic factors for HCC (Figures 9A, B). Therefore, with FPS and Stage information, we constructed nomogram for HCC prognosis (Figure 9C), and the calibration curve indicated that Nomogram predicted 1-year, 3-year, and 5-year survival and actual observations of HCC with high fitting tendency (Figure 9D). According to DCA, we observed that nomogram and FPS possessed excellent clinical observation value (Figure 9E). Finally, we integrated information of Age, Gender, TNM Stage, Stage, Grade, FPS, and nomogram to graph ROC curves for predicting 1-year, 3-year, and 5-year survival of HCC. We found that FPS and nomogram possessed the highest AUC values (Figures 9F–H), therefore, FPS and nomogram were trustworthy prognostic scoring tools for HCC.

FIGURE 9

FIGURE 9. Nomogram for Predicting HCC Survival. (A–B) Univariate and multivariate COX analysis of Age, Gender, Stage, Grade, and FPS in TCGA-LIHC. (C) Nomogram constructed from Stage and FPS information. (D) Calibration curves for 1 year, 3 years, 5 years survival. (E) Decision curve. (F–H) ROC curves for multiple factors predicting 1-year, 3-year, and 5-year survival in HCC.

Discussion

HCC remains a major cause of cancer death progressing rapidly (Llovet et al., 2021). Inhibitors of key genes in the Foxo signaling pathway were proposed to be an influential factor for novel targeted therapies for HCC (Sun et al., 2021). Therefore, identification of key regulators in the Foxo signaling pathway was essential to optimize HCC survival and treatment options.

Although a lot of progress in systemic therapy, for example, molecular targeting reagents, HCC is one of the worst prognostic tumors attributed to drug tolerance as well as frequent recurrence and metastasis (Oura et al., 2021). As research deepens, TME is receiving an increasing attention, and for this reason, immunosuppressive therapy has been launched (Llovet et al., 2022). A previous study showed that FOXO1 played a part in macrophages through transcriptionally controlling IRF-1/nitrio oxide (NO) axis and reduced the secretion of IL-6 from macrophages in TME indirectly (Cui et al., 2023). In order to Further explore the relationship between Foxo signaling pathway and TME, we also evaluated the immune landscapes in 3 FMSs molecular subtypes. Interestingly, the FMSs exhibited different immune landscapes. Tregs in the TME of FMS 3 patients exhibited increased infiltration abundance. Previous studies indicated that Tregs were significantly enriched and amplified in progressive HCC with the amplified Tregs leading to CD8⁺ T cells being depleted or their function being suppressed (Zheng et al., 2017). Another study revealed that in HCC, Tregs infiltration abundance was reduced in TME by inducing apoptosis to inhibit liver carcinogenesis (Zhou et al., 2021). Additionally, lower levels of immune checkpoint gene expression were discovered in FMS 1, and the highest level of immune checkpoint gene expression was demonstrated in FMS 3. Among these immune checkpoint genes, highly expressed CTLA4, LAG3 and TIGIT has been treated as diagnostic biomarkers for colorectal cancer (Sasidharan Nair et al., 2018). These findings altogether may account for the reason for excellent survival advantage in FMS 1 and provide promising therapy targets toward immune checkpoint gene expression discrepancy for HCC patients. The discovery of FMSs would contribute to refine HCC classification and further shed insights into the connections between Foxo signaling pathway, TME and targeted therapies for immune checkpoint gene.

FPS, a prognostic scoring system for HCC, was constrained based on the DEGs in FMSs. HCC patients with high TMB typically exhibited suboptimal prognostic status (Owada-Ozaki et al., 2018). The combined survival analysis of FPS and TMB showed that patients with low FPS and low TMB exhibited an excellent prognostic status. Lower TIDE score, IFNG score, Exclusion score, Dysfunction score, MDSC score were realized in low FPS, suggesting that low FPS are more sensitive to immunotherapy (Peng et al., 2021). In addition, Nomogram, constructed in accordance with FPS, was a potential tool for clinicians to accurately prognosticate patients based on their FPS characteristics.

9 key prognostic genes were screened out to construct an FPS model. Several of them were related to FOXO pathways. Numerous studies have shown that SIRT1, a member of the silent information regulator 2 (Sir2) family, plays an important role in the deacetylation of FOXO and regulation of autophagy in cells (Singh and Ubaid, 2020). Based on FOXO pathways genes, CENP-A was selected. A literature reported that SIRT7 could promote CENP-A assembly in nucleosome and restrain the tumorigenesis in gut (Liu et al., 2020). SIRT7 is also a member of the Sir2 family. Our findings may indicate that the Sir2 gene family has complex regulatory effects on FOXO pathways. Cyclin B2 (CCNB2) belongs to Foxo signaling pathway and is regarded as a promising biomarker for prognosis in LIHC (Li et al., 2021). Anti-silencing function 1b (ASF1b) was discovered as an oncogenic indicator for gastric cancer (Chen et al., 2022). Other study uncovered that ASF1B activated PI3K/AKT/mTOR signaling to promote cisplatin resistance in triple-negative breast cancer cells (Wang et al., 2022). Intestinally, GSEA analysis revealed that cell division cycle associated 8 (CDCA8) regulated gene sets relevant to PI3K/AKT/mTOR signaling (Bi et al., 2018). Collectively, we excavated some genes which were directly or indirectly related to Foxo signaling pathway or Foxo signaling pathway activity.

Notwithstanding original molecular subtype and prognostic tool based on Foxo signaling pathway signature was developed for HCC patients, some limitations need to be solved for future applying. First of all, more clinical statistics are necessary to calibrate the model. Secondly, 9 gens selected for FPS riskscore signature requires more studies on mechanism to confirm their roles as biomarkers for HCC prognosis.

Conclusion

Overall, we defined novel molecular subtypes in HCC in accordance with genes in the Foxo signaling pathway, providing new research perspectives in the study of HCC tumor heterogeneity. In this study, we also constructed the HCC prognostic scoring system for Foxo signaling pathway activity characteristics based on DEGs in FMSs, which demonstrated excellent robustness in assessing HCC prognosis, immune activity, and mutational status. Encouragingly, the FPS also exhibited excellent and promising outcome in guiding HCC drug selection. The clinical value of FPS was not further validated in a large sample of clinical cases. Further prospective trials are needed to investigate the prognostic value and therapeutic guidance of FPS in the future.

Data availability statement

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

Author contributions

All authors contributed to this present work: ST designed the study, YQ and ST acquired the data. YQ and ST drafted the manuscript, YQ and ST revised the manuscript. All authors contributed to the article and approved the submitted version.

Funding

This study was funded by National Key R&D Program of China (2021YFC2301800) and (2022YFC2304500).

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fphar.2023.1213506/full#supplementary-material

Abbreviations

HCC, Hepatocellular carcinoma; WHO, World Health Organization; EMT, epithelial mesenchymal transition; TGF-β, transforming growth factor-β; TCGA, The Cancer Genome Atlas; GEO; Gene expression omnibus; KEGG, Kyoto encyclopedia of genes and genomes; DEGs, differentially expressed genes; FMS, Foxo Molecular subtypes; AIC, Akaike information criterion; FPS, Foxo prognosis score; K-M, Kaplan-Meier; ROC, Receiver Operating Characteristic; CNV, copy number variation; SNV, Single nucleotide variant; TMB, Tumor mutation burden; TME, tumor microenvironment; ssGSEA, single sample gene set enrichment analysis; TIDE, Tumor Immune Dysfunction and Exclusion; IFNG, interferon gamma; MDSCs, myeloid-derived suppressor cells; GDSC, Genomics of Drug Sensitivity in Cancer; IC50, half maximal inhibitory concentration; DCA, Decision curve; Tregs, T cells regulatory.

References

Ahn, H., Kim, H., Abdul, R., Kim, Y., Sim, J., Choi, D., et al. (2018). Overexpression of forkhead box O3a and its association with aggressive phenotypes and poor prognosis in human hepatocellular carcinoma. Am. J. Clin. Pathol. 149 (2), 117–127. doi:10.1093/ajcp/aqx132

PubMed Abstract | CrossRef Full Text | Google Scholar

Auslander, N., Zhang, G., Lee, J. S., Frederick, D. T., Miao, B., Moll, T., et al. (2018). Robust prediction of response to immune checkpoint blockade therapy in metastatic melanoma. Nat. Med. 24 (10), 1545–1549. doi:10.1038/s41591-018-0157-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Barbie, D. A., Tamayo, P., Boehm, J. S., Kim, S. Y., Moody, S. E., Dunn, I. F., et al. (2009). Systematic RNA interference reveals that oncogenic KRAS-driven cancers require TBK1. Nature 462 (7269), 108–112. doi:10.1038/nature08460

PubMed Abstract | CrossRef Full Text | Google Scholar

Becht, E., Giraldo, N. A., Lacroix, L., Buttard, B., Elarouci, N., Petitprez, F., et al. (2016). Estimating the population abundance of tissue-infiltrating immune and stromal cell populations using gene expression. Genome Biol. 17 (1), 218. doi:10.1186/s13059-016-1070-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Bi, Y., Chen, S., Jiang, J., Yao, J., Wang, G., Zhou, Q., et al. (2018). CDCA8 expression and its clinical relevance in patients with bladder cancer. Med. Baltim. 97 (34), e11899. doi:10.1097/MD.0000000000011899

PubMed Abstract | CrossRef Full Text | Google Scholar

Blanche, P., Dartigues, J. F., and Jacqmin-Gadda, H. (2013). Estimating and comparing time-dependent areas under receiver operating characteristic curves for censored event times with competing risks. Stat. Med. 32 (30), 5381–5397. doi:10.1002/sim.5958

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, B., Khodadoust, M. S., Liu, C. L., Newman, A. M., and Alizadeh, A. A. (2018). Profiling tumor infiltrating immune cells with CIBERSORT. Methods Mol. Biol. 1711, 243–259. doi:10.1007/978-1-4939-7493-1_12

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, C., Bao, H., Lin, W., Chen, X., Huang, Y., Wang, H., et al. (2022). ASF1b is a novel prognostic predictor associated with cell cycle signaling pathway in gastric cancer. J. Cancer 13 (6), 1985–2000. doi:10.7150/jca.69544

PubMed Abstract | CrossRef Full Text | Google Scholar

Cui, X., Zhao, H., Wei, S., Du, Q., Dong, K., Yan, Y., et al. (2023). Hepatocellular carcinoma-derived FOXO1 inhibits tumor progression by suppressing IL-6 secretion from macrophages. Neoplasia 40, 100900. doi:10.1016/j.neo.2023.100900

PubMed Abstract | CrossRef Full Text | Google Scholar

Dong, T., Zhang, Y., Chen, Y., Liu, P., An, T., Zhang, J., et al. (2017). FOXO1 inhibits the invasion and metastasis of hepatocellular carcinoma by reversing ZEB2-induced epithelial-mesenchymal transition. Oncotarget 8 (1), 1703–1713. doi:10.18632/oncotarget.13786

PubMed Abstract | CrossRef Full Text | Google Scholar

European Association for the Study of the Liver (2018). EASL clinical practice guidelines: management of hepatocellular carcinoma. J. Hepatol. 69 (1), 182–236. doi:10.1016/j.jhep.2018.03.019

PubMed Abstract | CrossRef Full Text | Google Scholar

Geeleher, P., Cox, N., and Huang, R. S. (2014). pRRophetic: an r package for prediction of clinical chemotherapeutic response from tumor gene expression levels. PLoS One 9 (9), e107468. doi:10.1371/journal.pone.0107468

PubMed Abstract | CrossRef Full Text | Google Scholar

Grinchuk, O. V., Yenamandra, S. P., Iyer, R., Singh, M., Lee, H. K., Lim, K. H., et al. (2018). Tumor-adjacent tissue co-expression profile analysis reveals pro-oncogenic ribosomal gene signature for prognosis of resectable hepatocellular carcinoma. Mol. Oncol. 12 (1), 89–113. doi:10.1002/1878-0261.12153

PubMed Abstract | CrossRef Full Text | Google Scholar

Hanzelmann, S., Castelo, R., and Guinney, J. (2013). Gsva: gene set variation analysis for microarray and RNA-seq data. BMC Bioinforma. 14, 7. doi:10.1186/1471-2105-14-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Harrell, F. E., and Harrell, M. F. E. (2019). Package ‘hmisc’. CRAN2018 2019, 235–236.

Jiang, P., Gu, S., Pan, D., Fu, J., Sahu, A., Hu, X., et al. (2018). Signatures of T cell dysfunction and exclusion predict cancer immunotherapy response. Nat. Med. 24 (10), 1550–1558. doi:10.1038/s41591-018-0136-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, T., Fu, J., Zeng, Z., Cohen, D., Li, J., Chen, Q., et al. (2020). TIMER2.0 for analysis of tumor-infiltrating immune cells. Nucleic Acids Res. 48 (W1), W509–W14. doi:10.1093/nar/gkaa407

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, Y., Li, P., and Wang, N. (2021). Effect of let-7c on the PI3K/Akt/FoxO signaling pathway in hepatocellular carcinoma. Oncol. Lett. 21 (2), 96. doi:10.3892/ol.2020.12357

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, Z., Kwon, S. M., Li, D., Li, L., Peng, X., Zhang, J., et al. (2022). Human constitutive androstane receptor represses liver cancer development and hepatoma cell proliferation by inhibiting erythropoietin signaling. J. Biol. Chem. 298 (5), 101885. doi:10.1016/j.jbc.2022.101885

PubMed Abstract | CrossRef Full Text | Google Scholar

Lian, Q., Wang, S., Zhang, G., Wang, D., Luo, G., Tang, J., et al. (2018). Hccdb: a database of hepatocellular carcinoma expression Atlas. Genomics Proteomics Bioinforma. 16 (4), 269–275. doi:10.1016/j.gpb.2018.07.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, X., Li, C., Li, Q., Chang, H. C., and Tang, Y. C. (2020). SIRT7 facilitates CENP-A nucleosome assembly and suppresses intestinal tumorigenesis. iScience 23 (9), 101461. doi:10.1016/j.isci.2020.101461

PubMed Abstract | CrossRef Full Text | Google Scholar

Llovet, J. M., Castet, F., Heikenwalder, M., Maini, M. K., Mazzaferro, V., Pinato, D. J., et al. (2022). Immunotherapies for hepatocellular carcinoma. Nat. Rev. Clin. Oncol. 19 (3), 151–172. doi:10.1038/s41571-021-00573-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Llovet, J. M., Kelley, R. K., Villanueva, A., Singal, A. G., Pikarsky, E., Roayaie, S., et al. (2021). Hepatocellular carcinoma. Nat. Rev. Dis. Prim. 7 (1), 6. doi:10.1038/s41572-020-00240-3

CrossRef Full Text | Google Scholar

Mayakonda, A., Lin, D. C., Assenov, Y., Plass, C., and Koeffler, H. P. (2018). Maftools: efficient and comprehensive analysis of somatic variants in cancer. Genome Res. 28 (11), 1747–1756. doi:10.1101/gr.239244.118

PubMed Abstract | CrossRef Full Text | Google Scholar

Oura, K., Morishita, A., Tani, J., and Masaki, T. (2021). Tumor immune microenvironment and immunosuppressive therapy in hepatocellular carcinoma: a review. Int. J. Mol. Sci. 22 (11), 5801. doi:10.3390/ijms22115801

PubMed Abstract | CrossRef Full Text | Google Scholar

Owada-Ozaki, Y., Muto, S., Takagi, H., Inoue, T., Watanabe, Y., Fukuhara, M., et al. (2018). Prognostic impact of tumor mutation burden in patients with completely resected non-small cell lung cancer: brief report. J. Thorac. Oncol. 13 (8), 1217–1221. doi:10.1016/j.jtho.2018.04.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Peng, Y., Liu, C., Li, M., Li, W., Zhang, M., Jiang, X., et al. (2021). Identification of a prognostic and therapeutic immune signature associated with hepatocellular carcinoma. Cancer Cell Int. 21 (1), 98. doi:10.1186/s12935-021-01792-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Ritchie, M. E., Phipson, B., Wu, D., Hu, Y., Law, C. W., Shi, W., et al. (2015). Limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res. 43 (7), e47. doi:10.1093/nar/gkv007

PubMed Abstract | CrossRef Full Text | Google Scholar

Sasidharan Nair, V., Toor, S. M., Taha, R. Z., Shaath, H., and Elkord, E. (2018). DNA methylation and repressive histones in the promoters of PD-1, CTLA-4, TIM-3, LAG-3, TIGIT, PD-L1, and galectin-9 genes in human colorectal cancer. Clin. Epigenetics 10 (1), 104. doi:10.1186/s13148-018-0539-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Shannon, P., Markiel, A., Ozier, O., Baliga, N. S., Wang, J. T., Ramage, D., et al. (2003). Cytoscape: A software environment for integrated models of biomolecular interaction networks. Genome Res. 13 (11), 2498–2504. doi:10.1101/gr.1239303

PubMed Abstract | CrossRef Full Text | Google Scholar

Shen, W., Song, Z., Zhong, X., Huang, M., Shen, D., Gao, P., et al. (2022). Sangerbox: a comprehensive, interaction-friendly clinical bioinformatics analysis platform. iMeta 1 (3), e36. doi:10.1002/imt2.36

CrossRef Full Text | Google Scholar

Singh, V., and Ubaid, S. (2020). Role of silent information regulator 1 (SIRT1) in regulating oxidative stress and inflammation. Inflammation 43 (5), 1589–1598. doi:10.1007/s10753-020-01242-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Song, S. S., Ying, J. F., Zhang, Y. N., Pan, H. Y., He, X. L., Hu, Z. M., et al. (2020). High expression of FOXO3 is associated with poor prognosis in patients with hepatocellular carcinoma. Oncol. Lett. 19 (4), 3181–3188. doi:10.3892/ol.2020.11430

PubMed Abstract | CrossRef Full Text | Google Scholar

Sun, E. J., Wankell, M., Palamuthusingam, P., McFarlane, C., and Hebbard, L. (2021). Targeting the PI3K/Akt/mTOR pathway in hepatocellular carcinoma. Biomedicines 9 (11), 1639. doi:10.3390/biomedicines9111639

PubMed Abstract | CrossRef Full Text | Google Scholar

Szklarczyk, D., Gable, A. L., Nastou, K. C., Lyon, D., Kirsch, R., Pyysalo, S., et al. (2021). The STRING database in 2021: customizable protein-protein networks, and functional characterization of user-uploaded gene/measurement sets. Nucleic Acids Res. 49 (D1), D605–D612. doi:10.1093/nar/gkaa1074

PubMed Abstract | CrossRef Full Text | Google Scholar

Therneau, T. M., and Lumley, T. (2015). Package ‘survival. R. Top. Doc. 128 (10), 28–33.

Google Scholar

Thorsson, V., Gibbs, D. L., Brown, S. D., Wolf, D., Bortone, D. S., Ou Yang, T. H., et al. (2018). The immune landscape of cancer. Immunity 48 (4), 812–830.e14. doi:10.1016/j.immuni.2018.03.023

PubMed Abstract | CrossRef Full Text | Google Scholar

Vogel, A., Meyer, T., Sapisochin, G., Salem, R., and Saborowski, A. (2022). Hepatocellular carcinoma. Lancet 400 (10360), 1345–1362. doi:10.1016/S0140-6736(22)01200-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, L., Yang, X., Zhou, F., Sun, X., and Li, S. (2022). Circular RNA UBAP2 facilitates the cisplatin resistance of triple-negative breast cancer via microRNA-300/anti-silencing function 1B histone chaperone/PI3K/AKT/mTOR axis. Bioengineered 13 (3), 7197–7208. doi:10.1080/21655979.2022.2036894

PubMed Abstract | CrossRef Full Text | Google Scholar

Wilkerson, M. D., and Hayes, D. N. (2010). ConsensusClusterPlus: a class discovery tool with confidence assessments and item tracking. Bioinformatics 26 (12), 1572–1573. doi:10.1093/bioinformatics/btq170

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, S., Pang, L., Dai, W., Wu, S., Ren, T., Duan, Y., et al. (2021). Role of forkhead box O proteins in hepatocellular carcinoma biology and progression (review). Front. Oncol. 11, 667730. doi:10.3389/fonc.2021.667730

PubMed Abstract | CrossRef Full Text | Google Scholar

Zheng, C., Zheng, L., Yoo, J. K., Guo, H., Zhang, Y., Guo, X., et al. (2017). Landscape of infiltrating T cells in liver cancer revealed by single-cell sequencing. Cell 169 (7), 1342–1356. doi:10.1016/j.cell.2017.05.035

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhou, Y. F., Song, S. S., Tian, M. X., Tang, Z., Wang, H., Fang, Y., et al. (2021). Cystathionine beta-synthase mediated PRRX2/IL-6/STAT3 inactivation suppresses Tregs infiltration and induces apoptosis to inhibit HCC carcinogenesis. J. Immunother. Cancer 9 (8), e003031. doi:10.1136/jitc-2021-003031

PubMed Abstract | CrossRef Full Text | Google Scholar