PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of bioinfoLink to Publisher's site
 
Bioinformatics. 2013 July 1; 29(13): i274–i282.
Published online 2013 June 19. doi:  10.1093/bioinformatics/btt225
PMCID: PMC3694676

Automated target segmentation and real space fast alignment methods for high-throughput classification and averaging of crowded cryo-electron subtomograms

Abstract

Motivation: Cryo-electron tomography allows the imaging of macromolecular complexes in near living conditions. To enhance the nominal resolution of a structure it is necessary to align and average individual subtomograms each containing identical complexes. However, if the sample of complexes is heterogeneous, it is necessary to first classify subtomograms into groups of identical complexes. This task becomes challenging when tomograms contain mixtures of unknown complexes extracted from a crowded environment. Two main challenges must be overcomed: First, classification of subtomograms must be performed without knowledge of template structures. However, most alignment methods are too slow to perform reference-free classification of a large number of (e.g. tens of thousands) of subtomograms. Second, subtomograms extracted from crowded cellular environments, contain often fragments of other structures besides the target complex. However, alignment methods generally assume that each subtomogram only contains one complex. Automatic methods are needed to identify the target complexes in a subtomogram even when its shape is unknown.

Results: In this article, we propose an automatic and systematic method for the isolation and masking of target complexes in subtomograms extracted from crowded environments. Moreover, we also propose a fast alignment method using fast rotational matching in real space. Our experiments show that, compared with our previously proposed fast alignment method in reciprocal space, our new method significantly improves the alignment accuracy for highly distorted and especially crowded subtomograms. Such improvements are important for achieving successful and unbiased high-throughput reference-free structural classification of complexes inside whole-cell tomograms.

Contact: alber/at/usc.edu

Supplementary information: Supplementary data are available at Bioinformatics online.

1 INTRODUCTION

Cryo-electron tomography enables the 3D imaging of macromolecular complexes at nanometer-scale resolution in near native conditions (Best et al., 2007; Lucic et al., 2005; Medalia et al., 2002). Tomograms of individual cells are essentially 3D representations of the entire proteome providing a snapshot of the distributions of protein complexes (Beck et al., 2011). However, comprehensive detection of individual complexes in cell tomograms is challenging because of the inherent low signal-to-noise ratio (SNR), missing data, nonisotropic resolution and the fact that individual macromolecules are difficult to recognize in a highly crowded environment (Beck et al., 2009; Best et al., 2007; Böhm et al., 2000; Frangakis et al., 2002; Medalia et al., 2002; Nickell et al., 2006).

Most methods for detecting complexes in cell tomograms rely on a template structure, which serves as a reference in searching for a similar pattern in the tomogram (e.g. Beck et al., 2009). However, for an unbiased detection and classification of all the different complexes in a cellular tomogram, template-free methods are needed (Xu et al., 2011). Such an analysis is challenging and relies on three main steps: First, the locations of potential complexes are detected by using particle-picking methods and a subregion surrounding the potential complex is extracted (i.e. the subtomogram). Second, the subtomogram regions corresponding only to the target complex must be detected, which allows masking out all background regions, which in turn contain noise and in case of crowded subtomograms also fragments of surrounding structures. Third, reference-free structural classification of the masked subtomograms is performed, which is typically based on an iterative process of subtomogram alignments, classifications and averaging (e.g. Amat et al., 2010; Bartesaghi et al., 2008; Chen et al., 2013; Förster et al., 2008; Volkmann, 2010; Xu et al., 2012). Finally, the averaging of aligned subtomograms of the same complexes will enhance the nominal resolution of the resulting density maps.

In this article, we address two main challenges in reference-free subtomogram classifications. First, we propose an automatic method for adaptive masking of target regions in crowded subtomograms without the knowledge of the shape of the complex. Such a method is particularly important for subtomograms extracted from crowded environments, such as the cell cytoplasm. In such a case, subtomograms will also contain fragmental regions of other complexes owing to the high particle density in the tomogram. Masking out these regions is of great importance in the subsequent classification process. Unlike automatic masking methods based on voxel weighting (Xu et al., 2012) or dimension reduction methods such as Principle Component Analysis (Heumann et al., 2011), our method is highly scalable because it is independent of the classification process and does not involve iterative processing of a large number of subtomograms. Second, we propose a new method for fast rotational matching of subtomograms in real space. The 3D subtomogram alignment is computationally the most intensive step in the classification process. Currently, most alignment methods are based on maximizing the overlap similarity of two subtomograms through exhaustive search over all rigid transforms (rotation and translation) of one subtomogram with respect to the other (e.g. Amat et al., 2010; Förster et al., 2008). Such methods scan exhaustively through each rotational angle, and find the best corresponding translation using Fast Fourier Transform (FFT) (e.g. Amat et al., 2010; Förster et al., 2008). These methods are computationally intensive, which limits their applicability in reference-free classifications. To increase alignment accuracy, local search methods have been developed to refine initial alignments (e.g. Bartesaghi et al., 2008; Hrabe et al., 2012; Xu and Alber, 2012). To increase computational efficiencies, two approximate alignment methods have been developed, which separate the translational from the rotational subtomogram alignments (Bartesaghi et al., 2008; Xu et al., 2012). These methods are based on similarity scores defined in the reciprocal space (i.e. Fourier space). To separate translational and rotational search, these methods use an approximate score that introduces translational invariance by eliminating the phase information of the complex coefficients in the reciprocal space and use only the magnitude of the Fourier coefficients (Bartesaghi et al., 2008; Xu et al., 2012). By expressing the structural information in Spherical Harmonics (SH) Expansion, it is possible to formulate the rotational alignment as a fast rotational matching, which simultaneously calculates alignment scores of all rotations using FFT. This procedure significantly increases the alignment speed. For example, our previously published method (Xu et al., 2012) achieved between hundreds to thousands fold speedup (depending on subtomograms size and rotational angle interval) compared with the standard scanning-based alignment method (Förster et al., 2008) by using a translational invariance approximation of a similarity score equivalent to a popular score proposed by (Förster et al., 2008).

