Screening of co-pathogenic genes of non-alcoholic fatty liver disease and hepatocellular carcinoma

Background Non-alcoholic fatty liver disease (NAFLD) is a risk factor for hepatocellular carcinoma (HCC). However, its carcinogenic mechanism is still unclear, looking for both diseases’ transcriptome levels, the same changes as we are looking for NAFLD may provide a potential mechanism of action of HCC. Thus, our study aimed to discover the coexisting pathogenic genes of NAFLD and HCC. Methods We performed a variance analysis with public data for both diseases. At the same time, weighted gene correlation network analysis (WGCNA) was used to find highly correlated gene modules in both diseases. The darkturquoise gene module was found to be highly correlated with both diseases. Based on the diagnosis related module genes and the differential genes of the two diseases, we constructed diagnostic and prognostic models by logistic regression, univariate Cox regression, and LASSO regression. Public datasets verified the results. Meanwhile, we built a competing endogenous RNA (ceRNA) network based on the model genes and explored the related pathways and immune correlation involved in the two diseases by using Gene Ontology, Kyoto Encyclopedia of Genes and Genomes, and gene set enrichment analyses. Immunohistochemistry was used to verify the different expression of ABCC5 and TUBG1 among the normal liver, NAFLD, and HCC tissues. Sodium palmitate/sodium oleate was used to establish high-fat cell models, and Real Time Quantitative Polymerase Chain Reaction (RT-qPCR) was used to verify the messenger RNA (mRNA) expression of ABCC5 in lipidization cells. Results A total of 26 upregulated genes and 87 downregulated genes were found using limma package identification analysis. According to WGCNA, the darkturquoise gene module was highly correlated with the prognosis of both diseases. The coexisting genes acquired by the two groups were only three central genes, that is, ABCC5, DHODH and TUBG1. The results indicated that the diagnostic and prognostic models constructed by ABCC5 and TUBG1 genes had high accuracy in both diseases. The results of immunohistochemistry showed that ABCC5 and TUBG1 were significantly overexpressed in NAFLD and HCC tissues compared with normal liver tissues. The Oil Red O staining and triglyceride identified the successful construction of HepG2 and LO2 high-fat models using PA/OA. The results of RT-qPCR showed that the lipidization of LO2 and HepG2 increased the mRNA expression of ABCC5. Conclusions The gene model constructed by ABCC5 and TUBG1 has high sensibility and veracity in the diagnosis of NAFLD as well as the diagnosis and prognosis of HCC. ABCC5 and TUBG1 may play an important role in the development of NAFLD to HCC. In addition, lipidization could upregulate the mRNA expression of ABCC5 in HCC.


Introduction
Cancers like hepatocellular carcinoma (HCC) have one of the highest mortalities in the world. Pathological metabolic diseases are the primary causes of HCC, such as hepatitis C or B and alcoholic and non-alcoholic fatty liver disease (NAFLD) (1). Due to the rapid increase in obesity and NAFLD patients, the prevalence of NAFLD-related HCC has increased accordingly (2). As the harmful effects of NAFLD can ultimately cause cirrhosis and even cancer, the prevalence of HCC has also started to rise suddenly. Based on some research findings (3), NAFLD can act synergistically with other HCC risk factors to accelerate cancer development. Over the past decade, an analysis of US health insurance databases found (4), according to the report, that 59% of HCC patients enrolled in medical plans are believed to have NAFLD as their primary cause. In a meta-analysis by Younossi et al. (5), all studies from 1989 to 2015 involving NAFLD were screened and estimated to determine point estimates [95% confidence intervals (CIs)] of HCC incidence, mortality, and incidence rates using a random-effects model. In NAFLD patients, HCC incidence was found to be 0.44/1,000 person-years (range: 0.29-0.66). A retrospective cohort study (6) evaluated the cancer risk in patients with NAFLD; it is also required to follow up the successive prevalence of HCC. Tumor biomarkers can assist medical researchers with early screening, diagnosis, treatment evaluation, recurrence, and prognosis prediction of tumors. With the prevalence of HCC and the estimated continued increase in this tumor (7), it is suggested that early monitoring and screening of tumor markers can reduce the mortality rate. In fact, in most patients, even after a complete resection or ablation of HCC tumors, the remaining carcinogenic tissue will progress to incurable terminal disease (8), thus, early detection and prevention of HCC development are considered to be the most effective strategies worldwide for improving patient outcomes. Up to now, the exact pathogenesis of NAFLD transition to HCC has not been explained. Identifying new marker genes can be used to predict risk and target treatment for individual patients. Hence, the prediction of candidate genes based on NAFLD-HCC pathogenesis is also an option. Based on bioinformatics analysis, earlier studies have established that telomerase reverse transcriptase (hTERT) (9), hypoxia-inducible transcription factor-2a (HIF-2a) (10), S100 calcium-binding protein A11(S100A11) (11), squalene epoxidase (SQLE) (12), and deoxycytidine kinase (DCK) (13) were associated with NAFLD to HCC progression and prognosis and could act as a related marker. Therefore, in an era of precision medicine, the diagnostic markers of NAFLD progressing to HCC need to be explored further.
In the human genome, the ATP-binding cassette (ABC) transporter is one of the largest protein families. More than 48 genes encoding human ABC transporters have been identified and sequenced, and the largest member of the family is multidrug resistance-associated protein (MRP), including MRP1 (ABCC1), MRP5 (ABCC5), and MRP6 (ABCC6) (14). In addition to its role in chemotherapy-resistant tumors, as one of the major classes of membrane ATPase, it is also implicated in the treatment failure of HCC with a multidrug-resistant phenotype (15). Studies in NAFLD have shown that several ABC transporter family members, such as ABCA1 and ABCC5, are closely related to the progression of NAFLD (16).
g-Tubulin (TUBG1) is a member of the human tubulin family. There are five known tubulin subtypes in the human tubulin family based on the c-terminal region with the greatest variation: a-tubulin, b-tubulin, g-tubulin, d-tubulin, and ϵ-tubulin. Complex cortical malformations related to tubulin gene mutations, also known as tubulin diseases, are a heterogeneous group of diseases with a wide range of clinical severity (14). A family of GTP enzymes called tubulin is highly enriched in microtubules and centrosomes.
This study aimed to identify candidate biomarkers for predicting NAFLD and HCC diagnoses by using the full public repository for Gene Expression Omnibus (GEO) and The Cancer Genome Atlas (TCGA) databases. We screened differentially expressed genes (DEGs) from public datasets of NAFLD and HCC and also screened overlapping genes between DEGs and the constructed weighted gene correlation network analysis (WGCNA); constructed NAFLD and HCC models, functional and pathway enrichment analyses, established ceRNA network, and protein-protein Interaction (PPI) network, respectively, and finally identified the hub gene. Therefore, this study hopes to provide new candidate diagnostic markers during the progression of NAFLD to HCC, offer the best screening strategy for high-risk patients, and advance the accuracy of HCC prevention.

