1.

Introduction

Recent advances in molecular medicine have provided a greater opportunity to understand the genetic basis of disease, as well as the cellular mechanisms of disease, and to select appropriate treatments with the greatest likelihood of success. One such technique for molecular diagnosis is in situ hybridization in which labeled hybridizing agents [such as deoxyribose nucleic acid (DNA), ribonucleic acid, or single stranded or double stranded DNA probes] are exposed to intact tissue sections. The probes can be labeled by direct or indirect means. When fluorescent dyes are used as labels, the technique is referred to as fluorescence in situ hybridization (FISH).^{1 2 3 4} The probe hybridizes to a defined target nucleotide sequence of DNA in the cell, and the dye fluoresces to some particular color when excited by a mercury arc lamp or argon laser (in the case of a laser scanning microscopy), so that the labeled probe can be visually detected when the probed tissue is viewed through a microscope. Each chromosome containing the target DNA sequence will produce a fluorescent signal (spot) in every cell when the specimen is illuminated with suitable excitation. FISH is an excellent method for the detection of gene copy number alterations in cancer and other diseases.

Application of FISH technology has been hampered because fluorescent spot counting is tedious, inaccurate, often highly subjective, and subject to substantial intraobserver variability. It also requires a highly trained technician to recognize the cells or tissue to be analyzed, and who can recognize and count tiny fluorescent spots accurately. Finally, at most 100 or 200 cells are typically analyzed per specimen, and in the case of gene amplification, much less than that (such as 20 per specimen).

An instrumental impediment to accurate FISH spot counting is that the probes hybridize throughout a three-dimensional tissue sample, and therefore the use of a single focal plane can result in a low estimation of the number of spots. This can happen in two ways. First, a lower signal results for spots outside the focal plane, thereby exacerbating the confusion between true spots and noise, and resulting in a low estimation of the number of true spots. Second, if one spot lies below another, they produce a single signal relative to the focal plane, and again there is low spot estimation. Although automated FISH spot counting algorithms have been developed based on a single focal plane,^{5 6 7} they have been subject to both problems. Satisfactory results have been achieved by using the morphological top-hat transform in a Bayesian context, where prior distributions are assumed for spot and noise intensities.⁸ This approach takes into account the distribution of intensities resulting from a single focal plane; however, it requires accurate prior distributions and therefore is very sensitive to image acquisition, in particular, precise protocol implementation on the part of technicians.

New technology now allows the cost-efficient acquisition of a stack of images on multiple focal planes throughout a tissue sample. Each of these provides an intensity gradient at a particular focal plane, and thereby overcomes the two low biases mentioned previously. This paper presents an algorithm for three-dimensional spot estimation from the image stack by using the morphological top-hat transform on each slice image in the stack. Owing to the greater accuracy of the spot intensity readings, good results are obtainable without a Bayesian methodology.

A drawback to previous automated spot counting algorithms is their dependence on the detection of nuclear boundaries. This can be extremely difficult (or impossible) in many kinds of tissue sections. Moreover, overlapping nuclei require segmentation. Not only is segmentation difficult for irregularly shaped nuclei, or nuclei sitting within complex backgrounds, even when nuclei are successfully segmented, there is the problem of determining to which nuclei a spot belongs. The difficulty is illustrated in Figure 1(a). The figure illustrates two-color FISH, in which two probes are labeled with different dyes that fluoresce with different colors. The red label (black) may, for example, be attached to a probe that hybridizes to a gene of interest (such as a hormone receptor gene that may be amplified in certain tumors). The green label (gray) may be attached to a probe that hybridizes to a known chromosomal locus that is not expected to vary in disease states (such as the centromere of a chromosome on which the gene of interest is found). The gene of interest could be recognized by a probe on the X chromosome labeled with spectrum orange (to provide an orange-red spot signal), and the reference probe could be labeled with spectrum green (to provide a green spot signal) for the centromere of the X chromosome. A single green signal would therefore be observed in the nuclei of the schematic representation of Figure 1(a) from male cells (which have only one X chromosome), while two green signals would be seen in female cells. Amplification of the gene of interest would be noted in certain cells in the schematic representation of Figure 1(a) in which there is an increase in the ratio of red signals to green signals.