However, formulating the fast rotational matching in reciprocal space has some limitations. In reciprocal space, the majority of informative signals of a complex is usually contained in a relatively small amount of high magnitude Fourier coefficients (Amat et al., 2010). These coefficients are often occupying relatively few voxels in the reciprocal space located within a small region centered at the origin. Therefore, the SH expansion in the fast rotational matching may be hampered by interpolation errors, which reduces the accuracy in the alignment process. In contrast, in real space, the signal covers a wider area within a subtomogram, and the SH expansion is expected to be more accurate. It is therefore beneficial to formulate fast rotational matching in real space. In addition, real space fast rotational matching uses the full constrained alignment score instead of a translation invariant approximation when the method is expressed in reciprocal space (Bartesaghi et al., 2008; Xu et al., 2012).

Here, we devise an improved alignment method, which formulates the fast rotational matching of subtomograms in real space. To achieve this goal, we first identify a keypoint in the target complex whose relative location is invariant to rigid transformations and serves as the center of rotation in the fast rotational matching of the two subtomograms. A natural choice is the center of mass of a complex. However, detection of the center of mass is not trivial and cannot be approximated by the geometrical or mass density center of the subtomogram. The reason for this complication is that tomograms are images of the crowded cellular environment (Beck et al., 2009), and a subtomogram contains not necessarily mass density of only one single complex. When a potential complex is detected and its subtomogram extracted as a rectangular cube then this subtomogram typically contains also fragments of other surrounding structures that occupy parts of the subtomogram. Then, the center of mass estimation and subtomogram alignment are affected by the existence of the additional mass density in the subtomogram. It is necessary to focus only on the regions that contain the actual target complex. Even when a subtomogram does not contain any surrounding structures, its mean intensity is often close to background intensity owing to the suppression of low frequency signal owing to the Contrast Transfer Function (CTF) effect, which also makes the estimation of the center of mass of the complex inaccurate.

Applying our automatic target complex segmentation method, it is possible to estimate a center of mass based only on the target complex regions in the subtomogram (i.e. the constrained center of mass). Once the constrained center of mass is detected, we design a fast subtomogram alignment method that uses real space signals and takes into account missing wedge corrections using reciprocal space signals. Moreover, the alignment focuses exclusively on the target complex regions in each subtomogram, which significantly decreases the influence of background noise and surrounding structures to the alignment process.

Our experiments show that our new approach significantly increases the alignment accuracy compared with our previous proposed fast alignment method (Xu et al., 2012) for highly distorted subtomograms. Most importantly, the new method is highly robust to the influence of surrounding structures inside crowded subtomograms.

2 METHODS

This section begins by describing the automatic segmentation of the target complex and then focuses on the fast rotational alignment in real space.

2.1 Automatic segmentation of the target complex in a subtomogram

The automatic target complex segmentation method consists of several steps (Figs 1, ,22 and and7):7): (i) automatic scale selection is performed to maximally enhance the detection of structural boundaries in a subtomogram; (ii) segments are defined that describe the boundaries between structural elements; and (iii) a classification of the detected structural segments into candidate regions is performed by using a statistical model based clustering.

Fig. 1.
Flow chart for segmentation of target complex regions and the real space alignment of subtomogram using only target complex regions
Fig. 2.
Four typical examples for target complex segmentation based on simulated data. Each row is an example, whereas each subfigure is a slice through the x-z plane in the 3D image. (A) Ground truth. (B) Simulated subtomograms at SNR 0.005 and tilt angle range ...
Fig. 7.
Segmentation of target complex regions in Leptospira interrogans subtomograms. Each row shows the segmentation process for a ribosome instance. Each subfigure in a row consists of the consecutive slices through the y-axis of the 3D image. (A) Original ...

2.1.1 Automatic scale selection for optimal smoothing and boundary enhancement

In this article, we assume that in a subtomogram high image intensity corresponds to high electron density. When a subtomogram is renormalized so that its mean intensity is zero, the boundary of a complex tends to have a negative intensity surface (Fig. 3, yellow arrow). This property is due to the CTF effects in the electron tomography imaging process, which suppresses low frequency signals. We use this negative intensity surface as a main characteristic when identifying the segments at the boundary of structural elements, such as the target complex and those structures captured at the outer regions of the subtomogram. To reduce the influence of noise, smoothing of the subtmograms is needed. However, too much smoothing will weaken the negative intensity boundary. To find the optimal degree of smoothing, we use the scale-space representation method and propose an automatic scale selection method to find the optimal degree of smoothing that maximally enhances the boundary surface while minimizing the noise. Scale-space (Witkin, 1983) is a framework for multi-scale signal representation. An image is represented as a one-parameter family of smoothed images, the scale-space representation, parameterized by the size of the Gaussian smoothing kernel (Lindeberg, 1994):

equation image

Given a subtomogram, f, a Gaussian smoothing can be expressed as the convolution between f and An external file that holds a picture, illustration, etc.
Object name is btt225i4.jpg:

equation image

In the smoothed subtomogram An external file that holds a picture, illustration, etc.
Object name is btt225i5.jpg, all features with a size smaller than σ are filtered out, whereas others are preserved. Our objective is to determine the optimal scale σ to enhance the formation of a negative intensity surface while minimizing the influence of noise. To do so, we scan a range of feasible σ values. For a given smoothed subtomogram with scale σ, we find all local minima with negative intensities in An external file that holds a picture, illustration, etc.
Object name is btt225i6.jpg and denote this set as An external file that holds a picture, illustration, etc.
Object name is btt225i7.jpg. We then define a local surface score for each minima in An external file that holds a picture, illustration, etc.
Object name is btt225i8.jpg, to identify those minima that are likely part of a negative intensity surface. The surface score is calculated as follows: we calculate the Hessian matrix An external file that holds a picture, illustration, etc.
Object name is btt225i9.jpg, whose elements are second-order derivatives of An external file that holds a picture, illustration, etc.
Object name is btt225i10.jpg evaluated at An external file that holds a picture, illustration, etc.
Object name is btt225i11.jpg,

equation image

where An external file that holds a picture, illustration, etc.
Object name is btt225i12.jpg is the corresponding location of a voxel k.

Fig. 3.
A subtomogram’s image features. The boundary feature is a surface that has strong negative intensity (yellow arrow). A local maxima in An external file that holds a picture, illustration, etc.
Object name is btt225i2.jpg is highlighted by a red cross, which is inside the macromolecular complex. Its corresponding structural boundary ...

