PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of bioinfoLink to Publisher's site
 
Bioinformatics. 2010 June 15; 26(12): i29–i37.
Published online 2010 June 1. doi:  10.1093/bioinformatics/btq194
PMCID: PMC2881379

A spectral graph theoretic approach to quantification and calibration of collective morphological differences in cell images

Abstract

Motivation: High-throughput image-based assay technologies can rapidly produce a large number of cell images for drug screening, but data analysis is still a major bottleneck that limits their utility. Quantifying a wide variety of morphological differences observed in cell images under different drug influences is still a challenging task because the result can be highly sensitive to sampling and noise.

Results: We propose a graph-based approach to cell image analysis. We define graph transition energy to quantify morphological differences between image sets. A spectral graph theoretic regularization is applied to transform the feature space based on training examples of extremely different images to calibrate the quantification. Calibration is essential for a practical quantification method because we need to measure the confidence of the quantification. We applied our method to quantify the degree of partial fragmentation of mitochondria in collections of fluorescent cell images. We show that with transformation, the quantification can be more accurate and sensitive than that without transformation. We also show that our method outperforms competing methods, including neighbourhood component analysis and the multi-variate drug profiling method by Loo et al. We illustrate its utility with a study of Annonaceous acetogenins, a family of compounds with drug potential. Our result reveals that squamocin induces more fragmented mitochondria than muricin A.

Availability: Mitochondrial cell images, their corresponding feature sets (SSLF and WSLF) and the source code of our proposed method are available at http://aiia.iis.sinica.edu.tw/.

Contact: wt.ude.acinis.sii@nannuhc

Supplementary information: Supplementary data are available at Bioinformatics online.

1 INTRODUCTION

Recently, high-throughput image-based assay technologies, or high-content analysis, have become a useful tool for drug discovery (Jones et al., 2009; Lang et al., 2006), small molecule screen (Tanaka et al., 2005), subcellular localization (Huang and Murphy, 2004; Lin et al., 2007), etc. These technologies make it possible to visualize, trace and quantify cellular morphological changes in high resolution and play an increasingly crucial role to the understanding of biological processes.

Since the cost of acquiring a large number of microscopic cell images is decreasing, the analysis tools must be powerful. Consider a study aiming at correlating drug influences to morphological changes of mitochondria. This topic has become increasingly important mainly because of its relation to apoptosis and aging process. Figure 1 shows some representative cell images with different levels of fragmentation of mitochondria. As we can see, intact or completely fragmented mitochondria are relatively easy to recognize, but it is difficult to quantify partial fragmentation. In fact, in addition to fragmentation, mitochondria may become tubular, shortened, elongated or swollen. All these morphological changes are important biological indicators (Brooks et al., 2007; Lee et al., 2004; Tamai et al., 2008). However, off-the-shelf analysis software only distinguishes morphological patterns at cellular level. No partial change or mixture of patterns at subcellular level can be readily detected. Quantification requires human inspection, which is infeasible for high-throughput screening (Taguchi et al., 2007; Zhou and Wong, 2006).

Fig. 1.
Representative microscopic fluorescent images of single cells with different levels of fragmentation of mitochondria. Intact mitochondria forms filamentous networks (top row), completely fragmented mitochondria have round shape (bottom row) and partially ...

Our goal is to derive an objective and reliable quantification tool to correlate drug influences and morphological differences at cellular or subcellular levels. We assume that cell images are characterized by a feature set extracted from the images and can be considered as a data point in the feature space. A straightforward quantification is to measure the Euclidean distance between data points or apply other sophisticated similarity metrics. However, there are key challenges that must be addressed.

  • A wide variety of morphological categories must be considered.
  • Influence of a drug is a random variable, the difference must be estimated between pairs of cell-image collections. Averaging pairwise differences of single images is not appropriate here because of the presence of noise and outliers.
  • Whether a selected similarity metric can actually quantify the difference that matches the purposes of the study depends heavily on the selected features.
  • Calibration is essential for a practical quantification method because we need to measure the confidence of the quantification.

Our solution is a graph-based approach. Previously, graph-based approaches have been shown to be effective for clustering, semi-supervised learning and image segmentation (Belkin and Niyogi, 2003b; Zhu and Ghahramani, 2002; Zhu et al., 2003). In our approach, graph transition energy is defined to quantify the similarity between collections of images. To accommodate a wide variety and combinations of morphological differences to be quantified, we adopt a supervised paradigm where two sets of extremely different cell images are assumed to be given as training examples. For example, to quantify mitochondrial fragmentation, sets of cell images with intact and completely fragmented mitochondria are assumed to be given. Figure 2 illustrates this graph-based approach. By applying a spectral graph theoretic regularization (Chung, 1997), we developed a method to transform the feature space based on the training examples so that regularized graph energy between data points of extremely different morphology is minimized. In this way, calibration of the quantification can be achieved. Then we can quantify a new set of cell images by computing the graph transition energy between the set and training examples in the transformed feature space. Experimental results show that our method quantifies the morphological difference more accurately and sensitively than that without transformation. Results also show that our method outperforms competing methods, including neighbourhood component analysis (NCA; Goldberger et al., 2005) and the multi-variate drug profiling method by Loo et al. (2007). Finally, we illustrate the utility of our method with a study of Annonaceous acetogenins and their impact to mitochondrial fragmentation. Our result reveals that squamocin induces more fragmented mitochondria than muricin A.

