# Identification and Validation of a Prognostic Model Based on Three Autophagy-Related Genes in Hepatocellular Carcinoma

Authors: Fanbo Qin; Junyong Zhang; Jianping Gong; Wenfeng Zhang
Journal: BioMed Research International
 (2021)

---

## Abstract

Background. Accumulating studies have demonstrated that autophagy plays an important role in hepatocellular carcinoma (HCC). We aimed to construct a prognostic model based on autophagy-related genes (ARGs) to predict the survival of HCC patients. Methods. Differentially expressed ARGs were identified based on the expression data from The Cancer Genome Atlas and ARGs of the Human Autophagy Database. Univariate Cox regression analysis was used to identify the prognosis-related ARGs. Multivariate Cox regression analysis was performed to construct the prognostic model. Receiver operating characteristic (ROC), Kaplan-Meier curve, and multivariate Cox regression analyses were performed to test the prognostic value of the model. The prognostic value of the model was further confirmed by an independent data cohort obtained from the International Cancer Genome Consortium (ICGC) database. Results. A total of 34 prognosis-related ARGs were selected from 62 differentially expressed ARGs identified in HCC compared with noncancer tissues. After analysis, a novel prognostic model based on ARGs (PRKCD, BIRC5, and ATIC) was constructed. The risk score divided patients into high- or low-risk groups, which had significantly different survival rates. Multivariate Cox analysis indicated that the risk score was an independent risk factor for survival of HCC after adjusting for other conventional clinical parameters. ROC analysis showed that the predictive value of this model was better than that of other conventional clinical parameters. Moreover, the prognostic value of the model was further confirmed in an independent cohort from ICGC patients. Conclusion. The prognosis-related ARGs could provide new perspectives on HCC, and the model should be helpful for predicting the prognosis of HCC patients.


---

## Body

## 1. Introduction

Hepatocellular carcinoma (HCC) is the fifth most common cancer worldwide and the third most common cause of cancer mortality [1]. The main cause of HCC is nonalcoholic fatty liver disease (NAFLD) and hepatitis virus infection including chronic hepatitis B virus (HBV) infection and chronic hepatitis C virus (HCV) infection [2]. In America, although the 5-year survival rate of HCC patients increased by 90% between 2006 and 2012 compared to 1992-1999, the overall prognosis of HCC remains poor [3]. Therefore, establishing an effective prognostic model can provide new guidance for clinical management. However, conventional clinical parameters used to predict clinical outcome often have some limitations given the heterogeneity of HCC [4]. Therefore, the establishment of a novel effective prognostic model must be based on the molecular heterogeneity of HCC.Autophagy is a lysosomal degradation pathway for unwanted, damaged, and defective intracellular components and is essential for survival, differentiation, development, and homeostasis [5]. In recent decades, accumulating evidence has demonstrated that abnormal expression of autophagy-related genes (ARGs) is involved in the development of various cancers [6, 7]. However, the role of ARGs in cancer remains inconclusive. In general, at the early stage of liver cancer, autophagy protects cells from carcinogenesis. At an advanced stage, however, autophagy can promote cancer progression [8, 9]. For liver cancer, dysregulated autophagy is dependent on changes in ARG expression [10]. Thus, it is of special significance to identify biomarkers from ARGs that are associated with the prognosis of patients with cancer and to calculate their influence on the prognosis model. Recent reports have shown that prognostic models based on ARGs provide better prediction of clinical outcomes for patients with cancer, such as thyroid cancer and bladder cancer [11, 12]. However, studies on prognostic models based on ARGs in HCC are limited. Thus, our aims in this study were to construct a prognostic model based on ARGs to predict HCC patient survivals.In this study, we identified the 62 differentially expressed ARGs in HCC based on The Cancer Genome Atlas (TCGA) database. Enrichment analysis of differentially expressed ARGs may help to identify the potential molecular mechanisms in autophagy. In addition, we performed univariate Cox regression analysis to identify 34 prognosis-related ARGs and conducted multivariate Cox regression analysis to select 3 key genes (PRKCD, BIRC5, and ATIC) to construct a prognostic model. Univariate and multivariate Cox regression analyses confirmed that the risk score calculated by the model formula was an independent risk factor for patient survival. Survival analysis and ROC curve analysis demonstrated that the model had good performance in predicting prognosis. Moreover, we confirmed the reliability of this model by analyzing another independent data cohort obtained from the International Cancer Genome Consortium (ICGC) database. In summary, the prognostic model we constructed in this study is not only helpful in identifying the potential mechanism of autophagy but also important for developing an effective prediction tool for HCC patients.


## 2. Materials and Methods

### 2.1. Data Acquisition

