Signature maps for automatic identification of prostate cancer from colorimetric analysis of H&E- and IHC-stained histopathological specimens

Prostate cancer (PCa) is a major cause of cancer death among men. The histopathological examination of post-surgical prostate specimens and manual annotation of PCa not only allow for detailed assessment of disease characteristics and extent, but also supply the ground truth for developing of computer-aided diagnosis (CAD) systems for PCa detection before definitive treatment. As manual cancer annotation is tedious and subjective, there have been a number of publications describing methods for automating the procedure via the analysis of digitized whole-slide images (WSIs). However, these studies have focused only on the analysis of WSIs stained with hematoxylin and eosin (H&E), even though there is additional information that could be obtained from immunohistochemical (IHC) staining. In this work, we propose a framework for automating the annotation of PCa that is based on automated colorimetric analysis of both H&E and IHC WSIs stained with a triple-antibody cocktail against high-molecular weight cytokeratin (HMWCK), p63, and α-methylacyl CoA racemase (AMACR). The analysis outputs were then used to train a regression model to estimate the distribution of cancerous epithelium within slides. The approach yielded an AUC of 0.951, sensitivity of 87.1%, and specificity of 90.7% as compared to slide-level annotations, and generalized well to cancers of all grades.

patient cohort. A total of 184 prostate specimens were obtained from a cohort of 63 patients who underwent radical prostatectomy for definitive treatment of biopsy-proven prostate adenocarcinoma at our institution between November 2009 and January 2012. A summary of the patient characteristics is detailed in Table 1.
Histopathology processing and staining. The prostate specimens were fixed and paraffin-embedded using a previously described protocol and sliced into 4 µm-thick axial sections 8,12 . From each tissue block, two sections were selected from tissue levels no more than 100 µm apart, and were de-paraffinized and rehydrated using standard methods. H&E and IHC staining was performed on the two sections, respectively. H&E staining was performed in three batches using routine clinical protocols. IHC staining was performed using a Ventana BenchMark ULTRA automated immunostainer platform (Ventana Medical Systems, Tucson, AZ). Antigen retrieval and blocking were performed as previously described 31 . Slides were incubated for 32 minutes with the triple-antibody cocktail containing primary antibodies to the basal cocktail of HMWCK + p63 (monoclonal mouse; clones 34βE12 and 4A4 respectively; prediluted; Ventana, Tucson, AZ) and AMACR (monoclonal rabbit; clone 13H4; prediluted; Dako, Glostrup, Denmark). Detection was performed with the Ventana ultraView Universal DAB Detection Kit and ultraView Universal Alkaline Phosphatase Red Detection Kit according to manufacturer's instructions. This was followed by rinsing, counterstaining with hematoxylin, dehydrating, and coverslipping. In summary, HMWCK + p63 expression in benign basal epithelium was demonstrated as brown by 3,3-diaminobenzidine (DAB), AMACR expression in malignant epithelium was demonstrated as red by Fast Red chromogen, and stroma was demonstrated as blue by hematoxylin counterstain. slide digitization and slide-level annotations. Both H&E and IHC slides were digitized at 20x magnification (0.5 µm/pixel) using a whole slide scanner (Aperio ScanScope CS, Leica Biosystems, Buffalo Grove, IL). Digitized H&E WSIs were annotated at the slide-level for PCa by pathology trainees (B.M.B., A.D.J., N.P.R.) under the supervision of a board-certified pathologist (S.C.S.) using Aperio's ImageScope software (Leica Biosystems, Buffalo Grove, IL) and a pen-tablet screen (Wacom Cintiq 22HD, Saitama, Japan). The slide-level annotations were carried out by demarcating the borders of distinct regions of cancer and assigning a Gleason score (GS) to each region (Fig. 1c). Using the same tools, negative annotations, defined as regions containing artifacts of the histological processing (e.g., tissue folds, debris, irregular staining), were demarcated on the IHC WSIs by www.nature.com/scientificreports www.nature.com/scientificreports/ technologists (A.E.R., J.C.H.). Regions of negative annotations were ultimately excluded from analysis, and typically comprised no more than 5% of a given slide. Digitized WSIs and annotations were stored and managed as previously described 32 .
SigMap software was used to further process the digitized WSIs 12,33 . First, it was used to register the IHC WSI to the H&E WSI using a rigid transformation (Fig. 1a,b). Next, binary masks of the slide-level cancer annotations and the negative annotations were created to transfer the annotations between H&E WSIs and IHC WSIs (Fig. 1c). A virtual grid composed of analysis squares (dimensions 1,000 × 1,000 pixels, area of 0.25 mm 2 ) was then generated by SigMap and added to both WSIs (Fig. 1d). Analysis squares whose areas overlapped at least 75% with the cancer annotation mask were labeled as cancer and assigned the GS of the corresponding annotation (Fig. 1e). Analysis squares whose areas overlapped at least 75% with the negative annotation mask were excluded from further analysis.
Colorimetric image analysis algorithms. The following three quantitative image analysis algorithms (Aperio Brightfield Image Analysis Toolbox, Leica Biosystems, Buffalo Grove, IL) were configured by a technologist (J.C.H.), then applied to H&E and IHC WSIs in order to extract features for prediction of cancer.
The Positive Pixel Count (PPC) algorithm was applied to H&E WSIs. Briefly, the PPC algorithm counts the number of stained pixels within each analysis square that falls within and out of a specified range of hue-saturation-brightness (HSB) color values (positive and negative pixels, respectively). HSB values were sampled from three types of regions that predominantly contained a single histological feature of interest (nuclei, cytoplasm, or stroma). Fifteen of each type of region were manually identified on control H&E WSIs and sampled. Ranges of HSB values were calculated for each type of region and were manually adjusted to eliminate overlap between ranges. A separate PPC algorithm was configured for each type of region and its corresponding range of HSB values. The three configured PPC algorithms were then applied prospectively to analysis squares of H&E WSIs. The resulting numbers of positive pixels were converted to percentages of the analysis square occupied by nuclei, cytoplasm, and stroma (% nuclei, % cytoplasm, and % stroma, respectively), which were in turn used as predictive features. The unstained percentage of each analysis square was also calculated as % unstained = 100%− (% nuclei + % cytoplasm + % stroma), and analysis squares with % unstained >99% were excluded from further analysis on the basis that they are taken from regions outside of the tissue boundaries. To account for variations in H&E staining intensity across the three batches, a different set of PPC algorithms was configured and applied to each batch.
Color Deconvolution (CD) and Co-expression (CE) algorithms were applied to IHC WSIs to measure the colorimetric features of the IHC stain. Briefly, the CD algorithm isolates individual staining components of IHC WSIs for quantification, while the CE algorithm quantifies how often the staining components occur separately and together. These algorithms were first configured on control slides. Three control slides were cut, processed, and singly-stained with either DAB chromogen (brown), Fast Red chromogen (red), or hematoxylin counterstain (blue), using the same protocols as the triple-stained IHC slides described above. The average red-green-blue (RGB) optical density (OD) values of the three components were sampled from the corresponding WSIs of the control slides and were measured as Fast Red (R: 0.283, G: 0.949, B: 0.757), DAB (R: 0.461, G: 0.826, B: 1.0), and hematoxylin (R: 0.21, G: 0.276, B: 0.176), and intensity thresholds were manually configured for each component to define positively-stained pixels. The configured CD and CE algorithms were then applied prospectively to analysis squares of IHC WSIs, from which the percentage of each analysis square that was positively staining (%Pos) was calculated. As previously described, the OD quantifies the stain intensity, as it is linearly related to the amount of staining 13 www.nature.com/scientificreports www.nature.com/scientificreports/ Using the configured RGB OD and intensity threshold values, IHC WSIs were then separated into brown, red, and blue color channels corresponding to each staining component. The brown and red staining were separately quantified by the CD algorithm as previously described 32 . Specifically, the average OD and %Pos were measured by the CD algorithm for both brown and red components, and the products OD × %Pos were calculated and used as predictive features. The co-localization of brown and red staining was quantified by the CE algorithm, which was then used to calculate the percentage of the analysis square that was positively staining for only red or only brown, but not both (%Pos CE ). %Pos CE for red and brown components were used as predictive features.
In summary, seven features were extracted from each analysis square ( Table 2). The features derived from H&E WSIs were the percentages of nuclei, cytoplasm, and stroma, while the features derived from IHC WSIs were the percentages and stain intensities (quantified by the OD) of brown and red staining, which corresponded to the characteristics of the basal cell staining (HMWCK + p63) and the AMACR staining, respectively. training data and analysis square-level annotations. Ten of the 63 patients in our cohort were randomly selected, and one pair of WSIs was created from each for purposes of training the regression model. Forty analysis squares were randomly selected from each of the ten pairs of WSIs (400 analysis squares in total) and were manually annotated in much greater detail than usual (S.C.S.). The analysis square-level annotations were carried out by meticulously delineating the benign and cancerous epithelium, gland lumens, stroma, and regions of clear glass within the 1,000 × 1,000 pixel-area of each analysis square. The fractional areas of each of the aforementioned components were then summated for each annotated analysis square (Fig. 2b). The percentage of cancerous epithelium within each analysis square-level annotation was taken to be the ground truth. Details on the slides can be found in Table 3.

