Construction of a solid Cox model for AML patients based on multiomics bioinformatic analysis

Acute myeloid leukemia (AML) is a highly heterogeneous hematological malignancy. The bone marrow (BM) microenvironment in AML plays an important role in leukemogenesis, drug resistance and leukemia relapse. In this study, we aimed to identify reliable immune-related biomarkers for AML prognosis by multiomics analysis. We obtained expression profiles from The Cancer Genome Atlas (TCGA) database and constructed a LASSO-Cox regression model to predict the prognosis of AML using multiomics bioinformatic analysis data. This was followed by independent validation of the model in the GSE106291 (n=251) data set and mutated genes in clinical samples for predicting overall survival (OS). Molecular docking was performed to predict the most optimal ligands to the genes in prognostic model. The single-cell RNA sequence dataset GSE116256 was used to clarify the expression of the hub genes in different immune cell types. According to their significant differences in immune gene signatures and survival trends, we concluded that the immune infiltration-lacking subtype (IL type) is associated with better prognosis than the immune infiltration-rich subtype (IR type). Using the LASSO model, we built a classifier based on 5 hub genes to predict the prognosis of AML (risk score = -0.086×ADAMTS3 + 0.180×CD52 + 0.472×CLCN5 - 0.356×HAL + 0.368×ICAM3). In summary, we constructed a prognostic model of AML using integrated multiomics bioinformatic analysis that could serve as a therapeutic classifier.

Introduction

Acute myeloid leukemia (AML) is a highly heterogeneous group of hematological malignancies that are characterized by various cytogenetic and molecular heterogeneities (1, 2). Although substantial progress has been achieved with combinatorial therapies including radiation, chemotherapy, immunotherapy and/or targeted therapy, the cure rate of patients remains only 35%-40% in younger patients (age< 60 years) and 5%-15% in older patients (age< 60 years) (3). Relapse and refractory disease continue to be major obstacles in the treatment of AML, with 29% or fewer patients living beyond 5 years.

Several studies have shown that changes in the bone marrow (BM) microenvironment in AML largely promote distinct biological processes in leukemogenesis, drug resistance and leukemia relapse (4). Thus, insights into BM microenvironment action may provide better diagnosis and treatment strategies for AML patients.

The BM microenvironment in AML is comprised of leukemia cells, stromal cells, endothelial cells and distinct immune cell subsets. Among them, vascular endothelial cells (ECs) promote leukemia cell proliferation, drug resistance, and recurrence through paracrine vascular endothelial growth factor (VEGF), adhesion, and fusion with leukemia cells, resulting in poor prognosis (5–8). Stromal cells can promote chemotherapy resistance in leukemia cells through ligand-receptor interactions (9) such as SDF-1/CXCR4 (10) and VLA-4/VCAM-1 (11). Adipocytes promote the proliferation, growth and chemotherapeutic resistance of leukemia cells by breaking down stored triglycerides into free fatty acids (12) and secreting tumor-related proinflammatory cytokines (13).

The leukemia immune microenvironment presents with immune dysregulation and suppression, leading to an imbalance of suppressor T cells and effector T cells, T cell exhaustion and an increase in myeloid-derived suppressor cells (MDSCs) and leukemia-supporting macrophages compared to normal bone marrow tissue (14). Recent studies on the characterization of the leukemia immune microenvironment could aid in the search for novel prognostic biomarkers and potential therapeutic targets (15, 16). In addition, treatments such as chemotherapy, immunotherapy, and combination therapy to alter the immune microenvironment of AML have been widely used (17, 18). However, different immune cells have different effects in AML. Therefore, understanding the distribution and function of immune-related genes in different immune cells is of great significance to further explore the BM immune microenvironment of AML patients.

Here, we investigated the impact and potential mechanisms of immune-related genes on the prognosis of AML. We are the first group to screen for hub genes in this disease using multiomics analysis. We constructed a LASSO-Cox regression model to predict the prognosis of AML according to the characterization of the leukemia immune microenvironment. The distribution of hub genes in immune cells was revealed through single-cell sequencing data and may provide the potential for precise patient stratification and treatment.

Materials and methods

Datasets

The test cohort of AML was downloaded from The Cancer Genome Atlas (TCGA) database (https://www.cancer.gov/) and includes mRNA data from 151 cases (RNASeq V2), miRNA data from188 cases (Illumina HiSeq miRNAseq) and Illumina Human Methylation450 Bead Array data from 140 cases. Samples were selected for the study according to the following criteria: 1) acute myelocytic leukemia was pathologically diagnosed, 2) all three kinds of data were available for the patient, and 3) the clinical information was complete. Finally, 97 patients were selected for our following study.

The other datasets were obtained from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/). The independent validation cohort, GSE106291 dataset (251 samples), which was generated using the GPL18460 (Illumina HiSeq 1500, Homo sapiens) platform.

The single-cell RNA sequence dataset GSE116256, including 16 untreated samples (D0), was used to reveal the distribution of hub genes in immune cell types. The immune gene set, including 776 genes, was acquired from a previous study (19).

Screening of candidate genes and hierarchical clustering

Differential mRNA and miRNA expression were analyzed by the DESeq2 (20) function (P<0.05, |logFC|>1). For each probe in the methylation data, the value shown is the β value (β =U/(M+U+1)), where M is the methylated probe signal strength and U is the unmethylated probe signal value. The methylmix package (21) was used to analyze the correlation between the gene methylation level and mRNA expression value (Pearson correlation coefficient test, R>0.5, P<0.05). Unsupervised hierarchical clustering was performed (Euclidean distances and ward.D2 method) based on survival-related immune genes (SIGs) to establish an immunogenomic classification of TCGA-AML patients.

Immune infiltration analysis

The enrichment scores of 28 immune signatures in each AML sample were quantified using single-sample gene set enrichment analysis (ssGSEA) (22). Stromal, immune, and estimate scores were further calculated to evaluate tumor purity and immune cell infiltration in tumor tissues based on the mRNA expression data using the ESTIMATE (Estimation of STromal and Immune cells in MAlignant Tumor tissues using Expression data) algorithm (23). The ESTIMATE algorithm is based on the immune signature genes and stromal signature genes in the mRNA data, and then calculates the immune score and stromal score by ssGSEA.