Fig. 2.
Given two sets of cell images with extremely different morphology patterns as training examples; our method transforms the feature space so that regularized graph transition energy between the two sets is minimized. Then we can quantify the morphological ...

The remainder of this article is organized as follows. Section 2 reviews related work. Section 3 presents our method. Section 4 reports experimental evaluation of our method and an application to the study of Annonaceous acetogenins. The last section concludes and discusses our future work.

2 REVIEW OF PREVIOUS WORK

We briefly review previous work on image-based approaches for drug screening, mitochondrial fragmentation, graph-based approaches and feature space transformation.

2.1 Image-based approach for drug screening

Recently, high-throughput image-based approaches have received great attention for drug screening (Carpenter, 2007; Lang et al., 2006; Loo et al., 2007). Among them, Loo et al. (2007) proposed an image-based multivariate profiling method for drug screening. In their method, support vector machines (SVMs) are applied to establish a hyperplane in the feature space between cell images representing control and images of cells under different perturbation. Then the unit normal vector of the hyperplane is used as a multivariate profile to indicate the phenotypic direction of the greatest separation between the two cellular populations. By clustering these profiles, a new compound with unknown effect can be associated with compounds with known effects to infer the effects of the new compound. Clustering these profiles may also reveal new usage of an old drug if its profile is similar to another drug with target effects. Their method is interesting in that it can consider a large number of features and is not study-dependent. However, their approach fails when the images in the feature space are not linearly separable. Also, the images may distribute in the feature space so widely that the hyperplane may have a very high sampling variance. In this article, we report an experimental comparison of our graph-based method with their approach in Section 4.2.

2.2 Mitochondrial fragmentation

Mitochondria are essential organelles in eukaryotic cells, because they play a central role in energy metabolism and the integration of diverse apoptotic stimuli. Recent investigations reveal that dynamics of mitochondria morphology is important during apoptosis as the organelle undergoes massive fragmentation. Yaffe (1999) suggested that over fragmentation of mitochondria may result in breakage of mtDNA that causes excess free radicals due to impaired functions in respiration. Other reports show that mitochondrial fragmentation is an early step of apoptosis (Desagher and Martinou, 2000; Frank et al., 2001; Karbowski and Youle, 2003). However, mechanisms behind the relationship between mitochondrial fragmentation and apoptosis is still unclear (Jeong and Seol, 2008). High-content screening can play a key role in revealing their relationship and shed new light on the discovery of potential drugs for the treatment of related diseases, such as neurodegenerative diseases and diabetes.

2.3 Graph-based approach

Graph-based approaches have been successfully applied in various machine-learning problems including classification (Belkin and Niyogi, 2003b; Goldberger et al., 2005; Zhan et al., 2009; Zhu et al., 2003), spectral clustering (Azran and Ghahramani, 2006) and dimension reduction (Belkin and Niyogi, 2003a; Roweis and Saul, 2000; Tenenbaum et al., 2000). These methods create a graph whose nodes correspond to data points, while the edge weights encode the similarity between each pair of data points. It is crucial to choose a proper distance metric to estimate the similarity as the performance of graph-based models depends considerably on the similarity measure of the graph.

Zhu et al. (2003) uses a Gaussian random field model to construct a weighted graph representing labeled and unlabeled data for semi-supervised learning. For dimension reduction, Isomap (Tenenbaum et al., 2000) applies geodesic distance to create a graph that captures the manifold structure to recover the intrinsic dimension of the data. Though geodesic distance can serve as a useful quantification of a given type of morphological difference, the isomap method requires input data uniformly distributed in the feature space with a high density to accurately recover the manifold, a requirement that is infeasible in most of biological study because the distribution of cell images is usually very skew.

Other promising applications of the graph-based approaches in bioinformatics include immune cells detection (Chang et al., 2008), cell images segmentation (Isfahani et al., 2008) and protein function prediction (Borgwardt et al., 2005).

2.4 Feature space transformation

