Morphological Estimation of Cellularity on Neo-Adjuvant Treated Breast Cancer Histological Images

This paper describes a methodology that extracts key morphological features from histological breast cancer images in order to automatically assess Tumour Cellularity (TC) in Neo-Adjuvant treatment (NAT) patients. The response to NAT gives information on therapy efficacy and it is measured by the residual cancer burden index, which is composed of two metrics: TC and the assessment of lymph nodes. The data consist of whole slide images (WSIs) of breast tissue stained with Hematoxylin and Eosin (H&E) released in the 2019 SPIE Breast Challenge. The methodology proposed is based on traditional computer vision methods (K-means, watershed segmentation, Otsu’s binarisation, and morphological operations), implementing colour separation, segmentation, and feature extraction. Correlation between morphological features and the residual TC after a NAT treatment was examined. Linear regression and statistical methods were used and twenty-two key morphological parameters from the nuclei, epithelial region, and the full image were extracted. Subsequently, an automated TC assessment that was based on Machine Learning (ML) algorithms was implemented and trained with only selected key parameters. The methodology was validated with the score assigned by two pathologists through the intra-class correlation coefficient (ICC). The selection of key morphological parameters improved the results reported over other ML methodologies and it was very close to deep learning methodologies. These results are encouraging, as a traditionally-trained ML algorithm can be useful when limited training data are available preventing the use of deep learning approaches.


