We downloaded a scRNA-seq dataset, GSE116256 (van Galen et al., 2019), and two bulk gene expression datasets, GSE6891 (Verhaak et al., 2009) and TCGA-LAML (Cancer Genome Atlas Research et al., 2013), for AML from the Gene Expression Omnibus (GEO) data repository (RRID:SCR_005012) and Genomic Data Commons data portal (RRID:SCR_014514), respectively. The scRNA-seq dataset contains single-cell GEPs and cell annotations of 30,712 cells and 27,899 genes from the bone marrow aspirates of 16 AML patients. The cell annotations comprise information such as the number of unique molecular identifiers (UMIs), the number of expressed genes, and the inferred cell type for each cell. A total of 21 cell types were defined, including HSC, HSC-like, progenitor (Prog), Prog-like, granulocyte–monocyte–progenitor (GMP), GMP-like, promonocyte (ProMono), ProMono-like, monocyte (Mono), Mono-like, conventional dendritic cell (cDC), cDC-like, plasmacytoid dendritic cell (pDC), early erythroid progenitor (earlyEry), late erythroid progenitor (lateEry), progenitor B cell (proB), mature B cell (B), plasma cell (plasma), naïve T cell (T), cytotoxic T lymphocyte (CTL), and natural killer cell (NK). Details of the scRNA-seq dataset can be learned from elsewhere (van Galen et al., 2019). For the bulk gene expression dataset GSE6891, 537 GEPs of AML patients profiled by microarray were obtained. For TCGA-LAML, 151 GEPs with fragments per kilobase million normalization were downloaded. The corresponding clinical characteristics and survival information for each sample were downloaded from the cBioPortal database (RRID:SCR_014555).

The workflow of this study is illustrated in Figure 1. We first constructed the GES reference matrix of the 21 cell types required in CIBERSORT (Newman et al., 2015) (RRID:SCR_016955) using AML scRNA-seq profiles. The CTCs of patients in the bulk gene expression datasets of GSE6891 and TCGA-LAML were subsequently estimated. A CTC-based prognostic model was established, with GSE6891 as the training set, and was validated in TCGA-LAML subsequently.

For the scRNA-seq dataset GSE116256, we excluded cells derived from samples of AML314, AML371, AML722B, and AML997 due to the unconfident cell type annotations (van Galen et al., 2019). Next, we computed the ratio of UMI counts to the number of expressed genes for each cell, termed UTG ratio below. In each cell type, cells with outlier values of UTG ratio were prone to be low in quality. The threshold to filter such cells was determined to be the median UTG ratio plus–minus three times the median absolute deviation (Leys et al., 2013). A total of 27,023 cells remained (Supplementary Figure S1).

The bulk GEPs of GSE6891 were generated by Affymetrix Human Genome U133 Plus 2.0 Array (Verhaak et al., 2009). The raw CEL files were processed using affy (version 1.66.0) and normalized by the Gene Chip Robust Multi-array Average (Wu et al., 2004) algorithm using gcrma (version 2.60.0) Bioconductor R package. The probe set IDs were transformed to the corresponding gene symbols according to the chip definition file (GEO accession: GPL570). The probe sets that did not match any gene symbols or matched multiple gene symbols were filtered out. To retain enough genes for subsequent analysis, we computed the mean expression of probe sets that matched the same gene and chose the probe set with the highest average gene expression to represent that gene (Ng et al., 2016). Among the cases in GSE6891, we only retained de novo AML cases whose age at diagnosis were greater or equal to 18 with completed survival information.

For TCGA-LAML, the ensemble gene IDs of the downloaded GEPs were transformed to gene symbols according to the comprehensive gene annotation files of GENCODE release 38 (GRCh38.p13; RRID:SCR_014966) in gene transfer format. We filtered out the ensemble gene IDs matching the same gene symbol due to the difficulty in determining which ensemble gene ID to represent that gene. Among the cases in TCGA-LAML, we took the same filtering criteria as implemented for the GSE6891 dataset.

We constructed the cell type-specific GES reference matrix based on the AML scRNA-seq GEPs using Seurat (Stuart et al., 2019) (version 3.2.0) R package. First, the single-cell GEPs of AML patients were integrated and imported into a Seurat object. All cells were labeled as the cell type in the annotation file. Then, we normalized the UMI counts to counts per million (CPM) and performed natural-log transformation [log (CPM+1)]. Subsequently, we conducted the differential expression analysis using FindAllMarkers function to acquire highly expressed genes of each cell type by comparing the cells of 1 cell type against all others in turn. The tests of comparisons between groups used the “bimod” method, a likelihood-ratio test for single-cell gene expression (McDavid et al., 2013). The “min.pct” parameter was set to 0. Other parameters were set as default. The acquired highly expressed genes of each cell type with the adjusted p-values lower than 0.05 and the average natural-log fold-change (logFC) above 1 were retained (Supplementary Figure S2 and Supplementary Table S1). Notably, the highly expressed genes selected to build the GES reference matrix are the dominant influence factor for CTC estimations, thereby affecting subsequent modeling. Therefore, we extracted the top 25, 50, 100, and 150 most significantly highly expressed genes for each cell type and computed the mean expression by cell type to build 4 cell type-specific GES reference matrices (GES25, GES50, GES100, and GES150; Supplementary Figure S3, Supplementary Table S2, Supplementary Table S3, Supplementary Table S4, and Supplementary Table S5) (Donovan et al., 2020).

