Robust Method for Unsupervised Scoring of Inmunohistoquemical Staining

1. Introduction

Detection of biomarker expression in tissue images is a common technique both in research laboratories and in the clinic. This detection is carried out through the use of chromogens of different colors that bind to the antigens of interest using an antibody-antigen detection system, a technique known as Immunohistochemistry (IHC) [1]. This allows the presence of the chromogen and the protein to which it binds by observing the tissue using bright field microscopy. Capturing images of the tissue allows the subsequent preservation and analysis of the tissue expression pattern of analyzed proteins. Usually, the quantification of the presence, concentration and distribution of the chromogen is carried out manually by viewing and comparing the images by experts.

Two widely used chromogens are 3,3’-Diaminobenzidine (DAB) and Hematoxylin (H) [2,3,4,5,6]. While the presence of DAB is detected as a brown stain, H appears as a blue one. When combined, H stain marks cells nuclei, whereas DAB reveals the greater or lesser presence of the protein of interest. The IHC image evaluation and quantification procedure lacks objectivity, as it is a subjective method done by experts, exposed to observer variations. Indeed, even the same observer can make different evaluations of the same image at different times [7]. Furthermore, manual scoring is affected by visual distortion due to the perception of colors and their intensity depending on the surrounding colors [8]. This is particularly important in protocols that use two dyes. Finally, IHC visual scoring it is a time-consuming procedure for experts. For this reason, demand for automated scoring procedures is increasing in laboratories and in the clinic [9]. This type of protocols must include a stage for separating both stains to perform the scoring with greater reliability.

Scoring systems can be classified on those based on the perceived intensity of staining, and those based on the percentages of stained cells, while some others mix both criteria. Regarding the most used categories, they range from the simplest systems that divide tissue samples into positive or negative, depending on whether or not the observed staining exceeds certain thresholds; to those that use a set of scoring levels depending on the greater or lesser presence of stain, the so-called semi-quantitative scoring systems [8]. Although some platforms and applications for semiquantitative scoring of IHC images already exist, they require the intervention of researchers as well as reference images. Examples of these systems are QuPath [10], IHC Profiler [11] or DAB-quant [12].

The stage of stain separation is necessary in automated scoring systems, since it allows the obtention of image features based on each stain. The separation stage cannot be performed in the RGB color space, because the concentration of the light-absorbing material follows the Beer-Lambert law (although not exactly, due to some degree of dispersion). For this reason, a transform based on this law is applied to the RGB image, obtaining the coordinates of each pixel in a new space: the optical density (OD). Within this space, the intensity of each coordinate can be considered as a linear combination of the intensities due to each staining [13].

Separation methods in the literature can be classified into two categories: supervised and unsupervised (or blind) methods [14,15]. In the first case, the coordinates of each used stain are known in the OD space; otherwise, they require the intervention of the observer or the use of reference images to complete their calculation. In the second case, unsupervised methods use techniques like independent component analysis (ICA) [16,17], non-negative matrix factorization (NMF) [13,18] or non-negative least squares (NNLS). These techniques suffer from a strong dependence on initialization. Thus, these methods need a previous estimation of the stain vectors coordinates to achieve a correct separation, guiding the search towards the correct solution. Some approaches which are based on deep learning [4,19] need a large number of training images, so they cannot be considered as unsupervised techniques.

Some approaches have exploited singular value decomposition (SVD) in order to calculate the subspace where the OD data lie [20,21] in the context of supervised methods. Principal component analysis (PCA) works in the same direction as SVD, searching for the principal direction of the data: those that preserve the greatest quantity of variance or power (when uncentred PCA is performed) [22].

In this work, we propose an unsupervised method for comparative biomarker quantification that improves the method presented in our previous work [14]. Accordingly, this new methodology represents a simpler and more robust method for color separation, with very low dependence on the selected parameters, and a more accurate quantification without the need for reference or control images.

We circumvent the use of NMF for the stain separation stage, whose results are very dependent on the starting point. Instead, the proposed method exploits the concentration of DAB (brown) and H (blue) to estimate the vectors associated to the color of both stains. This stage has been provided with robustness, in order to achieve consistent results for the estimation of the stain vectors with no prior information, even avoiding the use of reference images that the method in [14] needs.