Introduction
Digital pathology has recently become a major player in Cancer research, disease detection, classification, and even in outcome prognosis [1][2][3][4]. Perhaps the most common imaging technique is Hematoxylin and Eosin (H&E), where H stains nuclei blue and E stains cytoplasm pink [5]. Additionally, other immunohistochemistry methods (IHC) that use antibodies to stain antigens or proteins in the tissue are more specific and they can complement H&E [6]. For instance, for TC assessment. A full hand-engineering method for this task was proposed by Peikari [45]. This methodology was based on the extraction of a vector with 125 parameters [45].
Separate to traditional ML techniques, another approach to estimate TC is by Deep Learning (DL) techniques. In the last decade, there has been increased interest in DL methods. One of the main architectures or models is based on the use of neural networks, and a particular case is a Convolutional Neural Network (CNN) [46]. The main advantage of DL techniques is that they do not require the extraction of features manually or by training, as the network learns a series of parameters and weights by itself. However, to achieve this, the network requires a relatively large number of training images with labels, and the number of training data can impact on the performance. In many histopathology applications, labelled training data are still limited as compared to other imaging applications, such as everyday images, like cats and dogs. With the spread use of whole slide digital scanners [47], numerous histopathology images have become available for research purposes. Some have been released to the research community in general in the form of challenges (e.g., https://grand-challenge.org/challenges/) in order to encourage research groups globally to work together, gather annotations, provide training, and testing data sets and benchmark algorithms. The present work follows the 2019 SPIE-AAPM-NCI BreastPathQ: Cancer Cellularity Challenge with the specific objective of "development of quantitative biomarkers for the determination of cancer cellularity from whole slide images (WSI) of breast cancer hematoxylin and eosin (H&E) stained pathological slides" ( https://breastpathq.grand-challenge.org/) and addresses the development of quantitative biomarkers for the challenge.
CNNs have been used for different Histopathology tasks, like segmentation [48,49], detection of a specific image properties [50] , and image grade classification [38,51]. Breast cancer tumour cellularity has also been addressed by deep techniques. For instance, Ziang Pei [52] implemented a direct method based on deep and transfer learning approach with the advantage of avoiding cell segmentation. Akbar [20] presented a traditional hand-engineering approach and a deep neural network.
The methodology in this study is based on a hand-engineering approach similar to the one that was described by [45] based on a 125 parameter vector size, but we selected 22 parameters after a correlation analysis of extracted features with TC. Additionally, the methodology described here is similar to Dong et al. parameters from nuclei; we also include parameters from the neighbourhood of the nuclei and from the full image patch. The study by Fondon [36] also analysed regions around the nuclei; however, our study includes morphology parameters from whole breast ducts region. Thus, the main contributions of this paper in the methodology are the segmentation algorithm that is based on the enhancement of the nuclei region, and an algorithm for breast ducts detection and the parameter derivation and selection and its correlation to tumour cellularity. In the results, the correlation analysis of the morphological parameters with cellularity revealed that stroma concentration has the strongest correlation with TC, which is in agreement with the results that were presented by [16]. Finally, the results were validated with the ICC and compared with similar studies [20,45] indicating an increase in ICC. The methodology described only requires a reduced set of training parameters. Therefore, the methodology described in this paper improves previous results and it may be useful in cases when large training data sets, which are normally required for deep learning approaches, are not available.

Materials
The data set used in this work consists of 2579 patches of tissue stained with H&E, which were extracted from 64 different patients under Neo-adjuvant Treatment. The size of each patch is 512 × 512 pixels. As a reference, a breast tissue is formed of a connective tissue, named stroma, as seen in pink. Lymphocytes can be seen as dark blue round objects and fat zones as white areas.
It also contains ducts,which are responsible for carrying milk, and sometimes arteries, which can be seen as regular clusters of darker blue nuclei grouped into the region (Figure 1). The residual cellularity of each patch was evaluated by 2 pathologists and this assessment was considered to be the ground truth (GT) for this study. The data was released as part of the challenge 2019 and it was collected at the Sunnybrook Health Sciences Centre, Toronto. It comprises a set of whole slide images (WSIs) that have been stained with H&E [45] from 64 patients with residual invasive breast cancer on re-section specimens following NAT therapy. The specimens were handled according to routine clinical protocols and WSIs were scanned at 20× magnification (0.5 µm/pixel).
The images were divided into a training and validation set of images. The cellularity value distribution for the whole data-set is uneven, i.e., most of the patches correspond to cellularity zero and fewer patches were available for higher cellularity. The training images were selected uniformly distributed from cellularity zero to one in order to have an even amount of benign and malignant nuclei. First, 212 images with selected cellularity values from 0 to 1 were processed and 4533 nuclei cells from those selected images were extracted and its corresponding features computed and used as a training set. The remaining 2367 patches were used for validation. Figure 2 shows the selected patches with different levels of cellularity. A third set of test data of 1119 images was available for the challenge contest. These three sets were used in this work.

Methodology
The proposed methodology consists of a sequence of traditional image processing routines, computer vision, and machine learning algorithms that automatically process the full validation set. A large amount of morphological parameters can be extracted and stored orderly in an output table file. A master control routine is responsible for selecting one by one the corresponding image patch to be processed. There are two operational modes. A manualmode useful to train machine learning algorithms. In this mode, nuclei from the selected image with a known classification assignment are fed to the algorithm and output features are saved in an output data file. Subsequently, an automated mode processes the full validation set and gives a TC estimation. This mode is able to process thousands of patch images. Figure 3 presents this process in a diagram. The methodology was implemented in Matlab R (The Mathworks TM , Natick, MA, USA) 2019b version, with functions from the digital image processing, statistical, and machine learning toolboxes. Additionally, QuPath [53], an open source software for digital pathology image analysis, was used to validate the segmentation results obtained by the methodology. Three regions of interest were selected from each image patch, and more than 150 morphological parameters were extracted at inner segmented cell region, neighbourhood around segmented cell, and the full image patch. Under two operational modes images are processed to extract features at nuclei, regional, cluster, and global image regions either to classify and assign Cellularity or to extract same features in order to an archive output.

Segmentation of Nuclei
At the nuclei region, cells from the image patch are segmented by colour separation and binarisation. The image is converted to HSV colour space and, by K-means clustering [54], three main colour images are extracted, say I p for pink, Ib for blue, and Iba for the remaining background component image. A special procedure enhances Ib and weakens background region is calculated as: where Io is the RGB image converted to gray levels and Im is gray image after a median filter of size 3 × 3 was applied. γ is the gamma correction Im of parameter 0.1 represented as γ(Im, 0.1). This value was selected as low to start from a lighter background. K 1 and K 2 are constants that control enhancement of nuclei and weaken background intensity, respectively, and they were experimentally adjusted. First, both of the constants were fixed to 0.5 and, as K 1 is increased and K 2 is reduced, the nuclei region is enhanced. The enhanced image In is binarised by Otsu's algorithm [55] combined with watershed [56] separation of touching nuclei. This segmentation procedure is seen in Figure 4 for three TC cases. A validation analysis of this procedure was done by means of Jaccard Index [57], as shown in (e). The reference images were obtained by a manual segmentation while using the QuPath platform. The results indicate an average Jaccard index of 0.73 after comparing cases of low, medium, and high TC. Although this result might be improved, Peikari and co-authors [45] suggest that cell segmentation accuracy does not have a significant effect in TC assessment.

Extraction of Morphological Parameters
Using the binary image, a set of morphological parameters has been obtained: area, perimeter, roundness, eccentricity, centroid X, centroid Y, major axis, minor axis, and orientation angle. Additionally, a sub region inside the segmented cell body was determined in order to compute the texture and mean HSV values of the nuclei cell body. The parameters were extracted from different window sizes surrounding the segmented cell: 30 × 30, 60 × 60, 90 × 90, and 120 × 120 pixels. The statistical analysis indicated that the best correlation was obtained from the smallest window size. The parameters estimated are the four bins histogram from HSV image and also the regional concentration R of I p, Iba, and Ib from its binary image, determined by the ratio of total white pixels T W by total pixels T P inside window neighbourhood region, R = T W /T P . Figure 5 shows a sample patch image, with its corresponding pink, blue, binary, and concentration image components. Brown areas indicate low concentration and white and yellow areas are high concentration. Breast normal tissue images present ducts and sometimes arteries, which are clearly seen as regular clusters of darker blue nuclei grouped into the region, see Figure 6. These are detected based on set theory: let U be the full image, D be the cluster or duct, D⊆U = ∅. Subsequently, D regions are obtained by subtracting background image Iba from original image U. The following morphological parameters can be extracted from these regions: total cluster area, roundness, number of cells inside cluster, and distance from cells centroid to cluster centroid.
Concentration parameters from the complete full image patch were computed using original RGB image and transformed to HSV. Additionally, a pink colour filter was implemented by means of the Matlab R colour thresholder, which extracts stroma region from the image and computes its mean value. A summary of the full set of features extracted is presented in Table 1; notice that regional parameters are computed at the different window sizes.

Correlation Analysis of Morphological Parameters to TC
The full training set was processed to determine all of the morphological parameters at the three region of interest and those corresponding to Cellularity values equal to (0, 0.1, 0.2, ..., 1.0) were selected, and their mean, standard deviation, maximum, and minimum values were computed. Subsequently, a linear regression and lasso analysis were determined from both analyses. Parameters of coefficient above 0.80 from linear regression and in concordance with lasso selected parameters after redundant removal yield 22 parameters that have the strongest correlation with TC. The plot in in Figure 2 illustrates three parameters with the strongest correlation with TC.
First, the morphological parameters related to segmented nuclei are: eccentricity, roundness, major axis, minor axis, and perimeter. Subsequently, parameters computed at neighbourhood region are: nuclei density concentration, Hue, Value, and Saturation histograms of the regional HSV colour map image. Finally, parameters that are computed from the full image are: HSV average values from HSV components, nuclei concentration, basement concentration, and stroma average concentration determined at output of the pink colour filter. Some of these parameters are graphically displayed in Figure 7.

Training of Machine Learning Algorithms
Machine learning algorithms were trained with the parameter vector of size 22 determined by the statistical correlation analysis. From the training set, 4533 segmented cells were processed and its corresponding extracted parameters were used to generate a prediction function that classifies nuclei cells between malignant and normal cells. Three algorithms were tested: Support Vector Machines, Nearest K-Network, and AdaBoost. The accuracy of training process for every selected method with the training data showed values up to 0.99 due to a high correlation selected features used for training. TPR achieved are 0.97 for SVM, 0.95 for AdaBoost, and 0.97 for KNN (Figure 8). Support Vector Machines (SVM) is a training algorithm for optimal margin classifier [41], and it is based on a determination of a decision function of pattern vectors x of dimension n classifying in either A or B, in our case benign or malignant cells. The input is a set of p examples of x i , i.e., the 22 strongest correlated features extracted. K-Nearest Neighbour method [42] was also selected, because it is one of the most well known algorithms within clustering and data classification, in our case between benign and malignant classes. AdaBoost [43] is a decision tree type learning algorithm that starts from observations of a certain item that is represented by branches and goes to conclusions about item target value or leaves. It has a best performance on binary classification problems.

Assessment of Tumour Cellularity
The methodology to estimate TC is illustrated in Figure 9, with three representative images with increasing TC from top to bottom in each row. In Figure 9a the nuclei are segmented and their corresponding parameters are extracted. Subsequently, the prediction function classifies every cell in either benign (green) or malignant (red), as illustrated in Figure 9b. Next, an estimation of full cell cytoplasm is done by morphological dilation drawn as the white circles around malignant cells (Figure 9c). The full cellularity that is detected region is shown in Figure 9d. TC is computed as the ratio of the area that is covered by cellularity (white in the figure) over the total area (white and black in figure).

Results
An automated estimation of TC was computed from two test data sets. Three prediction functions that were trained by machine learning algorithms were determined to be used with the automated processing software of breast cancer images that classifies cells and computes TC. The method was tested with a training set of 2579 images that were already classified by a pathologist with a TC value. Additionally, it was tested with the 1119 images for submission of SPIE Breast Challenge, with an unknown TC value. Figure 10 shows the statistical behaviour of the method's result for the training set as boxplots.
Dispersion plot indicates the method for approximating to the pathologist classification assignment. The results have a better approximation at higher cellularity values (TC > 0.70) and performs well with KNN algorithm. Additionally, around the middle region (0.4 < TC < 0.6) has a good approximation with AdaBoost. At low cellularity values (TC < 0.3), three methods present deviation, with its higher at cellularity zero, which correspond to images with only benign nuclei cells. According to Minimum Square Error (MSE), SVM performs better overall in the cellularity region.
This result can be validated by a visual inspection of boxplots of Figure 10 in the three cases there is a positive correlation between the actual cellularity (horizontal) and the estimated cellularity (vertical). However, SVM shows less dispersion, especially in the lower values of cellularity as compared with the other two techniques. In order to analyse the limitations of the algorithm, the TC assessment outcomes with the highest errors were analysed. Figure 11 presents three cases where the TC was incorrectly calculated. All three cases in the figure are of benignant tissue with no cellularity; this means that most of the segmented cells should be marked as benignant (green), but several of them are shown in red (b-d), which corresponds to false-positive cases. The expanded area around the cell (e) yields a high TC value instead of the correct value which should have been zero. We assume that this problem is because of the limited amount of cells used to train the algorithm. This problem also explains some of the outliers on the boxplots presented in Figure 10, which are mainly observed at low cellularity. Additionally, this suggests that other classification algorithms should be evaluated.  Statistical analysis of the training set revealed 22 key parameters that have a strong correlation to TC. The Stroma concentration (r = −0.9786), global Value of HSV component (r = −0.9728), regional histogram bins (r = −0.9659), and minor axis (r = 0.8939) from nuclei morphology were the strongest parameters, as shown in Figure 2.
This result revealed that the stroma region has a significant relation to TC, which is in agreement with the results of the hand engineering method by [45] that was trained by a 125-dimensional feature vector reported a 0.75 ICC (first column of Table 2). Lower and upper-bounds are shown in square brackets. Our methodology is also a hand engineering, but trained with only 22 key morphological parameters (second column of Table 2), indicates a 0.78 ICC. This result outperformed those that were obtained by the method of Peikari. The results based on Deep Learning techniques like a combined hand engineering Deep Neural Network reported by Akbar [20] is slightly above the proposed methodology (third column of Table 2). Finally, the methodology proposed was used to process the test set and the results were submitted to the challenge contest. The prediction probability result obtained from contest was P k = 0.76.

Discussion
A computer methodology that automatically processes H&E histopathology digital images based in the extraction of main morphological parameters at a cell, regional, and global level is presented in this paper. The methodology processed a training set of breast cancer images under NAT treatment and the results indicate 22 key morphological parameters are strongly correlated with cellularity. Interesting results were revealed from the correlation analysis of the morphological parameters. The strongest related parameter was stroma density, in agreement with Beck et al. [16], which is, the histology of stroma correlates with prognostic in breast cancer. Three different machine learning algorithms for cell classification were evaluated and compared in order to determine tumour regions. The best result was obtained with Support Vector Machines (SVM) algorithm. The relevance of this paper is a selection of a key parameters to train the algorithms, which results in a better performance of similar techniques; however, the reported deep learning algorithms outperform this result, which is a motivation to explore these techniques in the future.  Acknowledgments: The authors would like to thank Sunnybrook Health Sciences Centre, Toronto, for the breast cancer NAT images. We also would like to thank Kenny H. Chan for validating our test results which were submitted to SPIE Breast Challenge contest.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.