Data download and processing
Through the GEOquery package (17), NAFLD datasets GSE48452, GSE89632, and GSE37031 were downloaded from the GEO database. The GSE48452 dataset is from Homo sapiens, and the data platform is GPL11532, which has 73 samples in total, including 14 normal liver tissue control samples, 27 liver tissue samples from obese patients, and 32 NAFLD tissue samples (18).This study included 14 normal liver tissue control samples and 32 NAFLD tissue samples from the dataset. The GSE89632 dataset is from Homo sapiens, and the data platform is GPL14951, which contains 63 samples, including 24 normal liver tissue control samples and 39 NAFLD tissue samples. This study included all the samples (19). The GSE37031 dataset is from Homo sapiens, and the data platform is GPL14877, containing 15 samples, including 7 normal liver tissue control samples and 8 NAFLD tissue samples. In this study, all samples were included as a validation set (20). GSE48452 and GSE89632 datasets are combined, and batches are removed through the SVA package (21). The data were normalized and standardized through the limma package (22).
Through the TCGAbiolinks package (23), the hepatoma dataset was downloaded from the TCGA database, TCGA-LIHC (Liver Hepatocellular Carcinoma, n=445 cases). The data type Count was selected and converted to the TPM format. HCC data were taken from the International Cancer Genome Consortium database, and TCGA-LIHC (n = 419 cases) was used as a validation set.

Model building and validation
For the construction of diagnostic models, the minimization of absolute contraction and selection operator LASSO regression are commonly used machine learning algorithms. In curve fitting, regularization is used to solve the overfitting and improve model accuracy. The glmnet package (24) was used to construct the model, and the parameter was set as: seed (2), family = " binomial ".
The prognostic model of HCC patients was developed using univariate and multivariate COX regression.

Variance analysis
The limma package was used for the gene difference analysis of GEO chip data in the normal group and disease group, and R package DESeq2 was used for the gene difference analysis of the normal group and disease group (25). LogFC >0.5 and the adjp value <0.01 were fixed as the threshold of differential genes. LogFC >0.5 was considered to be upregulated in the high-risk group, while logFC <0.5 was downregulated. The outcomes of difference analysis were represented by R package pheatmap heatmap and GGplot2 volcano map (26).
Analysis of enrichment (gene ontology/ kyoto encyclopedia of genes and genomes/gene set enrichment analysis/ GSVA/ssGSEA) Gene Ontology (GO) is a common analysis for large-scale functional enrichment studies, such as biological processes (BPs), molecular functions (MFs), and cellular component (CC). The Kyoto Encyclopedia of Genes and Genomes (KEGG) is a widely used database keeping information about genomes, biological pathways, diseases, and drugs. The Cluster Profiler R software package was used to analyze the GO annotation and KEGG pathway enrichment analysis of differential genes, and the critical value of False discovery rate (FDR) < 0.05 was considered statistically significant (27).
The gene set enrichment analysis (GSEA) method analyzes whether a particular gene set is statistically different between two biological states and is commonly used to estimate changes in the pathway and biological process activities in expression dataset samples (27). For GSEA analysis, the "msigdb. v7.0. Symbols" gene set was downloaded from the MSigDB database. Furthermore, the enrichment fractions of related pathways in the MSigDB database were calculated according to the gene expression matrix of each sample using the Gene Set Variation Analysis (GSVA) method by R-packet (28). The limma package was used to screen the differences, and the relevant enrichment items with statistically significant differences were displayed by a heat map. Using the R-packet GSVA and single-sample Gene Set Enrichment Analysis (ssGSEA) methods combined with the 28 characteristic gene matrices of immune cells, we calculated the immune cell enrichment fraction of each sample.