Feature space transformation is widely applied to data visualization and learning algorithm enhancement. Here, we reviewed two methods, linear discriminant analysis LDA and NCA, that are closely related to our method. Fisher LDA uses the label information to find informative projection such that the separation of data of different classes can be maximized. To that end, LDA tries to maximize the inter-class scatter matrix and minimize the intra-class scatter matrix simultaneously. However, LDA suffers from a small sample size problem when dealing with high-dimensional data, the intra-class scatter matrix can be nearly singular (Chen et al., 2000). In our approach, we re-weight the feature space and make no attempt to compute projections as in LDA. Therefore, the problem is not an issue here because no inverse matrix is necessary. Also, LDA is designed for classification. Quantifying difference between cell images is more of a regression problem than classification. NCA (Goldberger et al., 2005) learns a linear transformation of the feature space to optimize leave-one-out (LOO) performance on training data for k-nearest neighbor (k-NN) classifiers. NCA is a non-parametric learning method and makes no assumptions about data distributions. Promising performance of NCA have been shown in many applications including face recognition (Butman and Goldberger, 2008) and hyperspectral classification (Weizman and Goldberger, 2007). Since our quantification method is based on a distance-weighted k-NN, we empirically compared our method with NCA and report the results in Section 4.1.2.naveenh

3 METHODS

Our task is to re-weight the feature space to separate the data points with different labels as far apart as possible. We consider a graph-based model to solve our problem.

3.1 Adjacency graph

Let {(xi, yi)}[set membership]X × Y, i=1,…, n be a set of data points where X[set membership]RD is an arbitrary feature space and Y a finite set of labels. In a graph-based model, a graph G := {V, E} is constructed for the data points, where V corresponds to data points and eij[set membership]E connects xi and xj. A non-negative weight is assigned to eij to build a weight matrix W by a radial basis function such that:

equation image
(1)

where σd[set membership]σ is the length scale for the d-th coordinate in the feature space. A common practice in graph-based methods is to consider weights of nearby data points only. Other links will be assigned zero weights. Given the weight matrix W, the graph energy (Zhu and Ghahramani, 2002) can be defined as a function to measure the stability of a set of data points:

equation image
(2)

where I(ϕ) denotes the indicator function whose value is 1 if statement ϕ is true and 0 otherwise.

3.2 Cheeger's constant

Our goal is to transform the feature space by re-weighting σ so that the Euclidean distance in the transformed space correlates with the difference of the labels of data points. The transformation requires a set of training data points with different labels. The idea is to re-scale σd such that the energy is minimized. To avoid overfitting, regularization is required. Well-known regularization includes graph Laplacian (Belkin and Niyogi, 2003b) and Cheeger's constant (Chang et al., 2008). It has been shown that Cheeger's constant effectively removes noise, like Laplacian, but preserves better sharp boundaries needed for classification (Chang and Moura, 2007; Chung, 1997). Therefore, we chose Cheeger regularization for our problem.

Let f : XY be a function that labels data points in G. Cheeger regularization is defined as

equation image
(3)

where D=diag(di) is a diagonal matrix with di=∑jwij, L=DW is the Laplacian of G, 1=(1,…, 1)T and β is the non-negative weight. We note that fTLf is the graph Laplacian regularization given in (Belkin and Niyogi, 2003b). In our problem, f is simply the known label assignment y of training data. Mapping two labels in the training data to {0, 1}, we have

equation image

To extend its applicability to multiple label problems, we modified the second term fTD1 in Equation (3) so that we only count if the labels are the same or not. That is, An external file that holds a picture, illustration, etc.
Object name is btq194i1.jpg Therefore, the Cheeger regularization is modified as

equation image
(4)

3.3 Graph transition energy

We embed the constraint that scaling factors sum to one for each i into Equation (4) by replacing the weight matrix W with a transition matrix T where An external file that holds a picture, illustration, etc.
Object name is btq194i2.jpg. With the normalization term, T is usually asymmetric. By doing so, we also need to modify the definition of graph energy function given in Equation (2) by replacing W with T, which leads to the Graph Transition Energy, defined as a function of σ:

equation image
(5)

A low E indicates that G is stable because data points with different labels are far apart, and vice versa. Therefore, given two sets of cell images as feature vectors, if their graph transition energy is low than they are quite different, and vice versa.

3.4 Feature space transformation

We continue our derivation of an objective function from Cheeger's constant. In addition to embedding the constraint to the objective function, we also re-parameterize An external file that holds a picture, illustration, etc.
Object name is btq194i3.jpg in Equation (1) to enforce another constraint so that σd≥0, d=1,…, D. Then

equation image

Let κ be a vector with κd's as its components. We complete the derivation of our objective function as a function of κ:

equation image
(6)

Minimizing J creates a transformed feature space with wide label separation and therefore calibrates the quantification with the graph transition energy. Any minimization method can be applied to minimize J. We simply used a steepest gradient descent method in our implementation.

The gradient for the objective function is

equation image
(7)

In Equation (7), the computation of An external file that holds a picture, illustration, etc.
Object name is btq194i4.jpg by using chain rule is

equation image
(8)

where An external file that holds a picture, illustration, etc.
Object name is btq194i5.jpg is

equation image
(9)

and An external file that holds a picture, illustration, etc.
Object name is btq194i6.jpg is

equation image
(9)