For the scoring stage, features whose values mostly increase with the concentration of the stains helped us to develop a procedure for the initialization of the clustering algorithm (k-means). While in our previous work a set of three parameters, together with the use of two reference images were needed, in the present work the initialization procedure provides excellent results for common values of the two parameters that we use and, again, avoiding the use of reference images. This makes totally unnecessary any intervention of the observer and allows a totally automated procedure.

2. Methods

We developed an automated and unsupervised procedure for the scoring of IHC images with two stains: Hematoxylin and DAB. The quantification of the set of images is done in three steps: stain separation, feature extraction and, finally, clustering and scoring. A semi-quantitative scoring system clusters the images in five groups and assigns scores from “1+” to “5+” for each image, being “1+” the lowest level of DAB expression.

2.1. Stain separation

The stain separation step is carried out within the color space given by the transform provided by the Beer-Lambert Law. Let ${\mathbf{I}}_m\in {\mathbb{R}}^{3\times N}$ be the matrix of RGB intensities of the mth analyzed image, where N is the total number of pixels of each color plane. Then the relative optical density matrix for this image, ${\mathbf{X}}_m$, is obtained as

$${\mathbf{X}}_m=-log\left(\frac{{\mathbf{I}}_m}{I_o}\right),$$

(1)

where $I_o$ is the maximum level of intensity at each pixel (255 for 8 bit images) and the function $log\left({\mathbf{I}}_m\right)$ is applied to each entry of the matrix ${\mathbf{I}}_m$. This provides a representation of each pixel in a ${\mathbb{R}}^3$ space where each color corresponds to a direction, being a higher darkness of the color represented by a higher norm of the vector. So, each color (including its darkness) can be represented by the coordinates of its corresponding vector within this space. Since the image was obtained by means of an IHC procedure, the color vector of each pixel lies, mainly, within the subspace generated by the color vectors associated to both stains: that due to H (blue) and that due to DAB (brown). This allows to estimate the matrix ${\mathbf{X}}_m$ as [13]

$${\widehat{\mathbf{X}}}_m={\mathbf{W}}_m{\mathbf{H}}_m,$$

(2)

where ${\widehat{\mathbf{X}}}_m$ is the matrix that estimates ${\mathbf{X}}_m$, ${\mathbf{W}}_m$ is a $3\times 2$ matrix whose columns are the color vectors associated to both stains for mth image and ${\mathbf{H}}_m$ is a $2\times N$ matrix whose rows are the intensities or activations of each column of ${\mathbf{W}}_m$ at each pixel of the mth image.

We exploit this fact in order to simplify the problem by means of dimension reduction. The subspace generated by the columns of ${\mathbf{W}}_m$ is found by means of the application of uncentred Principal Component Analysis (PCA) to the matrix ${\mathbf{X}}_m$. This provides the principal directions of the data, that lie within the above mentioned subspace (although the principal directions do not have to coincide with the columns of ${\mathbf{W}}_m$) as can be seen in Figure 1.

Let

$${\mathbf{Q}}_m=\left[\begin{array}{ccc}{{\mathbf{q}}_m}_1& {{\mathbf{q}}_m}_2& {{\mathbf{q}}_m}_3\end{array}\right]$$

(3)

be the unitary matrix whose columns are the eigenvectors of ${\mathbf{R}}_m={\mathbf{X}}_m{\mathbf{X}}_m^T/N$ and

$${\mathbf{L}}_m=\mathrm{diag}\{{L_m}_1,{L_m}_2,{L_m}_3\}$$

(4)

be the diagonal matrix whose nonzero entries are the eigenvalues of ${\mathbf{R}}_m$. The superscript T stands for transpose. Both the eigenvectors and the eigenvalues are in descending order. Since ${\mathbf{X}}_{\mathbf{m}}\approx {\widehat{\mathbf{X}}}_m$, the principal subspace is given by the two first columns of ${\mathbf{Q}}_m$, as the third eigenvalue, ${L_m}_3$, is very small compared to ${L_m}_1$ and ${L_m}_2$. The columns of ${\mathbf{X}}_m$ are projected onto the principal subspace (the plain defined by the first two principal direction), thus obtaining the $2\times N$ matrix of principal components

