Autophagy-related genes prognosis signature as potential predictive markers for immunotherapy in hepatocellular carcinoma

Datasets

The RNA-seq transcriptome profiling and corresponding clinicopathological information were obtained for 377 HCC patients from The Cancer Genome Atlas (TCGA) (https://portal.gdc.cancer.gov/) and 260 HCC patients from International Cancer Genome Consortium (ICGC) (https://icgc.org/). HCC patient protein expression information and normalized RPPA value of TCGA-LIHC-L3 datasets downloaded from The Cancer Proteome Atlas (TCPA) (https://tcpaportal.org/tcpa/index.html). A total of 232 ATGs were obtained from the Human Autophagy Database (http://www.autophagy.lu/index.html).

Bioinformatic analysis

Based on the expression of ATGs in normal tissues and tumor tissues of HCC based on the TCGA dataset, we screened 63 differentially expressed genes (logFC > 1 or logFC < −1, and adjusted P-value < 0.05) for later analysis. The Database for Annotation, Visualization and Integrated Discovery (David, V6.8) online website carries out Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) annotation and visualization on differential genes. Univariant Cox regression analyses were used to determine the prognostic value of ATGs in HCC. A total of 21 ATGs were screened and were closely related to HCC prognosis (P < 0.001). Finally, we selected seven genes from 21 prognosis-related genes to establish a prognostic risk signature according to the LASSO Cox regression algorithm using the “glmnet” and “survival” R-packages (Sauerbrei, Royston & Binder, 2007). The seven genes and corresponding coefficients are established by minimum partial likelihood deviance. The sum of the seven genes and the coefficient product is the risk score for each patient. Based on the median risk score, the TCGA and ICGC HCC patients were divided into high-risk and low-risk groups. Gene set enrichment analysis (GSEA) for high-expression genes in the high-risk group (Subramanian et al., 2005). The CIBERSORT method calculates the infiltration abundance of immune cells using the “e1071,” “BiocManager” and “parallel” R-packages, which calculates the cell composition according to the complex tissue gene expression profile (Newman et al., 2015). The deconvolution approach Tumor IMmune Estimation Resource (TIMER) was used to verify the results (Li et al., 2016). In this study, we only revealed macrophage infiltration results.

Statistical analysis

Unpaired Student’s t-test was used to compare differences in gene expression between the two groups. The test of correlation coefficients was used to compare the data correlation between the two groups. The TCGA and ICGC datasets HCC patients were divided into high and low risk groups according to the median value of the risk score. The chi-square test was used to compare the clinical-pathological features of HCC between the two risk groups. Univariate and multivariate Cox regression analyses were used to evaluate the prognostic value of risk score and clinicopathological features in patients with HCC. Receiver operating characteristics (ROC) curves test the efficiency of risk score and clinicopathological features in predicting the 5 year survival rate of HCC. The Kaplan–Meier curve with the log-rank test was performed to compare the overall survival (OS) of the HCC patients in two groups. All statistical analyses were conducted using R v3.6.0, SPSS 19.0 and Prism 7.

Screening and enrichment analysis of ATGs

We obtained 232 ATGs from Human Autophagy Database. Based on the expression in HCC, 63 differentially expressed genes in the TCGA dataset (logFC > 1 or logFC < −1, adjust P-value < 0.05) (Fig. 1A) were screened out, of which FOS, DIRAS3, NRG1, FOXO1 were down-regulated in HCC. The other genes were up-regulated in expression. GO enrichment analysis showed that these genes were mainly enriched in autophagy, the release of cytochrome c from mitochondria and regulation of protein localization to the membrane. KEGG enrichment analysis showed that these genes were mainly concentrated in the PI3K–AKT signaling pathway, Platinum drug resistance, Hepatitis B (Figs. 1B and 1C).

Prognostic value of ATGs and a risk signature built using seven selected ATGs

To investigate the prognostic role of ATGs in HCC. Univariate Cox regression analysis was performed to explore the expression levels in the TCGA dataset. The results showed that 21 genes are significantly correlated with OS (P < 0.001). All 21 genes are risky genes with Hazard ratio >1 (Fig. 2A). To better predict the clinical pathological features and prognosis of HCC with ATGs, the least absolute shrinkage and selection operator (LASSO) Cox regression algorithm was applied to the 21 prognosis-associated genes in the TCGA dataset, which was used as a training set. Seven ATGs were selected to build the risk signature based on the minimum partial Likelihood Deviance and the coefficients were used to calculate the risk score for both the TCGA and ICGC datasets. At the same time, we also analyzed the prognostic value and clinical pathology for seven ATGs which were used to build the risk signature in the TCGA dataset (Figs. 2B–2H; Figs. S1 and S2). These ATGs affect OS in HCC patients. Based on the median risk value, we divided the TCGA and ICGC dataset into a high-risk group and a low-risk group. The survival time of the high-risk group was significantly lower than that of the low-risk group (P < 0.001) (Figs. 2I and 2J). This indicates that the risk score has an important predictive value for the prognosis of HCC (Fig. S3). Seven ATGs were differentially expressed in the high and low risk group (P < 0.0001) (Fig. S4).

Prognostic risk score showed a correlation with clinicopathological features of HCC