For any An external file that holds a picture, illustration, etc.
Object name is btt225i13.jpg, let An external file that holds a picture, illustration, etc.
Object name is btt225i14.jpg be the ordered eigenvalues of An external file that holds a picture, illustration, etc.
Object name is btt225i15.jpg such that An external file that holds a picture, illustration, etc.
Object name is btt225i16.jpg. Then, we can construct a local surface score An external file that holds a picture, illustration, etc.
Object name is btt225i17.jpg similar to (Martinez-Sanchez et al., 2011) as follows:

equation image

We find the σ value that maximally enhances the surface scores of all the minima so that the most negative minima tends to have a strong surface score. To do so, we calculate an accumulative surface score from all the minima. First, the minima are ordered in ascending order in An external file that holds a picture, illustration, etc.
Object name is btt225i18.jpg such that

equation image

where An external file that holds a picture, illustration, etc.
Object name is btt225i19.jpg. Then the accumulative surface score An external file that holds a picture, illustration, etc.
Object name is btt225i20.jpg is defined as

equation image

Because An external file that holds a picture, illustration, etc.
Object name is btt225i21.jpg is in ascending order with respect to i, this score An external file that holds a picture, illustration, etc.
Object name is btt225i22.jpg forms a monotonically increasing function of An external file that holds a picture, illustration, etc.
Object name is btt225i23.jpg. We can measure the Area Under Curve, An external file that holds a picture, illustration, etc.
Object name is btt225i24.jpg, using An external file that holds a picture, illustration, etc.
Object name is btt225i25.jpg and An external file that holds a picture, illustration, etc.
Object name is btt225i26.jpg. An external file that holds a picture, illustration, etc.
Object name is btt225i27.jpg tends to be large when the most negative minima have large surface scores. To measure how much the AUC reflects true surface signals rather than random noise, we also randomly permutate all voxels in f and repeat the aforementioned procedure to calculate a permutated AUC, An external file that holds a picture, illustration, etc.
Object name is btt225i28.jpg.

To determine the optimal Gaussian filtering value σ, we scan through a range of different σs to find the An external file that holds a picture, illustration, etc.
Object name is btt225i29.jpg value that maximizes

equation image

For simplicity, in the following sections, we denote An external file that holds a picture, illustration, etc.
Object name is btt225i30.jpg (Figs 2C and and7B)7B) and also denote An external file that holds a picture, illustration, etc.
Object name is btt225i31.jpg for any voxel i.

2.1.2 Detecting structural, boundary and background segments in the subtomogram

After identifying the optimal scale An external file that holds a picture, illustration, etc.
Object name is btt225i32.jpg, we now determine those local minima that have strong surface scores and are likely to separate different structural elements (i.e. the target complex and surrounding structures). We define them as those local minima with both high surface scores and low intensity values. To do so, first we select the minima whose intensity is significantly smaller than the average intensity, i.e.

equation image

where An external file that holds a picture, illustration, etc.
Object name is btt225i33.jpg is a constant parameter to control the significance level, and An external file that holds a picture, illustration, etc.
Object name is btt225i34.jpg is the Median Absolute Deviation, which is a robust measure of variability.

In addition, we also collect those minima whose surface score is significantly larger than the average scores of all minima, i.e.

equation image

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

All minima with strong surface scores are defined as the intersection of both groups.

equation image

We denote these minima as the boundary surface minima.

We then determine all local intensity maxima An external file that holds a picture, illustration, etc.
Object name is btt225i36.jpg in An external file that holds a picture, illustration, etc.
Object name is btt225i37.jpg, which are likely to be inside structural elements. We first perform Watershed segmentation (Beucher and Lantuéjoul, 1979; Volkmann, 2002) using An external file that holds a picture, illustration, etc.
Object name is btt225i38.jpg as seeds (Figs 2D and and7C).7C). Each segment corresponds to one local maximum, and therefore the terms local maximum and segment are used interchangeably. We then select all the segments that border with local minima defined in An external file that holds a picture, illustration, etc.
Object name is btt225i39.jpg and denote the corresponding set of maxima as An external file that holds a picture, illustration, etc.
Object name is btt225i40.jpg. The actual separations between structural elements and background regions are always at the boundary between segments in An external file that holds a picture, illustration, etc.
Object name is btt225i41.jpg. The segments in An external file that holds a picture, illustration, etc.
Object name is btt225i42.jpg can be divided into two types: a structural boundary segment is located in the target complex or the surrounding structures (Fig. 3, inside red circle), and it contains a maximum with high intensity value (Fig. 3, highlighted by red cross); a background boundary segment is located in the regions that do not contain structural elements (Fig. 3, inside green circle), and which contain a maximum with low intensity value (Fig. 3, highlighted by green cross).

To separate these two types of boundary segments, we perform k-means clustering on An external file that holds a picture, illustration, etc.
Object name is btt225i43.jpg (with cluster number An external file that holds a picture, illustration, etc.
Object name is btt225i44.jpg), resulting in the two groups An external file that holds a picture, illustration, etc.
Object name is btt225i45.jpg and An external file that holds a picture, illustration, etc.
Object name is btt225i46.jpg. An external file that holds a picture, illustration, etc.
Object name is btt225i47.jpg are segments containing local maxima with higher intensity values and are defined as structural boundary segments (Figs 2E and and7D),7D), whereas An external file that holds a picture, illustration, etc.
Object name is btt225i48.jpg are defined as background boundary segments.

After focusing on boundary segments (i.e. those that border with surface minima), we now characterize the remaining segments as either part of structural elements or background regions. To do this, we first determine a cutoff intensity value An external file that holds a picture, illustration, etc.
Object name is btt225i49.jpg as the largest intensity value in all the background boundary segments An external file that holds a picture, illustration, etc.
Object name is btt225i50.jpg. All remaining segments that contain a local maxima with an intensity value smaller that An external file that holds a picture, illustration, etc.
Object name is btt225i51.jpg are assigned as background regions.