Figure 1

Two-color FISH.

Relative to Figure 1(a) it is conventional to count the number of signals in each nucleus of a large number of cells. This approach has been adopted because the amplification or deletion of a gene occurs in large populations of cells, and significant changes in the copy numbers of genes are often only detected by examining a large number of (for example, at least 200) cells. Because the amplification has been considered to be a nuclear event, a change in the copy number of a gene with respect to each nucleus has been counted, both manually and in automated systems.

Our approach here [Figure 1(b)] avoids nuclear segmentation by determining the probe ratios without reference to the cells (or the nucleus) containing the probes. Thus, Figure 1(b) shows the FISH spots of Figure 1(a) in a region of interest [such as the microscope field of view shown in Figure 1(a) but without reference to the nuclear contours]. We calculate the ratio of test probes (red) to reference probes (green). This ratio provides sufficient information over a sufficiently large number of cells in a region of interest to be informative about the relative amplification or deletion of a target gene. In some instances, a region of interest is a microscopic field of view at low magnifications (e.g., 100–200×). An entire microscope field of view can then be used for the image capturing and analysis at 400–1000× magnification (×40–100 objectives) for FISH analysis. The thickness of the tissue sections used for FISH analysis is the same as in sections routinely used for histopathological analyses, ranging from 4–10 μ in thickness.

We note that a different algorithm for multiple-focal-plane FISH spot counting has been proposed.⁹ Besides involving different image processing tools, that algorithm requires nuclei segmentation and assumes even illumination. By using the top-hat transform as the basis for spot identification, the algorithm introduced here is fairly insensitive to substantial and irregular illumination gradients.

The algorithm is presented as a morphological algorithm for counting three-dimensional grains via multiple slice images. It belongs to the large class of morphological algorithms used for counting grains, without reference to the particular FISH application, and is based on the morphological top-hat transform.^{10 11} Other morphological grain-counting algorithms depend on different morphological transformations, including those based on the watershed transformation^{9 12 13} and those involving granulometric measurements.^{14 15}

There exist commercial software packages for single-focal-plane spot counting and, more recently, for multiple focal planes. We do not wish to comment on these because we do not have access to the details or knowledge of the kinds of images for which they give satisfactory results.

2.

Morphological Algorithm for Spot Estimation

The algorithm has the following general structure. A morphological top-hat transform is applied to each of L gray-scale slice images to yield L outputs, each possessing brightness intensity spikes jutting above a fairly flat background. Each bright spike can correspond either to a signal or to noise. Each top-hat image is thresholded to produce a stack of L binary images showing spike locations. Morphological filters are applied to the binary images to eliminate noise, and touching spots are segmented. Binary spot markers occurring as vertical neighbors in the stack are grouped into one final spot located at a particular assigned stack level. Various algorithm parameters are set according to characteristics of the physical images. These include window size for the top-hat transform, threshold levels, and sizes of filter structuring elements.

For a detailed description of the algorithm, we begin with an intensity function, ξ(x,y,z), defined on three-dimensional (3D) space. To model the situation in which ξ results from a set of concentrated 3D intensity signals, assume it can be decomposed as

Figure 2

Diagram of the algorithm.

The intensity ξ is sampled by taking slices at L values of z, thereby yielding L gray-scale slice images, f₀(x,y),f₁(x,y),…,f_L−1(x,y). These slice images are the actual input to the algorithm (block 1). From them, we calculate the max-image, f _max (x,y), formed by taking pixel maxima over the slice images (block 2). The histograms of f _max and f _mid , the midlevel slice image, are calculated and smoothed by a moving average filter to form H _max (m) and H _mid (m), respectively (block 5).

Figure 3 shows five slice images, sampled from top to bottom from a total of 16 slices, and the max-image arising from normal glands. In our nomenclature, these slices form Stack 64-TRI (red dye). As expected for normal glands, there is an equal number of androgen receptor (AR) signals and reference probe signals, since normal glands do not have chromosomal aberrations. The red/green (FITC/TRITC) ratio is therefore 1. Image size is 436×345 pixels, magnification is ×100, and the stack step size is 0.35 μm.