Substituting Equation (9) and (10) into (8), we obtain

equation image
(10)

Recall that each data point in the graph only connects to its nearest neighbors. This approximation is sufficient because in a Gaussian energy model, links of points too far apart are nearly zero. But we cannot rule out the possibility that a pair of points becomes closer with a new κ after optimization. Therefore, we design an iteration method where each iteration applies a gradient descent update to obtain a new κ. Then the graph is rewired by connecting each point with its nearest neighbors after each iteration. The iteration repeats until convergence.

We can derive a lower bound of the objective function as follows.

Proposition 3.1. —

If there exists a κ such that J(κ)=−λ · n, where n is the data size, then the data in the transformed feature space can be separated completely. That is, the transition probability Pr(xixj)≥0 if yi=yj, otherwise Pr(xixj)=0

Proof. —

Because under perfect label separation, we have Tij≥0 if yi=yj and 0 otherwise. In that case J(κ)≥−λ · n. [filled square]

According to Proposition 3.1, the value of J(κ) reveals useful information to determine if the data with different labels or conditions can be well-separated or not given the set of features. If J(κ) is still far away from its lower bound after the re-weighting process, the feature set may not be sufficient or appropriate for representing the data points. The re-weighted κ also provides useful information to rank important of the features.

3.5 Quantifying differences

We summarize our method as follows:

  1. Define the label space Y. For example, we may define Y as different levels of mitochondrial fragmentation, days after drug applications or drug dosages.
  2. Acquire data points (i.e. cell images) with labels of the two extremely different cases (e.g. mitochondrial intact and complete fragmentation).
  3. Extract features for X and minimize Equation (6) to obtain a transformed feature space.
  4. Then we can acquire more data points and apply Equation (5) to quantify the difference.

3.6 Feature extraction

We extracted two types of features from each single cell image to serve as the input in our experiments. Strong detectors are knowledge-driven features that are supposed to provide strong hints, while weak detectors are randomly extracted patterns to allow the learning algorithm to consider subtle characteristics of a class. This set of features has been used in (Lin et al., 2007) on recognizing fluorescent protein-tagged subcellular organelles in cell images, including mitochondria.

3.6.1 Strong detectors

Our strong detectors consist of both geometric and texture features. Table S1 in the supplementary data gives the list of geometric features. The texture features are extracted based on the gray-level co-occurrence matrix (GLCM) method proposed by Haralick (1979). Twelve GLCM with distances 1, 2 and 10 and angles 0, 45, 90 and 135 are applied to the bi-leveled images. Then, various co-occurrence quantities including entropy, energy, contrast, homogeneity and correlation can be evaluated from the co-occurrence matrix to produce our texture features. The definition of these quantities are given in Supplementary Table S2.

The above method extracts a total of 155 geometric features and 500 texture features from a cell image. We conducted a backward stepwise discriminant analysis to select 132 most discriminant features as our final set of strong detectors.

3.6.2 Weak detectors

The weak detectors for each image are extracted in four steps as follows.

  1. Randomly pick five images of each class as templates.
  2. Randomly extract a set of eight fragments from each template. The fragments vary in size from 9 × 9 to 25 × 25.
  3. Convolve each fragment i with a set of four filters. This set includes the original image, x derivative, y derivative and a Laplacian filter.
  4. Then for a given image either for training or testing, apply normalized cross-correlation between the given image and the fragment i to find where the fragment i occurs and then record the maximum correlation as the i-th component in the vector of weak detectors for the given image.

In this study, we have two extreme classes of mitochondrial images: intact and fragmented completely. We can generate 5(template) × 8(fragment) × 4(filter) × 2(class)=320 weak detectors for each image. If we want to generate more features, we just increase the number of templates, fragments or filters.

3.7 Mito-Q

Because we applied our method to detect the mitochondrial fragmentation in microscopy cell images, we have developed an auxiliary tool called Mito-Q to assist human inspectors to visually quantify the percentage of fragmentation of mitochondria for each cell image as our golden standard. Mito-Q helps human inspectors to filter out noise and poorly focused light spots in an input image and identify which patterns correspond to mitochondria. A morphological filter (Serra, 1998), with disk structuring element (of size 2), and power-law transformation with the basic form s=crγ (Gonzalez and Woods, 2002) are used to enhance the mitochondria objects. A bi-level operator with an optimal threshold TH, selected by the human inspector, is applied to generate the template of mitochondrial, as shown in Figure 3. At last, Mito-Q will generate mitochondrial features for evaluating the percentage of fragmentation. Though Mito-Q is helpful, it takes a skillful human inspector to generate desirable output scores and is not designed to replace automated quantification methods.

Fig. 3.
The results of Mito-Q image processing. (a) The original image; (b) the image processed by the morphological filter; (c) the image after power-law transformation; (d) the locations of mitochondrial identified by Mito-Q.

4 RESULTS AND DISCUSSION

