We downloaded RNA-seq BAM files, mRNA expression matrixes, and clinical features of acute myeloid leukemia (TCGA-LAML) from the TCGA data portal (https://portal.gdc.cancer.gov/) (Weinstein et al., 2013). The validation datasets were downloaded from the Beat AML (http://www.vizome.org/) (Tyner et al., 2018a) and Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/, GSE37642 (Li et al., 2013), GSE12417 (Metzeler et al., 2008), and GSE10358 (Tomasson et al., 2008)). To ensure the stability of the prediction, patients with survival time <1 month were excluded. The PRO-seq BAM files, Histone ChIP-seq, and P300 ChIP-seq data of the Kasumi-1 cell line were downloaded from the GEO database (GSE83660 and GSE100446) (Zhao et al., 2016; Xu et al., 2019). From the Encyclopedia of DNA Elements (ENCODE) portal (https://www.encodeproject.org/), we downloaded the Hi-C interactions data of K562 and GM12878 cell lines with the following identifiers: ENCFF996XEO and ENCFF355OWW (Inoue et al., 2017). The Histone ChIP-seq data of the mononuclear cells of the bone marrow were collected from the Blueprint project (IHEC portal: https://epigenomesportal.ca/ihec/index.html, IHECRE00000277) (Adams et al., 2012; Stunnenberg et al., 2016). The ribosome-deleted RNA-seq data of two AML patients’ BM samples were obtained from GSE87285 (Loke et al., 2017). All detailed information on datasets mentioned above is listed in supplementary table 1.

Annotation of enhancer loci was downloaded from the FANTOM5 project (http://fantom.gsc.riken.jp/5/) (Andersson et al., 2014). We filtered out eRNA loci overlapping with known annotated elements and 1 kb extension from both the transcription start site (TSS) and transcription end site (TES) to eliminate introns, exons, promoters, rRNA repeats, and other transcribed elements. The 3 kb of the middle locus of the enhancer was designated as the potential eRNA region (Zhang et al., 2019). RNA-seq data of the TCGA-LAML cohort were mapped to eRNA regions by featureCounts and normalized to reads per million (RPM) (Liao et al., 2014). eRNA with RPM ≥1 in at least one sample was defined as an active eRNA (Zhang et al., 2019; Zhang et al., 2021). The enhancer loci from all the human tissues and cell lines in the database were included in the analysis. Liftover was used to convert genomic regions from various genomic versions to the hg19 genome version (Rhead et al., 2010). The annotation of enhancer loci was supported by R package annotatr (Cavalcante and Sartor, 2017). The histograms of promoters, enhancers, and random loci were created by HOMER (http://homer.ucsd.edu/homer/) (Heinz et al., 2010). Immune-related eRNAs (ir-eRNAs) were defined as at least correlated with an immune cell type (absolute Spearman’s correlation Rs > 0.3, p < 0.05). We used univariate Cox regressions (p < 0.01, hazard ratio (HR) > 1.1 or <0.9) to identify prognostic eRNA.

The pair of eRNAs and protein-coding genes was identified based on close distance (1 MB) and the corresponding expression level (absolute Spearman’s correlation Rs > 0.3 and FDR <0.05) concurrently (Zhang et al., 2019). We selected the defined eRNA loci, promoter of protein-coding genes, and random loci of the genome (0.5–4.5 kb away from the region of the defined enhancer and excluded any overlapping annotation) to exhibit the genome coverage. The Histone and P300 ChIP signal was visualized by Deeptools, and the HiC data were illustrated by HiCExplorer (Wolff et al., 2020).

To investigate distinct patterns, we used the R package ConsensusClusterPlus for unsupervised clustering analysis (Wilkerson and Hayes, 2010). The number of iterations was set at 1,000 to ensure clustering stability. The other parameters, clusterAlg, distance, and linkage, were set as pam, Pearson, and ward.D.

GSVA is a nonparametric, unsupervised method for assessing alternations in the activity of signaling pathways and biological processes using expression datasets. We compared signaling pathways across different ir-eRNA patterns using the R package GSVA (Hänzelmann et al., 2013). The gene set “c2. cp.kegg.v6.2. symbols” from the Molecular Signatures Database (MSigDB) was used as the predetermined dataset. Statistical significance was defined as an adjusted p value of less than 0.05. We used the R package clusterProfiler to annotate the ir-eRNA-related gene for Gene Ontology (GO) analysis by biological processes (BPs) (p < 0.05, FDR <0.05) and the differentially expressed genes (DEGs) for gene set enrichment analysis (GSEA) by the Kyoto Encyclopedia of Genes and Genomes (KEGG) with a default parameter.

By generating a gene set using cell-specific markers, single-sample gene set enrichment analysis (ssGSEA) allows us to assess the number of immune cells in AML BM samples. A previous study yielded gene sets for labeling each kind of TME-infiltrated immune cells, such as activated dendritic cells, activated CD8⁺ T cells, macrophages, active CD8 T cells, and regulatory T cells (Barbie et al., 2009; Charoentong et al., 2017). The relative abundance of TME cells in each sample was represented by the estimated abundance score in silico. In addition, the ESTIMATE method was also used to compute immune and stromal scores, as well as tumor purity, corroborating the ssGSEA approach.

Based on prognostic ir-eRNAs, AML patients were split into two different expression patterns. To determine pattern-related genes, we used the R package limma to compare the expression profiles between the two patterns. An empirical Bayesian algorithm using moderated t-tests was used to carry out this strategy. The adjusted p value of 0.05, determined with Benjamini–Hochberg correction, was set as the significant criteria for detecting DEGs (Ritchie et al., 2015).

We developed a scoring system to assess the immune pattern of individual AML patients. First, we standardized the expression of DEGs recognized from different immune subtypes in all TCGA samples. Patients were divided into two subgroups by DEGs through unsupervised clustering analysis. Next, we used a univariate Cox regression model for prognostic analysis of each gene in the signature (p < 0.01, HR > 1.1 or <0.9). We calculated the correlation between different subtypes and signature genes (Pearson correlation analysis). The positively correlated genes were defined as gene set A, and the negatively correlated genes were defined as gene set B. We then performed principal component analysis (PCA) to construct pattern-related signatures. The benefit of this strategy was that it concentrated the scores on the dominant network of gene signatures while downweighting contributions from genes that did not correlate with other set members. Principal component 1 of each gene set was extracted as PC1. A and PC1.B to calculate E-score, which was defined in an approach similar to that of the Genomic Grade Index (GGI) (Van 'T Veer et al., 2002; Zeng et al., 2019): E _ s c o r e = ∑ ( P C 1. A i + P C 1. B i )   , where i is the expression of a gene related to each signature gene set.

Spearman’s rank correlation analysis evaluated the correlation coefficient between the expression of ir-eRNAs and TME-infiltrated immune cells (Sokol et al., 2017). We utilized the R package survminer to identify the cutoff point between each E-score subgroup based on the connection between the E-score and the overall survival of the patients. To get the maximum rank statistic, all possible cutoff points were evaluated repeatedly. Patients were separated into high- and low-E-score groups using log-rank statistics to minimize batch effects for maximum selection (Kerschke et al., 2020). The Kaplan–Meier curves were used to create survival curves for prognostic analysis, and the significance of differences was determined using the log-rank test. The hazard ratios (HRs) and p values of ir-eRNAs and signature genes were calculated using a univariate Cox regression model. To find independent prognostic factors, we employed a multivariate Cox regression model. The R package forestplot was applied to visualize the results of univariate and multivariate Cox regression. Receiver operating characteristic (ROC) curves were used to assess the sensitivity and specificity of E-score, and the area under the curve (AUC) was determined using the R package survivalROC. The R package RCircos was used to illustrate the chromosomal location of the ir-eRNAs and target genes. The significance criterion for all statistical comparisons was 0.05, and all comparisons were two sided. R 4.0.2 was used to conduct all statistical analyses.

All expression profiles and clinical annotation files were downloaded from public databases. Our code is available on request.

The overview of the study is depicted in Figure 1. We started by gathering approximately 65,000 enhancer loci from the FANTOM5 project, a robust enhancer database including enhancers from most human tissues and cell lines (Andersson et al., 2014). Then, we identified 13,251 enhancer RNAs in the TCGA-LAML cohort, none of which overlapped with annotated genomic regions. Active enhancer loci were typically described as eRNA transcription and unique epigenomic modification patterns, such as high H3K4me1, high H3K27ac, low H3K4me3, and high P300 levels (Klemm et al., 2019). To verify the reliability of enhancer loci in AML, we then investigated the epigenomic modification patterns of 13,521 eRNA regions, which showed enhancer-like epigenomic patterns compared with random loci and protein-coding gene loci (Supplementary Figure S1). In addition, the distribution of RNA transcription abundance revealed obvious transcription peaks in a classic manner among enhancer regions, demonstrating that the expression of enhancer regions was not the transcription background noise (enhancer loci vs. random loci, Figure 2A). On average, the expression level of eRNAs was less than that of protein-coding genes (enhancer loci vs. random loci and transcription start site (TSS) of protein-coding gene vs. random loci, Figures 2A,B). Enhancers usually play roles as cis-elements (Klemm et al., 2019). Previous studies revealed that eRNAs were excellent markers for active enhancers, so eRNAs usually correlated with adjacent genes. To illustrate the characteristics of eRNAs, pairwise Pearson correlation analysis revealed that eRNAs were more correlated with neighboring genes in a distance-related way, which indicated that eRNA has a stronger association with the closest gene than distal ones (Figure 2C). The closest coding genes near active eRNAs (RPM >1) exhibited a higher expression level than those near nonactive eRNAs (p < 0.001, Wilcoxon test; Figure 2D). Similar to previous studies, we observed that most eRNAs were mapped to interact with one adjacent gene (51.20%, 1,171/3,459, Figure 2E), and most genes were only matched to a single eRNA as well (59.74%, 2,664/4,459, Figure 2F). In conclusion, we demonstrated the reliability of the defined enhancer loci and presented the landscape of eRNA in AML. Our results and recent studies suggested that eRNA levels were intimately associated with the expression of nearby genes (Andersson et al., 2014; Chen and Liang, 2020).

We then concentrated on eRNAs correlated with prognosis, which may play crucial roles in tumor progression and proliferation. We filtered out 258/13,251 prognostic eRNAs (Supplementary Table S2) by univariate Cox regression. To explore the potential functions of prognostic eRNAs, co-expression analyses were performed. Correlated mRNAs in the 1 Mb region were identified as potential target genes (Rs > 0.3, Spearman correlation analysis). Functional enrichment analysis revealed that the bitter taste-related pathway and several vital oncogenic signaling pathways such as PI3K-Akt-mTOR and MAPK were enriched in GO enrichment analysis (Figure 3A) (Nepstad et al., 2020; Kim et al., 1999). Surprisingly, prognostic eRNAs were most concentrated on immune-related pathways, such as myeloid leukocyte activation and granulocyte migration (Figure 3A). Therefore, we suggested that AML blasts shape immune microenvironment patterns through eRNAs and related genes. We evaluated the immune cell abundance of each sample by ssGSEA and filtered 72 eRNAs correlated with immune cell abundance (Supplementary Figure S2), which were defined as ir-eRNAs. Due to the diverse expression of ir-eRNAs in different samples, we explored whether ir-eRNAs exhibited distinct expression patterns in AML. Unsupervised clustering analysis was performed, and two clusters were identified (Supplementary Figures S3A–E, Figure 3B). One cluster defined as an immune-resistant subtype showed a higher infiltration level of suppressive cells, including dendritic cells, macrophages, regulatory T cells, and myeloid-derived suppressive cells (MDSCs). The other cluster defined as the immune-active subtype showed a lower level of suppressive cells and a trend of a higher level of activated T cells without statistically significant at the 5% level (Figure 3C). The stromal and immune scores determined by ESTIMATE were also higher in the immune-resistant subtype (stromal score, p < 0.001, Wilcoxon test; immune score, p < 0.001, Wilcoxon test; Figure 3D), while the immune-active subtype had a higher tumor purity (p < 0.001, Wilcoxon test; Figure 3E), in line with the work of Yan et al. (2019). Cell clustering revealed complicated correlation networks of 23 immune cells and their prognostic significance by the log-rank test (Figure 3F). In addition, Figure 3G shows the detailed hazard ratio of each immune cell by univariate Cox regression analysis. Plasmacytoid dendritic cell and macrophage, previously reported as cancer-promoting cells, had hazard ratio >1 (HR, 2.06 (1.314–3.228), p = 0.001; HR, 1.694 (1.083–2.648), p = 0.023, respectively) (Miari et al., 2021; Xiao et al., 2021). The Kaplan–Meier curves validated the poor prognosis of several immune cells, and the representative cell types are shown in Supplementary Figures S3F–P (log-rank test, cutoff determined by R package survminer). Furthermore, patients of the immune-resistant subtype had a significantly shorter overall survival time than patients of the immune-active subtype (log-rank test, p < 0.001, Figure 3H). In general, we discovered that two clinically significant immune subtypes characterized by distinct ir-eRNAs profiles existed in AML, suggesting that ir-eRNAs played vital roles in modulating the immune microenvironment of BM. The immune-resistant subtype was characterized by a high level of suppressive and cancer-promoting immune cells. The immune-active subtype was accompanied by a relatively high level of protective immune cells. Although the absolute abundances of activated T cells, including CD4 T cells and CD8 T cells, were not significant between the two types, fewer suppressive immune cells in the immune-active subtype provided an immune microenvironment with less immune evasion and T-cell dysfunction.

Based on the abovementioned findings, we applied the GSVA algorithm to explore pathways that differed between two immune subtypes. The result showed that the immune-resistant subtype was highly associated with innate immune response pathways, while the immune-active subtype was related to the T-cell receptor signaling pathways (Supplementary Figure S4A). The GSEA of immune-related pathways further confirmed that the NF-κB pathway and other innate immune response pathways were highly enriched in immune-resistant subtypes (Figure 4A). Surprisingly, the PDL1 expression and PD1 checkpoint pathway showed a close connection to the immune-resistant subtype, leading to the following analysis to predict the response to ICI therapy. To further explore the possible biological roles of each immune subtype, we identified 1,057 ir-eRNA-related differentially expressed genes (DEGs) between immune-resistant and immune-active subtypes by the limma package (Supplementary Table S3). To investigate whether the DEGs directly correlated with ir-eRNAs, we subsequently screened genes with correlation index >0.3 in the approximately 1 Mb adjacent regions of ir-eRNAs (Zhang et al., 2019; Zhang et al., 2021) and identified 148 potential target genes (Supplementary Table S4). The Venn diagram illustrated that 55 genes showed both different expression levels between two clusters and correlations with ir-eRNAs, which suggested that these genes might be directly activated by immune-related enhancers (Supplementary Figure S4B). Most of these genes were reported as immune-related genes in AML (Li et al., 2020; Brenner and Bruserud, 2018; Borrego, 2013). Among these genes, we focused on five genes (S100A4, S100A6, S100A8, S100A9, and S100A12) from the S100 protein family that was widely reported in AML. The five genes were correlated to a specific eRNA transcribed from a chromosomal region (chr1:153488459–153488748) (supplementary table 4, Figure 4B). We defined this enhancer and eRNA as S100-enhancer and S100-eRNA, respectively. In addition, Supplementary Figure S4C indicates the compact correlation among the five genes. Figure 4C reveals that the five genes exhibited consistent higher expression in the immune-resistant subtype (S100A4 p < 0.001, S100A6 p = 0.001, S100A8 p = 0.014, S100A9 p = 0.006, S100A12 p = 0.056, Wilcoxon test). Brenner and Bruserud (2018) summarized that the S100 protein family promoted tumor growth by enhancing a suppressive immune microenvironment in AML. Thus, we examined the relationship between S100-eRNA and immune cells and found a strong correlation with suppressive immune cells, including MDSC, macrophage, and regulatory T cell, as was expected (Supplementary Figure S4D). Figure 4D demonstrates that the S100 protein family tightly correlated with the immune checkpoint-related gene expression, suggesting the potential efficacy in predicting response to ICI treatment. Because S100 proteins were positively correlated with the immune-resistant subtype, a relatively poor prognosis would be observed. Supplementary Figures S4E–I demonstrate the hypothesis. We then expanded the S100 protein family to pan-cancer survival analysis. S100 proteins represented general risky factors in the pan-cancer scale, except for S100A6/A8/A9 in Sarcoma (SARC) and S100A4/A8 in Head and Neck Squamous Cell Carcinoma (HNSC) (univariate regression analysis, Figure 4E). These results suggested that the S100-enhancer and S100-eRNA generally facilitated a suppressive inflammatory microenvironment, resulting in a worse outcome via the S100 protein family and related pathway, as Figure 4F illustrates. The illustration of the genomic region revealed that the S100-enhancer and correlated genes were in the same transcription activation domain (TAD) in leukemia cells, which suggested S100-enhancer interacted with nearby genes and induced their transcription (Figure 5). Additionally, PRO-seq demonstrated an active transcription event in the enhancer region, and histone modifications confirmed the epigenomic characteristics of S100-enhancer.

To explore the roles of immune microenvironment patterns in predicting clinical prognosis and potential response to immunotherapy, we created a gene signature based on prognostic DEGs between immune-resistant and immune-active subtypes (246/1,057 prognostic DEGs, Supplementary Table S5). Random Forest regression, a machine learning-based algorithm, was used to get representative genes (66/246 after Random Forest regression, Supplementary Table S5). We then performed an unsupervised clustering analysis with the 66 signature genes and obtained two gene clusters (Supplementary Figures 5A–D). The immune infiltration and survival analysis based on both ir-eRNAs and signature genes were similar, thus supporting the robustness of immune subtypes (Supplementary Figures 5E–G). The result of PCA revealed the well-defined characteristics of immune subtypes as well (66 signature genes, Figure 6A).

The abovementioned results demonstrated that ir-eRNA alteration had a significant impact in generating distinct TME patterns. These findings, however, were limited to the population level and could not properly predict the immune profiles of AML at the individual level. Given the alternative heterogeneity and diversity of ir-eRNA-mediated immune profiles, we developed a scoring system based on these phenotype-related genes to quantify the TME pattern of individual patients in AML through the PCA algorithm, which was named eRNA-derived score (E-score) (see Methods, Figure 6B). E-score well distinguished the two immune subtypes as was expected (p < 0.001, Wilcoxon test; Figures 6B,C). The high-E-score group had a better prognosis than the low-E-score group (log-rank test, p < 0.001, Figure 6D). Nearly all immune-active subtypes were stratified into the high-E-score group (Figure 6C).

AML was commonly classified by the French–American–British (FAB) classification and the morphology, immunology, cytogenetics, and molecular biology (MICM) classification via genetic, immunophenotypic, pathological, and clinical features. CALGB cytogenetics-based risk classification and ELN 2017 risk classification provided widely used predictive models based on the integration of cytogenetics and mutational status (Döhner et al., 2017; Tallman et al., 2019). We examined the efficacy of E-score in distinguishing various clinical characteristics by the TCGA cohort and another independent cohort, Beat AML. It was found that the high- and low-E-score groups showed distinct clinical characteristics (Figure 6E). When favorable CALGB risk patients were almost concentrated in the high-E-score group, the low-E-score group contained a greater proportion of intermediate and poor CALGB risk patients (Figure 6F). The M3 subtype was categorized exclusively into the high-E-score group, while M4 and M5 subtypes were most concentrated in the low-E-score group (Figure 6G). Furthermore, PML-RARα translocations only existed in the high-E-score group (Figure 6H). On the contrary, CBFB-MYH11, GATA2-MECOM, and MLLT3-KMT2A only occurred in the low-E-score group (Figure 6H). Differences in other clinical features, including the response to chemotherapy, ELN2017 risk classification, and WHO classification, are shown in Supplementary Figures 5H–J. The group with both low E-score and intermediate/poor risk showed poor prognosis compared with the other groups, which indicated that a good proportion of patients with intermediate/poor risk showed relatively better outcomes (Figure 6I). We have mentioned that the immune-resistant subtype was enriched in the PD−L1 expression and PD−1 checkpoint pathway (Figure 4A). The immune-resistant TME with abundant suppressive immune microenvironment generally adopted a better response to immune checkpoint therapy, such as the anti-PD1/PDL1 agents, in solid tumors (Turan et al., 2018). We predicted the response to immunotherapy in the two groups. The Microsatellite Instability (MSI) score and T-cell dysfunction score of Tumor Immune Dysfunction and Exclusion (TIDE) were higher in the low-E-score group, suggesting a predictive better response to immunotherapy (Figure 6J). We also detected other non-redundant immune checkpoints and checkpoint-related gene expression between the two clusters (Figure 6K). Several checkpoint molecules, including IFNG, TNFRSF9, CTLA4, CD86, TIM-3, and PD-L2, had significantly increased expression in the low-E-score group, and their effects on T-cell dysfunction and immune escape have been verified in AML (Hobo et al., 2018; Kikushige, 2021). However, there was no significant difference, but only a trend, in the expression of PD1 and PDL1, implying that the TME modulation is complicated and understudied. On the other side, our result partly explained the reason why current clinical trials of anti-PD1 agents exhibited controversial results. A group of intermediate/poor-risk patients with low E-score might be the potential population benefiting from the ICI treatment, which requires clinical trials for further validation. Based on the aforementioned results, our E-score system firstly proposed a quantizable system to predict response to ICI therapy in AML.

To discover whether E-score could provide satisfactory predictive efficacy in AML, we proved that E-score was an independent prognostic factor in 4 validation cohorts (Beat AML, GSE37642, GSE12417, and GSE10358). Consistent survival differences were observed in all the validation cohorts mentioned above (Beat AML cohort, log-rank test, p = 0.021; GSE37642, log-rank test, p = 0.020; GSE12417, log-rank test, p = 0.003; and GSE10358, log-rank test, p = 0.019; Figures 7A–D). Both univariate cox regression and multivariate cox regression confirmed the independent predictive power of age, CAGLB risk classification, and E-score (TCGA-LAML cohort, Supplementary Figure 6A). Log-rank regression of these clinical features also revealed that age and CAGLB risk classification were independent predictors (TCGA-LAML cohort, Supplementary Figures 6B–F). E-score showed a high area under the curve (AUC) in a time sequence (AUC_1year = 0.741, AUC_3year = 0.730, AUC_5year = 0.810, Figure 7E). Additionally, compared with other clinical characteristics (gender, age, race, FAB subtype, and CAGLB risk classification), the E-score exhibited the best capability at predicting OS in the TCGA cohort (Figure 7F). Given that the E-score showed an excellent predict power in OS, we constructed a nomogram integrating E-score with other clinical features (TCGA-LAML cohort as the training cohort, Figure 7G). The calibration curves of the nomogram for 1-year, 3-year, and 5-year survival probability exhibited excellent consistency with the ideal performance, indicating the high accuracy of our nomogram (TCGA-LAML cohort, Figure 7H). The nomogram discrimination as evaluated by C index was 0.714 (95% CI = 0.651–0.777, p < 0.0001), which was superior to that of the E-score alone (0.674, 95% CI = 0.608–0.740, p < 0.0001), CAGLB risk classification (0.594, 95% CI = 0.526–0.663, p < 0.0001), and age (0.661, 95% CI = 0.596–0.726, p < 0.0001) in the TCGA-LAML cohort. The nomogram was externally confirmed in the validation cohort, the Beat AML cohort. The C index was 0.684 (95% CI = 0.614–0.755, p < 0.0001; Beat AML cohort), indicating a significant discriminative ability. The calibration curve presented a good agreement for 1-year OS rates in the validation cohort (Supplementary Figures 6G). However, calibration curves exhibited underestimated 3- and 5-year OS rates (Supplementary Figures 6H–I). In all, the E-score system had outstanding predict power in both OS and response to ICI therapy, which improved the current predictive models and made it potential for clinical practice.