Gene Expression Profiling of Advanced Stage Hepatocellular Carcinoma: A Bioinformatic Analysis.

Introduction
Hepatocellular carcinoma (HCC) ranks as the fourth leading cause of mortality worldwide. Globally, HCC is now the sixth most diagnosed cancer, and its incidence continues to rise, with projections estimating 1.4 million new cases by 2040 [1]. Key risk factors for HCC include chronic infection with hepatitis B virus or hepatitis C virus, a history of heavy alcohol consumption, diabetes, and nonalcoholic fatty liver disease [2].
The Barcelona Clinic Liver Cancer (BCLC) system is a widely used staging system in modern medical practice. The European Association for the Study of the Liver recommends the BCLC staging system for determining prognosis and guiding appropriate therapy. The BCLC system categorizes HCC into five stages, linking each stage to specific treatment options. It uses the Child-Pugh score, tumor size, and the Eastern Cooperative Oncology Group scale to assess performance status [2]. Although clinical and histopathological staging remains in widespread use, these techniques have failed to accurately predict treatment response or prognosis [3]. Therefore, more personalized strategies are indicated, such as precision medicine that considers tumor-specific molecular targets in the planning of therapy. Next-generation sequencing illustrates the molecular heterogeneity of HCC and possibly can be utilized to identify susceptibility to the tumor and monitor postoperative cancers [4]. According to a review conducted by the International Liver Cancer Association, a few biomarkers have been implemented into HCC clinical practice despite the need for enhanced patient stratification at several levels of clinical management [3]. Using molecular and genetic information, like combining genomic, transcriptomic, and immunologic data, has been proven to improve predictions about outcomes and how well HCC patients will respond to immunotherapy. Furthermore, these approaches facilitate the recognition of tumor molecular subtypes, thus allowing a more tailored and effective therapeutic approach [5].
Patients with the same HCC stage may show different therapeutic responses and clinical outcomes due to genetic, immunological and tumor microenvironmental variations. This emphasizes the need for molecular-based patient stratification. A study highlighted that the integration of molecular biomarkers in clinical systems can improve the accuracy of prognosis and therapy selection [6]. Bioinformatics studies represent a recent approach to biomarker discovery. These studies are significant because they enable researchers to conduct large-scale gene activity studies, discover biomarkers, and map how tumors interact with the immune system. Utilizing publicly accessible datasets, such as the Gene Expression Omnibus (GEO), these studies can identify molecular subtypes that traditional clinical methods are unable to discern. This, in turn, facilitates the development of more precise and targeted treatments [7].
Furthermore, immune-tumor system interactions play a pivotal role in patient course and treatment response in HCC. It is important to research the same to understand the body dynamics properly in patients, as seen in a model that examined the influence of lymphocyte counts on radiotherapy response. The results indicated that the shorter radiotherapy fractionation regimens may reduce the time of lymphocyte depletion and accelerate recovery after treatment. Understanding these interactions better could help find important biomarkers and create more tailored treatments, as well as enhance treatment methods based on molecular biology for HCC [8].
This study aims to analyze and elucidate genetic expression, crucial genetic markers, and immune-tumor interactions in HCC samples classified as BCLC-B and BCLC-C using datasets available in GEO.

Materials and Methods

Study Design and Setting
This study aims to analyze and elucidate genetic expression, crucial genetic markers, and immune-tumor interactions in HCC samples classified as BCLC-B and BCLC-C using datasets available in GEO. This study is a bioinformatic study that involved data collection and analysis using several database and software.

Data Collection
We obtained the gene expression datasets from the GEO database. We only considered datasets that were classified as BCLC B and BCLC C and had at least 100 significant genes in the volcano plot analysis, which compared the combined BCLC B and C data against non-BCLC B and C datasets. Data collection was carried out from October 21, 2024, to November 4, 2024. Subsequent bioinformatic analyses commenced until November 18, 2024.