Weighted gene correlation network analysis
The aim (29) is to identify coexisting gene modules, explore the relationship between gene networks and phenotypes, and study the core genes in the network. Soft threshold was calculated by pick soft threshold function, and 5 were considered to be the optimal soft threshold. Scale-free networks were subsequently constructed based on soft thresholds, topology matrices were created, and hierarchical clustering was performed. Considering 50 as the minimum number of genes in the module, the identification gene module was dynamically cut and eigengenes were calculated. Based on eigengenes, the correlation between modules was constructed, and hierarchical clustering was carried out. The modules with a correlation above 0.7 were merged again, and lastly, 13 modules were obtained. The correlation between modules and clinical features was understood through Pearson's correlation analysis.

Construction of ceRNA network
CeRNA discloses the interaction mechanism among lncRNA, microRNA, and mRNA. The possible upstream miRNAs of TUBG1 and ABCC5 were searched through the mirTarbase database using multiMiR package (30). The StarBase database was used to download lncRNA-miRNA data (starBaseV3_hg19_CLIP-seq_LncRNA_all). Using pancancerNum >10 and clipExpNum>4 as criteria, lncRNAs are screened and intersected with those with significant differences in HCC data, and finally, four lncRNAs and 30 miRNAs are obtained. We completed the visual analysis using Cytoscape.