The Human Autophagy Database (HADb,http://www.autophagy.lu/index.html) is the first human autophagy-dedicated database. It is a public repository containing information about the human genes involved in autophagy described to date [13]. All ARGs for subsequent analysis were downloaded from this database. The expression data (type: FPKM) and clinical information of HCC were obtained from The Cancer Genome Atlas (TCGA, https://portal.gdc.cancer.gov/) database to identify differentially expressed ARGs and construct the prognostic model, and independent datasets obtained from the International Cancer Genome Consortium (ICGC, https://icgc.org) were used to validate the model. Data were downloaded from the publicly available database; hence, additional ethical approval was not necessary for this study.


### 2.2. Differentially Expressed ARGs and the Enrichment Analysis

The Wilcoxon rank-sum test was used for the identification of differentially expressed ARGs in HCC samples compared with noncancer tissue samples with the criteria of adjustedp value < 0.05 and ∣log2FC∣>1. The Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses were performed using the R clusterProfiler package, and the GO and KEGG results were visualized using the GOplot package. A p value < 0.05 was set as the cutoff criterion for both GO and KEGG functional analyses.


### 2.3. Construction and Validation of the Prognostic Model

Univariate Cox regression analysis evaluated the association between the expression of differentially expressed ARGs and patient overall survival (OS). Genes with ap value less than 0.05 were considered prognosis-related ARGs and then entered into a stepwise multivariate Cox regression analysis tested by AIC (Akaike Information Criterion, assessing the goodness of fit of a statistical model) to identify the predictive model. Then, a prognostic model was constructed, and the formula of the risk score was as follows: riskscore=theexpressionofgene1×regressioncoefficientofgene1+theexpressionofgene2×regressioncoefficientofgene2+⋯+theexpressionofgenen×regressioncoefficientofgenen.Based on the risk score, HCC patients were classified into low- or high-risk groups. Survival curves were generated using the Kaplan-Meier method, and two-sided log-rank tests were employed to compare the differences in OS between the low- and high-risk groups. Univariate Cox regression analysis and multivariate Cox regression analysis were conducted to explore whether the risk score could be an independent indicator of OS in the TCGA data cohort of HCC patients. The sensitivity and specificity of the prognostic model to predict the clinical outcome of HCC patients were analyzed by calculating the area under the curve (AUC) of the receiver operating characteristic (ROC) curve. Moreover, an independent data cohort of the ICGC database was used to further confirm the predictive value of the model.


### 2.4. Statistical Analysis

R 3.6.1 (https://www.r-project.org) was utilized for plot and statistical analysis. The Wilcoxon rank-sum test was used to identify differentially expressed ARGs. Univariate Cox regression analysis was performed to estimate the prognosis-related ARGs. Multivariate Cox regression analysis was used to construct the prognostic model. An independent t-test was used to analyze the associations between risk score and conventional clinical parameters. ROC analysis was performed to test the sensitivity and specificity of the model. The survival curve was plotted using the survival and survminer packages of R. The forest maps were plotted by the forestplot package of R. The survivalROC package of R was used to generate ROC curves and AUC values to calculate according to ROC curves.


## 2.1. Data Acquisition

The Human Autophagy Database (HADb,http://www.autophagy.lu/index.html) is the first human autophagy-dedicated database. It is a public repository containing information about the human genes involved in autophagy described to date [13]. All ARGs for subsequent analysis were downloaded from this database. The expression data (type: FPKM) and clinical information of HCC were obtained from The Cancer Genome Atlas (TCGA, https://portal.gdc.cancer.gov/) database to identify differentially expressed ARGs and construct the prognostic model, and independent datasets obtained from the International Cancer Genome Consortium (ICGC, https://icgc.org) were used to validate the model. Data were downloaded from the publicly available database; hence, additional ethical approval was not necessary for this study.


## 2.2. Differentially Expressed ARGs and the Enrichment Analysis

The Wilcoxon rank-sum test was used for the identification of differentially expressed ARGs in HCC samples compared with noncancer tissue samples with the criteria of adjustedp value < 0.05 and ∣log2FC∣>1. The Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses were performed using the R clusterProfiler package, and the GO and KEGG results were visualized using the GOplot package. A p value < 0.05 was set as the cutoff criterion for both GO and KEGG functional analyses.


## 2.3. Construction and Validation of the Prognostic Model

Univariate Cox regression analysis evaluated the association between the expression of differentially expressed ARGs and patient overall survival (OS). Genes with ap value less than 0.05 were considered prognosis-related ARGs and then entered into a stepwise multivariate Cox regression analysis tested by AIC (Akaike Information Criterion, assessing the goodness of fit of a statistical model) to identify the predictive model. Then, a prognostic model was constructed, and the formula of the risk score was as follows: riskscore=theexpressionofgene1×regressioncoefficientofgene1+theexpressionofgene2×regressioncoefficientofgene2+⋯+theexpressionofgenen×regressioncoefficientofgenen.Based on the risk score, HCC patients were classified into low- or high-risk groups. Survival curves were generated using the Kaplan-Meier method, and two-sided log-rank tests were employed to compare the differences in OS between the low- and high-risk groups. Univariate Cox regression analysis and multivariate Cox regression analysis were conducted to explore whether the risk score could be an independent indicator of OS in the TCGA data cohort of HCC patients. The sensitivity and specificity of the prognostic model to predict the clinical outcome of HCC patients were analyzed by calculating the area under the curve (AUC) of the receiver operating characteristic (ROC) curve. Moreover, an independent data cohort of the ICGC database was used to further confirm the predictive value of the model.


## 2.4. Statistical Analysis

R 3.6.1 (https://www.r-project.org) was utilized for plot and statistical analysis. The Wilcoxon rank-sum test was used to identify differentially expressed ARGs. Univariate Cox regression analysis was performed to estimate the prognosis-related ARGs. Multivariate Cox regression analysis was used to construct the prognostic model. An independent t-test was used to analyze the associations between risk score and conventional clinical parameters. ROC analysis was performed to test the sensitivity and specificity of the model. The survival curve was plotted using the survival and survminer packages of R. The forest maps were plotted by the forestplot package of R. The survivalROC package of R was used to generate ROC curves and AUC values to calculate according to ROC curves.


## 3. Results

### 3.1. Differentially Expressed Autophagy-Related Genes between HCC and Normal Tissues

To date, HADb datasets include a total of 232 ARGs that have been described as being involved in autophagy. A total of 50 nontumor tissues and 374 tumor tissue samples with mRNA expression profiles and 235 HCC clinical follow-up data points were downloaded from the TCGA database (Table1). Compared with the noncancer tissue samples, 62 differentially expressed ARGs were identified in the HCC samples (Table 2). The results of hierarchical cluster analysis showed that the HCC samples could be clearly distinguished from the normal tissues based on the differentially expressed ARGs (Figure 1(a)), and the volcano plot showed that there were 4 downregulated genes and 58 upregulated genes among these ARGs (Figure 1(b)). Furthermore, a scatter plot was generated to visualize the expression levels of differentially expressed ARGs between HCC and normal tissue (Figure 1(c)).Table 1
General characteristics of HCC patients included in the present study (data downloaded from the TCGA database).
ParametersPatients (N=235)Number%Age≤65 y16469.79>65 y7130.21GenderMale16168.51Female7431.49GradeG1/213256.17G3/410343.83Pathologic stageI/II16369.36III/IV7230.64T stageT1/216369.36T3/47230.64N stageN023198.3N1-341.7M stageM023198.3M141.7Table 2
Characteristics of the differentially expressed autophagy-related genes by using the Wilcoxon rank-sum test (tumor vs. normal).
Gene symbollogFCRegulationp valueFDRBIRC54.8146253Up2.35E−284.82E−26PEA151.83411226Up3.57E−273.66E−25RAB241.70056392Up8.37E−274.29E−25CLN31.7043542Up7.73E−274.29E−25HSP90AB11.63411809Up3.53E−261.45E−24MAPK31.37850531Up5.73E−261.68E−24CAPN101.49940395Up4.96E−261.68E−24CDKN2A4.65431982Up1.87E−254.26E−24RHEB1.21520158Up2.98E−256.10E−24RPTOR1.30199062Up4.43E−258.26E−24TSC11.56728821Up2.41E−244.11E−23ATG4B1.05882354Up4.52E−247.13E−23ATIC1.26886226Up1.56E−232.18E−22ITGA62.05950447Up1.59E−232.18E−22FKBP1A1.164278Up2.89E−233.37E−22TP734.10294233Up2.96E−233.37E−22ULK31.25487092Up3.47E−233.75E−22MLST81.04550722Up4.72E−234.84E−22CAPNS11.10213996Up5.81E−235.42E−22PARP11.34948718Up6.82E−236.08E−22BAX1.52707657Up7.52E−236.42E−22ATG101.0226174Up1.11E−229.11E−22DIRAS3-2.2142941Down2.04E−221.49E−21HGS1.59409057Up3.22E−222.28E−21TSC21.1673426Up4.51E−223.08E−21CANX1.10228843Up2.00E−211.32E−20WDR45B1.03610977Up4.09E−212.62E−20FOS-2.4192498Down6.68E−214.15E−20DAPK22.15868832Up1.19E−207.20E−20NPC11.40873288Up1.77E−209.80E−20GAPDH1.24841467Up5.02E−202.64E−19PELP11.02773866Up5.19E−202.66E−19HSPA51.07141566Up1.04E−195.18E−19HDAC11.0260015Up6.77E−192.89E−18VMP11.03638398Up1.06E−184.45E−18RUBCN1.05799379Up1.65E−186.76E−18BAK11.61591549Up1.55E−175.67E−17CD461.12531986Up1.83E−176.60E−17CASP31.01180546Up1.98E−176.98E−17DDIT31.50655255Up3.18E−171.10E−16CAPN21.4191918Up3.60E−171.23E−16SQSTM11.73640938Up9.02E−172.89E−16PRKCD1.41302501Up8.94E−162.82E−15CASP81.04011948Up1.52E−154.65E−15FOXO1-1.0757149Down4.28E−151.25E−14SPNS11.31587141Up2.15E−145.96E−14RB1CC11.01420177Up4.73E−141.29E−13ATG16L21.02505335Up2.00E−135.40E−13NRG1-1.280118Down5.44E−121.25E−11DRAM11.2000343Up1.11E−112.54E−11IKBKE1.86880957Up1.30E−112.93E−11TMEM742.33376628Up3.59E−117.83E−11ITGB42.37261675Up1.85E−103.88E−10RGS191.11586497Up5.15E−101.02E−09NRG22.87978625Up1.43E−092.68E−09ATG9B1.05236942Up5.60E−068.32E−06SPHK12.84999763Up8.75E−051.21E−04ITGA32.09689798Up4.48E−045.89E−04TP631.90806327Up2.37E−032.93E−03IRGM2.45185036Up3.98E−034.86E−03NKX2-32.56776198Up4.50E−035.46E−03HSPB82.02223664Up9.57E−031.13E−02logFC: log(fold change); FDR: false discovery rate.Figure 1
The differentially expressed autophagy-related genes (ARGs) between normal and hepatocellular carcinoma (HCC) tissues. (a) The hierarchical clustering of differentially expressed ARG expression levels. (b) The volcano plot for the differentially expressed ARGs. Red indicates high expression, green indicates low expression, and black indicates no difference between HCC and normal tissue. (c) The boxplot for the expression patterns of 62 significant differentially expressed ARGs in HCC and normal tissue.
(a)(b)(c)


### 3.2. Enrichment Analysis of the Differentially Expressed ARGs

To further investigate the biological functions of the differentially expressed ARGs in HCC, we performed GO term function and KEGG pathway enrichment analyses for these genes (Table3). For biological processes (BPs), these ARGs participate in autophagy, processes utilizing autophagic mechanisms, macroautophagy, and the regulation of autophagy. In terms of cellular components (CC), the genes were related to autophagosomes, phagophore assembly sites, vacuolar membranes, and autophagosome membranes. Changes in molecular function (MF) were enriched in protein kinase regulator activity, kinase regulator activity, cysteine-type endopeptidase activity, and BH domain binding (Figure 2, Table 3). Furthermore, KEGG pathway analysis showed that these genes mostly participated in autophagy-animal, apoptosis, platinum drug resistance, and longevity regulating pathways (Figures 3(a) and 3(b); Table 4).Table 3
GO analysis of differentially expressed autophagy-related genes.
CategoryIDTermCountp valueBiological processGO:0006914Autophagy318.54E−33Biological processGO:0061919Process utilizing autophagic mechanism318.54E−33Biological processGO:0016236Macroautophagy232.05E−26Biological processGO:0010506Regulation of autophagy181.38E−17Biological processGO:0016241Regulation of macroautophagy143.91E−16Biological processGO:0070997Neuron death131.53E−10Biological processGO:2001233Regulation of apoptotic signaling pathway137.96E−10Biological processGO:0097193Intrinsic apoptotic signaling pathway112.51E−09Biological processGO:0071900Regulation of protein serine/threonine kinase activity137.65E−09Biological processGO:0042475Odontogenesis of dentin-containing tooth71.68E−08Cellular componentGO:0005776Autophagosome71.26E−08Cellular componentGO:0000407Phagophore assembly site51.38E−08Cellular componentGO:0005774Vacuolar membrane114.53E−08Cellular componentGO:0000421Autophagosome membrane42.83E−06Cellular componentGO:0098589Membrane region81.05E−05Cellular componentGO:0005770Late endosome71.43E−05Cellular componentGO:0005765Lysosomal membrane81.74E−05Cellular componentGO:0098852Lytic vacuole membrane81.74E−05Cellular componentGO:0101031Chaperone complex33.25E−05Cellular componentGO:0045121Membrane raft77.15E−05Molecular functionGO:0019887Protein kinase regulator activity74.93E−07Molecular functionGO:0019207Kinase regulator activity71.14E−06Molecular functionGO:0004197Cysteine-type endopeptidase activity51.55E−06Molecular functionGO:0051400BH domain binding34.96E−06Molecular functionGO:0008234Cysteine-type peptidase activity61.35E−05Molecular functionGO:0051087Chaperone binding51.72E−05Molecular functionGO:0031072Heat shock protein binding53.93E−05Molecular functionGO:0044389Ubiquitin-like protein ligase binding76.46E−05Molecular functionGO:0030295Protein kinase activator activity49.71E−05Molecular functionGO:0019209Kinase activator activity40.000135Figure 2
Gene Ontology (GO) enrichment analysis of the differentially expressed autophagy-related genes. The size of each dot represents the count of genes; the color represents thep.adjust.Figure 3
Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis of differentially expressed autophagy-related genes (ARGs). (a) The outer circle shows a scatter plot for each term of the logFC of the assigned gens. Red circles display upregulation, and blue ones display downregulation. (b) The heat map of the relationship between ARGs and enriched pathways. The color of each block depends on the logFC value.
(a)(b)Table 4
KEGG analysis of differentially expressed autophagy-related genes.
IDTermCountp valuehsa04140Autophagy-animal162.50E−16hsa04210Apoptosis102.02E−08hsa01524Platinum drug resistance75.52E−07hsa04211Longevity regulating pathway72.15E−06hsa04136Autophagy-other52.24E−06hsa04215Apoptosis-multiple species52.24E−06hsa04115p53 signaling pathway68.78E−06hsa04218Cellular senescence81.14E−05hsa04657IL-17 signaling pathway63.84E−05hsa04151PI3K-Akt signaling pathway101.19E−04Table 5
The three ARG-based prognostic models identified by multivariate Cox regression analysis.
Gene symbolCoefHRHR.95LHR.95Hp valuePRKCD0.31751.37371.00551.87670.0461BIRC50.24791.28141.01941.61080.0337ATIC0.43971.55231.00382.40060.0481


### 3.3. Identification of Prognostic-Related ARGs

To identify ARGs associated with the clinical outcomes of patients, univariate Cox regression analysis was applied with the criterion of ap value < 0.05. As shown in Figure 4, 34 genes were selected and significantly associated with the overall survival (OS) of HCC patients. The hazard ratio (HR) of each gene was calculated, and all 34 genes had an HR>1, indicating that all of these genes are risk factors for OS in HCC patients.Figure 4
Thirty-four differentially expressed autophagy-related genes with a prognostic value determined by univariate Cox regression.


### 3.4. Construction and Validation of the ARG-Based Prognostic Model

Based on the prognosis-related ARGs, multivariate Cox regression analysis was performed to construct the ARG-based prognostic model. Three genes, PRKCD, BIRC5, and ATIC, were finally selected to construct the model (Table5), and the risk score formula of the model was as follows: riskscore=0.3175×PRKCDexpression+0.4397×ATICexpression+0.2479×BIRC5expression.Figure 5
Classification of patients into high- and low-risk groups by median risk score. (a) Kaplan-Meier plot of patients in a low- or high-risk group (p=1.049e−04). (b) The number of patients in different risk groups.
(a)(b)To assess the performance of the prognostic model in predicting the clinical outcome of HCC patients, we calculated the risk score of each HCC patient and divided patients into high- or low-risk groups using the median risk score as the cutoff value. The survival curve indicated that the high-risk group (n=117) suffered significantly lower 3- and 5-year survival rates than the low-risk group (n=118) (Figure 5). The risk score distribution, survival status of each patient, and heat map of the three gene expression profiles in the TCGA database are shown in Figure 6. The results indicated that survival time decreased and the mortality rate increased as the risk scores increased. We further analyzed the associations between the expression of the three ARGs and clinical parameters in HCC patients. The results of the independent t-test are shown in Figure 7. We observed significant correlations between overexpression of PRKCD and advanced histological grade (p=0.001), advanced T stage (p=0.011), advanced pathologic stage (p=0.011), and female sex (p=0.032). High BIRC5 expression was closely related to advanced T stage (p=0.007), advanced pathologic stage (p=0.007), and advanced histological grade (p=3.352e−04). Overexpression of ATIC occurred in the advanced T stage (p=0.047), advanced pathologic stage (p=0.047), and advanced histological grade (p=0.002). In particular, we used the same method to analyze the correlations between the risk score calculated by the model and clinical parameters. As shown in Figures 7(k)–7(n), we found that a high-risk score was significantly related to younger patients (p=0.047), advanced T stage (p=0.002), advanced histological grade (p=0.001), and advanced pathologic stage (p=0.002).Figure 6
Prognostic risk score of the differentially expressed autophagy-related genes. (a) The risk score of each hepatocellular carcinoma patient. (b) The patient survival based on the risk score. (c) Heat map of the expression of three genes in high- and low-risk groups.
(a)(b)(c)Figure 7
Correlations between expression of three key genes and clinical parameters (the results ofp<0.05 are shown).
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)(k)(l)(m)(n)Then, we wanted to determine whether the risk scores or other conventional clinical parameters of HCC patients were independent risk factors for the OS of patients. In Figure8, univariate Cox analysis showed that advanced pathologic stage, advanced T stage, advanced M stage, and risk score were risk factors for OS. However, after adjusting for clinical parameters, only the risk score remained an independent prognostic indicator for HCC patients in multivariate analyses (HR=1.745, 95%CI=1.420‐2.145, p<0.001). Moreover, we plotted the ROC curve to compare the predictive value of the risk score with other clinical parameters, as is shown in Figure 9. The results indicate that the T stage (AUC=0.708) and pathological stage (AUC=0.703) have the highest predictive value among the conventional clinical parameters. However, the predictive value of the risk score (AUC=0.758) was higher than that of the T stage and pathological stage.Figure 8
Univariate and multivariate Cox regression analyses of OS in HCC patients. (a) Univariate Cox regression analysis. (b) Multivariate Cox regression analysis. Age: ≤65 vs. >65; gender: male vs. female; grade: G1/G2 vs. G3/G4; stage: I/II vs. III/IV; T stage: T1/T2 vs. T3/T4; M stage: M0 vs. M1; N stage: N0 vs. N1-3; risk score: high-risk score vs. low-risk score (median risk score as the cutoff value).
(a)(b)Figure 9
Receiver operating characteristic (ROC) curve analysis for the prognostic value of the prognostic model and other conventional clinical parameters. AUC: area under the curve.


### 3.5. Confirmation of the Prognostic Model Using an Additional Data Cohort

The International Cancer Genome Consortium (ICGC) database contains large-scale cancer genome studies in tumors from 50 different cancer types [14]. To further validate the predictive value of this model, we downloaded the data cohort (Liver Cancer-RIKEN, JP, https://dcc.icgc.org/releases/current/Projects/LIRI-JP), which contains expression profiles and clinical follow-up information from 232 HCC patients. We analyzed the correlation between the expression levels of three genes and the survival of HCC patients in both the TCGA and ICGC databases. Kaplan-Meier analysis results indicated that high expression of the three genes was significantly associated with inferior OS in HCC patients (Figure 10). Similarly, these patients were divided into low- or high-risk groups based on the risk score calculated by the prognostic model constructed based on the TCGA data. The Kaplan-Meier analysis of the two groups was significantly different (p=4.505e−06<0.001) and a similar trend was observed in TCGA (Figure 11(a)). In addition, the AUC calculated by the ROC curve was 0.739, which indicated that the model has good performance for predicting HCC patient survival (Figure 11(b)).Figure 10
The correlation between the expression level of three genes included in the prognostic model and survival of HCC: (a) ATIC expression level and patient survival (data from TCGA), (b) PRKCD expression level and patient survival (data from TCGA), (c) BIRC5 expression level and patient survival (data from TCGA), (d) ATIC expression level and patient survival (data from ICGC), (e) PRKCD expression level and patient survival (data from ICGC), and (f) BIRC5 expression level and patient survival (data from ICGC).
(a)(b)(c)(d)(e)(f)Figure 11
Kaplan-Meier and ROC curves for the ARG prognostic model in the cohort of the ICGC. (a) The differences between the low-risk (n=116) and high-risk (n=116) groups determined by the log-rank test (p=4.505e−06). (b) ROC curve analysis for survival prediction by the prognostic model.
(a)(b)


## 3.1. Differentially Expressed Autophagy-Related Genes between HCC and Normal Tissues

To date, HADb datasets include a total of 232 ARGs that have been described as being involved in autophagy. A total of 50 nontumor tissues and 374 tumor tissue samples with mRNA expression profiles and 235 HCC clinical follow-up data points were downloaded from the TCGA database (Table1). Compared with the noncancer tissue samples, 62 differentially expressed ARGs were identified in the HCC samples (Table 2). The results of hierarchical cluster analysis showed that the HCC samples could be clearly distinguished from the normal tissues based on the differentially expressed ARGs (Figure 1(a)), and the volcano plot showed that there were 4 downregulated genes and 58 upregulated genes among these ARGs (Figure 1(b)). Furthermore, a scatter plot was generated to visualize the expression levels of differentially expressed ARGs between HCC and normal tissue (Figure 1(c)).Table 1
General characteristics of HCC patients included in the present study (data downloaded from the TCGA database).
ParametersPatients (N=235)Number%Age≤65 y16469.79>65 y7130.21GenderMale16168.51Female7431.49GradeG1/213256.17G3/410343.83Pathologic stageI/II16369.36III/IV7230.64T stageT1/216369.36T3/47230.64N stageN023198.3N1-341.7M stageM023198.3M141.7Table 2
Characteristics of the differentially expressed autophagy-related genes by using the Wilcoxon rank-sum test (tumor vs. normal).
Gene symbollogFCRegulationp valueFDRBIRC54.8146253Up2.35E−284.82E−26PEA151.83411226Up3.57E−273.66E−25RAB241.70056392Up8.37E−274.29E−25CLN31.7043542Up7.73E−274.29E−25HSP90AB11.63411809Up3.53E−261.45E−24MAPK31.37850531Up5.73E−261.68E−24CAPN101.49940395Up4.96E−261.68E−24CDKN2A4.65431982Up1.87E−254.26E−24RHEB1.21520158Up2.98E−256.10E−24RPTOR1.30199062Up4.43E−258.26E−24TSC11.56728821Up2.41E−244.11E−23ATG4B1.05882354Up4.52E−247.13E−23ATIC1.26886226Up1.56E−232.18E−22ITGA62.05950447Up1.59E−232.18E−22FKBP1A1.164278Up2.89E−233.37E−22TP734.10294233Up2.96E−233.37E−22ULK31.25487092Up3.47E−233.75E−22MLST81.04550722Up4.72E−234.84E−22CAPNS11.10213996Up5.81E−235.42E−22PARP11.34948718Up6.82E−236.08E−22BAX1.52707657Up7.52E−236.42E−22ATG101.0226174Up1.11E−229.11E−22DIRAS3-2.2142941Down2.04E−221.49E−21HGS1.59409057Up3.22E−222.28E−21TSC21.1673426Up4.51E−223.08E−21CANX1.10228843Up2.00E−211.32E−20WDR45B1.03610977Up4.09E−212.62E−20FOS-2.4192498Down6.68E−214.15E−20DAPK22.15868832Up1.19E−207.20E−20NPC11.40873288Up1.77E−209.80E−20GAPDH1.24841467Up5.02E−202.64E−19PELP11.02773866Up5.19E−202.66E−19HSPA51.07141566Up1.04E−195.18E−19HDAC11.0260015Up6.77E−192.89E−18VMP11.03638398Up1.06E−184.45E−18RUBCN1.05799379Up1.65E−186.76E−18BAK11.61591549Up1.55E−175.67E−17CD461.12531986Up1.83E−176.60E−17CASP31.01180546Up1.98E−176.98E−17DDIT31.50655255Up3.18E−171.10E−16CAPN21.4191918Up3.60E−171.23E−16SQSTM11.73640938Up9.02E−172.89E−16PRKCD1.41302501Up8.94E−162.82E−15CASP81.04011948Up1.52E−154.65E−15FOXO1-1.0757149Down4.28E−151.25E−14SPNS11.31587141Up2.15E−145.96E−14RB1CC11.01420177Up4.73E−141.29E−13ATG16L21.02505335Up2.00E−135.40E−13NRG1-1.280118Down5.44E−121.25E−11DRAM11.2000343Up1.11E−112.54E−11IKBKE1.86880957Up1.30E−112.93E−11TMEM742.33376628Up3.59E−117.83E−11ITGB42.37261675Up1.85E−103.88E−10RGS191.11586497Up5.15E−101.02E−09NRG22.87978625Up1.43E−092.68E−09ATG9B1.05236942Up5.60E−068.32E−06SPHK12.84999763Up8.75E−051.21E−04ITGA32.09689798Up4.48E−045.89E−04TP631.90806327Up2.37E−032.93E−03IRGM2.45185036Up3.98E−034.86E−03NKX2-32.56776198Up4.50E−035.46E−03HSPB82.02223664Up9.57E−031.13E−02logFC: log(fold change); FDR: false discovery rate.Figure 1
The differentially expressed autophagy-related genes (ARGs) between normal and hepatocellular carcinoma (HCC) tissues. (a) The hierarchical clustering of differentially expressed ARG expression levels. (b) The volcano plot for the differentially expressed ARGs. Red indicates high expression, green indicates low expression, and black indicates no difference between HCC and normal tissue. (c) The boxplot for the expression patterns of 62 significant differentially expressed ARGs in HCC and normal tissue.
(a)(b)(c)


## 3.2. Enrichment Analysis of the Differentially Expressed ARGs

To further investigate the biological functions of the differentially expressed ARGs in HCC, we performed GO term function and KEGG pathway enrichment analyses for these genes (Table3). For biological processes (BPs), these ARGs participate in autophagy, processes utilizing autophagic mechanisms, macroautophagy, and the regulation of autophagy. In terms of cellular components (CC), the genes were related to autophagosomes, phagophore assembly sites, vacuolar membranes, and autophagosome membranes. Changes in molecular function (MF) were enriched in protein kinase regulator activity, kinase regulator activity, cysteine-type endopeptidase activity, and BH domain binding (Figure 2, Table 3). Furthermore, KEGG pathway analysis showed that these genes mostly participated in autophagy-animal, apoptosis, platinum drug resistance, and longevity regulating pathways (Figures 3(a) and 3(b); Table 4).Table 3
GO analysis of differentially expressed autophagy-related genes.
CategoryIDTermCountp valueBiological processGO:0006914Autophagy318.54E−33Biological processGO:0061919Process utilizing autophagic mechanism318.54E−33Biological processGO:0016236Macroautophagy232.05E−26Biological processGO:0010506Regulation of autophagy181.38E−17Biological processGO:0016241Regulation of macroautophagy143.91E−16Biological processGO:0070997Neuron death131.53E−10Biological processGO:2001233Regulation of apoptotic signaling pathway137.96E−10Biological processGO:0097193Intrinsic apoptotic signaling pathway112.51E−09Biological processGO:0071900Regulation of protein serine/threonine kinase activity137.65E−09Biological processGO:0042475Odontogenesis of dentin-containing tooth71.68E−08Cellular componentGO:0005776Autophagosome71.26E−08Cellular componentGO:0000407Phagophore assembly site51.38E−08Cellular componentGO:0005774Vacuolar membrane114.53E−08Cellular componentGO:0000421Autophagosome membrane42.83E−06Cellular componentGO:0098589Membrane region81.05E−05Cellular componentGO:0005770Late endosome71.43E−05Cellular componentGO:0005765Lysosomal membrane81.74E−05Cellular componentGO:0098852Lytic vacuole membrane81.74E−05Cellular componentGO:0101031Chaperone complex33.25E−05Cellular componentGO:0045121Membrane raft77.15E−05Molecular functionGO:0019887Protein kinase regulator activity74.93E−07Molecular functionGO:0019207Kinase regulator activity71.14E−06Molecular functionGO:0004197Cysteine-type endopeptidase activity51.55E−06Molecular functionGO:0051400BH domain binding34.96E−06Molecular functionGO:0008234Cysteine-type peptidase activity61.35E−05Molecular functionGO:0051087Chaperone binding51.72E−05Molecular functionGO:0031072Heat shock protein binding53.93E−05Molecular functionGO:0044389Ubiquitin-like protein ligase binding76.46E−05Molecular functionGO:0030295Protein kinase activator activity49.71E−05Molecular functionGO:0019209Kinase activator activity40.000135Figure 2
Gene Ontology (GO) enrichment analysis of the differentially expressed autophagy-related genes. The size of each dot represents the count of genes; the color represents thep.adjust.Figure 3
Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis of differentially expressed autophagy-related genes (ARGs). (a) The outer circle shows a scatter plot for each term of the logFC of the assigned gens. Red circles display upregulation, and blue ones display downregulation. (b) The heat map of the relationship between ARGs and enriched pathways. The color of each block depends on the logFC value.
(a)(b)Table 4
KEGG analysis of differentially expressed autophagy-related genes.
IDTermCountp valuehsa04140Autophagy-animal162.50E−16hsa04210Apoptosis102.02E−08hsa01524Platinum drug resistance75.52E−07hsa04211Longevity regulating pathway72.15E−06hsa04136Autophagy-other52.24E−06hsa04215Apoptosis-multiple species52.24E−06hsa04115p53 signaling pathway68.78E−06hsa04218Cellular senescence81.14E−05hsa04657IL-17 signaling pathway63.84E−05hsa04151PI3K-Akt signaling pathway101.19E−04Table 5
The three ARG-based prognostic models identified by multivariate Cox regression analysis.
Gene symbolCoefHRHR.95LHR.95Hp valuePRKCD0.31751.37371.00551.87670.0461BIRC50.24791.28141.01941.61080.0337ATIC0.43971.55231.00382.40060.0481


## 3.3. Identification of Prognostic-Related ARGs

To identify ARGs associated with the clinical outcomes of patients, univariate Cox regression analysis was applied with the criterion of ap value < 0.05. As shown in Figure 4, 34 genes were selected and significantly associated with the overall survival (OS) of HCC patients. The hazard ratio (HR) of each gene was calculated, and all 34 genes had an HR>1, indicating that all of these genes are risk factors for OS in HCC patients.Figure 4
Thirty-four differentially expressed autophagy-related genes with a prognostic value determined by univariate Cox regression.


## 3.4. Construction and Validation of the ARG-Based Prognostic Model

Based on the prognosis-related ARGs, multivariate Cox regression analysis was performed to construct the ARG-based prognostic model. Three genes, PRKCD, BIRC5, and ATIC, were finally selected to construct the model (Table5), and the risk score formula of the model was as follows: riskscore=0.3175×PRKCDexpression+0.4397×ATICexpression+0.2479×BIRC5expression.Figure 5
Classification of patients into high- and low-risk groups by median risk score. (a) Kaplan-Meier plot of patients in a low- or high-risk group (p=1.049e−04). (b) The number of patients in different risk groups.
(a)(b)To assess the performance of the prognostic model in predicting the clinical outcome of HCC patients, we calculated the risk score of each HCC patient and divided patients into high- or low-risk groups using the median risk score as the cutoff value. The survival curve indicated that the high-risk group (n=117) suffered significantly lower 3- and 5-year survival rates than the low-risk group (n=118) (Figure 5). The risk score distribution, survival status of each patient, and heat map of the three gene expression profiles in the TCGA database are shown in Figure 6. The results indicated that survival time decreased and the mortality rate increased as the risk scores increased. We further analyzed the associations between the expression of the three ARGs and clinical parameters in HCC patients. The results of the independent t-test are shown in Figure 7. We observed significant correlations between overexpression of PRKCD and advanced histological grade (p=0.001), advanced T stage (p=0.011), advanced pathologic stage (p=0.011), and female sex (p=0.032). High BIRC5 expression was closely related to advanced T stage (p=0.007), advanced pathologic stage (p=0.007), and advanced histological grade (p=3.352e−04). Overexpression of ATIC occurred in the advanced T stage (p=0.047), advanced pathologic stage (p=0.047), and advanced histological grade (p=0.002). In particular, we used the same method to analyze the correlations between the risk score calculated by the model and clinical parameters. As shown in Figures 7(k)–7(n), we found that a high-risk score was significantly related to younger patients (p=0.047), advanced T stage (p=0.002), advanced histological grade (p=0.001), and advanced pathologic stage (p=0.002).Figure 6
Prognostic risk score of the differentially expressed autophagy-related genes. (a) The risk score of each hepatocellular carcinoma patient. (b) The patient survival based on the risk score. (c) Heat map of the expression of three genes in high- and low-risk groups.
(a)(b)(c)Figure 7
Correlations between expression of three key genes and clinical parameters (the results ofp<0.05 are shown).
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)(k)(l)(m)(n)Then, we wanted to determine whether the risk scores or other conventional clinical parameters of HCC patients were independent risk factors for the OS of patients. In Figure8, univariate Cox analysis showed that advanced pathologic stage, advanced T stage, advanced M stage, and risk score were risk factors for OS. However, after adjusting for clinical parameters, only the risk score remained an independent prognostic indicator for HCC patients in multivariate analyses (HR=1.745, 95%CI=1.420‐2.145, p<0.001). Moreover, we plotted the ROC curve to compare the predictive value of the risk score with other clinical parameters, as is shown in Figure 9. The results indicate that the T stage (AUC=0.708) and pathological stage (AUC=0.703) have the highest predictive value among the conventional clinical parameters. However, the predictive value of the risk score (AUC=0.758) was higher than that of the T stage and pathological stage.Figure 8
Univariate and multivariate Cox regression analyses of OS in HCC patients. (a) Univariate Cox regression analysis. (b) Multivariate Cox regression analysis. Age: ≤65 vs. >65; gender: male vs. female; grade: G1/G2 vs. G3/G4; stage: I/II vs. III/IV; T stage: T1/T2 vs. T3/T4; M stage: M0 vs. M1; N stage: N0 vs. N1-3; risk score: high-risk score vs. low-risk score (median risk score as the cutoff value).
(a)(b)Figure 9
Receiver operating characteristic (ROC) curve analysis for the prognostic value of the prognostic model and other conventional clinical parameters. AUC: area under the curve.


## 3.5. Confirmation of the Prognostic Model Using an Additional Data Cohort

The International Cancer Genome Consortium (ICGC) database contains large-scale cancer genome studies in tumors from 50 different cancer types [14]. To further validate the predictive value of this model, we downloaded the data cohort (Liver Cancer-RIKEN, JP, https://dcc.icgc.org/releases/current/Projects/LIRI-JP), which contains expression profiles and clinical follow-up information from 232 HCC patients. We analyzed the correlation between the expression levels of three genes and the survival of HCC patients in both the TCGA and ICGC databases. Kaplan-Meier analysis results indicated that high expression of the three genes was significantly associated with inferior OS in HCC patients (Figure 10). Similarly, these patients were divided into low- or high-risk groups based on the risk score calculated by the prognostic model constructed based on the TCGA data. The Kaplan-Meier analysis of the two groups was significantly different (p=4.505e−06<0.001) and a similar trend was observed in TCGA (Figure 11(a)). In addition, the AUC calculated by the ROC curve was 0.739, which indicated that the model has good performance for predicting HCC patient survival (Figure 11(b)).Figure 10
The correlation between the expression level of three genes included in the prognostic model and survival of HCC: (a) ATIC expression level and patient survival (data from TCGA), (b) PRKCD expression level and patient survival (data from TCGA), (c) BIRC5 expression level and patient survival (data from TCGA), (d) ATIC expression level and patient survival (data from ICGC), (e) PRKCD expression level and patient survival (data from ICGC), and (f) BIRC5 expression level and patient survival (data from ICGC).
(a)(b)(c)(d)(e)(f)Figure 11
Kaplan-Meier and ROC curves for the ARG prognostic model in the cohort of the ICGC. (a) The differences between the low-risk (n=116) and high-risk (n=116) groups determined by the log-rank test (p=4.505e−06). (b) ROC curve analysis for survival prediction by the prognostic model.
(a)(b)


## 4. Discussion

The incidence of HCC is characterized by insidiousness, rapid progression, and a low early diagnosis rate. Most patients are already in advanced tumor stages when they receive treatment [15, 16]. Due to the heterogeneity of liver cancer, conventional clinical parameters such as age, sex, grade, and TNM stage often do not accurately predict clinical outcomes [17]. Therefore, there is an urgent need to develop new prognostic features to facilitate the prediction of clinical outcomes in HCC patients. In recent decades, many studies have focused on identifying novel biomarkers to promote the prediction of HCC patient survival [18–20]. Based on the advantages of this type of research, the combination of multiple prognosis-related genes with conventional clinical parameters to construct a prognostic model may have better predictive value for HCC patients.Autophagy is an intracellular self-digestion process that plays a fundamental role in cell homeostasis, and numerous studies have demonstrated that ARGs play important roles in tumorigenesis [21–23]. Therefore, identifying biomarkers from ARGs may provide new perspectives of the diagnosis or treatment for various cancers. Recently, accumulating reports have shown that using ARGs to build a prognostic model can provide better prediction of clinical outcomes for cancer patients. Lin et al. established a prognostic signature based on three ARGs in thyroid cancer [11]. Wang et al. analyzed TCGA and GEO datasets to construct and validated an autophagy-clinical prognostic model in bladder cancer [12]. Although numerous studies have confirmed the close relationship between the development of HCC and autophagy, a prognostic model based on ARGs in HCC had not been reported previously.In this study, we used high-throughput expression profile data from TCGA to construct an ARG prognostic model and validated this model using data from TCGA and ICGC. First, we identified 62 ARGs were differentially expressed in tumor tissues. GO analysis and KEGG pathway enrichment analysis were used to explore the potential biological functions of these genes. GO enrichment revealed that these genes were mostly enriched in autophagy, a process utilizing autophagic mechanisms, macroautophagy, and autophagosomes, indicating that the differentially expressed genes mainly affect the progression of liver cancer by affecting autophagy [24, 25]. KEGG pathway enrichment was mainly enriched in autophagy-animal, apoptosis, and platinum drug resistance. The results are consistent with previous studies showing that dysregulated apoptosis and platinum-based resistance are common features of many cancers, including HCC [26–28]. Then, univariate Cox regression analysis identified 34 genes closely related to the survival of HCC patient, and multivariate Cox regression analysis finally selected three key genes (PRKCD, BIRC5, and ATIC) to construct the prognostic model. We divided patients into high- or low-risk groups and found that the high-risk group had a lower 3- or 5-year survival rate. Then, we confirmed that the risk score calculated by the model formula was an independent prognostic indicator after adjusting for other clinical parameters. In addition, ROC curve analysis demonstrated that the risk score has a better predictive value than other clinical parameters. Furthermore, we also observed a similar trend of survival analysis and ROC curve analysis in an independent dataset from the ICGC database which further confirmed the reliability of this prognostic model.The three ARGs that we have selected to construct the prognostic model have been reported to be involved in the development of cancer in other studies. PRKCD is one of the PKC family members and is a family of serine- and threonine-specific protein kinases that can be activated by calcium and the second messenger diacylglycerol [29]. Previous studies suggested that PRKCD can show completely opposite effects on tumors in different cell types [30]. For example, PRKCD overexpression protected keratinocytes from UV-induced apoptosis and enhanced long-term survival which is a protective mechanism against skin carcinogenesis [31, 32]. LV et al. found that PRKCD can promote tumor progression in pancreatic cancer [33]. In liver cancer, Zhang et al. reported that dihydromyricetin inhibits the migration and invasion of hepatoma cells by reducing MMP-9 expression via a mechanism that is dependent on the upregulation of PRKCD [34]. Nambotin et al. confirmed that Frizzled-7 displays anticancer properties against HCC involving PRKCD activation [35]. In addition, Cai et al. implied that PRKCD is an independent gene involved in the progression of NAFLD to HCC [36]. Combined with our research, we suggest that high expression of PRKCD in HCC may affect tumorigenesis and serve as a biomarker for predicting patient survivals. BIRC5 (also known as survivin) is a member of the inhibitor of the apoptosis (IAP) gene family, which is essential for cell division and can inhibit apoptotic cell death [37, 38]. Ambrosini et al. first described this gene as an oncogene that is prominently expressed in all common human cancers of the lung, colon, pancreas, prostate, and breast [39]. Consistent with the results of previous studies, Chen et al. recently found that increased BIRC5 expressions are associated with histological grade, tumor size, and TNM stage in HCC patients which is consistent with our findings [40]. The protein encoded by ATIC is the last enzyme in the de novo purine biosynthetic pathway [41]. Previous studies have demonstrated that the purine synthesis pathway correlates with cancer cell proliferation [42, 43]. Li et al. recently found that ATIC is an oncogenic gene that promotes survival, proliferation, and migration by targeting AMPK-mTOR-S6K1 signaling in HCC [44]. Although many studies have reported that all three ARGs we selected to construct prognostic models were directly or indirectly involved in various cancers, this is the first study to combine those genes with clinical information to predict the prognosis of HCC.In summary, we constructed and validated a prognostic model based on three ARGs and this model could be a useful tool to predict the survival of HCC. To our knowledge, the three ARG-related prognostic models have not been reported previously, and the differentially expressed ARGs may provide a new perspective for the study of HCC molecular mechanisms. However, there are some limitations of our research that should be taken into consideration. First, we only focused on the mRNA levels of these genes, and protein levels should be further investigated. Second, the results are exclusively based on the TCGA and ICGC datasets, and additional independent cohorts are necessary to confirm the reliability of this model. Third, the results of our study are descriptive and the potential molecular mechanisms of the three genes warrant additional functional experiments.



---

*Source: 5564040-2021-03-12.xml*