Differentially Expressed Genes Screening
We employed the GEO2R tool (https://www.ncbi.nlm.nih.gov/geo/geo2r/) to compare multiple sample groups and detect genes that exhibited differential expression across the experimental conditions. We identified differentially expressed genes (DEGs) between BCLC stages B and C with other stage samples from the dataset, using a threshold of |log2 Fold change (FC)| > 1. We applied the Benjamini-Hochberg method for adjustment. We displayed the significant DEGs using a volcano plot. We identified common genes using a Venn diagram.

Protein-Protein Interaction Network Construction
We explored the interactions between proteins and pathways by using the STRING database (https://string-db.org/) to build a protein-protein network (PPI) for the significant DEGs. We then applied Cytoscape software to identify the key interacting genes.

Hub-Bottleneck Genes Identification
We analyzed the network using the Network Analyzer and CytoHubba plugins. Four key centrality measures—degree, betweenness, closeness, and stress—were employed to rank the nodes. Hub genes were defined as those ranking in the top 20 based on degree, while bottleneck nodes were identified as the top 20 based on betweenness. Likewise, the top 20 nodes were also highlighted according to their closeness and stress centrality scores. Crucial genes, referred to as hub-bottleneck genes, were determined based on the overlap between the hub and bottleneck categories. The regulatory behavior of these critical genes was subsequently validated using the log2 FC score.

Function Enrichment
To study the important genes’ biological roles, we used the ClueGO plugin within Cytoscape to perform enrichment analyses for Gene Ontology (GO), Kyoto Encyclopedia of Genes (KEGG), and Reactome pathways. A term should consist of a minimum of two genes. We determined a p-value of less than 0.05 as the statistical threshold.

Immune Infiltration
To understand tumor-immune interaction, we used TIMER2 to find the influence of the specific immune cells in liver cancer on gene expression levels. We created a scatter plot to visualize the result.

Results

Identification of Differentially Expressed Genes
We took two datasets according to the specified inclusion criteria, namely the GSE222334 and GSE56545. Dataset GSE222334 consists of 100 HCC patients, and dataset GSE56545 consists of 35 HCC patients. A total of 3314 DEGs were found between the two datasets. There were 753 DEGs in the dataset GSE56545, of which 372 genes were upregulated, and 382 genes were downregulated. There were 2560 DEGs in the dataset GSE22234, with 1728 upregulated and 832 downregulated genes. The volcano plot in Figure 1 depicted this differential expression pattern. The red dot represented upregulated genes, while the blue dot represented downregulated genes. Then, we used the Venn diagram to identify the common DEGs between the two datasets. There were 75 overlapped DEGs, which would be further analyzed. Figure 2 showed the Venn diagram.

Protein-Protein Interaction Network
We constructed a PPI network using DEGs identified between BCLC stage B–C and non-B–C stages. Among 75 common DEGs, 63 genes were recognized and included in the network. The resulting PPI network consisted of 63 nodes and 99 edges, as illustrated in Figure 3. The network showed a significant enrichment p-value of 0.0143, indicating that the observed interactions occurred more frequently than expected by chance. Genes without any interactions were subsequently excluded from further analysis.

Identification of Candidate Crucial Genes
TUBB2A had the highest degree score of 26, identifying it as the most connected gene within the network. Its closest competitors, thioredoxin (TXN) and HSPA6, each had 24 connections. TUBB2A also ranked highest in betweenness centrality, with a score of 0.407, highlighting its key role in bridging various components of the network. Additionally, TUBB2A demonstrated the highest closeness score of 0.527 and stress score of 1760, further confirming its central position in regulating network dynamics and managing information flow. A total of 13 crucial genes, classified as hub-bottleneck genes, are listed in Table 1. Figure 4 visualizes the PPI construction for these genes, which includes one main component and two genes located outside the primary network cluster.

Functional Enrichment Analysis
The analysis of GO showed that NQO1, SRXN1, TXN, TXNRD1, CRYAA, and HSPA6 are important for “antioxidant activity,” “controlling transport inside cells,” and “helping proteins fold correctly.” KEGG analysis showed that HDAC1, HSPA6, NQO1, and TXNRD1 are involved in pathways such as “longevity regulation pathways,” “ubiquinone and other terpenoid-quinone biosynthesis,” and “selenocompound metabolism.” Reactome analysis showed that HDAC1, TXN, TXNRD1, SRXN1, and NQO1 are part of processes like “FOXO-mediated transcription,” “the reaction of oxidized thioredoxin with NADPH and H to form reduced thioredoxin and NADP,” and “nuclear activities controlled by NFE2L2,” as illustrated in Figure 5. Two genes were involved in those three analyses, which were NQO1 and TXNRD1.

Tumor Infiltration Analysis

TXNRD1 and NQO1 showed the strongest positive correlation with neutrophils with partial core values of 0.322 and 0.201, respectively (p-value <0.001). The findings suggested that higher neutrophil infiltration is associated with increased crucial gene expression levels. As shown in Figure 6, The y-axis showed gene expression levels, and the x-axis indicated immune infiltration levels.

Discussion
Gene mutations play a critical role in initiating carcinogenesis. Among various malignancies, HCC exhibits one of the highest numbers of mutations per genome [9]. Molecular studies have identified several key genetic and epigenetic alterations in HCC, including mutations in the TERT promoter, TP53, and CTNNB1, as well as widespread epigenetic abnormalities. Epigenetics refers to heritable changes in gene expression that occur without alterations to the deoxyribonucleic acid (DNA) sequence [10].
Numerous transcriptomic studies have investigated gene expression in HCC, including large-scale initiatives such as The Cancer Genome Atlas. For instance, Sarathi and Palaniappan (2019) found new genes that are expressed differently at various stages of HCC, showing unique molecular patterns linked to tumor growth [11]. A previous study also demonstrated the oncogenic role of SRXN1 in promoting HCC cell proliferation, migration, and invasion through its involvement in the antioxidant response [12]. Other studies have found possible indicators like FKBP11, SCRIB, SLC38A2, SORBS2, and STAB2 that are connected to the early return of HCC, but these were not directly related to the stage of the disease [13].
Bioinformatic analysis identified two key genes TXNRD1 and NQO as being strongly associated with BCLC stage B and C HCC. The liver plays a vital role in detoxifying substances that generate reactive oxygen species (ROS). Disruption of the hepatic reduction–oxidation (redox) system leads to oxidative stress, which has been closely linked to liver cancer development [14]. Thioredoxin is essential in maintaining cellular redox homeostasis and includes several components: cytoplasmic TXN, mitochondrial TXN2, cytoplasmic TXNRD1, mitochondrial TXNRD2, and the TXN inhibitor TXN-interacting protein [15]. Research by Cho et al. found that mutations in TXN and TXNRD1 are significantly associated with worse overall survival and disease-free survival, with TXNRD1 expression being 1.8-fold higher in HCC than in controls [16].
TXNRD1 is a key redox regulator in mammalian cells, functioning to maintain thioredoxin 1 (Trx1) in its reduced, active state [17]. Trx1 serves as an antioxidant that protects cells from oxidative stress and plays a crucial role in regulating cell proliferation and survival [18]. Trx1 also works with the tumor suppressor PTEN, and when TXNRD1 disrupts this connection, it speeds up the breakdown of PTEN, which activates the Akt/mammalian target of rapamycin (mTOR) signaling pathway and helps liver cancer cells spread [17]. TXNRD1 is an important gene that plays a key role in changing how cells use energy and is believed to help control oxidative stress and fat damage in liver cancer [19].
In addition, TXN stabilizes the expression of BTB and CNC Homology 1 (BACH1) by inhibiting its ubiquitination and directly interacting with it. Once stabilized, BACH1 translocates into the nucleus and activates the AKT/mTOR signaling pathway, enhancing stemness and metastatic potential in HCC. In metastatic HCC, dual inhibition of TXN and treatment with lenvatinib has demonstrated a synergistic therapeutic effect, suggesting that targeting the TXN–BACH1 axis could be a promising strategy for treating metastatic HCC [20].
We also identified NAD(P)H Quinone Dehydrogenase 1 (NQO1) as a gene associated with HCC in BCLC stages B and C. NQO1 plays a critical role in cellular defense by detoxifying quinones and preventing oxidative damage [21]. However, its overexpression or mutation can paradoxically increase the production of ROS, contributing to oxidative stress. This elevated stress activates signaling pathways that promote DNA damage, chronic inflammation, and cell proliferation—processes that are central to the initiation of hepatocarcinogenesis [22].
Recent evidence also suggests that NQO1 overexpression destabilizes sirtuin 6 (SIRT6), making it more vulnerable to degradation by the 26S proteasome. SIRT6, a nicotinamide adenine dinucleotide-dependent deacetylase, exhibits context-dependent effects during HCC progression. In the early stages of tumorigenesis, SIRT6 acts as a tumor suppressor by maintaining genomic stability, inhibiting inflammation, and modulating DNA repair mechanisms, thereby preventing malignant transformation. However, in established tumors, particularly in the context of HCC, reduced levels of SIRT6 contribute to increased activation of the AKT signaling pathway. Specifically, decreased SIRT6-mediated deacetylation of AKT leads to enhanced AKT phosphorylation and activity, which in turn stabilizes the anti-apoptotic protein XIAP. This stabilization allows liver cancer cells to resist apoptosis and promotes tumor cell survival [22]. Thus, while SIRT6 may suppress tumor initiation, its downregulation in established tumors facilitates progression and resistance to cell death [23].
A profound understanding of the signaling pathways associated with HCC is also essential for cancer diagnosis and the discovery of potential therapeutic targets. There are vital pathways that are repeatedly disrupted in HCC, such as tyrosine kinase-dependent signaling pathways, including Ras/Raf/MEK/extracellular-signal-regulated kinase and phosphoinositide three kinase/AKT/mTOR, Wnt/β-catenin, and Janus-kinase (JAK)/signal transducer activator of transcription (STAT) [24-25]. Lenvatinib-resistant cells have shown overexpression of the NQO1 gene. NQO1 overexpression reduces lenvatinib-induced ROS production, consequently suppressing apoptosis. Combining lenvatinib with dicoumarol, an NQO1 inhibitor, can overcome this resistance [26]. TXNRD1 also counteracts intracellular ROS production in patients with HCC. Genetic inhibition of TXNRD1 has been shown to suppress cell proliferation and induce apoptosis in vitro. Accordingly, administration of auranofin, a TXNRD1 inhibitor, may suppress HCC progression [27]. The PI3K/AKT/mTOR pathway helps cells grow and spread tumors, and it is often turned on in HCC, making it a key reason why treatments may not work. This mTOR pathway is activated by the TXNRD1 and TXN genes, which help HCC cells spread. Bioinformatic studies show that TXNRD1 and TXN are found in HCC cases at BCLC stages B and C. The Wnt/β-catenin pathway is also wrongly activated in advanced HCC and affects how tumors move, invade, maintain stem cells, and other functions [25].
Infiltration analysis revealed correlation with immune cells, mainly neutrophils. Nowadays, the key point of HCC immunotherapy is intruding on the immune checkpoint. Neutrophils have emerged as a novel key pathogenesis of HCC by affecting immunosuppression, cell survival, and progression. A study discovered a novel therapeutic strategy that targets tumor-infiltrating neutrophils, which might alter the tumor microenvironment to be more beneficial to systemic immunotherapies. This implies the possibility of a more potent therapeutic approach in HCC [28]. In liver injury, TXNRD1 influences the infiltration of inflammatory cells, including neutrophils, macrophages, T-lymphocytes, and B-lymphocytes. Neutrophil infiltration is suppressed in hepatocytes expressing TXNRD1. Neutrophils are attracted to hepatocytes by certain signals called chemokines, like Cxcl1, Cxcl2, Cxcr2, and the increase of ICAM1. Research has found that in hepatocytes that express TXNRD1, the levels of Cxcl1, Cxcl2, Cxcr2, Egr1, and ICAM1 are much lower, which matches the decrease in neutrophil presence [29].
Despite providing valuable insights, this study has several limitations that should be acknowledged to offer a balanced perspective and guide future research. First, the analysis was based on two publicly available datasets (GSE222334 and GSE56545) with a total of 135 HCC patients. While sufficient for exploratory bioinformatics, this sample size may not fully represent the genetic diversity of advanced-stage HCC, and the datasets lacked comprehensive clinical annotations such as treatment response, comorbidities, and survival data. Second, all findings were derived from in silico analysis without experimental validation. Although bioinformatics tools are powerful for hypothesis generation, the expression, and biological roles of the identified genes particularly TXNRD1 and NQO1 require confirmation through laboratory techniques such as qRT-PCR or immunohistochemistry. Third, although both genes have been linked to chemoresistance in other cancers, their direct association with therapeutic resistance in HCC was not assessed. Exploring this relationship in the context of standard treatments such as sorafenib or immune checkpoint inhibitors would be an important step forward. Fourth, the analysis of tumor immune infiltration was limited to neutrophils, without a broader evaluation of other immune components like T-cells, macrophages, or dendritic cells, which play crucial roles in the tumor microenvironment and immunotherapy response. Lastly, the study’s dependence on specific bioinformatics platforms (GEO2R, STRING, Cytoscape, and CytoHubba) introduces potential algorithmic and selection biases. These tools, while widely accepted, operate on inherent assumptions that may influence the interpretation of network topology and enrichment results. Future studies should aim to validate these bioinformatic findings experimentally, expand the scope of immune profiling, integrate multi-omics data, and explore the clinical relevance of TXNRD1 and NQO1 in chemoresistance and prognosis for patients with advanced-stage HCC.
In conclusion, this study identifies several genes related to the advanced stage of HCC, of which the two most crucial genes are NQO1 and TXNRD1. These genes are also involved in biological processes, like longevity regulation pathways, antioxidant activity, and FOXO-mediated transcription. Both genes also have a strong correlation with neutrophil infiltration. These findings suggest its potential as a biomarker and therapeutic target candidate for advanced-stage HCC.

Author Contribution Statement
All authors contributed to the study conception and design. Study conception and design were made by Ajib Zaim Alamsyah and Putu Enrico Pramana Okaniawan. The acquisition of data and critical revision were performed by Putu Enrico Pramana Okaniawan and Kadek Mercu Narapati Pamungkas. Data analysis and interpretation were performed by Ni Nyoman Gita Kharisma Dewi and Ajib Zaim Alamsyah. Drafting of manuscript were written by Kadek Mercu Narapati Pamungkas and Ni Nyoman Gita Kharisma Dewi. Revision of manuscript, final approval of manuscript, and study supervision were performed by Dwijo Anargha Sindhughosa and I Ketut Mariadi.