Finally, the remaining segments whose maxima have intensity values larger than An external file that holds a picture, illustration, etc.
Object name is btt225i52.jpg are defined as structural segments and are therefore assumed to be either part of the target complex or other structural elements, i.e. An external file that holds a picture, illustration, etc.
Object name is btt225i53.jpg.

2.1.3 Combining structural segments into complexes using model based clustering

We now cluster structural segments in An external file that holds a picture, illustration, etc.
Object name is btt225i54.jpg according to the intensity values and location of their local maxima so that each cluster of segments will approximately correspond to one consecutive structural region. We choose a model-based clustering approach because it allows an automatic determination of the optimal set of clusters. To perform the clustering, we assume that the location and intensity of local maxima from clusters follows a probabilistic mixture model, where each component probability distribution corresponds to a cluster (Fraley and Raftery, 2002). We model the probability of independently observing both location and intensity of a local maxima An external file that holds a picture, illustration, etc.
Object name is btt225i55.jpg from a cluster k as:

equation image
(1)

where An external file that holds a picture, illustration, etc.
Object name is btt225i56.jpg is a probability proportional to the intensity An external file that holds a picture, illustration, etc.
Object name is btt225i57.jpg, and An external file that holds a picture, illustration, etc.
Object name is btt225i58.jpg is a probability density function in the form of a multivariate normal distribution with mean An external file that holds a picture, illustration, etc.
Object name is btt225i59.jpg and covariance An external file that holds a picture, illustration, etc.
Object name is btt225i60.jpg. For simplification, we assume that the clusters are spherical with different sizes, i.e. An external file that holds a picture, illustration, etc.
Object name is btt225i61.jpg. The main difference between our model and the standard model-based clustering (Fraley and Raftery, 2002) is the inclusion of the image intensity related term An external file that holds a picture, illustration, etc.
Object name is btt225i62.jpg in Equation (1). Our simulation experiments show that the inclusion of image intensity improves the clustering performance (Supplementary Document).

The likelihood for the mixture model with K clusters is:

equation image

where x represents the collection of observed local maxima, and An external file that holds a picture, illustration, etc.
Object name is btt225i63.jpg is the probability of cluster k (i.e. cluster mixing probability). A model-based clustering determines the parameters An external file that holds a picture, illustration, etc.
Object name is btt225i64.jpg, An external file that holds a picture, illustration, etc.
Object name is btt225i65.jpg, An external file that holds a picture, illustration, etc.
Object name is btt225i66.jpg that maximizes An external file that holds a picture, illustration, etc.
Object name is btt225i67.jpg.

We use a scale-space representation to estimate cluster centers An external file that holds a picture, illustration, etc.
Object name is btt225i68.jpg. Given a scale An external file that holds a picture, illustration, etc.
Object name is btt225i69.jpg, we obtain the set An external file that holds a picture, illustration, etc.
Object name is btt225i70.jpg of local maxima with positive intensities in An external file that holds a picture, illustration, etc.
Object name is btt225i71.jpg. The local maxima An external file that holds a picture, illustration, etc.
Object name is btt225i72.jpg are ordered with decreasing intensity values, i.e. An external file that holds a picture, illustration, etc.
Object name is btt225i73.jpg, where An external file that holds a picture, illustration, etc.
Object name is btt225i74.jpg is the corresponding location of local maxima jk on An external file that holds a picture, illustration, etc.
Object name is btt225i75.jpg. We choose cluster centers as the first K local maxima, i.e An external file that holds a picture, illustration, etc.
Object name is btt225i76.jpg

Similar to the standard model based clustering (Celeux and Govaert, 1995; Fraley and Raftery, 2002), an Expectation Maximization algorithm is used to estimate An external file that holds a picture, illustration, etc.
Object name is btt225i77.jpg and An external file that holds a picture, illustration, etc.
Object name is btt225i78.jpg. The expectation maximization is defined by two iterative steps:

E-step: update an estimated posterior probability An external file that holds a picture, illustration, etc.
Object name is btt225i79.jpg that the local maximum i belongs to cluster k:

equation image

M-step: update the cluster mixing probability An external file that holds a picture, illustration, etc.
Object name is btt225i80.jpg, and variances An external file that holds a picture, illustration, etc.
Object name is btt225i81.jpg, given the new parameters in the previous step:

equation image

where

equation image

and

equation image

where

equation image

The Expectation Maximization algorithm iterates until there is no significant change in the likelihood An external file that holds a picture, illustration, etc.
Object name is btt225i82.jpg. The clustering of local maxima leads to clusters of the corresponding segments (Figs 2F and and7E).7E). Using statistical modeling for clustering enables the determination of an optimal set of clusters through model selection. By varying the scale σ and cluster number K, we obtain the optimal values for scale and cluster set using the Bayesian Information Criterion (Fraley and Raftery, 2002):

equation image

The final target complex is then represented by combining multiple segments that fulfill the following criteria: (i) they belong to same cluster in the model based clustering; (ii) they are in consecutive contact to each other; and (iii) the maximum intensity values An external file that holds a picture, illustration, etc.
Object name is btt225i83.jpg of the voxels between two consecutive segments i and j is larger than a threshold. The threshold is computed as An external file that holds a picture, illustration, etc.
Object name is btt225i84.jpg, where An external file that holds a picture, illustration, etc.
Object name is btt225i85.jpg is a constant parameter.

The target complex region is then selected as the combined segments that contain at least one structural boundary segment (defined by An external file that holds a picture, illustration, etc.
Object name is btt225i86.jpg) and are located closest to the center of the subtomogram. To prevent that small features of the target complex are excluded, a dilation operation is performed to define the target complex region A (Figs 2G and and77F).

2.2 Real space fast subtomogram alignment using the target complex region and its center of mass

Fast rotational matching in real space relies on an initial alignment of the center of masses of the two subtomograms. Therefore, we describe first how the center of mass of each target complex is calculated before we describe the fast rotational matching approach.

2.2.1 Calculation of center of mass for target complex region

The determination of the center of mass in subtomograms is not trivial, primarily owing to high noise levels, the suppression of low frequency signals in the reciprocal space and the existence of surrounding structures inside a subtomogram.

