Digital Measurement of Myelofibrosis Associated Platelet Derived Growth Factor Receptor [Beta] (PDGFR [Beta]) Expression in Bone Marrow Biopsies

Journal of Computing and Information Technology - CIT 21, 2013, 1, 47–56 doi:10.2498/cit.1002109 47 Digital Measurement of Myelofibrosis Associated Platelet Derived Growth Factor Receptor ββ (PDGFR ββ ) Expression in Bone Marrow Biopsies Szilvia Szeghalmy1, Judit Bedekovics2,G ´ abor M´ ehes2and Attila Fazekas1 1Department of Computer Graphics and Image Processing, University of Debrecen, Hungary 2Department of Pathology, University of Debrecen, Hungary In daily routine the reticulin silver staining is used on bone marrow biopsy samples as a gold standard for the characterization of myelofibrosis, however this method does not provide information about the prefibrotic stage. Recently a specific immunohistochemical method was introduced which may overcome these weaknesses of reticulin staining. Activated fibroblasts responsible for stromal proliferation are highlighted by increased PDGFR β expression, which can be presented by immunohistochemistry in bone marrow samples. Using this staining the pre-fibrotic stage can become detectable and we have information about the disease activity. During development of new staining method it is important to prove its reliability and usability. In this paper we introduce a digital image processing method to measure paranchymal damage in digitalized histological slides that can aid correct interpretation of the staining. Keywords: myelofibrosis, PDGFR, image processing 1. Introduction Myelofibrosis (MF)is a lingering disease which replaces normal cells in the bone marrow to fibrotic tissue containing variable degrees of reticulin and, in the advanced phase, collagen fibres. The severity of the disease is usually characterized by the degree of fiber content highlighted by the reticulin silver staining (G¨ om¨ ori’s staining)in bone marrow biopsy samples. Recently, different grading systems have been introduced to assess bone marrow fibrosis, most of them derived from the Bauermeister scoring system [2]. A few years ago, Thiele et al. reached a consensus defining a grading system [17]. This scheme consists of qualitative and quantitative analysis of bone marrow fibrosis and distinguishes four increasing categories, ranging from MF0 to MF3. This wide-spread grading system allows making more accurate prognosis [6], describes precisely the fiber content of the bone marrow, but does not provide information about fiber producing cells which are called fibroblasts. Lately we identified PDGFR β as a novel biomarker of actived fibroblasts. PDGFR β controls stromal proliferations and turned out to be selective for fibroblasts in the bone marrow. The amount of PDGFR β is strongly dependent on their number and activity and in this context it could be applied as a measure of the fibrotic process. A new scoring system based on PDGFR β expression was also introduced. This system resulted in excellent agreement with the classical silver staining method, although different aspects of the same process are analysed (Table 1, Figure 1). As a special feature, the new method also covers the earliest, prefibrotic phase of the disease due to the presentation of all fibroblasts potentially participating in reticulin fiber synthesis. The method is simple and requires only a standard immunohistochemical staining for PDGFR β by one of the commercially available antibodies (anti-PDGFR β clone ab-32570, Abcam was used in this study). 48 Digital Measurement of Myelofibrosis Associated Platelet Derived Growth Factor Receptor β (PDGFR β )... PDGFR β immunopositivity can be assessed in a semi-quantitative way during visual inspection in the microscope. However, a much better and reproducible assay would be required for the comparison of different clinical samples or samples from different time points from the same patient. In order to perform measurements on PDGFR β related immunopositivity in bone marrow samples, digital image analysis was used and special algorithms were applied to detect the parenchyma (useful region)and the positive component. Grade Description 0No significant amounts of positive fibroblasts. 1 Isolated positive fibroblasts and/or their branches clearly recognizable in the intercellular space. 2 Many positive fibroblasts with long processes are present accompanied by intersections forming a loose network. 3 Masses of positive fibroblasts and processes are present accompanied by frequent intersections and bands forming a dense network. Table 1. Grading system based on PDGFR β expression of activated fibroblasts. Figure 1. Grading system based on PDGFR β expression of activated fibroblasts. 1.1. Related works Although manual methods are the gold standards for analysis of immunohistological staining, efforts are being made to apply automatic digital analysis as well. Most of the solutions are semi-automatic: manual segmentation is performed by any image processing tool (ImageJ, PhotoShop, Tmaj, etc.), and only the measurement is automatic [5]. There are some fully automatic characterisation algorithms for similar staining method to ours, for example the authors in [8]analyse lung breast cancer in mice. But the bone marrow requires specific analysing method because: •Not only the tissue to be examined can be found in the sample. •The disease replaces normal cells to fibrotic tissue. •Fibrotic tissue can deform normal cells. Our samples can include some kinds of haematopoietic bone marrow cells, fat tissue, bone trabecules, blood clots (Figure 2), muscular tissue, connective tissue (Figure 9)and noises. Figure 2. There are bone marrow nuclei (a), positive fibroblast (b), plasma of bone marrow cells (c)fat (d) and deformed nuclei (e)inside the black curve. (f)is a part of bones, (g)is a blood clot (outside the curve)and (h)shows a staining failure part. The classification of different, sometimes overlapping and crinkled tissues is not an easy task. In [12]the authors note, the traditional automatic classification algorithms, such as Kmeans and Fuzzy-C-means do not give correct result. They finally developed a semi-automatic analyser using Generalised Regression Neural Digital Measurement of Myelofibrosis Associated Platelet Derived Growth Factor Receptor β (PDGFR β )... 49 Networks. The user can set the applied features (Fourier transform of region of common GLCM texture features)and some other parameters. They did not mention any staining procedure, and there is no information whether pathological cases were investigated or not. Anyway, fibers do not appear as a separate class. In [7]authors proposed a method to segment bone trabeculae from bone marrow biopsy dye with haematology and eosyn. First, they use Watershed Transform and compute features of each segment. After that, they classify the features to four groups by K-means. Using bone regions as markers, they perform the Watershed Transform again to find the exact borders of the bone trabeculae. It is very important for us that the procedure can distinguish between the positive fibers and normal cells, thus we developed an own method taking advantage of the properties of immunohistochemical staining. 1.2. Materials and methods For research 41 bone marrow biopsy samples were available. Each sample was processed by immunohistochemical staining. In immunehistochemistry every cell nucleus is displayed by haematoxylin counter-staining which gives the nuclei a blue color. That specific protein which is in the focus of interest is bounded by a specific antibody, in this case this was a primary antibody against PDGFR β subunit. Then the primary antibody is bounded by a secondary antibody which is linked to an enzyme. If we add a substrate to the sample, a chemical interaction will proceed between the enzyme and its substrate which is accompanied by a brown precipitation formation. In other words, under the microscope the brown color indicates the localization of the analysed protein within the tissue. The score of PDGFR β was assigned to stained samples manually. The samples were obtained at 40x magnification using Light microscopy (Leica DM2500 microscope, DFC 420 camera and Leica Application Suite V3 software, Wetzlar, Germany). The Mirax SlideAC SDK was used to read the slides. The size of the whole slides was 256216 ×78336 pixels, that is why we performed the analysis with 1 : 4 magnification, tile by tile. Each tile size was 512 ×512 pixels. 2. Digital Analysis In this section we introduce our own method. Figure 3 presents the main steps. In this paper, the words nucleus and cytoplasm will denote the nucleus and cytoplasm which belong to the haematopoietic parenchyma (the bone marrow cells, except cells of bone, fat, and connective tissue).Theuseful region is the nuclei and the cytoplasms around them. Figure 4 shows two ideal parts of samples. Normal nuclei are purple, positive fibroblasts are dark brown. The background is clear white, there are no inappropriate things, only nuclei, cytoplasm and white, round fat cells. Figure 3. Main steps of the algorithm. 50 Digital Measurement of Myelofibrosis Associated Platelet Derived Growth Factor Receptor β (PDGFR β )... Figure 4. Ideal parts of samples. The useful region is inside the black lines. On the other hand, ideal case is a curiosity. The examination procedure implies many mistakes, since each step can cause several errors from the biopsy to the digitalization. That is why, manual selection of some appropriate region for the digital analysis is a standard. In our opinion, this step should be omitted because the samples can be inhomogeneous, so the selection may influence the result. 2.1. Creating mask of positive component The ratio of the positive fibroblasts and the useful region is one of the most important information to measure the state of disease. As stated above, the positive fibroblast is dark brown, so first, we need to define the brown mask, Bm=1,if H≤60 or 340 ≤H, 0,otherwise,(1) where Hdenotes hue channel of HSV tile. It is more difficult to find threshold level, which can separate fibroblast from the plasma of other cells and the background. It can also be important information whether the positive fibroblasts are forming a dense network, and if the threshold is high, the network will break up. And when the threshold is low, the detector frequently finds false networks because the discoloured background can connect separated parts. Using saturation channel seems to be suitable for partition, although the shapes of histogram are various. To detect interesting parts the B=1,if Bm=1andS>t, 0,otherwise, formula is used, where tis the level of threshold computed on the Sb=S,if Bm=1, 0,otherwise, where Sdenotes saturation and Bmis defined by (1). We tested some well-known automatic threshold algorithms to find suitable threshold value, such as Otsu (1997), three-level Otsu [14],Yen [18], Kittler [11]. We observed, that each tested method is sensitive to the discoloured background and plasma, except the Ramesh [15] and Arora [1]with adequate parameters that performed well also in the problematic cases. The latter procedure with the parameters, which worked well in noisy pictures, was too strict in cases of cleaner pictures, thus we chose the Ramesh method. The result of an ideal case can be seen in Figure 5b and Figure 5c – 5e gives example that discoloured plasma can mislead other threshold methods. (Only the results using Ramesh and the two Otsu thresholds are presented here.) Figure 5. Detected fibroblasts using various threshold methods. (a)The original image, (b)Ramesh method on (a),(c)Original image with discoloured plasma, (d) Otsu method on (c),(e)Higher threshold of three-level Otsu method on (c),(f)Ramesh method on (c). 2.2. Creating mask of useful region To detect useful region, the starting-point is the ideal case, when almost all the brown and purple pixels belong to it. After that, we eliminate the unwanted objects. Finally, we extend the mask of nuclei inside the remaining part. Digital Measurement of Myelofibrosis Associated Platelet Derived Growth Factor Receptor β (PDGFR β )... 51 2.2.1. Detecting useful candidate pixels Since the plasma belongs to the useful region the criteria to highlight useful candidate points are more permissive than when the aim is to detect positive fibroblasts. The mask of useful candidate from brown parts is determined by Bc and of that from purple parts is determined by Pc. Bc=1,if Bm=1andS>20, 0,otherwise,(2) Pc=1,if 160 ≤Hand H≤320 and S>20, 0,otherwise, (3) where His the hue, Sis the saturation channel of the HSV tile. The threshold of the saturation ensures that the white fat tissue is not included in the mask. 2.2.2. Detecting lobe-like objects In an ideal case the Bcand Pcform together the useful region, but the pixels of other tissue also can satisfy the formula Pcor Bc.The non-useful parts in purple region, except the folds, are commonly smoother than the useful segments. Whereas the middle of nucleus is bright coloured surrounded by pale cytoplasm or brown fibers, the lobe-like objects are homogeneous. Figure 6 displays typical kinds of non-useful regions. Figure 6. Typical kinds of uninteresting parts. To detect lobes, the homogeneous region was highlighted from the purple region by L0=1,if 160 ≤Hand H≤320 and Eg+Es<40, 0,otherwise, formula, where Eg,Esare the edge images computed by Sobel operator with 5 ×5maskonthe gray-level respectively saturation image. The formula causes the segments containing group of nucleus to fall apart, while the lobe-like objects remain whole. (Figure 7b) After connected component labelling [4]was used on L0to detect each blob, we copied all the blobs larger than the given limit to L.(In our case the limit was set to 500.)Since the edge of blobs does not belong to the L, we dilate L by the following mask: M= ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 0011100 0111110 1111111 1111111 1111111 0111110 0011100 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ (4) In the next step, we used hole filling on Lto eliminate larger hole in the blobs. Figure 7c presents the result image. Finally, we remove L from Pc.LetPdenote the new image. Figure 7. Steps of eliminating the lobes. (a)Original image, (b)Raw lobe mask (L0),(c)Detected lobe (L), (d)Purple nuclei (P). 2.2.3. Extending mask of nuclei Remember that Bccontains only positive fibroblasts, but we also need the plasma between them. Since the interesting plasma regions are very similar to lots of non-useful parts, the region grow algorithms cannot be used. Instead of them, we fill the small holes between the nuclei inside the candidate region, but those pixels which were removed from Pc, cannot become useful pixels. U0=Por B, U=(U0•M)and (Por Bc), 52 Digital Measurement of Myelofibrosis Associated Platelet Derived Growth Factor Receptor β (PDGFR β )... where •is the morphological close operator and Mis the Formula 4. Figure 8a presents result of this step on Figure 7a. Figure 8c is the detected fibroblast and the normal cells (U0)in Figure 8a. We extend this region applying M, which can eliminate the small holes between the fibers and nuclei (Figure 8d). In good case the extended part is plasma, but background or lobe-like object may be added to the result as well. So we keep only pixels falling in the useful candidate region (without lobes)(Figure 8e). This is ensured by the “and” operator between the two parts of the formula. We can observe in this example that (Por B)almost perfectly matches the useful region. It’s a bit misleading, because this part can contain also the discoloured background, if the stain has not been perfect. Figure 8. (a)Detected useful region on 7.a. (b)-(f)The detailed example of computing U:(b)an original image, (c)positive fibroblasts and normal cells (U0), (d)U0dilated by M,(e)useful candidate region without lobes (Por Bc),(f)the result. 2.3. Safety filter However, most of the noise was removed by the previous steps, some of uninteresting parts could remain in the U. It can happen that the useful region in the slide is much more smaller than the other (foreground)regions. Accordingly, we have to deal with this problem, else lots of small errors could influence the results. One reason for the error is if another sort of tissue gets into the sample (See Figure 9), or folds remain in the detected useful region. Figure 9. Detected bone marrow cells in muscle tissue (a)and in connective tissue (b). When the mask of useful region was created, most of these errors were filtered out. The greater the difference between the Pcand P, the more likely that the remaining parts of purple one are bugs. Based on our experience the false remaining parts are less than 15% of the removed part. The left side of Formula 6 estimates the ratio of the remaining lobes area and the nuclei area. We also have to check whether the tile contains enough data (Formula 5), because we will introduce features to characterize the samples in the next section, and some of them may give false result if useful area can be hardly found in the tile. a(U)>s 10 (5) 0.15 ·(a(Pc)−a(P)) a(B)+a(P)<e,(6) where a(X)the area of Xand sis the size of the tile, and eis the given limit of error. When features of the whole slide are computed, only tiles satisfying (5)and (6)are used. On average, 48% of the tiles (contains data)satisfies the condition. In many cases almost whole sample is bone, fat and connective tissue. This is consistent with the experience of other researchers [12]. We have processed 25870 Mirax tiles, of which 6819 contained any data (others only contained background). Table 2 presents the detailed result. Rejected by All Contains inappropriate obj. Formula 5. 2488 1681 Formula 6. 884 883 Table 2. Number of tiles filtered out by Formula 5 and Formula 6 with e=0.01. Digital Measurement of Myelofibrosis Associated Platelet Derived Growth Factor Receptor β (PDGFR β )... 53 To find out how efficiently we filter out the irrelevant parts, we implemented [7]and tested it on our samples (Figure 10). It is obvious that the comparison cannot be complete because the staining methods are different. Moreover, the goal of the whole procedure in [7]is to find bones, while we remove lobe-like objects (including the bone)in various steps. Considering only the bone detection, it is varying which algorithm is more efficient. Gonzalez’s algorithm uses fix number of clusters, thus if there is no bone in the picture, this algorithm inevitably fails. In general, the folds can also cause an error because K-means often classify the pixels of the simple and the double layer of bones to different classes (Figure 10b).Our lobes detector is also not perfect. If the lobes are discoloured to brown a bit, our algorithm does not find them correctly (Figure 10g), but other steps can eliminate these parts from the sample (Figure 10h)and the safety-filter can throw out the tile. 3. Characterisation of Samples After creating masks, some parameters were computed to characterise the samples. The ratio of positive fibroblasts area and useful region area is commonly used in reticulin silver analysing [17]. Based on our results, this ratio is not always appropriate to distinguish the stage of samples, because the area of thin loose network (PDGFR β 2)can be less than the area of isolated positive fibroblast (PDGFR β 1). Thus, other parameters are also calculated. Table 3 includes calculated parameters. Parameters Description PFibrA Area of positive fibroblasts regions. FtbCount Number of the positive fibroblast blobs (pfb). SumPerm Sum of the perimeter of pfb. WPerm Weighted sum of the perimeter of pfb. Weight is the number of cross and end points of the skeletons of pfb. Top50A Sum of the area of the 50 biggest pfb on whole slide. Top50P Sum of the perimeter of the 50 most length perimeter of pfb. Top50S Sum of the length of skeleton the 50 most length perimeter of pfb. Table 3. Computed features. The PFibrA,FtbCount, SumPerm and WPerm were normalised by the area of the useful region. Figure 10. (a),(e)Original images with annotation. (The bone trabecula is inside the black lines.). (b)Result of Gonzalez method K-means step. (c)Result of Gonzalez method on (a). (d)Detected lobe by our method on (a).(f)Results of Gonzalez method on (e),(g) Detected lobes (L)by our method on (e).(h)Detected purple nuclei (P)on (e). 54 Digital Measurement of Myelofibrosis Associated Platelet Derived Growth Factor Receptor β (PDGFR β )... Out of the 41 scanned samples 11 are PDGFR β 0, 9 are PDGFR β 1, 10 are PDGFR β 2, 11 are PDGFR β 3 by manual grading. Our observation is that, most parameters are directly proportional to the manual grades, but there are some exceptions. For example, the FtbCount is often larger in PDGFR β 2thanPDGFR β 3, hence when the pfbs connect together, the number of pfbs is decreased. SumPerm also can be almost equal in both cases, because when the pfbs close to each other connect, the perimeter of the new object can be less than the sum of the perimeters of the original pfbs. The weight in the WPerm can compensate that phenomenon. Figure 11 presents the WPerm parameters. The “Top” parameters are mainly introduced to detect emergence of positive fibroblast networks. Figure 11. Relationship of WPerm and PDGFR β grade. At present, only the FtbCount is normal distributed based on Shapiro-Wilk Test with 0.05 level, thus we chose the Kruskal-Wallis test to compare means of groups. The H0(the median of variable are the same across categories of PDGFR β )was rejected for each variable. The MANOVA analysis confirmed this observation. Figure 12 presents the CVA scatter plot [9].Of course, this result must be treated with caution, because of the sample size. Figure 12. CVA scatter plot: green, blue, purple and red denote the a 0, 1, 2, 3 grades in this order. 3.1. Reliability We defined some other features that can indicate the quality of sample or reliability of results (Table 4). Some cases contain only a few usable tiles, which is relevant information about the sample, since this value is between 50 and 200 in general. After using the safety filter, most samples are homogeneous by PFibrA.Occasionally, the positive fibroblast is stronger in a few tiles, that can be arisen from detection error, but it also can be true. The stain error is more typical, and we cannot filter out this type of error well because the wrong segment can be easily mistaken for healthy one. The texture is the same, only the error part is smoother a bit. Figure 13 is a scenic example because PDGFR β grade of this sample is 3. Figure 13. Thumbnail of whole slide. Only the middle of the sample was stained well. We create the histogram PFibrA with 0.1 bin. If everything goes well, the histogram is unimodal. Otherwise one of the errors (mentioned above)occurs. We determine the Otsu threshold, and compute how many tiles are smaller (NumLow)or larger (NumUsable-NumLow)than the threshold. Our experience shows that if there are only few outliers, these tiles contain big vessel or the tissue is discoloured near the bone. If there are almost equal tiles in both sides of threshold, the sample was probably badly stained. In this case the program will not compute the result automatically. First, an expert decides which group of tiles is suitable for calculating the features. Digital Measurement of Myelofibrosis Associated Platelet Derived Growth Factor Receptor β (PDGFR β )... 55 Parameters Description NumUsable Number of tiles satisfy (5)and (6). NumFiltered Number of tiles satisfy (5)but do not statisfy (6). Range Range of histogram. HistoType Is histogram unimodal or not. NumLow Number of tiles with lower PFibrA than threshold. Table 4. Reliability features. 3.2. Conclusion and future plans Immunohistochemichal assessment of PDGFR β on bone marrow biopsy samples is a typical example of semi-quantitative scoring systems which are frequently used in pathological practice. Biological samples usually show heterogeneity, the features of stage can be mixed, so it can be difficult to find the grade which is the best to describe whole sample. We have been developing an image processing method that enables us to analyse the bone marrow samples prepared by immunohistochemical staining and make a possibility of objective comparison of the samples. At present the method works well on ideal samples, and can eliminate from the samples the inordinate objects such as bone trabecules, connective tissue, muscle tissue fairly well. Our method can also indicate certain staining failure. The incipient result shows the computed parameters enable a clear distinction of early (MF0-1) from advanced (MF2-3)myelofibrosis and provide a good basis for future classification of the disease. However, because of the small sample size, further investigations are necessary. The research is going on. Now, we have been dealing with some natural phenomena caused by the staining method. For example positivity may appear besides of bones edge in healthy samples as well, and the inside of blood vessels are discoloured typically. In the immediate future, further samples will be available, that will give us the opportunity to test our analysing method on more varied samples and select the best features to determine the state of disease. References [1]S. ARORA,J.ACHARYA,A.VERMA,PRASANTA K. PANIGRAHI, Multilevel thresholding for image segmentation through a fast statistical recursive algorithm. Pattern Recognition Letters,29(2)(2008), 119–125. [2]D. E. BAUERMEISTER, Quantitation of bone marrow reticulin – a normal range. Am. J. Clin. Pathol., 56(1)(1971), 24–31. [3]J. BEDEKOVICS,SZ.SZEGHALMY,A.FAZEKAS,G. M´ EHES, Platelet derived growth factor receptor – an early marker of myelofibrosis. In Proc. of Modern sejtanalitikai m´ odszerek: a VII. Magyar Sejtanalitikai Konferencia, 31 May–2 September,(2012) Budapest, Hungary, 160., ISBN 978-963-9070-899(in Hungarian) [4]F. CHANG,C-J.CHEN,C-J.LU, A Linear-Time Component-Labeling Algorithm Using Contour Tracing Technique. Computer Vision and Image Understanding,93(2)(2004), 206–220. [5]T. C. CORNISH,A.M.DEMARZO, Tissue Microarrays in Cancer Research. Modern Molecular Biology Applied Bioinformatics and Biostatistics in Cancer Research,(2010), pp. 157–184 [6]U. GIANELLI ET AL., The European consensus on grading of bone marrow fibrosis allows a better prognostication of patients with primary myelofibrosis. Modern Pathology,25(2012), 1193–1202 [7]M. A. GONZALEZ,T.R.CUADRADO,V.L.BALLARIN, Comparing marker definition algorithms for Watershed segmentation in microscopy images. Journal of Computer Science & Technology,8(3) (2008), 151–157. [8]T. GOUDAS,I.MAGLOGIANNIS,Advanced Cancer Cell Characterization and Quantification of Microscopy Images.LNAI(7297), Springer-Verlag, Berlin, Heidelberg, (2012), 315-322, ISBN: 978-3642-30447-7 [9]O. HAMMER,D.A.T.HARPER,P.D.RYAN ,PAST: Paleontological statistics software package for education and data analysis. Palaeontologia Electronica,4(1)(2001). [10]R. HOFFMAN,D.RONDELLI, Biology and treatment of primary myelofibrosis. American Society of Hematology, 2007(1)(2007), 346–354. [11]J. KITTLER,J.ILLINGWORTH, Minimum error thresholding. Pattern Recognition,19(1986), 41–47. [12]G. MESCHINO,E.G.MOLNER,L.I.PASSONI,Semiautomated Segmentation of Bone Marrow Biopsies Images Based on Texture Features and Generalized Regression Neural Networks. In Proc of CACIC, 4-6 October 2004, Universidad Nacional de La Matanza, ISBN: 987-9495-58-6 [13]N. OTSU, A threshold selection method from graylevel histograms. IEEE Trans. Sys., Man. Cyber., 9(1)(1979), 62–66.