$${\mathbf{Z}}_m=\left[\begin{array}{c}\sqrt{{{L_m}_1}^{-1}}·{{\mathbf{q}}_m^T}_1\\ \\ \sqrt{{{L_m}_2}^{-1}}·{{\mathbf{q}}_m^T}_2\end{array}\right]·{\mathbf{X}}_m.$$

(5)

Since all columns of ${\mathbf{X}}_m$ lie at the first octant, the color vectors are relatively close at ${\mathbb{R}}^3$ Euclidean space. A simple projection onto the ${\mathbb{R}}^2$ subspace defined by ${{\mathbf{q}}_m}_1$ and ${{\mathbf{q}}_m}_2$, keeping the scales of the principal components would lead to a 2D scattering where it would be difficult to precisely distinguish the color vectors associated to each stain. This is the reason why the columns of ${\mathbf{X}}_m$ are projected onto the subspace defined by ${{\mathbf{q}}_m}_1\sqrt{{L_m}_1}$ and ${{\mathbf{q}}_m}_2\sqrt{{L_m}_2}$, i.e., PCA with unit-power principal components is performed [22]. This expands the scattering on the second principal component direction, as we show in Figure 2, and allows to recover the stain color vectors more accurately.

2.1.1. Stain color basis estimation

The next step of the procedure consists in the estimation of the two basis vectors for the stain colors, i.e., the columns of the matrix ${\mathbf{W}}_m$. Due to the projection in (5) the columns of this matrix are projected into the same subspace defined by the principal directions

$${\mathbf{B}}_m=\left[\begin{array}{c}\sqrt{{{L_m}_1}^{-1}}·{{\mathbf{q}}_m^T}_1\\ \\ \sqrt{{{L_m}_2}^{-1}}·{{\mathbf{q}}_m^T}_2\end{array}\right]·{\mathbf{W}}_m.$$

(6)

Let ${\theta }_m\left(n\right)$ and $r_m\left(n\right)$ be the angle and norm of the nth column of ${\mathbf{Z}}_m$. Since, in many pixels, only one stain is dominant, the values of $r_m\left(n\right)$ are higher around the basis vectors. So the columns of ${\mathbf{B}}_m$ are estimated by means of the angles where large concentrations of high levels of $r_m\left(n\right)$ are found. The range of values of ${\theta }_m\left(n\right)$ is divided into 1000 bins for a sufficient angle resolution. For the kth bin, whose center is given by ${\varphi }_m\left(k\right)$, the mean of the values of $r_m\left(n\right)$ such that ${\theta }_m\left(n\right)$ lies within the bin is computed. Then, this function is smoothed by a lowpass eighth-order Butterworth filter with digital cutoff frequency equal to $0.035$ to get the function $v_m\left(k\right)$, as illustrated in Figure 3.

Therefore, ${\mathbf{B}}_m$ is estimated by means of ${\widehat{\mathbf{B}}}_m$, defined as

$${\widehat{\mathbf{B}}}_m=\left[\begin{array}{cc}cos\left({{\varphi }_m}_1\right)& cos\left({{\varphi }_m}_2\right)\\ sin\left({{\varphi }_m}_1\right)& sin\left({{\varphi }_m}_2\right)\end{array}\right]$$

(7)

where ${{\varphi }_m}_1$ and ${{\varphi }_m}_2$ are the main peaks in the filtered version of $r_m\left(n\right)$. Since it is possible the existence of more than two peaks, in order to avoid a bad selection of the true angles of both basis vectors, a procedure was implemented. Due to the nature of the data, the basis vectors are always near the extrema of ${\varphi }_m\left(n\right)$ (in an ideal case, the points of the scattering are a linear combination of the basis vectors with nonnegative coefficients). So, the maximum of $v_m\left(k\right)$ is selected as the position for one the basis vectors. The other peak is selected from the angle bins that lies within the opposite half of the total range of ${\varphi }_m\left(k\right)$. This ensures that the selected angles correspond to a good estimate of the basis vectors.

It is easy to see that the estimate of ${\mathbf{W}}_m$ corresponding to ${\widehat{\mathbf{B}}}_m$ is given by the representation of the columns of this matrix in the whole 3D space as

