Prediction of early recurrence and response to adjuvant Sorafenib for hepatocellular carcinoma after resection

Datasets of HCC

We carefully screened four gene expression datasets (GSE14520, GSE87630, GSE54236, and GSE109211) with a total of 504 HCC patients from Gene Expression Omnibus (GEO, http://www.ncbi.ncbi,nlm.nih.gov/geo) and downloaded Liver Hepatocellular Carcinoma (LIHC) dataset containing 219 patients with relapse-free survival (RFS) status from The Cancer Genome Atlas (TCGA, https://portal.gdc.cancer.gov). In our research, the datasets from GEO were used as the discovery cohorts, while TCGA-LIHC was taken as the validation cohort. Among the discovery cohorts, three datasets (GSE14520, GSE87630, and GSE54236) were acquired on different sequencing platforms (Affymetrix, Illumina, and Agilent) and different sample backgrounds (America, Korea, and Italy). GSE109211 dataset included 140 patients with hepatic resection, of whom 67 patients received sorafenib as adjuvant therapy and 73 patients with Placebo. In addition, the complete clinical prognostic information of patients in GSE14520 and LIHC-TCGA were downloaded, respectively.

Identification of differentially expressed genes

The online analytical tool of GEO2R (https://www.ncbi.nlm.nih.gov/geo/geo2r) was used to identify Differentially Expressed Genes (DEGs) between tumor and normal tissues, and the volcano plots of genes were visualized using the R package ‘ggplot2’. Genes with fold change (FC) > 1.5 and adjusted p-value < 0.05 were reserved as DEGs (Ritchie et al., 2015). Furthermore, Venn diagram analyses were performed with R package ‘VennDiagram’ to obtain the overlapping DEGs, and the overlapping DEGs within all series were retained as credible DEGs. Besides, Functions Enrichment Analyses (GO-biological processes) and the clustering analysis of the credible DEGs were respectively performed using the R package ‘clusterProfiler’ and ‘pheatmap’.

Establishment and validation of LASSO Cox regression model

To identify DEGs related to early relapse (within 2 years after resection) in GSE14520, we used the Youden index method in receiver operating characteristic (ROC) analysis. According to the cutoff values, the expression level of each gene was converted into binary status (low expression status was equal to 0 and high expression status was equal to 1). Based on the binary status data of each gene, the univariate Cox regression analysis were conducted to further select the early relapse related DEGs. The predicted time was equal to 24 months, and the criteria were as follows: the area under the ROC curve (AUC) > 0.55 and p-value < 0.25, respectively. Finally, these selected DEGs with the binary status were further analyzed using the least absolute shrinkage and selection operator (LASSO) Cox regression which could be used to construct risk models by screening variables in high dimension data (Gui & Li, 2005; Shi et al., 2017). According to 10-fold likelihood cross-validation and ran up to 100 times, the most powerful prognostic genes were identified and the risk score model was constructed. The model = ${\sum }_{i=1}^{n}coeffcient\ast ExpressionstatusofDEGs\left(i\right)$. Based on this formula, patients were divided into high- or low-risk groups according to the optimal cutoff risk score value determined by ROC for predicting tumor early recurrence. Next, ROC curves were applied to evaluate the accuracy of the risk score model for predicting tumor early recurrence and response to sorafenib. Kaplan–Meier method and Log-Rank test were carried out to compare the recurrence and OS rate between two groups. Univariate and Multivariate Cox regression analyses (Zhang et al., 2018a) were performed for potential prognostic factors such as genders (female vs. male), metastasis (high vs. low), cirrhosis (yes vs. no), clinicopathological features (TNM, BCLC, CLIP stage) and risk score model. Furthermore, decision curve analysis (DCA), which was widely used to measure the clinical utility of specific models was performed (Vickers et al., 2008; Zhang et al., 2018b). The advantage of DCA is that it integrates the preferences of patients or policy-makers into analysis. If the net benefit of a model is greater than treating all and none patients, the model has clinical utility.

Statistical analysis

Kaplan–Meier plots and log-rank tests were performed using GraphPad Prism (version 7.0). And ROC curves, Univariable or multivariable Cox regression analyses, LASSO Cox regression, DCA plots and Calibration curve analysis were performed using R (version 3.6.1) packages ‘pROC’, ‘survival’, ‘glmnet’, ‘devtools’, ‘DecisionCurve’, ‘rmda’, respectively. The values of p < 0.05 in all tests were considered statistically significant.

Identification and analysis of differentially expressed genes in HCC

Detailed information of all eligible datasets as described in Table 1. After analyzing these datasets using GEO2R, 2859, 2399, and 3226 DEGs were obtained and were displayed by volcano plots (Figs. 1A–1C). Then the overlapping analysis of DEGs in three eligible datasets (GSE14520, GSE87630, and GSE54236) was conducted, and a commonly dysregulated gene set containing 520 genes in HCC was obtained (Fig. 1D). Similarly, in the GSE109211 dataset, 7034 and 6874 DEGs were obtained from the adjuvant treatment group with sorafenib and the placebo group, respectively (Figs. 1E and S1A). Meanwhile, the setdiff function analysis of DEGs in GSE109211 showed that 964 genes were associated with sorafenib adjuvant treatment but not with placebo (Fig. S1B). Finally, 34 overlapping DEGs between HCC-related dysregulate genes and sorafenib-related genes were identified, which were considered to be commonly dysregulated in HCC and were related to the response to adjuvant sorafenib in patients, who received liver resection (Fig. 1F). Clustering analysis of these genes was conducted by a heat map, including 14 up-regulated and 20 down-regulated genes (Fig. 1G).

To explore the potential biological roles of these 34 DEGs, Gene Ontology (GO) enrichment analysis were performed. As shown in Table S1, most genes were enriched in the biosynthetic process (such as carboxylic acid, organic acid, and monocarboxylic acid). These results suggested that the 34 genes were closely associated with HCC processing.

Construction of an early relapse model from the GSE14520 cohort

Based on these 34 candidate genes and complete clinical prognostic information of GSE14520, the optimal cutoff and AUC values for predicting tumor early recurrence of each gene were determined by time-dependent ROC curve, and the p-values were determined by univariate Cox regression analysis, which could independently predict the tumor early recurrence (The optimal cutoff values and p-values were shown in Table S2). Finally, 15 DEGs with AUC > 0.55 and p-value < 0.25 were selected to construct the risk model using LASSO Cox Regression analysis.

Finally, a multi-mRNA model containing 11 genes with nonzero coefficients was generated at the lambda.min value (0.03786) chosen by 10-fold cross-validation and 100 times (Fig. 2A). These mRNAs included ubiquitin associated protein 2 like (UBAP2L), cell division cycle 20 (CDC20), midkine neurite growth-promoting factor 2 (MDK), signal sequence receptor subunit 2 (SSR2), centromere protein M (CENPM), cytochrome P450 family 1 subfamily A member 2 (CYP1A2), vasoactive intestinal peptide receptor 1 (VIPR1), acetyl-CoA acyltransferase 2 (ACAA2), cytochrome P450 family 3 subfamily A member 4 (CYP3A4), cell death-inducing DFFA-like effector b (CIDEB), and N-terminal EF-hand calcium-binding protein 3 (NECAB3) (Fig. 2B). The equation of 11-mRNA risk score model was: The model = (0.4044 * status of UBAP2L) + (0.1919 * status of CDC20) + (0.1714 * status of MDK) + (0.1375 * status of SSR2) + (0.0566 * status of CENPM) – (0.0538 * status of CYP1A2) – (0.1371 * status of VIPR1) – (0.1488 * status of ACAA2) – (0.1812 * status of CYP3A4) – (0.3578 * status of CIDEB) – (0.5150 * status of NECAB3).

Evaluation and validation of the model in the discovery cohort

According to the multi-mRNA model, the score for each patient in GSE14520 and GSE109211 was calculated. The ROC curve showed that the AUC values of the score for predicting tumor early recurrence and response to sorafenib were 0.746 and 0.722, respectively (Figs. 3A and 3B), and the score of the sorafenib resistance was significantly higher than that of responders (p = 0.0024, Fig. 4B). The optimal cutoff value for relapse was –0.058, with a sensitivity of 0.659 and specificity of 0.796 (Fig. 3A). Based on this cutoff value (–0.058), we divided the patients in GSE14520 into two subgroups (high-risk group and low-risk group, Fig. S4A). Kaplan–Meier curves showed that patients with a high-risk score had worse recurrence rate (2-year: p < 0.0001, hazard ratio (HR): 4.658, confidence interval 95% CI [2.895–7.495]; 5-year: p < 0.0001, HR: 3.251, 95% CI [2.155–4.904]), and the median time to recurrence was 13.5 months (Fig. 3C and Fig. S2C). After adjusting by clinical factors (gender, metastasis, cirrhosis, TNM stage, BCLC stage, and CLIP stage) determined using univariate Cox regression analyses, the risk score model still had an excellent ability to independently predict the early recurrence of HCC, the risk model and other independent predictors were presented as a nomogram (Fig. 3D). Meanwhile, decision curve analysis (DCA) showed that if the threshold probability is higher than 0.2 and less than 0.5, the risk score model had a larger net benefit versus clinicopathologic nomogram (Fig. 3F). In addition, the DCA curve showed favorable agreement between risk model predicted recurrence probability and actual recurrence probability, and clinical impact curve also revealed that the patients could be significantly benefited from this model (Figs. 3E and 3G).

Furthermore, to explore the potential clinical value of the risk score model among patients with the same tumor stage category, Kaplan–Meier curves analysis was performed. For the GSE14520 cohort, when stratified by clinicopathologic features (BCLC staging and TNM staging), patients with a high-risk score displayed higher recurrence rates, especially those patients with tumors at the early stage (BCLC 0 and A, TNM I and II) (p < 0.001, Figs. S2A, S2B and S2D–S2F), and regardless the risk of tumor metastasis (all p < 0.0001, Figs. S2G and S2H). Even if the patients with advanced tumor stage (such as TNM Stage III) were excluded, the risk model still had high predictive value demonstrated by ROC (AUC value: 0.725, Fig. S2I). In addition, to verify the prognostic value of the risk model, we further conducted the time-dependent ROC curve for predicting survival and Kaplan–Meier curve analyses. It demonstrated that the AUC value was 0.715 (Fig. S3A) and patients with the high-risk score had remarkably worse survival rates (p < 0.0001, HR:3.811, 95% CI [2.365–6.142], Fig. S3B). Multivariate Cox regression analysis demonstrated that the risk model was still significantly related to survival (p < 0.001, Fig. S3C). Moreover, DCA revealed that the risk score model had a similar survival benefit to BCLC stage (Fig. S3D); and Calibration curve analysis showed the predicted probability survival value was consistent with actual survival (Fig. S3E). In conclusion, the risk score model was a reliable model with clinically and statistically significant for early-RFS patients and for predicting response to adjuvant sorafenib.

External validation of the risk score model in TCGA-LIHC

Since the risk score model can be used to predict HCC recurrence and the response of patients to adjuvant sorafenib in the GSE14520 and GSE109211 cohort, then we further validated the predictive ability of the model in another HCC cohorts (TCGA-LIHC). According to the same formula, the risk score for each patient was calculated and the time-dependent ROC curve analysis was conducted. The AUC values of the model for predicting the tumor recurrence and response to adjuvant sorafenib were 0.618 (95% CI [0.522–0.713]) and 0.708 (95% CI [0.467–0.949]), respectively (Figs. 4A and 4B). Similar to the GSE14520 cohort, patients with high risk had a significantly higher recurrence rate and shorter median recurrence time (18.3 months) (all p < 0.001, Figs. 4C, 4E and Fig. S4C). And these sorafenib responders have a lower risk score (p = 0.0411, Fig. S4D). Moreover, patients with low risk had much better survival rate (p = 0.0016, HR: 2.331, 95% CI [1.288–4.218], Fig. 4F). Multivariate Cox regression analysis also demonstrated that the risk model was significantly related to the tumor recurrence (p < 0.05, Fig. 4D). In addition, Kaplan–Meier curves analysis was further performed in subsets of patients with different TNM stages. Interestingly, patients with a high-risk score also displayed higher recurrence rates in the validation cohorts with the different TNM stages (Figs. 4G–4I). In summary, the novel risk score model constructed using bioinformatic methods could effectively predict HCC early recurrence and the benefits from adjuvant sorafenib in patients with surgical resection.