This section reports the empirical evaluation of our feature space transformation method.

4.1 Effectiveness of feature space transformation

We first evaluate whether our graph-based feature space transformation method actually results in a new space that characterizes the difference between data points better. To this end, we evaluate whether our method improves a semi-supervised learning task and a k-NN regression task, both for cell-image analysis.

4.1.1 Improving semi-supervised learning

We applied our feature space transformation method to improve the classification performance of Zhu et al.'s (2003) semi-supervised learning method for the task of recognizing subcellular organelle patterns in the HeLa cell image data set. We used the 2D HeLa image data set available from the Murphy lab (See http://murphylab.web.cmu.edu/data/2Dhela_images.html). This data set contains 862 images. Each image contains a single cell with exactly one of the 10 distinct subcellular organelles tagged by a fluorescent protein. The feature set that we used was the superset of SLF16, which consists of 47 features selected by stepwise discriminant analysis (SDA) from a set of 180 features (See http://murphylab.web.cmu.edu/data/2Dhela_images_download.html). Among them, six are DNA-channel features. We used this 180-feature set of the HeLa images in our first experiment.

We split the test images into two subsets. One of the subsets was used as the labeled data and the other as unlabeled data. When our method was applied, the labeled data were also used to transform the feature space. We must first determine k, the number of the nearest neighbors to construct a graph. We also wanted to see if our method is sensitive to k. We selected the k-value from {6, 8, 10, 12, 14} according to their classification performance. We repeated the experiments on each data set for 10 runs with random partitions of 60% training/40% test instances. It can be observed from Figure 4a that there is no significant difference for different k. Therefore, we arbitrarily assigned k=8, the best performing setting, to construct the graph in our experiments. We also found that our method is insensitive to choices of λ in Equation (6) over a wide range. In all of our experiments, we simply set λ=0.1.

Fig. 4.
(a) Influence of the k nearest neighbors. (b) The performance of semi-supervised learning before and after learning κ in the HeLa cell image data set.

The result of semi-supervised learning is given in Figure 4b, which shows that though both semi-supervised methods improve their performance as the number of labeled training examples increases, applying our transformation clearly improves the performance of the plain semi-supervised learning method. The best accuracy is about 87%, about as good as a supervised neural net reported in (Huang and Murphy, 2004) for the same set of images. We expect that with more unlabeled images, semi-supervised learning with transformed feature space has the potential to achieve an accuracy as good as the best results by supervised learning.

Moreover, we can compute graph transition energy between images of different subcellular organelles. The resulting energy table can serve as a measure of pairwise similarity. Then we can applied multi-dimensional scaling (MDS) to help us visualize the morphological similarity of the shape of 10 distinct subcellular organelles as shown in Figure 5, where we can see that two Golgi types (Gpp130 and Giantin) are particularly challenging to distinguish, so are mitochondria and tubulin.

Fig. 5.
Visualizing shape similarity of 10 distinct subcellular organelles with a MDS plot derived from pair-wise graph transition energy.

4.1.2 Estimating partial fragmentation of mitochondria

Now we show that our approach can also be applied to improve the estimation in the fragmentation level of mitochondria in a single cell without actually identifying each mitochondrion object in a cell image. Since estimating fragmentation of a mitochondrion is a regression problem, we applied distance weighted k-NN for this task. We compared our feature space transformation method with NCA (Goldberger et al., 2005), which is designed to transform the feature space to minimize the error of k-NN.

To apply our method, in addition to training data for obtaining an optimal κ, we need other labeled data in the feature space to approximate the manifold of the data. With the labeled data, we can quantify each test image by the weighted average of its nearest neighbors:

equation image
(12)

where fk is the label of the k-th labeled data and K is a pre-defined constant for the number of the nearest neighbors. In our test case, fk is simply the score by human inspectors with aid of Mito-Q and K=6. We used the same κ obtained from mostly impact (MI) and mostly fragmented completely (MC) and all 392 images as the labeled data. We then acquired an additional set of 43 new images to evaluate the above method.

Table 1 shows the performance comparison in terms of correlation coefficient with the human inspectors' score and mean-square-error (MSE). The result shows that our method outperforms NCA by achieving the highest correlation coefficient and the lowest MSE. The result also shows that the transformation by learning κ effectively improves the quality of the estimation.

Table 1.
Comparing the performance of estimating partial fragmentation of mitochondria in individual cell images

4.2 Quantifying collective morphological difference

We now consider the main task that our approach is designed for. That is to quantify morphological difference between collections of cell images. This task is particularly useful for drug screening, as different image sets represent the results of the treatment of different drugs. The same approach can be applied to study effects of the same drug under different dosages and time courses.