PPI network building (STRING)
The STRING database (https://cn.string-db.org/) searches known and predicted protein interactions across 2,031 species, comprising 9.6 million proteins and 1.38 million PPIs. It includes results obtained from experimental data, PubMed Abstracts in Chinese, and the synthesis of data from other databases, as well as results predicted using bioinformatics methods. We used the STRING database to construct the PPI network of genes with common differences between the two diseases and set the parameter as correlation coefficient 0.7. The results from PPI are exported from the STRING database and visualized through Cytoscape. Furthermore, the CytoHubba plug-in was used to analyze Hub genes in the PPI network.

Specimen collection
The study was approved by the Ethics Committee of the First Affiliated Hospital of Guangxi Medical University [NO.2022-KY-E-(115)]. NAFLD, HCC, and normal liver tissues were collected from tissue sections preserved in the Department of Pathology of our hospital. Informed consent was obtained.
High-fat LO2 cells and HepG2 cells (HF-LO2 group, HF-HepG2 group) were constructed by using sodium palmitate/ sodium oleate (PA/OA) to simulate a high-free fatty acid environment. After the growth density of HepG2 and LO2 cells reached 70%-80% according to the manufacturer's instructions, approximately 1 × 10 5 cells were inoculated in 6well plates and incubated for 24 h. The NAFLD model was established by a complete culture medium containing 500-mM sodium oleate and 250-mM sodium palmitate for 24 and 48 h. A complete culture medium supplemented with solvent was used as the control group (Control-LO2 group, Control-HepG2 group). The lipid accumulation level was detected by oil red O staining and triglyceride (TG). The TG assay kit (Jiancheng Bioengineering Institute, Nanjing, China) was utilized for extracting the contents of TG using the GPO-PAP method. According to the instruction, the cell homogenate and triglyceride assay reagents were added to the 96-well plate, mixed, and incubated at 37°C for 10 min, and the optical density (OD value) of each well was measured with the microplate reader (546 nm). The calculation formula was as follows: TG content = (sample OD value − blank OD value)/ (calibration sample OD value − blank OD value) * Calibration sample concentration/sample protein concentration.

RNA extraction and RT-QPCR
According to the manufacturer's instructions, total RNA was isolated from cultured cells using a TRIzol reagent (TAKERA, Tokyo, Japan), and cDNA was obtained by reverse transcription using RT SuperMix for qPCR (Vazyme). The extracted cDNA was stored at -80°C. The target genes were amplified and detected using SYBR qPCR Master Mix (Vazyme, Nanjing, China) in the Applied Biosystems 7500 Fast real-time PCR system (ABI, Waltham, Massachusetts, USA). Primers for ABCC5 and Glyceraldehyde-3-Phosphate Dehydrogenase (GAPDH) are as follows: ABCC5-F (5'-ATC ATG GCT TGA GTG CTC TGA-3') and ABCC5-R (5'-AGA CCA CAC GTC CAT TGA-3'); GAPDH-F (5'-AAT CAA GTG GGG CGA TGC TG-3') and GAPDH-R (5'-GCA AAT GAG CCC CAG CCT TC-3'). The setting parameters were as follows: the denaturation temperature was 95°C for 15 min, the annealing temperature was 60°C for 1 min, and the extension was 40 cycles. The relative mRNA expression of ABCC5 was calculated by the 2 -DDCT method (GAPDH as a housekeeping gene).

Immunohistochemistry
Samples were prepared into paraffin sections, followed by dewaxing and hydration, and high-pressure and hightemperature antigen recovery in sodium citrate buffer (pH 6.0). Subsequently, it was incubated with endogenous peroxidase blockers and then incubated with a primary antibody for 12 h at 4°C. Sections were then rinsed with phosphate-buffered saline (PBS), followed by washing with an enzyme-coupled goat anti-rabbit/mouse IgG polymer for another 20 min and rinsed with PBS three times. Then, the sections were incubated with a DAB staining solution for 5 min. After washing with water, the paraffin sections were immersed in a hematoxylin-staining solution for 15 s and finally identified, rinsed, and returned to the laboratory for staining. The antibodies used included ABCC5 (Cat No. 19503-1-AP, 1:100; Proteintech, Wuhan, China) (31) and TUBG1 (Cat No. 15176-1-AP, 1:100; Proteintech, Wuhan, China) (32). Two pathologists observed the results. The number of positive cells was divided into five grades (0-4 score) according to the ratio (<5%, 5%-25%, 26%-50%, 51%-75%, 75%-100%), and the intensity of staining was divided into four grades (0-3 score). The group is obtained by multiplying the two fractions (from 0 to 12). A score of >4 was suggested to be positive.

Statistical analysis
We performed all data calculations and statistical analysis using R programming (https://www.r-project.org/, version 4.0.2). The statistical significance of the normally distributed variables was calculated using an independent Student's t-test, and the differences between the non-normally distributed variables were analyzed using the Mann-Whitney U test (i.e., Wilcoxon rank-sum test). All statistical p-values were bilateral, and P <0.05 was considered statistically significant.

Results
Weighted gene correlation network analysis consensus module identification NAFLD can ultimately cause HCC, but the mechanism by which hepatitis leads to HCC remains unclear. We assumed that the same gene modules expressed in both diseases might play a role in the copathogenesis of both diseases. NAFLD datasets GSE48452 and GSE89632 were downloaded from the GEO database and HCC dataset from the TCGA database. WGCNA analysis was used to identify gene modules that were expressed uniformly in both NAFLD and HCC. We removed large outliers with hierarchical clustering analysis ( Figure 1A), scale-free networks and topological matrices were constructed by soft threshold ( Figure 1B), and then topological overlapping matrices were scaled to make them comparable between two disease sets ( Figure 1C). By extracting the minimum of parallel positions in two topological overlapping matrices [we now calculate the consensus topological overlap by taking the component-wise ["Parallel")], the minimum of the topological matrices in individual sets acted as the consensus topological overlap matrix. Finally, we obtained 13 gene modules with similar expression in both diseases ( Figure 1D). As clinical features, we explored modules related to the diagnosis of both diseases based on whether they were sick or not. The dark turquoise module had a high correlation and statistical significance with NAFLD and HCC among the 13 modules (Figures 2A-C). Then, we performed GO and KEGG function analysis on 58 genes in the dark turquoise module, and the results indicated that these genes were primarily involved in DNA-related enzyme activity, DNA shear, and the repair process ( Figures 2D, E).

Differences in gene
We performed difference analysis on the two datasets to understand the correlation of transcriptional level changes in the two diseases. Using principal component analysis (PCA), we found some differences in expression profiles between the NAFLD group and the normal group. We used the limma package for difference analysis to obtain DEGs in the two data groups, including 181 upregulated genes and 232 downregulated genes, which were visualized by the volcanic diagram and heat maps ( Figures 3A-C). Significant differences were found between HCC and the control group gene expression profiles. We used the Deseq2 package to analyze the differences of genes in different groups, comprising 2,490 upregulated genes and 1,077 downregulated genes, which were also used and visualized by the volcano map and heat map ( Figures 3D-F).
Subsequently, we searched for both upregulated and downregulated genes in the two diseases and demonstrated them through the Venn diagram. The results showed that 27 genes were upregulated in the two diseases compared with the control group, while 83 genes were downregulated in the two diseases compared with the control group ( Figure 4A). We performed the GO function analysis of genes to further understand the biological function of these genes with same trend. The results indicated that these genes were mainly involved in cell chemotaxis, neuroresponse, response to lipopolysaccharide, response to the molecule of bacterial origin, response to steroid hormone; a variety of lipid-related functional items are involved, and then we used the mesh diagram to visualize the main items ( Figure 4B; Supplement Table 1). Then, we used GSEA enrichment analysis to enrich the pathway-related dataset in MSigDb based on the above gene list. The results indicated that it was associated with multiple immune, HCC, and cancer-related pathways (Figures 4C, D; Supplement Table 2). Consequently, we constructed the PPI network through STRING (correlation coefficient 0.7) and visualized it using Cytoscape ( Figure 4E). The top three genes interacting with other proteins were EGR1, FOSB, and CCL2.  To construct a diagnostic model of nonalcoholic fatty liver disease The consensus module genes related to the diagnosis of the two diseases obtained by WGCNA were intersected with the different genes (DEG) with the same change trend in the two diseases to find efficient diagnostic and therapeutic markers ( Figure 5A), suggesting three coexisting genes (ABCC5, DHODH, and TUBG1). After univariate and multivariate logistic regression analyses, ABCC5 and TUBG1 genes were found to be independent risk factors for the diagnosis of NAFLD and HCC (OR>1; P <0.05), while DHODH was an independent protective factor in HCC (OR<1; P < 0, 05) ( Tables 1, Table 2). These findings propose that ABCC5 and TUBG1 may have certain roles in the occurrence and development of the two diseases. We visualized the expression of the two genes by a heat map in the group of NAFLD ( Figure 5B). The above three genes were reverified using minimal absolute contraction and selection operator LASSO regression. The results were in line with multivariate logistic regression. Two genes, ABCC5 and TUBG1, were involved in constructing a diagnostic model of NAFLD (Figures 5C, D). Its correlation coefficient was used to calculate the risk values of all samples, and a box plot was Identification of diagnostic-related modules. (A) Correlation analysis was used to find the correlation between 13 gene modules in NAFLD and whether they were sick or not. Red indicated positive correlation, blue indicated negative correlation, and P<0.05 was considered statistically significant. (B) Correlation analysis was used to find the correlation between 13 gene modules in HCC and whether they were sick or not. Red indicated positive correlation, blue indicated negative correlation, and P<0.05 was considered statistically significant. used to visualize the difference of risk values between the NAFLD group and the normal group, indicating that the risk values of NAFLD were significantly higher ( Figure 5E). According to the ROC curve results, the diagnostic model constructed by ABCC5 and TUBG1 genes had very good diagnostic efficiency (AUC = 0.893) ( Figure 5F). As further validation of our results, we downloaded the NAFLD dataset GSE37031 through the GEO database. The results showed that the risk value of the NAFLD group was significantly higher, and the TUBG1 expression level was also higher in the NAFLD group (P<0.01) (Figures 5G, H). While the box plot indicated that the ABCC5 expression level was high in the NAFLD group (P = 0.064), the p-value was insignificant as well as the AUC value of the ROC curve was 1, which may result from the small sample size. (Figure 5I, Supplement Figure 1). Results from these studies suggest that the high Variance analysis. (A) Volcanic diagram shows the differentially expressed genes (DEGs) of the NAFLD gene dataset, the horizontal axis is logFC value, the y-axis is -log10 (adjusted p-value), and the dotted line in the figure representslogFC> 0.5, adj-P < 0.05). Red represents upregulated genes, and blue represents downregulated genes. The first five gene labels are displayed in the adjusted p-value arrangement. (B) PCA principal component analysis of the NAFLD dataset: light green represents the samples of NAFLD, and red represents the samples of normal liver tissue (normal). (C) Heat map of differential genes. Red represents the sample of NAFLD, and dark blue represents the sample of normal liver tissue (normal). The top 20 genes are, respectively, visualized in adjusted p-value arrangement (yellow indicates high expression, and blue indicates low expression). (D) Volcanic diagram shows the DEGs of the HCC gene dataset; the horizontal axis is logFC value; the y-axis is -log10 (adjusted p-value), and the dotted line in the figure representslogFC> 1, adj. P < 0.05). Red represents upregulated genes, and blue represents downregulated genes. The first five gene labels are displayed in the adjusted p-value arrangement. (E) Principal component analysis (PCA) of HCC datasets: light green represents HCC tissue samples (cancer), and red represents normal liver tissue samples (normal). (F) Heat map of differential genes. Red represents HCC tissue samples (cancer), and dark blue represents normal liver tissue samples (normal). The top 20 genes are, respectively, visualized by adjusted p-value arrangement (red indicates high expression, and dark green indicates low expression). expression of ABCC5 and TUBG1 is an independent risk factor for NAFLD and can be used as a diagnostic indicator.

