Autophagy-related signature as Indicators for the Prognosis of Hepatocellular carcinoma

Hepatocellular carcinoma (HCC) is the most common and deadly type of liver cancer. Autophagy is the process of transporting damaged or aging cellular components into lysosomes for digestion and degradation. There is an accumulative evidence implies that autophagy is a key factor of the progression of cancer. The aim of this study was to determine a panel of a novel autophagy-related prognostic marker for liver cancer. We conducted a comprehensive analysis of ARGs expression proles and corresponding clinical information based on The Cancer Genome Atlas (TCGA) and International Cancer Genome Consortium (ICGC) database. The univariate Cox proportional regression model was used to screen candidate autophagy-related prognostic genes. In addition, the multivariate Cox proportional regression model were helped to prove ve key prognostic autophagy-related genes (ATIC, BAX, BIRC5, CAPNS1 and FKBP1A), which were used to construct prognostic signature.


Abstract
Background Hepatocellular carcinoma (HCC) is the most common and deadly type of liver cancer. Autophagy is the process of transporting damaged or aging cellular components into lysosomes for digestion and degradation. There is an accumulative evidence implies that autophagy is a key factor of the progression of cancer. The aim of this study was to determine a panel of a novel autophagy-related prognostic marker for liver cancer.

Methods
We conducted a comprehensive analysis of ARGs expression pro les and corresponding clinical information based on The Cancer Genome Atlas (TCGA) and International Cancer Genome Consortium (ICGC) database. The univariate Cox proportional regression model was used to screen candidate autophagy-related prognostic genes. In addition, the multivariate Cox proportional regression model were helped to prove ve key prognostic autophagy-related genes (ATIC, BAX, BIRC5, CAPNS1 and FKBP1A), which were used to construct prognostic signature.

Results
Based on the prognostic signature, liver cancer patients were signi cantly divided into high-risk and lowrisk groups in terms of overall survival (OS). Further multivariate Cox regression analysis indicated that the prognostic signature remained as an independent prognostic factor for OS. The prognostic signature in possession of a better Area Under Curves (AUC) has a better performance in predicting the survival of patients with HCC, compared with other clinical parameters.

Conclusion
This study provides a prospective biomarker for monitoring the outcomes in the patients with HCC.

Background:
Liver cancer is one of the most common malignant tumors in alimentary (digestive) system with poor clinical prognosis, and also the second leading cause of cancer-related mortality [1]. Hepatocellular carcinoma (HCC) accounts for 75%-85% of primary liver cancer cases, which make it the most common type of liver cancer 1 . Liver cancer patients were often diagnosed at a late stage so that the e cacy or result of operation was limited and dispirited. Thus, it is crucial to explore an effective therapeutic indicator to achieve better prognosis in clinic. The correlation between autophagy and HCC has been studied, showing a new direction and approach for clinical treatment of liver cancer. For example, it was demonstrated that autophagy can induce resistance of HCC cells to chemotherapeutic agents 2, 3 .
Autophagy is the process of transporting damaged, denatured or aging cellular components into lysosomes for digestion and degradation. Autophagy is an approach or process regulating cell metabolism and maintains homeostasis, and participating in a variety of pathophysiological processes including malignant tumors. However, the e cacy of autophagy in cancer still inconclusive, whereas it has been reported that autophagy acts a dual function in the development of tumors. On the one hand, autophagy degradation of cellular components can provide nutrient supply for tumor cell survival. On the other hand, it can inhibit tumorigenesis by clearing toxic cellular materials. To date, researchers on autophagy and liver cancer mainly focus on cancer progression and chemotherapy resistance 4-6 .
However, the role of autophagy in liver cancer prognosis has rarely been explored.
In this study, we established a signature involving ve autophagy-related genes (ARGs) to predict prognosis in the patients with liver cancer, and found that it was an independent prognostic index for liver patients. Through our study, it is suggested that ARGs play important roles in liver cancer and may show potential and valuable performance for predicting the prognosis of HCC patients.