To test our method for this task, we acquired a collection of 392 single cell images with different percentage of mitochondria fragmentation. We sorted the cell images according to the output scores of Mito-Q and labeled the top 36 images as MI and bottom 36 as MC to construct two extreme cases of mitochondria. The rest of the images (320) are labeled as mitochondrial fragmented partially (MP) in three equal sized subsets (MP1, MP2 and MP3), from the least fragmented to the most. For each cell image, we extracted 452 features that include 132 strong (e.g. morphology, texture, moment and intensity features) and 320 weak detectors, as described in Section 3.6.

We sampled with replacement from MI and MC to create 11 treatment sets with different mixture proportion. The k-th set was a mixture of (11−k) × 10% MI images with (k−1) × 10% MC. Then we regarded raw MI as the control set to compute the difference between each mixture set and control set. Note that the percentage here refers to the proportion of images sampled from MI or MC, not the degree of mitochondrial fragmentation. We compared our approach with the approach proposed by Loo et al. (2007). Recall that in their approach, SVM is applied to establish a hyperplane in the feature space to distinguish images of cells under different perturbation. Then the unit normal vector of the hyperplane is used as a profile and cosine similarity between profiles is used to quantify the difference. To duplicate their method, we obtain a profile between control and each treatment set. In our graph-based method, we use graph transition energy defined in Equation (5) to quantify the difference in the transformed feature space by minimizing Equation (6) with MI and MC as the training examples. Because of a small set of training data, we determine the number of the nearest neighbors from {2, 4, 6,…, 20,∞} via leave-one-out validation based on the classification performance of Zhu et al.'s (2003) semi-supervised learning method. It is observed from Figure 6 that the accuracy rate is robust to the number of the nearest neighbors. Therefore, we just use six nearest neighbors for this task.

Fig. 6.
Influence of the k-NN.

We repeated the above trial 1000 times to approximate the distribution of the similarity of profiles and the distribution of the energy values to determine which method has higher power of discrimination. In our experiment, we also compared the energy distribution obtained in the feature space without transformation to demonstrate the effectiveness of the re-weighting process. Figure 7 plots the distribution of the resulting 1000 quantification for the 1st (pure MI, 0% MC), 6th (mixture of MI and MC with equal proportions, 50% MC) and 11th (pure MC, 100% MC) treatment sets by different approaches. The result shows that the approach proposed by Loo et al. (2007) yields a large overlapping area among these distributions, implying that it is quite likely that their approach will fail to quantify that a sub-sample of MI is more similar to raw MI (control) than a mixture set of MI and MC. In contrast, our method yields a small overlapping area and will quantify the difference in a correct order. Though without transformation, the overlapping area is still small, but the mean between pure MI and pure MC is closer than that with transformation. Next, for all treatment sets, we estimate the KL divergence to measure the difference of the distributions of the quantification. For an ideal method, the KL divergence should increase gradually and then surge immediately after the proportion of MC is increased to more than 50%. Figure 8 plots the result, which reveals that our approach has the highest discrimination power and that the transformation of the feature space actually improves the quality of quantification substantially.

