Endoplasmic reticulum stress-related signature predicts prognosis and immune infiltration analysis in acute myeloid leukemia

ABSTRACT

Objectives

To construct an endoplasmic reticulum stress-related prognostic risk score (RS) model to predict prognosis and perform a preliminary analysis of immune infiltration in patients with acute myeloid leukemia (AML).

Methods

The whole-genome expression data for AML and endoplasmic reticulum stress (ER stress)-related genes were downloaded from the GEO and GSEA databases, respectively. The samples were divided into death and survival groups, combined with clinical prognosis information. LASSO regression was used to construct a prognostic RS model. The Kaplan–Meier curve method was used to evaluate the association between different risk groups and actual survival prognosis information. A cox regression analysis was used to screen for independent survival prognostic clinical factors and construct a nomogram. CIBERSORT and ssGSEA was used for immune-related analysis.

Results

Eighteen ER-stress related genes were identified and a comprehensive network was constructed. Further, 5 CC, 8 MF, 17 BP, and 2 KEGG pathways were enriched. Ten optimal DEGs were obtained and a prognostic risk model was constructed. Compared to the low RS group, the OS values of the high RS group were significantly lower. A significant correlation between the different risk groups and the actual prognosis was demonstrated. Ten immune cells with significantly different distributions in different risk groups were screened. KEGG enrichment analysis showed that there were 5 signaling pathways in the high-risk group.

Conclusions

The RS model can effectively predict the prognosis and has clinical implications for the prognosis of AML, combined with the correlation between different RS groups and the immune microenvironment.

KEYWORDS:

• Acute myeloid leukemia (AML)
• endoplasmic reticulum stress
• prognosis
• immune infiltration

1. Introduction

Acute myeloid leukemia (AML) is a malignant disease of myeloid stem/progenitor cells. The primary causes of AML are genetic and cytogenetic alterations [1]. Recent evidence indicates that the identification of novel AML biomarkers could contribute to a better understanding of the molecular basis of this disease [2]. However, some patients with relapsed and refractory AML and older patients have a poor prognosis, which is currently a challenging issue. The WHO also recently released the fifth edition of the Classification of Hematolymphoid Tumours [3], and the European Leukemia Net (ELN) recently published a revised system (ELN-2017) [4]. The authors used molecular genetic changes as risk stratification criteria, such as the presence of NPM1, FLT3-ITD, and CEBPA, to indicate a favorable prognosis and the presence of TP53, RUNX1, and ASX1, to indicate a poor prognosis. With the development of many molecularly targeted drugs, new mutated genes to modify the genetic risk stratification system and accommodate the progression of treatment strategies must be continued.

The endoplasmic reticulum (ER) is present in various eukaryotic cells and comprises the biofilm system. It is the basis of intracellular protein synthesis and an important site for post-translational modification, folding, and transport of proteins [5, 6]. Tumor cells are often in a microenvironment of ischemia, hypoxia, nutrient deficiency, and low pH during the growth process, which leads to the accumulation of a large number of unfolded or misfolded proteins, causing ER stress. Recent studies have confirmed that ER stress occurs in various tumor cells, leading to various pathological states. ER stress activates three major signaling pathways: the unfolded protein response (UPR), ER overload response (EOR), and sterol regulatory cascades. UPR activation also occurs during hematopoiesis, and increasing evidence supports the cytoprotective role of ER stress in the emergence and proliferation of leukemic cells. Previous studies have shown that ER stress plays an important role in AML pathogenesis [7]. The increased basal UPR is induced, at least in part, by adverse growth conditions in the bone marrow, such as hypoxia [8]. However, the exact role of ER-stress in AML development and progression remains unclear. Therefore, a comprehensive understanding of ER stress can contribute to a better diagnosis and help develop effective treatments for AML.