Protein–protein interaction network construction and gene ontology functional enrichment analysis

mRNA interaction data were obtained from the STRING database (https://string-db.org/) (24). The threshold of interaction score is 0.4. The PPI network was established using Cytoscape software (25). The Database for Annotation, Visualization and Integrated Discovery (DAVID, https://david.ncifcrf.gov/) (26, 27) was used for GO enrichment analysis (P < 0.05), and the results were plotted by the GO plot package (28).

Survival analysis and prognostic model construction

A Cox regression model was constructed to identify immune-related genes that significantly correlated with the OS (SIGs) of AML patients in TCGA. The genes that met the standard of P < 0.05 were used for subsequent study. Least absolute shrinkage selection operator (LASSO)-penalized Cox regression analysis (29, 30) was applied to distinguish the most important SIGs and to construct a prognostic model for AML based on a linear combination of the regression coefficients. Kaplan-Meier survival curves and receiver operating characteristic (ROC) curves were used to test and validate the performance of the classifier.

scRNA dataset analysis

We employed the Seurat (31, 32) and SingleR (33) packages to generate Uniform Manifold Approximation and Projection (UMAP) plots and reveal the distribution of hub genes in each immune cell type.

Patients and samples

A total of 55 patients with newly diagnosed AML were enrolled between Jan 2016 to Dec 2020 at the Xinqiao Hospital of the Army Medical University in China. Samples were selected for the study according to the following criteria: 1) acute myelocytic leukemia diagnoses were based on morphological findings, karyotype, and immunophenotypic features of leukemia cells by consultant hematologist, 2) patients were newly diagnosed, and remained untreated at the time of collection, 3) the clinical information was completed.

Molecular docking

