Identification of single-cell blasts in pediatric acute myeloid leukemia using an autoencoder

Study participants’ details

Cryopreserved bone marrow specimens from 20 pediatric patients with de novo AML were obtained from the Children’s Oncology Group (COG) Biopathology Center. All patients were enrolled on the phase III treatment protocol AAML1031 and were randomized to the standard arm (Arm A). The treatment regimens and outcomes for this trial have been reported (Aplenc et al, 2020). Sample selection criteria included the following: diagnosis/relapse pairs from patients enrolled on AAML1031, Arm A, who consented to banking and who had three or more vials of bone marrow (not peripheral blood) banked at both timepoints. The requirement for multiple vials is imposed by the NCI to avoid depletion of the holdings for any individual patient and timepoint. In addition, the vials should contain at least a certain number of cells.

All patients contributed a sample from diagnosis and relapse. Four of the patients also contributed a remission sample at the end of induction 1 chemotherapy (EOI1). Clinical flow cytometry reported residual blasts in only one of these samples (0.02%). One additional remission sample without residual AML at EOI1 was included from a non-COG patient who was not enrolled in the AAML1031 study but followed the same regimen. The National Cancer Institute’s Central Institutional Review Board (IRB) and IRB at each enrolling center approved the study, and patient families gave written informed consent, in accordance with the Declaration of Helsinki, to allow the bone marrow samples to be banked for future research. Bone marrow samples were enriched for mononuclear cells by density centrifugation and cryopreserved. Patient demographic and clinical information was provided by the COG Statistics and Data Center (Tables S1, S2, and S3). The cohort shows minor differences in sex, age, and mutational frequency compared with global pediatric AML frequencies. These differences are due to sample selection based on the above-mentioned criteria. The incidence in males and females is equal in global pediatric AML statistics (Hossain & Xie, 2015), whereas the cohort contains 62% males (13) and 38% females (8). In addition, we have an overrepresentation of KMT2A rearrangements (Conneely & Rau, 2020).

Spectral flow cytometry

Frozen bone marrow samples were thawed by centrifugation using cryo–thaw devices (Medax) and collected in cell culture medium (RPMI 1640, 10% FBS [Gibco], 1 × GlutaMAX [Gibco], and 1 × penicillin–streptomycin [Gibco]) supplemented with 5 U/ml Benzonase. Cells were spun down and resuspended in cell culture medium supplemented with 2 U/ml Benzonase.

For the staining with fluorophore-coupled antibodies, 8 × 10⁵ cells were transferred into a 96-well V-bottom plate, spun down, resuspended in 50 μl of Zombie NIR Fixable Viability Dye (1:500; BioLegend) for 15 min, and washed with PBS with 0.5% BSA. 50 μl of FcX Fc receptor block (1:200; BioLegend) was added for 15 min, followed by another washing step. Subsequently, cells were incubated with 50 μl of a flow cytometry antibody mixture for 15 min. After washing, samples were fixed using 50 μl of 2% PFA for 15 min, washed with PBS, and analyzed on Cytek Aurora 5L Spectral Analyzer. All incubation and washing steps during the staining procedure were performed at 4°C.

The initial flow cytometry panel consisted of 26 markers: CD45, CD19, CD3, CD4, CD8, CD56, CD34, CD45RA, HLA-DR, CD117, CD38, CD33, CD35, CD13, CD123, CD64, CD11b, CD11c, CD14, CD203c, CD71, CD42b, CD16, CD300e, CD112, and CD155. We ignored CD42b for the analysis because the signal for CD42b was due to platelet aggregates on cells and did not detect megakaryocytes specifically. We excluded CD16, CD300e, CD112, and CD155 entirely because of staining artifacts. We selected 21 markers from the flow cytometry panel and included CD45, CD19, CD3, CD4, CD8, and CD56 to distinguish the lymphoid and myeloid lineages. For a more fine-grained insight into the myeloid lineage, we included CD34, CD45RA, HLA-DR, CD117, CD38, CD33, CD35, CD13, CD123, CD64, CD11b, CD11c, and CD14 for granulocytes, monocytes, and dendritic cell development. We included CD203c for basophil identification, and CD71 for erythroblast identification. We also measured the forward and side scatter of the cells, which provided information about the cells' size and internal complexity. Finally, we used this 23-dimensional dataset for downstream analysis (antibodies are listed in Table S4). We excluded two samples from the analysis because of low cell counts (<10,000 cells).

Data processing

