Bioinformatics Advance Access originally published online on April 19, 2005
Bioinformatics 2005 21(13):2994-3000; doi:10.1093/bioinformatics/bti455
| ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Segmentation of cDNA Microarray Spots Using Markov Random Field Modeling
Department of Biostatistics, Epidemiology, and Scientific Computing King Faisal Specialist Hospital and Research Center MBC No. 03, PO Box 3354, Riyadh 11211, Saudi Arabia
*To whom correspondence should be addressed.
 |
Abstract |
|---|
 TOP Abstract 1 INTRODUCTION2 METHOD3 RESULTS4 DISCUSSIONREFERENCES |
|---|
Motivation: Spot segmentation is a critical step in microarray gene expression data analysis. Therefore, the performance of segmentation may substantially affect the results of subsequent stages of the analysis, such as the detection of differentially expressed genes. Several methods have been developed to segment microarray spots from the surrounding background. In this study, we have proposed a new approach based on Markov random field (MRF) modeling and tested its performance on simulated and real microarray images against a widely used segmentation method based on Mann–Whitney test adopted by QuantArray software (Boston, MA). Spot addressing was performed using QuantArray. We have also devised a simulation method to generate microarray images with realistic features. Such images can be used as gold standards for the purposes of testing and comparing different segmentation methods, and optimizing segmentation parameters.
Results: Experiments on simulated and 14 actual microarray image sets show that the proposed MRF-based segmentation method can detect spot areas and estimate spot intensities with higher accuracy.
Availability: The algorithms were implemented in MatlabTM (The Mathworks, Inc., Natick, MA) environment. The codes for MRF-based segmentation and image simulation methods are available upon request.
Contact: demirkaya{at}ieee.org
|
1 INTRODUCTION |
|---|
|
TOPAbstract1 INTRODUCTION2 METHOD3 RESULTS4 DISCUSSIONREFERENCES |
|---|
Microarray technology is a powerful and efficient means for measuring relative expression level of thousands of genes simultaneously. A comprehensive review of the biological and technological aspect of the microarray technology can be found in Nguyen et al. (2002). Despite the recent advances in microarray technology, e.g. introduction of high-density oligonucleotide arrays, such as Affymetrix' GeneChip®, custom-made or spotted cDNA microarrays are still widely used today. Scientists who have narrowed down their study into a small set of genes prefer custom-made cDNA microarrays since, they (1) are less costly compared to GeneChip® technology, (2) are easy to prepare and to analyze in-house and (3) allow co-hybridization, i.e. two samples can be compared in one experiment.
The gene expression information embedded in spots is obtained by scanning the hybridized slides. This information is extracted using the techniques of image processing and analysis techniques. In this data acquisition process, the segmentation of cDNA spots is the most challenging task and may significantly impact the gene expression analysis (Ahmed et al., 2004). Image segmentation is a process that divides an image into mutually exclusive regions. Each region is homogeneous with respect to a region property, such as gray-level intensity. Microarray images, in general, are difficult to segment, due to factors, such as highly varying image contrast from experiment to experiment, high background noise and image artifacts, just to name a few. Segmentation of spots in microarray images can further be complicated by nonuniform shape and surface-intensity distribution. It is therefore a formidable task to develop a segmentation method that would demonstrate an acceptable and robust performance under varying circumstances.
Several methods have been developed to segment microarray spots and have been incorporated into commercial microarray image analysis software packages. Gradient-based spot segmentation is used in Dapple (Buhler et al., 2000). This approach may not work well for the arrays that are produced in-house. The spot images produced in our laboratories tend to have a doughnut- or crescent-like patches, which may be attributed to mainly the problems in the spot printing process, that feature high-gradient boundaries. Histogram method in QuantArray (Boston, MA) uses lower and higher range of percentile values for the background and foreground (spot). By default, 80th and 95th percentiles are used for spot intensities.
Seeded region growing method adopted by SPOT (Yang et al., 2001) relies heavily on the selection of a seed point. Although the approximate locations of the spot centers are known a priori, the selection of a seed point is critical for the performance of the method, especially when spot surfaces exhibit nonuniform intensity patches.
The fixed and adaptive circle segmentation methods are also among the used methods. The former is implemented in ScanAlyze (Eisen, 1999 http://rana.lbl.gov/manuals/ScanAlyzeDoc.pdf) while the latter is implemented in GenePix software for the Axon scanner (GenePix4000, 1999 http://www.axon.com/GN_Genomics.html#software). Another approach, which is named as adaptive method by QuantArray software for the GSI Lumonics Scanner and DeArray by Scanalytics (Fairfax, VA), computes a threshold based on Mann–Whitney (MW) test (Chen et al., 1997). Intensity-based segmentation techniques, such as k-means clustering were also tried in an attempt to classify foreground and background pixels (Nagarajan, 2003).
Previously, Markov random field (MRF) modeling was proposed by (Katzer et al., 2003) for the purpose of the gridding or addressing. However, here we address specifically the classification of the foreground and background pixels in target regions, i.e. spot segmentation. Since both gridding and spot segmentation are challenging tasks, they have been generally dealt with separately in the literature. For instance, Jung and Cho (2002) in their study addressed only the gridding issue using k-nearest neighbors graph.
In this study, we assume that the addressing has already been done, i.e. the locations of spot centers are already known, and directly deal with the spot segmentation problem. To this end, we propose a unique segmentation method that utilizes the contextual or neighborhood information, along with the intensity information based on MRF modeling of the target regions. In addition, we propose a method to simulate microarray images for the purpose of testing and comparing spot segmentation methods. We used simulated images to optimize the performance of our approach as well. The proposed method was validated on the simulated (gold standard) images. We also compared our results with that of the adaptive method since it is widely used through QuantArray software. This method is referred to as MW in the remainder of this paper.
|
2 METHOD |
|---|
|
TOPAbstract1 INTRODUCTION2 METHOD3 RESULTS4 DISCUSSIONREFERENCES |
|---|
2.1 Image acquisition and preprocessing
The 14 microarrays used in this study were prepared in the Interferons and Cytokines laboratory at the Department of Biological and Medical Research, King Faisal Specialist Hospital and Research Center, Riyadh. A detailed description of the microarray preparation can be found in Khabar et al. (2004). Our typical microarray slides contain about 6000 spots arranged in multiple sub-arrays. The hybridized cDNA microarray slides were imaged at two different wavelengths (channels), one for Cy3-labeled (green) sample and the other for Cy5-labeled (red) sample. The images were acquired at 10 µm resolutions by GenePix 4000 scanner, and saved in 16-bit TIFF format. The images were analyzed using QuantArray (Boston, MA) software. The adaptive method, with the parameters P-value = 0.01 and Max Spot Diameter = 140 microns, was chosen. The P-value refers to the significance level at which the statistical test, whether or not foreground and background pixels come from the same population, is carried out. The output of this analysis includes the location of the spot centers in addition to the mean or median spot and background intensities. We have used this location information to define a target region which is a 20 x 20 pixel square box concentric with the spot. The region was large enough to include the spot and its surrounding background.
Then, we applied our proposed segmentation method to these target regions identified by QuantArray to classify pixels as spot or background. The proposed method consists of two stages, namely initial labeling and MRF-based segmentation. In the following two sections, these stages are described in detail.
2.2 Initial labeling of target pixels
In MRF-based segmentation, an initial classification of the pixels is necessary to initialize the parameters of the intensity distributions and to converge to the desired solution faster. This initial segmentation can be performed using manual or automated means. We have tested several automated thresholding algorithms (Kittler and Illingworth, 1986; Otsu, 1979) and did not find them robust enough, especially for low contrast spots (Fig. 1). Then we tested the percentile method to do initial labeling, and found out that it was robust within a reasonable percentile range. Target regions denoted by 
G and R from Cy3 and Cy5 images, respectively, were combined to obtain the initial labeled image as follows:
 | (1) |
G and R were thresholded using the kth percentile (Pk). The thresholded binary images denoted by 
G and R were combined with the logical OR operation 
to obtain the initially labeled 
in which the background and spot pixels were assigned 1 and 2, respectively. Depending on the hybridization, spots in one channel may have better contrast than those of the other. Combining two channels in the initial segmentation stage provides additional prior information about the possible locations of actual spot pixels.
|
To find the optimal percentile for initial labeling, we tested the proposed method on a simulated image at 75th, 80th, 85th and 90th percentiles, and computed the concordance correlation (
c), which is described in the performance evaluation section, for the spot areas and mean intensities. As seen in Figure 2, c between the computed and true spot areas peaks more noticeably. It should be noted that the c between the mean spot intensities has a wider region within which it remains >0.98.
|
Throughout this study, 80th percentile [Pk(G,R) = 0.8], which appears to be the optimal percentile for both the spot area and the mean spot intensity, was used as a threshold for initial labeling of the pixels.
2.3 Markov random field modeling of target regions
Segmentation methods based on MRF modeling have been extensively studied and have found widespread use in medical and nonmedical imaging applications (Geman and Geman, 1984; Held et al., 1997; Pappas, 1992; Zhang et al., 2001). It provides a convenient way to combine both observed intensity and the spatial information under Bayesian framework.
To describe the MRF modeling, first a neighborhood system and a clique need to be defined. Let S denote an M x N lattice indexing the pixels in a given target region, s be the lattice point (or pixel) and 
s denote the neighboring pixels of s. The neighborhood system s must be symmetric; r 
s 
s r. A clique is a set of points, c C, which are all neighbors of each other; 
s r c, r s. An 8-point neighborhood system around the center pixel s was used throughout this study. Let L = {1, 2, ..., k} denote the label set. The k is 2 (background and foreground) in our case. Let random variable X denote the labeling process of S such that xs L is the value of X at pixel s. According to the Hammersley–Clifford theorem (Besag, 1974), the density of X is given by the following Gibbs density:
 | (2) |
 | (3) |
indicate the estimate of x*, then the objective is to find given y. In a Bayesian framework, p(x) can be viewed as the prior distribution for the true image x*, then the posterior probability, by Bayes's theorem, is proportional to
 | (4) |
 | (5) |
, respectively. can be obtained by computing maximum a posteriori (MAP) estimate,
 | (6) |
 | (7) |
|
2.4 Simulation of microarray images
Evaluation of segmentation methods by quantitative means is necessary to test and compare their performances. In biological applications, it is often the case that a gold standard is not readily available. To validate the proposed method we devised a simple yet effective method to generate microarray images with realistic characteristics.
We started off with an original microarray image (Cy3 or Cy5). The entire image was thresholded using an empirical threshold of 250 to identify background (BG) and foreground (FG) pixels. The threshold was selected so that individual spots would be disconnected. In this binary mask image, we also identified individual FG regions using the spot centers that were obtained in the addressing stage. We then know the location of pixels that each FG and BG region comprises. We also computed the area (total number of pixels) of each FG region. This is referred to as the true area of the FG region.
Then, to generate a simulated image based on features extracted from a real cDNA image, intensities of each FG region were drawn from an exponential distribution using the mean FG intensities of the respective spots in the original image. The remaining BG pixels were drawn from an exponential distribution whose mean intensity was determined from the original image using the entire binary mask image. Note all the BG intensities were drawn from a single distribution while FG intensities of each target region were drawn from an exponential distribution whose mean was estimated separately from the original respective spot region. The expfit and exprnd functions of MatlabTM statistics toolbox (The Mathworks, Inc., Natick, MA) were used to estimate the mean intensities and to generate random samples from exponential distributions, respectively. The spots in our microarray images sometimes exhibited doughnut-like shapes. During the simulation, we assumed that doughnut holes, identified as BG during thresholding, have the same intensity distributions as the BG. In the remainder of this paper, by ‘simulated’ we refer to the simulated image set (Cy3 and Cy5). A more versatile cDNA microarray simulation method can be found in (Balagurunathan et al., 2002), if one needs to simulate various artifacts encountered in cDNA microarray images.
2.5 Performance evaluation
The performances of the algorithms were measured based on regression analysis, coefficient of determination (r2) and the concordance correlation (c). The c between two samples A and B is defined as:
 | (8) |
and 
are the sample means (Lin, 1989).
In general, the mean spot intensities and background intensities are used to obtain a measure of abundance of gene expression (Khabar et al., 2004). All the spot intensities were background-subtracted. In comparing results, we presented the c along with the r2, as the c takes both bias and variation into account. Maximization of c was used as the criterion to optimize the percentile used during the initial labeling.
2.6 Implementation
The algorithms described above were all implemented in Matlab programming environment. The image simulation was also done using the Matlab functions from its Statistics Toolbox. The code was run on an IBM PC with 2.8 GHz Pentium 4 processor and 1.75 GB of memory. The computational cost of the proposed method is reasonable. For instance, the segmentation of 5632 spots on one array took 
50 s, whereas finding the locations, i.e. addressing, of 6513 spots (on a different array) with QuantArray took 204 s on an IBM PC with 1.5 GHz Pentium 4 processor and 256 MB memory.
|
3 RESULTS |
|---|
|
TOPAbstract1 INTRODUCTION2 METHOD3 RESULTS4 DISCUSSIONREFERENCES |
|---|
In Figure 1, examples of a low contrast (left panel) and a higher contrast (right panel) spot from an actual microarray image are shown. Along with the results of the proposed method, binary images obtained by Otsu's thresholding algorithm (Otsu, 1979) are also presented. The spot on the left panel is more challenging than the one on the right.
In the left and right panels in Figure 4, five spots from the simulated Cy3 and Cy5 images, respectively, are shown. The initial labeling process produced the spots in the 2nd row. The spots in the 3rd and 4th rows were obtained by the proposed and the MW method, respectively. QuantArray does not allow access to the segmented images. It can, however, display the segmented spots on a separate window. The image of this window was captured and target regions were cropped and resized to match the original size. The two concentric circles around spots were added by QuantArray to indicate the maximum spot size and the background region. Therefore, while comparing spots in the 3rd and 4th rows, these circles should be disregarded.
|
In Figure 5 spot areas computed by the proposed (left panel) and the MW (right panel) method are plotted against the true areas. The regression lines (solid line) and the r2 and
c values are also shown. The dashed line indicates the unity line.
|
Figure 6 shows the performances of the proposed and MW methods on the simulated images. The true mean spot intensities were compared to the mean spot intensities computed by the proposed and MW methods.
|
The comparison of both methods on an actual microarray image set is demonstrated in Figure 7. Table 1 shows the results of the comparison performed on 14 microarray image sets. A total of 85 145 spots were processed. Saturated spots, i.e. with intensity values of 65 535, were excluded. In experiment No. 3, the spots were 100 µm in diameter different from that of the other arrays, which is
140 µm.
|
|
|
4 DISCUSSION |
|---|
|
TOPAbstract1 INTRODUCTION2 METHOD3 RESULTS4 DISCUSSIONREFERENCES |
|---|
In this study, we proposed a new segmentation method to identify microarray spots. The proposed segmentation method models the FG and BG intensities with exponential distributions. The FG and BG histograms shown in Figure 3 support this type of model. The graphical goodness-of-fit test for exponential distributions (Sturart and Ord, 1987) indicated a very good fit (r = 0.97) for both BG and FG intensities. Hence, our approach is unique in that it assumes appropriate distributions for the intensities. The MW method, for instance, uses a nonparametric statistical test that ascribes no distribution to the intensities. The other novelty of our approach is that it takes the spatial information, i.e. the labels of the pixels in the neighborhood, into consideration by modeling the region labeling process with MRF. When combined together, the intensity and spatial information result in a more robust and accurate pixel classification process.
We have also devised a simple but effective microarray image simulation method to test and compare the segmentation methods, and to optimize the parameters of our approach (i.e. selection of the optimal percentile for initial segmentation). Using the simulated images as gold standard, the performances of the proposed and the MW method were thoroughly investigated.
The initial labeling impacts the performance of the MRF segmentation. This impact may differ depending on the contrast of the spot. For low contrast spots, since there is a larger overlap between the intensity distributions, one may prefer a higher percentile to avoid the undesired effect of noise. Using the true spot areas and mean intensities, we were able to identify the range of percentiles within which the performance of the proposed method was optimal. This range was much wider for mean spot intensities indicated by the flatter line close to one (solid line in Fig. 2). We used and also recommend using 80th percentile as a good compromise for a wide range of contrast values.
The proposed method appears to be estimating true spot areas more accurately than the MW method as indicated by the c values and regression parameters (Fig. 5). This advantage will most probably translate into a better identification of differentially expressed genes. The MW method appears to be overestimating i.e. over-segmenting the areas when spot sizes are small, and underestimating when large. This trend is also visible in the proposed method for small spots but not as profound as in the MW. The proposed method agrees better for a wide range of spot sizes in the middle of the entire range. The downward trend for large areas visible in the MW method does not exist in the proposed method.
The proposed method performed similarly for different values of ß parameter of the MRF model controlling the size of clustering. The r2 and c values between the two runs (for ß = 1.5 and ß = 2) were both 1.0 for both channels. The effect of the number of iteration was not significant either. Again, the r2 and c values between the two runs (3 and 7 iterations) were 0.98 and 0.99, respectively. Since increasing the number of iteration from 3 to 7 increased the computation time by only 14%, one may iterate the ICM up to 7 times.
The original and segmented, by both methods, spots shown in Figure 4 demonstrate that the proposed method behaves more conservatively in separating spots from the background as compared to the MW method. The MW method was reported to be not robust with respect to the parameter Max Spot Diameter when it is set to a value larger than the spot size (Yang et al., 2000 http://www.stat.berkeley.edu/users/terry/zarray/HTml/image.html. We have also observed the same phenomenon.
Table 1 lists the results of the comparison performed on 14 microarray image sets consisting of 85 145 spots in total. The regression parameters indicate that the proposed method consistently produced larger spot intensities. This may be due to the over-segmentation by the MW method, which means that some BG pixels were incorrectly identified as spot pixels. The c value indicated a poor agreement for the experiment No. 3 in which spots were 100 µm in diameter. The over-segmentation problem was probably more severe in this case. The poor agreement in three cases was reflected by the low (<0.8) c values.
We have also observed that the inaccurate estimation of spot locations during the addressing stage may result in inexact overlap between the actual spot and the fixed circle placed over the spot. As a result, some of the foreground region may be left outside the circle while some background pixels may be included inside.
As seen in Figure 1, histogram thresholding algorithms alone, such as Otsu's, fall short of providing satisfactory results, especially for low-contrast spots.
Our current investigation aims at studying the effect of preprocessing using nonlinear diffusion filtering techniques on the performance of the segmentation method and the background intensity. We are also interested in using the simulated annealing approach (Geman and Geman, 1984) to obtain the MAP estimate. Although this approach is computationally slower than the ICM algorithm (Besag, 1974), it is worthwhile to investigate as it may provide a more robust and accurate segmentation. The investigators will soon be developing a fully independent microarray image analysis tool by adopting one of the existing addressing algorithms.
|
Acknowledgments |
|---|
We would like to thank Dr Khalid S. Abu Khabar for providing the microarray images and his valuable comments.
Received on December 20, 2004; revised on April 1, 2005; accepted on April 14, 2005
|
REFERENCES |
|---|
|
TOPAbstract1 INTRODUCTION2 METHOD3 RESULTS4 DISCUSSIONREFERENCES |
|---|
Ahmed, A.A., et al. (2004) Microarray segmentation methods significantly influence data precision. Nucleic Acids Res., 32, e50
Balagurunathan, Y., et al. (2002) Simulation of cDNA microarrays via a parameterized random signal model. J. Biomed. Opt., 7, 507–523[CrossRef][Web of Science][Medline].
Besag, J. (1986) On the statistical analysis of dirty images. J. R. Stat. Soc. B, 48, 259–302.
Besag, J. (1974) Spatial interaction and the statistical analysis of lattice systems. J. R. Stats. Soc. B, 26, 192–236.
Buhler, J., Ideker, T., Haynor, D. (2000) Improved techniques for finding spots on cDNA microarrays. , Seattle, WA August 2000 University of Washington.
Chen, Y., et al. (1997) Ratio-based decisions and the quantitative analysis of cDNA microarray images. J. Biomed. Opt., 4, 364–374[CrossRef].
Geman, S. and Geman, D. (1984) Stochastic relaxation, Gibbs distribution, and the Bayesian restoration of images. IEEE Trans. Pattern Anal. Mach. Intell., 6, 721–741.
GenePix4000. A User's Guide, (1999) Axon Instruments Inc.
Held, K., et al. (1997) Markov random field segmentation of brain MR images. IEEE Trans. Med. Imaging, 16, 878–886[CrossRef][Web of Science][Medline].
Jung, H.Y. and Cho, H.G. (2002) An automatic block and spot indexing with k-nearest neighbors graph for microarray image analysis. Bioinformatics, 18, S141–S151[Abstract].
Katzer, M., Kummert, F., Sagerer, G. (2003a) A Markov random field for microarray gridding. Proceedings of the ACM Symposium on Applied ComputingMelbourne, FL , pp. 72–77.
Katzer, M., Kummert, F., Sagerer, G. (2003b) Methods for automatic microarray imagesegmentation. IEEE Transactions on Nanobioscience, 2, 202–214[CrossRef][Web of Science][Medline].
Khabar, K.S.A., et al. (2004) Expressed gene clusters associated with cellular sensitivity and resistance towards anti-viral and anti-proliferative actions of interferon. J. Mol. Biol., 342, 833–846[CrossRef][Web of Science][Medline].
Kittler, J. and Illingworth, J. (1986) Minimum error thresholding. Pattern Recognition, 19, 41–47.
Lin, L.I. (1989) A concordance correlation coefficient to evaluate reproducibility. Biometrics, 45, 255–268[CrossRef][Web of Science][Medline].
Nagarajan, R. (2003) Intensity-based segmentation of microarray images. IEEE Trans. Med. Imaging, 22, 882–889[CrossRef][Web of Science][Medline].
Nguyen, D.V. (2002) DNA microarray experiments: biological and technological aspects. Biometrics, 58, 701–717[CrossRef][Web of Science][Medline].
Otsu, N. (1979) A threshold selection method from gray-level histogram. IEEE Trans. Syst. Man Cybern., 8, 62–66.
Pappas, T.N. (1992) An adaptive clustering algorithm for image segmentation. IEEE Trans. on Signal Processing, 40, 901–914[CrossRef].
Sturart, A. and Ord, J.K. In Sturart, A. and Ord, J.K. (Eds.). Kendall's Advanced Theory of Statistics, (1987) , London Charles Griffin & Company Limited.
Yang, Y.H., et al. (2002) Comparison of methods for image analysis on cDNA microarray. Journal of Computational and Graphical Statistics, 11, 108–136[CrossRef].
Yang, Y.H., et al. (2001) Analysis of cDNA microarray images. Brief. Bioinformatics, 2, 341–349
Zhang, Y., et al. (2001) Segmentation of brain MR images through a hidden Markov random field model and the expectation-maximization algorithm. IEEE Trans. Med. Imaging, 20, 45–57[CrossRef][Web of Science][Medline].
CiteULike 
Connotea 
Del.icio.us What's this?
This article has been cited by other articles:
|
|
 |
|
  A. Daskalakis, D. Cavouras, P. Bougioukos, S. Kostopoulos, D. Glotsos, I. Kalatzis, G. C. Kagadis, C. Argyropoulos, and G. Nikiforidis Improving gene quantification by adjustable spot-image restoration Bioinformatics, September 1, 2007; 23(17): 2265 - 2272. [Abstract] [Full Text] [PDF]
|
|
|
|
|
|
J. Baek, Y. S. Son, and G. J. McLachlan Segmentation and intensity estimation of microarray images using a gamma-t mixture model Bioinformatics, February 15, 2007; 23(4): 458 - 465. [Abstract] [Full Text] [PDF]
|
|
|
|
|
|
A. Lehmussola, P. Ruusuvuori, and O. Yli-Harja Evaluating the performance of microarray segmentation algorithms Bioinformatics, December 1, 2006; 22(23): 2910 - 2917. [Abstract] [Full Text] [PDF]
|
|
| ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||