Automated bone marrow cytology using deep learning to generate a histogram of cell types.

Background: Bone marrow cytology is required to make a hematological diagnosis, influencing critical clinical decision points in hematology. However, bone marrow cytology is tedious, limited to experienced reference centers and associated with inter-observer variability. This may lead to a delayed or incorrect diagnosis, leaving an unmet need for innovative supporting technologies.
Methods: We develop an end-to-end deep learning-based system for automated bone marrow cytology. Starting with a bone marrow aspirate digital whole slide image, our system rapidly and automatically detects suitable regions for cytology, and subsequently identifies and classifies all bone marrow cells in each region. This collective cytomorphological information is captured in a representation called Histogram of Cell Types (HCT) quantifying bone marrow cell class probability distribution and acting as a cytological patient fingerprint.
Results: Our system achieves high accuracy in region detection (0.97 accuracy and 0.99 ROC AUC), and cell detection and cell classification (0.75 mean average precision, 0.78 average F1-score, Log-average miss rate of 0.31). Conclusions: HCT has potential to eventually support more efficient and accurate diagnosis in hematology, supporting AI-enabled computational pathology.

Introduction

A bone marrow study is the foundation of making a hematological diagnosis, with an estimated 700 000 bone marrow studies performed annually in the US[1]. It is performed to investigate a clinically suspected hematological disorder, as part of lymphoma staging protocols and to assess bone marrow response to chemotherapy in acute leukemias[2]. Information is extracted by a hematopathologist from the multiple components that comprise a bone marrow study and then integrated with clinical information to make a final diagnostic interpretation[2]. Much of this interpretation relies on visual features of bone marrow cells and tissue viewed through a light microscope[2] or more recently, via high-resolution scanned digital whole slide images (WSIs) of pathology specimens, known as digital pathology[3,4]. One component of a bone marrow study, called the aspirate, consists of particles of bone marrow tissue that are smeared onto a glass slide to allow individual bone marrow cells to be analyzed for subtle and complex cellular features that represent the morphological semantics of the tissue, known as cytology[2,5]. As per international standards, aspirate cytology includes a nucleated differential cell count (NDC), where 300-500 individual bone marrow cells are manually identified, counted, and classified into one of many discrete categories by a highly experienced operator such as a hematopathologist[2]. Bone marrow cytology and the NDC are required for many critical clinical decision points in hematology. For example, the identification of leukemic blasts may lead to immediate initiation of flow cytometry, karyotype, and induction chemotherapy in acute myeloid leukemia (AML)[6,7]. Similarly, the identification of subtle cytological changes in bone marrow cells is necessary for the diagnosis and risk stratification in patients with a myelodysplastic syndrome (MDS)[8]. Failure to recognize and quantify abnormal cell populations in the aspirate in a timely and accurate manner may lead to delayed or incorrect diagnosis. In the context of a busy reference hematopathology lab, performing cytological review on every bone marrow aspirate specimen is tedious and subject to inter-observer variability[9-11]. At the same time, smaller community centers often lack sufficient technical expertise to correctly interpret bone marrow aspirate cytology[12]. One study estimated that up to 12% of MDS cases are misdiagnosed due to the inability to recognize morphological dysplasia in aspirate specimens in less experienced centers[11]. This leaves an unmet clinical need for innovative computational pathology tools that will support the aspirate review process.

Artificial Intelligence (AI) describes the aspiration to build machines, or computer software, with human-like intelligence[13,14]. One particular type of AI algorithm, called deep learning, has shown considerable success in digital image analysis and image classification tasks in many domains[15,16]. In the pathology domain, deep learning represents a computational pathology tool that has been successfully implemented in many non-hematopoietic pathology sub-specialties using WSIs of solid tissue pathology specimens, known as histopathology[17]. Numerous studies have demonstrated the ability of deep networks to perform tasks such as binary morphological classification, distinguishing tumor from normal tissue[18-23], as well as histomorphological tissue grading[24]. While these approaches generally deliver excellent classification results, they do not capture the nuances or complexity inherent in bone marrow aspirate cytology. Specifically, the vast majority of morphological analysis in the hematopoietic system is performed at the level of cellular resolution and represents non-binary classification based on subtle morphological features such as dysplasia in MDS. The application of deep learning to diagnostic hematopathology will therefore require unique solutions that are tailored to these distinct cytomorphological challenges.