ABCC5 and TUBG1 are independent risk factors for diagnosis of hepatocellular carcinoma
By performing univariate and multivariate logistic regression analysis on three coexisting genes (ABCC5, DHODH, and TUBG1) in TCGA HCC dataset, it was found that ABCC5 and TUBG1 were independent risk factors for HCC (OR>1; P <0.05), whereas DHODH was an independent protective factor in HCC (OR<1; < 0, 05) ( Table 2). To further confirm our conclusions, we used HCC dataset (LIRI-JP) from the ICGC database for validation and found ABCC5 and TUBG1 genes were independent risk factors for HCC (OR>1; P <0,05), while DHODH was not (Table 3). These outcomes showed that ABCC5 and TUBG1 were independent risk factors for NAFLD and HCC, suggesting that the two genes possessed potential diagnostic value in both diseases. Common differential genes. (A) Venn diagram shows the coupregulated and codownregulated genes in the two diseases. (B) Network diagram showing the GO function analysis of common differential genes. (C) Gene set enrichment analysis (GSEA) of common differential genes: the logFC value was used to arrange genes, and the GSEA algorithm was used to calculate the first two enrichment items and the last two enrichment items in descending order. (D) GSEA analysis: the top 10 enrichment pathway were visualized in the bubble diagram, arranged with the descending order of the absolute value of enrichment fraction (NES). (E) PPI network construction, from red to blue; continuous color change represents the logFC value from high to low.

