A clustering approach to multireference alignment of single-particle projections in electron microscopy

2.1. Similarity between two images: correntropy

The cross correntropy between two random variables (X, Y) is defined as

where E{·} is the expectation operator, and κ_σ(x) is a one-dimensional symmetric (κ_σ(x) = κ_σ(−x)), non-negative kernel. In practice, the true distribution of X − Y is unknown and the cross correntropy can be approximated by its empirical estimate

where x_i and y_i are actual measurements of the X and Y variables.

Before comparing images X and Y, both images should be aligned. For doing this, we first translationally align image X with respect to Y using the correlation property of the Fourier transform in cartesian coordinates producing a new image X′. Then, we rotationally align X′ with respect to Y using the correlation property of the Fourier transform in polar coordinates producing a new image X″. We repeat this process twice before computing the correntropy between both images. Since the Shift + Rotation sequence is not guaranteed to reach the correct alignment between the two images (the global minimum), we also perform two iterations of the Rotation + Shift sequence, and pick the aligned image that maximizes the correntropy between the two images.

Interesting properties of the correntropy are that it is symmetric (V_σ(X, Y) = V_σ(Y, X)), positive, bounded, it is the argument of Renyi’s quadratic entropy (from which the name correntropy comes from), it is related to robust M-estimates (Liu et al., 2007), and it involves all the even moments of the variable X − Y (not just its second order moment, as is the case of correlation):

where the coefficients a_2n depend on the Taylor expansion of the kernel (note that the Taylor expansion has only even terms because the kernel is even symmetric and by definition cannot have odd terms in its Taylor expansion). It is this latter property which makes it specially suited to non-Gaussian signal processing and our application. Let us assume that X and Y are two identical images except for the noise (assumed to be white and Gaussian). It can be easily seen that the two images are correctly aligned when the variance of the difference image, E{(X − Y)²}, is minimum (if it is not minimum it means that there is still some misalignment between X and Y yielding a difference in the signals themselves, not only the noise). At the point of minimum variance, the correntropy is also minimum. However, thanks to the higher order terms in the Taylor expansion of the correntropy, these differences in X − Y due to the misalignment contribute not only to the 2nd power (as in the variance), but also to the 4th, 6th, …(each term with less and less weight, a_2n, so that the whole series is convergent; note also that misalignments make the distribution of the difference image to be non Gaussian). The result is that if we plot the correntropy (or variance) landscape around its optimal alignment, the landscape of the correntropy has a much sharper optimum than that of the variance, meaning that it is more sensitive to small misalignments. The same reasoning applies when there are small differences between X and Y, not only because of the noise and misalignments, but because X and Y are slightly different images (they belong to different classes).

In our algorithm we use the unnormalized Gaussian kernel

which is bounded between 0 (no similarity between X and Y) and 1 (X and Y are identical).

A key point of kernel algorithms is the estimation of the kernel width σ. In our algorithm we first estimate the noise variance ${\sigma }_N^2$ from the background of the experimental images (defined as the region outside the maximum circle inscribed in the projection image). This estimation is performed over the whole dataset. Let us assume that X is an experimental image and Y its corresponding true class which has been calculated without noise. Then the variance of the variable X − Y is ${\sigma }_N^2$. However, in practice, Y is never calculated without noise. Let us assume that Y has been calculated from N_Y experimental images. Then, the noise in Y has a variance of $\frac{{\sigma }_N^2}{N_Y}$, and the variance of image difference X − Y is ${\sigma }_N^2+\frac{{\sigma }_N^2}{N_Y}$. Therefore, we set for our kernel

2.2. Divisive vector quantization

For completeness we present here the standard multireference classification in the framework of vector quantization (or clustering). Then, we present the standard divisive vector quantization algorithm, which is known to produce much more robust results in vector quantization, and show how this divisive algorithm can be applied to multireference classification in electron microscopy.

Considering each experimental image, X_i, as a vector in a P dimensional space (being P the total number of pixels), multireference classification with K class centers aims at finding the vectors X̂_k (k = 0, 1,…,K − 1) such that Σ_i∥X_i − X̂_k(i)∥² is minimized, being k(i) the index of the class center assigned to the ith experimental image. This is exactly the same goal as the one of vector quantization with K code vectors.

