Development and validation of a novel survival model for acute myeloid leukemia based on autophagy-related genes

Database

The TCGA-LAML dataset (n = 200) was obtained from the TGCA database (https://portal.gdc.cancer.gov/). After deleting data with imperfect clinical information, we included the remaining 140 patients in the study. The GSE37642 dataset was obtained from the GEO database (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE37642), and we specifically used the two datasets GSE37642- GPL96 and GSE37642- GPL570. After merging (n = 562), we used “sva” R package to eliminate any batch effects (Varma, 2020; Leek & Storey, 2007; Leek et al., 2012). The TCGA-LAML cohorts were the training group, the GSE37642 cohorts were the verification group. The autophagy gene set (Table S1) was obtained from the autophagy database (http://www.autophagy.lu/).

Autophagy signature construction and validation

Autophagy-related genes were extracted from TCGA-LAML, and univariate Cox analysis was used, with p < 0.05 considered significant. Next, we performed LASSO analysis and multivariate Cox to obtain the most critical prognostic genes, and then construct an autophagy model. The LASSO coefficients (β) as follows:

Risk Score = (βmRNA1 ×expression level of mRNA1) + (βmRNA2 ×expression level of mRNA2) + ⋯ + (βmRNAn ×expression level of mRNAn) (Livingston et al., 2016; Apfel et al., 1999; Toulopoulou et al., 2019).

The β in this formula refers to the regression coefficient. The GSE37642 data set was used as a validation 1 cohort. In addition, we further verified the reliability of the prognostic gene signature by randomly dividing the training set (TCGA-LAML) into a verification 2 cohort and a verification 3 cohort. The autophagy risk score of each patient was calculated according to the uniform formula determined in the training cohort. We determine the best autophagy risk scoring standard through the “survminer” software package (Walter, Sánchez-Cabo & Ricote, 2015), and then divide the patients into high- and low-risk groups. In addition, we also constructed a prognostic nomogram.

Estimation of immune cell type fractions

The CIBERSORT algorithm is used to estimate the immune cell types of TCGA data (Alaa et al., 2019; Gentles et al., 2015; Newman et al., 2019; Chen et al., 2018).

Generation of immunescore and stromalscore

The ESTIMATE package (Yoshihara et al., 2013) was used to estimate the ratio of immune-stromal components in each sample in the tumor microenvironment in the form of two kinds of scores: Immune Score, and Stromal Score, which positively correlate with the ratio of immune and stroma, respectively. Meaning the higher the respective score, the larger the ratio of the corresponding component in the tumor microenvironment.

Functional enrichment analysis

The Kyoto Encyclopedia of Genes and Genomes (KEGG) and Gene Ontology (GO) analysis of all differentially expressed genes (DEGs) by R software with p < 0.01 set as the threshold. Gene Set Enrichment Analysis (GSEA software, version 4.0.1) was used to investigate the pathways enriched in the high-risk subgroups. The number of random sample permutations was set at 10.

Statistical analysis

LASSO analysis was performed using the “glmnet” package (Engebretsen & Bohlin, 2019; Blanco et al., 2018). The number of folds used in cross-validation was 10. The Time-dependent receiver operating characteristic (ROC) curve was used to evaluate the predictive performance of 10-gene features. The area under the ROC curve (AUC) was calculated by using the “survivalROC” package (Le et al., 2020; Do & Le, 2020; Li et al., 2021; Le et al., 2021). The decision curve analysis was carried out using the “rmda” software package. The “rms” software package was used for nomogram and calibration diagrams. We use one-way ANOVA to analyze multiple sets of normalized data. All statistical analyses were performed using R software (version 3.5.1) and GraphPad Software (version 7.00). p < 0.05 is considered statistically significant.

Establishing an autophagy-related model and functional enrichment analysis

Thirty-five autophagy genes were related to prognosis in TCGA (Fig. 1A), and LASSO regression analysis narrowed down the list (Figs. 1B, 1C), to include 10 autophagy genes (BAG3, BNIP3, CANX, CDKN2A, DIRAS3, NRG2, PARP1, PRKCD, VAMP3, WDFY3) for prognostic model construction (Fig. 1D).

The GO results indicated that 10 autophagy genes were significantly enriched in the biological process (BP) and cellular components (CC) categories (Fig. 1E), such as positive regulation of protein localization to nucleus, regulation of muscle cell apoptotic process, muscle cell apoptotic process, regulation of protein localization to nucleus, negative regulation of organelle organization, positive regulation of muscle cell apoptotic process, positive regulation of protein import into nucleus, positive regulation of protein import, intrinsic apoptotic signaling pathway in response to oxidative stress, protein localization to nucleus, intrinsic apoptotic signaling pathway, regulation of striated muscle cell apoptotic process, regulation of protein import into nucleus, negative regulation of mitochondrion organization, striated muscle cell apoptotic process, selective autophagy, regulation of protein import, inclusion body, integral component of organelle membrane, intrinsic component of organelle membrane, and nuclear envelope. In addition, it is worth noting that the results of the KEGG analysis did not enrich for obvious pathways.

Evaluation of autophagy risk score

After dividing patients into high-risk and low-risk subgroups, we found an important result that the high-risk group was significantly associated with poor prognosis in the TCGA-LAML cohort (P = 6.975e−09; Fig. 2A). The AUC of the one-, three-, and five-year overall survival (OS) in the TCGA-LAML cohort were 0.819, 0.846, and 0.887, respectively (Fig. 2B). Compared with the other six signatures (Chen et al., 2020), our signature showed a higher C-index (0.7240) and AUCs for one-, three-, and five-year OS predictions (Figs. 2C, 2D).

In order to verify the predictive value of the 10-gene signature, we calculated the risk scores of patients in the GSE37642 cohort (validation 1 set). We found that the results of the GSE37642 cohort were consistent with the results in the TCGA cohort, and the OS of the high-risk group was significantly lower than that of the low-risk group (P < 0.001). The AUCs for one-, three-, and five-year OS were 0.638, 0.553, and 0.532, respectively (Fig. 2F). In addition, we further verified the reliability of the model. We randomly dividing the training set into a verification 2 set (Figs. S1A–S1D) and a verification 3 set (Figs. S1E–S1H), the signature had reliable predictive ability (Fig. S1). Taking this together, the 10-gene signature was capable of predicting OS in AML. The clinical information of the patients was shown in Table S2.

Clinical correlation analysis

Univariate and multivariate COX analysis of clinically relevant factorsshowed that age (p < 0.001) and riskScore (p < 0.001) were independent prognostic indicators in the TCGA-LAML cohort (Figs. 3A, 3B), and that age (p < 0.001), runx1-mutation (p < 0.001), and riskScore (p = 0.019) were independent prognostic indicators in the GSE37642 cohort (Figs. 3C, 3D).

Nomogram analysis results of TCGA-LAML cohort and GSE37642 cohort

In order to better evaluate the relationship between genes and prognosis in the model, we used a nomogram to analyze it. The results show that in the TCGA-LAML cohort, BNIP3, CANX, and WDFY3 have a positive correlation with OS, and BAG3, CDKN2A, DIRAS3, NRG2, PARP1, PRKCD, and VAMP3 have a negative correlation with OS (Fig. 4A). In addition, in the GSE37642 cohort, CANX, CDKN2A, NRG2, and VAMP3 have a positive correlation with OS, and BAG3, BNIP3, DIRAS3, PARP1, PRKCD, and WDFY3 have a negative correlation with OS (Fig. 5A). The calibration plots showed that the nomogram could accurately predict the one-, three-, and five-year OS (Figs. 4B–4D, Figs. 5B–5D) with a harmonious consistency (TCGA-LAML, C-index = 0.72; GSE37642, C-index = 0.66) between the predicted and observed survival.

Significant differences between high- and low-risk subgroups

The patients were scored by autophagy-related gene models, and the patients were divided into high- and low-risk groups based on the optimal score. Principal components analysis (PCA) supports the classification of AML patients into two subgroups (Fig. 6A). In order to further analyze the difference between the high-risk and low-risk subgroups, the ESTIMATE algorithm was used to analyze the TCGA-LAML tumor microenvironment. The results showed that high ImmuneScore was significantly associated with poor survival (Fig. 6B). Another important finding was that ImmuneScore and StromalScore were higher in the high-risk group (Fig. 6C). In addition, age was significantly correlated with both Immune Score and Stromal Score (Fig. 6D).

In order to explore the differences in immune infiltrating cells in the high- and low-risk subgroups, we used the CIBERSORT algorithm to analyze the composition of 22 immune cells in the TCGA-LAML cohort (Fig. S2) and analyzed the correlation between different immune infiltrating cells (Fig. S3). In addition, the difference in immune infiltrating cells between high and low-risk subgroups is shown in Fig. 6E. Further analysis showed that the high expression of mast cells resting was associated with a better prognosis and the NK cells activated with high expression was associated with a poor prognosis (Fig. 6F).

PDL1 (CD274) and CTLA4 play a very important role in the immunotherapy of AML. We found that the high-risk group had higher expression levels of PDL1 and CTLA4 (Fig. 6G). GSEA analysis results showed that KEGG CHEMOKINE SIGNALING PATHWAY, KEGG CELL ADHESION MOLECULES CAMS, KEGG CYTOKINE CYTOKINE RECEPTOR INTERACTION, KEGG HEMATOPOIETIC CELL LINEAGE, and KEGG INTESTINAL IMMUNE NETWORK FOR IGA PROC were enriched in the high-risk group (Fig. 6H). The results of drug sensitivity analysis showed that there are significant differences between 24 chemotherapy drugs between high-risk and low-risk patients, which may provide help for personalized treatment of AML patients (Fig. 7).