While there are several commercial computational pathology workflow support tools developed for analysis of peripheral blood cytology[25], there are currently no clinical-grade solutions available for bone marrow cytology. In comparison to blood film cytology, bone marrow aspirates are complex cytological specimens. Aspirates contain only a small number of regions suitable for cytology, significant non-cellular debris and many different cell types that are often aggregated or overlapping[2,5]. This has rendered bone cytology as a relatively challenging computational pathology problem. Aspirate cytology can be roughly modeled into three distinct computational steps to reflect real-world hematopathology practice. The first problem is region of interest (ROI) detection, where a small number of regions or tiles suitable for cytology must be selected from large WSI prior to cell detection and classification. ROI selection has previously been accomplished in bone marrow aspirates by a human operator manually selecting and cropping the appropriate tiles in aspirate WSIs[26,27]. Second, there is the problem of object detection, where individual bone marrow cells or non-cellular objects must be identified in aspirate WSI as both distinct and separate from background. Prior approaches have employed deep learning for object detection such as regional CNN (R-CNN), fast and Faster R-CNN[28,29]. These approaches utilize region proposals for object detection followed by a separate method such as object classification, which renders them complex to train and hence computationally inefficient[26,30,31]. Third and finally there is the problem of object classification, where individual bone marrow cells or non-cellular objects must be assigned to one of numerous discrete classes based on nuanced and complex cytological features. This complexity increases in MDS, where morphological dysplasia creates subtle cytological changes.

One study attempted to address the second and third problems using fine-tuning of Faster R-CNN and the VGG16 convolutional network[26]. However, this approach proved operationally slow and is not likely scalable to a clinical diagnostic workflow. Therefore, novel, efficient and scalable computational pathology approaches are needed to support bone marrow aspirate cytology; specifically approaches that add full end-to-end automation, i.e., from unprocessed WSI to bone marrow cell counts and classification.

Recently, a deep learning model called You Only Look Once (YOLO) was developed for real-time object detection to specifically address the detection and classification problems in complex image analysis domains[31]. YOLO uniquely allows for object detection and classification to occur in a single step, where all objects in an image are simultaneously identified and localized by a “bounding box” and then assigned a class probability by the same deep network[31]. The YOLO model outputs a set of real numbers that captures both object localization in an image and an object class probability, therefore solving both object detection and classification problems simultaneously in a regression approach[31]. In addition, the most recent version of YOLO, YOLOV4, has been optimized for small object detection and uses complete intersection-over-union loss (CIoU), which results in faster convergence and better accuracy for bounding box prediction[32]. These factors collectively lead to increased computational efficiency and speed compared to previous methods[28-32]. YOLO can perform object detection and classification on multiple image objects which are complex and overlapping in virtual real-time (milliseconds)[31], and consequently has been applied in several real-world problems including autonomous driving[31,33-35]. Recently, YOLO has been applied in some medical domain problem such as pathology. For example in ref. [36], YOLO has been applied to assess the cell types in bone marrow smears. However, only 7 cell types have been considered in that study. Moreover, the tiles need to be selected manually by the user.

In this work, we demonstrate the first automated end-to-end AI architecture for bone marrow aspirate cytology. We first employ and implement a fine-tuned DenseNet model to rapidly and automatically select appropriate ROI tiles from a WSI for bone marrow aspirate cytology. Subsequently, we implement a YOLO model trained from scratch to detect and assign class probabilities to all cellular and non-cellular objects in bone marrow aspirate digital WSI. Collective cytological information for each patient is then summarized as a Histogram of Cell Types (HCT), which is a novel information summary quantifying the class probability distribution of bone marrow cell types, acting as a cytological fingerprint. A histogram is generally a representation of a distribution, a very old graphical technique to count discrete values[37]. Our approach shows cross-validation accuracy of 0.97 and precision of 0.90 in ROI detection (selecting appropriate tiles), and mAP (mean Average Precision) of 0.75 and average F1-score of 0.78 for detecting and classifying 16 key cellular and non-cellular objects in aspirate WSIs. Our approach has potential to fundamentally change the process of bone marrow aspirate cytology, leading to more efficient, more consistent and automated diagnostic workflows, and providing a foundation for computational pathology driven augmented diagnostics and precision medicine in hematology.

