Hierarchical transformation-guided image registration for two-dimensional gel electrophoresis image

Abstract

Many proteomics studies have been devoted to investigate the gel-based methods for separation and fractionation of complex protein mixture extracted from tissues, cells, and other biological specimens. In order to fractionate, identify, and quantify proteins when coupled with mass spectrometric identification or other immunological tests, accurate image alignment method for possible significant protein diversity under different experimental conditions is of high demand. Although lots of protein spot matching and alignment methods have been proposed for two-dimensional gel electrophoresis image, very few methods work well because of dynamic appearance changes and structure changes for gel electrophoresis images. To address these difficulties, we propose a hierarchical transformation-guided registration method to not only align the global shape changes by affine transformation but also alleviate the local structure changes between two gel electrophoresis images by no-rigid deformation registration. The non-rigid registration method is used to obtain the local displacement with shape changes. Since such initial affine transformation can alleviate the global transformation difference between two images, it becomes efficient to use non-rigid registration to refine the registration result. The visualized and quantified experimental results show that our method is more accurate and robust to identify the variability between individual proteins.

Introduction

Two-dimensional polyacrylamide gel electrophoresis (2DE) is considered as a powerful tool for separation and fractionation of complex protein mixtures from tissues, cells, or other biological samples. Many proteomics and protein studies are dedicated to understand 2DE development for further analysis of the differentially regulated protein spots using mass spectrometry and western blotting. It provides the first step for further analysis of the differentially regulated protein spots using mass spectrometry and western blotting to disclose both physiological mechanisms and proteins associated with clinical pathologies that is used to discover biomarkers in many cases.^1–3

In order to fractionate, identify, and quantify across individual proteins, 2DE image registration is of great importance to align the geometric distortion and to find the correct matched spots to identify the differential protein. The image registration methods for 2DE images can be classified into three main categories: (1) landmark-based method, (2) intensity-based method, and (3) combination of landmark-based method and the intensity-based method.

The landmark-based method is used to detect the landmarks (e.g. spot, scale-invariant feature transform [SIFT], and Harris) as feature points and then finds the correspondences for estimating the transformation between two images.^4–6 The landmark-based method can be utilized to make up the non-uniform distribution of the gel spots from different gels. Rogers and Graham ⁴ proposed an iterative method to match the 2D points and combine the transform estimation. However, these methods need to extract the corresponding spots that are computationally expensive and hard to guarantee accuracy. And generally, the landmark is chosen by the transformation estimation. Also, the optimal choice of the spots for guiding the matching remains an open question that is not necessarily in conformity with the physical reality. ⁷

Intensity-based method uses some criterions (e.g. sum of squared difference, mutual information, and cross-correlation) to drive the registration.^8–11 The software Z3 using intensity-based frame is the first system for 2D gel image registration. ¹⁰ These methods can make full use of pixel information. But it is also computationally expensive.

The third frame is the landmark combing with the intensity information.^8,12 The extracted landmark constraints can compensate for appearance changes between images. And the intensity information constraints can focus on the shape changes between the spots in two images. However, the inaccuracy of the landmark and the insufficient corresponding pairs lead to deteriorate the performance for 2D gel image registration.

Although lots of methods have been proposed, very few works well for gel comparisons due to complex physical and chemical processes. The locations of proteins generally vary in different images. Figure 1 shows the two pairs of 2D gel images. (1) The red contour of ventricle region in the top images denotes the different shapes in each pair of images. (2) The location of protein spots is slightly different in images (c) and (g). And the image in (c) is the directly aligned image from (a) and (b). The image in (g) is the directly aligned image from (e) and (f). (3) It shows that the intensity range is different from the histograms (d) and (h). Not only whole gel expands but also protein spot shapes as well as the image intensities slight change in the local regions. Hence, image registration algorithm is required to have the capacity to deal with the challenge of both global and local deformation of shape and intensity for 2DE image.

Figure 1.

Examples showing the differences between gel electrophoresis images. The red contour of ventricle region in the top images denotes the different shapes in each pair of images. The second row is the difference of location and intensity range. There are two image pairs (a) and (b), (e) and (f). The location of protein spots is slightly different in image (c) and (g). And the (c) is the directly aligned image from (a) and (b). The (g) is the directly aligned image from (e) and (f). The intensity range is different from the histogram (d) and (h), and (d) is for image pairs (a) and (b), (h) is for image pairs (e) and (f).