$${\widehat{\mathbf{W}}}_m=\left[\begin{array}{cc}{{\mathbf{q}}_m}_1\sqrt{{L_m}_1}& {{\mathbf{q}}_m}_2\sqrt{{L_m}_2}\end{array}\right]·{\widehat{\mathbf{B}}}_m.$$

(8)

2.1.2. Color deconvolution

From (2) we can write the decomposition of ${\mathbf{X}}_m$ in terms of ${\widehat{\mathbf{W}}}_m$ and ${\mathbf{H}}_m$ as

$${\mathbf{X}}_m={\widehat{\mathbf{W}}}_m·{\mathbf{H}}_m.$$

(9)

This describes, for each column of ${\mathbf{H}}_m$, an overdetermined system of equations as ${\mathbf{W}}_m$ and ${\mathbf{X}}_m$ are known. Therefore, the stain concentration matrix for the mth is estimated by means of a linear least squares problem as

$${\widehat{\mathbf{H}}}_m={\widehat{\mathbf{W}}}_m^{+}{\mathbf{X}}_m$$

(10)

where ${\widehat{\mathbf{W}}}_m^{+}={\left({\widehat{\mathbf{W}}}_m^T{\widehat{\mathbf{W}}}_m\right)}^{-1}·{\widehat{\mathbf{W}}}_m$ is the Moore-Penrose pseudoinverse of the full-rank matrix ${\widehat{\mathbf{W}}}_m$ [23]. This procedure is known as color deconvolution [24]. Since the concentration of each stain at every pixel must be nonnegative, the values ${\left[{\widehat{\mathbf{H}}}_m\right]}_{ij}<0$ are set to 0.

The estimate is improved by means of the following procedure. The values ${\left[{\widehat{\mathbf{H}}}_m\right]}_{ij}$ for columns of ${\widehat{\mathbf{H}}}_m$ that have only one nonzero entry are recalculated by using the corresponding column of ${\widehat{\mathbf{W}}}_m$. After this, a new estimate of the basis vectors is computed as the mean of the columns of ${\mathbf{Z}}_m$ for which one stain is dominant, i.e., when the value for a concentration is, at least, 10 times that for the other. A new estimate of ${\widehat{\mathbf{H}}}_m$ is computed following the procedure described from equation (10).

2.1.3. Average basis vectors

Once the previous procedure is carried out, the average of the matrices ${\widehat{\mathbf{W}}}_m$ is computed

$$\overline{\mathbf{W}}={〈{\widehat{\mathbf{W}}}_m〉}_m.$$

(11)

Then the concentration matrix is estimated for each image as ${\overline{\mathbf{H}}}_m={\overline{\mathbf{W}}}^{+}{\mathbf{X}}_m$. Again, the values ${\left[{\overline{\mathbf{H}}}_m\right]}_{ij}<0$ are set to 0 and then the values ${\left[{\overline{\mathbf{H}}}_m\right]}_{ij}$ for columns of ${\overline{\mathbf{H}}}_m$ that have only one nonzero entry are recalculated.

2.2. Feature extraction

A set of four features is computed for each image. Since the columns of ${\widehat{\mathbf{W}}}_m$ and $\overline{\mathbf{W}}$ were not normalized, we take into account them together with ${\widehat{\mathbf{H}}}_m$ and ${\overline{\mathbf{H}}}_m$. The jth column of ${\widehat{\mathbf{W}}}_m$ (resp. $\overline{\mathbf{W}}$) is denoted as ${\widehat{\mathbf{W}}}_{m_{:,j}}$ (resp. ${\overline{\mathbf{W}}}_{:,j}$) and the ith row of ${\widehat{\mathbf{H}}}_m$ (resp. ${\overline{\mathbf{H}}}_m$) is denoted as ${\widehat{\mathbf{H}}}_{m_{i,:}}$ (resp. ${\overline{\mathbf{H}}}_{m_{i,:}}$). Let $M_1\left(m\right)$ be the mean of the the intensity average of all the pixels of DAB stain for the mth image according to the results for the average basis vectors