FCS files were loaded into FlowJo (FlowJo Software, 2023) (Tree Star), and compensation of the flow cytometry data was manually adjusted. Two low-quality samples (Section 5.2) were excluded because of low cell counts (FlowJo Software, 2023). After exclusion of cell doublets and dead cells (for this gating strategy, see Fig S1A), 10,000 random cells per file were exported. The data were transformed with an inverse hyperbolic sine function (cofactor range 500–9,030, Table S5), typically used to manage the wide dynamic range of flow cytometry data. These cofactors were manually determined for each marker, to achieve a distribution with two clear peaks with the negative peak being around zero (den Braanker et al, 2021). Cofactors remained constant for all data subsets. Second, each marker was scaled between 0 and 1 using MinMaxScaling (Section 5.12). The data before these preprocessing steps are considered raw data. After the preprocessing, the data were clustered using FlowSOM (Van Gassen et al, 2015) into 80 clusters. Based on the expression of specific markers in the clusters, we merged some of them into the main lymphocyte and myeloid lineages and subsequently excluded B cells, T cells, NK cells, erythroblasts, plasmacytoid dendritic cells, and basophils. Neutrophils were not observed as the samples were enriched for mononuclear cells before cryopreservation. All remaining clusters expressed either CD45, CD34, CD33, CD64, HLA-DR, or CD71 and comprised either healthy HSPCs, healthy myeloid cells, or malignant AML blasts. These 350’250 cells are collectively referred to as HSPCs/myeloid cells from this point forward and used for downstream analyses.

Individuality

In our preliminary analyses, we noticed that patients and timepoints had quite distinct marker expression profiles. To quantify the degree of similarity across samples, patients, or timepoints, we calculated an individuality score (Wagner et al, 2019). This metric evaluates the class distribution of a cell’s nearest neighbors. If most of the neighbors fall in the same class as the cell investigated, then the cell resembles other cells of its own class. In this case, the class is self-contained. Conversely, if the neighboring classes are evenly represented, there is class mixing and the cell exhibits phenotypic similarities with multiple classes, reducing its class-specific uniqueness. We used the kNN (Cover & Hart, 1967) method from sklearn in combination with posterior probability analyses to determine how self-contained the classes are. The posterior probability that a cell belongs to a class i given its k-nearest neighbors is given byp(class_i|k_NN)=n_NNc=in_class_i∑i=1C1n_class_i.

In the above equation, k_NN are the k-nearest neighbors of a cell, and n_NN_c=i are the number of nearest neighbors belonging to class i. n_class_i is the number of cells that belong to class i in the data. C is the total number of classes. This probability was calculated for each class and each cell.

Single-cell annotation

We split the HSPC/myeloid cell raw data into subsets: remission cells and per-patient non-remission cells (comprising diagnosis and relapse cells). The MinMaxScaling was applied to each subset independently. Each subset was clustered with PhenoGraph (Levine et al, 2015) using different numbers of nearest neighbors (k parameter). For each k value, the clustering was repeated 10 times with different random initializations. The adjusted rand index (ARI [Rand, 1971]) between all pairs of clustering was calculated, and the clustering with the highest average ARI was selected. 19 clusters for the remission samples were obtained and annotated based on established markers for hematopoietic differentiation such as CD33, CD34, CD117, CD64, CD11b, CD11c, and CD14 (Porwit & Béné, 2018). This allowed us to merge the 19 clusters into four subsets, namely, early HSPC, monoblasts/myeloblast, promonocyte, or monocyte. The low expression signal for CD19 expression on monocytes was due to background staining of the antibody and not considered aberrant. Three clusters were considered as cell debris because of expression profiles, which are not present in healthy bone marrow. For the annotation of non-remission HSPC/myeloid cells, all patients were clustered separately to achieve a high resolution for blast annotation. Clusters with aberrant expression profiles of normal myeloid differentiation with at least one of the following characteristics were considered blasts: (a) the aberrant expression of the NK cell marker CD56, (b) an abnormal low expression of a canonical myeloid (CD33) or leukocyte marker (CD45), and (c) the overexpression or asynchronous marker expression for myeloid cells, for example, a cluster with the overexpression of a precursor marker CD117 with an absence of CD34, but the overexpression of CD4 (Porwit & Béné, 2018).

Autoencoder training

The raw HSPC/myeloid remission cells were split into a test and training set. The training set was further split into five cross-validation (CV) folds, yielding five training and validation sets. For each trained model, the training data were preprocessed as described in “Data processing.” The test and validation data were transformed with the same cofactors and min-max–scaled with the minimum and maximum of the training data. The test data were not used for the autoencoder training but for the training of the classifier (Section 5.7).

Autoencoders were implemented in Python using the PyTorch (Paszke et al, 2019) and PyTorch Lightning (Falcon et al, 2020) libraries. We tested both autoencoder and variational autoencoder architectures with different hyperparameter choices (Table S6). We randomly selected 100 hyperparameter configurations and trained both model architectures on each of the five cross-validation folds, resulting in a total of 1,000 tested models. All models were trained to reconstruct the 23 flow cytometry inputs (Section 5.2) using a mean-squared error (MSE) loss. The variational autoencoder used an additional Kullback–Leibler loss (Kingma & Welling, 2013 Preprint). The model with the lowest average MSE over all 5 CV folds was retrained on the full training set and used for further analyses.

Classification of malignant cells