In this study, we analyzed ER-related genes in AML, explored the functions of related genes, and constructed an ER stress-related prognostic risk model. Further, the prognostic risk model was integrated with the clinical parameters, and a survival curve and normality were determined to verify the clinical predictive ability of the model. A preliminary immune infiltration analysis was also performed. The flow of the analysis is shown in Figure 1. Overall, a comprehensive understanding of ER stress can contribute to a better diagnosis and prognosis of AML.

Figure 1. Fow chart of analysing.

2. Materials and methods

2.1. Research objects

GSE106291[9, 10] and GSE37642 [11–13] genome-wide expression level datasets, including 250 AML bone marrow samples and 417 AML bone marrow samples with clinical survival prognosis, were downloaded from the NCBI GEO database (http://www.ncbi.nlm.nih.gov/geo/). GSE106291 and GSE37642 were used as training and validation datasets, respectively. ER stress-related genes were downloaded from the MSigDB of the Gene Set Enrichment Analysis database (http://software.broadinstitute.org/gsea/downloads.jsp) [14]. The GPL18460I Illumina HiSeq 1500 (Homo sapiens) and GPL96 Affymetrix Human Genome U133A Array were used for the analysis.

2.2. Differentially expressed ER stress-related genes search

Following the screening of genes related to ER stress, the expression levels of the corresponding genes in the samples were extracted. The samples combined with the clinical prognosis information were divided into the death group and survival group [15–17]. T-test in R3.6.1 was used to compare ER stress-related genes between the two groups, and P < 0.05 was considered statistically significant for the screening threshold.

2.3. Network construction and GO function and KEGG enrichment analysis

The STRING [18] database (Version:11.0, http://string-db.org/) was used to search for interactions between DEGs and their protein products, and an interaction score higher than 0.4 was selected as the screening threshold. The cor function (http://77.66.12.57/R-help/cor.test.html) was calculated using the cor coefficient (PCC). Relationship pairs with a P < 0.05 and an absolute PCC > 0.3 were selected as significant correlation co-expression relationship pairs. A comprehensive network was built using Cytoscape Version 3.6.1(http://www.cytoscape.org/) [19]. Finally, GO function and KEGG signaling pathway enrichment analysis of DEGs was performed using DAVID Version 6.8（The Database for Annotation, Visualization and Integrated Discovery (https://david.ncifcrf.gov/) [20, 21].

2.4. Prognostic model construction and validation

The DEGs and clinical survival prognosis information were combined to screen overall survival prognosis using Cox regression analysis (http://bioconductor.org/packages/survivalr/) [22]. Survival regression analysis was performed to determine the optimal combination of ER stress-related genes (https://cran.r-project.org/web/packages/lars/index.html) [23]. According to the LASSO prognostic coefficient of each element in the optimized ER stress-related gene combination and the expression level of the target gene in the dataset, the risk score (RS) model was constructed to calculate the RS as follows: $\mathrm{R}\mathrm{i}\mathrm{s}\mathrm{k}\mathrm{s}\mathrm{c}\mathrm{o}\mathrm{r}\mathrm{e}\left(\mathrm{R}\mathrm{S}\right)=\sum \mathrm{C}\mathrm{o}\mathrm{e}{\mathrm{f}}_{\mathrm{g}\mathrm{e}\mathrm{n}\mathrm{e}\mathrm{s}}\times \mathrm{E}\mathrm{x}{\mathrm{p}}_{\mathrm{g}\mathrm{e}\mathrm{n}\mathrm{e}\mathrm{s}}$where Coef_genes represents the LASSO prognostic coefficient of the target gene and Exp_genes. represents the expression level of the target gene.

The RS values of the training and validation datasets were calculated, and using the RS median value as the dividing line, the samples were divided into high (RS score ≥ RS median value) and low (RS < RS median value) groups. The association of the actual survival prognosis information between the groups with high and low scores was evaluated using the Kaplan–Meier (KM) curve.

2.5. Relationship between risk scoring and clinical factors and construction of nomogram survival model

Analysis of the distribution of clinical information in samples in different risk groups in the training and validation datasets and screening of independent prognostic clinical factors were performed using Cox regression [22]. A nomogram was constructed to predict one-, three-, and five-year survival rates combined with independent prognostic factors and risk information (https://cran.r-project.org/web/packages/rms/index.html) [24, 25] and to calculate the C-index coefficient (C-index can be regarded as the score of all individual pairs correctly ranked based on Harrell C statistics) [25] for predicting survival time[ 26]. (http://www.bioconductor.org/packages/release/bioc/html/survcomp.html) [27].

2.6. Immunological correlation evaluation

CIBERSORT (https://cibersort.stanford.edu/index.php) [28] was used to evaluate immune cell type ratios and compare differences among the immune cell distributions of various risk groups. The immune score, stromal score, and tumor purity were calculated and compared among different risk groups (http://127.0.0.1:29606/library/estimate/html/estimateScore.html) [29]. Finally, the correlation between the immune cells, the four estimated scores, and the expression levels of the 10 DEGs was calculated using the cor function.

All KEGG information for related genes was acquired from the GSEA database [14]. According to genome-wide expression data (http://www.bioconductor.org/packages/release/bioc/html/GSVA.html) [30], quantitative analysis was performed for each KEGG signaling pathway, the distribution difference of each pathway was compared, and the correlation between significantly different signaling pathways was analyzed.

2.7. Statistical analysis

All statistical analyses were performed using R3.6.1. A Student’s t-test was used to compare the two groups. Univariate and multivariate Cox regression analyses were used for screening ER stress-related genes and independent prognostic factors for overall survival. The Wilcoxon rank-sum test was used to compare two paired groups. The cor function was calculated for correlation analysis. P < 0.05 was considered statistically significant.

3. Results

3.1. Identification of significantly DEGs

In total, 295 ER stress-related genes were identified and the expression levels of the corresponding genes in the samples were extracted from the GSE106291 expression profile data. The samples were divided into death (n = 100) and survival (n = 150) groups based on clinical prognosis information. Further, 18 ER stress-related genes were identified, including 8 downregulated genes and 10 downregulated genes (Figure 2A).

Figure 2. Differential analysis ER-stress-related genes and network construction. (A) Violin plot of 18 ER-stress-related genes between survival and death groups; (B) Significantly differentially expressed ER-stress-related genes integrated network; (C) Histogram of BP, CC, MF and KEGG signaling pathways with significant gene correlation in the comprehensive network.

3.2. Network construction and GO function and KEGG enrichment analysis

Thirty-eight linked pairs of gene–protein relationships were obtained via network construction. The PCC of ER stress-related genes was then calculated, and 67 linked pairs of co-expression relationships were obtained and integrated to construct a comprehensive network (Figure 2B).

GO function and KEGG signaling pathway enrichment showed that a total of 5 cellular components (CC), 8 molecular functions (MF), 17 biological processes (BP), and 2 KEGG pathways were screened (Figure 2C). The most enriched signaling pathways in CC were ‘endoplasmic reticulum,’ ‘integral component of membrane,’ and ‘endoplasmic reticulum membrane;’ those in MF were ‘protein kinase bindin,’ ‘ubiquitin protein ligase activity,’ ‘ubiquitin protein ligase binding,’ and ‘enzyme;’ those in BP were ‘endoplasmic reticulum unfolded protein response,’ ‘ubiquitin-dependent ERAD pathway,’ and ‘cellular response to glucose starvation.’ ‘Protein processing in endoplasmic reticulum’ and ‘p53 signaling pathway’ were enriched in KEGG.

3.3. ER stress-related genes prognostic model construction

Fifteen DEGs were screened based on clinical survival prognosis information (Figure 3). To prevent overfitting and eliminate highly correlated ER stress-related genes, LASSO regression analysis was performed. Ten optimal DEGs were obtained using the LASSO algorithm (Figure 4A). These parameters are shown in Figure 4(B). Based on the median expression of the 10 genes, the samples were divided into high- and low-expression groups. A KM curve analysis was used to investigate the correlation between expression levels and survival prognosis (Figure 4C). High expression of FBXO44, TRIB3, RNF139, BRSK2, TRIM25, and FKBP14 was associated with low OS values, indicating that these six genes are risk factors, whereas high expression of MBTPS1, TBL2, VAPB, and PDIA4 was associated with high OS values, indicating that these four genes are protective factors. Based on the LASSO regression coefficients of the 10 genes and their expression levels, the following RS calculation formula was constructed: $\begin{array}{rl}\mathrm{R}\mathrm{S}& =\left(-0.055701696\right)\ast \mathrm{E}\mathrm{x}{\mathrm{p}}_{\mathrm{M}\mathrm{B}\mathrm{T}\mathrm{P}\mathrm{S}1}\\ & +\left(-0.01367718\right)\ast \mathrm{E}\mathrm{x}{\mathrm{p}}_{\mathrm{V}\mathrm{A}\mathrm{P}\mathrm{B}}\\ & +\left(-0.00820906\right)\ast \mathrm{E}\mathrm{x}{\mathrm{p}}_{\mathrm{P}\mathrm{D}\mathrm{I}\mathrm{A}4}\\ & +\left(0.006963321\right)\ast \mathrm{E}\mathrm{x}{\mathrm{p}}_{\mathrm{F}\mathrm{K}\mathrm{B}\mathrm{P}14}\\ & +\left(0.008961944\right)\ast \mathrm{E}\mathrm{x}{\mathrm{p}}_{\mathrm{T}\mathrm{R}\mathrm{I}\mathrm{M}25}\\ & +\left(0.025160051\right)\ast \mathrm{E}\mathrm{x}{\mathrm{p}}_{\mathrm{R}\mathrm{N}\mathrm{F}139}\\ & +\left(0.046016902\right)\ast \mathrm{E}\mathrm{x}{\mathrm{p}}_{\mathrm{F}\mathrm{B}\mathrm{X}\mathrm{O}44}\\ & +\left(0.02998434\right)\ast \mathrm{E}\mathrm{x}{\mathrm{p}}_{\mathrm{T}\mathrm{R}\mathrm{I}\mathrm{B}3}\\ & +\left(-0.018042718\right)\ast \mathrm{E}\mathrm{x}{\mathrm{p}}_{\mathrm{T}\mathrm{B}\mathrm{L}2}\\ & +\left(0.018018557\right)\ast \mathrm{E}\mathrm{x}{\mathrm{p}}_{\mathrm{B}\mathrm{R}\mathrm{S}\mathrm{K}2}\end{array}$

Figure 3. Forest map of ERS-related genes significantly associated with survival outcomes.

Figure 4. (A) LASSO parameter diagram, the intersection of horizontal and vertical black lines in the following figure is the location of the optimal parameter (−4.501,0.218), and the corresponding gene is the optimal gene combination; (B) LASSO coefficients of 10 genes; (C) 10 Prognostic KM curves of the optimized ERS-related genes.

3.4. Efficacy evaluation and comparison of RS survival prognostic risk prediction model

The scatter plot of the model in the GSE106291 training set shows the samples in high- and low-risk groups. The RS value distribution and survival time distribution for each group are shown in Figure 5(B). The KM curves showed that the OS of the high-RS group was significantly lower than that of the low-RS group (Figure 5A). The ROC curves for the training set model yielded an acceptable AUC of 0.905 (Figure 5C). Although the results of our GSE37642 validation dataset remained consistent (Figures 5D,E), ROC curve analysis of the CGGA validation set model yielded an acceptable AUC value of 0.885 (Figures 5F).

Figure 5. GSE106291 training set: (A) KM curves based on 10 optimized ER-stress related gene prognostic models; (B) RS distribution and survival time status; (C) The ROC curves of 10 optimal ER-stress related gene prognostic models were plotted, and the numbers in parentheses represent the corresponding ROC curves. GSE37642 validation set: (D) KM curves based on 10 optimized ER-stress related gene prognostic models; (E) RS distribution and survival time status; (F) The ROC curves of 10 optimal ER-stress related gene prognostic models were plotted, and the numbers in parentheses represent the corresponding ROC curves.

3.5. Relationship between risk score and clinical information

The differences among various RS groups with different clinical characteristics were analyzed using the R language. The results showed that the RS differed significantly according to age, runx1-runx1t1fusion, and treatment response (Figure 6A–E), with patients aged >60 years having a higher RS with runx1-runx1t1_fusion and those with sensitivity to treatment response having a lower RS.

Figure 6. (A) Differential relationship between age clinical traits in risk group; (B) Differential relationship between sex clinical traits in risk group; (C) Differential relationship between runx1 mutation in risk group; (D) Differential relationship between runx1-runx1t1_fusion in risk group; (E) differential relationship between treatment response in risk group.

3.6. Construction of a nomogram survival model with independent prognostic factors

The screening results for independent prognostic clinical factors are shown in Table 1. Age, treatment response, and the RS model were significant independent prognostic factors (Figure 7A). The nomogram survival model is shown in Figure 7(B). The ‘Total points’ axis was used to predict the survival of the sample based on various clinical indicators. Calibration plots showed that the nomogram graphs were consistent with the ideal model runs by comparing the consistency between the one-, three-, and five-year survival rates (Figure 7C).

Figure 7. (A) Independent prognostic Mori diagram of clinical prognostic factors; (B) Nomogram of nomogram for survival prediction model with independent prognostic factors; (C) Nomogram (left to right) Consistency plots of predicted and actual survival rates at one、three、and five years. The horizontal axis shows the predicted and the vertical axis shows the actual survival rates;

Table 1. Prognostic relevance clinical prognostic factors.

Clinical characteristics	Uni-variables cox			Multi-variables cox
	HR	95%CI	P	HR	95%CI	P
Age(years,mean ± SD)	1.039	1.024–1.054	1.05E-07	1.023	1.006–1.040	6.57E-03
Gender(Male/Female)	1.188	0.863–1.636	2.91E-01	–	–	–
Treatment response(Resistant/Sensitive)	0.189	0.134–0.269	2.00E-16	0.272	0.187–0.396	1.15E-11
runx1-runx1t1_fusion(Yes/No)	0.606	0.193–1.901	3.85E-01	–	–	–
runx1_mutation(Yes/No)	1.759	1.202–2.573	3.62E-03	1.181	0.785–1.774	4.26E-01
RS model status(High/ Low)	3.249	2.299–4.593	1.86E-12	2.722	1.873–3.956	1.50E-07

3.7. Immune infiltration analysis

Ten immune cells with significantly different distributions in the different risk groups were screened (Figure 8A). Figure 8(B) shows that the distribution of Stromal, Immune, and ESTIMATE Scores in the high-risk group was higher than that in the low-risk group, but TumorPurity in the low-risk group was higher than that in the high-risk group. Correlation analysis showed that TRIM25 was positively and negatively correlated with Immune Score and TumorPurity, respectively. Macrophage M2 was positively and negatively correlated with PDIA4 and T cell CD4 + memory activation, respectively (Figure 8C).

Figure 8. (A) Distribution of 10 immune cells with significantly different distribution in different risk groups; (B) Distribution of the estimated scores in different risk groups; (C) Heat map of correlation between expression levels of 10 DEGs with 10 immune cells, 4 ESTIMATE scores;Scatter plot of the best positive and negative correlations on the side.

3.8. GSVA quantitative KEGG analysis based on genome-wide expression levels

We identified 21 KEGG signaling pathways with significant differences between the different risk groups (Figure 9A). Twelve signaling pathways were enriched in the low-risk group, including ‘valine leucine and isoleucine biosynthesis,’ ‘aminoacyl tRNA biosynthesis,’ ‘non homologous and joining,’ ‘homologous recombination,’ ‘propanoate metabolism,’ ‘mismatch repair,’ ‘valine leucine and isoleucine degradation,’ ‘selenoamino acid metabolism,’ ‘sulfur metabolism,’ ‘cysteine and methionine metabolism,’ ‘histidine metabolism,’ and ‘DNA replication’. Nine signaling pathways were enriched in the high-risk group, including ‘cytokine cytokine receptor interaction,’ ‘thyroid cancer,’ ‘natural killer cell-mediated cytotoxicity,’ ‘folate biosynthesis,’ ‘type I diabetes mellitus,’ ‘allograft rejection,’ ‘glycosaminoglycan biosynthesis keratan sulfate,’ ‘kegg graft versus host disease,’ and ‘kegg circadian rhythm mammal’. The correlation among KEGG signaling pathways with significant differences in distribution showed that the PDIA4 gene was positively correlated with the valine, leucine, and isoleucine degradation signaling pathways and negatively correlated with the mammalian circadian rhythm signaling pathway (Figure 9B).

Figure 9. (A) Heatmap of KEGG signaling pathways with significant differences between risk groups; (B) Heat map of the correlation between the quantified KEGG signaling pathways and the expression levels of the 10 genes that constitute the RS model. Scatter plot of the best positive and negative correlations on the side.

3.10. KEGG enrichment analysis related to risk groups

Based on the genome-wide expression level of the GSE106291 tumor samples, KEGG signaling pathways significantly related to risk grouping were screened. The result showed that 20 significantly related KEGG signaling pathways were screened out, of which 5 and 15 pathways were enriched in the high-risk and low-risk groups, respectively (Figure 10). The five pathways in the high-risk group were ‘Allograft rejection,’ ‘glycosaminoglycan biosynthesis keratan sulfate,’ ‘graft versus host disease,’ ‘prostate cancer,’ and ‘type I diabetes mellitus.’

Figure 10. Twenty KEGG signaling pathways were significantly associated with different risk groups.

4. Discussion

Recent studies have shown that ER stress is involved in AML pathogenesis [7]. Simultaneously, the ER stress response in myeloid leukemia cells has been associated with advanced disease stages and cell resistance [31, 32]. Although the role of ER stress in the development of AML is well-established, the extent to which ER stress affects the prognosis of AML remains unclear.

In this study, ER stress-related genes were analyzed in AML. First, 18 ER stress-related genes were identified, including 8 downregulated genes and 10 downregulated genes. A comprehensive network was constructed, including 38 pairs of gene and protein relationships, 67 pairs of co-expression relationships A total of 5 CC, 8 MF, 17 BP, and 2 KEGG pathways were enriched. The most enriched signaling pathways in CC were ‘endoplasmic reticulum,’ ‘integral component of membrane,’ and ‘endoplasmic reticulum membrane;’ those in MF were ‘protein kinase binding,’ ‘ubiquitin protein ligase activity,’ ‘ubiquitin protein ligase binding,’ and ‘enzyme;’ those in BP were ‘endoplasmic reticulum unfolded protein response,’ ‘ubiquitin-dependent ERAD pathway,’ and ‘cellular response to glucose starvation.’ ‘Protein processing in endoplasmic reticulum’ and ‘p53 signaling pathway’ were enriched in KEGG. Further, 10 ER stress-related genes were optimized, among which MBTPS1, TBL2, VAPB, and PDIA4 were protective factors, and FBXO44, TRIB3, RNF139, BRSK2, TRIM25, and FKBP14 were risk factors. We then combined the DEGs with their corresponding clinical information to establish an RS model. In the prognostic prediction mode, MBTPS1 and FBXO44 had the greatest impact on prognosis. The higher the expression level of MBTPS1, the better the prognosis. The MBTPS1 gene, also known as the membrane-bound transcription factor peptidase site 1 gene, encodes Site-1 protease (S1P). S1P regulates lipogenesis, ER function, lysosomes, and biogenesis in mice and cultured cells. Studies have found [33] that deletion of MBTPS1 can impair ER function. The investigators also developed a predictive model that included six membrane protein genes in patients with AML without cytogenetic abnormalities, among which MBTPS1 was a risk factor. However, the higher the expression level of FBXO44, the worse the prognosis. FBXO44 belongs to the F-box protein family (FBA family) and shares a conserved G domain for substrate binding. In addition to sequence homology, all members of the FBA family, except FBXO44, can recognize the glycosylation substrate FBXO44, which has a substrate-recognition function during ubiquitin-mediated proteolysis [34–36]. A misfolded protein response occurs when the ER lumen misfolds or when excessive protein synthesis occurs. Studies [37] have shown that FBXO44 is an E3 ubiquitin ligase responsible for BRCA1 degradation. FBXO44 promotes the silencing of DNA replication-coupled repetitive elements in cancer cells [38].

In addition, we evaluated the performance of the risk prediction models. The KM curves showed that the OS of the high RS group was significantly lower than that of the low RS group. Meanwhile, we explored the relationship between the RS and clinical information, which revealed that the RS differed significantly according to age, runx1-runx1t1 fusion, and treatment response. Compared with patients with low RS scores, patients with high RS scores had worse prognosis in the training and validation groups, which further proved that ER stress-related gene markers have high prediction accuracy and robustness and that the RS is an independent prognostic factor. Furthermore, we constructed a quantitative and objective nomogram based on a multivariate analysis. We analyzed the consistency between the one-, three-, and five-year survival rates predicted by the nomogram survival model and the actual survival rates, verifying that our nomogram survival model has high predictive value and is suitable for clinical applications. The ROC curve analysis of the training set and CGGA validation set models yielded acceptable AUC values of 0.905 and 0.885, respectively. AUC values of 0.8–0.9 indicate good sensitivity and specificity. The results of the ROC analysis indicated that the metrics of our nomogram graphs are steadily predictive. Therefore, ER stress-related prognostic markers can be used as potential prognostic biomarkers for AML.

Immune escape has become an important factor affecting the long-term disease-free survival of patients with AML [39]. Therefore, we conducted a preliminary study on the immune infiltration of AML, and 10 immune cells with significantly different distributions in different risk groups were screened. In the high-risk group, CD8+ T cells, CD4 + memory activated T cells, regulatory T cells (Tregs), resting natural killer (NK) cells, and Macrophage M0 immune cells infiltrated in a high proportion. Whereas, in the low-risk group, B cell plasma, CD4 + naïve T cells, T cell gamma delta, macrophage M2, and mast cells infiltrated in a high proportion. A previous study found that the T-cell and Tregs population in the bone marrow of patients with AML were higher than those in healthy donors. Moreover, targeting immune receptors provides a theoretical basis for T-cell immunotherapy in AML [40]. In AML, T cells are not functionally impaired; increased PD-1 expression is observed only at relapse and is associated with shifts to memory T-cell compartments [41]. Tregs play an important role in tumor immune escape and angiogenesis [42]. NK cells play a key role in inducing an effective adaptive anti-tumor immune response. However, in patients with AML, significantly impaired NK cell function can prevent immune surveillance and affect prognosis [43]. Various immune escape mechanisms have been reported in patients with AML, including decreased expression of NK cytotoxic receptors, downregulation of NK ligands, and secretion of soluble NK inhibitory factors [44]. Monocytes and macrophages are important immune cell types involved in tumor immune escape. Monocytes have been reported to differentiate into macrophages. Stimulation with cytokines and chemokines leads to the differentiation of macrophages into classically activated macrophages (M1) and alternately activated macrophages (M2). Tumor-associated macrophages (TAMs), which are important components of tumor-infiltrating immune cells, have been associated with the growth, angiogenesis, and metastasis of a variety of cancers, most likely through TAMs polarization to M2 (an alternative phenotype). Clinical and experimental evidence has shown that highly invasive tumor tissues with TAMs are associated with poor prognosis and drug resistance. Therefore, targeting TAMs is considered a promising immunotherapeutic strategy for tumors [45]. Some researchers have reported that the degree of mast cell infiltration positively correlates with tumor size, grade, and metastasis [46].

Subsequently, we correlated 10 immune cells with significantly different distributions, 4 ESTIMATE scores, and the expression levels of 10 genes constituting the RS model. The results showed that the distribution of stromal, immune, and ESTIMATE scores in the high-risk group was higher than that in the low-risk group. However, TumorPurity in the low-risk group was higher than that in the high-risk group. TRIM25 positively correlated with Immune Score and negatively correlated with TumorPurity. Macrophage M2 positively correlated with PDIA4 but negatively correlated with CD4 + T-cell memory activation. Based on the GSVA quantitative KEGG analysis of genome-wide expression levels, 21 KEGG signaling pathways with significant differences among the different risk groups were identified. The results showed that PDIA4 was positively correlated with valine, leucine, and isoleucine degradation signaling pathways. There was a negative correlation with mammalian circadian rhythm signaling pathway. Valine, leucine, and isoleucine degradation play a key role in the regulation of energy homeostasis, nutrient metabolism, gut health, immunity, and disease in humans and animals. Genetic research related to the mammalian circadian rhythm signaling pathway has revealed that the circadian system is closely linked to the processes that control sleep and metabolism.

TRIM25, a member of the TRIM (tripartite motif) protein subfamily (a subfamily of Ring-type E3 ubiquitin ligases), plays an important regulatory role in carcinogenesis. Dysregulation of the ubiquitin-mediated degradation of oncogene products or tumor suppressors may be related to the etiology of cancer and leukemia [47]. Studies have found that TRIM25 is the most significant inducer gene in ER-stress [48]. As a feedback mechanism, TRIM25 is required for tumor cell defense against ER stress. TRIM25 promotes tumor cell survival. PDIA4 is a protein disulfide isomerase that catalyzes the formation of disulfide bonds that can activate platelets, promote thrombosis, and respond to ER stress. In addition, current studies have shown that activated platelets and ER stress responses influence tumor progression, and PDIA4 is involved in tumor development by affecting apoptosis and DNA repair mechanisms. Therefore, PDIA4 inhibitors may have great potential for inhibiting tumor progression [49].

KEGG signaling pathways that were significantly related to risk groups were identified. The result showed that 20 significantly related KEGG signaling pathways were screened out, of which 5 pathways were enriched in the high-risk group. The five pathways in the high-risk group were ‘allograft rejection,’ ‘glycosaminoglycan biosynthesis keratan sulfate,’ ‘graft versus host disease,’ ‘prostate cancer,’ and ‘type I diabetes mellitus.’

Our study had some limitations. First, the results were analyzed using bioinformatics and need to be experimentally verified further. Second, the robustness of this prognostic risk model must be verified in prospective clinical studies with large sample sizes.

5. Conclusions

We developed and validated an ER stress-associated gene marker for the prognostic stratification of patients with AML. In addition, this biomarker can be used to predict patient survival, providing a new strategy for the prognostic prediction of AML. Survival outcomes and immune cell infiltration levels differed among the different RS groups. By combining prognostic risk models and clinical parameters to construct nomogram plots, our results provide new insights into the identification of prognostic biomarkers.

Acknowledgements

Author contributions: Lu Dong: research design and performance, data collection, data analysis and manuscript-writing; Zefeng Miao: data analysis; Haili Wang, Yanhui Yu, Dongzheng Gai, Guoxiang Zhang and Li Ge: data recheck; Xuliang Shen: research design and manuscript review.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Data availability statement

Our data and related clinical information were retrieved from NCBI Gene Expression Omnibus database (GEO http://www.ncbi.nlm.nih.gov/geo/).

Additional information

Funding

This work was supported by the Shanxi Province ‘1331’ Project Key Discipline Construction Project under Grant (Jin Jiao-Research [2017] No 3); the Subject of Health Commission of Shanxi Province under [grant number2021001]; the 2019 Science and Technology Innovation Team Project of Changzhi Medical College under [grant number CX201901]; and the 2022 Shanxi Higher Education Science and Technology Innovation Project [grant number 2022L367].