$$M_1\left(m\right)=\frac{1}{3N}{∥{\overline{\mathbf{W}}}_{:,2}·{\overline{\mathbf{H}}}_{m_{2,:}}∥}_1=\frac{1}{3N}{∥{\overline{\mathbf{W}}}_{:,2}∥}_1·{∥{\overline{\mathbf{H}}}_{m_{2,:}}∥}_1=\frac{1}{3N}\sum _{i=1}^3{\left[\overline{\mathbf{W}}\right]}_{i2}\sum _{j=1}^N{\left[{\overline{\mathbf{H}}}_m\right]}_{2j}$$

(12)

where the 1-norm of the matrix is computed entry-wise, i.e., it is the sum of the entries of the matrix (they are all nonnegative). For row vectors we computed the norm as for column vectors. Then the first feature for the mth image is defined as the average of the DAB intensity in this image, normalized by the maximum value that this quantity reaches among the entire set of images

$$f_1\left(m\right)=\frac{M_1\left(m\right)}{{\displaystyle \underset{m}{max}}\left\{M_1\left(m\right)\right\}}.$$

(13)

The second feature is based on the relation between the robust maxima (99th percentile) of both stains: H and DAB. Let $R_1\left(m\right)$ (resp. $R_2\left(m\right)$) be the robust maximum of the intensities matrix for H stain (resp. DAB stain), according to the results for the average basis vectors, i.e., the robust maximum of the matrix ${\overline{\mathbf{W}}}_{:,1}·{\overline{\mathbf{H}}}_{m_{1,:}}$ (resp. ${\overline{\mathbf{W}}}_{:,2}·{\overline{\mathbf{H}}}_{m_{2,:}}$). Then, we define

$${\overline{R}}_2\left(m\right)=\frac{R_2\left(m\right)}{{\displaystyle \underset{m}{max}}{R}_2\left(m\right)}$$

(14)

and

$${\overline{R}}_1\left(m\right)=\frac{R_1\left(m\right)}{R_1(m=arg{max}_m{R}_2\left(m\right))}.$$

(15)

The second feature is then defined as

$$f_2\left(m\right)=\frac{{\overline{R}}_2\left(m\right)}{{\overline{R}}_1\left(m\right)}.$$

(16)

The third feature is analogue to the first one, but using 2-norm, i.e., based on the average power of the intensity. Let $M_2\left(m\right)$ be

$$\begin{array}{cc}\hfill M_3\left(m\right)=\frac{1}{\sqrt{3N}}{∥{\overline{\mathbf{W}}}_{:,2}·{\overline{\mathbf{H}}}_{m_{2,:}}∥}_F=\frac{1}{\sqrt{3N}}{∥{\overline{\mathbf{W}}}_{:,2}∥}_2·{∥{\overline{\mathbf{H}}}_{m_{2,:}}∥}_2\hfill & \\ & \hfill =\frac{1}{\sqrt{3N}}\sqrt{{\displaystyle \sum _{i=1}^3}{\left[\overline{\mathbf{W}}\right]}_{i2}^2}\sqrt{{\displaystyle \sum _{j=1}^N}{\left[{\overline{\mathbf{H}}}_m\right]}_{2j}^2}\end{array}$$

(17)

where ${∥·∥}_F$ is Frobenius norm of the matrix. Then we define the third feature, $f_3\left(m\right)$ as

$$f_3\left(m\right)=\frac{M_3\left(m\right)}{{\displaystyle \underset{m}{max}}\left\{M_3\left(m\right)\right\}}.$$

(18)

Although very goods results are achieved when using the above defined features, the inclusion of a fourth one improves the behavior of the clustering procedure. The feature $f_4\left(m\right)$ is defined as the square root of the average power of the intensity of the DAB stain at every pixel according to the results for the basis vectors estimated for each image,

$$f_4\left(m\right)=\frac{1}{\sqrt{3N}}{∥{\widehat{\mathbf{W}}}_{m_{:,2}}·{\widehat{\mathbf{H}}}_{m_{2,:}}∥}_F=\frac{1}{\sqrt{3N}}{∥{\widehat{\mathbf{W}}}_{m_{:,2}}∥}_2·{∥{\widehat{\mathbf{H}}}_{m_{2,:}}∥}_2$$