To address the above difficulties, we propose a hierarchical transformation-guided registration framework to accurately align the two images. Since the transformation not only includes the global affine transformation but also local deformation, we estimate two mappings in a hybrid manner: (1) affine transformation model that characterizes the mapping by spatial mapping consistency to iteratively recover the corresponding point and estimates the motion distortion transformation $F$ between two images. (2) No-rigid deformation model that encodes the displacement from each point in moving image to the corresponding point in the fixing image by free-form deformation.

In our experiments, we have comprehensively evaluated the registration performance for 2DE images with comparison to the landmark-based or intensity-based registration method, including the Demons^13,14 and landmark-based method with affine transformation. Through quantitative measurements and visual inspection, our proposed alignment method outperforms other methods.

Methods

The goal of our 2DE registration is to align any two images with possible significant physical and chemical gap for further analysis of the differentially regulated protein spots. We register a moving image $M$ to a fixed image $F$ in two steps that consist of affine transformation estimation and no-rigid deformation estimation stages:

Affine transformation estimation. We first extract the SIFT feature in the moving image and fixed image, which can sparsely cover the whole image. The correspondence detection scheme is established for global consistence alignment. Further refinement of the initial correspondences is performed by using iterative spatial vector consistency mapping. Then, we can remove mismatches via transformation-guided with corresponding spatial vector consistency mapping, which can alleviate the motion distortion and make global alignment much more efficient.

No-rigid deformation estimation. After we approximately transform moving image $M$ from its native space ${\mathrm{\Omega }}_M$ to the fixed image space ${\mathrm{\Omega }}_F$ by affine transformation, we further use the no-rigid deformational model to estimate the displacement for each pixel. Since the global deformation and local deformation are almost applied to the moving image for aligning the fixed image, we employ the logic to complete the detection of different spots between the moving image and fixed image.

Affine transformation estimation stage

Considering the global alignment efficiency, we select a set of landmarks automatically with any structure features that can sparsely cover the whole image. Then, an iterative correspondence detection scheme is established for motion distortion. Thanks for completing and complementary imaging information, we can deploy the state-of-the-art landmark detector method to accurately extract the stability landmarks and obtain underlying correspondence in both images. The establishment of correspondence between images is to find the point pairs that represent the same structures in two images. Correspondence matching can be performed either on point intensities or on higher level feature descriptors. The SIFT ¹⁵ algorithm is used in our schemes. Since the approximate affine transformation highly depends on the ratio of inliers and the distribution of inliers has been established, we obtain the well correspondences set by the following two steps.

The correspondences are estimated coarsely

We employ the Best-Bin-First (BBF) algorithm ¹⁵ to estimate the correspondences between two images. And the BBF algorithm cannot guarantee that all matches are correct. Fortunately, the SIFT descriptors are invariant to image scaling and rotation. The main incorrect matches can be removed by the SIFT keypoints’ orientations. We measure the orientation difference by them. If the difference of orientations is more than a threshold, then the correspondence is discarded. Experiments show that most incorrect matches are excluded by this criterion. We calculate the distance from each mapped point to the correspondent point in pre-defined reference image to discard the mismatches. The distance is used to measure whether thecorrespondence is along the estimated spatial deformation. ¹⁶

The correspondences are refined finely using L₂E estimator

Since the transformation between gel electrophoresis images is non-rigid, we can apply the optimal transformation with ${\mathrm{L}}_2$ -minimizing estimate of transformation to discard the mismatching for refining the correspondences and estimating the transformation parameters. The ${\mathrm{L}}_2$ -minimizing estimate is robust for a significant fraction of outliers that can estimate the transform $\mathrm{ }\mathrm{f}$ and fit the inliers, which is used to build the correspondences at same time. Assuming that the noise of the inliers is given of the inliers by a normal distribution to obtain the following functional criterion^17–19

$${\text{L}}_2\text{E}\left(\text{f},{\text{σ}}^2\right)=\frac{1}{2^{\text{d}}{\left(\text{πσ}\right)}^{\text{d}/2}}-\frac{2}{\text{n}}{\displaystyle \sum }_{\text{i}}^{\text{n}}\text{Φ}\left({\text{y}}_{\text{i}}-\text{f}\left({\text{x}}_{\text{i}}\right)|0,{\text{σ}}^2\text{I}\right)$$

(1)

where $\mathrm{\Phi }\mathrm{ }$ is a multi-dimensional Gaussian function, $\mathrm{f}\mathrm{ }$ is an optimal transformation model. The transformation $\mathrm{f}\mathrm{ }$ can be obtained by minimizing the L₂E functional (1) then takes the form ¹⁹

$$\text{f}\left(\text{x}\right)={\displaystyle \sum }_{\text{i}=1}^{\text{m}}\text{Γ}\left(\text{x},{\tilde{\text{x}}}_{\text{i}}\right){\text{c}}_{\text{i}}$$

