Patients and samples
This study was reviewed and approved by the regional ethics review board in Stockholm (2008/1330-31/2 and 2017/2085-31/2) and the institutional ethics committees of Kyoto University (G608 and G1110) and participating institutions (Hyogo Prefectural Amagasaki General Medical Center, Chugoku Central Hospital, Dokkyo Medical University Saitama Medical Center, Gifu University Hospital, Gifu Municipal Hospital, Hyogo Medical University, Hokkaido University, Japanese Red Cross Kyoto Daini Hospital, Kobe City Medical Center General Hospital, Kurashiki Central Hospital, Kitano Hospital, Kyoto Medical Center, Kyoto City Hospital, Matsushita Memorial Hospital, Japanese Red Cross Nagano Hospital, National Cancer Center Hospital, NTT Medical Center Tokyo, Osaka International Cancer Institute, Japanese Red Cross Osaka Hospital, University of Osaka, Otsu Red Cross Hospital, Shizuoka City Shizuoka Hospital, Shinko Hospital, Shiga General Hospital, Sumitomo Hospital, Takeda General Hospital, Takatsuki Red Cross Hospital, Uji Tokushukai Medical Center and Japanese Red Cross Society Wakayama Medical Center) (no. G608), and was performed in accordance with the Declaration of Helsinki. Informed consent was obtained from all participants at participating institutions.
A total of 1,563 patients diagnosed with AML and related neoplasms were consecutively enrolled from the participating institutes and hospitals between 1997 and 2022 in Sweden and between 2011 and 2022 in Japan. Although no preselection criteria were applied with respect to clinical characteristics, sample inclusion was contingent on the availability of viable cryopreserved tumour cells suitable for ATAC-seq for consecutively diagnosed patients with AML. Diagnoses were made at each participating institution and hospital in accordance with the WHO classification in use at the time. In Sweden, patients were registered through the national AML registry, which prospectively collects clinical and genomic data on all newly diagnosed cases. In Japan, patients were diagnosed and treated at participating hospitals, with biospecimens and clinical information sent to and managed by the Kyoto University biobank, where clinical annotations were updated annually. Treatment was administered according to institutional standards of care and, for Swedish individuals, according to the national treatment guidelines for AML, with a subset of patients participating in clinical trials. All samples analysed in this study were collected at the time of initial diagnosis before treatment. Detailed clinical annotations, including diagnosis, demographic and laboratory data, treatment regimens and outcomes, were extracted from electronic medical records and the Swedish AML registry. No statistical methods were used to predetermine sample size. Patient characteristics and diagnoses are summarized in Supplementary Tables 1 and 2.
Tumour samples, such as bone marrow or peripheral blood, and matched control buccal samples, were obtained from patients. We also analysed normal bone-marrow samples from 25 individuals without haematological malignancies who underwent hip joint replacement surgery, serving as controls for the ATAC-seq analysis. Bone marrow and peripheral blood cells were isolated, subjected to erythrolysis or mononuclear cell isolation by Ficoll gradient centrifugation, resuspended in CELLBANKER 1 solution (Nippon Zenyaku Kogyo) or in 10% dimethyl sulfoxide (DMSO) (Merck KGaA) with 90% fetal calf serum (FCS; Thermo Fisher Scientific) and cryopreserved in liquid nitrogen. Genomic DNA was extracted using the QIAamp DNA Mini Kit (QIAGEN), the Gentra PureGene Kit (QIAGEN) or the Maxwell RSC Genomic DNA Kit (Promega). RNA was extracted using the RNeasy Mini Kit (QIAGEN) or the Maxwell RSC simplyRNA Tissue Kit (Promega). A summary of this cohort, including diagnosis, ATAC subgroup, AML classifications, driver genes and available multiomics data, is provided in Supplementary Table 3.
Cell lines
K562, KG-1, THP-1, SKM-1, KY821 and NOMO-1 cells were obtained from the RIKEN BioResource Center Cell Bank (Tsukuba, Japan); Kasumi-1, HL-60, MOLM-13 and TF-1 cells from the American Type Culture Collection (ATCC); and OCI-AML3 cells from the German Collection of Microorganisms and Cell Cultures (DSMZ). These cell lines were authenticated by short tandem repeat profiling and tested for mycoplasma by the providing cell banks.
Targeted-capture sequencing
Targeted-capture sequencing was performed using the SureSelect custom kit (Agilent Technologies), with an in-house gene panel including 331 known AML driver genes (Supplementary Table 4) and an additional 1,158–1,317 probes for copy-number detection. Captured targets were sequenced using the NovaSeq 6000 (Illumina) or DNBSEQ-G400 (MGI) with a 150-bp paired-end read protocol.
Sequence alignment and mutation calling were performed using the hg19 reference genome and Genomon pipeline25,51,52,53 (https://github.com/Genomon-Project). Unless otherwise stated, sequencing data were aligned to the hg19 reference genome. The called variants were further filtered by assessing the oncogenicity of variants on the basis of an in-house curation program25,51,52,53 that uses the COSMIC database (v.96), an in-house blacklist of error calls and public SNP databases, including the 1000 Genomes Project (October 2014 release), NCBI dbSNP build 138, National Heart, Lung, and Blood Institute (NHLBI) Exome Sequencing Project (ESP) 6500, the Human Genetic Variation Database (HGVD) and our in-house dataset.
For copy number analysis, SNP probes included in the target bait for targeted-capture sequencing were used to allow the detection of copy-number changes and allelic imbalances25,51,52,53. This program (CNACS) is available at https://github.com/papaemmelab/toil_cnacs. A total copy number (TCN) of 2.22 or higher was defined as gain, and a TCN lower than 1.88 was defined as loss. Copy-number-neutral loss of heterozygosity was called with a B-allele frequency lower than 0.90 and a TCN between 1.88 and 2.22. Arm-level changes were called for regions with a total length greater than 1 Mb. Detected copy-number changes were manually curated.
Structural variations (SVs) were detected using the Genomon SV pipeline54, which uses both breakpoint-containing junction read pairs and improperly aligned read pairs. Detected putative SVs were filtered by removing (i) those with fewer than four supporting tumour reads and fewer than ten supporting tumour/normal reads and (ii) those present in control normal samples, whose breakpoints were manually inspected using the Integrative Genomics Viewer (IGV). For KMT2A and MECOM rearrangements, cases other than t(9;11)(p21.3;q23.3) (KMT2A::MLLT3) and inv(3)(q21.3q26.2) or t(3;3)(q21.3;q26.2) (GATA2::MECOM) were classified as SVs with rare partners and annotated with a distinct colour, reflecting their separate classification from KMT2A::MLLT3 and GATA2::MECOM rearrangements in the ICC classification.
WGS
WGS data were obtained as part of the Japanese national cancer genomics initiative (Genome Research in Cancers and Rare Diseases; G-CARD). A sequencing library was generated using the Illumina DNA PCR-Free Prep Tagmentation Kit according to the manufacturer’s protocol, and sequenced using the NovaSeq 6000 (Illumina) with a 150-bp paired-end read protocol at a target depth of 100× for tumours and 30× for normal controls.
The reads were aligned to the human hg38 reference genome by Parabricks. Somatic variants were called in tumour–normal paired mode using Mutect2 (GATK v.4.5.0.0), retaining variants with tumour log-odds (TLOD) ≥ 20 and variant allele frequency ≤ 0.3 in the matched normal. Variants listed in germline databases, such as gnomAD (v.3.1.2), the 1000 Genomes Project and Tohoku Medical Megabank Organization (ToMMo), with a minor allele frequency of 0.001 or higher, were removed. SVs were called using GRIDSS v.2.12.0 and Genomon SV v.0.8.0. Filtered outputs from GRIDSS (‘FILTER = PASS, QUAL ≥ 500, AS > 0, or RAS > 0’) and Genomon SV (‘overhang ≥ 150’) were combined. CNAs were inferred from tumour–normal paired WGS data using Battenberg (v.3.0.0) with default parameters, based on log R ratios and B-allele frequencies of germline heterozygous SNPs.
Bulk RNA-seq experiments and analysis
Libraries for RNA-seq were prepared using the NEBNext Single Cell/Low Input RNA Library Prep Kit for Illumina (New England BioLabs) and were subjected to sequencing using the NovaSeq 6000 (Illumina) with a paired-end protocol.
The sequencing reads were preprocessed by fastp55 with ‘–detect_adapter_for_pe -q 15 -n 10 -u 40’ parameters, and aligned to the reference genome using STAR56. Reads on each gene defined in the University of California, Santa Cruz (UCSC) hg19 gene annotation were counted with featureCounts57. The quality of sequencing data was assessed using mapping and count statistics. Samples were excluded from the analysis if any of the following criteria were met: ‘uniquely_mapped_percent’ (STAR) < 30%, ‘percent_assigned’ (featureCounts) < 30% or ‘assigned’ (featureCounts) < 3 million. The edgeR package58 was used to normalize read counts and calculate counts per million (CPM) values of genes. Genes with CPM > 1 in at least two samples and located on the autosomal chromosomes were kept and used for downstream analysis. To adjust for cohort-specific technical effects while preserving biological variation, batch correction between the Swedish and Japanese cohorts was performed using linear modelling in limma59. Differentially expressed genes (DEGs) for each subgroup were identified using the eBayes test in limma through a one-versus-rest comparison, with thresholds of FDR < 0.05 and |log2-transformed fold change| > 0.5. The top 3,000 DEGs in each subgroup are provided in Supplementary Table 8. GSEA analysis was performed on genes ranked by fold change using the GSEA function in the clusterProfiler package60 (minGSSize = 20, pAdjustMethod = “BH”) and a curated set of genes associated with haematopoietic cells and AML derived from the Human Molecular Signatures Database (MSigDB)61 (Supplementary Table 10). Gene fusions were detected using the Genomon fusion pipeline (https://github.com/Genomon-Project/fusionfusion) and filtered for known drivers of AML.
Bulk ATAC-seq experiments
ATAC-seq experiments were performed using the Fast-ATAC protocol21,62. Cryopreserved tumour cells were thawed, and 50,000 cells were pelleted. Fifty microlitres of transposase mixture (comprising 25 µl 2× TD buffer, 2.5 µl TDE1 (Illumina, FC-121-1031), 0.5 µl 1% digitonin (Sigma-Aldrich, D141) and 22 µl nuclease-free water) was added to the cells. After transposition reactions at 37 °C for 30 min, transposed DNA was purified using the QIAGEN MinElute Reaction Cleanup Kit or Sera-Mag Select (Cytiva) magnetic beads, and PCR-amplified using the NEBNEXT Q5 Hot Start HiFi PCR Master Mix and custom primers63. The TapeStation 4200 (Agilent) with High Sensitivity D5000 ScreenTape was used to assess fragment size and confirm nucleosomal periodicity, a characteristic feature of ATAC-seq libraries20,64. The resulting library was sequenced using the NovaSeq 6000 (Illumina) with a paired-end protocol.
Bulk ATAC-seq analysis
Reads were aligned to the reference genome using Bowtie265 with ‘-X 2000 –no-mixed –very-sensitive’ parameters, after adapter trimming using Skewer66. The quality of sequencing data was assessed using ataqv67 and ATACseqQC68. Samples were excluded from the analysis if any of the following criteria were met: fewer than six million uniquely mapped reads; mapping rate lower than 40%; more than 30% of reads on mitochondrial chromosome; or transcription start site (TSS) enrichment scores lower than 2 in both 1,000 and 2,000 bp windows. After removing duplicates and reads on the mitochondrial genome or blacklisted regions (ENCODE), peaks were called using HMMRATAC69 with the ‘–window 2500000’ parameter and peaks were fixed to a width of 501 bp centred on the peak summit64. When extended or merged peaks overlapped, the region with the highest peak score was retained. For each of the six datasets—including our own datasets (Swedish AML, Japanese AML, AML cell lines and normal bone marrow) as well as public datasets21 (sorted normal blood and AML cells)—we merged peaks from all samples in each dataset. Peaks that overlapped in two or more samples were considered the recurrent peak set for each dataset. For samples from patients with AML, Swedish and Japanese peak sets were further merged and filtered for those recurrently identified in both the Swedish and the Japanese samples. We then merged five peak sets generated as above (samples from patients with AML (eCHROMA), AML cell lines, normal bone marrow, public normal blood and public AML) to generate a union peak set fixed to a width of 501 bp, which was used for downstream analysis (Extended Data Fig. 1b). Annotation of peaks was done using the annotatePeaks function in HOMER70. Reads on peaks were counted and normalized to calculate CPM values using edgeR58. Differentially expressed ATAC peaks (DEPs) for each subgroup were identified using the eBayes test in limma through a one-versus-rest comparison, with thresholds of FDR < 0.05 and |log2-transformed fold change| > 0.5. The top 3,000 DEPs in each subgroup are provided in Supplementary Table 9. Inference of cellular contribution to each ATAC-seq data was evaluated by the CIBERSORTx27 with CPM values of our AML data and public data for 13 normal blood cell types21 as input files.
Clustering analysis using bulk ATAC-seq
ATAC log2 CPM were quantile normalized using preprocessCore (https://github.com/bmbolstad/preprocessCore), followed by the exclusion of peaks on chromosomes X and Y. Batch correction was applied as described above. To focus on leukaemia-intrinsic chromatin variation and minimize contamination from non-malignant cells, peak variance was calculated using samples with at least 75% bone-marrow blasts. The 3,000 most variable peaks were selected to capture dominant inter-sample chromatin heterogeneity while reducing noise from low-variance regions. The first 50 principal components were determined by prcomp in R, a nearest-neighbour graph was generated using buildSNNGraph in the scran package71 with a ‘k = 7’ parameter and Leiden clustering72 was performed using ‘cluster_leiden’ from the igraph package with parameters of ‘resolution = 0.2’ and ‘n_iterations = 100’. Clusters with less than 1% of samples (fewer than 16 out of 1,563 cases) were reassigned to the cluster of the nearest sample by calculating the squared Euclidean distance between samples in the two-dimensional UMAP space. ATAC subgroups were annotated with representative genetic alterations and/or differentiation states (for example, A, PML::RARA) to aid interpretability. Labels were assigned on the basis of enriched and relatively specific features within each subgroup, including mutation patterns and differentiation signatures, but do not indicate that subgroup identity is defined solely by these alterations. For high-resolution clustering, we used a resolution of 0.3 with all other parameters kept the same.
Clustering stability
To evaluate the robustness of the Leiden clustering results, we performed a resampling-based stability analysis. First, we repeatedly subsampled 90% of the samples and re-applied the same clustering parameters (resolution = 0.2, k = 7 neighbours) for 100 iterations, calculating the adjusted Rand index between the original and each resampled clustering to quantify overall agreement; adjusted Rand index values close to 1 indicate strong stability. Second, we computed a sample-wise consistency score for each sample, defined as the fraction of iterations in which the sample clustered together with its original cluster members, conditional on co-occurrence, thus measuring the stability of cluster membership at the individual sample level. Third, we constructed a cluster–cluster similarity matrix that captures the probability that samples from different clusters co-cluster across all resampling runs, providing a summary of within-cluster consistency and potential cross-cluster mixing.
Decision-tree modelling based on genomic alterations
To test whether ATAC-defined subgroups could be explained by simple combinations of genomic alterations, we constructed one-versus-rest decision-tree classifiers using gene mutations, CNAs and SVs as binary input features. For each of the 16 ATAC subgroups, a binary classifier was trained using a CART framework implemented in the rpart package in R. Rare alterations (present in fewer than four positive samples) were excluded from model training. Samples were randomly split into training (80%) and test (20%) sets. Tree depth was fixed at a maximum depth of 3 to constrain model complexity and emphasize interpretability. The splitting criterion was Gini impurity. Class predictions were derived using a fixed probability threshold of 0.5. Model performance was evaluated on the independent test set using sensitivity, positive predictive value (PPV), F1 score and balanced accuracy. Model complexity was quantified as the number of internal splits in the final tree. To examine how performance changed with increasing model complexity, we also trained models across a range of maximum tree depths (2–12), while keeping all other parameters fixed.
Bulk-RNA-seq-based prediction of ATAC subgroups
Normalized gene-count data for four external adult AML cohorts were obtained from a previous study10. The TARGET cohort was excluded from the analysis because it consisted mainly of paediatric AML. The count matrix was quantile normalized and merged with our RNA-seq count matrix, followed by batch correction using ‘removeBatchEffect’. Prediction models for ATAC subgroups were generated using the ClaNC algorithm, a nearest-centroid classifier that ranks genes by standard t-statistics and selects subgroup-specific genes30. Fifteen genes per subgroup were selected for subgroup-specific markers and incorporated into the prediction model. Gene lists used for each prediction model are summarized in Supplementary Table 11. The accuracy of prediction models was calculated by fivefold cross-validation, separating 80% of samples for training and the remaining 20% for validation.
DNA methylation experiments and analysis
DNA methylation profiles were analysed using the Infinium MethylationEPIC BeadChip Kit according to the manufacturer’s protocol. The raw files were processed using the ChAMP package in R, filtering probes and generating the normalized β-value matrix. Common probes across samples were used for downstream analysis. To correct batch effects between different cohorts and probe sets, the ‘removeBatchEffect’ command from the limma package59 was used.
For comparison between chromatin accessibility and DNA methylation, ATAC peak regions overlapping with methylation probes were analysed. DNA methylation levels in a given region were calculated as the average β-values of all probes in the region. For each ATAC subgroup, average DNA methylation levels were computed, and the top 3,000 regions with the most variable DNA methylation levels across subgroups were visualized in heat maps.
GRN analysis by bulk ATAC-seq and RNA-seq
GRNs were inferred using the ANANSE software35. In brief, this software integrates chromatin accessibility and gene-expression data to infer enhancer-based GRNs. Input data included the consensus ATAC peak set, a merged ATAC-seq BAM file, mean RNA expression (log2 CPM) and the default motif database (GimmeMotifs73).
For each TF, ANANSE estimates genome-wide binding potential, predicts target genes and incorporates expression levels of the TF and its targets to model regulatory interactions. The regulatory importance of each TF is quantified by two key metrics: out-degree (the number and strength of predicted regulatory connections) and link score (the strength of each individual connection). GRNs were constructed for each ATAC subgroup and for the entire AML cohort. To identify subgroup-specific regulators, we compared each subgroup GRN with the global AML GRN and calculated TF influence scores, which measure how much a TF explains differences in expression between the two groups (Extended Data Fig. 6a). The top 20 differential TFs per subgroup were visualized.
ChIP–seq experiments and analysis
ChIP–seq experiments were performed according to the SimpleChIP Plus Sonication Chromatin IP Kit (Cell Signaling Technology)62 with minor modifications. Cryopreserved cells were thawed, and more than one million cells were fixed with 1% formaldehyde (Thermo Fisher Scientific) in phosphate-buffered saline (PBS) for 10 min at room temperature with gentle mixing. The reaction was stopped by adding glycine solution (10×) (Cell Signaling Technology) and incubated for 5 min at room temperature, and the cells were washed twice in cold PBS. The cells were then processed with the SimpleChIP Plus Sonication Chromatin IP Kit (Cell Signaling Technology) and Covaris E220 (Covaris) according to the manufacturer’s protocol. The antibodies used for ChIP were as follows: SMC1 (Abcam, ab9262), CTCF (Cell Signaling Technology, D31H2), RPB1 (Cell Signaling Technology, D8L4Y), H3K27ac (Cell Signaling Technology, D5E4) and H3K27me3 (Cell Signaling Technology, C36B11). After purification of the precipitated DNA, libraries were constructed using the ThruPLEX DNA-seq Kit (Takara) as per the manufacturer’s protocol, and subjected to sequencing using the NovaSeq 6000 (Illumina). ChIP–seq experiments were performed with input controls. The sequencing reads were aligned to the reference genome using Bowtie74, after adapter trimming with Skewer66 and read tail trimming to a total length of 50 bp using Cutadapt75. The quality of sequencing data was assessed using ‘plotFingerprint’ in deepTools76. Samples were excluded from the analysis if any of the following criteria were met: ‘X-intercept’ (deepTools) > 0.85 or ‘Synthetic JS Distance’ (deepTools) < 0.225. After removing duplicates and reads on blacklisted regions (ENCODE)77, peaks were called using MACS278 and a P value threshold of 1 × 10–3 with an input control for each sample. For each ChIP–seq (SMC1, CTCF, RPB1, H3K27ac and H3K27me3), peaks were merged for all AML samples, and recurrently identified peaks were regarded as a consensus peak set.
SE analysis
To identify SEs, recurrent enhancers were first identified in all AML samples using H3K27ac ChIP–seq data. Identified enhancers were stitched and ranked with H3K27ac ChIP–seq and input data, using ROSE34 with a ‘-t 2500’ parameter. Mean signals were used to calculate enhancer ranks across all AML samples. To identify SEs for each ATAC subgroup, mean signals were calculated to identify the top 750 enhancers in each ATAC subgroup, which were regarded as SEs. SEs were separated into three categories on the basis of their distribution across subgroups: common SEs (present in 12 or more subgroups), partially shared SEs (present in 2 to 11 subgroups) and unique SEs (specific to a single subgroup). Subgroup-specific SEs were defined as those found in each subgroup and categorized as partially shared or unique SEs. Known SEs for various cell types and cancers were obtained from previous reports34. Annotation of SEs was done using annotatePeaks in HOMER70, filtering for genes expressed in our AML cohort (logCPM > 1 in more than 1% of patients). AML driver genes and TF genes were determined using databases, including the Catalogue of Somatic Mutations in Cancer (COSMIC) Cancer Gene Census (CGC) (as of 6 November 2024)79, the database of leukaemia gene literature (dbLGL)80, a list of TFs from a previous study81 and the JASPAR21,82 database, and were manually selected. A list of SEs identified in this study is provided in Supplementary Table 6.
SE-regulated gene signatures
Enriched gene ontologies in SE-associated genes were identified using enricher in the clusterProfiler package60 (pAdjustMethod = “BH”, qvalueCutoff = 0.25) and ontology geneset (‘hallmark’, ‘c2.cp.reactome’ and ‘c5.go.bp’) from MSigDB (v.2024.1)61. Expression levels of SE-associated genes in each haematopoietic cell type were calculated using public gene-expression data from DMAP83, and the average expression levels were determined for each cell type.
SE-based TF network analysis
To analyse SE-regulated TF networks, we applied the Coltron software36 with ROSE34 outputs generated from the merged H3K27ac BAM files for each ATAC subgroup. This software computes the inward binding (in-degree) of other SE-associated TFs to a given SE-associated TF, as well as the outward binding (out-degree) of the TF to other SEs. The Coltron score for each TF was determined as the sum of its in-degree and out-degree (Extended Data Fig. 6c).
scRNA/ATAC-seq experiments
Single-cell matched RNA-seq and ATAC-seq experiments were performed using the Next GEM Single Cell Multiome ATAC + Gene Expression Reagent Kit (10x Genomics), according to the manufacturer’s protocols (CG000365 for nuclei isolation and CG000338 for library generation). Cryopreserved cells were thawed, dead cells were stained with DAPI and live mononuclear cells were sorted using FACS Aria III (BD Biosciences). The samples were resuspended in lysis buffer (10 mM Tris-HCl (pH 7.4), 10 mM NaCl, 3 mM MgCl2, 1% bovine serum albumin (BSA; Miltenyi Biotec, 130-091-376), 0.1% Tween-20 (Bio-Rad, 1610781), 0.1% IGEPAL (Sigma-Aldrich, i8896), 0.01% digitonin (Thermo Fisher Scientific, BN2006) and 1 mM DTT (Sigma-Aldrich, 646563), plus 1 U μl−1 RNase inhibitor (Thermo Fisher Scientific, 10777019)), incubated on ice for 3 min, washed three times with wash buffer (10 mM Tris-HCl (pH 7.4), 10 mM NaCl, 3 mM MgCl2, 1% BSA, 0.1% Tween-20 and 1 mM DTT, plus 1 U μl−1 RNase inhibitor) and passed through a 40-μm Flowmi Cell Strainer (Bel-Art). After microscopy inspection and counting of nuclei using a Countess II FL Automated Cell Counter (Thermo Fisher Scientific), nucleus suspensions were prepared in a concentration targeting a maximum of 10,000 nuclei recovery, and incubated with transposase to add adapter sequences to the DNA fragments. The suspensions containing transposed nuclei were subjected to gel bead in emulsion (GEM) generation, incubation, and clean-up, using the Chromium Next GEM Chip J Single Cell Kit (10x Genomics) and Chromium Controller. The resulting suspensions contained ATAC fragments and cDNA with the same cell barcodes. Pre-amplification of cDNA was performed, and the amplified product was split and used as the input for ATAC and gene-expression library construction. Libraries were generated using the 10x Genomics Single Index N Set for ATAC and the 10x Genomics Dual Index TT Set A for RNA. The scATAC libraries were sequenced using the DNBSEQ-G400 (MGI) with a custom protocol (read 1: 50 cycles, read 2: 49 cycles, i5 Index: 24 cycles, i7 Index: 8 cycles). The scRNA libraries were sequenced using the DNBSEQ-G400 (MGI) with a custom protocol (read 1: 28 cycles, read 2: 90 cycles).
scRNA/ATAC-seq analysis
Cell Ranger ARC (10x Genomics)84 was used for data processing to generate BAM files and count matrices for scRNA and scATAC with the reference genome and GENCODE human v.19 annotation85. The report generated by Cell Ranger was manually evaluated and samples were filtered to remain with no errors or with only warnings. Basic quality-control reports including analysed cell numbers from Cell Ranger for each sample are summarized in Supplementary Table 7. ArchR86 was used to filter for cells with the following thresholds: number of unique molecular identifiers > 100 and proportion of mitochondria genes < 0.05 for scRNA; TSS enrichment score > 4 and number of fragments > 1,000 for scATAC. Doublet cells were removed by the ‘filterDoublets’ function in ArchR with a filterRatio of 1. The gene-expression count matrix was stored in the Seurat87 object, and mitochondrial genes were excluded from the following analysis. The consensus peak set generated from bulk ATAC-seq was filtered for peaks on autosomal chromosomes and used to count reads on peaks in the scATAC analysis, using the FeatureMatrix function from the Signac package88.
Merging individual data, normalization and clustering for scRNA/ATAC-seq
Raw count matrix data for single-cell gene expression and chromatin accessibility was written to disk using the BPCells package (https://github.com/bnprks/BPCells) to enable high-throughput data processing and merged across samples. The merged gene-count matrix was normalized for sequencing depth using the SCTransform function in Seurat87. scRNA-based clustering was done using the first 50 principal component analysis (PCA) dimensions with the FindNeighbors and FindClusters functions (resolution = 0.15, algorithm = 1). scRNA UMAP was computed using the first 50 PCA dimensions with RunUMAP (n.neighbours = 30, min.dist = 0.1). The merged count matrix for scATAC was normalized using the term frequency inverse document frequency (TF-IDF) normalization function in BPCells. Cells were clustered on the basis of scATAC using the PCA dimensions 2 to 20 with knn_hnsw (ef = 3000), knn_to_snn_graph and cluster_graph_louvain (resolution = 0.15) functions in BPCells. Each scATAC cluster was classified as AML-predominant if it met both of the following criteria: (i) the average proportion of cells derived from remission samples was less than 10%; and (ii) at least one AML subgroup contributed 25% or more of its cells to the cluster. Clusters not satisfying these criteria were classified as normal-mixed. scATAC UMAP was computed using the PCA dimensions 2 to 20 with the umap function in the uwot package (n.neighbours = 30, min.dist = 0.1).
Estimation of differentiation, pseudotime and LSC scores using scRNA-seq
The BoneMarrowMap package in R was used to project AML cells onto reference scRNA-seq data of human bone-marrow haematopoiesis41. Cells with low mapping quality were excluded from the analysis according to the following criteria: (1) mean absolute deviation of mapping error scores ≧ 2; (2) assigned to ‘orthochromatic erythroblast’ and AUC score of haemoglobin genes < 0.2, calculated using the AUCell package. Pseudotime analysis was performed using the ‘predict_Pseudotime’ function in BoneMarrowMap. The LSC score was computed as the AUC for gene set enrichment of LSC signatures, defined by DEGs specific to sorted LSC+ fractions41,42.
Single-cell GRN analysis
TF activity was inferred using SCENIC+43, which integrates matched scRNA-seq and scATAC-seq data to construct GRNs at single-cell resolution. The analysis followed the standard SCENIC+ pipeline with minor modifications (Extended Data Fig. 8a). We used a multiome mode with 5 cells per metacell, and empirically set the number of topics to 100. The search space for regulatory elements was defined as 0–500 kb from each gene. Established eRegulons were filtered with the following parameters: ‘rho_threshold’ = 0.03, ‘min_regions_per_gene’ = 0 and ‘min_target_genes’ = 10. For group-level comparisons, cells were downsampled to 2,000 per ATAC subgroup when calculating region-based and gene-based specificity scores. For pseudotime analysis, cells were downsampled to 300 per pseudotime bin in each ATAC subgroup if more than 300 cells were present in that bin. Direct positive eRegulons (those with positive TF-to-gene and region-to-gene links) were used for downstream analyses. Extended eRegulons were included only when direct ones were unavailable.
Drug sensitivity screening in samples from patients with AML
The procedure for drug sensitivity and resistance testing in the samples used for the study has been described in detail previously46. Biobanked mononuclear cells were thawed and added to pre-spotted drug plates (FIMM HTB)89. The mononuclear cells were incubated in HS-5 conditioned (12.5%; ATCC) complete RPMI (10% FBS (Thermo Fisher Scientific), 2 mM l-glutamine (Sigma-Aldrich), 100 IU ml−1 penicillin and 0.1 mg ml−1 streptomycin (Pen-Strep; Sigma-Aldrich)), for 72 h (37 °C, 5% CO2). Cell viability was measured by CellTiter-Glo (CTG; Promega) on an EnSight plate reader (PerkinElmer). Drug sensitivity scores (DSSs) were calculated with Breeze90. Selective DSS (sDSS) values were calculated by subtracting the DSSs of healthy bone-marrow control samples from the DSSs of samples from patients with AML.
Survival analysis
Survival analysis was performed for patients who were treated with standard intensive chemotherapy, and observations were censored at the last follow-up. For clinical parameters, the median values were used as the threshold unless otherwise specified. Overall survival was estimated using the Kaplan–Meier method, and differences between groups were evaluated by the log-rank test.
To assess the prognostic relevance of deconvolution patterns, we performed feature selection using LASSO-penalized Cox proportional hazards regression (glmnet R package, alpha = 1), incorporating the estimated proportions of 13 haematopoietic cell populations derived from ATAC-seq–based deconvolution as explanatory variables, with overall survival as the outcome. The optimal penalty parameter (lambda.min) was selected by cross-validation. The resulting coefficients were used to compute a LASSO-based risk score for each patient, defined as the weighted sum of the corresponding cell fractions.
To determine whether ATAC subgroups provided extra prognostic information beyond established genomic risk categories, multivariable Cox proportional hazards models were fitted separately within each ELN 2022 risk group. These models included age, sex, WBC, haemoglobin level, platelet count, blast percentage and ATAC subgroup membership as covariates. Subgroups that were significantly associated with worse overall survival within each ELN stratum were designated as high-risk ATAC subgroups. We stratified patients in each ELN category into two groups according to the presence or absence of these risk ATAC subgroups and compared overall survival in Kaplan–Meier analyses. To evaluate the additive prognostic value of ATAC subgroups, we built a Cox model for ELN 2022 risk classification and compared it with a combined model of ELN and ATAC risk groups (Fig. 5b). Model performance was assessed by calculating the C-index using 1,000 bootstrap replicates in both the training (Swedish) and the validation (Japanese) cohort. The incremental prognostic value of ATAC subgroups was quantified as the difference in C-index between models, assessed across the same bootstrap replicates.
Statistical analysis
Experimental replication was not performed. The robustness of the findings is supported by the large cohort size (n = 1,563) and the consistency of results across several independent analytical approaches, data types and independent cohorts. Statistical analyses were performed in R. Comparisons between groups were based on the two-sided Wilcoxon rank-sum test for continuous data and the Fisher’s exact test for categorical data, unless otherwise specified. Explained variances in several features by ICC, WHO and ATAC classifications were evaluated as R2 values calculated by the adonis function in R, with Euclidean distances. For gene-expression profiles, the first 50 principal components were used as input.
Reporting summary
Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.
