Identification of a novel autophagy signature for predicting survival in patients with lung adenocarcinoma

The flowchart and data acquisition

The workflow of this study is shown in Fig. 1. We downloaded LUAD mRNA sequencing data (level 3) and their corresponding clinical patient data from TGGA (https://cancergenome.nih.gov/). Only the samples with complete clinicopathological information and more than 30 days of overall survival (OS) were included in this study. We used GSE50081 from the GEO database as the external validation dataset. Overall, 490 patients were randomly assigned into a training cohort (n = 245) and a testing cohort (n = 245) to satisfy the following criteria: (1) samples were randomly divided into training and testing datasets; and (2) gender, age, and clinical stage distributions between the two groups looked similar (Table 1). In addition, another 127 samples from GSE50081 were used as the validation dataset. This study was conducted in accordance with TCGA publication guidelines (http://cancergenome.nih.gov/publications/publicationguidelines).

ARG set curation

ARGs were curated from the Human Autophagy Database (HADb, http://www.autophagy.lu/index.html), REACTOME AUTOPHAGY in the Molecular Signatures Database v6.2 (MSigDB, http://software.broadinstitute.org/gsea/msigdb), and genes with “autophagy” term relevance scores > =7 from the GeneCards website (https://www.genecards.org/). After eliminating the overlapping genes, these three gene sets were combined and integrated into an autophagy-associated gene set. The ARG list was comprised of 366 genes when finally constructed.

ARG differential expression analysis

We downloaded all the genes of the LUAD samples from TCGA database in FPKM format, which we then converted to TPM format using the formula: TPM = 10⁶ ×FPKM/sum(FPKM). All expression profiles were converted to [log2(TPM + 1)]. We used the “Limma” package (Ritchie et al., 2015) in R software to identify the differentially expressed genes (DEGs) between tumor and normal tissue. The Benjamini–Hochberg method was used to adjust p values, and we considered adjusted P < 0.05 and fold change (FC) > 1.5 as the cutoff criterion for DEG identification. The intersection of DEGs and ARGs was considered the set of significant DEARGs for further analysis. Additionally, we performed volcano plot and heatmap analysis to screen the common DEARGs across the datasets.

Kyoto Encyclopedia of Genes and Genomes and gene ontology analysis

We analyzed the function of significant DEARGs using the Kyoto Encyclopedia of Genes and Genomes (KEGG), gene ontology (GO) functional enrichment analyses (Ashburner et al., 2000), and the Database for Annotation, Visualization, and Integrated Discovery (DAVID, http://david.ncifcrf.gov/; Jiao et al., 2012). A P value of <0.05 was considered statistically significant.

Prognostic model construction and performance assessment

We first conducted univariate Cox proportional hazard regression to identify the DEARGs that were significantly associated with overall survival (P value < 0.05) in the training cohort using the survival package (http://bioconductor.org/packages/survival/) in R. The LASSO Cox regression method (Sauerbrei, Royston & Binder, 2007) was then employed to select optimal gene combination variables, and then the “glmnet” package in R was used to construct the risk signature. Only genes with non-zero coefficients in the LASSO model were put into the multivariate Cox regression model to calculate the risk score. The prognostic model formula was as follows: Risk score = (expr gene₁ ×β₁) + (expr gene₂ ×β₂) + ⋯ + (expr gene_n × β_n), where “expr” represents the gene_i expression value and “ β_i” represents the estimated regression gene_i coefficient.

Using the median risk score value as a cutoff, we divided the LUAD patients into low-risk and high-risk groups. We employed the Kaplan–Meier (K-M) survival curves to show the OS differences between the high-risk and low-risk groups. We used the area under the curve (AUC) of the time-dependent receiver operating characteristic (ROC) curve and the R package survivalROC to evaluate the risk signature’s efficiency. These analyses were conducted in the TCGA testing dataset and GEO datasets.

GSEA

GSEA (version 3.0, http://www.broadinstitute.org/gsea/index.jsp) was used to evaluate the biological pathways or gene sets that differed significantly between the high-risk and low-risk groups. The parameters were as follows: max gene set size of 500, min size of 15, number of permutations of 1,000, and enriched gene sets with a nominal P value < 0.05 were considered significant.

Immune cell analysis

CIBERSORT (Newman et al., 2015), a deconvolution algorithm, was employed to estimate the relative abundance of immune-infiltrating cell composition in tissues based on their expression profiles. We submitted the LUAD gene expression dataset to the CIBERSORT website (http://cibersort.stanford.edu/) and used LM22 (22 immune cell types) as the signature gene file. The program was implemented with 1,000 permutations. Next, the output values generated by CIBERSORT were defined as immune cell infiltration fractions per sample. The output results were used to compare immune cell infiltration fractions across low-risk and high-risk patients.

Nomogram construction and validation

A prognostic nomogram was constructed to combine ARG signatures and other clinicopathological factors using the “rms” package (https://cran.r-project.org/web/packages/rms/index.html) in R. To evaluate the accuracy of the nomogram, we applied calibration curves and K-M analysis to compare the concordance between the predicted survival and observed survival.

Statistical analysis

All statistical analyses in our study were conducted using R language (version 3.5.1, https://www.r-project.org/). Boxplots and violin plots were generated using the “ggplot2” package in R language. P < 0.05 was considered significant.

Identification of autophagy-related risk signature in the LUAD training cohort

We analyzed the expression of 366 ARGs in 490 LUAD and 59 normal lung tissue samples using the “limma” package in R software. In total, we identified 83 DEARGs from the LUAD samples, including 33 upregulated DEARGs and 50 downregulated DEARGs with a cutoff criteria of adjusted P < 0.05, |FC| > 1.5. The volcano plot and heatmap for these 83 DEARGs in normal and tumor tissues are displayed in Fig. 2. Additionally, to better explore the biological interpretation of these DEARGs, we performed functional enrichment pathway analyses. According to the GO enrichment analysis results, we found that these DEARGs were primarily involved in autophagy, protein binding, and autophagosome related to biological process (BP), molecular function (MF), and cellular component (CC) terms (Figs. 3A–3C). Moreover, KEGG enrichment analysis also indicated that these genes primarily participated in cancer pathways, protein processing in the endoplasmic reticulum, as well as the TNF signaling pathway (Fig. 3D).

Identifying prognostic risk DEARGs in the LUAD training set

Next, we performed univariate Cox proportional hazards regression analysis to identify prognostic DEARGs in the LUAD training set using the coxph function of the survival package in R (Zhang, 2002). We found a total of 20 DEARGs that were significantly associated with OS (P <  = 0.05, Table 2). Moreover, LASSO Cox regression was subsequently used to avoid overfitting problems in the risk signature. Ten key autophagy-related genes (BAK1, DAPK2, ERO1A, GAPDH, IL1B, ITGA6, NLRC4, NUPR1, SERPINA1, and TOMM40) were retained when the optimal lambda value was achieved (Figs. 4A–4B). Finally, an autophagy-related signature was established using multivariate Cox regression and the following risk score formula for each patient was as follows: ${\displaystyle \text{Risk score}=0.49915\times \left(\text{expression value of BAK1}\right)+\left(-0.10340\right)\times \left(\text{expression value of DAPK2}\right)+0.36101\times \left(\text{expression value of ERO1A}\right)+0.15103\times \left(\text{expression value of GAPDH}\right)+\left(-0.26998\right)\times \left(\text{expression value of IL1B}\right)+0.02068\times \left(\text{expression value of ITGA6}\right)+\left(-0.08339\right)\times \left(\text{expression value of NLRC4}\right)+\left(-0.19307\right)\times \left(\text{expression value of NUPR1}\right)+\left(-0.05327\right)\times \left(\text{expression value of SERPINA1}\right)+0.02443\times \left(\text{expression value of TOMM40}\right).}$

Risk scores for each patient were calculated, and the patients in the training set were divided into high-risk (n = 122) and low-risk groups (n = 123), according to the median risk score cutoff. In the training set, we determined the risk score distribution, OS status, and the corresponding expression profiles of 10 ARGs (Figs. 5A–5C). The heatmap showed that patients in the high-risk group tended to have higher expression patterns of risky ARGs (ERO1A, TOMM40, GAPDH, ITGA6, and BAK1). On the other hand, patients in the low-risk group tended to have higher expression patterns of protective autophagy genes (NLRC4, IL1B, DAPK2, SERPINA1, and NUPR1) (Fig. 5C). Moreover, the K-M survival curve and the log-rank test exhibited that patients in the high-risk group had a significantly shorter OS time than those in the low-risk group (median time = 2.15 years vs. 2.91 years, respectively, p < 0.001; Fig. 5D). Additionally, we evaluated the predictive performance of the risk signature model with these prognostic biomarkers using time-dependent ROC curves. The area under the ROC curves for one, three, and five-year OS predictions of the risk scores were 0.705, 0.715, and 0.778 (Fig. 5E), respectively.

Validation of the autophagy signature in TCGA and GEO datasets

To confirm our findings, we performed additional testing and used external validation datasets to assess the predictive performance of the 10-gene autophagy signature. First, we validated our autophagy-related signature using a TCGA testing set as an internal validation series. A total of 245 LUAD samples were collected and used to assess the risk signature’s performance. Using the same risk score cutoff, we classified the patients into high-risk (n = 126) and low-risk (n = 119) groups in the internal testing set. In accordance with our previous findings, the distribution of risk score, OS status, and ARG expression were similar as those in the training set (Figs. 6A–6C). Moreover, patients with higher risk scores had significantly shorter median OS than those with lower risk scores (log-rank test P < 0.001; Fig. 6D). The AUCs for one, three, and five-year OS predictions for the risk scores were 0.747, 0.739, and 0.634, respectively, which results were similar to the training set (Fig. 6E).

We further validated our autophagy signature using another independent data set obtained from GSE50081. The distribution of risk score, OS status, and ARG expression in the testing dataset are shown in Figs. 7A–7C. The results confirmed our model’s ability to predict survival. The 10-autophagy-related signature model could effectively predict the OS in patients from the GSE50081 dataset (log-rank test P =  < 0.0001; Fig. 7D). The AUC values for the one, three, and five-year OS models were 0.732, 0.736, and 0.762, respectively (Fig. 7E). These results confirmed that the autophagy-related signature could accurately predict the OS of LUAD patients.

GSEA of high-risk and low-risk LUAD patient characteristics

We carried out GSEA to explore the high-risk and low-risk groups’ biological processes and signaling pathways associated with the autophagy signature. We compared the gene expression profiles of high-risk and low-risk LUAD patients that were classified by the 10-autophagy-related gene signature in both the training set and testing set. The GSEA results revealed that the genes in the low-risk group were closely associated with several metabolism and immune-related pathways, including arachidonic acid metabolism (NES = 1.65, P = 0.014), surfactant metabolism (NES = 1.63, P = 0.028), and the CD22-mediated BCR regulation pathway (NES = 1.58, P = 0.010). Genes in the high risk-group were enriched in several tumor progression pathways, including cell cycle (NES = −2.16, P = 0), spliceosomes (NES = −2.10, P = 0), and nucleotide excision repair (NES =- 2.09, P = 0). The GSEA results are shown in Fig. 8.

Distinct immune phenotype characterization of high-risk and low-risk LUAD patients

To further analyze the association between the ARGs and the tumor immune microenvironment, we used CIBERSORT software to estimate the infiltration fraction across the 22 distinct immune cell types in LUAD patients. The distribution of the 22 immune cell types in each individual are shown in Fig. 9A. The relative proportions of the 22 immune cell types were found to be weakly to moderately correlated (Fig. 9B). Additionally, we intensively investigated the potential differences between the low-risk and high-risk groups. In the high-risk group, we observed that the relative fraction of M0- and M1- macrophages and T cell CD4+ memory activated were significantly increased, while the relative fraction of Myeloid dendritic cells, Mast cells activated, and T cell CD4+ memory resting were significantly decreased (Fig. 9C). However, we also found no significant differences in CD8 T cell infiltration between these two groups (Fig. 9C). To investigate the immune status of the LUAD tumors, we selected immune checkpoints (PD-L1) and other immune-related genes (including CD4, CD47, CD244, CSF1R, and IL1RN) to explore the differences between the high-risk and low-risk groups. Compared to the high-risk group, we found that CD4, CD244, PD-L1, CSF1R, and CD47 were significantly overexpressed in low-risk patients (P < 0.05, Figs. 10A–10F). Moreover, classic immune checkpoints such as PDCD1(PD-1), CTLA-4, HAVCR2(Tim-3), LAG3, and TIGIT were also compared between the high-risk and low-risk groups, but no significant differences were found (Fig. S1). Taken together, these data implied that the autophagy signature may serve as an indicator of LUAD immune status.

Nomogram construction and validation

A prognostic nomogram can quantitatively predict an individual’s risk by integrating autophagy signature risk scores and clinicopathologic features. The training and testing dataset, both from TCGA, were combined to construct nomograms for validation. We constructed a nomogram to predict OS by incorporating the risk scores with age, gender, and tumor, node, metastasis (TNM) stage. Each variable was assigned points in proportion to its risk contribution to survival, and the C-index to evaluate the OS of the model was 0.721(Fig. 11A). The calibration curves suggested agreement between the actual and predicted OS (Figs. 11B–11D).