Image processing in digital pathology: an opportunity to solve inter-batch variability of immunohistochemical staining

Immunohistochemistry (IHC) is a widely used technique in pathology to evidence protein expression in tissue samples. However, this staining technique is known for presenting inter-batch variations. Whole slide imaging in digital pathology offers a possibility to overcome this problem by means of image normalisation techniques. In the present paper we propose a methodology to objectively evaluate the need of image normalisation and to identify the best way to perform it. This methodology uses tissue microarray (TMA) materials and statistical analyses to evidence the possible variations occurring at colour and intensity levels as well as to evaluate the efficiency of image normalisation methods in correcting them. We applied our methodology to test different methods of image normalisation based on blind colour deconvolution that we adapted for IHC staining. These tests were carried out for different IHC experiments on different tissue types and targeting different proteins with different subcellular localisations. Our methodology enabled us to establish and to validate inter-batch normalization transforms which correct the non-relevant IHC staining variations. The normalised image series were then processed to extract coherent quantitative features characterising the IHC staining patterns.

Methods 2 and 3, respectively labelled "Li-init" and "Li-init-NMF", adapted from 11 . The initialisation step, "Li-init", is done by converting the RGB space into HSV. A saturation-weighted hue histogram is then computed and K-means clustering is used (with ! = 2) to find 2 cluster centres that correspond to the hue ( $ ) values of the stain colours. For each $ value so extracted from the histogram, corresponding % and & values are determined to provide HSV vectors by computing a saturationweighted mean of both the % and & values on a set of pixels close to the cluster centre. The two so extracted HSV vectors are converted back to RGB and then to SDA in the IHC context. To extract the final vectors, an NMF algorithm is then used to decompose non-negative matrix ' into two other nonnegative matrices, ( and ) (see equation (S1)). The complete process is labelled "Li-init-NMF". min ||' − ()|| / 0 with ( and ) ≥ 0 (S1) Methods 4 and 5, respectively labelled "SNMF" and "Li-init-SNMF", adapted from 16 . SNMF uses an objective function in which regularisation and sparsity terms are introduced (see equation (S2)). Xu et al. 16 propose to use 2.5 as a common value for B 7 , B 0 and C on the basis of an experimental and qualitative analysis. In our application to IHC images setting this constant resulted in null stain vectors with the two initialisation techniques mentioned in the main text (cf. subsection "Blind colour deconvolution"). After experimental tests we set the 3 parameters to 0.2 and tested this method with the two initialisation techniques, providing us two variants labelled "SNMF" and "Li-init-SNMF".

SDA distribution fitting
We tested two methods proposed in the literature 8,15 . In the first method 8 deconvoluted OD distributions are fitted by simply multiplying each value by the ratio of the 99 th percentiles (D99), i.e. D99 F,G D99 H,G , where I is the target image, J is the image to normalise with regard to I and K is the deconvoluted channel. Khan et al. 15 use B-Splines to fit intensity distributions (without OD transform) for each deconvoluted channel. The splines are fitted on 5 values, including 3 distribution parameters (5 th percentile, mean and 95 th percentile) and the theoretical minimum and maximum values (e.g. 0 and 255 in intensity values) to impose the identity at these extreme values. To apply this method on SDA directly we set the identity at 0 and −log 1 255. In our application to IHC images we use cubic smoothing splines to implement this method. We also generalise these two methods by either applying linear regression or B-spline regression on quantile-quantile (Q-Q) plots constituted of regularly distributed quantiles sampled from the %QR STG,H,G value distributions. We use the 100 percentiles for fitting linear regression whereas we fit cubic smoothing splines on 20, 30 and 100 quantiles regularly spaced and always beginning at the 1 st percentile.

Software and hardware specifications
NDPI images were generated using NDP.Scan 2.5, the proprietary software of Hamamatsu. Tiles were extracted using NDP.Read 1.1 (Hamamatsu) and were written to the PNG format using Python 2.7 along with scikit-image 0.11 or in the proprietary NDPI format using NDP.Write 1.1.2 (Hamamatsu). Like previously done for the method developed to identify each core in a TMA slide image 19 , the code of the proposed methodology was written in Python using numpy 1.10, scipy 0.17, scikit-learn 0.17, nimfa 1.2. The machine was a laptop running a 64-bits version of windows 10 with a Core I7-6700hq CPU, 16GB of ram and an SSD (solid-state drive) with a minimal overhead. In these conditions the computation time usually took 20 minutes to process a batch of 8 whole slides (scanned at 20x). The code could be optimised by rewriting it in another language or using the GPU. The overhead was totally manageable in our IHC image analysis workflow.

Supplementary information on immunohistochemistry
Technical, fixative and negative IHC controls were prepared on TMA slides using haematoxylin-eosin, anti-vimentin and (anti-mouse or anti-rabbit) IgG staining, respectively, as previously detailed 21 .
Immunostains were performed on Ventana discovery XT (Ventana, Roche Diagnostics, Vilvoorde, Belgium) using a Discovery DABMap Detection Kit (P53, PDGFRA and IGF2R) or a Discovery ChromoMap DAB Kit (ERG), both from Ventana. Briefly, deparaffinized and rehydrated (5µm-thick) tissue sections were incubated in EDTA buffer (pH 8.4) and endogenous peroxidase activity was inhibited using Endogenous Biotin Blocking Kit. After primary antibody incubation the slides were washed and incubated with a secondary antibody against mouse (for P53) or rabbit (for the other targets). The HRP binding secondary antibody was detected by incubation with diaminobenzidine (DAB) and hydrogen peroxide. Sections were counterstained with haematoxylin (HEM) and mounted. Additional details on the experimental conditions can be found in Table S1. Except for the CD3 experiment, the different IHC batches were carried out with the same lots for the antibodies and the kits.

Abcam, ab27478
Conditions used with the associated primary antibody Table S1: Details on the IHC assays Figure S1: Examples of IHC staining using DAB and HEM on a TMA core (partially shown): a) IGF2R in colon, b) ERG in liver, c) PDGFRa in glioblastoma, d) P53 in glioblastoma. e) Tiled image characterising an IHC experiment targeting IGF2R and extracted from a TMA slide composed of tissue samples from various origins. Figure S2: Application of the normalisation methodology, using "Macenko" for colour matching and a regression with 100 splines for SDA fitting, to reduce the variations between the CD3 batches. Illustration of CD3 expression evidenced in cores from a-b) the batch of reference, c-d) the batch to normalise and the corresponding images obtained after normalisation (colour and SDA fitting) based on a selection of 21 cores (e-f) or 43 cores (g-h).