Figure 3

The max image and five slices of the stack.

For each slice image f_k, the top-hat transform is calculated by

Figure 4

Top-hat transform: signal, opening, top-hat, threshold.

Each top-hat image g_k needs to be thresholded to obtain a binary image whose components mark spikes in the slice image f_k. The algorithm provides two methods to find the threshold T (block 6). The first method, a variant of a standard approach, estimates the minimum between the peaks of the histogram representing the high-intensity spikes and the lower intensity background. The program calculates the point with the first maximum value of H _mid . This point gives the first peak of the graph of H _mid . We would like to use the derivative of H _mid to find the minimum between the first and second peaks of H _mid , but even with smoothing, H _mid is not sufficiently smooth. To avoid the small changes in the graph, we instead compute the differential operators

Because the spatial areas of the spikes in the slice images comprise only a small portion of slice-image area, and since the gray-levels of the nonspike region in g _max tend to be small, the mass of H _top is concentrated mainly at low values. Moreover, when H _top drops off permanently to small values, these small values result from high intensities in g _max . Hence, T can be chosen as a domain point of H _top at the place where this final falloff is commenced. Specifically, T is the point at which the derivative of H _top falls below a very small threshold.

T having been determined, each top-hat transform g_k is thresholded by T, thereby yielding in a binary image having value 1 at (x,y) if g_k(x,y)>T and value 0 if g_k(x,y)⩽T (block 7). Figure 5 shows the top-hat image from the middle slice shown in Figure 3 from Stack 64-TRI, H _top , and the resulting binarization by T. Each binary image is composed of disjoint maximally connected components. These components are segmented (identified) and labeled using a rapid procedure whose details we omit (block 8). A component with too few pixels, below a specified threshold τ, is deleted under the assumption that it either represents noise or a small tail of a spot whose main concentration is in other slices (τ may be set to 1) (block 9). For processing efficiency, the filtering is performed inside the segmentation routine. Following the deletions, each binary image is of the form

Figure 5

Binarization of slice 8.

A binary maximum image, X _max , is obtained by taking the maximum of the slice signals (block 10), the segmentation routine is applied to X _max , and the same data are stored for X _max (block 12). Figure 6 shows examples of X _max components and binary slice images beneath X _max . The numbers at the top of the figure correspond to the enumeration of the bounding rectangles for X _max . The numbers at the bottom in brackets give the numbers of pixels in the spots. In rectangle No. 76 the algorithm performs a horizontal segmentation to find two spots; in No. 85 it performs a vertical segmentation; and in No. 124 it performs a segmentation in which there is vertical and horizontal interaction. A threshold filter is applied to omit from further processing all components in X _max possessing less than a required number of pixels (block 14). Deleted components are not discarded, but stored in case one desires to redo the analysis with a lower threshold. For each bounding rectangle in X _max , we elicit from previous steps all slice signals lying spatially within the rectangle (block 11). The slice signal with the largest cardinality is called the main slice signal and others are tagged as being above or below the main slice signal for the rectangle. This ordering is used for visualization and for subsequent processing (block 13). One needs to be cognizant of the possibility that there may be more than one spot associated with a rectangle. The main slice signal may represent one of a number of spots lying in the stack beneath the rectangle. Graphical tools can be used to display a stack of images in which the main slice signal is shown, along with slices above and below the main slice.

Figure 6

Decomposition of eight spots by slices.

The next stage of the algorithm analyzes the slice signals associated with each rectangle and estimates the number of spots associated with the rectangle. If taken alone, a spot α_i produces an intensity function ξ_i=ϕ+α_i and slice signals f_i0,f_i1,…,f_i,L−1. Applying the top-hat transform to the slice signals yields top-hat signals g_i0,g_i1,…,g_i,L−1, and thresholding yields the binary slice images X_i0,X_i1,…,X_i,L−1. If, momentarily, we assume that α_i is a radial function and a relatively high threshold is chosen, then most of the binary images will be null, and a few consist of a single component about the (x, y) center of α_i. One of these components, say C_k, will include the others as subsets, and C_k can be associated with α_i. If the spots and noise functions composing the intensity function in Eq. (1) are noninterfering and the threshold can be picked to yield null binary slice images for each β_j and at least one non-null binary slice for each α_i, then the number of spots is found exactly by counting the number of maximal components resulting from the labels. This ideal situation is not typical. First, when the labels are densely packed, they will agglomerate to some degree; indeed, even if the nuclei are sparse in FISH images, it is possible for labels to be interfering. Second, the label and noise slice images are usually not sufficiently different to exactly separate them by a thresholded top-hat transform.