Univariate and multivariate Cox regression analyses were used to observe whether risk signature is an independent prognostic indicator for HCC. Univariate analysis showed that stage and risk score was associated with OS in patients with HCC (P < 0.001, Fig. 3A). We put the clinicopathological features and risk score into multivariate Cox regression. The results showed that the stage and risk score were still closely related to the OS of patients (P < 0.001, Fig. 3A). In the validation of the ICGC dataset, we get the same result. Stage and risk scores are closely related to overall patient survival in univariate and multivariate Cox regression analyses. Besides, we also found that gender is associated with OS in the ICGC dataset (P < 0.05, Fig. 3B). In summary, through the above analysis, we fully determined that the risk score from seven ATGs can independently predict the prognosis of HCC.

The ROC curve showed that the risk score (AUC = 0.755) was a perfect predictor of the 5 year survival rate of HCC patients in TCGA dataset (Fig. 3C), and the predicted effect was higher than grade (AUC = 0.505), stage (AUC = 0.663), age (AUC = 0.535), gender (AUC = 0.513). In the validation of the ICGC dataset, we found that the risk score can also predict the patient’s 5 year survival rate (AUC = 0.698) (Figs. 3C and 3D).

In the front, univariate and multivariate Cox regression analyses showed that stage and risk scores were independent prognostic factors affecting the OS of patients with HCC. Further analysis showed that the risk score was closely related to clinicopathological features of HCC patients. As the stage state of HCC patients increases, the risk score increases accordingly in TCGA and ICGC datasets (Figs. 3E and 3F; Fig. S5). The above results fully prove that the prognostic risk signature is closely related to the malignant degree of tumor. The higher the degree of malignancy, the stronger the autophagy. Finally, we performed a gene set enrichment analysis (GSEA) on the high-risk group to show that immune-related hallmarks were enriched (Figs. 3G–3J). Including PI3K–AKT–mTOR signaling, TGF-BETA signaling, IL6-JAK-STAT3 signaling, IL2-STAT5 signaling. The results of the above enrichment analysis suggest that the prognosis signature of ATGs may be closely related to the immune response.

Prognostic signatures of ATGs are closely related to PI3K–AKT–mTOR signaling

In the previous, we found through GSEA that high-expression genes in the high-risk group were mainly enriched in PI3K–AKT–mTOR signaling (NES = 2.102, Nom P < 0.001, FDR q = 0.004). Therefore, here we analyzed the expression of PI3K–AKT–mTOR signaling-related proteins in different risk groups through the TCPA and TCGA databases, and the results showed that AKT and mTOR protein (Figs. 4A–4F) and mRNA (Figs. 4G–4I) expression were up-regulated, but phosphorylated AKT (AKT–pS473, P = 0.003) (Fig. 4B) and mTOR (mTOR–pS2448, P = 0.0116) (Fig. 4D) expression down-regulated, indicating inhibition of AKT and mTOR protein phosphorylation in the high-risk group, resulting in inhibition of the PI3K–AKT–mTOR signaling pathway. Meanwhile, we found that mTOR target molecule P70S6K1 protein and mRNA expression were up-regulated, while P70S6K–pT389 (P < 0.0001) expression was down-regulated. It is well demonstrated that ATGs promote HCC progression by inhibiting PI3K–AKT–mTOR signaling. Finally, we analyzed the OS of mTOR–pS2448 in different risk groups, and the survival time of high expression of mTOR–pS2448 in the high-risk group was significantly higher than that in the low-expression group (P = 0.045). However, in the low-risk group, the mTOR–pS2448 expression had no significant effect on the prognosis of patients (Figs. 4J and 4K). Also, mTOR mRNA shows the same phenomenon (Fig. S6). The phosphorylation of mTOR is inhibited during high-intensity autophagy, so the survival time of mTOR–pS2448 with low expression is shortened. Combined with the effect of mTOR–pS2448 on the OS of HCC patients, we found that autophagy affects malignant tumor progression and is closely related to PI3K–AKT–mTOR signaling.

Differences of biomarkers predicted by immune cell infiltration and immunotherapy among different risk groups

To explore the differences in immune cell infiltration between different risk groups, we used the CIBERSORT and TIMER methods to assess the abundance of each sample based on the TCGA dataset. The results showed that the infiltration level of macrophages in the high-risk group was higher than that in the low-risk group (Figs. 5A–5D). Using CIBERSORT for the ICGC dataset we come to the same conclusion (Fig. S7). At the same time, CD11B mRNA and F4/80 mRNA expression were found to be significantly different in different risk groups based on the TCGA dataset. CD11B and F4/80 mRNA are characteristic molecular markers on the surface of macrophages, which can indirectly predict macrophage abundance (Maruyama et al., 2005; Schledzewski et al., 2006).

PD-L1 and CTLA-4 have been shown to be targeted for a variety of tumor immunotherapy (Chae et al., 2018; Ott, Hodi & Robert, 2013), In HCC, Phase I clinical trials have shown that PD-L1 and CTLA-4 blockers have achieved good clinical results (Kudo, 2017). In this study, we observed significant differences in PD-L1 and CTLA-4 mRNA expression between different risk groups (Figs. 5E, 5F, 5H and 5I). In the high-risk group, PD-L1 and CTLA-4 mRNA expression were up-regulated in the TCGA and ICGC datasets. At the same time, we also noticed that the PD-L1 protein was also significantly elevated in the high-risk group based on the TCPA database (Fig. 5G). In addition, PD-L1 and CTLA-4 mRNA expression levels showed a significant positive correlation with the risk score in the TCGA dataset (Figs. 5J and 5K). The above results indicate that the risk signature can be used as a predictive marker for immunotherapy response.