Clinical data and RNA sequencing of HCC were downloaded from the LIHC project of The Cancer Genome Atlas (TCGA) database (https://cancergenome.nih.gov), and a total of 50 corresponding paracancerous tissues and 365 HCC tumor tissues were incorporated into the analysis. A set of HCC chip data numbered GSE14520 were collected from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/). Another set of data were downloaded from the HCCDB database (http://lifeome.net/database/hccdb/download.html), which provides expression profiles of HCC samples. Meanwhile, pan-cancer cell line drug sensitivity and genomic data resources were acquired from the Genomics of Drug Sensitivity in Cancer (GDSC) database (https://www.cancerrxgene.org/) (Yang et al., 2013).

WGCNA was performed for analyzing the expression profiles of cancer cell lines downloaded from GDSC2. Samples were clustered and used to construct a gene co-expression network, from which modules were identified and then related to external data. Key drivers in the “WGCNA” package (Langfelder and Horvath, 2008) of R were analyzed based on the relationships among the module. Under the selected parameters, the sampleTree function provided by “WGCNA” was used to cluster the HCC cell lines and present the outliers. The gene expression matrices of the retained samples were extracted, and Pearson correlation was computed to calculate the correlation between twisted genes. Under different power values, we selected the optimal soft threshold β by analyzing the scale independence and average connectivity of modules using the “pickSoftThreshold” function provided by the “WGCNA” package. The “scaleFreePlot” function was adopted to evaluate whether the topology of the network was scale-free. The hierarchical clustering of genes was implemented by using the “hclust” function. The distance clustering threshold (height = 0.25, deepSplit = 3) was set by the cutreeDynamic function in the dynamicTreeCut package, and the minimum number of genes was 30 in each module. The automatic module merge step was performed using the mergeCloseModules function. The interested modules were the most relevant to 5-fluorouracil; therefore, the IC₅₀ value of 5-fluorouracil for module–trait relationships was analyzed.

The “limma” package (Ritchie et al., 2015) was employed to discriminate differentially expressed genes (DEGs) meeting FDR<0.05 and log2(Fold Change) > 1. DEGs were then screened by overlapping analysis with gene modules associated with 5-fluorouracil. The “ConsensusClusterPlus” package was applied to run the consensus clustering on the TCGA-LIHC data matrix (Wilkerson and Hayes, 2010). The initial step was to subsample 80% items and features. Each subsample was then partitioned into k groups. Afterward, consensus values were calculated and stored in a consensus matrix for each k-value. The output graphical plots included the consensus matrix plot and the empirical cumulative distribution function (CDF) plot.

Genomic variation includes small insertions or deletions (indels), single-nucleotide variants, and CNVs. SNVs and CNVs, all belong to the category of genomic variation. After reading the MAF file from TCGA-LIHC, the generated MAF object was passed to the “maftools” package (Mayakonda et al., 2018) for SNV analysis and oncoplot drawing. GISTIC 2.0, which calculates a statistic involving the occurrence frequency and distortion amplitude, was employed to analyze CNV data. The characteristic of this method is to identify the regions of the genome, where anomalies occur more frequently than accidentally expected and gives more weight to high-amplitude events (high-level copy number gain or homozygous deletion) that are unlikely to represent random distortions. For each important region, the method defines a “peak region” with a maximum aberration frequency and amplitude (Beroukhim et al., 2007).

The ESTIMATE algorithm, which leverages the properties of the TCGA-LIHC transcriptional profiles to infer the degrees of stromal and immune cell infiltration, was applied to determine the ESTIMATE score (Yoshihara et al., 2013). Different methods for assessing the level of immune infiltration, including CIBERSORT, ssGSEA, MCPcounter, and TIMER, were properly applied. CIBERSORT was used to measure the intra-sample (within-leukocyte) proportions of immune cell populations (Newman et al., 2015). Different from CIBERSORT, MCPcounter outputs the estimated abundance of each cell population, thereby enabling a comparison between samples to be expressed in arbitrary units (Becht et al., 2016). TIMER takes tissue specificity into account when estimating immune cell populations (Li et al., 2017), and this method helps identify associations between six types of immune cell infiltration and clustering in the TCGA-LIHC cohort.

Univariate COX regression analysis identified prognostic genes from the intersection of DEGs and 5-fluorouracil-related gene modules, and introduced five different machine learning algorithms to complete the task of variable selection, including gradient boosting machine (GBM), least absolute shrinkage and selection operator (LASSO) regression, support vector machines (SVM), Decision Trees, and Random Forest. The intersection of genes selected by each machine learning algorithm was used for stepwise regression analysis in multiple linear regression to generate a fitting regression model to evaluate the risk of samples in different HCC cohorts.

This study integrated age, gender, T stage, stage, grade, and RiskScore information, and performed univariate COX and multivariate COX analyses to determine independent prognostic factors which influenced the prognosis of HCC. Based on these independent prognostic factors, we developed a nomogram for predicting HCC survival. Based on the actual and predicted survival outcomes, we developed calibration curves to validate the predictive power of the nomogram. In addition, we also graphed decision curves to determine the prognostic guidance value of the nomogram and RiskScore.

HCC cell lines including Hep3B2.1-7 and Huh-7 were obtained from COBIOER (Nanjing, China). Hep3B2.1-7 and Huh-7 cells were cultured in DMEM F12 with 10% FBS (Gibco, Thermo Fisher, USA). Human liver epithelial cells (THLE-3) were obtained from ATCC (Manassas, VA, USA) and stored in BEGM (Lonza, Walkersville). Cells were grown at 37°C in a humidified environment containing 5% CO₂.

TRIzol reagent (Thermo Fisher, USA) was used to extract total RNA from the Hep3B2.1-7, Huh-7, and THLE-3 cell lines. Using FastStart Universal SYBR Green Master (Roche, Switzerland), quantitative reverse transcription polymerase chain reaction (qRT-PCR) was performed on the RNA extracted from each sample (2 μg) on a LightCycler 480 PCR System (Roche, USA). The cDNA was utilized as a template with a reaction volume of 20 μl (2 μl of cDNA template, 10 μl of PCR mixture, 0.5 μl of forward and reverse primers, and an appropriate water volume). The following procedures were utilized for the PCR reactions: cycling conditions started with an initial DNA denaturation phase at 95°C for 30 s, followed by 45 cycles at 94°C for 15 s, 56°C for 30 s, and 72°C for 20 s. Three separate analyses were performed on each sample. Based on the 2^−ΔΔCT method, data from the threshold cycle (CT) were obtained and standardized to the levels of GAPDH in each sample. The expression levels of mRNA were compared to controls obtained from normal tissues. Sequence lists of primer pairs for the target genes are summarized in Supplementary Table S1.

All statistical analysis and verification were conducted in the R code. A chi-squared test was adopted to detect differences in clinical characteristics between subtypes. The survival of the samples was presented by the Kaplan–Meier curve. Statistical survival difference was analyzed using the log-rank test. The time-dependent receiver operating characteristic (ROC) curve and the area under the curve (AUC) of the risk layering tool were generated and calculated using the timeROC package. p < 0.05 meant that the difference was statistically significant. For the results of the statistical analyses, ns indicated no significance, * indicated p < 0.05, ** indicated p < 0.01, *** indicated p < 0.001, and **** indicated p < 0.0001.

We initially examined the sensitivity of different HCC cell lines to 5-fluorouracil. The IC₅₀ value of 5-fluorouracil was the lowest in Hep3B2-1-7 cells and the highest in HuH-7 cells, meaning that Hep3B2-1-7 cells were the most sensitive to 5-fluorouracil and HuH-7 cells had the strongest resistance to 5-fluorouracil (Figure 1A). All cell samples from GDSC2 were clustered (Figure 1B). The soft-threshold power satisfying the scale-free topology of the network was 6, the corresponding R2 value was 0.86, and the average connectivity was very close to 0 (Figure 1C). Next, all genes were clustered into 19 interacting modules (Figure 1D). Among the 19 clustering modules, their correlations with 5-fluorouracil resistance were analyzed. The result showed that midnight blue was the module with the highest significant negative correlation with 5-fluorouracil sensitivity (Figure 1E). We analyzed GO and KEGG annotation of genes within the midnight blue module. Biological process was annotated to regulation of mRNA processing, regulation of mRNA polyadenylation, and positive regulation of telomere capping in GO terms. The protein products of these genes might be the components of the transcription elongation factor complex, DNA polymerase III complex (Figure 1F).

The differences between normal tissues and HCC tissues of log2 (TPM+1) of TCGA-LIHC were analyzed. A total of 2,356 genes with log2 (Fold Change) > 1 and FDR <0.05, as well as 462 genes with log2 (Fold Change) <-1 and FDR <0.05, were identified (Figure 2A), and 35 upregulated and 8 downregulated DEGs were also detected from the midnight blue module (Figure 2B). The samples of TCGA-LIHC were clustered according to the expression of the aforementioned 43 genes. The CDF plot helps in finding the k-value that reached the approximate maximum value of 3 (Figure 2C). The consensus matrix showed the clustering partition was k = 3 (Figure 2D). Complete separation of survival curves and overall survival (OS) of the three clusters in the detected TCGA-LIHC, HCCB18, and GSE145203 cohorts had significant differences among subtypes. Specifically, C3 had the most satisfactory survival outcome when compared with C1 and C2, while C1 had the shortest OS (Figures 2E–G). Next, the expression of 43 genes was also shown in a heatmap, which demonstrated that most genes were highly expressed in C1 than C3 and C2 (Figure 2H).

Here, in the C1, C2, and C3 subtypes, we analyzed the status of genomic variations. The TP53 mutation rate of C1 was also significantly higher than that of C2 and C3. The mutation rate of CTNNB1 in C2 was the highest, which was significantly higher than that of C1 and C3. The mutation rate of TTN in C3 was the highest, which was significantly higher than that of C1 and C2 (Figure 3A). From the Manhattan plot, we observed the CNV of each subtype at the chromosomal level, the number of high-level DNA copies amplified, and deleted in C2 was significantly less than that in C1 and C3 (Figure 3B). Comparison of the clinical characteristics of the three subtypes showed that there were more male subject samples than female subject samples in each subtype. Differences in sex, age distribution and grade, and T stage and stage characteristics were statistically significant among the three subtypes. Compared with the other two clusters, C2 had the highest proportion of male subjects and samples aged over 60 years old. C1 with the shortest OS was characterized by a high proportion of clinical grade and later T stage and stage (Figure 3C).

By running ESTIMATE, the ESTIMATE score, stromal score, and immune score of each cluster were calculated, which showed significant differences among the three clusters, and the level was the lowest in C2 (Figure 4A). Of the 22 immune cells provided by CIBERSORT, 15 showed significant differences in infiltration among three subtypes (Figure 4B). Memory cells, immunosuppressive cells (regulatory T cells (Treg) and myeloid-derived suppressor cells (MDSC)), and cytotoxic cells (CD8 T cells, natural killer (NK) cells, and NK T cells) identified in 28 TIL subpopulations showed differential infiltration among the three subtypes, and almost all of them had the least infiltration in C2 (Figure 4C). Combining the results of MCPcounter and TIMER analysis, the infiltration of CD4 T cell, T cells, B cells, macrophage, neutrophils, CD8 T cells, endothelial cells, and dendritic cells, and fibroblasts in C2 was significantly lower than that in C1 and C3 (Figures 4D, E). The differences in enrichment of pathways in clusters were archived in a heatmap, from which we could observe that the enrichment level of most pathways relevant to metabolism decreased in C1, such as linoleic acid metabolism, tyrosine metabolism, phenylalanine metabolism, and pyruvate metabolism (Figure 4F).

In the midnight blue module, a total of 43 genes were identified as hub genes, and 25 HCC prognostic genes were screened from these hub genes by univariate Cox regression analysis (Figure 5A). Machine learning models of these genes were established based on machine learning methods, including LASSO regression, GBM, SVM, Random Forest, and Decision Tree. A total of 21 genes belonged to the intersection of five machine learning models (Figure 5B). The stepwise regression method screened five genes from the 21 genes suitable for constructing a risk model, including TOMM40L, SNRPA, ILF3, CPSF6 and NUP205. The risk coefficient of each gene was obtained from multivariate Cox regression analysis, and a fitted regression model was generated: RiskScore = 0.293*TOMM40L+0.558*SNRPA-0.823*ILF3+0.493*CPSF6+0.464*NUP205. Regression models were used to calculate risk scores in the test set TCGA-LIHC and two independent verification sets HCCDB18 and GSE14520. A significant negative correlation with the sample OS in their cohorts was found, with the patients of a high-risk score showing a shorter survival time (Figures 5C–E). Meanwhile, the regression model showed stability and availability in predicting 1–5 year(s) OS of cases in TCGA-LIHC and GSE14520 (Figure 5F).

The heatmap of clinical characteristics corresponding to risk score was drawn. Differences in molecular subtypes, survival status, T stage and stage, and grade proportion were statistically significant between low-risk and high-risk groups. The high-risk group was characterized by a high C1 ratio, high mortality rate, later T stage and stage, and grade distribution, while these clinical traits in the low-risk group were significantly weaker (Figure 6A). The actual univariate Cox regression analysis showed that the risk score and T stage and stage were significantly correlated with HCC prognosis. The risk score was identified as an independent prognostic index of HCC by Multivariate Cox regression analysis (Figures 6B, C). Synthesizing information on T Stage, Stage, and RiskScore, we constructed the nomogram for assessing clinical outcomes of HCC patients at 1, 3, and 5 years (Supplementary Figure S1A). The calibration curves showed that the predicted clinical outcomes fit well with the actual observed clinical outcomes and that the nomogram had a good predictive value (Supplementary Figure S1B). In addition, the decision curve also showed that there is an excellent applicability of the nomogram and RiskScore in assessing clinical outcomes in HCC (Supplementary Figure S1C). The risk score of the T3–T4 stage, stage Ⅲ–Ⅳ, and G3–G4 samples was significantly higher than that of the T1–T2 stage, stage Ⅰ–Ⅱ, and G1–G2 samples, respectively (Figures 6D–F).

The degree of immune cell infiltration was evaluated according to the risk score. We detected that the cells with the most different degrees of infiltration in the high-risk and low-risk groups were type 2 T helper cells, central memory CD4 T cells, type 1 T helper cells, plasmacytoid dendritic cells, effector memory CD4 T cells, activated CD4 T cells, activated CD8 T cells, eosinophils, natural killer T cells, CD56 dim natural killer cells, activated dendritic cells, and effector memory CD8 T cells (Figure 7A). The ssGSEA score correlation analysis of the risk score and immune cells showed that the degree of correlation between the 12 cells had infiltration and risk score differences between the high-risk and low-risk groups. The correlation between risk score and CD56 dim natural killer cells was almost negligible, and the other 11 cells showed a significant correlation with risk score. Among them, effector memory CD4 T cells, activated CD4 T cells, natural killer T cells, plasmacytoid dendritic cells, type 2 T helper cells, central memory CD4 T cells, and activated dendritic cells showed a significant positive correlation with the risk score, while activated CD8 T cells, effector memory CD8 T cells and type 1 T helper cells, and eosinophils were significantly negatively correlated with the risk sore (Figure 7B).

To verify the reliability of the RiskScore, we detected the expression of the five genes by PCR. The results of PCR corroborated the reliability of the RiskScore, and we found that TOMM40L, SNRPA, ILF3, CPSF6, and NUP205 were significantly upregulated in the HCC cell lines Hep3B2.1-7 and Huh-7 compared to human normal liver epithelial cells THLE-3 (Figures 8A–E).