Regression model training and evaluation.
Elastic net regression models were trained on these data using 10-fold cross validation 35 , with each fold containing the 40 analysis squares from a single pair of WSIs. The elastic net is a generalized linear regression model with both L1 and L2 regularization, and its corresponding objective function to be minimized is    For each model, the coefficients of the two regularization terms (α and ρ) were treated as hyperparameters and selected by cross-validation to minimize the mean value of the objective function (averaged across the ten folds). Trained models were then applied to the analysis squares of the other 174 pairs of slides to produce predicted maps of cancerous epithelium. Model outputs were compared to the slide-level annotations on a per-analysis   www.nature.com/scientificreports www.nature.com/scientificreports/ square level using receiver operating characteristic (ROC) curve analysis. Sensitivities and specificities were calculated using the optimum cut-off points for the ROC curves that corresponded to the maxima of the Youden indices. The 95% confidence intervals (CIs) were calculated using a bootstrap procedure that resampled both WSIs and analysis squares from the training set, and only WSIs from the test set. Two-sided p-values were found by inverting the 95% bootstrap CIs.  Figure 2g,h show sample outputs of the CD and CE algorithms, respectively. As the degree of co-localization between the brown and red components of the deconvolved IHC slides is generally small (typically <5% of the analysis square), the overlap is difficult to visualize in Fig. 2i. Quantitative evaluation of model performance. Cross-validation performance for the four models was evaluated by plotting cumulative scatterplots of predicted % cancer epithelium vs. the actual % cancer epithelium across the 10 cross-validation folds (Fig. 3). The cross-validation root mean square error, median absolute error, and maximum absolute error for each model are shown in Table 4.