A divisive clustering algorithm starts with N₀ class centers (in our algorithm we have chosen averages of random subsets from the experimental set of images). Then, the algorithm assigns each image to the closest class center, and recomputes the class centers as the average of the images assigned (see Fig. 1). This process is iterated until the number of images that change their assignment from iteration to iteration is smaller than a certain threshold (we used 0.5% of the total number of images in the experimental dataset) or a maximum number of iterations is reached. Once the N₀ class averages have been computed, the divisive clustering algorithm chooses the class center with the largest number of images assigned to it, and splits this subset of images into two new classes. The split method computes the similarity (either through correntropy or correlation) of each image in the cluster with the class center of that cluster. Then, it splits the whole set of images into two halves: the one with the 50% highest similarities, and the one with the 50% worse similarities. With this initial split, we compute two new class centers and apply the same clustering methodology as for the whole dataset (of the two clustering criteria discussed in the next section, for the split in two classes we found more useful the newly proposed robust criterion). After splitting the largest class, the algorithm proceeds to split the second largest class, and this process is iterated until N₀ splits have been performed. Then, the clustering algorithm is run again with 2N₀ class centers. Split phases and clustering phases are alternated until the desired number of class centers is reached.

Alternatively, we could have started directly with the final number of class centers and let the clustering optimize the initial class centers (this is exactly how multireference alignment is usually applied in electron microscopy). However, this algorithm is known to get more easily trapped into local minima. Divisive clustering is less prone to local minima (Gray, 1984), although it is not guaranteed that the global minimum will be reached (this is a general drawback of K-means algorithms).

Our algorithm is an adaptation of this standard clustering algorithm to the cryo-EM setting. Since our main goal is recognizing the most common projections in the dataset, if at any moment a class represents less than a user given percentage of the total number of images (in the order of $0.2\frac{N_{\mathit{img}}}{K}$, being N_img the total number of images and K the current number of classes), we remove that class from the quantization and split the class with the largest number of images assigned. Note that we do not remove the corresponding experimental images from the dataset. Instead, in the next quantization iteration, these images will be assigned to one of the remaining clusters.

2.3. Clustering criterion

Instead of using the L₂ norm of the error to measure the similarity between an experimental image and a class representative, we propose to use the correntropy introduced earlier. Correntropies closer to 1 mean a larger agreement between the experimental image and the class representative.

Unfortunately, simply assigning an image to the class with highest correntropy does not result in a good classification of the images. The problem is that the absolute value of the correntropy is meaningless when comparing an experimental image to several classes (a similar problem has already been reported in electron microscopy using the correlation index (Sorzano et al., 2004a)). The reason for this is that at low Signal-to-Noise Ratios, the differences between two images may be absolutely masked by the noise, and the wrong assignment may be done simply because one class average is less noisy than another one. It could be said in this situation that the cleanest class “attracts” many experimental images even if they belong to some other class.

Let us prove this statement with a simple example. Let us presume that our experimental data comes from only two underlying classes A₀ and A₁, and that, after alignment, the experimental images are simply noisy versions of these two classes: X_i = A_j + N_i, where X_i is the ith experimental image, A_j is either A₀ or A₁, and N_i is a white noise image added during the experimental acquisition. Let us assume that at a given moment we have perfectly aligned and classified the experimental images in these two classes and that we have M₀ images in the class generated by A₀ and M₁ images in the class generated by A₁. Note that A₀ and A₁ are unknown (in fact, the multireference classification problem consists in estimating them), and all we have are estimates obtained by averaging the images assigned to each class, i.e., Â_j = A_j + N̂_j, where Â_j are our estimates of the true underlying class averages and N̂_j is some noise image obtained as the average of all the noise images of the experimental images assigned to this class. If the variance of the experimental images is σ², then the variance of N̂_j is $\frac{{\sigma }^2}{M_j}$.

The classical multireference clustering criterion is “Choose for X_i the class A_j that minimizes

where P is the total number of pixels, while X_ip and Â_jp denote the pth pixel of the experimental image and the class average, respectively. Let us assume that X_i really belongs to class A₁, and let us see under which conditions it would be correctly assigned to class 1:

In other words, the experimental image will fail to be correctly classified unless the SNR of the difference between the two underlying classes is larger than a certain threshold (in electron microscopy images this is not always the case, as we try to capture very subtle differences in very noisy images). Let us further assume that there are many more images in A₀ than in A₁, i.e., $\frac{1}{M_0}$ is “negligible” versus $\frac{1}{M_1}$, moreover M₁ tends to be a small number and $\frac{1}{M_1}$ is particularly high. In other words, if one of the classes is large (M₀ ≫ M₁), X_i will fail to be assigned to its correct class and the class with more images assigned, A₀, “attracts” many images from the other class.

Although we have presented the problem with the variance, this problem also occurs when using correntropy or maximum likelihood (Scheres et al., 2005). We propose an alternative clustering criterion that we will refer to as the robust criterion. Instead of looking only at the energy of the difference between the experimental image and the two estimates of the class averages, we propose to compare this energy to the energy of the rest of images compared to each class average. For instance, consider an experimental image whose correntropy with the representative of class A is 0.8, and with the representative of class B is 0.79. Because of the noise, these numbers cannot be blindly trusted and simply assign the experimental image to class A. Let us consider the correntropies of all images assigned to class A and of all images assigned to class B. If most correntropies of images in class A is in the order of 0.9, our image would be a very bad member of class A. However, if most correntropies of images in class B is in the order of 0.7, our image would be a very good member of class B (despite the fact that in absolute terms, the correntropy to class B is a little bit smaller than that to class A).

How good a correntropy is with respect to each one of these sets can be measured through the distribution functions

The first one measures the probability of having a correntropy smaller than a given value υ within the set of images assigned to class k. The second one measures the probability of having a correntropy smaller than a given value υ within the set of images not assigned to class k. The class to which an image X_i is assigned should be the one maximally fulfilling both goals

The classical clustering criterion is good for large differences between the class averages, while the robust clustering criterion is specially designed for subtle differences embedded in very noisy images. That is why we prefer to use this second criterion when splitting the nodes during the divisive clustering algorithm.

2.4. Final algorithm: CL2D

With the pieces introduced so far (correntropy, divisive clustering, and robust clustering criterion) we propose to use the following clustering algorithm that we will refer to as CL2D.

Algorithm 1. CL2D

Input:

S: Set of images

K0: Number of classes in the first iteration

KF: Number of classes in the last iteration

Output:

CKF : Set of KF class representatives

begin

• // Initialize the classes

• k = K0

• Ck = Randomly split S into k classes

• Update the distributions of Eq. (7)

• σ2 = Measure the noise power in the background of the experimental images

• // Refine the classes

• Ck = Refine Ck with the data in S

• while k < KF do

  • Cmin(2k, KF) = Split the largest min(k, KF – k) classes of Ck

  • k = min(2k, KF)

  • Ck = Refine Ck with the data in S

• end

end

The refinement of the current classes is performed according to the following algorithm.

Algorithm 2. CL2D refinement of the current class representatives (C_k)

Input:

S: Set of images

Ck: Set of k class representatives

σ2: Noise power in the background of experimental images

Itermax: Maximum number of iterations

Nmin: Minimum number of images in a class

Output:

Ck: Refined set of k class representatives

begin

• Iter = 0

• repeat

  • // Assign all images to a class

  • foreach image ∈ S do

    • foreach classAverage ∈ Ck do

      • image =Align image to classAverage

      • Compute argument of Eq. (8)

    • end

    • Assign image to the class maximizing Eq. (8)

  • end

  • // Remove small classes and split the largest ones

  • while there are classes with less than Nmin images do

    • Remove the class with the smallest number of images assigned

    • Split the largest class of Ck

  • end

  • Update the distributions of Eq. (7)

  • Iter = Iter + 1

• until No: of Changes < 0:5% and Iter < Itermax;

end

The split of a class is performed in a way very similar to that of the refinement of the current classes. The main difference is how the two subclasses are initially calculated. While in the standard algorithm these are computed by a random split, in the CL2D split we separate the images with the largest correntropies to the class average from the images with the smallest correntropies:

Algorithm 3. CL2D split of a class

Input:

S(c): Set of images assigned to class c

c: Class representative of class c

Output:

c1, c2: Two subclasses derived initially from class c

begin

• // Split the images into two initial subclasses

• c1 = Average of all images assigned to c with the largest 50%

• correntropies

• c2 = Average of all images assigned to c with the smallest 50%

• correntropies

• Update the distributions of Eq. (7) for these two subclasses

• // Refine the two subclasses

• Refine {c1, c2} with the data in S(c)

end

3.1. Simulated data: bacteriorhodopsin

To test our algorithm we generated 10,000 projections randomly distributed in all projection directions from the bacteriorhodopsin monomer (PDB entry: 1BRD) and added white Gaussian noise with a SNR of 1/3 and 1/30 (see Fig. 2). We then applied our algorithm to compute 256 class averages with a maximum iteration count of 20. We found that this parameter is not critical as long as it is not too small so that the algorithm is stopped when the classes are not well established yet.

In our first experiment we tried to assess the effectiveness of each one of the three main modifications to the clustering algorithm (correntropy vs. correlation, divisive clustering vs. multireference clustering, classical clustering criterion vs. robust clustering criterion). For doing so, we conducted four experiments with the dataset with SNR = 1/3: the first one with the three new features (see its results in Supp. Fig. 1), and other three in which one of the new features was missing (correntropy, divisive clustering, or robust clustering). For each one of the experiments we evaluated each cluster by computing the angle between the projection directions of any pair of images assigned to that cluster. Finally, we computed the histogram of these angular distances in order to compare the different choices of the algorithm. Fig. 3 shows these histograms. It can be seen that the best combination is the one using two of the three new ideas introduced in this paper: correntropy and divisive clustering. As has already been mentioned before, this simply means that the SNR of the differences between classes is large enough to overcome the attraction effect (in the next subsection we present a case in which the robust clustering criterion performs better than the classical one because the input images have much lower SNR).

We repeated this experiment using ML2D (see Supp. Fig. 2). The algorithm could not produce more than 20 different classes that accounted for 97.1% of the particles (the algorithm was run with 50 classes, and the remaining 30 classes accounted only for 58 images, existing even empty classes). This result illustrates the “attractive” effect described during the introduction of the robust clustering criterion. Moreover, there was at least one class average (with 129 images assigned) in the ML2D approach that did not really represent any of the true projection images (see Fig. 4). In fact, if we analyze the 129 images assigned to this node with our algorithm, it can be seen that the ML2D class is actually a mixture of at least four different classes with 29, 35, 35, and 30 images, respectively.

