Identification and validation of inferior prognostic genes associated with immune signatures and chemotherapy outcome in acute myeloid leukemia

Identification of differentially expressed genes related to acute myeloid leukemia

This study obtained an expression profile data matrix that includes 67758 genes in 214 samples, the R software to process the GSE114868 transcriptome data. The workflow of the analysis procedure in this study was shown in Figure 1. Based on the screening criteria as |FC| > 2 and P < 0.05, 3893 differentially expressed genes (DEGs) were identified in total. These DEGs include 1657 up-regulated genes and 2236 down-regulated genes. We used R software to draw the volcano plot and heatmap of differentially expressed genes (Figure 2).

Figure 1. The overall analysis process of the present study.

Figure 2. The volcano and heatmap of the differentially expressed genes. (A) In the DEGs' volcano plot, red dots are up-regulated genes, blue is down-regulated, and grey is no different; (B) DEGs clustering heatmap, pink indicates healthy, and blue indicates AML. DEGs: differentially expressed genes.

Weighted co-expression network construction and identification of the critical module

To find the co-expressed gene modules (Module), explore the relationship between the gene network and the phenotype of interest, and the core genes in the network, we used the WGCNA method to construct a weighted gene co-expression network. In the process of network construction, we selected the soft threshold β = 8 (R² = 0.85) to establish a topological matrix whose gene distribution conforms to the scale-free network. The topological matrix was clustered with co-expressed genes through the coefficient of dissimilarity between genes. The gene co-expression network was divided into nine modules, and the topological overlap of adjacent modules was used to select characteristic functional modules. Finally, it determined that the Blue module has the highest module importance and containing 75 genes (Figure 3).

Figure 3. Identification of modules related to clinical features of acute myeloid leukemia. (A) Analysis of the soft threshold (β) through the scale-free fitting index and mean connectivity; (B) Clustering dendrogram of the DEGs through dissimilarity coefficient, which shows nine gene co-expression modules AML. Gray modules indicate no co-expression between genes; (C) The correlation heat map of WGCNA adjacent modules. The rectangles in each row and each column represent a module characteristic gene. Light blue represents low adjacency, and red represents high adjacency; (D) The TOM visualized the gene co-expression network's heat map in the module. In the TOM map, light colors indicate topological overlap. Dark colors indicate a higher degree of topological overlap. The gene tree diagram and corresponding modules are displayed on the upper left of the TOM diagram. The intersection of the two rectangles indicates the topological overlap in the Blue module. DEGs: differentially expressed genes; TOM: topological overlap matrix.

Gene ontology (GO) functional annotation and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis

We performed the GO functional annotations on 75 characteristic genes in the Blue module, including biological process (BP), Cellular Component (CC), and Molecular Function (MF) analysis. BP analysis showed that the co-expressed genes were mainly annotated in neutrophil degranulation, neutrophil activation, and neutrophils participating in the immune activation response. CC analysis showed that the co-expressed genes were significantly annotated in the cell components such as secretory granule membrane, ficolin granule, secretory granule lumen, and cytoplasmic cyst. MF analysis shows that the co-expressed genes are mainly related to carbon-nitrogen bond hydrolase activity, SH3-domain binding, and RAGE receptor binding. KEGG pathway enrichment analysis results show that the co-expressed genes mainly involve neutrophil migration, tumor carbon metabolism, iron death, FcγR-mediated phagocytosis, HIF-1, and NF-κB signaling pathway. These results indicate that the co-expressed genes in the necessary modules are mainly involved in cellular immune and inflammatory processes (Figure 4).

Figure 4. GO functional enrichment analysis and KEGG pathway analysis of characteristic genes in the Blue module. (A) Top 20 enriched biological process terms; (B) Top 20 enriched cell component enrichment; (C) The results of Molecular function enrichment analysis; (D) The results of the KEGG pathway enrichment analysis. Abbreviations: GO: gene ontology; KEGG: Kyoto Encyclopedia of Genes and Genomes.

Recognition of module driver genes based on PPI network