outputs of colorimetric image analysis algorithms.
Performance on the test set for the four models was evaluated by plotting the ROC curves (Fig. 4). The area under the ROC curves (AUCs), as well as the sensitivities and specificities calculated at the respective maxima of the Youden indices, are shown in Table 5. The AUC for the full model was significantly higher than that of the H&E model (p = 0.026), while the specificity for the full model was significantly higher than those of the H&E and full -CE models (p < 0.001 for both). The AUC and specificity for the full model were not significantly different than those of the IHC model (p = 0.542 and p = 0.108, respectively). The sensitivity of the full model was also not significantly different than those of the H&E, IHC, and full -CE models (p = 0.134, p = 0.748, and p = 0.939, respectively). The CIs of these summary statistics for the full model were notably narrower than those of the other three models, suggesting that the performance of the full model will likely be closer to what is reported here when it is applied prospectively.     36 , the latter of which may better reflect cancer aggressiveness (Table 6). Using the convention that GG ≤ 2 (GS 3 + 3 or 3 + 4) is low to intermediate-grade and GG ≥ 3 (GS 4 + 3, 4 + 4, 4 + 5, or 5 + 4) is high-grade, the sensitivity of detecting low to intermediate-grade cancers was 0.884, while it was 0.864 for high-grade cancers; this difference was found to be not significant (p = 0.107). Figure 5 shows representative H&E and IHC slides with slide-level, manually-annotated cancer by pathologists compared with maps generated by the full model. Note the high degree of correlation between the annotated cancer, distribution of AMACR staining (red on IHC slides), and predicted distribution of malignant epithelium.

