A High-Throughput Analysis Method to Detect Regions of Interest and Quantify Zebrafish Embryo Images

Introduction

Zebrafish is a widely used vertebrate model in neurobiological research, including neural development, gene expression and function, and drug development. ^1,2 Compared with other vertebrate models, zebrafish, especially its embryo, has unique features that make it well suited for high-throughput drug screening. For example, the zebrafish embryo is transparent to light and has a short development stage; thus, we can observe its morphogenesis over time using a microscope. For this reason, zebrafish has been widely used not only in neural developmental research ³ but also in vasculature studies. ⁴ For example, zebrafish has been explored for modeling neurodegenerative diseases such as Alzheimer disease. ^5,6

The small size of a zebrafish embryo makes it easy to be placed in a well plate for time-lapse high-throughput imaging. Phenotype-based small-molecule screening in zebrafish has been described in a number of studies. ⁷ The above features, in combination with the modern microscope, mean that one can obtain a large number of images of zebrafish embryos over a short time period. Correspondingly, challenges arise from the requirement to quickly analyze the large number of images effectively and efficiently.

Computerized analysis has been extensively explored to quantify zebrafish images such as their neurons. ^8,9 One can differentiate neurons close to each other by varying their contour level intensity as shown by Kamali et al. ¹⁰ Other methods to segment touching nuclei in zebrafish images include gradient flow tracking and gradient vector diffusion. ¹¹ An artificial intelligence–based method has also been developed to analyze the phenotypes of the zebrafish embryo. Vogt et al. ¹² adopted cognition network technology from satellite image analysis to process zebrafish embryo images acquired in well plates. In addition to analyzing cellular structure in zebrafish, researchers have also monitored the behavior of zebrafish and their embryos in neurological studies.

Affecting manual analysis and inspection are the problems of inter- and intraobserver variations that may impair experimental reproducibility. Reducing variation would come through automatic quantification of pigments in a region of interest (ROI) in a zebrafish embryo that will also improve the efficiency of processing experiment results. Such a method would become the foundation for a high-throughput screening (HTS) tool of zebrafish images. The pigmentation of different groups will reflect how the zebrafish embryos respond to various environmental cues. Although image analysis methods have been adapted to process zebrafish images in the past, automatically identifying the appropriate ROI at first and then correctly segmenting pigmentations in the ROI in a high-throughput manner has been challenging. We developed a pipeline of imaging processing techniques to automatically detect the ROIs and then quantify the pigments contained in the ROIs.

Methods

Image-processing algorithm

Preprocessing

The pipeline of computerized zebrafish image analysis is outlined in Figure 1A . It consists of 3 major steps: segmentation of the zebrafish embryo from a background, detection of the ROI, and quantitative measurement of pigments in the ROI.

OPEN IN VIEWER

Fig. 1.

(A) Pipeline of computerized zebrafish image analysis. (B) Illustration of a pyramid algorithm to reduce figure size. (C) An embryo image is divided into 4 quadrants. After rotating and flipping the embryo, its head will be quadrant 1 with its back pointing upward and abdomen facing downward.

The segmentation procedure includes downsampling, edge detection, Wiener filtering, median filtering, and morphological operations. For the purpose of quickly processing a large number of images in a short time, a pyramid algorithm ¹³ is used to downsample the original image of 2048x1536 pixels to the size of 512 × 384. Compared with bilinear algorithm, a pyramid algorithm, a series of images in which one image is derived from a low-pass filter and the downsampling of its predecessor image, can better preserve details of the image. ¹³ Given an original image g₀ with R × C pixels, after first application of low-pass filtering and downsampling, we obtain g₁, which has a smaller size and lower resolution than g₀. Repeating this process, we can obtain a series of images g₂,. . .,g _k . Here the low-pass filtering is equivalent to convolving g_k with a 2D weight function that resembles a Gaussian distribution (i.e., Gaussian pyramid) such that

g k ( i , j ) = ∑ m = − δ δ ∑ n = − δ δ w ( m , n ) g k − 1 ( 2 i + m , 2 j + n ) ,

(1)

where (i,j) denotes a pixel in g _k , k = 1,. . .,K and w is the Gaussian weight function. In this work, the δ is set to 2. For a pyramid of N levels, the row and column numbers of g_k are given by

R k = 2 − k ( R − 1 ) , C k = 2 − k ( C − 1 ) + 1 ,

(2)

respectively, where k = 0,. . . ,N − 1 ( Fig. 1B ). The reason for using the Gaussian pyramid algorithm to reduce the image size is its high-speed in computation. Because we are interested in the number of pigments of the zebrafish embryo, it is reasonable to use images of reduced size under the condition that the pyramid algorithm does not remove any valid pigments.