Fig. 7.
Distributions of the similarity of profiles as proposed by Loo et al. (2007) (a) and distributions of energy estimated by the graph-based method with (c) and without (b) feature space transformation for the 1st (pure MI, 0% MC), 6th (mixture of MI and ...
Fig. 8.
KL divergence of the distributions of the quantification of the difference between control (pure MI) and the treatment sets with increasing proportion of MC.

Next, we evaluated whether the quantification by our approach actually reflects the morphological difference. We computed the energy between each partially fragmented subset and MI/MC in the transformed feature space. The result is shown in Table 2. As expected, the energy values are in the same order as their score of fragmentation measured by human inspectors with aid of Mito-Q. In our approach, the value of the objective function J(κ) indicates how well the data are separated in the transformed feature space. If J(κ) approaches the lower bound of the objective function given in Proposition 3.1, then the data of extremely different classes can be separated well, which implies that the feature set characterizes the difference well, and that the confidence of our quantification is high. In this case, the lower bound of J is equal to −7.2 (λ=0.1 and n=72). With the optimal κ learned from MI and MC sets, J(κ) can reach −6.040 in the transformed space, suggesting a high confidence for the quantification.

Table 2.
Graph transition energy between the treatment sets with different levels of mitochondrial fragmentation

4.3 Analysis of Annonaceous acetogenins

We applied our method to the study of Annonaceous acetogenins. More than 250 compounds in the Annonaceous acetogenin family have been isolated from Annonaceae plants. In this article, we report our results with muricin A and squamocin, two members in the family with the potential of triggering mitochondria fragmentation and apoptosis in tumor cells. We measured the effects of the compounds at different time points by quantifying mitochondrial morphology changes on Chinese hamster ovary (CHO) cells for 3 days. Besides, we consider the effect of DMSO as the control. For each group (a compound–day combination), we acquired images, segmented them into 20–70 single cells, and then applied our method to analyze the results. We found that the values of graph transition energy measured against DMSO for both muricin A and squamocin decrease over time as shown in Figure 9. The lower the energy the higher the morphological difference compared with the DMSO groups. That is, from Figure 9, we can conclude that squamocin induces more fragmented mitochondria than muricin A in CHO cells. This conclusion can be verified by human inspection. Figure 10 shows the representative images in each group. We can clearly see that on day 2, mitochondria of squamocin-treated cells have large round shape, indicating a high degree of fragmentation, while on the same day, mitochondria of muricin A-treated cells are much smaller. Our approach allows us to automatically quantify to what extent the difference is. This is crucial for large scale studies where human inspection becomes increasingly time-consuming and error-prone.

Fig. 9.
Comparing the effect of mitochondrial fragmentation by squamocin and muricin A in CHO cells.
Fig. 10.
Representative images of the cells after treated by DMSO (a), muricin A (b) and squamocin (c) for a certain number of days.

Again, we can compute graph transition energy between each pair of image groups of muricin A and squamocin and apply MDS to plot all image groups in a 2D space as shown in Figure 11. Arrows in the plot indicate the time course of the morphological change by these compounds. The plot clearly shows a drastic change of mitochondrial morphology from day 1 to day 2 in squamocin-treated cells.

Fig. 11.
Morphological similarity of mitochondria in different image groups by MDS according to the pairwise graph transition energy values.

5 CONCLUSION

In the end, we summarize our contributions as follows:

  • We defined graph transition energy to quantify morphological difference of two cell-image collections under different perturbations.
  • We applied a spectral graph theoretic regularization to re-weight the feature space according to training examples of extremely different cases. Experimental results show that our feature space transformation method actually improves the quality of classification and quantification. The method also allows for calibration by providing a confidence score of the quantification.
  • We demonstrated how to quantify morphological difference of collections of cell images in our approach and showed that our approach has a higher discrimination power than a competing approach by Loo et al. (2007).
  • We showed that after applying our feature space transformation method, the performance of a semi-supervised learning method by Zhu et al. (2003) for subcellular localization can be improved significantly.
  • Coupling with a k-NN regression method, we demonstrated that our transformation method improves the accuracy of the estimation of partial fragmentation of mitochondria in cell images, and that our method is more effective than NCA.
  • In Section 4.3, we demonstrated a case study of two Annonaceous acetogenin compounds and their effects to mitochondrial fragmentation.

Our approach is not specific to any particular morphological differences and can potentially be applied to quantify any difference with an appropriate set of features. Therefore, it has the potential to rapidly screen many drug candidates and understand their effects. In our future work, we will apply this approach to continue our study of Annonaceous acetogenins and their various effects to mitochondria.

Supplementary Material

[Supplementary Data]

ACKNOWLEDGEMENTS

We wish to thank Dr F-R. Chang and Dr Y-C. Wu of Kaohsiung Medical University for providing muricin A and squamocin compounds, C-G. Yeh for analyzing 350 images with aid of Mito-Q, and Y-S. Wong for acquiring all the cell images.

Funding: (Advanced Bioinformatics Core); NSC97-3112-B-001-024. National Research Program in Genomic Medicine (NRPGM); National Science Council, Taiwan.

Conflict of interest: none declared.

REFERENCES

  • Azran A, Ghahramani Z. ICML '06: Proceedings of the 23rd international conference on Machine learning. New York, NY, USA: ACM; 2006. A new approach to data driven clustering; pp. 57–64.
  • Belkin M, Niyogi P. Laplacian eigenmaps for dimensionality reduction and data representation. Neur. Comput. 2003a;15:1373–1396.
  • Belkin M, Niyogi P. Advances in Neural Information Processing Systems 15 (NIPS*2002). Cambridge, MA: MIT Press; 2003b. Using manifold stucture for partially labeled classification; pp. 929–936.
  • Borgwardt KM, et al. Protein function prediction via graph kernels. Bioinformatics. 2005;21(suppl. 1):47–56. [PubMed]
  • Brooks C, et al. Bak regulates mitochondrial morphology and pathology during apoptosis by interacting with mitofusins. Proc. Natl Acad. Sci. USA. 2007;104:11649–11654. [PubMed]
  • Butman M, Goldberger J. Face recognition using classification based linear projections. EURASIP J. Adv. Signal Proc. 2008;8:1–7.
  • Carpenter AE. Image-based chemical screening. Nat. Chem. Biol. 2007;3:461–465. [PubMed]
  • Chang H.-H, Moura J.MF. Proceedings of IEEE International Conference on Image Processing (ICIP-2007). II. Los Alamitos, CA: IEEE Publications Office; 2007. Classification by Cheeger constant regularization; pp. II209–II212.
  • Chang H.-H, et al. Automatic detection of regional heart rejection in USPIO-enhanced mri. IEEE Trans. Med. Imag. 2008;27:1095–1106. [PMC free article] [PubMed]
  • Chen L.-F, et al. A new LDA-based face recognition system which can solve the small sample size problem. Pattern Recogn. 2000;33:1713–1726.
  • Chung F.RK. CBMS Regional Conference Series in Mathematics. Providence, RI: American Mathematical Society; 1997. Spectral Graph Theory.
  • Desagher S, Martinou J.-C. Mitochondria as the central control point of apoptosis. Trends Cell Biol. 2000;10:369–377. [PubMed]
  • Frank S, et al. The role of dynamin-related protein 1, a mediator of mitochondrial fission, in apoptosis. Dev. Cell. 2001;1:515–525. [PubMed]
  • Goldberger J, et al. Advances in Neural Information Processing Systems 17 (NIPS*2004). Cambridge, MA: MIT Press; 2005. Neighbourhood components analysis; pp. 513–520.
  • Gonzalez RC, Woods RE. Digital Image Processing. 2. Upper Saddle River, NJ: Prentice Hall; 2002.
  • Haralick RM. Statistical and structural approaches to texture. Proceedings of IEEE. 1979;67:786–804.
  • Huang K, Murphy RF. Boosting accuracy of automated classification of fluorescence microscope images for location proteomics. BMC Bioinformatics. 2004;5:78. [PMC free article] [PubMed]
  • Isfahani SN, et al. International Conference on BioMedical Engineering and Informatics. Washington, DC: IEEE Computer Society; 2008. A new approach for touching cells segmentation; pp. 816–820.
  • Jeong S.-Y, Seol D.-W. The role of mitochondria in apoptosis. BMB Rep. 2008;41:11–22. [PubMed]
  • Jones TR, et al. Scoring diverse cellular morphologies in image-based screens with iterative feedback and machine learning. Proc. Natl Acad. Sci. USA. 2009;106:1826–1831. [PubMed]
  • Karbowski M, Youle RJ. Dynamics of mitochondrial morphology in healthy cells and during apoptosis. Cell Death Differ. 2003;10:870–880. [PubMed]
  • Lang P, et al. Cellular imaging in drug discovery. Nat. Rev. Drug Discove. 2006;5:343–356. [PubMed]
  • Lee Y.-J, et al. Roles of the mammalian mitochondrial fission and fusion mediators Fis1, Drp1, and Opa1 in apoptosis. Mol. Biol. Cell. 2004;15:5001–5011. [PMC free article] [PubMed]
  • Lin C.-C, et al. Boosting multiclass learning with repeating codes and weak detectors for protein subcellular localization. Bioinformatics. 2007;23:3374–3381. [PubMed]
  • Loo L.-H, et al. Image-based multivariate profiling of drug responses from single cells. Nat. Methods. 2007;4:445–453. [PubMed]
  • Roweis ST, Saul LK. Nonlinear dimensionality reduction by locally linear embedding. Science. 2000;290:2323–2326. [PubMed]
  • Serra J. Image Analysis and Mathematical Morphology. Orlando, FL: Academic Press; 1998.
  • Taguchi N, et al. Mitotic phosphorylation of dynamin-related gtpase drp1 participates in mitochondrial fission. J. Biol. Chem. 2007;282:11521–11529. [PubMed]
  • Tamai S, et al. Characterization of the mitochondrial protein letm1, which maintains the mitochondrial tubular shapes and interacts with the aaa-atpase bcs1l. J. Cell Sci. 2008;121:2588–2600. [PubMed]
  • Tanaka M, et al. An unbiased cell morphology-based screen for new, biologically active small molecules. PLoS biology. 2005;3 [PMC free article] [PubMed]
  • Tenenbaum JB, et al. A global geometric framework for nonlinear dimensionality reduction. Science. 2000;290:2319–2323. [PubMed]
  • Weizman L, Goldberger J. Geoscience and Remote Sensing Symposium, 2007. IGARSS 2007. IEEE International. Los Alamitos, CA: IEEE Publications Office; 2007. A classification-based linear projection of labeled hyperspectral data; pp. 3202–3205.
  • Yaffe MP. Dynamic mitochondria. Nat. Cell Biol. 1999;1:149–150. [PubMed]
  • Zhan D.-C, et al. ICML '09: Proceedings of the 26th Annual International Conference on Machine Learning. New York, NY, USA: ACM; 2009. Learning instance specific distances using metric propagation; pp. 1225–1232.
  • Zhou X, Wong ST. Informatics challenges of high-throughput microscopy. IEEE Signal Proc. Mag. 2006;23:63–72.
  • Zhu X, Ghahramani Z. Technical Report CMU-CALD-02-107. Pittsburg, PA: Carnegie Mellon University; 2002. Learning from labeled and unlabeled data with label propagation.
  • Zhu X, et al. ICML '03: Proceeding of the 12th International Conference on Machine Learning. Menlo Park, CA: AAAI Press; 2003. Semi-supervised learning using gaussian fields and harmonic functions; pp. 912–919.

Articles from Bioinformatics are provided here courtesy of Oxford University Press