Data acquisition
A total of 232 ARGs was obtained from the Human Autophagy Database (HADb, http://autophagy.lu/clustering/index.html). The transcription factors (TFs) information, gene-expression pro les and the clinical data of liver cancer cohort were obtained from TCGA portal (https://portal.gdc.cancer.gov/) and ICGC portal (http://dcc.icgc.org/releases/current/Projects) respectively. At the same time, the samples acquired from the TCGA dataset were listed as the training group and those from the ICGC dataset as validation cohort.

Identi cation of differentially expressed ARGs
It was conducted that the differential ARGs expression analysis between HCC and normal liver tissue by using the "limma" package of R. Genes with false discovery rate (FDR) < 0.05 and |log2 fold change (FC)| > 1 was considered to be the signi cantly differentially expressed ARGs. Then Venn diagram was used to detect overlapping target genes.

Differentially expressed ARGs enrichment analysis
Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway genome and gene ontology (GO) genome were selected to conduct a series of gene functional enrichment analysis to nd the major biological attributes. In addition, we employed the GOplot package for visualing the enrichment terms.

Construction of OS risk prognostic model based on ARGs
Initially, the potential prognostic ARGs related to liver cancer patients' OS were selected by univariate Cox regression analyses. Next, using multivariate Cox regression analyses to explore whether these candidate prognostic genes can be independent predictors of prognosis monitoring in TCGA cohorts. Finally, several prognostic ARGs were selected (out or excluded) and Cox proportional hazards regression was employed to build the OS risk prognostic model. The autophagy-related risk score for each patient was calculated by relative expression value multiplied regression coe cients. X-tile analysis was used to select optimal cutoff value in the training group. Simultaneously, the HCC patients were separated into low-risk and high-risk groups according to the value. The overall survival was assessed by using Kaplan-Meier curve.
Then, the differences in survival were estimated into high-or low-lrisk groups based on log-rank test.The Cox regression model was employed to analyze the factors affecting the OS in patients with HCC. The risk score, age, gender, cancer grade, cancer stages were used as covariates.

Statistical analysis
The cut-off point for risk score was determined by X-tile 3.6.1 software, based on the minimal P value. Statistical analysis was performed with R software. Receiver operating characteristic (ROC) curve was used to estimate prognostic value of the risk score. AUC value equal to or greater than 0.7 were regarded a signi cant predictive value. P value < 0.05 was considered as signi cant.
Following differentially expressed TF genes were screened in TCGA and ICGC database respectively, Venn diagram was performed to nd the intersection (Fig. 3a-c). Then a co-expression network was established with these differentially expressed TF genes and ARGs (Fig. 3d). The results showed that differentially expressed ARGs were co-expressed interactively with ERG1, ERG2 and FOXM1, etc.

Functional enrichment evaluation of the differentially expressed ARGs
The functional enrichment analysis of those differentially expressed ARGs was performed by GO and KEGG pathway, which elucidated the biological functions (Figs. 4 and 5). The top enriched GO annotation were related to peptidyl-serine modi cation, regulation of protein binding and regulation of binding during the biological process. Cellular components included pseudopodium, pore complex, and integrin complex. And for molecular function, the terms of BH domain binding, SMAD binding, and chaperone binding were the top three. The KEGG enrichment pathways ware notably associated with Apoptosis, Colorectal cancer, Shigellosis, Platinum drug resistance, Protein processing in endoplasmic reticulum, and so on. Most of the Z-score of enriched pathways were less than zero, suggesting that most of them would be decreased. For ARGs differential express genes, the KEGG pathways including the IL-17 signaling pathway and the AGE − RAGE signaling pathway were displayed by heat-map.

Identi cation of prognostic ARGs
To reveal the relationship between the prognosis of patients and differential ARGs pro les, we evaluated the data obtained from TCGA by univariate Cox regression analysis. 14 ARGs were found signi cantly associated with liver cancer patients' prognosis concurrently (Fig. 6a). Then further multivariate Cox regression analysis were performed (Fig. 6b). Finally, ve genes including ATIC, BIRC5, BAX, CAPNS1 and FKBP1A were identi ed as risk genes in OS and to develop prognostic signature. Moreover, the Kaplan-Meier analysis was performed to focus on the prognostic ability of each ARG. The results from the TCGA database indicated that overexpression of ATIC, BIRC5, BAX, CAPNS1 and FKBP1A was closely related to inferior OS of liver cancer patients (Fig. 7a). And similar results were obtained from the ICGC database.
However, there was no signi cant difference in OS with regard to BIRC5 and CAPNS1 expression (Fig. 7b).

Establishment and de nition of the prognostic risk model
The risk score formulas for OS-based prognosis were constructed as follows: rick score = (0.5554 × expression value of ATIC) + (-0.3096 × expression value of BAX) + (0.2345 × expression value of BIRC5) + (-0.3780 × expression value of CAPNS1) + (0.5091 × expression value of FKBP1A). Subsequently, the X-tile analysis was performed to determine 1.9 as the optimal cutoff points of risk score in the TCGA set ( Fig. 8a-c). Using the cutoff point calculated from the risk score of each patient from the TCGA cohorts, we categorized them into two groups, including a high-risk group and a low-risk group. The risk scores distribution and survival status of patients in TCGA dataset were shown in Fig. 8d. Next, the heat-maps showed the expression pro les of the ve risk ARGs (Fig. 8e). We observed that high-risk patients had a higher level of risky genes (ATIC, BIRC5, BAX, CAPNS1 and FKBP1A). The results from ICGC dataset were similar (Fig. 8g).
The correlations between the risk model and clinical parameters in patients with HCC are summarized in Table 1. Furthermore, the results showed that histological grade (P = 0.029), pathological T stage (P < 0.001) and pathological stage (P < 0.001) were closely correlated with the risk score.
Prognostic risk model is an independent predictor of liver cancer To compare clinical parameters and the independent predictive value of the autophagy-related risk score model, the univariate and multivariate Cox regression models were utilized. In TCGA dataset, the univariate Cox analysis revealed that risk score, tumor stage, T and M stages were correlated with OS of HCC patients (Fig. 9a). And the ve-gene risk score remarkably correlated with survival in multivariate Cox regression analysis, whereas other factors were not signi cant (Fig. 9b). Both univariate and multivariate Cox regression analysis demonstrated that the risk characteristics, tumor gender and stage of patients with liver cancer were markedly correlated with OS in the ICGC dataset (Fig. 9e, f). These results con rmed that the autophagy-related prognostic model could independently predicted liver cancer patient's survival.
The risk score of each patient in the TCGA set was calculated with the ve-gene risk score formula. As expected, the risk score classi ed liver cancer patients into high-and low-risk group based on optimal cutoff point with a signi cantly different prognosis (Fig. 9c). And the high-risk patients had lower OS than the low-risk patients. The ICGC dataset was treated as an independent external validation and was processed in the same way. The Kaplan-Meier analysis demonstrated that the prognostic ability of the signature once again (Fig. 9g).
With the help of the OS-based ROC curves, the predictive performance of the risk score were assessed. For the risk score of OS in TCGA dataset, the AUC values was 0.746, which was apparently higher than that of age (AUC = 0.524), gender (AUC = 0.504), tumor grade (AUC = 0.501), tumor stage (AUC = 0.678), T stage (AUC = 0.679), M stage (AUC = 0.508), and N stage (AUC = 0.508) (Fig. 9d). Afterwards in the IGCG dataset, the AUC of the signature was also as high as 0.744 (Fig. 9h). These data suggested that the risk score was competitive in the survival prediction of liver cancer patients.

Discussion
HCC is one of the most lethal human cancers. With the development of clinical management of liver cancer, some prognostic factors including tumor volume, grade and stage and the number of lesions are characterized 7,8 . But, an effective molecular biomarker for HCC prognosis monitoring is still urgently needed. There is an increasingly emerging evidence suggest that autophagy is closely implicated in tumorigenesis and progression 9 . The exploration of mechanism about autophagy opens up a new prospect approaching for liver cancer. In present, high-throughput biological technologies have been widely applied in the early diagnosis of cancer. Therefore, using large-scale databases will help to explore the expression pattern of ARGs and reveal the prognosis of HCC patients.
In the current study, based on the high-throughput expression data from two datasets (TCGC and ICGC), we aimed to screen key ARGs having intimate connection with the prognosis of patients with HCC. Firstly, 22 ARGs were found having differentially expression between liver tumor and normal tissues, including 3 up-regulated and 19 down-regulated genes. The results of GO terms and KEGG pathways analysis showed the major enrichment in tumor biological process and molecular function. It was observed that the KEGG pathways enriched by the differentially expressed ARGs was involved in several types of cancer. The results suggested that speci c autophagy pattern may played a role in occurrence and progression of liver cancer. Then we observed 14 risk ARG related to OS in the TCGA database by dimensionality reduction in univariate survival analysis. And further multivariate survival analysis decided ve key prognostic ARGs (ATIC, BAX, BIRC5, CAPNS1, FKBP1A) to construct prognostic risk models, which could provide accurate prognosis for patients with liver malignant tumor. Multivariate Cox regression analysis of prognostic risk model and clinical parameters showed that the ve-gene risk score is conducive to independently predict the prognosis of patients with liver cancer. Meanwhile the results of ROC curve analysis suggested that the prognostic risk model could distinguished healthy people from liver cancer patients accurately.
ATIC, a bi-functional protein enzyme, catalyzes the nal two steps of the de novo purine biosynthetic pathway. Recent studies have shown that ATIC is expressed at high levels in lung cancer and related to poor patient prognosis 10 . There are few reports on its function in liver cancer. It has reported that upregulated of ATIC in HCC was correlated with poor survival and it supported propagation of HCC cells by regulating AMPK-mTOR-S6 K1 signature 11 .
Bax has been proved to be one of the most widely-characterized participants in autophagy. It is a member of the BCL2 protein family and functions as an apoptotic activator by forming a heterodimer with BCL2. Bax participates in the mitochondria-initiated intrinsic apoptotic signaling and the extrinsic apoptotic pathway that is triggered by transmembrane death receptors 12,13 . Also, Bax has the effect in the crosstalk between endoplasmic reticulum (ER) signaling pathway and mitochondrial pathways 14 . Considering Bax be important role in apoptosis, it is not unexpected that Bax participates in multiple anticancer drugs that induce apoptosis of cancer cells. SKA3, a component of the spindle and kinetochore-related complexes, in uences cell apoptosis by regulation of BAX/Bcl-2 expression in HCC cells 15 . The modulation of MT1G (a low-molecular weight protein with high a nity for zinc ions) on p53 resulted in upregulation of Bax, which leads HCC cells apoptosis 16 .
BIRC5 belongs to the inhibitor of apoptosis (IAP) family, which can prevent apoptotic cell death. It has been reported that BIRC5 exerted its in uence on HCC cells in promoting proliferation 17 . Recent studies have shown that BIRC5 could be regarded as potential biomarkers for molecular diagnosis as well as therapeutic intervention of HCC and with signi cant in uence on prognosis of HCC patients 18,19 .
CAPNS1, a member of the calpain small subunit family, operate as heterodimers and essential for the calpain activity and function. Several studies have shown that CAPNS1 is associated with biological function to tumorigenesis. Indeed, the expression of CAPNS1 has been detected in various cancers, including hepatocellular carcinoma 20, ovarian carcinoma 21 , colorectal cancer 22 , nasopharyngeal carcinoma 23 , and other tumors. In vitro experiments found that CAPNS1 enhances the growth and metastasis of HCC by activating FAK-Src signaling pathway and MMP2 24 .
FKBP1A is a cis-trans prolyl isomerase that binds to FK506 and rapamycin, and inhibits calcineurin and the activity of mTOR 25 . FKBP1A was validated to have overexpressed in HCC and predicted the poor prognosis 26 . However, the speci c mechanism of FKBP1A in HCC is still unclear, which deserve to be further investigated.
In summary, there is a close relationship between autophagy and the prognosis of liver cancer patients. In addition, this study suggests that risk ARGs be potential candidates for prognostic biomarkers in HCC with the value to guide making decisions for the choice of clinical treatment. Further researches are needed to provide more speci c information about these results. At the same time, based on the discovery of autophagy-related patterns which could impact the prognosis of patients with tumor, our work conducted a risk model to bene t distinguish high-risk liver cancer patients. However, the main limitation of our study is that we used already available data from two public databases, and above results are required to be further investigated in prospective studies.         The autophagy-related risk score analysis of HCC patients. X-tile analysis of the prognostic risk score in TCGA training cohort. (a) The X-tile plots of training group. The black circles highlight the optimal cutoff values. (b) Histogram of the entire cohort divided into low-risk score and high-risk score subgroups according to the cutoff value of 1.9. Blue bars represent the low-risk score group, and grey bars represent the high-risk score group. (c) Kaplan-Meier plot of overall survival (OS) in groups strati ed using the optimal cut-off value of risk score. Blue curves represent the low-risk score group, and grey curves represent the high-risk score group. (d) The risk score distribution and OS of patients in TCGA database.
(e) The heat-map of the ve genes expression pro les in TCGA database. (f) The risk score distribution and OS of patients in ICGC database. (g) The heatmap of the ve genes expression pro les in ICGC database.

Figure 9
The prognostic value of the ve-gene risk score of liver cancer patients. (a, b) The univariate and multivariate Cox regression analyses reveal the independent value of the autophagy-related signature for