Construction of prognostic model in hepatocellular carcinoma
Univariate and multivariate Cox regression analyses were performed to explore whether three coexisting genes (ABCC5, DHODH, and TUBG1) were associated with the prognosis of HCC, and the results showed that ABCC5 and TUBG1 were independent risk factors for the prognosis of HCC (HR>1, P<0.05) ( Figure 6A). Subsequently, we calculated the risk value of each sample based on the correlation coefficient. There were  more deaths, shorter survival time, and higher expression levels of ABCC5 and TUBG1 in the high-risk group (Figures 6B, C). Afterward, we verified our prognostic model using ROC curves with survival time nodes of 1, 2, and 3 years. The results presented the AUC value of 1-year, 2-year, and 3-year ROC curves to be 0.80, 0.68, and 0.68, respectively ( Figure 6D). As shown in Figure 6E, the risk value was higher in the death group (P<0.05). Using the nomogram, we found that when combined with age, gender, stage, and the depth of invasion, the risk score significantly contributed to survival risk ( Figure 7A). We  Figure 7B). Next, we tested the predictive accuracy of the model using the calibration curve, and the results showed that the prediction accuracy of 1-year, 2-year and 3-year survival rates was very high ( Figure 7C). The survival curve showed that both ABCC5 and TUBG1 were significantly correlated with poor prognosis in HCC and the survival rate of patients in the highrisk group was significantly lower than that in the low-risk group ( Figures 7D, E). So as to further verify our prognostic model, HCC-related data (LIRI-JP) from the ICGC database were used for confirmation. The ABCC5 and TUBG1 expressions were significantly increased in the HCC group. In the meantime, based on the model, it was determined that the HCC group had a significantly higher risk score than the normal group ( Figures 8A, B). Then, we verified our prognostic model by using ROC curves with survival time nodes of 1, 2, and 3 years. The results disclosed that the AUC of 1-year, 2-year and 3-year ROC curves was 0.70, 0.66, and 0.74, respectively ( Figure 8C). Meanwhile, the survival curve showed that high risk score was significantly correlated with poorPrognosis in LIRI-JP ( Figure 8D).

Changes of immune microenvironment and construction of ceRNA network
We used the GSVA algorithm to analyze pathway-related datasets in MSigDb in the two diseases and the limma package for difference analysis. LogFC = 0.5, ADj. P<0.05 were used as thresholds to screen related pathways with significant differences between the disease and normal groups. The results showed 46 significantly upregulated and 122 significantly downregulated pathways in NAFLD and 1,418 significantly upregulated and 1,669 significantly downregulated pathways in HCC. We visualized the top 10 pathways with significant changes using heat maps (Figures 9A, B) and then obtained both upregulated and downregulated pathways in the two diseases ( Figure 9C, Table 4). The outcomes showed that NAFLD and HCC had different fat metabolism, immunity, tumor, and other related pathways as compared to the normal group. These pathways may disclose the principal mechanisms by which NAFLD progresses to HCC. We conducted enrichment analysis of 28 immune cell components in the two datasets ssGSEA to understand the immune infiltration of the two diseases and found that the immune microenvironment of the two diseases had significantly changed. Furthermore, effector memory CD4 T cells, eosinophils, mast cells, neutrophils, and type 1 T helper cells showed significant differences and the same trend in the high-and low-risk groups of the two diseases ( Figure 9D). CeRNA discloses the interaction mechanism among lncRNA, microRNA, and mRNA. The upstream miRNAs and lncRNAs of ABCC5 and TUBG1 were predicted through the ceRNA network, and visualization analysis was executed using  Cytoscape ( Figure 9E). According to the ceRNA mechanism, lncRNAs can upregulate the expression of their target mRNA. Hence, we used the upregulated lncRNAs expressed in HCC data based on TCGA and intersected with the predicted lncRNAs to find four lncRNAs with high confidence and consistent expression.