We trained binary classifiers using the autoencoder’s reconstruction error for each marker; that is, the input of the classifier was a 23-dimensional input vector containing each marker’s error. As described in Section 5.6, the remission test data (i.e., data not used to train the autoencoder) were used as training data for the classifier. These data were combined with the manually annotated malignant cells from non-remission samples to create leave-one-out cross-validation splits. In each fold, annotated malignant cells (blasts) of one patient were set aside for validation. The size of the excluded malignant cluster was subsampled to match the size of the validation fold of the remission cells (1/20 of the remission test set). The remaining malignant cells served as the training set for that fold and were also subsampled to match the size of the remission training set (19/20 of the remission test set). This resulted in 20 balanced training and validation splits. For each split, the data were preprocessed with the same scaling factors used in the best-performing autoencoder (Section 5.3).

We trained 100 models using these cross-validation splits and selected the model with the highest average validation accuracy over all folds. The selected model was then retrained on all folds, maintaining the balance in the data. The models and parameters explored are listed in Table S7. All models were implemented with the Scikit-learn (sklearn; Pedregosa et al, 2011) library in Python. For any parameter not listed, we used the default sklearn values.

Cell-type prediction

We used the sklearn (Pedregosa et al, 2011) k-nearest neighbor (Cover & Hart, 1967) classifier to predict cell types of the non-remission samples based on their latent space encoding. The kNN classifier was trained using the latent encoding of the remission cells and used to classify the non-remission cells. We used the default sklearn parameters (n_nearest_neighbors = 5) and only considered the remission cells as potential neighbors in the classification task.

Surviving cell analysis

We used a Gaussian mixture model (GMM) to investigate whether diagnosis cells that survive therapy might give rise to new relapse phenotypes; namely, using sklearn we fitted two separate GMMs with a single component each to the encoded diagnosis and relapse cells (i.e., diagnosis-GMM and relapse-GMM, respectively) (Pedregosa et al, 2011). For each cell, we evaluated the log-likelihood under both GMMs. Cells present at diagnosis that had higher log-likelihood under the relapse-GMM were labeled “Diagnosis like relapse.” Similarly, relapse cells with higher log-likelihood under the diagnosis-GMM were labeled “Relapse-like diagnosis.”

Changes in the fraction of positive cells

To investigate the changes in marker expression between diagnosis and relapse blasts, we measured the fraction of cells positively stained for each marker. For the cells predicted by our classifier to be blasts (Section 5.7), we used marker-specific thresholds to determine which cells were positive for each marker. These thresholds were based on naturally occurring internal control cells with negative expression in the remission cells (Table S5). We calculated the fraction of predicted blasts that were positive for each marker and compared this fraction between diagnosis and relapse samples of the same patient (18 patients).

Additional models trained directly on expression

We compared our computational workflow with models trained directly on the marker intensities. We trained these models following the same cross-validation scheme and hyperparameter search as described in Section 5.7 with three major differences. (1) The input was the marker expression/intensities instead of the reconstruction error. (2) Because there was no autoencoder used, we did not use the scaler used in the autoencoder training. Instead, we used the minimum and maximum of each train fold of the cross-validation to scale both the training and validation sets of that cross-validation fold with min-max scaling. (3) We trained extra multilayer perceptron (MLP) models, because we want to compare the latent space of the MLP with that of the autoencoder. However, the initial architectures of the MLPs had hidden layer dimensions ranging from 10 to 50 and the autoencoder has a latent dimension of 4. Therefore, we trained 50 extra models with parameter “hidden_layer_sizes” = [10, 4], resulting in a 4-dimensional layer comparably to the autoencoder. We randomly sampled the hyperparameters “activation function” and “solver.”

We evaluated the accuracy of these models in single-cell blast classification on the remission cells and annotated blasts (Section 5.5, Fig S3A and B). For the comparison of the latent space, we investigated the latent space of the MLP with the highest accuracy overall and the highest accuracy of the MLPs with the “hidden_layer_sizes” = [10, 4]. For the overall best model, the hidden layers were as follows: [40, 60, 40, 10] and we used the UMAP to map the activation of each of these layers to a 2D space. We observed complete separation between the remission and malignant cells. The remission cells did form a developmental trajectory, but the malignant cells formed separate clusters outside of this trajectory. For the smaller MLP, we looked at the activations of the last, four-dimensional layer. Fig S4A shows the first two dimensions. We also investigated the third and fourth dimensions, as well as an UMAP of all four dimensions. In all cases, we observed a separation of the remission and malignant cells.

MinMaxScaling

Min-max scaling scales all values between 0 and 1. The equation is as follows:Xmin_max_scaled=X-minXmaxX-minX

We performed this scaling for each marker separately. Thus, to scale for n cells the marker m, min(X) is the minimum value among the n cells for marker m. Max(X) is the maximum value of marker m among the n cells. X is the value of marker m in a given cell, and this formula is applied to each cell and each marker.