A signature of 13 autophagy‑related gene pairs predicts prognosis in hepatocellular carcinoma

ABSTRACT Growing evidences suggest that autophagy plays a momentous part in the tumorigenesis and development of hepatocellular carcinoma (HCC). However, there are not many researches to predict the prognosis of HCC using autophagy-related genes. Therefore, based on the clinical information and RNA-Seq data of The Cancer Genome Atlas data portal (TCGA), 13 autophagy‑related gene pairs were screened to build the autophagy‑related signature to predict the prognosis by least absolute shrinkage and selection operator (LASSO) regression analysis. Besides, the International Cancer Genome Consortium (ICGC) cohort was further applied to verify the autophagy‑related prognostic signature. Gene ontology (GO), Kyoto Encyclopedia of Genes and Genomes (KEGG) and gene set enrichment analysis (GSEA) were also used to predict the relevant function of the autophagy-related gene pairs signature. As shown in the results, the autophagy-related gene pairs were mainly involved in process utilizing autophagic mechanism, autophagy, macroautophagy and cellular response to oxidative stress. The immune cell levels in the high-risk and low-risk group were explored, which showed that the three immune cells were obviously increased in the high-risk group, while the five immune cells were obviously increased in the low-risk group. In conclusion, an autophagy‑related prognostic signature was established to predict the prognosis of HCC patients with great accuracy and we found that autophagy‑related prognostic signature was related to infiltrating immune cells.


Introduction
According to GLOBAOCAN, liver cancer was one of the most deadly cancers around the world in 2018 [1]. The most frequent category of primary malignant liver cancers was hepatocellular carcinoma (HCC) [2]. Because of the high mortality of HCC, researching effective and stable biomarkers of HCC was a significant task [3]. Nevertheless, the molecular pathological mechanisms of the development of HCC are still unclear. The accurate prediction of HCC prognosis could contribute to individual treatment and reduce waste of resources. At present, the assessments of prognosis mainly include the tumor stage and the degree of liver function damage [4]. These methods could not accurately reflect the prognosis of HCC. Hence, it is urgent to develop and validate a novel prognostic assessment system.
Many researches have proposed signatures based on gene expression for outcome stratification in HCC patients. 11 genes, including TNRC6B, VAV3 andRNF213, with a range of mutational prevalence from 1% to 3%, could predict the prognosis of HCC patients [5]. Jing-Xian Gu et al. have established a signature of six-lncRNA to predict the prognosis of HCC [6]. The 4-gene signature, which was strongly correlated with the age, grade, and T stage of HCC patients, acted as a candidate prognostic factor [7]. However, the above signatures have not been widely used in clinical applications because of the lack of abundant validation and overfitting on small development datasets [8]. The normalization of gene expression level was insufficient in previous studies, which posed a significant challenge of the data analysis with high heterogeneity among diverse datasets and sequencing platforms differences [9]. Recently, a novel way based on corresponding rank relative to the gene expression values has been put forward to conquer the defects of the gene expression levels normalization and scaling, and has obtained reliable consequences in many researches [10,11].
The word autophagy date from the Greek word, denoting 'self-eating', and refers to the collection of various processes that deliver cellular materials to lysosomes for digestion. In mammals, microautophagy, macro-autophagy [12], and chaperonemediated autophagy [13] are currently known. Three main types of autophagy, of which macroautophagy is the most studied form. Because the impact of autophagy was complex in cancer, many studies have looked at the relationship between autophagy and HCC. For example, ATG14, an autophagy-related protein, was abnormally increased in HCC, which is related with the poor prognosis of patients with HCC [14]. In addition, autophagy was connected with pattern recognition receptor (PRR), cell death pathways and inflammatory in tumor cells, and therefore changed the immunogenicity of the tumor microenvironment (TME) and the anti-tumor immune response development [15]. Nevertheless, the prognosis effect and the relationship with immune cells of autophagy have not been clarified in HCC.
In the study, we established the signature of autophagy-related gene pairs through analyzing the gene expression profile of HCC. the prognosis of HCC patients could be predicted using the signature. Moreover, the association of signature of autophagy-related gene pairs and tumor microenvironment was further analyzed.