Expression of ABCC5 and TUBG1 at tissue and cell levels
Subsequently, the expression of ABCC5 and TUBG1 in normal liver tissue, NAFLD tissue, and HCC tissue was detected by immunohistochemistry. The results showed that compared with normal liver tissue, the protein expression levels of ABCC5 and TUBG1 in NAFLD and HCC tissues were significantly increased (P <0.05) (Figures 10A, B).
In order to further explore the expression of ABCC5 and TUBG1 in NAFLD and HCC. We simulated a high-fatty-acid environment by mixing PA and OA (33). HF-LO2 group and HF-HepG2 group models were constructed (34-36). Oil red O staining results showed that OA significantly increased lipid accumulation in HepG2 cells and LO2 cells compared with the control group (Figures 10C, D). The results of the quantitative analysis of triglyceride were consistent with the above results, and the triglyceride content of cells in HF-LO2 and HF-HepG2 Construction of hepatocellular carcinoma prognostic model. (A) Univariate Cox analysis suggested that ABCC5 and TUBG1 were correlated with the prognosis of HCC. Multivariate Cox analysis suggested that ABCC5 and TUBG1 were independent risk factors (red represents HR>1; dark blue represents HR<1, P<0.05 is statistically significant). (B) Both ABCC5 and TUBG were upregulated in the high-risk group (red represents high-risk group; dark blue represents low-risk group). (C) The risk value was calculated by multivariate Cox regression. With the increase of the risk value, the survival time decreased, and the number of patients who died increased. (D) With risk score as the observed value, the time ROC curve confirmed that it had good diagnostic value at 1-year, 2-year, and 3-year survival nodes (the AUC of a 1-year ROC curve was 0.8, the 2year ROC curve was 0.68, and the 3-year ROC curve was 0.68). (E) The box plot shows that the death group has higher risk scores (0 represents alive group; 1 represents dead group).
groups was significantly higher than that in the corresponding control group ( Figure 10E). The above results indicated that the adipose cell model was successfully constructed.
The transcriptional levels of ABCC5 and TUBG1 in HepG2 cells and LO2 cells were detected by RT-QPCR. The results showed that ABCC5 expression was significantly increased in HepG2 cells compared with LO2 cells (P <0.05). Similarly, RT-QPCR results showed that ABCC5 expression was significantly increased in both the HF-LO2 group and the HF-HepG2 group after aliphication, and the increased level increased with the extension of aliphication time (P <0.05) ( Figure 10F). However, the transcription level of TUBG1 was not significantly changed in HepG2 and LO2 cell lines and before and after adipose treatment.

