Gene expression data and corresponding prognostic and clinicopathological data of HCC were downloaded from the International Cancer Genome Consortium (ICGC) (24), The Cancer Genome Atlas (TCGA) (25), and Gene Expression Omnibus (26). We obtained GSE104580 cohorts (100 TACE response tissues and 100 TACE nonresponse tissues) and TCGA Liver hepatocellular carcinoma (LIHC) cohorts (374 HCC samples and 50 control samples) for subsequent analyses. The “normalizeBetweenArrays” function in “limma” R package was performed to background adjustment and quantile normalization (27). Besides, a total of 500 HCC samples in GSE14520 and ICGC cohorts, containing gene expression array and prognostic and clinical information, were collected as validation sets. The genes with lower expression and samples with no prognostic information were excluded.

In TCGA and GSE104580, differentially expressed genes (DEGs) between different subgroups were respectively identified utilizing the “limma” package with a fold change of 1.5 and an adjusted p < 0.05. Heatmaps for expression of top 50 up-/down-regulated DEGs were drawn by “pheatmap” R package. Venn diagram of the common DEGs between two datasets was visualized by “venndiag” package. We regard these common DEGs as TRGs and further explored their potential functions and enriched pathways using the “clusterprofiler” package in R and p < 0.05 was considered as significant (28).

All HCC patients in TCGA were randomly categorized into training (n = 223) and testing (n = 151) sets according to a ratio of 6:4. The TRGs were subjected to univariate Cox regression analysis for identifying those related to OS. Next, the LASSO regression algorithm was conducted to obtain the optimal prognostic genes with 10-fold cross-validation using the “glmnet” R package (29). LASSO regression can improve the accuracy and interpretability of the model and also reduce the risk of overfitting between independent variables (30). Based on these candidate genes, a TACE refractoriness–related prognostic score (TRP score) was established using the multivariate Cox regression. The formula is as follows: TRP score = Σ (Expi * coefi), where Expi represents gene expression level, and coefi denotes risk coefficient. Receiver operating characteristic (ROC) curve analysis was used to illustrate the diagnostic veracity of TRP score in different HCC cohorts (TCGA, ICGC, and GSE14520). A nomogram including candidate genes was constructed using the “survival” and “rms” packages, and its prediction accuracy was evaluated by calibration curves. Patients in the training set were divided into the high-risk group and the low-risk group based on the median TRP score and subsequently carried out to Kaplan–Meier survival analysis. Principal component analysis (PCA) was performed to assess the differences in the distribution of the low- and high-risk groups, and ROC curves were applied to verify the exactness of the score model in TCGA. Simultaneously, to further verify the accuracy of the score model, the testing and external validation sets (ICGC and GSE14520) were also subjected to survival analysis and ROC curve analysis.

Samples in GSE104580 were randomly divided into training (n = 120) and testing (n = 80) sets. In the training set, identical candidate genes were performed to establish TACE refractoriness diagnostic score (TRD score) using a logical regression algorithm. The relevance between TRP score and TRD score was assessed by Pearson approach. A diagnostic nomogram was constructed to predict the probability of occurrence of TACE refractoriness in HCC patients. Diagnostic ROC curves and calibration plots analyses on account of training sets, test sets, and whole datasets were used to evaluate the efficiency and accuracy of the diagnostic model.

To further understand the clinical relevance of the TRP score, we explore the impact of the TRP score on other clinical characteristics (age, sex, Child grade, histologic grade, TNM stage, cancer type, and radiation treatment). The difference in TRP score between clinical subgroups was compared by independent-samples t test procedure. In addition, the univariate Cox regression analysis was used to evaluate the correlation between prognosis and clinical characteristics, and multivariate Cox regression analysis was applied to analyze the independent prognostic ability of the risk factors.

Normalized enrichment score of pathways and functional annotation was calculated by the gene set variation analysis (GSVA) method using the “GSVA” R package, and p < 0.05 and false discovery rate (FDR) < 0.25 were considered to be statistically significant (31). The pathways significantly distributed between the high- and low-risk groups were selected by t tests according to the criteria of p < 0.05. The DEGs between two risk stratifications were identified by the “limma” R package (p < 0.05, logFC >1.5), and the Gene Ontology (GO) analysis was used to identify the characteristic biological attributes in DEGs. The difference in immune response between two risk subtypes was estimated by gene set enrichment analysis (GSEA) algorithm (32).