We performed a simulation analysis to examine the accuracy of CIBERSORT using the custom cell type-specific GES reference matrices. Specifically, we first generated 100 artificial samples using scRNA-seq profiles. For each sample, we selected a random number of cells for each cell type from at least 50 to the maximum number of cells for that cell type through the cell barcodes (Donovan et al., 2020). The normalized GEPs of these cells were summed to create the artificial sample with known cell type compositions. Subsequently, we ran CIBERSORT on these artificial samples using different GES matrices. Additionally, two other deconvolution methods, MuSiC (Wang et al., 2019) and MOMF (Sun et al., 2019), were also used for comparisons. The Pearson correlation coefficients of the real proportions and the estimated proportions were computed by each cell type as the metric of accuracy.

The simulation results showed that the performances of CIBERSORT, MuSiC, and MOMF were similar (Supplementary Figure S4). However, we noticed that MuSiC and MOMF took a much longer running time and much more memory consumptions (data not shown). Accordingly, we chose CIBEROSRT to estimate the relative proportions of 21 AML cell types for the bulk gene expression datasets GSE6891 and TCGA-LAML, setting 100 permutations and disabling the quantile normalization option.

After estimating the CTCs (Supplementary Table S6, Supplementary Table S7, Supplementary Table S8, and Supplementary Table S9), we found that the estimated proportions of some cell types were almost 0 for most of the samples, probably due to estimation error. To reduce the influence on subsequent modeling, we converted the cell types whose mean proportions were lower than 0.05 or median proportions were equal to 0 to dichotomous variables, with 0 as the cutoff value. Cell types converted to dichotomous or that remained continuous in both datasets and whose Pearson correlation coefficient was r > 0.8 in the simulations were used to train and validate the prognostic model.

The bulk gene expression dataset GSE6891 was set as the training set, and TCGA-LAML was set as the validation set to establish and validate a novel prognostic model for AML based on CTCs. With OS as the survival outcome, we performed the least absolute shrinkage and selection operator (LASSO) Cox regression (Simon et al., 2011) and 10-fold cross-validation using glmnet (version 4.1–1) R package. To obtain a robust model, we repeated this process 100 times using different random seeds, and cell types with non-zero coefficients in at least 95 fittings were retained. The coefficients of 100 fitting processes for the retained cell types were averaged as the final coefficient (Elsayed et al., 2020). The linear combination of the selected cell types in the LASSO Cox regression model weighted by the coefficients served as the prognostic marker for AML, called CTC score. For better interpretation and visualization, we partitioned all patients into low- and high-CTC-score groups by median.

The established CTC score was validated in TCGA-LAML. We computed the CTC scores for patients in TCGA-LAML based on the linear equation above (Supplementary Table S10). We likewise partitioned the patients in the validation set into low- and high-CTC-score groups based on the median. Kaplan–Meier curves were used to display the different prognoses between low- and high-CTC-score groups.

We considered displaying the CTC score established on the CTCs estimated with GES100 as the reference matrix to be the main results. Other prognostic models based on the CTCs estimated using GES25, GES50, and GES150 were considered as sensitivity analysis and could be accessible in Supplementary Figure S5. The Harrell’s concordance index (C-index) was used to compare the performance of these models (Harrell et al., 1996).

We found that GMP-like has a great weight when computing the CTC score (see results part). It has been reported that GMP-like is highly associated with two abnormal karyotypes (i.e., PML-RARA and RUNX1-RUNX1T1), both of which indicate a favorable prognosis (Appelbaum et al., 2006; Wang and Chen., 2008; van Galen et al., 2019). Thus, it is crucial to verify whether the prognostic significance of CTC score was dominantly captured by existing prognostic factors such as karyotypes and cytogenetic risk classifications. To verify this point, we first implemented univariable Cox regressions for clinical characteristics. The clinical characteristics significant in both training and validation dataset and CTC score were introduced to multivariable Cox regressions using survival (version 3.2–7) R package.

We further evaluated the performance of CTC score by comparing it with the LSC17 score and APS. The LSC17 score was constructed by the expression of 17 genes highly expressed in LSCs, while the APS was constructed by the expression of 16 genes acquired by LASSO Cox regression (Ng et al., 2016; Docking et al., 2021). The LSC17 score and APS for patients in the validation set TCGA-LAML were computed in compliance with the data processing flow and calculation equation according to the original articles (Supplementary Table S10) (Ng et al., 2016; Docking et al., 2021). Considering the comparability, all three prognostic scores were not converted to dichotomous variables. We implemented the time-dependent receiver operating characteristic (ROC) curve analysis to evaluate and compare the predictive accuracy using area under the ROC curve (AUC) as the indicator. The predictive sensitivities and specificities of CTC score, LSC17 score, and APS at 1-, 2-, 3-, and 5-years timepoints were computed and compared using timeROC (Blanche et al., 2013) (version 0.4) R package.