Methods

This work proposes a new end-to-end AI architecture for bone marrow aspirate NDC based on machine learning and deep learning algorithms (Fig. 1).

Fig. 1

End-to-end AI architecture for bone marrow aspirate cytology.

In this architecture, initially, our Region of Interest (ROI) detection model is run on unprocessed bone marrow aspirate WSI. A grid is created on an original Whole-Slide Image (WSI) and ROI tiles are selected using ROI detection model. Subsequently, a You-Only-Look-Once (YOLO)-based object detection and classification is run to localize and classify cells in the selected tiles and generate the Integrated Histogram of Cell Types (IHCT).

Dataset

This study was approved by the Hamilton Integrated Research Ethics Board (HiREB), study protocol 7766-C. As this study protocol was retrospective, it was approved with waiver of patient consent. Digital whole slide images (WSI) were acquired retrospectively and de-identified and annotated with only a diagnosis, spanning a period of 1-year and 1247 patients. This starting dataset represented the complete breadth of diagnoses over this period in a major hematology reference center. WSI were then sampled from this dataset for model development and validation as described in Table 1 and Supplementary Table S3. These images were scanned with either an Aperio Scanscope AT Turbo or a Huron TissueScope at 40X and acquired as SVS and tif file format.

Table 1

Diagnostic tags and the number of patient WSI used for training and test-validation in each category for the ROI detection model.

Diagnostic tags	Used in training	Used in test-validation	Number of patients
Normal	80	18	98
Myelodysplastic syndrome (MDS)	15	3	18
Acute leukemia	23	5	28
Lymphoproliferative disorder	28	7	35
Plasma cell neoplasm	19	4	23
Hypercellular	5	1	6
Erythroid hyperplasia	3	0	3
Myeloproliferative neoplasm (MPN)	4	1	5
Inadequate	11	3	14
Hypocellular	6	2	8
MPN/MDS	2	0	2
MPN	3	1	4
Necrosis	2	1	3
Carcinoma	3	0	3
Total	204	46	250

Data annotation and augmentation strategy

ROI tiles and individual bone marrow cell types included were annotated by expert hematopathologists as the ground truth or reference standard in WSI images used for model training and test-validation as described below. This follows ICSH guidelines, where expert pathologists are considered the reference standard for bone marrow aspirate cytology in clinical diagnosis[2]. Data augmentation was applied to increase the diversity of the input image types. Generally, there are two categories for pixel-wise adjustments augmentation, photometric distortion, which includes hue, contrast, brightness, saturation adjustment, and adding noise; and geometric distortion, which includes flipping, rotating, cropping, and scaling. As we had an imbalanced class distribution within our dataset for the ROI detection model (70,250 inappropriate and 4,750 appropriate tiles), it was necessary to apply one of the over-sampling or under-sampling methods to prevent misclassification. To address this, a number of the above augmentation techniques were applied to the training data during the learning process to over-sample the appropriate ROI tiles and train the model correctly. Subsequently, after applying augmentation, the dataset in this phase contained 98,750 annotated images for training, including 70,250 inappropriate ROI tiles and 28,500 appropriate ROI tiles (Supplementary Table S1). For the cell detection and classification model, after annotating the objects inside ROI tiles by using LabelImg tool (Supplementary Fig. S1)[38], in addition to the above augmentation categories, other techniques were also applied, like cutmix[39] which mixes 2 input images, and mosaic, which mixes 4 different training images. Accordingly, after applying augmentation, the dataset in this phase contained 1,178,408 annotated cells for training, including 119,416 neutrophils, 44,748 metamyelocytes, 52,756 myelocytes, 17,996 promyelocytes, 173,800 blasts, 117,392 erythroblasts, 1,012 megakaryocyte nuclei, 57,420 lymphocytes, 25,036 monocytes, 7,744 plasma cells, 10,956 eosinophils, 308 basophils, 4,664 megakaryocytes, 246,532 debris, 8,404 histiocytes, 1,452 mast cells, 174,724 platelets, 25,740 platelet clumps, and 88,308 Other cell types (Supplementary Table S2). To enhance generalization, the augmentation was only applied on the training set in each fold of the cross-validation.