More generally, suppose spot α_i produces a set of r(i) non-null binary slice images, where the non-null binary slice images appear in a sequence of adjacent slices and the foreground of each binary slice image corresponding to α_i is connected. The slice images need to be grouped together to form a single spot for counting. The problem with counting spots is that the slice images from a particular spot are not necessarily separated from those of a different spot. This means that a slice signal may be composed of slice images from more than a single spot. Moreover, if one label is above another, there may not be a slice separating their slice images; indeed, owing to thresholding sensitivity, it may not be that all binary slices corresponding to a single spot are connected. Therefore to count spots, the slice signals must be segmented and then the resulting images combined in such a way as to provide estimates of the spots.

There are various ways in which spot contiguity manifests itself in the slice images. Many are handled by the algorithm. Others are eliminated by the deletion of very small components prior to formation of the slice signals. For a situation that is too complicated for the algorithm (not recognized by the algorithm), the data (all channels) in that portion of the stack images are not used. The data correspond to a single component of the max image. This situation has proven to be very rare in practice. To explain the kinds of contiguities treated by the algorithm, we refer to Figure 7. Each part of the figure shows a two-dimensional (2D) idealization of the spot situation, the components resulting from the spots, and the resulting slice signals after deletion of too-small components by the threshold τ. In each case, the topmost line represents X _max . Segmentation need only be applied locally under connected components of X _max .

Figure 7

Examples of spot view by Z axis.

In Figure 7(a), the labels are vertically contiguous, but are separated upon application of the threshold τ. Two spots are formed: the upper and lower sets of vertically contiguous slices. Note that we are implicitly defining a spot to be a collection of binary components. The situation is more complicated in Figure 7(b), where application of τ does not separate the slice signals of the two spots. This case is recognized by two conditions: (1) each slice beneath the maximum component contains at most a single slice signal; (2) running downward, the pixel-count sequence has two local maxima and a local minimum between the two maxima that has at least κ pixels less than the smaller of the two maxima, where κ is some preassigned parameter. The second condition insures that there are two spots, not a single large one. To form two spots as in the case of Figure 7(a), the local minimum between the maxima is deleted. If κ is not obtained, then segmentation will not occur and the algorithm will estimate there to be a single spot. A horizontal touching case is shown in Figure 7(c). After the application of the threshold τ, there is one large segment slice, two signal slices above the large slice, and two signal slices below the large segment slice. Upon the segmentation of the large slice, two spots are formed. Finally, Figure 7(d) shows another situation where there are two spots corresponding to one signal slice in the max-image. These can be isolated into two spots because they do not touch. These segmentation and grouping techniques are implemented in block 16. Graphical tools are available for visualization of reconstructed spots. In Figure 8, the algorithm identifies three spots in the red channel.

Figure 8

Complex particle No. 81 with 21 pieces on eight slices consists of three spots.

At this point, there are N spots, S₁,S₂,…,S_N, each composed of one or more binary slice signals. No two slice signals from separate spots contain common binary pixels. The identified spots are visually shown on the max-image (even though there may be more than one because of vertical occlusion). Each spot corresponds to an authentic spot or to noise. Prior to counting it is necessary to select the authentic spots (block 17). To filter out small spots, we define a size measure φ and a threshold ρ such that spot S is classified as authentic if and only if φ(S)⩾ρ. A straightforward approach is to define φ(S) to be the maximum of the pixel counts in the binary components composing S. This approach works fairly well. A problem with it is that narrow spots may be missed if their elongated axis lies vertically so that the maximum of the component pixel counts is not large, whereas the volume of the spot is relatively large. This problem is addressed by defining φ(S) to be the sum of the pixel counts in the components composing S. We have found that this latter approach works best. Not only is it necessary to eliminate small noise spots; huge spots are also not likely to be authentic. These are located on the max-image by utilizing the minimum bounding rectangles. For each rectangle R the algorithm takes the product of the maximum intensity within R and the area of R. If this exceeds a preset threshold, then the spot is not counted. Figure 9 shows the marked max-images for red-channel image (stack64TRI) and the green channel (stack64FIT) after filtering.