For the clinical characteristics of patients in the bulk gene expression datasets GSE6891 and TCGA-LAML, continuous variables were described by medians and ranges, and categorical variables were described by frequencies and proportions. We used the Wilcoxon test or Kruskal–Wallis test for group comparisons of continuous variables and the chi-square test or Fisher’s exact test for that of categorical variables. All statistical tests were two-tailed, and p-values lower than 0.05 were considered statistically significant. All the analyses were performed in R-4.0.2.

For the bulk gene expression dataset GSE6891, 11 patients whose age at diagnosis was lower than 18, 17 patients of myelodysplastic syndrome, and four patients with missing survival information were filtered out. Eventually, 429 patients were eligible, whereas all patients in TCGA-LAML passed the filtering criteria. The descriptive characteristics of patients in these two datasets are shown in Table 1. Patients in GSE6891 were younger than those in TCGA-LAML (p < 0.0001). FAB classification (p = 0.0010) and cytogenetic risk (p = 0.0313) were also different between GSE6891 and TCGA-LAML. Patients in GSE6891 comprises more FAB-M5 subtype (23.3% in GSE6891 vs 9.9% in TCGA-LAML) and less poor cytogenetic risk strata (19.3% in GSE6891 vs 23.8% in TCGA-LAML). The CTCs for patients in the GSE6891 and TCGA-LAML datasets estimated with GES100 as the reference matrix are shown in Supplementary Figure S6.

The median follow-up time of patients in the bulk gene expression datasets GSE6891 and TCGA-LAML was 20.11 months [interquartile range (IQR), 7.89–92.78 months] and 19 months (IQR, 6.45–42.1 months), respectively. We fitted a LASSO Cox regression model and defined the CTC score computed by the following equation:

CTC score = (−1.7016 × GMP-like) + (0.2015 × HSC-like) + (−0.293 × T), where HSC-like and T were dichotomous. The negative coefficient of GMP-like indicated that lower relative proportions of GMP-like at diagnosis would predict worse survival outcomes. The estimated HSC-like greater than 0 and T equal to 0 would predict worse prognoses.

Comparing with the low-CTC-score group, the high-CTC-score group showed a 1.57-fold (95% CI, 1.23 to 2.00; p = 0.0002) higher overall mortality risk in the training set GSE6891 and 2.32-fold (95% CI, 1.53 to 3.51; p < 0.0001) in the validation set TCGA-LAML (Table 2). The 5-years OS rate for GSE6891 was 47.7% (95% CI, 41.4–54.9%) in the low-CTC-score group and 31.1% (95% CI, 25.5–37.9%) in the high-CTC-score group. For TCGA-LAML, the 5-years OS rate was 41.2% (95% CI, 29.7–57.1%) and 17.7% (95% CI, 10.2–30.7%) in the low-CTC-score group and high-CTC-score group, respectively (Figure 2).

The individual-level results of CTCs estimated using GES25, GES50, GSE100, and GES150 could be obtained in Supplementary Table S6, Supplementary Table S7, Supplementary Table S8, and Supplementary Table S9. As displayed in Supplementary Figure S5, the CTC-based scores established by reference matrices with different GES matrices were robustly associated with the OS of AML in the validation set, with C-index ranging from 0.64 (95% CI, 0.58–0.70) to 0.67 (95% CI, 0.61–0.73).

We performed univariable and multivariable Cox regressions in both the training and validation sets to test whether the CTC score is an independent factor associated with the OS for AML in adults. Among the clinical characteristics, age at diagnosis, cytogenetics risk, and karyotype were significantly associated with OS in both datasets (Table 2). The multivariable Cox regression results showed that CTC score remained statistically significant in GSE6891 (HR = 2.25; 95% CI, 1.20 to 4.24; p = 0.0119) and TCGA-LAML (HR = 7.97; 95% CI, 2.95 to 21.56; p < 0.0001) when adjusting for age at diagnosis, cytogenetic risk, and karyotype (Figure 3). These results suggested that CTC score can predict the prognosis of AML independent of age at diagnosis, cytogenetic risk, and karyotype.

In TCGA-LAML, we evaluated the predictive accuracy of 1-, 2-, 3-, and 5-years OS using ROC curves. The corresponding AUCs and 95% CIs for CTC score, LSC17 score, and APS were computed as shown in Figure 4. The differences in AUCs of CTC score versus LSC17 score and CTC score versus APS at four time points were not statistically significant (Supplementary Table S11), suggesting that CTC score can achieve a similar predictive accuracy compared with LSC17 score and APS. Additionally, we simultaneously included CTC score, LSC17 score, and APS into the multivariable Cox regression (Figure 5). CTC score (HR = 3.65; 95% CI, 1.37 to 9.7; p = 0.0095) and APS (HR = 1.84; 95% CI, 1.06 to 3.18; p = 0.0297) remained statistically significant, suggesting that both CTC score and APS could capture additional prognostic information compared with LSC17 score. Furthermore, the additional prognostic information captured by the CTC score was different from that captured by APS.