Region of interest (ROI) detection method

The first phase in the proposed architecture is ROI detection. The ROI detection was applied to extract tiles from a WSI and examine if that tile was suitable for diagnostic cytology. To accomplish this, a deep neural network was built, fine-tuned and evaluated on aspirate digital WSI tiles. In the ROI detection method, initially 98,750 tiles (including augmented data) in 512 × 512-pixel size in high resolution are extracted and acquired from 250 WSI. To choose the tiles, a grid of 15 rows and 20 columns was created on each digital WSI and tiles were selected from the center of each grid cell, ensuring all tiles have been sampled from the WSI evenly. Appropriate and inappropriate ROI tiles were annotated by an expert hematopathologist; appropriate ROI tiles needed to be well spread, thin and free of red cell agglutination, overstaining and debris; and contain at least one segmentable cell or non-cellular object as outlined above. Then, a deep neural network based on DenseNet121 architecture[39] was fine-tuned to extract features from each tile. A binary classifier was added in the last layer of the model to classify appropriate and inappropriate tiles. This network was trained using a cross entropy loss function and AdamW optimizer with learning rate 1e-4 and weight decay 5.0e-4. Also, a pretrained DenseNet121 was applied to initialize all weights in the network prior to fine-tuning. The entire network was fine-tuned for 20 epochs with 32 batch size.

We applied patient-level 5-folds cross-validation to train and test the model. Hence, the dataset (98,750 tiles) was split into two main partitions in each fold, training and test-validation, 80% (204 WSIs including 80,250 tiles) and 20% (46 WSIs including 18,500 tiles), respectively. The test-validation was also been split into two main partitions, 70% validation and 30% test. To ensure that enough data for each class was chosen in our dataset, the above split ratios were enforced on appropriate and inappropriate tiles separately. The dataset was split into training, validation and test sets at patient level, such that each set has a patient WSI that does not come in the other sets to prevent data leakage. In each fold, the best model was picked by running on the validation partition after the training and then evaluated on unseen patients in the test dataset. Extracting ROI tiles for further processing for the cell detection and classification model was the primary aim of using the ROI detection model. To this end, the ROI detection model should be able to minimize false positives in the result. Therefore, the precision has been considered as a key performance metric to select the best model.

Cell detection and classification

The next phase was cell detection and classification applied on ROI tiles of 512 × 512 pixels in high resolution. To accomplish this, the YOLOv4 model was customized, trained and evaluated to predict bounding boxes of bone marrow cellular objects (white blood cells) and non-cellular objects inside the input ROI tile and classify them into 19 different classes. In this architecture, CSPDarknet53[40] was used as the backbone of the network to extract features, SPP[41] and PAN[42] were used as the neck of the network to enhance feature expressiveness and robustness, and YOLOv3[43] as the head. As bag of specials (BOS) for the backbone, Mish activation function[44], cross-stage partial connection (CSP) and multi input weighted residual connection (MiWRC) were used. For the detector, Mish activation function, SPP-block, SAM-block, PAN path-aggregation block, and DIoU-NMS[45] were used. As bag of freebies (BoF) for the backbone, CutMix and Mosaic data augmentations, DropBlock regularization[46], and class label smoothing were used. For the detector, complete IoU loss (CIoU-loss)[45], cross mini-Batch Normalization (CmBN), DropBlock regularization, Mosaic data augmentation, self-adversarial training, eliminate grid sensitivity, using multiple anchors for single ground truth, Cosine annealing scheduler[47], optimal hyperparameters and random training shapes were used. In addition, the hyperparameters for bone marrow cell detection and classification were used as follows: max-batches is 130,000; the training steps are 104,000 and 117,000; batch size 64 with subdivision 16; the polynomial decay learning rate scheduling strategy is applied with an initial learning rate of 0.001; the momentum and weight decay are set as 0.949 and 0.0005 respectively; warmup step is 1,000; YOLO network size set to 512 in both height and width; anchor size set to 13, 14, 19, 18, 29, 30, 19, 64, 62, 20, 41, 39, 35, 59, 50, 49, 74, 35, 56, 62, 68, 53, 46, 87, 70, 70, 95, 65, 79, 85, 101, 95, 87,129, 139,121, 216, 223.