After being processed by the pyramid algorithm, the resulting image, g_p, undergoes a series of steps to find the outer boundary of an embryo. The Sobel operator is applied to find edges in g _p , where the embryo is considered the foreground. Due to the presence of noise, the edge map, g _E , contains clusters, broken boundaries, and the internal structure of the embryo. To improve upon g _E , we next use a low-pass filter to smooth out the artifacts. The filtered image g _F tends to have a well-connected boundary outlining the shape of the embryo that, in turn, allows better segmentation. However, g _F may still contain holes inside the embryo and some clusters outside the embryo, so we apply morphological operators of dilation followed by erosion to fill the holes and connect adjacent foreground pixels. After morphological operations, the image g _M usually contains a well-connected embryo object. In some cases, clusters located outside the embryo in g _M are then removed by median filtering on the binary image, g _Bin . Performing a pixel-wise multiplication of the binary image with g _P , we obtain the resulting gray-scale image of embryo g _Emb that is then used to detect the ROI.

Detect ROI and segment pigmentations

As the orientation of the embryo in each image varies, we next rotate the image so that the embryo’s head is located in the upper right quadrant with its back pointing upward and abdomen downward. The reason for rotating the image is to facilitate ROI detection in which we explore a priori anatomic information of the embryo so that the ROI is always located between the back and abdomen. The rotational operation first divides g _Emb into 4 regions of same size ( Fig. 1C ). Then, on the basis of knowledge that the head region, especially the eyes, of an embryo in g _Emb is much darker than its other regions, we can easily detect the quadrant that contains the largest part of the head by using a global threshold. If the head region is determined to be in quadrant 1, no rotation is performed. If it is in quadrant 2, the image will be rotated 90 degrees clockwise. Similar rotations are performed if the head region is in quadrant 3 or 4. To differentiate the back and abdomen of the embryo, we compare the curvatures of the 2 boundaries using g _Bin where the curvature is defined as the change rate of the unit length of the arc on the curve. For a pixel P at location (x,y) on a curve C = f(x,y), the curvature at the point is given by

K ( x , y ) = | f ″ | ( 1 + f ′ 2 ) 2 / 3 ,

(3)

where f′ and f″ are the first- and second-order derivative of C at point (x,y), respectively. Assuming a discrete curve has N points P ₁,. . .,P_i ,. . .,P_N , according to the formula for calculating the continuous curve, the formula for calculating the discrete curvature k_i at point P_i = (x_i ,y_i ) is given by

k i = x ′ i y ″ i − y ′ i x ″ i ( x ′ i 2 + y ′ i 2 ) 3 / 2 ,

(4)

where x_i ′ = x_i+1 − x_i-1 , x_i ′ = x_i+1 + x_i-1 − 2x _i and similarly for y′ and y″. Then the average curvature of a discrete curve in discrete format is

k ¯ = ∑ 2 N − 1 | k i | N − 2 ,

(5)

Based on the a priori information that the back of an embryo usually has a lower curvature than the abdomen, the boundary with a smaller k– is considered to be the embryo’s backside. By rotating and flipping embryos to ensure that the head locates to the upper right position and the back to the upper direction, we obtain consistent positioning of embryos in all images.

Within the embryo of g _Emb , a watershed method is used to separate the head region from the torso that contains the pigmentation. Specifically, the pigmentation is contained in a region approximately halfway between the back and abdominal boundaries of the torso.

Pigments are determined by finding the connected domain tracking enclosed in the ROI. All the connected domains in the ROI are labeled after using Otsu’s method for segmentation. ¹⁴ If they satisfy the following two conditions, the connected domains are considered valid pigments.

The size of the connected domain should be within a certain range. We set the threshold of size from 10 to 70 pixels based on the magnification of the microscope and pyramid algorithm that is used to reduce the image size.

The Euclidean distance between the geometric center of a connected domain and the central line of the embryo is within a certain threshold as the pigments of interests are located close to the central line of the embryo. In this work, the threshold is set to 16 pixels.

In addition to pigments, we calculate the curvature of the centerline, defined as the ratio of the Euclidean distance between its two endpoints over the length of the curve,

L = ∑ k = 1 n − 1 l ( k ) ≡ ∑ k = 1 n − 1 ( ( x k + 1 − x k ) 2 + ( y k + 1 − y k ) 2 ) 1 / 2

(6)

C = ( ( x n − x 1 ) 2 + ( y n − y 1 ) 2 ) 1 / 2 L ,

(7)

where x_k ,y_k ,k = 1,. . .,n are the pixels on the centerline. By definition, C is between 0 and 1, where a large C means the curve is close to a straight line.