The ESTIMATE score was calculated using R package “ESTIMATE” to evaluate the composition of the stromal cells and the immune cells (33). CIBERSORT and IMMUNCELL AI algorithms were used to quantify infiltration levels for different immune cell types between high- and low-risk groups using gene expression profiles from the TCGA database (34, 35). The distinction in immune infiltrating level and ESTIMATE score between high- and low-risk groups were analyzed by the Wilcoxon test. To further validate the immune phenotype of different prognostic subtypes, we investigated the distribution of immune molecules, which negatively regulated the antitumor immune response between high- and low-score groups. These immunosuppressive genes were downloaded from the Tracking Tumor Immunophenotype website (http://biocc.hrbmu.edu.cn/) (36, 37).

In order to explore the somatic mutations between high- and low-risk groups, the “maftools” R package (37) was applied to analyze the mutation annotation format from TCGA. The somatic mutation data were applied to calculate tumor mutation burden (TMB), an effective immunotherapy predictor (38). Meanwhile, other accepted biomarkers of immunotherapies, such as mismatch repair deficiency (39, 40), HLA genes (41), and immune checkpoint genes (ICPs) (42), were used to reflect the difference in immunotherapeutic sensitivity between high- and low-risk groups. In addition, a series of predicted scores to the checkpoint immunotherapy [cytolytic activity (CYT) score, Tumor Immune Dysfunction and Exclusion (TIDE) score, TME score, and immunophenoscore (IPS)] were conducted to forecast the response to ICP in tumors. CYT score was defined as the geometric mean of PRF1 and GZMA expression (43). Then, we uploaded the HCC transcriptome profiles to the TIDE website (http://tide.dfci.harvard.edu/) and then obtained corresponding TIDE scores (44). TME score, a novel approach to evaluate the efficacy of immune checkpoint blockade, was computed by “TMEscore” R packages (45). The patient’s IPS was obtained from The Cancer Immune Group Atlas (https://tcia.at/home), and a higher IPS is positively correlated to the increased immunogenicity (46). All the aforementioned scores were carried out based on gene expression data from TCGA.

The response to chemotherapy drugs of HCC patients was predicted based on the public pharmacogenomics database Genomics of Drug Sensitivity in Cancer (47). To compare the therapeutic effects of chemotherapeutic drugs in the two different prognostic groups, we measured the semi-inhibitory concentration (IC₅₀) values of commonly used chemotherapeutic drugs for LIHC by the “pRRophetic” package (48).

In accordance with the outcome of the independent prognosis analysis, we developed a predictive nomogram by integrating diverse clinicopathological information, including TNM stage, cancer type, and TRP score, to comprehensively estimate patient survival. Time-dependent concordance index (C-index) and time-dependent ROC analysis was conducted to compare the accuracy of nomogram and different variables by R package “pec” and “riskRegression.” Calibration plots illustrated the consistency between the predicted 1-, 3-, and 5-year endpoint events and the authentic outcomes. Simultaneously, the decision curve analysis of 1, 3, and 5 years was performed to evaluate the clinical utility of the nomogram. Moreover, the external validation of the combined model was carried out based on the ICGC dataset.

All statistical analyses were conducted via R software (version 3.6.7). The Student t test was used for statistical comparisons. Spearman correlation was applied for the analysis of the correlation. The Benjamini–Hochberg FDR method was used for p value adjustment. Fisher test was used to identify the significant GO terms. p < 0.05 was regarded as statistically significant. The cutoff values of continuous variables, such as gene expression and immune infiltration level, were the median.

The general analysis flow of our study is shown in Figure 1 . To explore potential TRGs in HCC, we identified 1,183 DEGs from TCGA set and 347 DEGs from GSE104580 set using R package “limma” ( Figure 2A ). Expression patterns of the top 50 up-/down-regulated DEGs of the two datasets are shown in Figure 2B . By the intersection of DEGs in two datasets, a total of 487 common DEGs were screened out and named as TRGs ( Figure 2C ). According to the GO annotations, the TRGs were mainly associated with the regulation of cell cycle, chromosome segregation, mitotic nuclear division, regulation of inflammatory response and immune effector process, response to drugs, and organelle fission ( Figure 2D ).

Univariate Cox regression analysis found 25 TRGs significantly associated with OS in the training set. Subsequently, four TRGs (TTK, EPO, SLC7A11, and PON1) were further confirmed as the optimal combination parameters of prognosis prediction according to LASSO–Cox regression ( Figures 3A – C ). Ultimately, the TRP score was calculated according to the following equation: TRP score = (0.264 ∗ TTK) + (0.147 ∗ EPO) + (0.226 ∗ SLC7A11) + (−0.091 ∗ PON1). We constructed a nomogram incorporating the four markers for calculating TRP score easily and predicting the 1-, 3-, and 5-year OS of tumors roughly ( Figure 3D ). The C-index of TRP score ranked first among all the variables ( Figure 3E ). The calibration curves revealed that the 1-, 3-, and 5-year OS predicted by the nomogram were close to the ideal performance ( Figure 3F ). ROC analysis proved that TRP score exhibited a great diagnostic accuracy for HCC in diverse independent verification datasets, among which area under the ROC curve (AUC) = 0.785 in GSE104580, AUC = 0.825 in TCGA, AUC = 0.825 in TCGA, AUC = 0.79 in GSE14520, and AUC = 0.928 in ICGC ( Figure 3G ). The PCA results demonstrated that patients in the high- and low-risk groups were clearly distributed on both sides ( Figure 3H ).

These prognostic signatures were also involved in logical regression to evaluate their diagnostic property on TACE refractoriness based on the training set. The formula of diagnostic score is as follows: TACE refractoriness diagnosis score (TRD score) = (0.315 ∗ TTK) + (0.411 ∗ EPO) + (0.151 ∗ SLC7A11) + (−0.262 ∗ PON1). A diagnostic nomogram for predicting the likelihood of HCC patients suffering TACE refractoriness is shown in Figure 4A . There is a strong correlation between the TRP score and TRD score (p = 2.2e-16, cor=0.8) ( Figure 4B ), which implies that TACE refractoriness may affect the HCC patients’ outcomes. Diagnostic ROC curves demonstrate that TRD score had excellent prediction performance in training, test, and whole sets ( Figure 4C ). The calibration curves of the diagnostic model presented that predicted values approach actual performance ( Figure 4D ).

Based on the training set, whole TCGA set, ICGC set, and GSE14520 set, patients were divided into high- and low-risk groups by the median TRP score. The distribution of the risk score and survival status in these datasets was visualized in Figure 5A and S1A . The expression patterns of four prognostic TRGs were obviously different in the two risk stratifications ( Figure 5A and S1A ). Kaplan–Meier survival curves showed a longer OS in patients with low-risk scores than those with high-risk scores ( Figure 5B and S1B ). Furthermore, time-dependent ROC curves of risk score displayed high prediction accuracy for OS at 1, 3, and 5 years in several datasets ( Figure 5C and S1C ). These results indicate the reliable ability of the TRP score in distinguishing differences in tumors’ outcomes.

To evaluate the independent prognostic ability of TRP score, we combined a variety of clinical features with the TRP score to perform univariate and multivariate COX regression in TCGA cohorts. As shown in Figures 6A, B , TRP score, TNM stage, Child grade, and cancer type showed significant differences, which meant that these indicators could be independent prognostic predictors. Subsequently, we explored the clinical relevance between TRP score and other clinical parameters. The results suggest that TRP score was significantly related to histology grade (p < 0.001), Child grade (p < 0.001), T stage (p < 0.001), and TNM stage (p < 0.001), and it gradually elevated with the development of these clinical factors ( Figures 6C, D ).

Through GSVA algorithm to estimate the variation of pathway activity in the TCGA sample population, we found that the high-risk group was significantly enriched for the cell cycle, such as homologous recombination, mitotic spindle and G2/M checkpoint, and those related to the DNA repair, such as mismatch repair, whereas the activity of pathways related to metabolism, such as fatty acid metabolism, glycine, serine, and threonine metabolism, enriched significantly in the low-score group ( Figures 7A, B ). The up-regulated DEGs in the high-risk group significantly enriched in those terms related to mitotic nuclear division, chromosome segregation, mitotic spindle organization, and nuclear division, which is in keeping with the result of GSVA and implied that the molecular signature of the high-score group was closely associated with the cell cycle ( Figures 7C, D ). It is worth noting that the immune responses were more active in the low-risk group rather than the high-risk group, as revealed by GSEA ( Figure 7E ).

The infiltration of immune cells in HCC was estimated by three independent methods (CIBERSORT, IMMUNECELL AI, and ESTIMATE) to evaluate the interaction of prognostic score with tumor immune microenvironment. The immune score, stromal score, and ESTIMATE score of the high-risk group were lower than those in the low-risk group, whereas tumor purity was higher in the high-risk group ( Figure 8A ). After examining the results of CIBERSORT and IMMUNECELL AI, we found that regulatory T cells (Tregs) and macrophages were both significantly elevated in high-risk group ( Figures 8B, C ). In addition, based on 22 immune cell subpopulations from CIBERSORT, patients with high- and low-risk could be clearly divided into two discrete groups by PCA analysis ( Figure 8D ). As shown in Figure 8E , the immunosuppressive genes were essentially up-regulated in high-risk group. On the other hand, with the help of CIBERSORT and IMMUNECELL AI algorithms, we found that macrophages increased in the TACE nonresponse group ( Figures S2A, B ). This evidence indicates that patients with high TRP score and TACE refractoriness had low activities of antitumor immune processes.

Through the mutation landscape of HCC patients, we discovered a higher TP53 mutation rate existed in the high-score group rate than the low-score group (45% vs. 16%, Figures 9A, B ). The TMB slightly increased in the high-score group ( Figures 9C ), and high TMB suggested a bad prognosis ( Figures 9D ). Meanwhile, preferential expression of MMR-related genes (MLH1, MSH2, MSH6, and PMS2) and HLA genes were found in the low-risk group ( Figures 9E, F ). Six immune checkpoint molecules (CTLA4, PD1, PD-L1, LAG3, HAVCR2, and TIGIT) were significantly up-regulated in the high-risk group ( Figures 10A, B ). To evaluate the impacts of TRP score on the tumors’ response to ICP therapy, we analyzed the correlation between TRP score and multiple immunotherapeutic indices. The results demonstrated that the TRP score was actively associated with the CYT score and TME score, while negatively associated with the TIDE score ( Figures 10C – E ). Furthermore, IPS analysis indicated the IPS, IPS-CTLA4, IPS-PD1, and IPS-PD1-CTLA4 scores were higher in the high-score group ( Figure 10F ). As for TACE refractoriness, higher ICP (PD1, HAVCR2, and TIGIT) expression and lower TIDE score mainly concentrated in the TACE nonresponse group ( Figure S2C ). Thus, these results revealed that TRP score may serve as targets for immunotherapy, and patients with high-TRP and TACE refractoriness seem more sensitive to ICPs.

Based on TCGA cohorts, we calculated the IC₅₀ of common chemotherapeutic drugs recommended for liver cancer treatment by the American Joint Committee on Cancer guidelines to evaluate the sensitivities of patients in the low- and high-risk groups to these drugs. We discover that the patients in the high TRP score group had lower IC₅₀ value for sorafenib, doxorubicin, vinblastine, and cisplatin ( Figure 10G ). Meanwhile, patients with TACE refractoriness showed lower IC₅₀ value for sorafenib ( Figure S2D ). These results indicated that TRP score and TACE refractoriness can significantly impact patients’ sensitivity to some drugs.

For synergizing clinical efficacy of PRG score in the prognosis, we established a nomogram containing the PRG score and other clinicopathological characteristics to predict the 1-, 3-, and 5-year OS with better precision ( Figure 11A ). Time-dependent C-index analysis demonstrated that the nomogram exhibited superior predictive performance about survival rate compared with other clinical parameters ( Figure 11B ). The calibration curves of the nomogram of 1-, 3-, and 5-year OS prediction were approached to the ideal model ( Figure 11C ). Meanwhile, time-dependent ROC curves at 1-, 3-, and 5-year OS indicated that the nomogram possessed the best predicting precision compared with TRP score and TNM stage ( Figure 11D ). DCA at 1-, 3-, and 5-year OS proved the nomogram had great clinical utility ( Figure 11E ). As for validation analyses in the external dataset (ICGC-LIHC), the nomogram showed the same prognostic value through comprehensive evaluation with time-dependent C-index, calibration curve, and decision curve ( Figure S3 ).