Similar to the ROI detection method above, patient-level 5-folds cross-validation was applied to train the model here. Therefore, each fold is divided into training and test-validation partitions, 80% and 20% respectively. The test-validation data portion was split into two main partitions (70% validation and 30% test). Additionally, to ensure that enough data for each class was chosen in our dataset, the mentioned portions were enforced on each object class type individually. In each fold, the best model was picked by running it on the validation partition and then evaluation on the test (unseen) dataset was performed using the mean average precision (mAP).

After training and applying the cell detection and classification model on each tile, the Chi-square distance (Eq. (1)) was applied to determine when the IHCT converges.If IHCT converged, the bone marrow NDC is completed and represented by the IHCT, otherwise, another tile is extracted, and the previous process applied again iteratively until it converges.

To calculate Chi-square distance, the number of following cellular objects, as well as BMME ratio (Eq. (2)), were utilized: “neutrophil”, “metamyelocyte”, “myelocyte”, “promyelocyte”, “blast”, “erythroblast”, “lymphocyte”, “monocyte”, “plasma cell”, “eosinophil”, “basophil”, “megakaryocyte”.

Cell types were chosen to include all bone marrow cell types traditionally included in the NDC, as well as several additional cell or object types that have diagnostic relevance in hematology (“megakaryocytes”, “megakaryocyte nuclei”, “platelets”, “platelet clumps” and “histiocytes”).

Evaluation

To evaluate the ROI detection model in predicting appropriate and inappropriate tiles, we calculated common performance measures such as accuracy, precision (PPV-positive predictive value), recall (sensitivity), specificity, and NPV (negative predictive value), as shown by the following equations:

For the ROI detection model, Tp, Tn, Fp and Fn in the above equations are defined as:

Tp (True Positive): The number of appropriate tiles which predicted correctly

Tn (True Negative): The number of inappropriate tiles which predicted correctly

Fp (False Positive): The number of inappropriate tiles which predicted as appropriate tiles

Fn (False Negative): The number of appropriate tiles which predicted as inappropriate tiles

To assess the performance of the proposed cell detection and classification method, Average Precision (AP) was used with 11-point interpolation (Eq. (8)). Also at the end, the mean Average Precision (mAP) [48] was calculated for all the AP values (Eq. (10)). The value of recall was divided from 0 to 1.0 points and the average of maximum precision value was calculated for these 11 values. It is worth mentioning that the value of 0.5 was considered for Intersection over Union (IoU) in AP for each object detection and >0.75 has been used for class probability. In addition, Precision, Recall, F1-score (Eq. (12)), average IoU (Eq. (11)) and log-average miss rate (Eq. (13)) have been calculated here for each object type.where

For the cell detection and classification model, Tp, Fp and Fn in Eq. (4) and Eq. (5) are defined as:

Tp (True Positive): The number of all cellular and non-cellular objects which predicted correctly.

Fp (False Positive) and Fn (False Negative): The number of all cellular and non-cellular objects which not predicted correctly

The Log-average miss rate[49] is calculated by averaging miss rates at 9 evenly spaced FPPI points between 10−2 and 100 in log-space.Where a1, a2,a9 are positive values corresponding the miss rates at 9 evenly spaced FPPI points in log-space, between 10−2 and 100.

Table 2

Evaluation of the ROI detection using 5-fold cross-validation to calculate accuracy, precision (PPV-Positive Predictive Value), recall (Sensitivity), specificity, and NPV (Negative Predictive Value).

Metrics	%
Average Cross-validation Accuracy	0.97
Average Cross-validation Precision (PPV)	0.90
Average Cross-validation Specificity	0.99
Average Cross-validation Recall (Sensitivity)	0.78
Average Cross-validation NPV	0.99

All these metrics were computed in each test (unseen) fold separately and then the average was calculated.

Table 3

Performance result of the proposed cell detection and classification model.