We also repeated this experiment with SVD/MSA (classes not shown but we did not observe the “attraction” effect), PCA/Diday, PCA/Hierarchical, and PCA/K-means (with both dataset SNR = 1/3 and SNR = 1/30). As suggested in the multivariate data analysis web page of Spider (http://www.wadsworth.org/spider_doc/spider/docs/techs/MSA), we filtered our images with a Butterworth filter between the digital frequencies 0.25 and 0.33 (these frequencies are normalized in such a way that 1 corresponds to the Nyquist frequency). We also prealigned the images and reduced their dimensionality using a PCA with nine components. This preprocessing was applied to the data before using the three algorithms of Spider. In Fig. 5 we show the different estimates of the probability density function of the angular distance between the projection directions of images assigned to the same cluster for the six algorithms. It can be clearly seen that CL2D produces the tightest clusters in terms of projection directions.

The execution times in a single CPU (Intel Xeon 2.6 GHz) for PCA/Diday, PCA/Hierarchical and PCA/K-means was about 0.3 h, 50 h for SVD/MSA, 100 h for CL2D, and 74 h for ML2D (with only 64 classes; CL2D took 41 h in computing 64 classes). However, CL2D has been parallelized with MPI and has been run in parallel with 32 nodes, reducing the computation time to 3.5 h.

The experiment was repeated in a more realistic setup by introducing the effect of the microscope Contrast Transfer Function (CTF) with an acceleration voltage of 200 kV, a defocus of −5.5 μm, and a spherical aberration of 2.26 mm. Results were similar to the previous case (results not shown).

3.2. Simulated data: E. coli ribosome

For this test we used the public dataset of simulated ribosomes available at the Electron Microscopy Data Bank (http://www.ebi.ac.uk/pdbe/emdb/singleParticledir/SPIDER_FRANK_data, (Baxter et al., 2009)). This dataset contains 5000 projections from random directions of a ribosome bound with three tRNAs at A, P, and E sites, and other 5000 projections from random directions of an EF-G(GDPNP)-bound ribosome with a deacylated tRNA bound in the hybrid P/E position. Besides these differences in the ligands, these two ribosomes also have different ribosomal conformations: the first ribosome is in the normal conformation, while the second ribosome is in a ratcheted configuration (30S subunit rotated counter-clockwise relative to the 50S subunit). Fig. 6 shows some sample images from this dataset.

The goal of this experiment is to characterize the capability of the new algorithm to separate the two conformational states, i.e., to produce classes in which the majority of the images assigned comes from the same conformational state.

With our previous dataset we already showed the superiority of the correntropy over correlation, and of the divisive clustering over a classical multireference clustering. In the previous dataset it also appeared that the classical clustering criterion was superior to the robust clustering criterion. However, as was already discussed in the Methods Section, this holds for sufficiently high SNRs. This dataset has a much lower SNR than the previous one. We ran our algorithm to partition the input data into 256 classes using the classical clustering criterion and the new robust criterion. For each class we tested the hypothesis that the two ribosomes were equally represented (the proportion of images of each type was 0.5). The classical clustering criterion only found 67 classes (out of 256) where the proportion of the majority of the images were significantly different (with a 95% of confidence) from 50%. The amount of images assigned to these classes was 2822. However, the new robust clustering criterion was able to identify 98 classes where the majority of images significantly came from a single ribosome type. The amount of images assigned to these 98 classes was 4130. The difference between the proportions of classes with a majority of ribosomes of one type using the classical and the robust clustering criteria was significantly different with a confidence of 95%.

We repeated the same experiment with SVD/MSA, ML2D, PCA/Diday’s clustering, PCA/Hierarchical clustering and PCA/K-means. SVD/MSA identified 107 classes with a statistically significant majority of images of one type (4534 images were assigned to these classes). However, the difference between the proportions of SVD/MSA and CL2D was not statistically significant with a confidence level of 95%. ML2D identified 9 classes with a statistically significant majority of images of one type (1439 images were assigned to these classes). Images were preprocessed in exactly the same way as in the previous experiment before applying any of Spider’s methods. PCA/Diday method produced 11 classes with a statistically significant majority of images of one type (899 images were assigned to these classes). The PCA/Hierarchical clustering produced 74 classes with a statistically significant majority of images of one type (with 3049 images assigned to them). The PCA/K-means produced 53 classes with a statistically significant majority of images of one type (with 2114 images assigned to them).

3.3. Experimental data: p53

We applied our algorithm to a dataset that contained 7600 projections of a mutant of p53 without the C-terminal domain bound to the GADD45 DNA sequence (Ma et al., 2007; Tidow et al., 2007). The sample was negatively stained. The pixel size was 4.2 Å/pixel and the image size was 60 × 60 pixels (see Fig. 7). We applied our algorithm to generate 256 class averages (20 iterations per Hierarchical level). Supp. Fig. 3 shows the resulting classes. In an experimental setting it is difficult to validate the 2D classes produced by an algorithm. In this experiment at least we tested that the classes produced captured enough information from the original dataset. For doing so, we built a reference volume based on the common lines of the classes produced by CL2D (see Fig. 8a). This reference volume was refined in two different ways: first, using the CL2D classes themselves as the only projections in the dataset (the refined volume is shown in Fig. 8b); second, using the whole dataset of projections (the refined volume is shown in Fig. 8c). Interestingly, the volumes in Fig. 8b and c are consistent with each other up to 32 Å (at this frequency the Fourier Shell Correlation (Harauz and van Heel, 1986) between the two volumes drops below 0.5). Considering that the resolution of the volume in Fig. 8c is around 27 Å, the approximation up to 32 Å of the volume reconstructed using only the CL2D classes is not a bad approximation. The same projection dataset was refined using a reference volume obtained independently through Random Conical Tilt (Scheres et al., 2009).