We downloaded the raw data from TCGA (https://portal.gdc.cancer.gov/), ICGC (https://icgc.org/), and GEO (https://www.ncbi.nlm.nih.gov/geo/) website, processed the raw data, and finally obtained the expression matrix of hepatocellular carcinoma tissues (HCC). Datasets GSE22058, GSE25097, GSE36376, and GSE39791 containing data on HCC and para-cancerous tissues were used to screen differentially expressed genes (DEGs). All datasets were normalized by Robust Multi-Array Average (RMA), and different probes were converted into gene symbols. The dbEMT2 database (http://dbemt.bioinfo-minzhao.org/download.cgi) was used to download EMT-related genes. Perl language was used to extract the expression matrix of EMT-related genes from the TCGA, ICGC, and GSE14520 HCC datasets in batches.

The MetaDE package was used to screen DEGs between HCC and para-cancerous tissues (Ramasamy et al., 2008; Wang et al., 2012). Before performing the microarray meta-analysis, we extracted common genes between different datasets and eliminated the top 30% genes with average signal intensity and standard deviation. The four methods Fisher, maxP, roP, and AW were used to analyze differentially expressed genes (DEGs). A corrected p-value of less than 0.05 was statistically significant. Given the heterogeneity among different statistical methods, the R package VennDiagram was used to identify common differentially expressed EMT-correlated genes.

The clusterProfiler package was used to identify the biological functions and pathways of differentially expressed EMT-related genes (Yu et al., 2012). Gene Ontology (GO) annotation was conducted for three domains: cellular component (CC), molecular function (MF), and biological process (BP). A bubble plot was used to visually display the results of GO and KEGG function analyses. Gene set enrichment analysis (GSEA) was used to evaluate the distribution of genes in a predefined gene set. Results were displayed in a gene table with genes ranked by their relevance to phenotype, reflecting their contribution to the phenotype. The input data consisted of two parts, a gene set with known functions and an expression profile matrix, to explore significantly enriched pathways between the high and the low EMT risk score group. An FDR and p-value of less than 0.05 were statistically significant.

After excluding HCC patients who survived less than 30 days, 342, 221, and 230 patients remained in the TCGA, GSE14520, and ICGC datasets, respectively. The prognostic-related EMT genes were determined in the TCGA and GSE14520 datasets, yielding 43 common EMT-related prognostic genes based on the univariate Cox regression model. Then, we used the least absolute shrinkage and selection operator (LASSO) cox regression model to further screen the EMT-related risk genes in the TCGA dataset and finally determined 9 EMT gene signatures could predict hepatocellular carcinoma prognosis. The EMT risk score is equal to the sum of the product of EMT gene expression and its coefficient. We calculated the EMT risk score of each patient in the TCGA, GSE14520, and ICGC datasets and divided HCC patients into a high and a low EMT risk score group based on the median value in each dataset. The TCGA dataset was used to construct the EMT gene signature, while the GSE14520 and ICGC datasets were used to validate the constructed signature. The survival curve between high and low-risk groups of HCC patients was generated using the survminer package, and the difference in survival times was assessed by the log-rank test. The survival ROC package assessed the predictive performance of the EMT risk score for survival outcomes.

The nomogram was based on multivariate regression analysis, which integrates multiple predictors, and scaled line segments were drawn on the same plane according to a certain proportion to ensure that the survival probability of patients can be predicted objectively. In the TCGA dataset, we included patients with exact age, sex, tumor differentiation, tumor stage, and survival information. Univariate and multivariate cox regression models were used to screen variables that affect survival outcomes (including OS and RFS). Therefore, the risk score and tumor stage data from the TCGA dataset were included to build a predictive nomogram model for OS and RFS. Datasets GSE14520 and ICGC were used for external validation. The time-dependent ROC curve was used to assess the predictive ability, and the consistency of the nomogram model was visually represented by the 1-year, 3-year, and 5-year calibration curve. The larger the area under the time-dependent ROC curve, the better the discriminative ability of the nomogram model. The higher the degree of coincidence between the actual 1-year, 3-year, and 5-year calibration curves and the lines at the 45-degree angle, the better the consistency.

The Human Protein Atlas (http://www.proteinatlas.org/) database is a well-recognized online database that provides tissue and cell distribution information of all 24,000 human proteins. It can provide the expression of a specific protein in tumor tissues and normal tissues. In the present study, we analyzed differences in protein expression of these 9 EMT risk genes between HCC and normal liver tissues.

The TIMER (Tumor Immune Estimation Resource) database (https://cistrome.shinyapps.io/timer/) is an online database that provides information on the infiltration status of the six immune cells (B cells, CD4⁺ T cells, CD8⁺ T cells, Neutrophils, Macrophages, and Dendritic cells) in the TCGA database (Li et al., 2017). Besides, we downloaded the immune score, stromal score, and ESTIMATE score of HCC tissues of the TCGA dataset from the ESTIMATE database (https://bioinformatics.mdanderson.org/estimate). Correlation analysis was used to explore the correlation between EMT risk score and immune-related indexes (including six immune cell types, immune score, stromal score, and ESTIMATE score). Immune checkpoint genes included CD274 (PD-L1), PDCD1(PD-1), TIGIT, CTLA4, HAVCR2, and ICOS. Finally, we analyzed the differences in mRNA expression levels of immune checkpoint genes between high and low-risk groups.

All statistical analyses were conducted using R3.6.1 and SPSS26.0 software. Chi-square analysis was used to analyze the correlation between the EMT risk score group and clinicopathological characteristics of HCC patients. EMT genes with a p-value less than 0.05 during univariate analysis were included in the LASSO Cox regression model. STATA software was used to explore the combined effect between EMT risk score and survival outcome. A p-value less than 0.05 was statistically significant.

4 HCC cohorts containing HCC and para-cancerous tissues were integrated. The details of the four cohorts are shown in Table 1. The detailed process of this research is displayed in Figure 1. The AW, Fisher, maxP, and roP methods yielded 442, 442, 318, and 318 EMT-related differentially expressed genes, respectively (Figure 2A). Figure 2B shows the expression matrix heat map of the EMT related-genes with FDR values less than 1E-19 in the four HCC cohorts using the maxP method. Finally, we identified 317 differentially common expressed EMT-related genes (Figure 2C).

The clusterProfiler package was used to explore the functions and pathways of differentially expressed EMT-related genes. The top 5 significantly enriched GO terms in the BP domain included gland development, epithelial cell proliferation, regulation of vasculature development, neuron death, and response to transforming growth factor-beta (Figure 3A). The top 5 significant GO CC terms included focal adhesion, cell-substrate adherens junction, cell-substrate junction, collagen-containing extracellular matrix, and cell leading edge (Figure 3A). Moreover, the top 5 significantly enriched MF GO terms were RNA polymerase II transcription factor binding, cell adhesion molecule binding, beta-catenin binding, DNA-binding transcription activator activity, and repressing transcription factor binding (Figure 3A). Finally, the differentially expressed EMT-related genes were significantly enriched in KEGG pathways, including MicroRNAs in cancer, Human papillomavirus infection, Human cytomegalovirus infection, Shigellosis, and Chronic myeloid leukemia (top 5) (Figure 3B).

We used univariate Cox regression models to analyze the correlation between EMT-related genes and overall survival in the TCGA and GSE14520 datasets to establish a robust EMT-related gene signature. In the TCGA and GSE14520 datasets, we identified 131 and 68 EMT prognostic-related genes, respectively, and finally got 43 common prognostic EMT-related genes. Based on the TCGA dataset, the Lasso-Cox regression model identified a 9 EMT-related gene signature through 10-fold cross-validation (Figures 4A,B). The coefficients of the screened EMT genes are shown in Figure 4C. The formula for calculating the EMT risk score is as follows: 0.3084* HDAC2 + 0.0364* SPP1−0.0890* PPARGC1A+ 0.0407* LGALS3 + 0.1885* ENO1+ 0.0322* LMNB1+ 0.0119* CKS2+ 0.0143* SERPINE1+ 0.1062* HDAC1. According to the median EMT risk scores of HCC patients in the TCGA, GSE14520, and ICGC datasets, we divided HCC patients into low-risk and high-risk groups. In the TCGA dataset, the AUC values of the EMT risk score for predicting the overall survival of HCC patients at 1, 3, and 5 years were 0.792, 0.714, and 0.702, respectively (Figure 5A). HCC patients in the low EMT risk score group exhibited a better overall survival rate than HCC patients in the high EMT risk score group (p < 0.001) (Figure 5D). Similarly, in the GSE14520 dataset, the predictive performance of the EMT risk score for the overall survival of HCC patients at 1, 3, and 5 years was 0.703, 0.682, and 0.706 (Figure 5B). HCC patients with low EMT risk scores were associated with a longer overall survival time than HCC patients with high EMT risk scores (p < 0.001) (Figure 5E). In the ICGC validation set, the EMT risk score exhibited a better predictive power for overall survival (Figure 5C), and HCC patients with a high EMT risk score were associated with worse overall survival than patients with low EMT risk scores (Figure 5F). Finally, we explored the predictive power of the EMT risk score in the TCGA dataset and the GSE14520 dataset for HCC recurrence. The results of the time-dependent ROC curve showed that the EMT risk score has a moderate diagnostic value in the TCGA and GSE14520 datasets for predicting RFS (Figures 6A,B). The survival analysis results showed that HCC patients with high EMT risk scores had shorter recurrence-free time than HCC patients with low EMT risk scores (Figures 6C,D). Supplementary Figure S1 shows the distribution of EMT risk scores (A, B, C), patient survival status (D, E, F), and the expression heat map of the selected EMT genes (G, H, I) in the TCGA, GSE14520, and ICGC datasets.

Clinicopathological characteristics associated with the TCGA, GSE14520, and ICGC datasets, including age, gender, ALT, history of malignancy, Histopathological grade, tumor stage, and survival status, were extracted. In the TCGA dataset, a significant correlation was found between the EMT risk score and histopathological grade (p < 0.001), tumor stage (p < 0.001), and survival status (p < 0.001). In the GSE14520 dataset, the EMT risk score was significantly correlated to age (p = 0.043), TNM stage (p < 0.001), and survival status (p < 0.001). In the ICGC dataset, the EMT risk score was significantly correlated to the TNM stage (p < 0.001) and survival status (p = 0.001) (Table 2).

In the TCGA dataset, the EMT risk score (OS: p < 0.001; RFS: p = 0.005) and tumor stage (OS: p < 0.001; RFS: p < 0.001) were independent risk factors that affected the overall survival and recurrence-free survival of HCC patients (Supplementary Table S1, S2). In the GSE14520 dataset, the EMT risk score (p < 0.001) and tumor stage (p < 0.001) were independent risk factors that affected the overall survival of HCC patients (Supplementary Table S3). Moreover, the tumor stage (p = 0.001) was an independent predictor of recurrence-free survival in HCC patients (Supplementary Table S4). In the ICGC dataset, the EMT risk score (p = 0.002) was an independent predictor of overall survival in HCC patients (Supplementary Table S5).

In the present study, 3 OS and 2 RFS datasets were analyzed. We integrated the hazard ratio in the multivariate Cox regression analysis of the EMT risk score in the OS and RFS datasets. The results of the meta-analysis indicated that a high EMT risk score was a risk factor for overall survival (HR: 2.57; 95% CI: 1.94–3.41, p < 0.001) and recurrence-free survival (HR: 1.47; 95% CI: 1.15–1.88, p = 0.002) of HCC patients (Supplementary Figure S2).

Nomograms were established to objectively predict the survival outcome (including OS and RFS) of HCC patients at 1, 3, and 5 years. In the TCGA dataset, the EMT risk score and tumor stage were independent risk factors that affected the overall survival and recurrence-free survival of HCC patients. Therefore, OS and RFS nomogram models were constructed based on these two variables (Figure 7). The results of the time-dependent ROC curve showed that the OS and RFS datasets-based nomogram model yielded the best predictive performance compared with other clinicopathological variables (Figure 8). The 1-year, 3-year, and 5-year calibration curves showed that the constructed nomogram model had good accuracy in predicting OS (Figure 9) and RFS (Figure 10).

We analyzed the mRNA expression of 9 EMT-related genes in the TCGA, GSE14520, and ICGC datasets and found that these nine genes were dysregulated in these three HCC datasets (Supplementary Figures S3A–C). The HPA database was used to explore the protein expression levels of EMT-related genes. No corresponding data was found for PPARGC1A in the database. The differences in expression of the remaining eight genes between HCC tissue and normal liver tissue is shown in Supplementary Figure S3D. The results showed that HDAC2, HDAC1, SPP1, CKS2, and LGALS3 were significantly increased in HCC tissue compared with normal liver tissue. No significant difference in ENO1, SERPINE1 LMNB1 expression was found between normal liver and HCC tissues.

GSEA was used to explore significantly enriched KEGG pathways associated with the EMT gene signature in the TCGA dataset. The high EMT risk group was significantly enriched in the spliceosome, RNA degradation, oocyte meiosis, pyrimidine metabolism, and cell cycle (Top 5) (Figure 11A). In contrast, the low EMT risk group was significantly enriched in Drug metabolism cytochrome P450, fatty acid metabolism, Glycine serine and threonine metabolism, retinol metabolism, and Primary bile acid Biosynthesis (Top 5) (Figure 11B).

We downloaded the expression of tumor-infiltrating immune cells, tumor immune score, stroma score, and ESTIMATE score from the TIMER and ESTIMATE databases. Correlation analysis showed a significant correlation between the EMT risk score and six tumor immune cells in the TCGA dataset (Supplementary Figure S4). Moreover, a significant correlation was found between the EMT risk score and immune score (Supplementary Figure S5), further substantiating the relationship between the EMT risk score and immune infiltration levels. Immunotherapy has gained significant momentum in recent years, especially in tumor therapy. Mounting evidence suggests that immunotherapy has significantly prolonged the survival of many cancer patients. In our study, the differences in gene expression in immune checkpoints and T cell exhaustion were analyzed between the high-risk and low-risk groups. The results showed that HCC patients in the high-risk group had higher mRNA expression of immune-related genes (Figure 12A). Moreover, the ESTIMATE score, immune score, B cell, CD4T cell, CD8T cell, Dendritic cell, Macrophage cell, Neutrophil were significantly higher in the high-risk group than in the low-risk group (Figures 12B,C).