To further study the function of the characteristic genes in the Blue module at the protein level after the WGCNA analysis, we used the STRING database (https://string-db.org/) to screen 75 of the Blue modules co-expression of genes constructs a protein-protein interaction (PPI) network (Figure 5A).

Figure 5. The module driver genes identified in the PPI network. (A) The PPI network consists of 45 nodes and 66 edges, and 75 of the Blue module's co-expression of genes constructs a PPI network; (B) Clusters of driver genes in the Blue module, and the squares marked in yellow to red indicate the top 10 module driver genes in sub-function cluster 1, and the blue squares represented the other related genes in cluster 1. Abbreviation: PPI: protein-protein interaction.

The PPI network consists of 45 nodes and 66 edges, and we used the "MCODE" program in Cytoscape software to further identify the module driver genes (MDGs) in the PPI network. The cut-off value is set as K-core = 2. Finally, two functional clusters were screened out. Cluster 1 contains 10 MDGs (S100A9, S100A8, NCF2, CYBB, ITGAM, FCER1G, FPR1, FCN1, HK3, VNN2) (Figure 5B), and cluster 2 includes 4 MDGs (HCAR3, HCAR2, C5AR1, P2RY13) (Table 1).

Identification of core differentially expressed mRNAs (DEmRNAs) in the whole transcriptome

To further verify the regulatory function of the obtained MDGs in the whole transcriptome, we also identified 745 DEmRNAs in the transcriptome data with more stringent screening conditions (|logFC ≥ 2|, P < 0.05), and using the STRING database constructed a PPI network, which contains 392 nodes and 2357 edges (Figure 6A). Besides, we also use the "cytoHubba" program in Cytoscape software to calculate the degree of contribution (Degree) in the PPI network, then obtain the ten core DEmRNAs (FPR2, PPBP, ITGB2, ITGAM, APP, CEACAM8, STOM, CKAP4, ORM2, FPR1) with the highest degree of contribution (Degree) in the global network (Figure 6B). Then the obtained top 10 DEmRNAs are compared with the MDGs in the module (Blue). The results show that ITGAM and FPR1 are common genes, and these two common genes show vital global regulatory functions. They may play an essential role in the pathogenesis of AML.

Figure 6. The core DEmRNAs identified in the PPI network. (A) The darker the color (red) of the genes, the higher degree of contribution in the PPI network; (B) The darker color (red) of the DEmRNA represents the gene with a higher degree. Abbreviations: DEmRNAs, differentially expressed lncRNAs. PPI: protein-protein interaction.

Screening and correlation analysis of differentially expressed lncRNAs (DElncRNAs)

To further determine the potential regulatory relationship between differentially expressed lncRNAs (DElncRNAs) and core DEmRNAs in AML and clarify the regulatory function of DElncRNAs. We screened 474 DElncRNAs from the transcriptome data and displayed the volcano plot of DElncRNAs by the R software (Figure 7A).

Figure 7. Identification of DElncRNAs and the co-relationship of core DEmRNAs and DElncRNAs. (A) In the volcano plot of the DElncRNAs, red dots are up-regulated lncRNAs, blue is down-regulated, and grey is no different; (B) The expression relevance of DEmRNAs and DElncRNAs. The darker the color (red) of the circle, the stronger the correlation. Abbreviations: DEmRNAs: differentially expressed mRNAs; DElncRNAs: differentially expressed lncRNAs.

After that, we used the PCC method to evaluate the correlation between core DEmRNAs and DElncRNAs, then screened the noteworthy DElncRNA-mRNA pairs (PCC > 0.8, P < 0.05). The results showed that the expression of 7 lncRNAs (SIRPAP1, ITGB2-AS1, FAM157B, FAM157A, CASP1P2, SMIM25, IL10RB-DT) significantly correlated with the expression of core mRNAs (Figure 7B). Among them, the expression of DElncRNAs includes FAM157A and SIRPAP1 that are mainly related to the expression of ITGAM. Furthermore, DElncRNAs as FAM157A, FAM157B, SMIM25, ITGB2-AS1, and SIRPAP1 are significantly associated with the expression of FPR1.

Prognostic evaluation of critical mRNAs and co-expressed lncRNA

To verify the prognostic value of the screened essential genes and co-expressed DElncRNAs, we conducted prognostic verification on the 18 important genes (including two common genes) and 7 lncRNAs selected based on the TCGA-LAML clinical data in the GEPIA online database (http://gepia2.cancer-pku.cn).

The univariate Cox regression analysis showed that both the critical genes as S100A9, S100A8, NCF2, ITGAM, HK3, VNN2, PPBP, ITGB2 (Figure 8A), and co-expressed DElncRNA as ITGB2-AS1 (Figure 8B) are the risk factors affecting the prognosis of AML patients (P < 0.05). The results of K-M survival analysis showed that AML patients with high expression of ITGAM, PPBP, and ITGB2-AS1 had a poor prognosis (P < 0.05). In contrast, the common gene expression as FPR1 did not significantly affect AML patients' prognosis. It is worth noting that, combining the results of the two prognostic analyses, we found that ITGAM, PPBP, and ITGB2-AS1 both significantly impact AML patients' prognosis (Figure 8C). Finally, The Nomogram (Figure 8D) shows the probability of these three genes influencing the prognostic outcome.

Figure 8. Prognostic analysis results of essential mRNAs and co-expressed DElncRNAs. (A) The result of univariate Cox regression analysis showed that the key genes such as S100A8 (HR:1.119, 95% CI:1.01–1.16), S100A9 (HR:1.08, 95% CI:1.00–1.16), NCF2 (HR:1.19, 95% CI:1.05–1.35), ITGAM (HR:1.19, 95% CI:1.05–1.36), HK3 (HR:1.09, 95% CI:1.01–1.08), VNN2 (HR:1.12, 95% CI:1.01–1.24), PPBP (HR:1.10, 95% CI:1.02–1.19), and ITGB2 (HR:1.33, 95% CI:1.14–1.56) both have significant impact on the prognosis of AML patients (P < 0.05). (B) The development of the research showed that the expression of DElnRNAs as ITGB2-AS1 (HR:1.24, 95% CI:1.10–1.41) has a significant impact on the prognosis of AML patients (P < 0.05). (C) The results of K–M survival analysis showed that AML patients with high expression of ITGAM, PPBP, and ITGB2-AS1 had a poor prognosis (P < 0.05). (D) The Nomogram was established based on the clinical information of TCGA-LAML. The points for 11 factors (gender, cytogenetics risk category, age, leukocyte, hemoglobin, monocyte, platelet, FAB classification, and the expression level of ITGAM, PPBP, or ITGB2-AS1) were listed in the Nomogram. The score for each factor in the Nomogram was read out by drawing a straight line from the predictor to the point axis, and then the survival rates of 1, 3, and 5 years could be estimated by adding the points corresponding to each factor in the bottom scale. Abbreviations: DElncRNAs: differentially expressed lncRNAs; HR: hazard ratio; CI: confidence interval.

Quantitative real-time PCR (qRT-PCR) validated the expression level of critical genes in initial diagnosis patients with AML

We used 20 clinical samples to verify the different expression levels of essential genes between AML (15 bone marrow samples from initial diagnosis of AML) and health group (15 peripheral blood samples from anonymized healthy individuals). The validation result showed that ITGB2-AS1,ITGAM, and PPBP were significantly higher expressed in the AML group (Figure 9).

Figure 9. The differential expression of critical genes in clinical samples between AML and healthy individuals. The result of clinical validation showed that ITGB2-AS, ITGAM, and PPBP are significantly higher in the initial diagnosed AML. ^***P < 0.001; ^*P < 0.05.

The ROC curve analysis of the diagnostic value between MDGs, core DEmRNAs and co-expressed DElncRNAs

To further verify the diagnostic value between MDGs, core DEmRNAs and co-expressed DElncRNAs in AML, we used receiver operating characteristic (ROC) curve analysis to evaluate these genes. The study results showed that the area under the curve (AUC) of MDGs, core DEmRNAs and co-expressed DElncRNAs were both ≥ 0.7 (Figure 10), indicating that the critical genes that were screened in this study are good diagnostic value for AML.

Figure 10. Weights of the key mRNAs and co-expressed DElncRNAs were determined by the AUC values of the ROC curves in AML. (A) ROC analysis revealed that the AUC for 18 mRNAs was ≥ 0.7; (B) ROC curve analysis shows that the AUC for 7 co-expressed DElncRNAs was also ≥ 0.7. Abbreviations: ROC: receiver operating characteristic; AUC: area under the curve; DElncRNAs: differentially expressed lncRNAs.

Retrospective research about the relationship between the expression of ITGAM and AML patient's characteristics

To understand the correlation between the expression level of ITGAM and the clinical characteristics of AML patients, we also investigated the clinical data of 179 AML patients and analyzing the impact of ITGAM expression on patient's features through a retrospective study. We have reported the general information of AML patients with differentially expressed ITGAM in Table 2. The study results show that the positive expression of ITGAM is substantially related to AML classification (P < 0.001) and a higher count of white blood cells (P < 0.01). It is noteworthy that the complete response (CR) rate of AML patients, received primary chemotherapy with negative expression of ITGAM was significantly higher than the positive expression of ITGAM (P < 0.05). Still, there was no remarkable difference in the no remission (NR) rate or the partial remission (PR) between the two groups.

Pathway analysis based on grouping ITGAM as high and low expression in TCGA-LAML

Based on the substantial relationship between positive ITGAM expression and AML classification and clinical features, we further analyzed the pathway enrichment in the high versus low ITGAM expression group using TCGA-LAML expression profile data. The results show that highly expressed ITGAM is mainly enriched in cytokine-cytokine receptor interactions, cell adhesion molecules (CAMs), lysosomes, and hematopoietic cell lineages, while a large number of highly expressed ITGAM group genes are involved in immune infiltration and inflammation-related signaling pathways, such as Intestinal immune network for IgA production, antigen processing and presentation, Fc gamma R-mediated phagocytosis, Toll-like receptor signaling pathway, NOD-like receptor signaling pathway, Natural killer cell-mediated cytotoxicity, MAPK signaling pathway, Fc epsilon RI signaling pathway, and TGF-beta signaling pathway, etc. (Figure 11).

Figure 11. The pathway enrichment in the high versus low ITGAM expression group is based on the TCGA-LAML expression profiles. The significant enriched KEGG pathways were confirmed as an enrichment when FDR < 0.05.

Association of prognostic genes expression levels with immune score and stroma score in TCGA-LAML cohort

Since the highly expressed ITGAM, PPBP, and LncRNA ITGB2-AS1 are highly expressed in our internal cohort, we investigated the association of these three prognostic genes expression levels with immune score and stroma score. Interestingly, we found that the expression level (log2 transformation) of LncRNA ITGB2-AS1 was strongly correlated with the immune score (Spearman's correlation test, R = 0.81, P < 0.001) and moderately correlated with stroma score (Spearman's correlation test, R = 0.49, P < 0.001) (Figure 12A and 12B). Besides, the expression level (log2 transformation) of ITGAM was also strongly correlated with the immune score (Spearman's correlation test, R = 0.82, P < 0.001) and moderately correlated with the stroma score (Spearman's correlation test, R = 0.64, P < 0.001) (Figure 12C and 12D). However, we did not find significant correlations between the expression level of PPBP and immune scores and stromal scores in the TCGA-LAML cohort (R < 0.30, P < 0.01). This result indicated that the expression levels of ITGB2-AS1 and ITGAM are associated with the modulation of immune and stromal activity in LAML.

Figure 12. Association of prognostic genes expression levels with the tumor microenvironment (TME). (A) Strong positive correlation between ITGAM expression (log2 transformation) and immune score. (B) A moderate correlation between ITGAM expression (log2 transformation) and stroma score. (C) Strong positive correlation between LncRNA ITGB2-AS1 expression (log2 transformation) and immune score. (D) A moderate correlation between LncRNA ITGB2-AS1 expression (log2 transformation) and stroma score. R, Spearman’s correlation coefficient; ^***P < 0.001.

Association of ITGAM and LncRNA ITGB2-AS1 expression levels with immune infiltrations in the TCGA-LAML cohort

Since ITGAM and LncRNA ITGB2-AS1 expression levels are correlated with immune and stromal scores in the TCGA-LAML cohort, we investigated the correlations of ITGAM and LncRNA ITGB2-AS1 with the various immune signatures. We revealed significant positive correlations of the ITGAM expression levels with the enrichment levels (ssGSEA scores) of immune inhibitory cells, including Tregs, Th17, MDSC, macrophages, TAM, and M2 macrophages (Spearman's correlation test, P < 0.001) (Figure 13A) but not correlated with other immune cells including B cell, CD8+ T cells, CD4+ regulatory T cells, NK cells, and CAFs. In addition, LncRNA ITGB2-AS1 expression levels are correlated with the enrichment levels (ssGSEA scores) of CD4 regulatory T cells, MDSC, macrophages, TAM, and M2 macrophages (Figure 13B).

Figure 13. Association of ITGAM and LncRNA ITGB2-AS1 expression level with immune signature in AML. (A) The expression of ITGAM exhibit a significant positive correlation with six immune cells (M2 Macrophages, Macrophages, Treg, MDSC, TAM, and Thr17). The Spearman's correlation test P values are shown; (B) The expression of LncRNA ITGB2-AS1 exhibit a significant positive correlation with five immune cells (M2 Macrophages, Macrophages, MDSC, TAM, and CD4 regulatory T cells). The Spearman's correlation test P values are shown; (C) High infiltration levels (ssGSEA scores) of M2 macrophages associated with shorter survival time in LAML patients. (D) the ratios of pro-/anti-inflammatory cytokines are significantly lower in AML with highly expressing of ITGAM and LncRNA ITGB2-AS1 (expression levels > median) than in those lowly expressing of ITGAM and LncRNA ITGB2-AS1 (expression levels < median). The mean expression (log2 transformed) ratio of the marker genes levels was defined as pro-inflammatory cytokines representing immune-stimulatory signature with marker genes as IFN-γ, IL-1, and IL-2. The anti-inflammatory cytokines represent with immune-inhibitory signature with marker genes as TGFB1, IL-10, IL-4, and IL-11. Abbreviations: Treg: The regulatory T cells, TAM: Tumour-associated macrophages, MDSC: Myeloid-derived suppressor cells, TGFB: transforming growth factor–β1. ^***P < 0.001.

Furthermore, we investigated the survival differences between high and low infiltrated cellular groups of all positively correlated immune signatures. We found that the increased infiltration of M2 macrophage groups had a more insufficient survival time than the lower infiltration group (Figure 13C). In contrast, Tregs, Th17, MDSC, macrophages, TAM, and CD4+ regulatory T cells are not correlated with the survival differences.

Finally, the ratios of pro-/anti-inflammatory cytokines are significantly lower in AML with highly expressing of ITGAM and lncRNA ITGB2-AS1 (expression levels > median) than in those lowly expressing of ITGAM and lncRNA ITGB2-AS1 (expression levels < median) (P < 0.001) (Figure 13D).

Altogether, these results suggest that higher expression of ITGAM and lncRNA ITGB2-AS1 are related to the prognosis and immune inhibition in AML. It also could explain that ITGB2-AS1 and ITGAM play a specific role in immune inhibition in AML.

Identification of candidate drugs targeting hub genes

We screened these two essential genes (ITGAM and PPBP) and co-expressed lncRNA (ITGB2-AS1) for drug-gene interactions by using DGIdb 4.0 database [20]. The prediction results showed that we obtained many candidate drugs that potentially target two essential genes. Still, most of the interaction types between target genes and drugs are uncertain (Table 3). Some of the immunosuppressive drugs as rovelizumab and clarithromycin had been identified as a potential medicine for AML treatment. However, as far as we know, the copper target's inhibitory effects on PPBP have not been tested for the treatment of AML. Our data suggest that these genes may be promising targets for developing anticancer drugs to treat patients with AML.

Acquisition and process the multi-omics data from different databases

The transcriptome profiling of GSE114868 contributed by Cheng et al. [57] was download from the GEO database (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE114868). GSE114868 was a gene expression profile based on the GPL17586 platform (Affymetrix Human Transcriptome Array 2.0) containing 194 AML samples and 20 normal samples from the National Taiwan University Hospital. "Limma" [58] program in R software (version 4.0.3) was performed to background correction, normalization, and probe annotation. The DEG screening criteria were set as FC > 2 and P < 0.05. "Pheatmap" [59] and "ggrepel" [60] programs in R software (version 4.0.3) were used to draw DEGs volcano plot and heatmap.

Also, we downloaded the gene expression profiles and clinical data of 151 AML patients from the UCSC Xena [61] database (http://xena.ucsc.edu/) for further analysis. GEPIA [9] database was used to verify selected genes' effect on TCGA-LAML patients' OS. The TCGA-LAML lncRNA data was downloaded from The Cancer Genome Atlas (TCGA) database (https://portal.gdc.cancer.gov/). Furthermore, we validated the expression level of selected genes in an independent cohort based on our 20 clinical samples. We also investigated the relationship between patients' clinical characters and the expression level of a crucial gene in retrospective research based on 179 cases. Finally, we predicted the candidate drugs that target hub genes in the DGIdb 4.0 database (https://dgidb.org/).

Construction of co-expression gene network

After we obtained the DEGs in the GSE114868 transcriptome data, R software (version 4.0.3) "WGCNA" program [62] was further used to perform a weighted gene co-expression network analysis (WGCNA). We adopted a soft threshold β to ensure the scale-free network's stability during the co-expression gene network construction. In the scale-free co-expression network, genes with high correlation were clustered in the same module. The hierarchical clustering of the weighted coefficient adjacency matrix was used to identify the functional modules. The topological overlap matrix (TOM) has been calculated further.

The PPI network analysis and identification of crucial module driver genes

After performing WGCNA on the transcriptome data, we established a total of 8 AML related dysfunctional modules. The R software (version 4.0.3) "ClusterProfiler" program [63] was used to annotate the GO function of DEGs in the Blue module. At the same time, we used Cytoscape software (version 3.8.0) "ClueGO" program [64] to perform KEGG pathway enrichment analysis in Blue module genes. Besides, we also used the STRING database [65] to analyze the PPI in Blue module genes; the interaction network was constructed by interaction relationship score > 0.9. Finally, we used Cytoscape software (version 3.8.0) "MCODE" program [66] to screen for the MDGs, and the top 10 module genes with the highest scores were defined as MDGs. The PPI network and MDRs were both visualized in Cytoscape software (version 3.8.0).

Screening of core differentially expressed mRNAs in whole transcriptome and identification of co-expressed differentially expressed lncRNAs

To further verify the DEmRNAs and co-expressed lncRNAs in the transcriptome. We firstly screened the DElncRNAs and using R software (version 4.0.3) "pheatmap" [59] programs to draw the heatmap of DElncRNAs. Secondly, we used more stringent screening criteria (|logFC ≥ 2|, P < 0.05) to identify DEmRNAs, and perform PPI analysis in the STRING [65] database (score > 0.9). Then, we used Cytoscape software (version 3.8.0) "cytoHubba" [67] program screened the core DEmRNAs, and the top 10 DEmRNAs with the highest contribution (Degree) were defined as the core DEmRNAs in the whole transcriptome. The Cytoscape (version 3.8.0) displays the PPI network and the regulatory network of core DEmRNAs across the entire transcriptome. Finally, the Person correlation coefficient (PCC) [68] was used to analyze the correlation between core DEmRNAs and DElncRNAs. The DElncRNAs-mRNAs pair with a strong correlation (PCC > 0.8, P < 0.05) were selected for the further two-way clustering, the R software (version 4.0.3) "corrplot" [69] program performed to draws a heatmap of the expression correlation between the core DEmRNAs and DElncRNAs.

Prognostic analysis of critical genes in TCGA-LAML

To further verify the prognostic value of both the MDGs, core DEmRNAs, and co-expressed DElncRNAs in AML, we used the GEPIA database to perform clinical verification of TCGA-LAML patients' overall survival (OS) for these critical genes. Moreover, by using the downloaded clinical data of TCGA-LAML from the UCSC Xena [61] database and the "survival" [70] and "nomogramEx" [71] programs in R software (version 4.0.3), we have analyzed the univariate Cox risk proportional model and the Nomogram respectively. Survival analysis uses the Log-rank test for hypothesis testing. A Cox proportional hazard model was used to estimate the essential genes' hazard ratio and 95% confidence interval. P < 0.05 is considered statistically significant.

The ROC curve analysis of the diagnostic value in core DEmRNAs and co-expressed DElncRNAs

Only to verify the diagnostic value of core DEmRNAs and co-expressed DElncRNAs in AML. ROC curve analysis was used for these genes to determine the AUC. AUC ≤ 0.50 indicates no, 0.8 ≤ AUC ≤ 1.0 shows good diagnostic value. Statistical analysis and visualization are based on R software (version 4.0.3) "ROCR" [72] program.

qRT-PCR validated the expression of critical genes in initial diagnosis patients with AML

To further increase the research results' reliability, we verified the expression level of three essential genes in bone marrow mononuclear cells of newly diagnosed AML patients (WHO-AML criteria confirmed the diagnosis [73], cases with AML-M3 were excluded) by qPCR. The bone marrow samples were collected in heparinized tubes before treatment and shipped to the laboratory within 24 to 36 hours. The leukemic cells were isolated by density gradient centrifugation using 1.077 g/mL Ficoll-Isopaque (Pharmacia). The proportion of leukemic cells was estimated using May-Grünwald- Giemsa-stained cytocentrifuge preparations and light microscopy. The cell samples selected for analysis contained at least 90% blasts after separation. Pellets of 2 to 10 million cells were stored in TRIzol (Invitrogen, Carlsbad, CA, USA) and immediately frozen at –80°C. Peripheral blood mononuclear cells from 15 anonymized healthy volunteers were included as control samples.

cDNA was generated using the Reverse Transcription Kit (Foregene, Chengdu, China). The expression levels of ITGAM, PPBP, and ITGB2-AS1 were quantified using SYBR Green Master Mix (SYBR GREEN, Beijing, China), and the housekeeping gene GAPDH was used as an internal control. The following primers used in Table 4. qRT-PCR performed on ViiATM 7 System software (Thermo Fisher Scientific, ABI7500, USA). The results were normalized to the expression of GAPDH and presented as the fold change (2^−ΔΔCT).

Evaluation of the impact of critical genes on patient's clinical characteristics through a retrospective study

Also, to further study the effect of ITGAM on the patient's clinical characteristics with AML, we have retrospectively investigated 179 cases of newly diagnosed with AML admitted to the First Affiliated Hospital of Xinjiang Medical University from March 2019 to September 2020. The patients included in this study all received Immunophenotyping tests according to WHO diagnostic criteria. The test process followed the manufacturer's indication as these steps: we obtained 1–2 ml bone marrow from patients with initial diagnosed AML and collecting in an anticoagulated tube (1–10 × 10⁵ cells/ml). Mixed the bone marrow (50μL) with fluorescein (FITC, PE, PerCP, and APC; BD Bioscience, USA) labeled antibodies (5–20 μL). Incubate for 15 minutes without light at room temperature. Add 1ml hemolysin at room temperature for 10 minutes without light. Centrifuge at 1500 rpm for 5 minutes, then discarded the supernatant. After that, the precipitation was washed twice with PBS, and the cells were analyzed using a flow cytometer (FACSCalibur, BD Bioscience, USA).

Finally, we analyzed the relationship between the expression of these essential genes and the clinical characteristics and chemotherapy efficacy of AML patients. The outcomes of chemotherapy presented as CR, PR, NR, and death. The statistical analyses used GraphPad Prism 8.0 software (CA, USA) as appropriate. Comparing patient characteristics and chemotherapy outcomes between ITGAM positive and negative groups adopts Student t and Chi-square test. The features with P < 0.05 were considered a significant difference between the two groups.

Pathway analysis based on grouping ITGAM with high and low expression in TCGA-LAML

To explore the pathway enrichment after grouping ITGAM with high and low expression, we downloaded the expression profile files of AML samples (n = 151) from the TCGA database and performed DEGs analysis after grouping them with the median ITGAM expression. Then the GSEA [74] analysis was used to identify the enriched KEGG pathways, and the significant pathways were confirmed as an enrichment when FDR < 0.05.

Evaluation of the ESTIMATE scores and immune signature enrichment levels

ESTIMATE is an algorithmic tool based on the R package for predicting tumor purity, which uses the gene expression profiles of 141 immune genes and 141 stromal genes to generate ESTIMATE scores, immune scores, and stromal scores [18]. The presence of infiltrated immune cells and stromal cells in tumor tissues were calculated using related gene expression matrix data, represented by immune and stromal scores.

We also identified the immune signature's enrichment level in the TCGA-LAML sample as the single-sample gene-set enrichment analysis (ssGSEA) score [75]. The gene set contains the collection of all marker genes of an immune signature. We included 11 immune signatures: B cell, CD8+ T cells, CD4+ regulatory T cells, NK cells, Tregs, MDSC, TAM, macrophages, M2 macrophages, CAFs, and Th17. The threshold is the absolute value of R is not less than 0.30 and P < 0.01. The marker genes set of individual immune signatures displayed in Supplementary Table 1.

Furthermore, we used the Wilcoxon test to compare the ratio of pro-inflammatory cytokines and anti-inflammatory cytokines between the TCGA-LAML with high expression levels of the ITGAM and lncRNA ITGB2-AS1 (expression levels > median) and the TCGA-LAML with low expression levels of the ITGAM and lncRNA ITGB2-AS1 (expression levels < median). Besides, the ratio of mean expression (log2 transformed) of the marker genes levels defined as the ratio of pro-inflammatory cytokines (IFN-γ, IL-1, IL-2) and anti-inflammatory cytokines (TGFB1, IL-10, IL-4, IL-11) in AML [76].

In addition, to investigate the survival differences between high and low infiltrated cellular groups of all positively correlated immune signatures, the immune cells infiltrate data, and the TCGA-LAML survival data (n = 151) were used in the K-M survival analysis, "survival" programs [70] in R software (version 4.0.3) were used to analyze and visualize the results.

Prediction of candidate drugs target on essential genes

We identified the drugs that target the hub genes using DGIdb 4.0 [20]. DGIdb collects drug-gene interaction data from 30 disparate sources, including ChEMBL, DrugBank, Ensembl, NCBI Entrez, PharmGKB, PubChem, Clinical Trial Databases, and literature in NCBI PubMed. The drug-gene interactions supported by at least one database or PubMed reference were identified. From the identified drug-gene interactions, we selected the drugs that have been shown in the DGIdb database.

Ethics approval and consent to participate

All research was conducted with integrity and in line with generally accepted ethical principles and approved by the Research Ethics Committee of the First Affiliated Hospital of Xinjiang Medical University (reference number: 220200320). The study was performed following the principles of the Declaration of Helsinki. All of the participants provided written informed consent to participate in this study. The ethics mentioned above committee approved this consent procedure. All personal information relating to the samples involved in the study was anonymized.