(2)

where ${\mathrm{c}}_{\mathrm{i}}$ is a $1\times \mathrm{d}$ dimensional vector. In order to avoid significant complexity problems, a sparse approximation is applied by randomly picking only a subset of size input points to have nonzero coefficients in the expansion of the solution. Adding a regularization term for imposing smooth constraint on the transformation, the L₂E functional has the following form

$${\text{L}}_2\text{E}\left(\text{C},{\text{σ}}^2\right)=\frac{1}{2^{\text{d}}{\left(\text{πσ}\right)}^{\text{d}/2}}-\frac{2}{\text{n}}{\displaystyle \sum _i^n\frac{1}{2^{\text{d}}{\left(\text{πσ}\right)}^{\text{d}/2}}{\text{e}}^{-\frac{\Vert {\text{y}}_{\text{i}}^{\text{T}}-{\text{U}}_{\text{i}},{\text{C}}^2\Vert }{2{\text{σ}}^2}}+\text{λ}{\Vert \text{f}\Vert }_{\text{Γ}}^2}$$

(3)

where $\text{λ}> 0$ controls the strength of regularization, the stabilizer ${\Vert \text{f}\Vert }_{\text{Γ}}^2$ is defined by an inner product. By choosing a diagonal decomposable kernel, ¹⁹ with determining the width of the range of interaction between samples, the L₂E functional can be expressed in the following matrix form

$$\begin{array}{l}{\text{L}}_2\text{E}\left(\text{C},{\text{σ}}^2\right)=\frac{1}{2^{\text{d}}{\left(\text{πσ}\right)}^{\text{d}/2}}-\frac{2}{\text{n}}{\displaystyle \sum _i^n\frac{1}{2^{\text{d}}{\left(\text{πσ}\right)}^{\text{d}/2}}{\text{e}}^{-\frac{\Vert {\text{y}}_{\text{i}}^{\text{T}}-{\text{U}}_{\text{i}},{\text{C}}^2\Vert }{2{\text{σ}}^2}}}\\ +\text{λtr}\left({\text{C}}^{\text{T}}\text{ΓC}\right)\end{array}$$

(4)

The derivative of the L₂E cost function is used to calculate the minimized L₂E functional

$$\frac{\partial {\text{L}}_2\text{E}}{\partial \text{C}}=\frac{2{\text{U}}^{\text{T}}\left[\text{P}\circ \left(\text{Q}\otimes 1_{1\times \text{d}}\right)\right]}{{\text{nσ}}^2{\left(2{\text{πσ}}^2\right)}^{\text{d}/2}}+2\text{λΓC}$$

(5)

where $\mathrm{P}=\mathrm{UC}-\mathrm{Y}$ and $\mathrm{Y}={\left({\mathrm{y}}_1,\dots ,{\mathrm{y}}_{\mathrm{n}}\right)}^{\mathrm{T}}$ are matrices of size $\mathrm{n}\times \mathrm{d}$ , diag $\left(·\right)$ is the diagonal of matrix, $\mathrm{Q}=\exp \left\{\mathrm{diag}\left({\mathrm{PP}}^{\mathrm{T}}\right)∕2{\mathrm{\sigma }}^2\right\}$ is an $1\times \mathrm{d}$ dimensional vector, $1_{1\times \mathrm{d}}$ is an $1\times \mathrm{d}$ dimensional row vector of all ones, “ $\circ \mathrm{”}$ denotes the Hadamard product, and $\mathrm{“}⨂\mathrm{”}$ denotes the Kronecker product.

Although the robust estimation algorithm can calculate the global transformation matrix correctly with less mismatches, the alignment accuracy cannot be guaranteed, especially, when shape changes in the local region. We take the further procedure for improving alignment accuracy after L₂E estimator removing the mismatches.

No-rigid deformation estimation stage

In order to relieve the intensity difference in two images, we utilize the histogram matching to normalize two images with the same local illumination in whole image before no-rigid deformation. Given two images, moving image $M$ and fixed image $F$ , first their histograms are calculated. Then the cumulative distribution functions of the two image histograms are calculated: $M_h$ for the moving image and $F_h$ for the fixed image. And the gray level $G_2$ is needed with $G_1\in [0,\mathrm{ }255]$ , and $G_2$ can be obtained from the following format

$$M_h\left(G_1\right)=F_h\left(G_2\right)$$

(6)

The goal of this is to obtain the histogram matching function

$$H\left(G_1\right)=G_2$$

(7)

Finally, the objective function $H$ is used on each pixel of the fixed image.