Discussion
The annual global increase in the HCC incidence is inseparable from NAFLD, an independent risk factor. In the wake of increasing living standards, the incidence of NAFLD is increasing at an alarming speed. The prevalence of NAFLD Risk model testing. (A) Nomogram recommends that risk score would contribute the most than other clinical characteristics (age, gender, stage, and the depth of invasion). (B) Multivariate Cox regression analysis of the risk score combined with clinical characteristics (age, gender, stage, and the depth of invasion) suggested that the risk score was still an independent risk factor. (C) The results of 1-year, 2-year, and 3-year calibration curves show that the prediction accuracy of risk models is high. (D) The survival curve suggested that patients in both ABCC5 and TUBG1 high-expression groups had a poor prognosis (red represents a high-expression group; yellow represents low-expression group). (E) The survival curve suggested that patients in the high-risk group had a poor prognosis (red represents a high-risk group; yellow represents low-risk group, P<0.05 is statistically significant).
is increasing at an alarming rate and has been considered an ind uci ble c au se of HCC ( 3 7 ), a n d n on-a lc oholic steatohepatitis (NASH), as a manifestation of NAFLD, can progress to cirrhosis or even HCC (38). Elucidating the mechanism of NAFLD to HCC will enable us to have the opportunity to take measures at the early stage, so as to achieve early prevention and reduce the incidence of HCC, which will be of great significance. Currently, some excellent studies have tried to explore the underlying mechanism of the progression of NAFLD to HCC. Dabin Liu (39)  Risk model was tested with the ICGC-HCC dataset. (A) The box plot shows that the expression of ABCC5 and TUBG in the ICGC-HCC dataset weas significantly higher in the HCC group. (B) The box plot shows that in the ICGC-HCC dataset, the risk score in the HCC group was significantly higher than that in the normal group. (C) With the risk score as the observed value, the timeROC curve confirmed that our prognostic model still had a high diagnostic value (the AUC of the 1-year ROC curve was 0.7, the 2-year ROC curve was 0.66, and the 3-year ROC curve was 0.74). (D) The survival curve suggested that patients in the high-risk group had a poor prognosis (red represents high-risk group; yellow represents low-risk group, P<0.05 is statistically significant).
targeting SQLE. However, the pathogenesis of HCC induced by NAFLD still stays uncleared. Therefore, we explored the potential mechanism of NAFLD to HCC and found two key genes. In this study, 26 upregulated and 87 downregulated genes were obtained by limma package analysis using the NAFLD dataset (GSE48452, GSE89632, GSE37031) and HCC dataset TCGA-LIHC (n = 445 cases). By using the limma package to analyze the differences in the data of the two diseases, 26 coupregulated genes and 87 codownregulated genes were obtained. Then, we used 13 consensus modules screened by WGCNA from the two diseases, among which the darkturquoise module was highly correlated with the diagnosis of both diseases. The intersection of the dark turquoise module and common downregulated genes was obtained, and the coexpressed genes ABCC5, DHODH, and TUBG1, which were expressed uniformly in both diseases and were related to diagnosis, were finally obtained. After the establishment of diagnostic and prognostic models and the verification process of prediction ability, it was concluded Immune infiltration and ceRNA. (A) Heat map showing pathways that differ significantly between the normal and NAFLD groups (red for NAFLD group and dark blue for normal group). (B) Heat map shows pathways with significant differences between the normal group and the HCC group (red represents the cancer group; dark blue represents the normal group). that the expression patterns of ABCC5 and TUBG1 in NAFLD and HCC were consistent, which were common diagnostic markers and highly correlated with the prognosis of HCC. Both may be closely related to NAFLD-induced HCC and may become new biomarkers for the biological diagnosis and prognosis of HCC. Current studies indicate that ABCC5 plays an oncogene role in HCC. Borel et al. (40) studied the expression of 15 ABC transporters related to multidrug resistance in 19 samples of HCC patients (16 untreated and 3 treated with chemotherapy); 12 ABC transporters were found to be upregulated in HCC as compared to healthy livers, and one of them is ABCC5. Chen has found that (41) ABCC5 levels are significantly increased in HCC, and it is considered to be an independent risk factor for HCC progression and prognosis mainly because of immune cell infiltration. Due to its drug resistance characteristics, Huang et al. (42) proved that the regulation of the PI3K/AKT/NRF2 pathway activates and upregulates ABCC5 to promote the resistance of sorafenib (a first-line molecular targeted drug for advanced HCC) acquired in human HCC cells. This result suggests that ABCC5 is an important regulator and a promising therapeutic target of sorafenib resistance acquired in human HCC cells. ABCC5 has also been found to be able to regulate and reduce GPx4 consumption in sorafenib-resistant HCC cells to inhibit lipid peroxidation. Additionally, ABC transporters mainly mediate the transmembrane transport of cyclic nucleotides or nucleotide analogs, able to transport organic anionic compounds, including glucoside acids, and the sulfate conjugates of steroids (43). Further studies (44) have found that it actively transports bile acids made by cholesterol in liver cells. Lyu (45), Liu (46), and Akbulut (47)  The results of our bioinformatics analysis are consistent with the above studies. ABCC5 is a diagnostic marker in both diseases and is associated with the prognosis of HCC. In addition, it was found for the first time in vitro that the expression level of ABCC5 was significantly increased in hepatocellular carcinoma after adipose transformation, and the increased level was positively correlated with the degree of adipose transformation.
TUBG1 is a member of the human tubulin family and is involved in the development of various tumors. In the study of Hubert (49) and Hsu (50) et al. it is found the breast Cancer 1 protein (BRCA1), an inhibitor of breast and ovarian tumorigenesis, maintained a high fidelity of cell division after centrosome binding to TUBG1 produced a ubiquitination reaction. Furthermore, its association with tumor suppressor P53 and centrosome was explained by Morris (51) and Kanai (52) et al. It was also revealed by Horějsǐ́et al. that high expression of TUBG1 antagonized its inhibitory effect on DNA damage (53). In a hospital-based case-control study (54), variations in TUBG1 were also associated with breast cancer risk. Blanco et al. (55) further explored that TUBG1 interacts with associated proteins in sporadic breast tumors to regulate the mammary epithelial polarization and affect patient survival. However, the relationship between TUBG1 and the progression of HCC and NAFLD patients is still unclear, and few studies have been conducted in this context. Through bioinformatics multiple analysis, we believe that this gene may be a common carcinogenic factor in the pathogenesis of NAFLD and HCC. Our immunohistochemical results were consistent with the results of the bioassay, and the TUBG1 expression level was significantly increased in both NAFLD and HCC tissues. However, there was no significant change in the transcription level after fatty treatment. We speculated that the increase of the TUBG1 protein level may not be caused by its own increase but rather due to changes in the microenvironment of liver cells under a high-lipid environment. We believe that TUBG1 is a very promising research site. Although we verified the above results by using a variety of methods, there are still limitations. The regulatory mechanism of ABCC5 and TUBG1 in the development of HCC in the context of NAFLD is not detailed yet, and whether the differences of the last two genes in different stages of NAFLD-HCC progression also needs further study and confirmation.
Our results showed that ABCC5 and TUBG1 were increased in both NAFLD and HCC and were independent risk factors. The diagnostic model of NAFLD and the prognostic model of