Figure 9

The maximum images of the filtered stacks.

In addition, if two spots appear the same in both the red and green channels, then we conclude that they result from fluorescent noise, since it is very unlikely that they would appear together if they were legitimate (block 15). Spots are considered the same in both channels according to a threshold υ : if the number of pixels in the symmetric difference of the spots does not exceed υ, then the spots are considered to be the same. In practice, υ has been chosen to be very small, no greater than 2. Figure 10 shows the max-image containing both red and green signals, with spots appearing in both channels being shown in yellow.

Figure 10

Green and red spots on stack64.

3.

Experimental Results

To illustrate the basics of the algorithm, we have used the red-dye stack64TRI from normal glands. The corresponding green-dye stack is stack64FIT. Here we summarize the results. Thresholding the top-hat images for stack64TRI yields 1166 slice signals across all 16 slices. Using the threshold τ=1, no slice signals are eliminated. With the spot-size threshold φ(S)=3, there are 133 spots (3 of these being hidden in the max image). A 5×5 top-hat transform has also been used for the green-dye stack stack64FIT. For this stack, 1233 slice signals have been found, and following the threshold φ(S)=3, there are 137 spots (six hidden). This gives a ratio of 133/137=0.97 prior to the final filtering of large spots and spots showing identically on both channels. There are six too-large spots for stack64TRI, ten too-large spots for stack64FIT, and 16 identical spots. This gives final counts of 111 for both stacks, a ratio of 1.

For a second example, we consider stack87, in which there is a large amplification, and in which there are ten slices. Here there are clusters in the red channel and the red (TRITC) signals greatly outnumber the green (FITC) signals. In our notation, the two stacks are stack87FIT and stack87TRI. The max images are shown in Figure 11. Following application of the top-hat transform, there are 4083 slice signals in the red channel, and this is reduced to 1732 by the threshold τ=3. After grouping the algorithm yields 469 (seven hidden) spots. In this case, setting φ(S)=3 means there is no further filtering. For the green channel, there are 366 slice signals and none are filtered out by the threshold τ=3. Grouping yields 88 spots and the threshold φ(S)=3 reduces this to 71 (three hidden). In this case, no spots are removed because they are too large or cross channels. Hence, the final ratio is 469/71=6.6. This estimate is below the manually counted ratio of 10. The problem appears to be too few slices. In fact, when the same samples are examined using 15 slices, the algorithm gives a ratio of 9.5.

Figure 11

The maximum images of the stack87.  

The results from five other stacks are given in Table 1. The columns give the specimen, FISH probes, spot counts found by the algorithm, the ratios corresponding to the algorithm spot counts, and the ratios from manual counting. The abbreviations for the FISH probes are: AR=androgen receptor, Xcen=X centromere, Her2=Her-2 gene, 17cen=17 centromere, and 17q23 is the chromosome location of the Her-2 gene.

Table 1

4.

Simulation Software

To study problems such as the number of slices required for good estimation, we have built simulation software for the model of Eq. (1) under the assumption that the background is flat. This assumption means that the top-hat transform has performed well, since the top-hat transform performs a local thresholding to eliminate background effects. In practice, the top-hat transform has worked very well, so that this assumption is warranted. The simulation is in the framework of the MATLAB-based graphical user interface (GUI). It takes a geometrically simple form of Eq. (1) in which the spots are intensity functions defined on three-dimensional balls distributed in a 3D cube of dimensions N×M×H. In effect, we are dealing with a four-dimensional image. Each ball has a brightness function α(x,y,z), with the maximum intensity at the center and the intensity falling radially to the outside. The GUI allows the visualization of the balls as well as their 2D projections onto slices, which themselves appear as ordinary images. The dynamic interface provides the manipulation of various parameters, such as the box size, number of balls, ball radii, and intensity functions. It also allows one to choose a desired number of slices of the modeled four-dimensional (4D) image, along with storage of slices and data derived from each slice. The GUI has many graphical tools, including displaying and printing data in the form of 2D and 3D images, plotting histograms, and displaying density functions.