The virtual screening of molecular docking was performed using AutoDock Vina 1.1.2 (34) to predict the most likely optimal ligands. The three-dimensional structure of hub genes was retrieved from the Protein Data Bank (https://www.rcsb.org/). A library of 2115 FDA approved compounds were extracted from ZINC15 druglike database (http://zinc.docking.org/). The visualization of active interactions between proteins and compounds was performed by Discovery Studio Visualizer v4.5.0 (BIOVIA).

Results

Classification of AML based on immune-related genes that significantly affect patient prognosis

For a more extensive study of immune genes in AML, we retrieved transcriptome, microRNA, and DNA methylation profile data and integrated clinical information for 97 samples from TCGA database (Table S1). And the flowchart of analysis is shown in Figure S1. A Cox proportional hazard regression model was employed to analyze 776 immune-related genes (19) in the mRNA expression data of 97 samples (P<0.05), and 98 survival-related immune genes (SIGs) that significantly affected the survival of AML patients were identified (Table S2).

Using unsupervised clustering analysis (Euclidean statistics and complete linkage method) of 98 SIGs, those 97 samples were clustered into three distinct immune subtypes (Im1: immune cluster 1, Im2: immune cluster 2, Im3: immune cluster 3) based on the gene expression signature (Figure 1A). As shown in immune gene heatmaps, most of the SIGs were highly expressed in the Im1 and Im3 clusters but expressed at low levels in the Im2 cluster (Figure 1B). Kaplan-Meier survival analysis revealed that the prognosis of the Im2 cohort was significantly better than that of the Im1 and Im3 cohorts (Figure 1C, log-rank test, P=0.0008).

FIGURE 1

Figure 1 Unsupervised clustering analysis of AML patients based on 98 survival-related immune genes. (A) All 97 TCGA-AML patients were divided into 3 clusters (green: Im1 cluster, red: Im2 cluster, blue: Im3 cluster). (B) Heatmap of 98 survival-related immune genes in different AML clusters. (C) The Kaplan-Meier survival analyses along with the Log-rank test were used to compare the overall survival (OS) of the Im1, Im2 and Im3 clusters.

The immune infiltration was significantly different in different cluster

As the immune microenvironment was significantly correlated with the occurrence and development of AML, a single-sample gene set enrichment (ssGSEA) algorithm was utilized to explore differences in the immune microenvironment among the three immune clusters. The results showed that the Im2 cluster had fewer infiltrating immune cells than the Im1 and Im3 clusters (Figure 2A), and the specific infiltration of immune cells in different clusters is shown in Figure S2. Consistent findings demonstrated that the immune scores were significantly lower in the Im2 cluster (Figure 2B, unpaired t test, P < 0.001), while tumor purity was significantly higher in the Im2 cluster but significantly lower in the Im1 and Im3 clusters (Figure 2C, unpaired t test, P < 0.001). Generally, we can conclude that patients with less immune infiltration and lower immune scores may have a better prognosis than those with more immune infiltration and higher immune scores.

FIGURE 2

Figure 2 Immune functional characters of 3 AML patients clusters. (A) Heatmap of the Im1, Im2 and Im3 cohorts of 97 TCGA-AML patients using ssGSEA scores from 28 immune cell types. Violin plots depict the immune score (B) and tumor purity (C) of the Im1, Im2 and Im3 cohorts (***: P<0.001).

42 DEG-SIGs were screened by mRNA expression data analyzing

Based on the significant differences in immune infiltration and survival trends between the Im2 cluster and Im1/3 cluster, we defined Im2 as the immune infiltration-lacking subtype (IL type) and Im1/3 as the immune infiltration-rich subtype (IR type). To reveal the potential mechanisms of different prognoses between IL and IR subtypes, an elaborate analysis of the mRNA expression profiles of the two types of AML patients was implemented. We performed differentially expressed gene analyses and identified 1936 differentially expressed genes (DEGs) with significant differences between IL and IR subtypes. There were 42 SIG-DEGs which were common members of 1936 DEGs and 98 SIGs (P < 0.05, |Fold Change|>2 or <0.5) (Table 1) (Figures 3A, B).

TABLE 1

Table 1 Common 42 intersecting genes of differentially expressed genes (DEGs) and survival-related immune genes (SIGs) between immune infiltration-lacking subtype (IL type) and immune infiltration-rich subtype (IR type) (P < 0.05, | Fold Change |>2 or <0.5).

FIGURE 3

Figure 3 Difference analysis of the mRNA expression dataset from TCGA-AML patients. (A) Volcano plot of differentially expressed genes between immune infiltration-lacking subtype (IL type) and immune infiltration-rich subtype (IR type). (B) Venn plot of the intersecting genes between differentially expressed genes (DEGs) and survival-related immune genes (SIGs). (C) PPI network of 42 overlap genes DEG-SIGs. (D) The bubble plot represents the GO functional enrichment analysis of 42 DEG-SIGs.

To elucidate the mechanism of prognosis difference between IL and IR subtype we obtained the interaction data of the 42 DEG-SIGs from STRING website (https://string-db.org/) (interaction score > 0.4), and then constructed protein-protein interaction (PPI) network by Cytoscape.(Figure 3C). Gene ontology (GO) functional enrichment analysis distinguished some enriched terms in three subontologies: biological processes (BP), cellular component (CC), and molecular function (MF) (Figure 3D). For BP, 42 DEG-SIGs were enriched in defense response, inflammatory response, and immune system process. With regard to CC, 42 DEG-SIGs were enriched in integrin complex, external side of plasma membrane, and cell surface. For MF, 42 DEG-SIGs were enriched in cell part, tertiary granule, and whole membrane. These results may partially illustrate the potential mechanisms of 42 DEG-SIGs affecting the prognosis of AML patients.

19 hub genes were screened by integrated analysis of mRNA expression data, miRNA expression data and methylation data

Considering the complex mechanism of leukemogenesis and progression, we next conducted integrated multiomics analysis to identify hub genes that were associated with prognosis. Comparing the miRNA expression profiles of patients between IL and IR subtypes, we revealed 93 miRNAs that were significantly differentially expressed (P<0.05, |FC|≥2) (Figure 4A). A total of 7294 target miRNA genes (TDEmiRs) were identified using DIANO TOOLS/microT-CDS (threshold=0.9). Through integrated bioinformatics analysis, we selected 15 commonly differentially expressed genes from 42 DEG-SIGs and 7294 TDEmiRs between the IL and IR subtypes (Figure 4C).

FIGURE 4

Figure 4 Multiomics analysis of 97 TCGA-AML patients. (A) Volcano plot of differentially expressed miRNAs between IL and IR types. (B) Correlation between mRNA expression and DNA methylation level of 6 DEG & MethylCor genes. (C) Venn plot of the intersection of DEGs, SIGs, targets of DEmiRs and MethylCor gene set.

Combined analysis of mRNA and methylation profiles indicated that there were significantly negative correlations between mRNA expression level and degree of methylation for 355 genes (R < - 0.5, p< 0.05). When these 355 methylation correlation genes (MethylCor) were cross-referenced with the 42 DEG-SIGs, we identified 6 common genes associated with immune infiltration and differential expression, methylation and prognosis between IL and IR subtypes (Figures 4B, C).

A prognostic model based on 5 hub genes was constructed

Having observed significant differences in immune infiltration, gene expression and clinical behavior between IL and IR types, we next developed a LASSO-Cox proportional hazards regression model based on 19 immune-correlated DEGs by combining microRNA and epigenetic regulation data. Using the LASSO model, we built a classifier based on the 5 hub genes to predict the prognosis of AML (risk score = -0.086×ADAMTS3 + 0.180×CD52 + 0.472×CLCN5 - 0.356×HAL + 0.368×ICAM3) (Figures 5A, B). Kaplan-Meier plots displayed OS differences between patients in various subtypes (P=3.931×10^-06) (Figure 5C), and the ROC curve suggested that the model can effectively predict the 1-, 3- and 5-year prognosis of AML (AUC=0.82, 0.83, 0.99, respectively) (Figure 5D). Consistent with earlier analysis, we found similar predictive performance for 151 mRNA samples of the TCGA-AML profile (P=3.369×10^-06, AUC=0.63, 0.74, 0.83) (Figures 5E, F).

FIGURE 5

Figure 5 Construction of the COX regression model. (A) LASSO coefficient profiles of 19 candidate SIGs. (B) Tuning parameter (λ) selection cross‐validation error curve. The vertical lines were drawn at the optimal values determined by the minimum criteria and the 1‐SE criteria. (C, E, G) OS in patients with high vs. low risk scores depicted by Kaplan-Meier plots in the TCGA-AML-97, TCGA-AML-151 and GSE106291 cohorts. (D, F, H) ROC curves depicting the accuracy of the Cox regression model in identifying AML subtypes with poor prognosis in the TCGA-AML-97, TCGA-AML-151 and GSE106291 cohorts.

To test this model further, validation cohorts were obtained from the GEO database. Kaplan-Meier plots and ROC curves at 1, 3 and 5 years confirmed the prognostic value of the 5-hub-gene-based model: GSE106291 (log-rank test, P=3.311×10^-06, AUC=0.66, 0.66, 0.61) (Figures 5G, H). After stratification by disease classification, the results showed that the risk score of the IL type was significantly lower than that of the IR type (P=4.39×10^-07) (Figure S3). These evaluations demonstrated that the 5-hub-gene-based model could identify a group of high-risk patients within conventionally assigned risk groups and may guide clinical practice.

The better efficacy of the prognostic model for the prognosis of AML patients was further verified by clinical samples

For verifying the prognostic value in the 5-hub-gene-based model, we collected 6575 genes mutation detected in 200 newly diagnosed AML patients (TCGA.LAML.mutect.somatic.maf, https://portal.gdc.cancer.gov/files/27f42413-6d8f-401f-9d07-d019def8939e) and 38 genes mutation detected in 55 newly diagnosed AML patients (Xinqiao Hostpital). The common mutated genes were DNMT3A, IDH1, NRAS, RUNX1 and TET2. In this model classification, the high risk was significantly associated with the mutation of RUNX1 (p=0.015) and TET2 (p=0.054) considered by chi-square test (Table 2). Kaplan-Meier analyses of 55 patients with prognostic information indicated that patients with mutation of RUNX1 (Figure 6A, p=0.0001) and TET2 (Figure 6B, p=0.2257) were correlated with a poor prognosis and had a shorter median survival duration.

TABLE 2

Table 2 Results of Chi-square test to 5 common mutated genes based on regrouping LASSO model in TCGA-AML database.

FIGURE 6

Figure 6 Analysis of hub genes and mutated genes in AML. (A) Kaplan-Meier estimates of the OS according to RUNX1 mutation status. (B) Kaplan-Meier estimates of the OS according to TET2 mutation status. (wt, wild type; mut, mutation).

The diverse distribution of hub genes in immune cells of AML patients

To explore the value of these five hub genes in AML pathogenesis, we further identified the single-cell sequencing dataset GSE116256 to describe the distribution of the 5 hub genes in immune cells using the Seurat package for clustering and the SingleR package for annotation (Figure 7A). As shown in the scatter plot (Figure 7B) and violin plot (Figure 7C), CD52, ICAM3 and CLCN5 were widely expressed in granulocytes, monocytes, T lymphocytes, B lymphocytes, dendritic cells and NK cells, whereas ADAMTS3 was rarely expressed in those cells. HAL is highly expressed in granulocytes and monocytes but rarely expressed in other immune cells. Accordingly, we hypothesized that these hub genes play various roles through the regulation of gene expression in specific cells. The hub gene expression of blood cells in the Protein Atlas database (https://www.proteinatlas.org/) further confirmed this result (Figure S4).

FIGURE 7

Figure 7 scRNA analysis of hub genes. (A) Clustering analysis of the UMAP plot, color coded based on cell types. (B) Overlaying gene expression on UMAP clusters to illustrate the distribution of hub genes in each cell type. (C) Violin plots of hub genes in each cell type.

Investigation of best-fitting compounds on hub genes

To investigate best-fitting compounds, we performed virtual screening of molecular docking using the three-dimensional structure of CD52 (PDB ID: 6OBD), CLCN5(PDB ID: 2J9L), ICAM3(PDB ID: 1T0P) and 2115 FDA approved compounds in ZINC15 database. The predicted binding affinities of the top 2 hit compounds against respective targets are ranked from highest to lowest. Binding energy (Kcal/mol) for interaction of proteins and compounds are as follows: CD52 with ZINC164528615 (Glecaprevir) (−6.4 Kcal/mol), CD52 with ZINC3938684 (Toposar) (−6.3 Kcal/mol); ICAM3 with ZINC52955754 (Ergotamine) (−8.3 Kcal/mol), ICAM3 with ZINC1612996 (Irinotecan) (−8.2 Kcal/mol); CLCN5 with ZINC3978005 (Dihydroergotamine) (−11.8 Kcal/mol), CLCN5 with ZINC52955754 (Ergotamine) (−11.5 Kcal/mol). 2D visualization of the most probable interactions of these proteins and candidate compounds are represented in Figure 8.

FIGURE 8

Figure 8 Shows 2D interaction representations of the best pose of (A) CD52 with ZINC164528615 (Glecaprevir), (B) CD52 with ZINC3938684 (Toposar); (C) ICAM3 with ZINC52955754 (Ergotamine), (D) ICAM3 with ZINC1612996 (Irinotecan); (E) CLCN5 with ZINC3978005 (Dihydroergotamine), (F) CLCN5 with ZINC52955754 (Ergotamine).

Discussion

AML is a highly heterogeneous malignant tumor with poor prognosis. A large number of studies have demonstrated that the immune microenvironment of AML patients is altered significantly to promote leukemogenesis in AML (35, 36). In this study, we classified AML patients according to the difference in the SIG expression signature. The results showed that the prognosis of AML patients with immune infiltration deficiency (IL subtype) was significantly better than that of patients with immune infiltration enrichment (IR subtype), which was contrary to the effect of the immune infiltration degree in solid tumors (37). Recent studies have shown that the heterogeneity of immune infiltration may be affected by the components of cytokines in the microenvironment (38, 39) and the expression of tumor driver genes (40). These heterogeneities were significantly associated with immunotherapeutic effects.

We further compared the gene expression profiles of IL and IR patients and found that the functions of genes with differential expression between the two subtypes were mainly enriched in defense response, inflammatory response and immune system process (BP); integrin complex, plasma membrane and cell surface (CC); and cellular partial, tertiary granules, and whole membrane (MF). Consistent with our results, several studies have confirmed that the biological functional diversity of the immune system is significantly altered in the BM immune microenvironment.

Previous studies have constructed many prognostic models on the basis of diverse omics data. Hu et al. constructed a DNA methylation-based prognostic model for AML patients (41), but the accuracy of the model (AUC=0.67-0.75) was lower than that of our model (AUC=0.82-0.99). A three-miRNA signature was built for non-M3 AML patients by Xue et al., but due to the lack of validation by any independent data, its clinical application value may need further verification (42). A four-gene-based prognostic model from Huang et al. (43) was also less accurate (AUC=0.66-0.71) than our model.

These conventional prognostic analyses generally consider only mRNA, miRNA, or methylation data. To the best of our knowledge, this is the only study to date to construct a prognostic model through the integrated analysis of multiple omics datasets (mRNA, miRNA and methylation data). The TCGA-AML test dataset and an independent GEO-AML validation dataset confirm that the model is very effective in predicting the prognosis of AML patients. Further, the prognostic model was verified by mutation data of AML patients in Xinqiao Hospital. Although the data were limited, it also indirectly proved that the model could effectively predict the prognosis of AML patients. The model consists of five hub genes: ADAMTS3, CD52, CLCN5, ICAM3 and HAL.

ADAMTS3 is a member of the ADAMTS family and plays an important role in the genesis and development of a variety of tumors (44, 45). Previous studies demonstrated that ADAMTS3 is expressed in mast cells (19), and mast cells have both tumor-promoting and tumor-inhibiting effects (46), thus leading to different prognoses in different tumors (47–49). Our study confirmed that ADAMTS3 expression was significantly negatively correlated with the prognosis of AML and is rarely expressed in other immune cells. Therefore, we will pursue further in-depth analysis of how ADAMTS3 exerts its antitumor effect through mast cells to explore its potential therapeutic value.

CLCN5 (chloride voltage gated channel 5) is a member of the chloride channel family. It is widely expressed in a variety of tumor cells (50) and can enhance the chemotherapy resistance of chronic lymphocytic leukemia and multiple myeloma (51, 52). Hsa-let-7c-5p and hsa-mir-495-3p, two miRNAs that target CLCN5, were found to be involved in the occurrence and development of a variety of tumors (53, 54). Our study confirmed that CLCN5 plays a role in promoting tumorigenesis in AML, but its specific mechanism needs deeper investigation.

HAL (histidine ammonia lyase) is the rate-limiting enzyme of histidine catabolism (55). Kanarek et al. found that tumor cells with higher HAL expression levels possess higher sensitivity to methotrexate. Acute lymphocyte leukemia (ALL) patients with higher HAL expression were more likely to benefit from methotrexate treatment (56). In our study, we found that the expression level of HAL was significantly correlated with the degree of gene methylation. However, as a demethylation agent, hypomethylating agent (HMA) alone has difficulty achieving the intended effect in the clinical treatment of AML (57). We speculated that the combination of methotrexate and HMA may achieve better efficacy when changes in HAL expression are detected. Surprisingly, we found that HAL was only expressed in granulocytes and monocytes. Consistent with AML cells, granulocytes and monocytes originate from myeloid precursor cells (14). Therefore, insight into the methylation mechanism of HAL and its effect on the innate immune response and AML cells will be conducive to improving the therapeutic effect of demethylation agents. We also assessed the interaction between HAL and hsa-miR-582-3p, which has been previously indicated as a tumor suppressor miRNA in AML (58). Thus, deeper investigation of the relationship between HAL and miR-582-3p will be helpful to further understand the potential mechanism of malignant progression of AML.

Our study found that CD52 was widely expressed in leukemia cells and all types of immune cells. Evidence has revealed the high expression of CD52 in CD34+ stem cells of AML (5q-) patients and a significantly negative correlation between the CD52 expression value and the prognosis of AML (59), which is consistent with our results. ICAM-3 (intercellular adhesion molecule 3, CD50) belongs to the ICAM (intercellular adhesion molecule) immunoglobulin superfamily and plays important roles in immune response and tumor development (60, 61). Consistent with our scRNA sequencing results, ICAM-3 is expressed in granulocytes, monocytes and lymphocytes (62). In vivo and in vitro experiments confirmed that ICAM-3 participates in the proliferation, stemness and radiotherapy resistance of various tumors through the FAK pathway, PI3K/Akt pathway or other mechanisms (63–65). Interestingly, we found that the expression levels of CD52 and ICAM-3 were both significantly correlated with the degree of methylation but had opposite impacts on the prognosis of AML patients. Therefore, the roles of CD52 and ICAM-3 should be considered in HMA reagent treatment.

In conclusion, using a multiomics analysis and validation approach, we constructed and validated a novel, 5-hub-gene-based model that allows robust risk stratification and facilitates the identification of prognosis in AML. The distribution of the 5 hub genes in immune cells was revealed through scRNA sequencing analysis. Furthermore, we conducted virtual screening of three genes (CD52, CLCN5 and ICAM3) with known protein structure, and found the compounds with the lowest binding energy with them, which provided ideas for further searching for targeted inhibitors. Of course, this study is limited to bioinformatics analysis, and the proposed approaches need to be further tested in the clinic.

Data availability statement

The original contributions presented in the study are included in the article/Supplementary Material. Further inquiries can be directed to the corresponding authors.

Author contributions

YS and LG conceived and designed the experiments. FL and JC performed experiments. FL analyzed data. S-CY and XZ provided bioinformatics support. JL characterized patient samples. LG supervised the project. FL and JC wrote the manuscript. All authors contributed to the article and approved the submitted version.

Funding

This work was supported by the Chinese National Natural Science Foundation (Grant No. 82170161), Chongqing National Natural Science Key Foundation (Grant No. cstc2019jcyj-zdxmX0023), the Basic and Frontier Research Project of General Hospital of the Chinese People’s Liberation Army Western Theater (Grant No. 2021XZYG-C44).

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fonc.2022.925615/full#supplementary-material

References

1. Dohner H, Weisdorf DJ, Bloomfield CD. Acute myeloid leukemia. New Engl J Med (2015) 373(12):1136–52. doi: 10.1056/NEJMra1406184

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Short NJ, Rytting ME, Cortes JE. Acute myeloid leukaemia. Lancet (2018) 392(10147):593–606. doi: 10.1016/S0140-6736(18)31041-9

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Döhner H, Estey EH, Amadori S, Appelbaum FR, Büchner T, Burnett AK, et al. Diagnosis and management of acute myeloid leukemia in adults: Recommendations from an international expert panel, on behalf of the European leukemianet. Blood (2010) 115(3):453–74. doi: 10.1182/blood-2009-07-235358

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Ghobrial IM, Detappe A, Anderson KC, Steensma DP. The bone-marrow niche in mds and mgus: Implications for aml and mm. Nat Rev Clin Oncol (2018) 15(4):219–33. doi: 10.1038/nrclinonc.2017.197

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Cogle CR, Bosse RC, Brewer T, Migdady Y, Shirzad R, Kampen KR, et al. Acute myeloid leukemia in the vascular niche. Cancer Lett (2016) 380(2):552–60. doi: 10.1016/j.canlet.2015.05.007

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Padró T, Bieker R, Ruiz S, Steins M, Retzlaff S, Bürger H, et al. Overexpression of vascular endothelial growth factor (Vegf) and its cellular receptor kdr (Vegfr-2) in the bone marrow of patients with acute myeloid leukemia. Leukemia (2002) 16(7):1302–10. doi: 10.1038/sj.leu.2402534

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Aguayo A, Estey E, Kantarjian H, Mansouri T, Gidel C, Keating M, et al. Cellular vascular endothelial growth factor is a predictor of outcome in patients with acute myeloid leukemia. Blood (1999) 94(11):3717–21. doi: 10.1182/blood.V94.11.3717

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Wegiel B, Ekberg J, Talasila KM, Jalili S, Persson JL. The role of vegf and a functional link between vegf and P27kip1 in acute myeloid leukemia. Leukemia (2009) 23(2):251–61. doi: 10.1038/leu.2008.300

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Hazlehurst LA, Dalton WS. Mechanisms associated with cell adhesion mediated drug resistance (Cam-Dr) in hematopoietic malignancies. Cancer Metastasis Rev (2001) 20(1-2):43–50. doi: 10.1023/a:1013156407224

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Tabe Y, Jin L, Tsutsumi-Ishii Y, Xu Y, McQueen T, Priebe W, et al. Activation of integrin-linked kinase is a critical prosurvival pathway induced in leukemic cells by bone marrow-derived stromal cells. Cancer Res (2007) 67(2):684–94. doi: 10.1158/0008-5472.Can-06-3166

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Jacamo R, Chen Y, Wang Z, Ma W, Zhang M, Spaeth EL, et al. Reciprocal leukemia-stroma vcam-1/Vla-4-Dependent activation of nf-κb mediates chemoresistance. Blood (2014) 123(17):2691–702. doi: 10.1182/blood-2013-06-511527

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Shafat MS, Oellerich T, Mohr S, Robinson SD, Edwards DR, Marlein CR, et al. Leukemic blasts program bone marrow adipocytes to generate a protumoral microenvironment. Blood (2017) 129(10):1320–32. doi: 10.1182/blood-2016-08-734798

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Ye H, Adane B, Khan N, Sullivan T, Minhajuddin M, Gasparetto M, et al. Leukemic stem cells evade chemotherapy by metabolic adaptation to an adipose tissue niche. Cell Stem Cell (2016) 19(1):23–37. doi: 10.1016/j.stem.2016.06.001

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Mendez LM, Posey RR, Pandolfi PP. The interplay between the genetic and immune landscapes of aml: Mechanisms and implications for risk stratification and therapy. Front Oncol (2019) 9:1162. doi: 10.3389/fonc.2019.01162

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Vadakekolathu J, Minden MD, Hood T, Church SE, Reeder S, Altmann H, et al. Immune landscapes predict chemotherapy resistance and immunotherapy response in acute myeloid leukemia. Sci Trans Med (2020) 12(546):eaaz0463. doi: 10.1126/scitranslmed.aaz0463

CrossRef Full Text | Google Scholar

16. Yan H, Qu J, Cao W, Liu Y, Zheng G, Zhang E, et al. Identification of prognostic genes in the acute myeloid leukemia immune microenvironment based on tcga data analysis. Cancer Immunol Immunother (2019) 68(12):1971–8. doi: 10.1007/s00262-019-02408-7

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Allie SR, Zhang W, Tsai CY, Noelle RJ, Usherwood EJ. Critical role for all-trans retinoic acid for optimal effector and effector memory Cd8 T cell differentiation. J Immunol (Baltimore Md 1950) (2013) 190(5):2178–87. doi: 10.4049/jimmunol.1201945

CrossRef Full Text | Google Scholar

18. Assi R, Kantarjian H, Ravandi F, Daver N. Immune therapies in acute myeloid leukemia: A focus on monoclonal antibodies and immune checkpoint inhibitors. Curr Opin Hematol (2018) 25(2):136–45. doi: 10.1097/moh.0000000000000401

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Charoentong P, Finotello F, Angelova M, Mayer C, Efremova M, Rieder D, et al. Pan-cancer immunogenomic analyses reveal genotype-immunophenotype relationships and predictors of response to checkpoint blockade. Cell Rep (2017) 18(1):248–62. doi: 10.1016/j.celrep.2016.12.019

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for rna-seq data with Deseq2. Genome Biol (2014) 15(12):550. doi: 10.1186/s13059-014-0550-8

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Gevaert O. Methylmix: An r package for identifying DNA methylation-driven genes. Bioinformatics (2015) 31(11):1839–41. doi: 10.1093/bioinformatics/btv020

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Subramanian A, Tamayo P, Mootha VK, Mukherjee S, Ebert BL, Gillette MA, et al. Gene set enrichment analysis: A knowledge-based approach for interpreting genome-wide expression profiles. Proc Natl Acad Sci USA (2005) 102(43):15545–50. doi: 10.1073/pnas.0506580102

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Yoshihara K, Shahmoradgoli M, Martínez E, Vegesna R, Kim H, Torres-Garcia W, et al. Inferring tumour purity and stromal and immune cell admixture from expression data. Nat Commun (2013) 4:2612. doi: 10.1038/ncomms3612

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Szklarczyk D, Gable AL, Lyon D, Junge A, Wyder S, Huerta-Cepas J, et al. String V11: Protein-protein association networks with increased coverage, supporting functional discovery in genome-wide experimental datasets. Nucleic Acids Res (2019) 47(D1):D607–13. doi: 10.1093/nar/gky1131

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Shannon P, Markiel A, Ozier O, Baliga NS, Wang JT, Ramage D, et al. Cytoscape: A software environment for integrated models of biomolecular interaction networks. Genome Res (2003) 13(11):2498–504. doi: 10.1101/gr.1239303

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Huang da W, Sherman BT, Lempicki RA. Systematic and integrative analysis of Large gene lists using David bioinformatics resources. Nat Protoc (2009) 4(1):44–57. doi: 10.1038/nprot.2008.211

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Huang da W, Sherman BT, Lempicki RA. Bioinformatics enrichment tools: Paths toward the comprehensive functional analysis of Large gene lists. Nucleic Acids Res (2009) 37(1):1–13. doi: 10.1093/nar/gkn923

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Walter W, Sánchez-Cabo F, Ricote M. Goplot: An r package for visually combining expression data with functional analysis. Bioinformatics (2015) 31(17):2912–4. doi: 10.1093/bioinformatics/btv300

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Simon N, Friedman J, Hastie T, Tibshirani R. Regularization paths for cox's proportional hazards model Via coordinate descent. J Stat Software (2011) 39(5):1–13. doi: 10.18637/jss.v039.i05

CrossRef Full Text | Google Scholar

30. Friedman J, Hastie T, Tibshirani R. Regularization paths for generalized linear models Via coordinate descent. J Stat software (2010) 33(1):1–22. doi: 10.18637/jss.v033.i01

CrossRef Full Text | Google Scholar

31. Butler A, Hoffman P, Smibert P, Papalexi E, Satija R. Integrating single-cell transcriptomic data across different conditions, technologies, and species. Nat Biotechnol (2018) 36(5):411–20. doi: 10.1038/nbt.4096

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Stuart T, Butler A, Hoffman P, Hafemeister C, Papalexi E, Mauck WM III, et al. Comprehensive integration of single-cell data. Cell (2019) 177(7):1888–902.e21. doi: 10.1016/j.cell.2019.05.031

PubMed Abstract | CrossRef Full Text | Google Scholar

33. Aran D, Looney AP, Liu L, Wu E, Fong V, Hsu A, et al. Reference-based analysis of lung single-cell sequencing reveals a transitional profibrotic macrophage. Nat Immunol (2019) 20(2):163–72. doi: 10.1038/s41590-018-0276-y

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Trott O, Olson AJ. Autodock vina: Improving the speed and accuracy of docking with a new scoring function, efficient optimization, and multithreading. J Comput Chem (2010) 31(2):455–61. doi: 10.1002/jcc.21334

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Mittal D, Gubin MM, Schreiber RD, Smyth MJ. New insights into cancer immunoediting and its three component phases–elimination, equilibrium and escape. Curr Opin Immunol (2014) 27:16–25. doi: 10.1016/j.coi.2014.01.004

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Teague RM, Kline J. Immune evasion in acute myeloid leukemia: Current concepts and future directions. J immunotherapy Cancer (2013) 1(13):1. doi: 10.1186/2051-1426-1-13

CrossRef Full Text | Google Scholar

37. Thorsson V, Gibbs DL, Brown SD, Wolf D, Bortone DS, Ou Yang TH, et al. The immune landscape of cancer. Immunity (2018) 48(4):812–30 e14. doi: 10.1016/j.immuni.2018.03.023

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Benci JL, Xu B, Qiu Y, Wu TJ, Dada H, Twyman-Saint Victor C, et al. Tumor interferon signaling regulates a multigenic resistance program to immune checkpoint blockade. Cell (2016) 167(6):1540–54.e12. doi: 10.1016/j.cell.2016.11.022

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Li J, Byrne KT, Yan F, Yamazoe T, Chen Z, Baslan T, et al. Tumor cell-intrinsic factors underlie heterogeneity of immune cell infiltration and response to immunotherapy. Immunity (2018) 49(1):178–93.e7. doi: 10.1016/j.immuni.2018.06.006

PubMed Abstract | CrossRef Full Text | Google Scholar

40. Bezzi M, Seitzer N, Ishikawa T, Reschke M, Chen M, Wang G, et al. Diverse genetic-driven immune landscapes dictate tumor progression through distinct mechanisms. Nat Med (2018) 24(2):165–75. doi: 10.1038/nm.4463

PubMed Abstract | CrossRef Full Text | Google Scholar

41. Hu L, Gao Y, Shi Z, Liu Y, Zhao J, Xiao Z, et al. DNA Methylation-based prognostic biomarkers of acute myeloid leukemia patients. Ann Trans Med (2019) 7(23):737. doi: 10.21037/atm.2019.11.122

CrossRef Full Text | Google Scholar

42. Xue Y, Ge Y, Kang M, Wu C, Wang Y, Rong L, et al. Selection of three mirna signatures with prognostic value in non-M3 acute myeloid leukemia. BMC Cancer (2019) 19(1):109–. doi: 10.1186/s12885-019-5315-z

PubMed Abstract | CrossRef Full Text | Google Scholar

43. Huang R, Liao X, Li Q. Identification and validation of potential prognostic gene biomarkers for predicting survival in patients with acute myeloid leukemia. OncoTargets Ther (2017) 10:5243–54. doi: 10.2147/OTT.S147717

CrossRef Full Text | Google Scholar

44. Rocks N, Paulissen G, El Hour M, Quesada F, Crahay C, Gueders M, et al. Emerging roles of Adam and adamts metalloproteinases in cancer. Biochimie (2008) 90(2):369–79. doi: 10.1016/j.biochi.2007.08.008

PubMed Abstract | CrossRef Full Text | Google Scholar

45. Gomis-Rüth FX. Catalytic domain architecture of metzincin metalloproteases. J Biol Chem (2009) 284(23):15353–7. doi: 10.1074/jbc.R800069200

PubMed Abstract | CrossRef Full Text | Google Scholar

46. Oldford SA, Marshall JS. Mast cells as targets for immunotherapy of solid tumors. Mol Immunol (2015) 63(1):113–24. doi: 10.1016/j.molimm.2014.02.020

PubMed Abstract | CrossRef Full Text | Google Scholar

47. Rajput AB, Turbin DA, Cheang MC, Voduc DK, Leung S, Gelmon KA, et al. Stromal mast cells in invasive breast cancer are a marker of favourable prognosis: A study of 4,444 cases. Breast Cancer Res Treat (2008) 107(2):249–57. doi: 10.1007/s10549-007-9546-3

PubMed Abstract | CrossRef Full Text | Google Scholar

48. Fleischmann A, Schlomm T, Köllermann J, Sekulic N, Huland H, Mirlacher M, et al. Immunological microenvironment in prostate cancer: High mast cell densities are associated with favorable tumor characteristics and good prognosis. Prostate (2009) 69(9):976–81. doi: 10.1002/pros.20948

PubMed Abstract | CrossRef Full Text | Google Scholar

49. Groot Kormelink T, Abudukelimu A, Redegeld FA. Mast cells as target in cancer therapy. Curr Pharm Design (2009) 15(16):1868–78. doi: 10.2174/138161209788453284

CrossRef Full Text | Google Scholar

50. Ernest NJ, Weaver AK, Van Duyn LB, Sontheimer HW. Relative contribution of chloride channels and transporters to regulatory volume decrease in human glioma cells. Am J Physiol Cell Physiol (2005) 288(6):C1451–60. doi: 10.1152/ajpcell.00503.2004

PubMed Abstract | CrossRef Full Text | Google Scholar

51. Zhang H, Pang Y, Ma C, Li J, Wang H, Shao Z. Clc5 decreases the sensitivity of multiple myeloma cells to bortezomib Via promoting prosurvival autophagy. Oncol Res (2018) 26(3):421–9. doi: 10.3727/096504017X15049221237147

PubMed Abstract | CrossRef Full Text | Google Scholar

52. Ruiz-Lafuente N, Alcaraz-Garcia MJ, Sebastian-Ruiz S, Garcia-Serna AM, Gomez-Espuch J, Moraleda JM, et al. Il-4 up-regulates mir-21 and the mirnas hosted in the Clcn5 gene in chronic lymphocytic leukemia. PloS One (2015) 10(4):e0124936. doi: 10.1371/journal.pone.0124936

PubMed Abstract | CrossRef Full Text | Google Scholar

53. Wu Q, Liu P, Lao G, Liu Y, Zhang W, Ma C. Comprehensive analysis of circrna-Mirna-Mrna network in cervical squamous cell carcinoma by integrated analysis. OncoTargets Ther (2020) 13:8641–50. doi: 10.2147/ott.S254323

CrossRef Full Text | Google Scholar

54. Liao K, Qian Z, Zhang S, Chen B, Li Z, Huang R, et al. The lgmn pseudogene promotes tumor progression by acting as a mir-495-3p sponge in glioblastoma. Cancer Lett (2020) 490:111–23. doi: 10.1016/j.canlet.2020.07.012

PubMed Abstract | CrossRef Full Text | Google Scholar

55. Brand Lm Fau - Harper AE, Harper AE. Studies on the production and assessment of experimental histidinemia in the rat. Biochim Biophys Acta (1976) 444(1):294–306. doi: 10.1016/0304-4165(76)90246-4

PubMed Abstract | CrossRef Full Text | Google Scholar

56. Kanarek N, Keys HR, Cantor JR, Lewis CA, Chan SH, Kunchok T, et al. Histidine catabolism is a major determinant of methotrexate sensitivity. Nature (2018) 559(7715):632–6. doi: 10.1038/s41586-018-0316-7

PubMed Abstract | CrossRef Full Text | Google Scholar

57. Duchmann M, Itzykson R. Clinical update on hypomethylating agents. Int J Hematol (2019) 110(2):161–9. doi: 10.1007/s12185-019-02651-9

PubMed Abstract | CrossRef Full Text | Google Scholar

58. Li H, Tian X, Wang P, Huang M, Xu R, Nie T. Microrna-582-3p negatively regulates cell proliferation and cell cycle progression in acute myeloid leukemia by targeting cyclin B2. Cell Mol Biol Lett (2019) 24:66. doi: 10.1186/s11658-019-0184-7

PubMed Abstract | CrossRef Full Text | Google Scholar

59. Blatt K, Herrmann H, Hoermann G, Willmann M, Cerny-Reiterer S, Sadovnik I, et al. Identification of campath-1 (Cd52) as novel drug target in neoplastic stem cells in 5q-patients with mds and aml. Clin Cancer Res (2014) 20(13):3589–602. doi: 10.1158/1078-0432.Ccr-13-2811

PubMed Abstract | CrossRef Full Text | Google Scholar

60. Cavallaro U, Christofori G. Cell adhesion and signalling by cadherins and ig-cams in cancer. Nat Rev Cancer (2004) 4(2):118–32. doi: 10.1038/nrc1276

PubMed Abstract | CrossRef Full Text | Google Scholar

61. Xiao X, Mruk DD, Cheng CY. Intercellular adhesion molecules (Icams) and spermatogenesis. Hum Reprod Update (2013) 19(2):167–86. doi: 10.1093/humupd/dms049

PubMed Abstract | CrossRef Full Text | Google Scholar

62. de Fougerolles AR, Springer TA. Intercellular adhesion molecule 3, a third adhesion counter-receptor for lymphocyte function-associated molecule 1 on resting lymphocytes. J Exp Med (1992) 175(1):185–90. doi: 10.1084/jem.175.1.185

PubMed Abstract | CrossRef Full Text | Google Scholar

63. Chung YM, Kim BG, Park CS, Huh SJ, Kim J, Park JK, et al. Increased expression of icam-3 is associated with radiation resistance in cervical cancer. Int J Cancer (2005) 117(2):194–201. doi: 10.1002/ijc.21180

PubMed Abstract | CrossRef Full Text | Google Scholar

64. Kim YG, Kim MJ, Lim JS, Lee MS, Kim JS, Yoo YD. Icam-3-Induced cancer cell proliferation through the Pi3k/Akt pathway. Cancer Lett (2006) 239(1):103–10. doi: 10.1016/j.canlet.2005.07.023

PubMed Abstract | CrossRef Full Text | Google Scholar

65. Shen W, Xie J, Zhao S, Du R, Luo X, He H. Icam3 Mediates Inflammatory Signaling to Promote Cancer Cell Stemness. Cancer Letters (2018) 422:29–43. doi: 10.1016/j.canlet.2018.02.034

PubMed Abstract | CrossRef Full Text | Google Scholar