Next step is to estimate dense deformation for each pixel. Since the affine transformation is done by above stage, the next stage is to finish an optimization that aims at calculating the displacement of each pixel to obtain an intelligent alignment of the images.

Let $\mathrm{f}$ be the intensity of moving image $M$ and $\mathrm{m}$ be the intensity of fixed image $F$ . Using the edges in the static image as passive internal force, Wang et al. ²⁰ add an equation with the image edge forces of the moving image that improves the registration convergence speed and stability. The estimated displacement (velocity) of $\mathrm{u}$ required for point to match the corresponding point in $M$ is given as^20,21

$$\text{u}=\frac{\left(\text{m}-\text{f}\right)\nabla \text{f}}{{\left|\nabla \text{f}\right|}^2+{\text{α}}^2{\left(\text{m}-\text{f}\right)}^2}+\frac{\left(\text{m}-\text{f}\right)\nabla \text{m}}{{\left|\nabla \text{m}\right|}^2+{\text{α}}^2{\left(\text{m}-\text{f}\right)}^2}$$

(8)

where $\mathrm{u}=\left({\mathrm{u}}_{\mathrm{x}},{\mathrm{u}}_{\mathrm{y}}\right)$ is the 2D displacement in point $\left(\mathrm{x},\mathrm{y}\right)$ , $\nabla \mathrm{f}$ and $\nabla \mathrm{m}$ are the gradient of the fixed and moving images. The forces include internal edge-based force $\nabla \mathrm{f}$ , the external force $\left(\mathrm{m}-\mathrm{f}\right)$ , and ${\left(\mathrm{m}-\mathrm{f}\right)}^2$ . However, the forces $\left(\mathrm{m}-\mathrm{f}\right)$ and ${\left(\mathrm{m}-\mathrm{f}\right)}^2$ just indicate the local intensity differences (each pixel) between the two images.

Results

Dataset

In total, there are three group protein images dataset: KO-WHY1, KO1-3, WT, which are obtained from Fujian Provincial Key Laboratory of Plant Functional Biology. Total proteins were separated by sodium dodecyl sulfate-polyacrylamide gel electrophoresis (SDS-PAGE) after isoelectric focus (IEF) and then the gels were stained by Comassie Blue R250. And 56 image pairs from three group protein images dataset are chosen in the following experiments, where each pair has 2D protein separation patterns visualized.

Setting of parameters

Demons and our improved Demons

The iteration number is 20. For L₂E, since it is non-convexity that the L₂-minimizing estimate with the multi-dimensional Gaussian function in our method, the deterministic annealing is employed by controlling the annealing rate γ. The initialization of γ, λ, and σ² are 0.5, 0.1, and 0.05, respectively. ¹⁸

Registering 2D gel image and results

To evaluate the proposed framework, first the number of feature matching with different frames is compared. The initial matches are obtained by SIFT, the refined matches by the SIFT keypoints’ orientations, and the final matches are shown in Figure 2. In order to visualization of results, the matches are drawn by one-eighth sparse. Actually, the number of initial matches is 987 in (a) and it contains many mismatches, the number of refined correspondences is 936 in (b), the number of final correspondences is 823 in (c). The mismatches have been removed to a large extent. It is apparent in Figure 2 that the SIFT keypoints’ orientations can be used to delete the apparent mismatches. And the L₂E frame does well to obtain the correspondences without mismatches in this example.

Figure 2.

The correspondences with different steps for our frame: (a) the correspondences (987 matches) are obtained by BBF frame, (b) the correspondences (936 matches) are refined by the SIFT keypoints’ orientations, and (c) the final correspondences (823 matches) are refined by the L₂E frame.

We compare the correspondences obtained by different frames for eight pairs of image data. For example, RANSAC, L₂E, and the proposed framework are used for correspondence estimation between the images. And all original matches are from the SIFT. Table 1 shows (1) the obtained matched correspondences (O) by calculating Euclidean distance (ED) between the matched point and registration point and (2) the final correspondence fraction of correct matches (F (%)) with different estimators, which are measured by the ED between the correspondences. If the distance D is less than a threshold T (the value of T is 0.5 in our experiment), the correspondence is a correct match. The best result for each dataset is shown in bold.

Table 1.

Comparison of the number of obtained (O) matched correspondences as well as fraction of correct matches (F (%)) with different estimators for intra subjects.

Data	Original matches	RANSAC		L2E		Our method
		O	F	O	F	O	F