Two distributions must be specified for the random model. A number k of balls is selected, and the ball locations are uniformly distributed in the box. An interval [1,R] is set for ball radii. The radii are beta distributed in the interval, and the chosen beta distribution B(α_r,β_r) can be displayed. The default beta distribution is the uniform distribution over the interval. Once a ball has been chosen, a radial intensity function must be defined over the ball. As presently constructed, these intensity functions are symmetric beta densities, and the maximum intensity (at the center of a ball) satisfies a beta distribution B(α_i,β_i).

Once a set of balls is generated according to the random model, a number n of desired slices is chosen, a set of equations describing model geometry is automatically solved, and the slices are displayed (and stored). The morphological spot-counting algorithm is applied to the model by applying it to the slices. The GUI allows the number of slices to be changed. Figure 12 shows the rendering of an intensity ball and a slice of the ball.

Figure 12

Random 4D ball of radius 12.

To illustrate the random model, we use 128 balls in the box 512×640×64. The radii interval is [1, 16] and the radii possess a beta distribution with α_r=1.5 and β_r=4.0. The intensity interval is [0, 255] and the intensity distribution is set with α_i=5.0 and β_i=2.5. Figure 13 shows the balls in the box, absent their intensity functions (which, as shown in Figure 12, can be seen quite well on the monitor). Using 22 slices, a total of 312 slice signals result. With the threshold τ=2, the algorithm gives the correct number of spots, 128 (11 hidden). If we use less than 22 slices, the algorithm does not perform as well for the model with the given settings. For n<22, the sampling rate (number of slices) is too low. For n⩾22, the algorithm does very well, with at most two errors for 22⩽n⩽40, the range tested. The error curve, with the percentage error as a function of the number of slices is shown in Figure 14.

Figure 13

128 random balls.

Figure 14

Error of counting of 128 balls.

Besides the effect of the number of slices, the model can be used to check many other effects on the algorithm. The model can have both spots and noise, and the distributions of these can be varied to check the effects of model parameters on separating spots from noise. The ball locations can be made nonuniform, in particular, they can be made to cluster to check the effect of the algorithm on clustered spots. The spots can be made smaller or larger, and dimmer or brighter. More peaked or flatter intensity functions can be studied. The model provides a toolbox to study the spot counting algorithm, or derivatives of the algorithm that might be developed in the future for different imaging environments.

5.

Conclusion

The algorithm proposed herein provides global spot counting in stacked three-dimensional slice FISH images without the necessity of nuclei segmentation. It has been designed to work in complex backgrounds, when there are agglomerated nuclei, and in the presence of irregular and substantial illumination gradients.

Filters are employed to remove noise and these work best when the spots are not too small. In the extreme, it is impossible to separate spots from noise if spots consist of only a single pixel. When authentic spots are large, it is easier to segment them and there is less chance they will be confused with noise. Thus, centromeric probes tend to produce more accurate results than those that attach to unique DNA sequences. Imaging at sufficient resolution and using a sufficient number of slices to discriminate the sequence probes from noise alleviates this problem.

Except for parameter setting the algorithm is fully automated. In practice, parametric settings are stable for a consistent imaging environment, and for fixed parameters the results have been robust relative to modest changes in imaging environment. The salient situation in which the algorithm is sensitive to parametric settings is when there is large amplification, say, in the neighborhood of 10–1. At that point, there are many tight clusters in which the spots are very small and often contiguous. This can result in an undercount, and the algorithm can report a ratio as low as 5–1. Still, a high amplification is detected. Moreover, the algorithm has a facility that identifies clusters and reports their number, thereby alerting the user to the cluster problem.