Results and Discussion

An example of using the method to segment pigments is shown in Figure 2 . Figure 2A shows the result of a Gaussian pyramid method, in which an original image of 2048 × 1536 pixels is reduced to 512 × 384 pixels. Applying the Sobel operator on Figure 2A yields the edge map of Figure 2B , from which we observe many broken edges, holes, and clusters that will be removed in the next steps. Using the Wiener and median filter on Figure 2B to remove the artifacts creates a smoothed result ( Fig. 2C ), although the embryo area still contains some holes. Next, morphological operations of dilation and erosion fill the holes to produce a better binarized image ( Fig. 2D ). Applying the median filter on the binary image again and retaining the largest connected domain of the image gives Figure 2E , which accurately depicts the foreground of the embryo. Overlaying that with Figure 2A yields the segmented embryo of Figure 2F . From the example, we observe that preprocessing overcomes the uneven illumination and the inhomogeneous structure of the embryo to extract it from the background. After aligning the embryo to approximately the same direction, we proceed with ROI detection and pigment quantification. Figure 2G is the result of using the watershed algorithm to segment the torso and head of Figure 2F . The torso is processed by Otsu’s method to determine the largest connected domain, as seen in Figure 2H . The upper boundary of the largest connected domain of Figure 2H is detected and marked in yellow. Its linear curve-fitting result is the red line in Figure 2I . Similarly, we found the lower boundary of the largest connected domain and mark it in blue in Figure 2J so that the region contained between the yellow and blue lines is the valid ROI ( Fig. 2K ), in which the centerline is shown in green. The segmented pigments are marked in white boundaries in Figure 2L .

OPEN IN VIEWER

Fig. 2.

(A) A representative image after pyramid algorithm, g _P . (B) Edge map given by Sobel operator, g _E . (C) Low-pass filtering fills some holes, g _F . (D) Morphological operation fills most of the remaining holes, g _M . (E) After applying the median filter and retaining the largest connected domain, we obtain the binarized embryo, g _Bin . (F) Using E as a mask, we can easily segment the embryo, g _Emb . (G) Watershed method is used to remove the head region. (H) Otsu’s method is applied on the torso for segmentation. (I) The upper boundary of the region of interest (ROI) is depicted in the yellow line, and the fitted curve is in red. (J) Lower boundary is the blue line. (K) The final ROI is contained between the 2 boundaries. The centerline of the ROI is shown in green. (L) Segmented pigments are marked by white boundaries.

Additional examples are shown in Figure 3 , where in the first example, the segmentation and detection of the largest connected domain keep the embryo in the center while it removes the second embryo, Figure 3A-D . Another example is shown in Figure 3E-H , where we find that the preprocessing steps can segment the embryo with a large curvature. From Figure 3H , we note that the rotation is robust in the sense that although the image only contains the torso of an embryo, the algorithm detects the correct rotational direction so the head region points upper right for further processing.

OPEN IN VIEWER

Fig. 3.

Two examples of the preprocessing steps. (A) An image contains 2 embryos. (B) The corresponding g _Emb of A, where the smaller embryo is removed. (C) Rotated embryo of B. (D) Segmented pigments that are contained in the ROI. (E-H) Another example where only the torso of the embryo is captured.

We applied the image-processing pipeline to 18 images of zebrafish embryos that were treated with DMSO or GSI18 at various concentrations. Among the tested embryos, 7 were treated with DMSO and 11 with GSI18. The average lengths of centerlines and curvatures of the two groups are shown in Figure 4A . Although there is no statistically significant difference between the average lengths, a t-test shows that their average curvatures are statistically different (p < 0.05). Using the same data set, we performed segmentation and compared the results with manual analysis ( Table 1 ). Using manual analysis as the ground truth, we plot the average number of pigmentation of the two groups in Figure 4B . Student’s t-test shows no statistical difference between the two methods of measurements in both groups (p > 0.05).

OPEN IN VIEWER

Fig. 4.

(A) The average lengths (left) and curvatures (right) of 2 groups of embryos, one treated with DMSO and the other with GSI18. The 2 groups have approximately the same length, but their curvatures are statistically different, whereas the GSI18-treated group has more curved centerlines (*p < 0.05). (B) Segmented pigments by the image-processing algorithm and manual analysis. The automatic algorithm has a similar performance with manual analysis.

OPEN IN VIEWER

Table 1.

Comparing Segmented Pigments by the Proposed Method and Manual Analysis

				FP		FN
	Image	Automatic Method	Manual Analysis	#	%	#	%
