Level 3 mRNA expression, somatic mutation data, and clinicopathological data were obtained from the TCGA website (https://portal.gdc.cancer.gov/repository). A segment of copy number for the TCGA-LIHC cohort was accessed from the GDAC FireBrowse (http://firebrowse.org/). The raw count data were transferred to transcripts per kilobase of exon model per million mapped reads (TPM) data which would represent the expression of mRNA in the TCGA-LIHC cohort. After filtering mRNAs with low median absolute deviation (mad ≤0.5) across all samples and removing the samples without complete survival information, a total of 351 HCC samples were enrolled in this study. RNA-seq data, somatic mutation data, and clinical data with 240 tumor samples were downloaded from the International Cancer Genome Consortium (ICGC) portal (https://dcc.icgc.org/projects/LIRI-JP). Raw read count values were transformed into TPM values for subsequent analysis. The expression data and detailed clinical information of GSE14520 (including 219 HCC samples based on the GPL3921 platform) were downloaded from the Gene Expression Omnibus (GEO) (http://www.ncbi.nlm.nih.gov/geo/). Additionally, a LICH-CN cohort with 159 Chinese HCC patients was downloaded for somatic mutation, clinical outcome, and transcriptome expression FPKM value from the literature (Gao et al., 2019).

DDR-gene list including 557 genes was assembled from relevant gene lists, including MSigDB from the Broad Institute (http://www.broad.mit.edu/gsea/msigdb/) or literature (Pearl et al., 2015; Knijnenburg et al., 2018; Chen J. et al., 2021). Finally, 384 DDR genes detected in all datasets were analyzed in this study (Supplementary Table S1). Then, each gene pair was calculated via their gene expression level in each HCC sample. According to the pairwise comparison, the calculated score was 0 when the first expression level of the DDR gene was higher than that of the following DDR gene; otherwise, the calculated score was 1. DRGP scoring 0 or 1 in more than 90% of the samples were removed because they could not provide discriminative patients with different survival. The remaining DRGPs were considered as initial candidate DRGPs.

Patients were randomly divided into training and testing sets at cut-of 7:3 in the TCGC cohort. Then, a univariate Cox regression analysis was performed to identify the significant DRGPs related to overall survival (OS) if the FDR p-value was less than 0.05. Next, candidate DRGPs were submitted to LASSO-penalized Cox regression analysis based on package “glmnet” in R to construct an optimal prognostic signature in TCGA training datasets (Friedman et al., 2010). A DDR-related gene pair riskscore of each sample was calculated based on the lasso Cox regression model coefficients (β) multiplied with its DRGP score, as follows: R i s k s c o r e = ∑ i = 1 n ( β i × S c o r e i ) where Score _i is the relative expression of DRGPs for patient j in each cohort and β _i is the LASSO Cox coefficient of the DRGPs _i . Then, all patients were separated into low- (LRisk) or high-risk (HRisk) groups at the median cut-off. Kaplan–Meier survival curves were plotted for prediction of the clinical outcomes in the two groups via the “survival” package in R. The differences in survival were evaluated via the log-rank test. Time-dependent receiver operating characteristic (ROC) analysis curves were built, and the area under the curves (AUCs) for 1-, 3-, and 5-year overall survival (OS) were calculated utilizing the “survivalROC” package in R (Heagerty et al., 2000). The same method was further investigated in the TCGA testing cohort, TCGA whole cohort, GSE14520 cohort, and ICGC-LIRI cohort.

To explore the diagnostic capability of the DRGP prognostic signature in different levels of other clinical prognostic parameters, HCC samples in TCGA sets were stratified into different subgroups based on age (≥60 and <60), gender (female and male), TNM stage (I and II + III + IV), grade (G1+G2 and G3+G4), and TP53 (mutation and wild). Then, cancer samples in each subgroup were clustered into HRisk and LRisk groups. The differences in prognosis between the two groups were assessed via Kaplan–Meier OS analysis, followed by a log-rank test.

To elucidate whether the prognostic model for OS is independent of other prognostic factors, we presented univariate Cox regression analysis and multivariate Cox regression survival analysis to predict the clinical outcomes of HCC patients, which was visualized via package “forestplot” in R. Hazard ratio (HR), 95% confidence interval (CI), and p-value were calculated, respectively.

A nomogram was constructed based on all independent prognostic parameters screened by univariate and multivariate Cox proportional hazards regression analysis to predict the probability of 1-, 3-, and 5-year OS using the “rms” package of R software. Then, we used a calibration curve to visualize the performance of the nomogram with the observed rates of the TCGA whole set at corresponding time points by a bootstrap method with 1000 resamples. Furthermore, decision curve analysis (DCA) and calibration curves were detected to check the reliability of our nomogram (Kerr et al., 2016).

To investigate the difference in biological process between HRisk and LRisk, we performed some enrichment analysis using “GSVA” and “clusterprofiler” R packages (Yu et al., 2012; Hänzelmann et al., 2013). The infiltrating score of 24 microenvironment cell types was calculated with single-sample gene set enrichment analysis (ssGSEA) in the “GSVA” R package. Gene set enrichment analysis (GSEA) was conducted between HRisk and LRisk by using the R “clusterProfiler” package. A signature of eleven oncogenic pathways and a DDR gene list, which include eight core DDR pathways, were obtained from the literature (Pearl et al., 2015; Sanchez-Vega et al., 2018; Lu et al., 2021). Then, we used the gene set variation analysis (GSVA) method to generate enrichment scores for each cohort using the R package “GSVA”. The KEGG gene sets (c2.cp.kegg.v7.4.symbols.gmt) was selected as the reference datasets, which was obtained from the MSigDB database.

For immunotherapy, the tumor immune dysfunction and exclusion (TIDE) algorithm (http://tide.dfci.harvard.edu/) was applied to predict potential clinical response to immune checkpoint inhibitors (Jiang et al., 2018). Based on Genomics of Drug Sensitivity 2016 (GDSC 2016; https://www.cancerrxgene.org/), the R package “pRRophetic” was applied to estimate the chemotherapeutic sensitivity by the half-maximal inhibitory concentration (IC₅₀) of each HCC sample in four cohorts. Therefore, we could investigate the different sensitivity of common liver cancer chemotherapy drugs between HRisk and LRisk. In addition, in order to identify potential drugs in HCC samples, we performed a two-step analysis to find candidate compounds as described previously (Yang et al., 2021). First, differential drug response analysis between top decile riskscore samples and bottom decile riskscore samples was conducted to verify drugs with significantly different estimated IC₅₀ in two riskgroups (|log2FC| > 0.2). Next, Spearman correlation analysis was utilized to calculate the correlation coefficients between riskscore and IC₅₀ of each candidate drug (|Spearman correlation coefficient| > 0.4).

HBV oncoproteins were quantified for expression according to the previous study: HBVgp2_S, HBVgp3_X, HBVgp4_c, and HBVgp2_pre-S1/S2 (Xue et al., 2021). The four HBV oncoprotein expressions were identified and presented as FPKM values. To comprehensively explain the original expression level of HBV oncoproteins, we established a variable that was calculated by principal component analysis (PCA) as the previous study described (Lu et al., 2019). HBV _pca was derived from the first and second principal components that represented 76.73 and 17.10% of the variation in the original data, respectively. The coefficients of four HBV oncoproteins to the first and second principal components are shown in Table 1.

Mathematically, let E _ij denotes the log₂(FPKM +1) value of specific oncoprotein j in sample i, and C _jk represents the corresponding coefficient of HBV oncoprotein (HBV _j ; j∈{1,2,3,4}) for principal component k (k∈{1,2}). The HBV _pca can be calculated as follows: H B V p c a =   [ E 11 ⋯ E 1 j ⋮ ⋱ ⋮ E i 1 ⋯ E i j ] [ C 11 ⋯ E 1 k ⋮ ⋱ ⋮ C j 1 ⋯ E j k ] [ 0.7673 0.1710 ]

We then investigated the genomic variation between HRisk and LRisk groups. The mutation landscape was analyzed by the R package “maftools” with the initial removal of 100 FLAGS (Mayakonda et al., 2018). The data CNV segments were detected by Genomic Identification of Significant Targets In Cancer 2.0 (GISTIC 2.0) analysis. In the process of GISTIC 2.0 analysis, except for the refgene file which was “Human_Hg19.mat”, parameters were set to the default parameters. The individual fraction of genome altered (FGA), fraction of genome lost (FGL) and fraction of genome gained (FGG) for the HCC in the TCGA cohort were calculated as the study described (Lu et al., 2021). We also obtained GISTIC calls comprising −2 (deletion), −1 (loss), 0 (diploid), 1 (gain), and 2 (amplification) from GISTIC2.0 (Wu et al., 2020).

All statistical analyses were performed with R software (version 4.1.1: http://www.r-project.org) and R Bioconductor packages in this study. Kaplan–Meier analysis with the log-rank test was used to detect differences of OS between different groups through the package “survminer” in R. Time-dependent ROC was utilized to evaluate the predictive accuracy of the DRGP riskscore through package “survivalROC” in R. Cox proportional hazards regression for estimating the hazard ratios (HRs) and 95% confidence interval (CI). Comparison of a continuous variable in two groups was performed using Wilcoxon rank-sum test. Correlation between two continuous variables was measured by Spearman’s rank-order correlation. Differences in proportions were compared by the Chi-squared test or Fisher’s exact test.

The clinical features of HCC samples in the training and validation sets are listed in Supplementary Table S2. In the training datasets of 246 patients, 85 patients died during the follow-up. As tested by the univariate Cox regression OS analysis, 459 DRGPs had significant associations with OS of HCC in the training set (all FDR p-value < 0.05). Based on LASSO-Cox regression analysis, 23 DRGPs were independently related to the prognosis of HCC (Figures 1A,B). Among them, 11 gene pairs were risk factors for HCC prognosis (HR > 1). Supplementary Table S3 lists the 23 selected gene pairs and their coefficients. The regression coefficients and DRGP score of these 23 gene pairs in each sample were used to calculate the riskscore in each HCC cohort.

The patients in the training datasets were divided into HRisk or LRisk groups at median cut-off. Kaplan-Meier survival analysis indicated that the HRisk group have a poorer OS than the LRisk group (hazard ratio (HR) = 3.99, 95% CI = 2.48–6.44, p-value < 0.001, Figure 1C). The values of AUC are 0.76, 0.83, and 0.86 at 1-, 3-, and 5-year follow-up, respectively (Figure 1I), showing that the signature displays good sensitivity and specificity. The C-index of the DRGP model is 0.80. To determine the predictive ability of this prognostic model, we calculated individual riskscore with the aforementioned method and classified the patients in TCGA-testing set, TCGA whole set and other validation sets into HRisk and LRisk groups. Similarly, we validated the prediction of signature in these datasets. Consistent with the above findings, the HRisk patients in all cohorts have a markedly shorter OS than those in the LRisk group (Figures 1D–H). The AUCs of ROC curves for 1-,3-, and 5- year OS are shown in Figures 1J–N.

To further explore the clinical potentiality of the prognosis model in HCC, stratified analysis based on these clinical characteristics was conducted. As shown in Figures 2A–J and Supplementary Figure S1, Kaplan–Meier OS curves also showed that HRisk patients had considerably worse OS than LRisk patients, which further indicated the excellent prediction of the DRGP model. We further analyzed whether the riskscore was an independent prognostic predictor for OS. In univariate Cox regression analyses, high riskscore was significantly associated with shorter OS in TCGA cohort (HR = 3.22, 95% CI = 2.53–4.13, p-value < 0.001, Figure 3A) According to the multivariate Cox regression analysis results, we considered the TNM stage (HR = 1.27, 95% CI = 1.00–1.61, p-value = 0.0515) and riskscore (HR = 2.96, 95% CI = 2.29–3.83, p-value < 0.001] are both independent prognosis factors for TCGA (Figure 3A). The independence of the DRGP signature for HCC prognosis was also confirmed in GSE14520, ICGC, and LIHC-CN cohorts (Figures 3A,B). Collectively, the signature was an independent prognostic factor for HCC.

We constructed a nomogram based on multivariate Cox regression analysis for prediction of the 1-, 3-, and 5-year survival probability in TCGA datasets (Figure 3C). As shown in the calibration chart (Figure 3D), the nomogram could robustly predict OS for HCC patients. Moreover, the DCA curve suggested that riskscore was more beneficial when compared with the TNM stage alone (Figure 3E). The DCA curve demonstrated that the net benefit of the combined model was comparable to the riskscore. These results showed that the nomogram built with the combined model might help clinical management.

A previous study indicated that the high expression of HBV16 E6/E7 was significantly linked to a favorable prognosis because the inflammatory/immune response of the host may be stimulated (Lu et al., 2019). To investigate whether HBV oncoproteins were differentially expressed between two risk groups in HCC, we calculated the HBV _pca by principal component analysis (PCA) based on the FPKM value of four HBV oncoproteins in the TCGA cohort. In our study, 96 HVB-infected HCC patients were identified, the first and second principal components were used since they covered almost the variations (93.83%; Figure 4A). We found a significant difference in HBV virus expression between HRisk and LRisk (Wilcoxon test, p = 0.018, Figure 4B), and HBV _pca has also observed a mildly negative correlation with riskscore (R = −0.280, p = 0.006, Figure 4C). It means a high level of HBV _pca corresponds to low risk in HCC patients.

To evaluate the tumor immune microenvironment in different groups, a ssGSEA method was used to estimate the infiltration levels of the 24 types of immune cells. The ssGSEA score and immune cell types which were differentially infiltrated between LRisk and HRisk groups in the TCGA set are presented in Figure 4D. The proportion of 24 immune cells in each group is shown in a bar plot. The results revealed that the level of infiltration of dendritic.cells.resting, macrophages.M0, mast.cells.activated, and T.cells.CD4.memory.activated in the HRisk group was significantly higher than that in the LRisk group, while the level of endothelial cells, mast.cells.resting, and NK cells.activated in the LRisk were higher than that in the HRisk group. For further investigating the immune landscape of different risk groups reflected by the DRGP signature, validation cohorts GSE14520, ICGC, and LIHC-CN were also calculated by ssGSEA to verify the differences in risk groups at the immune level (Supplementary Figure S2).

We also detected and compared the expression levels of several immune checkpoints between the LRisk and HRisk groups. Results showed that the mRNA expression levels of CTLA4, PDCD1, and TIGIT were consistently overexpressed in the HRisk in TCGA and ICGC datasets (Figure 5A), but we could not observe a significant difference in GSE14520 and LIHC-CN cohorts (both, p > 0.05, not shown). These results suggested that the HRisk group may contribute to tumor immune dysfunction and immune exclusion in HCC. To investigate whether LRisk responds to immune checkpoint inhibits, we harnessed the TIDE algorithm to predict the potential response to immunotherapy in different groups. The higher TIDE score represented less promising treatment for response to immunotherapy. In our results, the HRisk group contained lower proportion of TIDE-predicted responders than LRisk group in all four cohorts (TCGA (p < 0.001), GSE14520 (p = 0.23), ICGC (p = 0.012), and LIHC-CN (p = 0.016); Figure 5B). These results suggest that HCC patients in LRisk might be more beneficial from immune checkpoint inhibitors.

Considering that chemotherapy is a common way to treat liver cancer, we first verify whether the DDR patterns groups may affect the sensitivity of chemotherapeutic drugs commonly used for treating liver cancer (including cisplatin, 5-fluorouracil, gemcitabine, oxaliplatin, doxorubicin, mitoxantrone, gefitinib, and sorafenib). We found that HRisk groups of all four HCC cohorts were highly sensitive to 5-fluorouracil (all, p < 0.05; Figure 5C, Supplementary Figure S3), and four LRisk groups presented with a significantly higher response to gefitinib and gemcitabine (all, p < 0.05; Figure 5C, Supplementary Figure S3). Next, we performed a two-step analysis to find potential therapeutic compounds. Eventually, the analysis obtained seven compounds (including CDK9_5576, CDK9_5038, bleomycin, midostaurin, SNX-2112, BMS-754807, and podophyllotoxin bromide) that had lower IC₅₀ in HRisk and a negative correlation with riskscore (Figure 5D, Supplementary Figure S4) across four datasets and five compounds (including trichostatin A, gefitinib, afatinib, selumetinib, and EphB4_9721) were observed to present a significant response to LRisk and a positive correlation with riskscore (Figure 5E, Supplementary Figure S4).

To better characterize the two HCC riskgroups, differential analyses were performed. Gene set enrichment analysis (GSEA) was conducted using the “clusterProfiler” package, and enrichment differences of pathways were significant if the FDR p-value < 0.15 and |NES| >1 in all four cohorts. The results indicated that 28 metabolism-relevant pathways were significantly upregulated in LRisk, while HRisk enriched in cell cycle, DNA replication, spliceosome, and DDR-relevant pathways (Figure 6A). Thus, the HRisk group presents upregulated cell cycle procession and DDR pathways which might contribute to the hyperproliferation and development of tumor cells. Pathway with significant differences in enrichment in all four cohorts was considered subclass specific pathway. GSVA was conducted to quantify and visualize the enrichment of 28 metabolism-related pathways which were classified into four specific metabolism signatures, including lipid metabolism relevant pathway, drug metabolism relevant pathway, carbohydrate metabolism relevant pathway, and amino acids metabolism relevant pathway (Figure 6A). Results confirmed that the LRisk group has significant upregulation of metabolism signatures, consistent with the results from GSEA.

To further investigate the activation of oncogenic pathways among HRisk and LRisk (Figure 6C), We found that the cell cycle oncogenic pathway was significantly activated in HRisk, while LRisk had a higher score of angiogenesis and Wnt activation-relevant pathways than HRisk (Figure 6D). Considering that the risk groups were divided based on DDR-relevant genes signature, we then decided to further explore whether different characteristics exist in distinct DDR pathways. Eight DDR core pathways were quantified using the GSVA algorithm. HRisk group exhibited higher expression for all eight pathways than LRisk. The same results were demonstrated in cohorts GSE14520, ICGC, and LIHC-CN (Supplementary Figures S5−S7).

We finally investigated the genomic variations between two different risk groups in the TCGA cohort. To analyze whether differences exist in the somatic variations (10%) of HCC between two riskgroups, R package “maftools” was used (Figure 7A). The results showed that the HRisk group had significantly high TP53 mutations and less CTNNB1 mutations than the LRisk group (TP53: 46.3 vs. 12.5%, p < 0.001; CTNNB1: 17.7 vs. 31.8%, p = 0.007; Figure 7B), and we obtained consistent results in ICGC cohort (TP53: 51.7 vs. 18.8%, p < 0.001; CTNNB1:35.0 vs. 41.9%, p = 0.2883; Figure 7B) and LIHC-CN cohort (TP53: 73.4 vs. 43.8%, p < 0.001; CTNNB1: 17.7 vs. 21.2%, p = 0.5439; Figure 7B), respectively. Mutation in TP53 is the most common genetic change in HCC, and patients with mutated TP53 have shorter OS than those with wild-type TP53. Therefore, we conducted special subgroup analyses stratifying samples according to the combination of TP53 mutation status and riskgroups. We found that some patients with mutant type TP53 in the HRisk group had significantly shorter OS than those with mutant wild TP53 in the HRisk group (p = 0.89 in TCGA, p = 0.58 in ICGC, p = 0.34 in LIHC-CN; Figure 7C), while patients with mutant type in the LRisk group had longer OS than those with wild type in the HRisk group (p = 0.0028 in TCGA; Figure 7C), but no significance could be calculated in ICGC and LIHC-CN datasets (p = 0.58 in ICGC, p = 0.40 in LIHC-CN; Figure 7C). These results confirm again that the prognosis model is robust and superior. Unfortunately, we could not find the mutation data of GSE14520.

Copy number variations were a common form of genomic structural change, and an amount of research has demonstrated chromosomal abnormalities play key roles in HCC. GISTIC 2.0 was used to analyze the copy number of HRisk and LRisk in TCGA-HCC samples. We first calculate the FGA, FGL, and FGG scores to evaluate differences in chromosomal instability between the two risk groups. We found the LRisk group had significantly lower copy number loss or gain than HRisk, so it was obvious that LRisk had better chromosomal stability than HRisk (all, p < 0.001; Figure 8B). Next, we analyzed the copy number in different specific regions in LRisk and HRisk groups. The most frequent arm-level aberrations in the HRisk group identified were 13q, 11q, 4q, etc., for copy number loss, and significantly amplified regions in the HRisk were 1q, 8q, 3q, etc. (Figure 8A). Therefore, we decided to further explore the relationship between the copy number variations of specific DDR genes. The relationship between genes in eight core DDR pathways and their copy number alterations were calculated by Spearman analyses. The GISTIC calls of 38 DDR genes (R < −0.2 or R > 0.2) are shown in Figure 8C, and the corresponding expression of DDR genes are shown in Figure 8D. These results suggest that HRisk existed higher levels of copy number alterations, and both were associated with overexpression of DDR genes.