Here, the center of mass is calculated only based on the regions in the target complex A. We set a thresholded region An external file that holds a picture, illustration, etc.
Object name is btt225i87.jpg on the filtered subtomogram An external file that holds a picture, illustration, etc.
Object name is btt225i88.jpg, i.e. An external file that holds a picture, illustration, etc.
Object name is btt225i89.jpg. Then, we use An external file that holds a picture, illustration, etc.
Object name is btt225i90.jpg to calculate the center of mass. We define another intensity function

equation image

We then calculate the constrained center of mass An external file that holds a picture, illustration, etc.
Object name is btt225i91.jpg of An external file that holds a picture, illustration, etc.
Object name is btt225i92.jpg.

2.2.2 Fast rotation alignment

We represent two subtomograms f and g as two integrable functions An external file that holds a picture, illustration, etc.
Object name is btt225i93.jpg. They have been transformed such that (i) their mean intensities are adjusted to zero, (ii) the voxel intensities outside the target complex region A is set to zero, and (iii) they are translated according to An external file that holds a picture, illustration, etc.
Object name is btt225i94.jpg so that their constrained centers of mass are aligned. The constrained center of mass is used as rotational centers in our fast rotational alignment. To correct for the missing wedge distortions, we only keep the Fourier coefficients within the overlap of nonmissing wedge regions of both subtomograms. Following Bartesaghi et al. (2008), Förster et al. (2008) and Xu et al. (2012), we use a binary missing wedge mask function as An external file that holds a picture, illustration, etc.
Object name is btt225i95.jpg, which defines valid and missing Fourier coefficients in reciprocal space. For example, with a tilt angle range An external file that holds a picture, illustration, etc.
Object name is btt225i96.jpg, the missing wedge mask function can be defined as An external file that holds a picture, illustration, etc.
Object name is btt225i97.jpg, where An external file that holds a picture, illustration, etc.
Object name is btt225i98.jpg is the indicator function. Given An external file that holds a picture, illustration, etc.
Object name is btt225i99.jpg, the real space subtomogram that excludes any coefficients inside of any of the two missing wedge regions is defined as,

equation image

equation image

where An external file that holds a picture, illustration, etc.
Object name is btt225i100.jpg denotes the real part of a complex function; An external file that holds a picture, illustration, etc.
Object name is btt225i101.jpg is the Fourier transform operator; and An external file that holds a picture, illustration, etc.
Object name is btt225i102.jpg and An external file that holds a picture, illustration, etc.
Object name is btt225i103.jpg are missing wedge masks for f and g, respectively.

Given a 3D rotation R, we can calculate a correlation

equation image
(2)

where An external file that holds a picture, illustration, etc.
Object name is btt225i104.jpg is the rotation operator such that An external file that holds a picture, illustration, etc.
Object name is btt225i105.jpg for any function An external file that holds a picture, illustration, etc.
Object name is btt225i106.jpg, and An external file that holds a picture, illustration, etc.
Object name is btt225i107.jpg is the rotation matrix corresponding to 3D rotation R. As in Xu et al. (2012), the correlation measure An external file that holds a picture, illustration, etc.
Object name is btt225i108.jpg in Equation (2) is equivalent (up to a constant factor.) to a popular constrained correlation score with missing wedge correction (Förster et al., 2008).

The fast rotational alignment is based on sampling of An external file that holds a picture, illustration, etc.
Object name is btt225i109.jpg simultaneously over all rotations R. To do so, we apply a fast 3D volumetric rotational matching (Kovacs and Wriggers, 2002; Xu et al., 2012). It can be seen that Equation (2) can be formulated as being composed of three rotational correlation functions of the form An external file that holds a picture, illustration, etc.
Object name is btt225i110.jpg, where p and q are component functions; An external file that holds a picture, illustration, etc.
Object name is btt225i111.jpg denotes the complex conjugate (when e is a real valued function, the complex conjugate An external file that holds a picture, illustration, etc.
Object name is btt225i112.jpg.) of a function e. Specifically, we can represent An external file that holds a picture, illustration, etc.
Object name is btt225i113.jpg as

equation image
(3)

where An external file that holds a picture, illustration, etc.
Object name is btt225i114.jpg, An external file that holds a picture, illustration, etc.
Object name is btt225i115.jpg, An external file that holds a picture, illustration, etc.
Object name is btt225i116.jpg, An external file that holds a picture, illustration, etc.
Object name is btt225i117.jpg, An external file that holds a picture, illustration, etc.
Object name is btt225i118.jpg, and An external file that holds a picture, illustration, etc.
Object name is btt225i119.jpg.

When represented in spherical coordinates, the p and q components of these three functions An external file that holds a picture, illustration, etc.
Object name is btt225i120.jpg can be approximated by an SH expansion. Following Garzón et al. (2007), Kovacs and Wriggers (2002) and Xu et al., (2012), the p and q components in Equation (3) can be expressed as:

equation image

where B is the bandwidth, and An external file that holds a picture, illustration, etc.
Object name is btt225i121.jpg are the coefficients associated with the complex-valued spherical harmonic function An external file that holds a picture, illustration, etc.
Object name is btt225i122.jpg with degree l and order m, where r, β and λ are the radial distance, co-latitude and longitude, respectively, which define the position in spherical coordinate.

When a suitable parameterization of the 3D rotational group is achieved, the rotational correlation function An external file that holds a picture, illustration, etc.
Object name is btt225i123.jpg of all sampled rotations R can be represented as an inverse FFT of a 3D array of the sum of integrals as follows (Garzón et al., 2007; Kovacs and Wriggers, 2002; Xu et al., 2012):

equation image

where An external file that holds a picture, illustration, etc.
Object name is btt225i124.jpg are real coefficients defining the elements of the Wigner small d-matrix evaluated at An external file that holds a picture, illustration, etc.
Object name is btt225i125.jpg.

As a consequence of the aforementioned formulation, the cross-correlation functions can be efficiently sampled by using FFT simultaneously over all sampled rotations (Kovacs and Wriggers, 2002; Xu et al., 2012). The sampling is given as twice the bandwidth B used in the SH transformations. Therefore, An external file that holds a picture, illustration, etc.
Object name is btt225i126.jpg can be efficiently computed over all sampled rotations R. Similar to Bartesaghi et al. (2008) and Xu et al. (2012), the set of candidate rotations are then obtained by identifying the local maxima of ρ with respect to the rotational degrees of freedom R. To obtain the optimal translation, a fast translational search is performed for each candidate rotation over the full correlation function An external file that holds a picture, illustration, etc.
Object name is btt225i127.jpg by using FFT. Finally, the best combination of rotation and translation is chosen.