(19)

2.3. Clustering and scoring

As a general behavior the values of the four defined features increase when the concentration of DAB stain increases. Consequently, a simple, automated an efficient procedure was used to initialize the algorithm k-means for clustering. Also, this property allows us to assign the scores in a increasing order of the centroids norms. Let $\mathbf{f}\left(m\right)$ be the vector of features whose ith entry is $f_i\left(m\right)$. In order to initialize the clustering algorithm we computed the p-norm of the vector

$${∥\mathbf{f}\left(m\right)∥}_p={\left(\sum _{i=1}^T{\left|f_i\left(m\right)\right|}^p\right)}^{1/p}$$

(20)

where T is the length of the vector of features $\mathbf{f}\left(m\right)$. A certain range of values can lead to good initialization and results for the clustering algorithm. In Section 3 we discuss about these values.

When the values of ${∥\mathbf{f}\left(m\right)∥}_p$ are sorted in ascending order the result is not a linear function of m, but rather a power function. For an automated initialization of the clustering algorithm, the range of values of ${∥\mathbf{f}\left(m\right)∥}_p^q$ are divided into 5 intervals with the same length, which provides a good starting point. Therefore, the kth interval is $I_k=\left[E_k^{\left(i\right)},E_k^{\left(s\right)}\right]$ where

$$\begin{array}{ccc}\hfill E_k^{\left(i\right)}& =& \hfill \frac{k-1}{5}\left(max{∥\mathbf{f}\left(m\right)∥}_p^q-min{∥\mathbf{f}\left(m\right)∥}_p^q\right)+min{∥\mathbf{f}\left(m\right)∥}_p^q\hfill \end{array}$$

(21)

$$\begin{array}{ccc}\hfill E_k^{\left(s\right)}& =& \hfill \frac{k}{5}\left(max{∥\mathbf{f}\left(m\right)∥}_p^q-min{∥\mathbf{f}\left(m\right)∥}_p^q\right)+min{∥\mathbf{f}\left(m\right)∥}_p^q.\hfill \end{array}$$

(22)

Although the value of q could change for different sets of features, square root is a good selection for the four proposed. We followed the criterion of fitting the first and second half of the sorted values of ${∥\mathbf{f}\left(m\right)∥}_p$ to a power function and set q as the average of both resulting powers. For the selected 4 features, the value is close to 2. In Section 3 we discuss about the possible values of this parameter. The initial clusters are assigned based on the norm of the features vector; if this norm lies in the kth interval, then the corresponding vector belongs to the kth cluster. So the initial kth cluster, ${\mathsf{\Omega }}_k^{\left(0\right)}$, is set as

$${\mathsf{\Omega }}_k^{\left(0\right)}=\left\{m:{∥\mathbf{f}\left(m\right)∥}_p^q\in I_k\right\}\forall k=1,\dots ,5.$$

(23)

The initial kth centroid is set as the mean of the features vectors that belong to the kth cluster

$${\mathbf{c}}_k=\sum _{m\in {\mathsf{\Omega }}_k}\mathbf{f}\left(m\right)\forall k=1,\dots ,5.$$

(24)

These initial centroids were used for k-means algorithm equiped with a measure of distance based on p-norm. Only 3 iterations of the algorithm are needed to achieve the scores of the images.

3. Results and discussion