Discussion
In this work, we show that a predictive model that uses features derived from colorimetric analysis of both digitized H&E and IHC slides is able to detect and delineate PCa on WSIs with accuracy comparable to pathologists' slide-level annotations. The performance of the full model was found to be superior to those of the other three models, individually (Table 5), though this difference was only significant in comparison to the H&E model. Furthermore, despite the relatively small amount of training data that included predominantly low and intermediate-grade cancers, the model performed well across a large number of test set slides. Its sensitivity was also largely consistent across cancers with different Gleason scores, and was only slightly worse for high-grade cancers (0.864 vs. 0.871 for all cancers, Table 6).
In contrast to most published works, the regression model described here uses a compact set of seven features that were calculated from the outputs of standard image analysis algorithms applied to H&E WSIs (% nuclei, % cytoplasm, % stroma) and IHC WSIs (OD × %Pos and %Pos CE , for brown and red). These features were ultimately chosen for their simplicity and interpretability. Although %Pos from the CD algorithm and %Pos CE from the CE algorithm appeared to be redundant, the results demonstrate that excluding the two %Pos CE features resulted in worse specificity of cancer detection (full −CE model vs. full model, Table 5). This is most likely due the fact that some slides contained a larger fraction of non-cellular components (e.g., intraglandular cellular debris, corpora amylacea) that stained both brown and red with IHC staining. This would cause %Pos (red) as calculated by the CD algorithm to be falsely elevated, but would not affect %Pos CE (red) as calculated by the CE algorithm as the CE algorithm excludes regions that stained both brown and red. Therefore, inclusion of the two %Pos CE features increased the specificity of cancer detection for slides containing a significant fraction of such regions, and in turn increased the overall specificity of cancer detection.
Another notable aspect of the work is that the full model was trained to predict the percentage of cancerous epithelium, which was made possible by the unique ground truth obtained from the meticulous annotation of individual analysis squares. For purposes of model training, this ground truth is superior to the slide-level annotations, as those are known to have finite accuracy and precision 9,10 .
The trained models were evaluated against the slide-level annotations on a per-analysis square level using ROC curve analysis. However, despite the high AUC achieved by the full model, it is difficult to assess the true performance of the model due again to the limitations of the slide-level annotations. Accurate assessment would require analysis square-level annotations of all the slides in the test set, which would be prohibitive. Qualitatively, visual comparison of the model-generated maps of cancerous epithelium with the slide-level annotations shows generally good concordance (Fig. 5). Sources of disagreements between the two can be divided into four categories, which are illustrated in Fig. 6: 1. Cancer missed by the model (Fig. 6a). This was most often due to cancer with poor AMACR staining; while AMACR is a sensitive positive marker of PCa, it is well-documented that some variants of PCa do not exhibit increased expression of AMACR 29,30,37,38 . Alternatively, inconsistencies in the staining procedure may have caused variabilities in AMACR staining, and these variabilities would be amplified in regions of cancer. www.nature.com/scientificreports www.nature.com/scientificreports/ 2. Cancer incorrectly annotated by the pathologist (Fig. 6b). This was most often due to large regions of glass (e.g., cystic areas, luminal areas of malignant glands) that are looped in with the slide-level annotations. More rarely, benign glands were incorrectly annotated as cancer; usually, these were examples of PIN with low AMACR expression. 3. Cancer incorrectly labeled by the model (Fig. 6c). This was most often due to PIN with high AMACR expression. 4. Cancer missed by the pathologist (Fig. 6d). This was most often due to small, isolated regions of cancer that were not annotated. More rarely, HGPIN and/or glands with IDC-P were missed by the pathologist but identified as cancer by the model due to high AMACR expression.
In summary, accurately distinguishing the different possible presentations of PCa from PIN is a challenge for both pathologists and the full predictive model. Although in theory glands with PIN are characterized by the presence of an intact basal cell layer, the basal cell layer may be quite fragmented, which would make it difficult to assess by either visual inspection or by quantitative assessment of brown staining intensity.
There are three major limitations of the features. First, since they are calculated from colorimetric analysis of stained WSIs, their consistency is highly-dependent on the reproducibility of the staining procedure and digitization process. As noted in previous works, the use of different histology protocols and/or different slide scanners can cause large variations in the morphological features of WSIs, which in turn degrade the predictive performance of trained models 20 . In our work, the three separate batches of H&E staining presented a major source of potential variability in the calculated H&E features, as the stain intensities were visibly different between H&E WSIs of different batches. To compensate, a different set of PPC algorithms was configured for each batch, though this was not ideal. In order to minimize batch effects in the future, H&E staining will also be performed on an automated platform using a standardized protocol, like what was done for the IHC staining. Additionally, algorithms like the PPC algorithm that rely on the analysis of intensity values (e.g., RGB or HSB values) are naturally quite prone to variations in stain intensity. Therefore, it would be worth extending the use of color deconvolution algorithms for the analysis of H&E WSIs, as proposed in previous works 39 . To expand this work to larger datasets, further methods for normalization of WSIs and/or derived features may also be investigated 20,39 . www.nature.com/scientificreports www.nature.com/scientificreports/ Another limitation of the features is that they are derived from WSIs of different tissue levels of the tissue block. Due to technological constraints, H&E and IHC staining were not performed on the same tissue section, and thus two sections were taken from each specimen. Although the sections were spatially adjacent (separated by ≤100 µm), there were sometimes noticeable differences between the two when digitized and viewed at the magnification level of individual analysis squares (Figs 2a,f, 6b). However, these differences were relatively minor and unlikely to significantly affect the calculated features, and can further be minimized in the future by always selecting serial levels of the tissue block for staining. Methods for visualizing multiple stains on the same tissue section may also be considered 40, 41 .
Lastly, the features only characterize the composition within each analysis square, and not the arrangement (i.e., cellular architecture) of the components. Therefore, the predictive model has difficulties distinguishing between analysis squares containing PIN and those containing a mixture of benign and malignant glands. This limitation may be addressed in future work by identifying additional IHC markers that are differentially expressed in cancer and PIN; for example, IDC-P is characterized by decreased expression of PTEN, which can be used to distinguish HGPIN from IDC-P 42 . An alternative approach could be to develop custom algorithms for object detection and segmentation (e.g., for identification of whole prostate glands). A more straightforward approach could be to supplement the training data with examples of PIN or PIN-like entities. Augmenting the training data may also allow the use of deep learning approaches such as convolutional neural networks [15][16][17][18][19] that can learn features that account for differences in the glandular architecture within analysis squares.
In summary, the methods introduced in this work can be modularly integrated into digital pathology frameworks for detection of prostate cancer on whole-slide images of histopathology slides. The unique aspect of this work is that it incorporates information from slides with conventional H&E staining as well as those with IHC staining, and as demonstrated in this work, the combination of both allows for more accurate identification of prostate cancer. Given the number of previously identified and characterized genetic markers in other types of cancers, the methods presented here may be extended naturally to other types of cancer as well.

Data Availability
The processed digitized pathology data, the outputs of the colorimetric analysis algorithms and the regression model, and the statistical analyses for the current study are available from the corresponding author upon reasonable request. Figure 6. Illustrative examples of analysis squares with discrepancies between the slide-level annotation and the model output. The top row shows H&E analysis squares, and the bottom row shows the corresponding IHC analysis squares. (a) Analysis square with predominantly malignant glands that have poor AMACR staining. 100% overlap with the slide-level annotation, but incorrectly labeled as non-PCa by the model. (b) Analysis square with PIN that has poor AMACR staining (blue arrow) and a large cystic region (bottom half). Correctly labeled non-PCa by the model, but had 100% overlap with slide-level annotation. (c) Analysis square with PIN that has strong AMACR staining. Did not overlap with the slide-level annotation, but was incorrectly labeled as PCa by the model. (d) Analysis square with small malignant glands (green arrows), including an example of HGPIN/IDC-P (yellow arrows). Correctly labeled PCa by the model, but did not overlap with the slide-level annotation.