We obtained an RNA-seq dataset (n = 377), which was used as a training set to build the model, and the corresponding clinical information of HCC patients from The Cancer Genome Atlas (TCGA) using the UCSC Xena browser (https://xenabrowser.net/). The validation set was based on a second RNA-seq dataset (n = 243) downloaded from the International Cancer Genome Consortium (ICGC) (https://dcc.icgc.org) and a microarray dataset (GSE14520, n = 225) from GEO database (http://www.ncbi.nlm.nih.gov/geo). We extracted 232 autophagy-related genes from the Human Autophagy Database (HADb, http://www.autophagy.lu/index.html) and 394 from the GO_AUTOPHAGY gene set in the Molecular Signatures Database v7.1 (MSigDB, http://software.broadinstitute.org/gsea/msigdb). Our autophagy-related gene set was formed by the integration of these two gene sets.

Specimens of HCC patients who survived less than one month or whose clinical data were incomplete were not included. We removed the data of normal tissue samples and only kept the data of the primary tumor. When multiple specimens were taken from the same patient, the average gene expression value was used to represent the patient’s gene expression level. Only the sequencing data of autophagy-related genes were retained. When the same gene was matched by multiple probes, we used the average expression value of multiple probes to indicate the expression level of the gene. For the RNA-seq data from TCGA database, we excluded HCC samples in which more than half of the gene probe expression values were zero. The expression profiles of common autophagy-related genes were screened from the three data sets.

We compared the expression values of autophagy-related genes in each sample to obtain the score of each ARGP. The expression values of autophagy-related genes in each sample were compared in pairs to calculate the score of each ARGP. In a pairwise comparison, if the previous value is greater than the next value, the output is 1, and if it is not, the output is 0. We excluded ARGPs that scored 0 or 1 in more than 90% of the samples in each dataset, and the remaining ARGPs were used to establish the prognostic model for HCC. First, we performed univariate Cox regression analysis using the R package “survival” to select gene pairs that are related to the overall survival of HCC patients in TCGA. Differences with p < 0.001 were considered statistically significant. To minimize the risk of overfitting, we used “glmnet” R package to conduct LASSO penalized Cox regression (3,000 iterations) to calculate the frequency of the models. The gene pair model with the highest frequency among the iterations was used to establish a prognostic model. Stepwise multivariate Cox regression analysis was performed.

After calculating the risk score in every dataset, the patients were classified into high-risk and low-risk groups according to the median value of the risk score. Kaplan–Meier survival analysis (p < 0.05) was used to analyze the over survival (OS) of the high-risk and low-risk groups. After drawing ROC curves for 1, 3 and 5 years, we used the area under curve (AUC) value to verify the accuracy and sensitivity of this model. The closer the AUC value is to 1, the better the predictive effect of the prognostic model. To perform multivariate Cox regression analysis, available clinical and pathological data were integrated with the ARGP signature. Tumor stage, grade, age, and sex were regarded as continuous variables. In addition, we selected three representative prognostic gene models for HCC. Our model was compared with these existing models using the 5-year multiple ROC curves. The respective AUC values were used to estimate the prognostic accuracy of each signature.

In order to further reveal the biological mechanisms through which the identified autophagy-related genes contribute to the development of HCC, we used the MSigDB hallmark gene set (h.all.v7.1.symbols.gmt) to run gene set enrichment analysis (GSEA). We used an FDR value < 0.25, a nominal (NOM) p < 0.05, and | NES | > 1 as the screening criteria to identified signaling pathways that are highly related to the model genes.

Independent prognostic factors that are highly correlated with OS in HCC patients (p < 0.05) were screened out using univariate and multivariate Cox regression analyses. We then integrated these independent prognostic factors using the R package “RMS” and constructed the predictive nomogram and corresponding calibration diagram for 1, 3, and 5 years. The calibration maps were verified by calibration and discrimination. The expected possibility of collinearity was plotted graphically as an observable indicator to assess the alignment of the nomogram. The closer the calibration curve was to the reference line (diagonal line), the better the predictive effect of the nomogram.

All statistical analyses were performed using R software (version 3.6.3, https://www.r-project.org/). The OS of the HCC patients in the low- and high-risk groups was compared using the log-rank test, and the Kaplan–Meier survival curves were drawn using the R package “survminer” (version: 0.4.6). The gene pair prognostic signature was established based on the LASSO Cox regression algorithm using the R package “glmnet” (version: 3.0.2). ROC curves and multiple ROC curves were drawn using the R packages “survivalROC” and “timeROC”, respectively.

After eliminating duplicate genes from HADb and MSigDB, an autophagy-related gene set comprising 527 genes was obtained. As shown in Table 1, this study included 808 HCC patients from three cohorts. ARGPs were constructed using a total of 269 autophagy-related genes that are represented in all three data sets.

We removed ARGPs with a score of 0 or 1 in more than 90% of the samples in all datasets, resulting in 1885 ARGPs. We used univariate Cox regression analysis to screen 117 prognostic ARGPs that were significantly associated with overall survival (p < 0.001), and established a prognostic gene model of ARGP using Lasso penalty score Cox regression in the TCGA dataset. After multivariate Cox regression analysis, 9 ARGPs were selected to construct the most stable prognostic signatures, and the corresponding coefficients were used to calculate the risk score for our datasets. Details of the 9-ARGP prognostic model are listed in Table 2.

The 9-ARGP signature in the TCGA dataset reflected the postoperative prognosis of patients very well (Figure 1). We calculated the risk score of each HCC patient in TCGA according to the prognostic characteristics of autophagy, and then divided the 346 cases into high- and low-risk groups based on the median score. As shown in Figure 1A, the low-risk group had a significantly lower mortality rate than the high-risk group (95%CI: 8.09–21.07, p < 0.0001). We evaluated the specificity and sensitivity of the prognostic model using time-dependent ROC curve analysis. The AUC values for 1, 3, and 5 years after surgery reached 0.812, 0.786, and 0.828, respectively, which demonstrated that our ARGP prognostic signature has a promising predictive ability (Figure 1B). The distribution of the autophagy-related prognostic model for patients in the TCGA data set is shown in Figure 1C,D.

In order to further verify its predictive power, we applied the prognostic signature to the ICGC database (containing 243 HCC cases) and the GSE14520 database (containing 225 HCC cases) for analysis. According to the median risk value calculated using the 9-ARGP prognostic signature, HCC patients in the two databases were assigned into high- and low autophagy-based risk groups, respectively. Consistent with the conclusions obtained using the training set, the OS of the high-risk groups in the two validation datasets was significantly lower than that of the low-risk group (p < 0.05) (Figure 2A,E). In the ICGC cohort, the AUC values of the prognostic model were 0.842 at 1 year, 0.725 at 2 years, and 0.727 at 3 years (Figure 2B), while in the GSE14520 cohort they were 0.625 at 1 year, 0.667 at 2 years, and 0.644 at 3 years (Figure 2F). Figure 2C,D show the distribution of risk scores corresponding to gene expression levels in the ICGC cohort, while Figure 2G,H shows the corresponding data for the GSE14520 cohort. In univariate Cox regression analysis, TNM staging and the ARGP signature risk score were significantly related to the OS in the three cohorts (HR > 1.00, p < 0.05). After correcting for age, gender, grade and TMN staging in multivariate Cox regression analysis, the ARGP signature risk score was still significantly associated with the OS as an independent prognostic factor in the TCGA dataset (HR: 1.25, 95%CI: 1.19–1.31, p < 0.0001), ICGC dataset (HR: 1.11, 95%CI: 1.06–1.17, p< 0.001) and GSE14520 dataset (HR: 1.03, 95%CI: 1–1.17, p = 0.07) (Table 2).

The gene pairs with the largest coefficients were STAM/TP53, PLOD2/CDKN1B and NTHL1/BLCR5, since the large coefficients indicate that they have the greatest influence on the model. We examined the expression levels of these genes in a normal liver cell line and in several common HCC cell lines. As shown in Figure 3A, the qPCR results showed that the values of STAM/TP53 and PLOD2/CDKN1B in the HCC cell lines HepG2 and Hep3B were significantly higher than in the normal liver cell line HL7702, which further supported the significance of our model. The coefficient of the gene pair NTHL1/BLCR5 was −0.52, and was significantly lower in the HCC cell lines HepG2 and Hep3B than in the normal liver cell line HL7702, which is also consistent with this conclusion. At the protein level, we found that the expression of PLOD2 in liver cancer cell lines was higher than in normal liver cells, and it was also higher in non-invasive liver cancer cells (Figure 3B).

Further experiments were performed on the representative gene pair Plod2/CDKN1B. We overexpressed these genes in MHCC97-H cells, and the protein levels of PLOD2/CDKN1B were confirmed to be increased after transfection with pcDNA-PLOD2 (oe-PLOD2) and pcDNA-CDKN1B (oe-CDKN1B) (Figure 3C). The results of the transwell assay showed that PLOD2 can promote HCC migration, while CDKN1B had the opposite effect. Since their coefficient is greater than zero, this experimental result was consistent with the previous conclusion (Figure 3D).

We verified the expression levels of the representative gene pair PLOD2 /CDKN1B in 45 liver cancer samples, and combined with the OS and RFS of the corresponding patients, we found that PLOD2 /CDKN1B can better predict the prognosis of the patients. The results are shown in Figure 7. Samples were collected from surgical biopsies of patients who underwent radical resection of liver cancer without preoperative treatment at Tongji Hospital in Wuhan, China, between 2015 and 2018. The Ethics Committee of Wuhan Tongji Hospital authorized this study on patient tissues with written informed consent of the patients.

Our ARGP prognostic marker was compared with three published representative gene prognostic markers (Lin et al., 2018; Long et al., 2018; Wang et al., 2020) using ROC curves for 5-year OS. All the data for validation were derived from TCGA. As shown in Figure 4, the AUC value was 0.828 for our prognostic signature, which was obviously more predictive and accurate than the existing autophagy-related signature (AUC = 0.628), the four-gene prognostic model (AUC = 0.688) and the nine immune-related gene model (AUC = 0.595).

Multivariate Cox regression analysis showed that only TNM stage and the ARGP signature were significant independent prognostic factors for OS (Figure 5A). We attempted to provide a method to more intuitively and accurately predict the survival of HCC patient, which can aid individual clinical decision-making and selection of treatment options. Therefore, a predictive nomogram was constructed based on multivariate Cox regression analysis and combined with two independent prognostic factors (Figure 5B). The scores of each independent prognostic factor were calculated according to the different degree of influence of each independent prognostic factor on the clinical outcomes of the patients, after which the scores were summed up to obtain the total score. Finally, the 1, 3, and 5-year survival rates were predicted based on the functional relationship between the total score and the survival rate. According to the calibration curves of the 1-, 3-, and 5-year nomograms, which were all close to the optimal prediction curve, the predicted OS rate was highly consistent with the actual observed values (Figure 5C–E).

We performed GSEA in the high- and low-risk groups from the TCGA cohort, divided according to the median risk value. A total of 12 cancer hallmark gene sets were identified in the high-risk group (Figure 6A). Some of these pathways are “MYC_TARGETS”, “GLYCOLYSIS” and “DNA_REPAIR”, indicating that these signaling pathways are closely related to the progression of HCC. To make the results of enrichment analysis more intuitive, we visualized the significance, the number of included genes and the enrichment score in a bubble chart with different colors, sizes and locations (Figure 6B).