No. 1	874	644	81.1	667	95.8	663	100
No. 2	670	642	89.0	452	95.3	448	100
No. 3	871	651	67.5	708	95.2	700	100
No. 4	750	593	84.1	555	98.3	552	100
No. 5	1378	1201	5.6	1223	100	1205	100
No. 6	674	468	17.9	472	100	466	100
No. 7	470	249	6.7	278	100	273	100
No. 8	590	434	14.3	403	97.3	398	97.1

Note: The best result for each row is shown in bold.

And to validate the proposed registration framework, we compare the accuracy of registering with Demons, landmark-based method with affine transformation, and the proposed hierarchical transformation-guided registration for 2DE image. The preprocessing is same that the affine transformation is done by our proposed framework. The result of landmark-based method is obtained using the refined correspondences as control points. The registration accuracy is measured through the following four qualities of alignment image: mutual information (MI)

$$\text{MI}=\text{h}\left(\text{f}\right)+\text{h}\left(\text{m}\right)-\text{h}\left(\text{f},\text{m}\right)$$

(9)

where h(f) and h(m) are the entropies for images f and m, respectively. h(f, m) is the joint entropy for images f and m; peak signal-to-noise ratio (PSNR)

$$\text{PSNR}=10{\text{log}}_{10}\left(\frac{{\text{Q}}^2}{\frac{1}{\text{C}}{\Vert \text{f}-\widehat{\text{m}}\Vert }_2^2}\right)$$

(10)

where $\widehat{\mathrm{m}}$ is the deformable image from $\mathrm{m}$ , $\mathrm{Q}$ is the maximal intensity value of $\mathrm{f}$ and $\widehat{\mathrm{m}}$ , $\mathrm{C}$ is the number of pixels in the image; and root mean square error (RMSE)

$$\text{RMSE}=\sqrt{\frac{{\left(\text{f}-\widehat{\text{m}}\right)}^2}{\text{C}}}$$

(11)

In Table 2, the registration accuracy is listed with different registration methods, which is measured by MI, PSNR, and RMSE for all image pairs. The values for those criterions are obtained from calculating the measure between the registration point and the matching point (the best result is that the two locations overlap). The values of MI and PSNR are as higher as better, and the higher value indicates that the results have more relevant items. The value of RMSE is as lower as better, and the lower value indicates that the result is more accurate. The best result for each row is shown in bold. It is apparent that the proposed registration framework achieves more accurate registration results than other methods.

Table 2.

The mean and standard deviation of MI, PSNR, and RMSE with different registration methods.

Method	AFF	Demons	Our method
MI	1.6598 ± 0.0377	1.7247 ± 0.0444	1.7303 ± 0.0480
PSNR	18.6209 ± 3.1783	20.0601 ± 3.2563	20.2006 ± 3.3514
RMSE	0.1582 ± 0.0122	0.1344 ± 0.0142	0.1324 ± 0.0155

Note: The best result for each row is shown in bold. MI: mutual information; PSNR: peak signal-to-noise ratio; RMSE: root mean square error.

As shown in Figure 3, the results are shown for four pairs of images, and each pair of images is put in a blue rectangle. And in each blue rectangle, the first row is the two original images and second row is the directly aligned images (the left of the second row) and registered images (the right of the second row). In the second row, the images are aligned by single channel image. In the aligned images, the image colored in red is from the fixed image, the image colored in green is from the moving image, and the image colored in yellow is the overlapped spots. If the aligned images show many yellow spots, the spot is aligned well and the image pairs have greater similarity. The spots can be aligned well in the registered image pairs. The different spots and misaligned spots can be identified easily, which is especially in the regions indicated by colored circles, where No. 4 image pair has an apparent misaligned point pair as shown by a green circle. It is clear that the proposed registration method achieves more accurate registration results for the further analysis of the differentially regulated protein spots.

Figure 3.

Registration results of datasets No. 1 (a), No. 2 (b), No. 3 (c), and No. 4 (d), respectively. Top row shows the original images, bottom row shows the unregistered overlap image and registered overlap image. In the second row, the images are aligned by single channel image. In the aligned images, the image colored in red is from the fixed image, the image colored in green is from the moving image, and the image colored in yellow is the overlapped spots.

Conclusion

In this paper, we have presented a hierarchical transformation-guided registration method for 2DE image. In order to address the changes of protein locations generally and shape development, we deploy affine transformation model to obtain the global alignment and use the deformable model to alleviate the local shape differences. Since there is no rigid deformation between the gel electrophoresis images, we combine the conventional non-parametric registration method – demons and the entropy of image to complete the deformable registration to alleviate the appearance of differences. We have intensively evaluated the registration accuracy on 2DE image and compared with landmark-based method and intensity-based method, where our proposed method achieves highest registration accuracy compared with other counterpart registration methods.