Our dataset was the same as in [14] and consisted of 94 images taken from stained xenograph tumors, using HCT116 cells implanted in mice treated with either vehicle or with a cell-permeable $\alpha$-ketoglutarate derivative, dimethyl$\alpha$-KG (DMKG), followed by treatment with the mTORC1 inhibitor temsirolimus (TEM) or with metformin (MET). Samples were processed as described in [25]. All procedures were approved by the corresponding institutional organizations (APAFIS# 10090 2017052409402562 v2). Omission of the primary antibody in the immunostaining procedure was used as a negative control. Images were acquired in TIFF format with a Leica DM6000B microscope using $\times 20$ or $\times 40$ objective lenses and a Leica DFC500 digital camera. All images were independently scored from “1+” to “5+” by four expert observers. The agreement among observers measured as the percentage of images that have been annotated with the same score by all experts is $70.21\%$. This makes that the maximum mean agreement of a scoring with all observers be $90.96\%$.

3.1. Results for stain separation step

The stain separation procedure achieves satisfactory results, providing a good decomposition of DAB and H planes, showing consistent solutions along the whole set of images. In Figure 4 the estimated color stain vectors for each image, together with the mean of these color vectors, are depicted. Note that the direction of the vectors is very similar for all the images. In Figure 5 and Figure 6 the results for the separation procedure are show for two example images from the set. In the first case the processed image has a high concentration of both (DAB and H) stains, whereas the second image has a very low concentration of DAB satin. In both figures two results are shown: that for individual analysis without considering the rest of images, based on the estimate from (8); and the estimation of stain concentrations from the average stain vectors estimated in (11). In all cases the color of H stain is correctly estimated. Furthermore, even if the color of DAB stain for the second image is not very accurate at the first analysis, the calculated concentrations of this stain are similar to that obtained from the average stain vectors.

3.2. Performance of the scores prediction

In [14] we showed that 1-norm of the DAB concentration stain correlates better than average threshold method (ATM) score and pixel-wise H-score to the scorings of the observers. There we used two other features (related to that denoted as $f_3\left(m\right)$ in the present work, but not the same) that improved the correlation. Here we show how the correlation for 1-norm increases when we add the features $f_2\left(m\right)$, $f_3\left(m\right)$ and $f_4\left(m\right)$. All the features that we propose increase (in general) when the score increases, which allows us to use the norm of the feature vector for an automated initialization of the clustering algorithm, due to the high correlation between this norm and the score of each image, as it can be seen in Figure 7.

Furthermore, the use of only the first three features allows to achieve excellent results by using the above mentioned automated procedure for the initialization. Figure 8 shows a comparison of the prediction scores and the median of the observers’ scores for each triple of feature values. Despite small differences the automated scoring achieves a high level of similarity with the reference truth.

Although the procedure was conceived to use four features, we also tested its behavior when only two or three features are exploited. Table 1 shows the main results for 2, 3 and 4 features in terms of percentage of coincidences with each observer and with the median of the observers’ scores. Features $f_1$ and $f_2$ were the couple used for testing our method with 2 features. After, we added $f_3$ and, finally, $f_4$. A high level of coincidence is achieved for 3 and 4 features, whereas the level of coincidences for two features was acceptable.

A wide range of values are suitable for the parameters p (p-norm for initialization and distance in clustering algorithm) and q (the power for reshaping the curve of values of the norm of the features vector). Figure 9 shows the percentage of coincidences of the scoring assigned by the proposed algorithm with the median of the scores assigned by the observers. For values of more than $90\%$ low correlation between p and q can be observed. This leads to consider the proposed method quite robust, in terms of the selection of the parameters, which enables an automated procedure for the scoring of the IHC images.

4. Conclusions

In recent years, automated scoring methods have been developed for the evaluation of IHC images. However, these methods suffer from problems related to the need for researchers to intervene at a certain point in the process, the use of control or training images or the dependence of the results on the choice of certain parameters.

In the present work we have proposed a completely unsupervised method that allows scoring IHC images without the need for reference images and without dependence on the set up of several parameters. The procedure consists of a first stage of separation of the DAB and H stains. In this stage, PCA allows working in a lower dimensional space (a plane), which allows finding the vectors associated with both stains in a simple and efficient way. Once these vectors and the matrix of intensities of both stains are obtained, four characteristics are used to cluster the images, although it has been proven that the use of three of them provides very good results and, even, the use of only two achieves acceptable results. The clustering method used (k-means) is based on automated initialization, without the need for reference images or parameter setup, making it a robust unsupervised method.

The proposed method has been tested with a database of 94 images and the reference truth was taken from four expert observers. The scoring was made up of five categories and the algorithm reached a mean percentage of correct predictions of 87.23% for a maximum of 90.96% (due to the experts’ opinion percentage of coincidences, which limited the average maximum percentage of correct predictions). Adapting the method to a different number of categories would be immediate, as in our previous work.

As a future line of work, we consider it interesting to find simple, fast and effective methods to locate nuclei and to develop new features related to the DAB stain of nuclei, since this could help to improve the results.