Object class	Precision	Recall	F1 score	Log-average miss rate	AP@0.5
Neutrophil	0.84	0.91	0.87	0.21	0.90
Metamyelocyte	0.68	0.79	0.73	0.37	0.77
Myelocyte	0.80	0.82	0.81	0.34	0.80
Promyelocyte	0.60	0.67	0.64	0.53	0.62
Blast	0.87	0.90	0.88	0.34	0.84
Erythroblast	0.86	0.92	0.89	0.17	0.92
Megakaryocyte nucleus	0.80	0.57	0.67	0.18	0.60
Lymphocyte	0.73	0.65	0.69	0.49	0.66
Monocyte	0.84	0.71	0.77	0.36	0.72
Plasma cell	0.75	0.69	0.72	0.33	0.72
Eosinophil	0.93	0.94	0.93	0.06	0.97
Megakaryocyte	1.00	0.79	0.88	0.19	0.82
Debris	0.85	0.80	0.82	0.34	0.79
Histiocyte	0.90	0.53	0.67	0.5	0.54
Platelet	0.84	0.64	0.73	0.33	0.64
Platelet clump	0.93	0.61	0.73	0.41	0.62
Average	0.83	0.75	0.78	0.32	mAP@0.5 =0.75

Table 4

Performance result of using active learning.

	Iteration 1		Iteration 2		Iteration 3		Iteration 4		Iteration 5		Iteration 6		Iteration 7		Iteration 8
Object class	Count	AP	Count	AP	Count	AP	Count	AP	Count	AP	Count	AP	Count	AP	Count	AP
Neutrophil	680	0.75	1256	0.82	1568	0.83	1756	0.85	1895	0.86	2050	0.89	2398	0.91	2714	0.90
Metamyelocyte	480	0.60	605	0.66	752	0.69	785	0.72	856	0.76	925	0.75	986	0.76	1017	0.77
Myelocyte	390	0.53	589	0.55	665	0.59	720	0.62	869	0.70	950	0.78	1015	0.79	1199	0.80
Promyelocyte	65	0.44	102	0.46	256	0.52	285	0.54	320	0.59	326	0.62	360	0.64	409	0.62
Blast	1050	0.69	1785	0.76	2029	0.78	2590	0.81	2896	0.80	3268	0.83	3526	0.84	3950	0.84
Erythroblast	620	0.72	1150	0.78	1390	0.80	1580	0.82	2028	0.89	2295	0.90	2480	0.92	2668	0.92
Megakaryocyte nucleus	5	0.32	7	0.35	18	0.52	19	0.55	19	0.55	23	0.60	23	0.59	23	0.60
Lymphocyte	390	0.47	530	0.48	689	0.50	706	0.51	780	0.52	1015	0.59	1150	0.62	1305	0.66
Monocyte	62	0.47	98	0.51	295	0.57	368	0.61	423	0.62	485	0.65	520	0.68	569	0.72
Plasma cell	29	0.57	45	0.59	50	0.61	82	0.63	105	0.67	135	0.68	158	0.71	176	0.72
Eosinophil	31	0.59	38	0.63	135	0.83	172	0.86	185	0.88	221	0.95	228	0.95	249	0.97
Megakaryocyte	25	0.49	30	0.52	90	0.77	90	0.77	92	0.78	95	0.80	100	0.81	106	0.82
Debris	1380	0.58	2680	0.62	3450	0.65	3920	0.68	4490	0.73	4901	0.77	5260	0.77	5603	0.79
Histiocyte	38	0.34	72	0.42	147	0.48	163	0.48	168	0.51	174	0.52	182	0.54	191	0.54
Platelet	790	0.41	1680	0.46	2150	0.48	2560	0.52	2890	0.58	3250	0.65	3680	0.65	3971	0.64
Platelet clump	93	0.37	146	0.41	320	0.54	409	0.56	475	0.57	536	0.58	563	0.61	585	0.62
Average	6128	0.52	10813	0.56	14004	0.64	16205	0.66	18491	0.69	20649	0.72	22629	0.74	24735	0.75

Model training started with a small dataset at the first and second iteration, and then is improved (especially on rare cellular objects) in the subsequent iterations by using active learning.