Acquisition of HCC datasets and autophagy-related genes
The RNA-seq expression matrix and corresponding clinic information of 374 HCC specimens were gained from The Cancer Genome Atlas (TCGA) database (https://portal.gdc.cancer.gov) [16], which was acted as the discovery dataset. The International Cancer Genome Consortium (ICGC) database (https://icgc.org/) including 243 samples was served as the validation dataset [17]. 232 autophagy-associated genes were further obtained from the Human Autophagy Database (HADb, http://autophagy.lu/clustering/index. html) [18]. The expression data of autophagyassociated genes were extracted from the discovery dataset and the validation dataset. Subsequently, the autophagy-related genes with uneven distribution (determined by median absolute deviation < 0.5) and relatively small change were removed [19]. Consequently, we obtained 127 autophagyassociated genes, which were shared between the two cohorts.

Establishment of autophagy-related gene pairs signature
The score for each autophagy-related gene pair was computed by the paired comparison of the gene expression levels in each patient. Based on the suggested method, if the expression level of preceding gene was much higher than the latter gene expression level in a specific autophagyrelated gene pair, the score of this autophagyrelated gene pair was 1, otherwise, the score was 0. The common autophagy-related gene pairs in two cohorts were screened out as the candidate gene pairs to reflect the prognosis of HCC. The Kaplan-Meier analysis coupled with the log-rank test was applied for the univariate analysis, which was used to select autophagy-related gene pairs related with prognosis. The least absolute shrinkage and selection operator (LASSO) analysis was applied to establish the final model [20]. The autophagy-related gene pairs signature was computed as following: Prognostic score = (Exprgenepair-1 × Coefgenepair-1) + (Exprgenepair-2 × Coefgenepair-2) + . . . + (Exprgenepair-n × Coefgenepair-n) [21]. Exprgenepair was the gene pair expression value and Coefgenepair was the LASSO regression coefficient. The cut off of the autophagy-related gene pairs signature was determined by time-dependent receiver operating characteristic (ROC) curve for the overall survival (OS) in the TCGA cohort [22].

Validation of autophagy-related gene pairs signature validation
In order to validate the prognostic signature, the prognostic score was calculated by the autophagyrelated gene pairs signature in the ICGC cohort. Furthermore, the HCC patients were also separated into high risk and low risk group based on the identical cut-off value in the ICGC cohort. The log-rank test was performed to compared survival curves between two groups. Meanwhile, patients with different tumor stage and grade were also performed Kaplan-Meier analysis. The ROC curve analysis was applied to demonstrate the specificity and sensitivity and of autophagy-relevant gene pairs signature in outcome prediction. Univariate/multivariate cox-regression models were applied to estimate risk elements that impacted OS.

Functional analysis of autophagy-related gene pairs signature
To explore the potential functions of the autophagyrelated gene pairs signature, GO analysis and KEGG analysis were performed by the 'clusterProfiler' package. Biological processes were collected from the Gene Set Enrichment Analysis (MSigDB C2 and C5 databases) (https://www.gsea-msigdb.org/gsea/ index.jsp) [23]. The false discovery rate (FDR) adjusted p-value below 0.05 was considered significant.

Profiling of immune cells infiltration
CIBERSORT algorithm, a novel computational method, uses 547 reference gene expression values for estimating enrichment of various immunocyte compositions [24]. In our research, CIBERSORT was performed to detect the abundance of 22 infiltrating immune cells in the high risk and low risk group. According to the Monte-carlo sampled, the deconvolution p-values of all patients were computed to offer reliability in the assessment. The expression profiling of gene of TCGA was put on the online platform CIBERSORT website (http://cibersort.stanford.edu/) for analysis with the default signature matrix at 1000 permutations.

The analysis of correlation
The heat map of correlation was performed by the 'corplot' package, which was aimed to explore the correlation of autophagy-related gene pairs. The relationships between infiltrating immune cells and risk score of patients were also performed.

Statistical analysis
The statistical analyses were processed by R version 3.6.0 software (https://www.r-proje ct. org/) and Perl version 5.32.0.1 software (https:// www.perl.org/). The differences among the two groups were compared with student's t-tests. The 'survival' package was used to perform the survival analysis. The ROC curves were performed by 'survivalROC' package. For all studies, P-value < 0.05 was served as statistically significantly different. Statistical significance was displayed as *P < 0.05, **P < 0.01, ***P < 0.001.

Results
In this study, we aimed to establish autophagy-related gene pairs signature for predicting the prognosis in patients with HCC. We filtered out autophagy-related gene to make the gene pairs and these autophagy-related gene pairs were applied to construct a prognostic signature. Function enrichment analysis was applied to analyze the potential functions of autophagy-related gene pairs. In addition, according to the autophagy-related gene pairs signature, we further explored the infiltration levels of 22 immune cells in the high-risk and lowrisk group to explore the correlation between and immune cells and autophagy-related gene pairs signature.

Validation and assessment of the autophagy-related gene pairs signature
The risk scores of all HCC patients were calculated by the above signature in the ICGC dataset. Subsequently, based on the cut-off value of −0.05, all patients were assigned into two groups: low-risk and high-risk and survival rate were also calculated. Consistent with the result acquired from the TCGA dataset, the low-risk patients had the better outcome than the high-risk patients (P < 0.001; HR = 2.3 (1.5, 3.5)) ( Figure 2d-Figure 2f). In addition, we validated prognostic effect of the autophagy-relevant gene pairs signature in diverse tumor stage and grade. The result showed that the autophagy-related gene pairs signature could significantly predict survival outcome in stage I (P < 0.001; HR = 3.7 (2.2, 6.1)), stage II (P < 0.001; HR = 4.3 (2.2, 8.5)), stage 2162 (P < 0.001; HR = 3.2 (2.0, 5.4)), grade I (P < 0.001; HR = 4.6 (2.2, 9.7)), grade II (P < 0.001; HR = 3.3 (2.2, 5.1)) and grade III (P < 0.001; HR = 4.1 (2.5, 6.7)) ( Figure S3B-S3G).In the ICGC dataset, univariate Cox analysis indicated that the clinical information including gender and tumor stage showed a prognostic ability and multivariate Cox analysis demonstrated that the autophagy-related gene pairs signature remained an independent prognosis factor ( Figure  3b,Figure 3d).

Functional assessment of the autophagy-related gene pairs signature
To assess the potential function of the autophagyrelated gene pairs signature, we performed the GO, KEGG, and GSEA analysis. The GO analysis consisting of biological processes (BP), molecular functions (MF) and cellular components (CC), and suggested that the autophagy-related gene pairs signature genes were mostly enriched in autophagy, transcription factor complex and ubiquitin-like protein ligase binding (Figure 4a). The KEGG enrichment analyses showed that signature genes only enriched in One carbon pool by folate (Figure 4b). GSEA analysis was carried out to analyze the pathways which were obviously altered between the high-risk group and low-risk group.  Figure S4).

The immune cells infiltration between two groups
According to the CIBERSORT, we estimated the relative proportion of the 22 immune cells for each HCC patient within the high-risk class and the low-risk class. The comparison of different immune cells in the low-risk class and high-risk class using CIBERSORT was displayed in Figure 5. B cells memory, T cells follicular helper and Macrophages M0 were obviously increased in the high-risk class than the low-risk class; however, the expression levels of B cells naive, Macrophages M1, T cells CD4 memory resting, Monocytes and Mast cells resting were obviously lower in the high-risk class (Figure. S5). Moreover, the potential correlations between the autophagyrelated gene pairs signature score and the eight immune cells were calculated. Our results demonstrated that the autophagy-related gene pairs signature was remarkably positively correlated to the T cells follicular helper, B cells memory and Macrophages M0, and negatively correlated to the B cells naive, Macrophages M1, Mast cells resting and T cells CD4 memory resting. However, the autophagy-related gene pairs signature score had no association with the Monocytes (Figure. S6).

Discussion
Despite of significant advances in the prevention, diagnosis, and therapy of HCC, it remained the primary reason of cancer death around the world [25][26][27]. According to exome sequencing, genomic characterization studies and transcriptome analysis, HCC had various molecular characteristics and clinical outcomes [28,29]. Because of the heterogeneity of cancer, patients of HCC should select different therapeutic approaches to achieve more significant therapeutic effects. Autophagy played different roles in different stages of HCC and markedly affected the efficacy of multiple therapeutic methods [30]. The analysis of autophagyrelated proteins and pathways in HCC could be helpful for finding biomarkers of specific treatments and predicting the prognosis more accurately. In HCC, the molecular mechanism and clinical significance of autophagy remain unclear. In this research, we obtained a novel prognosis signature from 13 autophagy-related gene pairs of HCC in TCGA cohort and verified in ICGC cohort. The autophagy signature could effectively separate patients into different subgroups according to their survival outcomes. Based on multivariate analysis and Kaplan-Meier method, autophagy-related gene pairs signature could act as the independent predictive factor of HCC patients. The above results were successfully verified in the ICGC cohort. In a word, autophagy-related genes may be closely associated with prognosis of HCC patients.
The traditional prognosis signatures required to preprocess the gene expression levels, which was the vital reason to restrain the wide use of models. The autophagy-related gene pairs signature was obtained through the relative sequencing of the gene expression levels, and paired comparison of weights was only carried out within the patient gene expression profiles, which was a huge advantage. Our prognosis signature could not only eliminate batch effects among different platforms but also did not need the normalization and scaling of data. The formula for the risk score and truncation value could be used in various datasets, which made it easier to be applied in clinical. This method has been performed robustly in cancer-relevant studies [31], and it was the primary superiority in our research.
Autophagy was associated with a wide range of physiological processes and human diseases, and it has been linked to both tumorigenesis and oncotherapy in recent years. Hence, it was important to build the autophagy-related signature for predict the survival outcome for HCC patients. Our 13 autophagy-related gene pairs signature contained 26 genes, which were mainly participated in the functions of autophagy and tumor occurrence and development. For example, the Tp53, as the gene of autophagy-related prognosis signature, could contribute to the chemoresistance and cell survival in HCC under conditions of nutritional deficiency by modulating the activation of autophagy [32]. In addition, the study has demonstrated that HDAC6 could regulate the tumor progression by suppressing the signaling pathway of canonical Wnt/β-catenin in HCC [33]. More and more evidences have shown that the tumor microenvironment (TME) was intricately related to cancer initiation, promotion, and progression. And autophagy and immunity took great participation in TME. Some studies have demonstrated that autophagy played a momentous role in the immunity regulation. Autophagy regulated cellular differentiation, development and homeostasis of multiple immune cells like monocytes, T-cells, and B-cells, thus was deeply involved in the innate immune system as well as lymphocyte activation and survival [34]. For instance, the Bcl-2/Bax ratio was significantly associated with both the absolute number and the rate of CD4(+) T cells and natural killer cell activity in head and neck cancer patients [35]. CX3CL1, served as the autophagy-related gene, was significant associated with NK cellmediated cytotoxicity, T cell activation, and leukocyte migration in lung adenocarcinoma, and a study demonstrated that the differential expression of CX3CL1 could affect the migration activity of CD8(+) T cells and CD56(+) NK cells in hepatoma cells, which had hepatitis B virus replication [36]. All of the above-mentioned autophagyrelated genes, the changes of immune cells, and related pathways and cytokines stated clearly that autophagy, immunity and cancer were all linked.
Nevertheless, there are limitations in this research. Firstly, we adopted a retrospective method to analyze the data, and to validate the results, a prospective cohort needs to be done. Secondly, more details about the biological functions and molecular mechanisms of 24 autophagy-related genes in HCC patients required further assessment. RT-PCR could be used for further clinical applications.

Conclusion
In summary, we established the autophagy-related gene pairs signature with data from the TCGA, ICGC and HADb datasets, which could accurately predict the HCC patient prognosis. The results demonstrated that the low-risk group had a significantly better OS than the high-risk. Moreover, we discovered that the immune cell infiltration levels were obviously different between highrisk group and low-risk group, which indicated that the autophagy-related gene pairs signature was relevant to infiltrated immune cells. The autophagyrelated gene pairs signature will act as a crucially predictive marker to offer insights into the therapeutic strategies.

Highlights
(1) Signature of 13 autophagy-related gene pairs could predict the HCC prognosis. (2) Autophagy-related gene pairs were involved in autophagy and oxidative stress. (3) The immune cell levels in the high-risk and low-risk group were different.