1	0.5 % DMSO	5	5	0	0	0	0
2	0.5 % DMSO	6	8	0	0	2	25.0
3	0.5 % DMSO	5	6	0	0	1	16.7
4	0.5 % DMSO	6	5	1	16.7	0	0
5	0.1 % DMSO	9	9	2	22.2	2	22.2
6	0.1 % DMSO	10	10	0	0	0	0
7	0.1 % DMSO	11	10	1	10.0	0	0
8	5 % GSI18	12	11	3	27.3	2	18.2
9	5 % GSI18	5	8	0	0	3	37.5
10	5 % GSI18	14	16	0	0	2	12.5
11	5 % GSI18	6	5	1	20.0	0	0
12	5 % GSI18	5	6	0	0	1	16.7
13	5 % GSI18	9	7	2	28.6	0	0
14	10 % GSI18	6	5	1	20.0	0	0
15	10 % GSI18	8	10	0	0	2	20.0
16	50 % GSI18	8	10	0	0	2	20.0
17	50 % GSI18	0	0	0	0	0	0
18	50 % GSI18	0	0	0	0	0	0

FP, false positive; FN, false negative. A t-test on the number of pigments detected by the automatic method and manual analysis at the significance level of 0.05 has a p-value of 0.3313, indicating that there is no statistically significant difference between the 2 results.

In summary, as zebrafish has become a widely used vertebrate model in neuroscience, a need has arisen for improved image analysis of that model. With the fast development of the automated microscope, researchers can now set a microscope to automatically image zebrafish embryos in well plates over a short period of time. ¹⁵ However, accompanying the capability to acquire a large number of zebrafish images is a lack of corresponding image analysis tools to quickly process and quantify the images. Although efforts have been dedicated to develop algorithms to analyze features of zebrafish images, existing work assumes that a predefined ROI is available. Here we present an automatic ROI detection and analysis technique that is well suited for HTS of zebrafish morphology. The results demonstrate that the fully automatic technique can achieve satisfactory quantification results without human interference, which is a key requirement of HTS.

The pipeline of image analysis algorithms that we developed automatically searches for an ROI and segments pigments. Our testing shows that, when compared with manual analysis, image-processing algorithms can obtain good results in terms of robustness and accuracy. Moreover, automatic analysis is repeatable and more objective than manual analysis. In our work, after detecting the ROI, the algorithm found the pigments by segmenting the ROI using Otsu’s method that maximizes the interclass variation between the foreground and background. In more challenging cases in which the ROI may have an inhomogeneous background, other methods may be used to segment the pigments. However, care needs to be taken to ensure low rates of false negatives and false positives. The thresholds used in the work are selected based on the prior information about zebrafish embryo and were shown to generate satisfactory results. As an HTS method, it is important to achieve a balance between the accuracy and computational complexity of the algorithm. For example, in this work, the pyramid algorithm reduced the image size for fast computations in the following steps at the cost of performing the pyramid computation itself. As computers become more powerful with a large memory size, it may be beneficial in the future to skip the pyramid computation. Similarly, although one can apply morphological operations only to fill the holes in the edge map, it is faster to use a low-pass filter to quickly fill a majority of the holes and then use morphological dilation and erosion to finalize the results. The automated method in its current form can analyze a single embryo located in the center of the image. With appropriate modification in segmentation and object identification, the automated method can be adapted to analyze multiple embryos in an image. The images of this work were acquired from a 24-well plate. Zebrafish embryo can also be assayed in 96- or 384-well plates, and in these cases, whole-well imaging will allow us to capture the entire embryo and extract additional features such as its full length. As the number of images increases dramatically in 96- or 384-well plate imaging, the requirements will be different as the size and resolution of the image change and the features of interest may also change. In this case, corresponding modification of the algorithms may be necessary.

Like all HTS algorithms, errors are unavoidable and may in some cases require user intervention. Acquiring images of high quality in the first place is often the most cost-effective step as it reduces the chances of errors in all steps of image processing. To some extent, a large sample size will smooth out image-processing errors. Therefore, when possible, it is beneficial to use more images to achieve a higher statistical power. In statistical analysis, often the variation in the quantification results from a group of zebrafish can be indicative of outliers in image quality and/or accuracy of HTS algorithms. In other words, a large variation may prompt a researcher to examine the original image quality and/or resulting images produced from the HTS algorithm to determine the variation’s cause. Another limitation of HTS algorithms is that they can analyze a large number of images with high efficiency but cannot replace a sound design of an HTS experiment in the first place. An HTS experiment must take many factors into consideration, including but not limited to time, number and cost of the subjects, number and cost of the agents to be screened, capability of the hardware, and image analysis tools available to the researchers.