3 RESULTS

We assess the performance of our method by using simulated subtomograms of a phantom model. In addition, we test our method by experimentally determined structures of five complexes as well as experimental ribosome subtomograms extracted from a whole-cell tomogram.

3.1 Generation of simulated tomograms from phantom and experimental structures

We follow a previously applied methodology for generating subtomograms from initial structures by simulating the tomographic imaging process as realistically as possible, allowing for the inclusion of noise, tomographic distortions due to missing wedge and electron optical factors such as CTF and Modulation Transfer Function (MTF) (Beck et al., 2009; Förster et al., 2008; Nickell et al., 2005; Xu and Alber, 2012; Xu et al., 2011, 2012).

A density map of the complexes (Fig. 4) are generated by applying a low pass filter at 4 nm resolution to the atomic structures using the PDB2VOL program of the Situs 2.0 package (Wriggers et al., 1999) with voxel length of 1 nm. In addition, a density map of a phantom model is generated (Fig. 4A). These density maps are used as base maps, both of size 643 voxels, and they are randomly rotated and translated to generate density map instances in uncrowded conditions. To generate a crowded environment, two additional structures (randomly chosen from the phantom model or the ribosome complex) are randomly oriented and placed into the density map ~32 voxels away from the center of the map.

Fig. 4.
Structures used for simulating the tomographic image process. (A) Top: a phantom model that consists of four elliptical Gaussian functions as branches and one spherical Gaussian function to connect the elliptical functions. Bottom: Ribosome complex (PDB ...

The density maps of the target complex and neighboring structures are used as input for simulating electron micrograph images at different tilt angles. Our simulated subtomograms therefore contain a wedge-shaped region in reciprocal space for which no structure factors have been measured (i.e. the missing wedge effect), resulting in distortions in the density map as observed in the experimental measurements. To generate realistic micrographs, noise is added to the images according to a given SNR level (ranging between 0.05 and 0.001), defined as the ratio between the variances of the signal and noise (Förster et al., 2008) (Fig. 4A).

The resulting image is convoluted with a CTF, which describes the imaging in the transmission electron microscope in a linear approximation (Frank, 2006; Nickell et al., 2005). We also convolute the density map with the corresponding MTF of a typical detector used in tomography. Typical experimental acquisition parameters (Beck et al., 2009) were used: voxel grid length An external file that holds a picture, illustration, etc.
Object name is btt225i128.jpg nm, the spherical aberration An external file that holds a picture, illustration, etc.
Object name is btt225i129.jpgm, the defocus value An external file that holds a picture, illustration, etc.
Object name is btt225i130.jpgm, the MTF corresponded to a realistic electron detector (McMullan et al., 2009), defined as An external file that holds a picture, illustration, etc.
Object name is btt225i131.jpg where ω is the fraction of the Nyquist frequency. Finally, we use a backprojection algorithm (Nickell et al., 2005) to generate a tomogram from the individual 2D micrographs generated at the various tilt angles (Beck et al., 2009; Xu et al., 2011, 2012; Xu and Alber, 2012).

3.2 Segmentation of the target complex in simulated subtomograms

To determine the optimal scale for surface enhancement, we vary σ from 1.0 nm to 3.0 nm. At high noise levels (e.g. An external file that holds a picture, illustration, etc.
Object name is btt225i132.jpg), the optimal An external file that holds a picture, illustration, etc.
Object name is btt225i133.jpg is usually obtained at ~2.0 nm.

For the detection of boundary segments, we set the significance levels An external file that holds a picture, illustration, etc.
Object name is btt225i134.jpg. For determining the target complex, we choose An external file that holds a picture, illustration, etc.
Object name is btt225i135.jpg, ignoring all boundary segment maxima with negative intensities in An external file that holds a picture, illustration, etc.
Object name is btt225i136.jpg. For combining segments, we scan σ from An external file that holds a picture, illustration, etc.
Object name is btt225i137.jpg to 10 nm and scan the cluster number K from 1 to 5.

To assess the target complex segmentation, we calculate the overlap between the determined target complex A and the ground truth An external file that holds a picture, illustration, etc.
Object name is btt225i138.jpg. We also calculate the overlap between A and regions occupied by the surrounding structures in the crowded ground truth density map An external file that holds a picture, illustration, etc.
Object name is btt225i139.jpg. Then, we calculate the median of the following three overlap scores across 100 simulated subtomograms at each distortion level and for each of the benchmark structures:

equation image

where Ac denotes the complementary set of A. Our results show that for all distortion levels, the median of An external file that holds a picture, illustration, etc.
Object name is btt225i140.jpg is at least 0.9, and An external file that holds a picture, illustration, etc.
Object name is btt225i141.jpg is at most 0.004, demonstrating a good performance in segmenting the areas of the corresponding target structures. In contrast, the median of An external file that holds a picture, illustration, etc.
Object name is btt225i142.jpg is close to zero for most of the distortion levels, indicating that our target complex segmentation method can successfully exclude surrounding structures in the subtomograms. Only at high distortion levels, the median of An external file that holds a picture, illustration, etc.
Object name is btt225i143.jpg starts to deviate from zero (Supplementary Tables S7–S12).

3.3 Assessment of constrained center of mass detection

The constrained center of mass is expected to be invariant to rigid transforms of the corresponding structures. The assessment of the rotational invariance of our constrained center of mass is performed as follows. We randomly select 100 pairs of subtomograms for each of the benchmark structures. For each pair of subtomograms i and j, we calculate

equation image

as a measure of the degree of invariance, where An external file that holds a picture, illustration, etc.
Object name is btt225i144.jpg and An external file that holds a picture, illustration, etc.
Object name is btt225i145.jpg are the corresponding locations of the constrained center of masses in the original complex that is used to generate the instances by random rigid transforms. The smaller An external file that holds a picture, illustration, etc.
Object name is btt225i146.jpg, the higher is the degree of invariance, indicating a good performance.

The performance is expressed as the median of An external file that holds a picture, illustration, etc.
Object name is btt225i147.jpg, which is listed in Figure 5 and Supplementary Tables S1, S4, S13 and S16. We can see that the estimation error is generally smaller than 3, indicating high accuracy even at high distortion levels. As can be seen in the following section, the degree of invariance is sufficiently small to allow successful pairwise fast rotational alignments in real space.

Fig. 5.
The performance of the center of mass invariance, expressed as the median of An external file that holds a picture, illustration, etc.
Object name is btt225i148.jpg across 100 pairwise comparisons of subtomograms simulated at different distortion levels. The surface region color is rescaled hue saturation value (HSV) color and proportional ...

3.4 Fast rotational alignment of subtomograms

Our fast alignment method relies on fast rotational matching. In the following two subsections, we calculate the rotational alignment errors for the pairs of complexes, both for the crowded as well as uncrowded cases. The rotational alignment error is calculated as follows: Suppose An external file that holds a picture, illustration, etc.
Object name is btt225i149.jpg is the rotation matrix calculated for an alignment, and An external file that holds a picture, illustration, etc.
Object name is btt225i150.jpg is the true rotation matrix. The rotation alignment error can then be measured as

equation image

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

3.4.1 Uncrowded subtomograms

First, we test our alignment method with uncrowded subtomograms, which contain one complex without surrounding structures. We compare the alignment errors between our new alignment method presented here and our previous method that relied on fast rotational matching in reciprocal space (Xu et al., 2012) (left column of Fig. 6 and Supplementary Tables S2 and S5). It is evident that our new method improves the alignment accuracy on highly distorted subtomograms.

Fig. 6.
Rotational alignment error for the phantom model and the ribosome complex inside uncrowded and crowded subtomograms simulated at different distortion levels. The error is measured in terms of the median An external file that holds a picture, illustration, etc.
Object name is btt225i154.jpg across 100 pairwise alignments. The nodes on the ...

3.4.2 Crowded subtomograms

Next, we compare the two alignment methods (Xu et al., 2012) with crowded subtomograms, which also contain fragments of additional complexes (right column of Fig. 6 and Supplementary Tables S14 and S17). It can be clearly seen that when considering crowded subtomograms without excluding surrounding structures, our previous method (Xu et al., 2012) generally fails in finding the correct alignment. By contrast, our new alignment method aided by target complex segmentation is only marginally affected by crowding and the existence of the surrounding structures. Our new method clearly outperforms our previous method when dealing with subtomograms from crowded cellular environments. Similar improvement in performance is also observed by testing the crowded subtomograms of four additional benchmark complexes (Supplementary Tables S19–S22) (Fig. 4B) used in (Xu et al., 2012).

We have also measured the correlation between An external file that holds a picture, illustration, etc.
Object name is btt225i156.jpg and An external file that holds a picture, illustration, etc.
Object name is btt225i157.jpg (Supplementary Tables S3, S6, S15 and S18) and tested the performance with respect to the levels of crowdedness (Supplementary Figs S1–S3). At low to median distortion levels, we observe moderate correlation between the aforementioned two quantities. At high distortion levels, the correlations are small, showing that the alignment is more affected by high levels of distortions.

3.5 Segmenting target complexes in subtomograms extracted from experimental whole-cell tomogram of the human pathogen Leptospira interrogans

We tested our method also on experimental subtomograms extracted from a whole-cell tomogram of Leptospira interrogans in undisturbed conditions (Beck et al., 2009). Segmenting of complexes in such subtomograms is significantly more challenging because the macromolecular crowding is high and because the SNR levels and resolution of these subtomograms is typically low. Beck et al. (2009) have identified several complexes in the tomogram using template matching (Best et al., 2007). The templates were generated by using density maps of structures in the Protein Databank and convoluted with CTF and MTF. As a test case, we use two subtomograms with high template matching scores to the ribosome template. For each of these matches, we extracted a subtomogram of size 493 voxels that was large enough to contain an instance of the complex and certain amount of surrounding structure fragments.

The results shown in Figure 7 demonstrate that our method can successfully segment the target ribosome complex while surrounding structures are clearly excluded. The ribosome positions agree with those detected in the template matching. We further extended the study to the subtomograms corresponding to the top 20 template matching scores. Overall, 18 of them showed successful segmentation (Supplementary Fig. S4). We also analyzed the target complex segmentation performance with respect to the level of crowdedness. The crowdedness of the Leptospira interrogans subtomograms ranges from 0 to 60%. The true-positive and false-positive rates plotted in Supplementary Figure S4 shows that, in most situations, our method can correctly identify the target complex region. Only at high crowdedness level, (i.e. >50%), our method fails to correctly identify the target complex regions.

3.6 Analysis of computational costs

We implemented the methods in MATLAB and carried out our tests on a computer cluster consisting of 2.3–2.7 GHz computer nodes. On average, the target complex detection step takes 14.8 s. In this step, the major proportion of computation time was used for the scale selection (8.6 s on average). In practice, we can determine a fixed scale parameter σ from a set of training subtomograms and directly apply it to all subtomograms. Then, on average, only 6.2 s would be used for detecting the target complex. The computer memory consumption is small. Because a subtomogram is typically small (a subtomogram of size 643 corresponds to at most 2 M memory), our detection step uses no more memory than 20 times the size of a single subtomogram. The most computational intensive step in the classification process is the alignment. The computational cost of our new method is essentially the same as our previous method (Xu et al., 2012) and therefore allows high-throughput subtomogram classifications.

4 CONCLUSION

Cryo-electron tomograms provide useful information for simultaneously detecting the native structures and their spatial cellular locations of a large number of macromolecular complexes. However, the high level of macromolecular crowding and distortions in the extracted subtomograms makes such analysis challenging. In addition, the availability of a large number of subtomograms containing potential complexes requires high-throughput analysis techniques. In this article, we propose an automatic segmentation technique that isolates the target complex inside a crowded subtomogram. We further detect the rigid transform invariant center of mass of the target complex and propose an improved fast alignment method using fast rotational matching in real space. The new method shows good performance in segmenting target complexes in crowded subtomograms and performing fast rotational alignments using both simulated and experimental subtomograms. This performance is a necessary condition for high-throughput structural classifications of complexes in whole-cell tomograms. This work represents a step toward high-throughput and reference-free Visual Proteomics analysis of highly crowded whole-cell cryo-electron tomograms.

Supplementary Material

Supplementary Data:

ACKNOWLEDGEMENTS

The authors thank Dr Martin Beck for providing the experimental tomograms for testing. They thank the anonymous reviewers for their constructive suggestions.

Funding: Human Frontier Science Program RGY0079/2009-C, the Arnold and Mabel Beckman foundation; NIH [R01GM096089 and U54RR022220] and NSF CAREER [1150287] (to F.A.). F.A. is a Pew Scholar in Biomedical Sciences, supported by the Pew Charitable Trusts.

Conflict of Interest: none declared.

REFERENCES

  • Amat F, et al. Subtomogram alignment by adaptive Fourier coefficient thresholding. J. Struct. Biol. 2010;171:332–344. [PubMed]
  • Bartesaghi A, et al. Classification and 3D averaging with missing wedge correction in biological electron tomography. J. Struct. Biol. 2008;162:436–450. [PMC free article] [PubMed]
  • Beck M, et al. Visual proteomics of the human pathogen Leptospira interrogans. Nat. Methods. 2009;6:817–823. [PMC free article] [PubMed]
  • Beck M, et al. Exploring the spatial and temporal organization of a Cell’s proteome. J. Struct. Biol. 2011;173:483–496. [PubMed]
  • Best C, et al. Localization of protein complexes by pattern recognition. Methods Cell Biol. 2007;79:615–638. [PubMed]
  • Beucher S, Lantuéjoul C. International Workshop on Image Processing, Real-Time Edge and Motion Detection/Estimation, CCETT/IRISA, Rennes, France. 1979. Use of watersheds in contour detection.
  • Böhm J, et al. Toward detecting and identifying macromolecules in a cellular context: template matching applied to electron tomograms. Proc. Natl Acad. Sci. USA. 2000;97:14245–14250. [PubMed]
  • Celeux G, Govaert G. Gaussian parsimonious clustering models. Pattern Recognit. 1995;28:781–793.
  • Chen Y, et al. Fast and accurate reference-free alignment of subtomograms. J. Struct. Biol. 2013 [Epub ahead of print, doi: 10.1016/j.jsb.2013.03.002, March 22, 2013] [PubMed]
  • Förster F, et al. Classification of cryo-electron sub-tomograms using constrained correlation. J. Struct. Biol. 2008;161:276–286. [PubMed]
  • Fraley C, Raftery A. Model-based clustering, discriminant analysis, and density estimation. J. Am. Stat. Assoc. 2002;97:611–631.
  • Frangakis A, et al. Identification of macromolecular complexes in cryoelectron tomograms of phantom cells. Proc. Natl Acad. Sci. USA. 2002;99:14153–14158. [PubMed]
  • Frank J. Three-dimensional electron microscopy of macromolecular assemblies. New York: Oxford University Press; 2006.
  • Garzón J, et al. Adp_em: fast exhaustive multi-resolution docking for high-throughput coverage. Bioinformatics. 2007;23:427–433. [PubMed]
  • Heumann J, et al. Clustering and variance maps for cryo-electron tomography using wedge-masked differences. J. Struct. Biol. 2011;175:288–299. [PMC free article] [PubMed]
  • Hrabe T, et al. Pytom: a python-based toolbox for localization of macromolecules in cryo-electron tomograms and subtomogram analysis. J. Struct. Biol. 2012;178:177–188. [PubMed]
  • Kovacs J, Wriggers W. Fast rotational matching. Acta. Crystallogr. D. Biol. Crystallogr. 2002;58:1282–1286. [PubMed]
  • Lindeberg T. Scale-Space Theory in Computer Vision. The Kluwer International Series in Engineering and Computer Science. Dordrecht, The Netherlands: Kluwer Academic Publishers; 1994.
  • Lucic V, et al. Structural studies by electron tomography: from cells to molecules. Annu. Rev. Biochem. 2005;74:833–865. [PubMed]
  • Martinez-Sanchez A, et al. A differential structure approach to membrane segmentation in electron tomography. J. Struct. Biol. 2011;175:372–383. [PubMed]
  • McMullan G, et al. Detective quantum efficiency of electron area detectors in electron microscopy. Ultramicroscopy. 2009;109:1126–1143. [PMC free article] [PubMed]
  • Medalia O, et al. Macromolecular architecture in eukaryotic cells visualized by cryoelectron tomography. Science. 2002;298:1209–1213. [PubMed]
  • Nickell S, et al. TOM software toolbox: acquisition and analysis for electron tomography. J. Struct. Biol. 2005;149:227–234. [PubMed]
  • Nickell S, et al. A visual approach to proteomics. Nat. Rev. Mol. Cell. Bio. 2006;7:225–230. [PubMed]
  • Volkmann N. A novel three-dimensional variant of the watershed transform for segmentation of electron density maps. J. Struct. Biol. 2002;138:123–129. [PubMed]
  • Volkmann N. Methods for segmentation and interpretation of electron tomographic reconstructions. Methods Enzymol. 2010;483:31–46. [PubMed]
  • Witkin AP. International Joint Conference on Artificial Intelligence, Karlsruhe, Germany. 1983. Scale-space filtering; pp. 1019–1022.
  • Wriggers W, et al. Situs: a package for docking crystal structures into low-resolution maps from electron microscopy. J. Struct. Biol. 1999;125:185–195. [PubMed]
  • Xu M, Alber F. High precision alignment of cryo-electron subtomograms through gradient-based parallel optimization. BMC. Syst. Biol. 2012;6(Suppl. 1):S18. [PMC free article] [PubMed]
  • Xu M, et al. Template-free detection of macromolecular complexes in cryo electron tomograms. Bioinformatics. 2011;27:i69–i76. [PMC free article] [PubMed]
  • Xu M, et al. High-throughput subtomogram alignment and classification by Fourier space constrained fast volumetric matching. J. Struct. Biol. 2012;178:152–164. [PubMed]

Articles from Bioinformatics are provided here courtesy of Oxford University Press