We study the problem of computing a compressed representation of the covariance matrix of 2-D cryo-EM images for the purpose of performing principal component analysis (PCA). More precisely, we consider an image formation model where the measurement is defined by
where is a radial function, denotes convolution, is an image, and the noise term.
We are motivated by single particle cryo-electron microscopy (cryo-EM) imaging, which is an important technique for determining the 3-D structure of macromolecules. In particular, the single particle reconstruction (SPR) problem asks to recover the 3-D structure of a macromolecule from noisy 2-D images of its tomographic projections along unknown viewing angles. In cryo-EM, the mathematical model is a special case of (1) and is of the form
where are 3-D spatial coordinates, is the point spread function, is the electrostatic potential of a molecule, is a 3-D rotation, and
the noise term. In computational microscopy, it is typical to work with the Fourier transform of the point spread function, which is known as the contrast transfer function (CTF). In the simplest case, each measurement could correspond to a single fixed molecule potential function; however, in general, we may assume that each
could be a random variable representing a mixture of molecules, conformational heterogeneity, cases where the images are not perfectly centered, or other measurement imperfections[20, 22].
In general, each measurement can be associated with a different point spread function; however, in practice, a group of measurements, called a defocus group, can share a common point spread function. We assume that the measurements are grouped into defocus groups. Given and for , our goal is to estimate the 2-D covariance function of the images
where is a random variable from the same distribution as the images , and . We assume that the distribution of the images is invariant to in-plane rotations (which is typically the case in cryo-EM).
In cryo-EM, the random variable is of the form , where the random variable is an unknown viewing angle, and the random variable is a molecule potential. There is generally no physical reason for a molecule to prefer one in-plane rotation to another so distributions of random variables of this form are generally invariant to in-plane rotations. In the field of cryo-EM processing, the covariance function is simply referred to as the 2-D covariance.
The 2D-covariance is an essential component of a number of computational techniques in cryo-EM; we survey a few of these below.
First, we are motivated by PCA, which is a ubiquitous technique in statistics, data science, and computational mathematics and has applications to dimensionality reduction, denoising, visualization, among others. The principal components (that is, the top eigenvectors of the digitized covariance matrix) have a number of uses in the computational cryo-EM pipeline. The subspace corresponding to the top eigenvectors of the covariance matrix identifies salient features of the dataset which enables, for instance, improved methods for image classification and visualization, such as Multivariate Statistical Analysis[27, 26, 28, 29]
. These techniques improve computational speed, since clustering becomes computationally easier in a space of reduced dimension, as well as accuracy, since dimensionality reduction by PCA amplifies the effective signal-to-noise ratio because many coordinates for which noise dominates the signal are eliminated[32, Table 3].
Second, the covariance matrix has applications in the method of moments, a classical statistical inference method, applied to cryo-EM. In this method, the 2-D covariance is used to compute the similarly defined autocorrelation function of the underlying 3-D structure. Under further assumptions such as sufficient non-uniformity of the distribution of the viewing angles  or sufficient sparsity of the molecular structure , this autocorrelation function determines the 3-D density map either up to a finite list of possible structures or uniquely, respectively. This has been further developed into principled methods for ab initio estimation of cryo-EM structures [24, 2], with reduced risk of user-induced model bias in the initial model. Alternatively, when additional information is available, for instance one  or two  noiseless projection images, or the 3-D structure of a related, homologous structure [3, 5], the 3-D density map is uniquely determined by the autocorrelation, without requiring any structural assumptions.
Third, the covariance matrix has applications to denoising and CTF-correcting projection images. Covariance Wiener Filtering (CWF)  is an approach which uses the classic Wiener filtering framework with the estimated covariance matrix to solve the image deconvolution and denoising problem. The technique represents images in a lower dimensional subspace that is formed from PCA using the estimated covariance matrix. The method then applies Wiener filtering to correct the CTFs and denoise the images in this reduced subspace.
Compared to the standard PCA problem, the cryo-EM setting exhibits further computational challenges, since the estimation method also has to account for convolution with the point spread function, which destroys information of the resulting convolved function; see §2.2 for a more precise statement. On the other hand, the problem has additional symmetries making fast algorithms possible. In this paper, we present a new fast algorithm for estimating the covariance matrix that improves upon past approaches (especially when there are a large number of defocus groups) in terms of time and space complexity.
1.2. Main contribution
The main contribution of this paper is a new computational method for estimating the covariance (3) from measurements of the form (1) encoded by digitized images. The presented fast method has time complexity independent of the number of defocus groups. This is in contrast to past methods, where this complexity scales poorly with and involves operations , where is the number of iterations needed in a Conjugate Gradient step. Many modern cryo-EM experimental datasets fall into the computationally challenging regime where scales with .
Our fast method hinges on a new fast and accurate method for expanding images into the Fourier-Bessel basis, which provides a convenient way to handle convolution of radial functions (such as point spread functions) with images: namely, convolution with radial functions can be expanded as a diagonal operator operating on the basis coefficients .
The Fourier-Bessel basis functions are harmonics
on the disk: the standing waves associated with the resonant frequencies of a disk shaped drum with a fixed boundary. More precisely, the harmonics on the disk are eigenfunctions of the Laplacian on the unit disk that satisfy Dirichlet boundary conditions. In computational mathematics, this basis is referred to as the Fourier-Bessel basis, since the basis functions can be expressed as a product of a Bessel function and a complex exponential; see (5) for a definition.
Because of this simple structure, the covariance matrix of clean images can be estimated by a simple closed-form solution, without using the (computationally expensive) conjugate gradient method from previous approaches. Simultaneously, the covariance matrix retains its block diagonal structure, meaning that its diagonal blocks can be estimated separately and independently, which altogether makes PCA fast.
We present numerical results of covariance estimation on synthetic and experimental data. Additionally, we show how the estimated covariance matrix can be used to denoise images using CWF, and perform PCA to visualize eigenimages from experimental data. Code implementing the method is publicly available online111Code is available at https://github.com/yunpeng-shi/fast-cryoEM-PCA. Moreover, our approach has the potential to generalize to settings beyond cryo-EM, where PCA is used for signals estimated under more general group actions .
For two -matrices and , we denote their Hadamard (or entrywise) product by , the Hadamard division of and by and the th Hadamard power of by . These operations are defined elementwise by
respectively. If is an
-dimensional vector, thendenotes the matrix with along its diagonal, i.e., , and zeros elsewhere. If is a radial function, we write to mean that can be expressed as a function only of the magnitude of .
2.2. Technical details
We make the following assumptions
We assume that the point spread functions in the model (1) are radial functions; this implies that their Fourier transform (the CTFs) are also radial. In systems where astigmatism is present and the point spread function deviates slightly from a radial function, our approach can be used as an initial approximation that could be refined using the Conjugate Gradient method.
We assume that the underlying images in the model (1) are i.i.d. random variables whose distribution is invariant to in-plane rotations.
Informally speaking, assumption 3 can be interpreted as saying that for any pair of frequencies and there is a point spread function whose Fourier transform does not vanish at or . Assumption 1 implies that is a radial function. Figure 1 shows the values of for each pair of radial frequency in log scale, for 1081 distinct CTF images of size , whose defocus values range from 0.81 to 3.87 , for an experimental dataset; see §4 for more details.
This assumption is much weaker than assuming that for all that does not vanish at any frequency. In the latter case, we could just use to invert each equation to get access to the underlying functions . If we had direct access to the underlying functions , then we could approximate the covariance function by the sample covariance
where is the sample mean. Indeed,
by the law of large numbers. To clarify why it is useful fornot to vanish, note that in the Fourier domain our measurement model can be expressed as
where , and denote the Fourier transforms of and , respectively. Since taking the Fourier transform changes convolution to point-wise multiplication, if each Fourier transform for some , then we could estimate by first estimating from for , and then using the sample covariance matrix. However, in practice, the CTF is approximately a radial function with many zero-crossings, which means that multiplication by destroys information in the corresponding frequencies, making the restoration from a single image ill-posed.
The main ingredient in our fast covariance estimation is a fast transform into a convenient and computationally advantageous basis, known as the Fourier-Bessel basis (which consists of the harmonics on the disk: eigenfunctions of the Laplacian on the disk that obey Dirichlet boundary conditions). This specific choice of basis has a number of beneficial properties:
it is orthonormal,
it is ordered by frequency,
it is steerable i.e., images can be rotated by applying a diagonal transform to the basis coefficients,
it is easy to convolve with radial functions i.e., images can be convolved with radial functions by applying a diagonal transform to the basis coefficients.
These properties have made the Fourier-Bessel basis a natural choice in a number of imaging applications [25, 4, 23, 32, 31] and will be central to the development of our fast covariance estimation method. In polar coordinates in the unit disk the Fourier-Bessel basis functions are defined by
where is a normalization constant, is the -th order Bessel function of the first kind [6, §10.2], and is its -th smallest positive root; the indices run over .
Recent work  has devised a new fast algorithm to expand -images into Fourier-Bessel basis functions. Informally speaking, given basis coefficients, the algorithm can evaluate the function on an grid in operations; the adjoint can be computed in the same number of operations, which makes iterative methods fast.
2.4. Key property of Fourier-Bessel basis
A key property of the Fourier-Bessel basis is that convolution with radial functions are diagonal transformations in any truncated basis expansion. More precisely, the following result holds [21, Lemma 2.3]: suppose that for some index set , and is a radial function. Then,
where denotes orthogonal projection onto the span of , is the Fourier transform of , and is the -th positive root of .
We emphasize that the weights of the diagonal transform in (6) are not the coefficients of in the disk harmonic expansion. Indeed, since is radial, it follows from (5) that the coefficients of in the basis satisfy when . Computing the weights from the coefficients of in the basis, would involve computing weighted sums of the Fourier transforms of the basis functions: , for some index set .
As an alternative to the disk harmonic basis expansion, one can consider simply taking the Fourier-transform of (1), which also leads to a diagonal representation of the convolution operator. However, the discrete Fourier transform does not have the steerability property, which is essential for the covariance estimation. Another attempt could be to use the polar Fourier transform. However, this representation is not invariant to arbitrary in-plane rotations, but only to finitely many rotations as determined by the discretization spacing of the grid. These expansions are therefore unsuitable for the goal of this article and we instead use expansions into the Fourier-Bessel basis, although other steerable bases could be considered [14, 15]. Table 1 summarizes the considerations that make the Fourier-Bessel basis a natural choice of basis.
|Basis||Orth.||Cont. steerable||Radial convolution diag.||Fast expansion222The basis expansion from Cartesian grid representation can be completed within operations up to log factors.|
|2-D discrete Fourier||✓||✗||✓||✓|
|Polar Fourier||✗||✗ (discrete)||✓||✓|
|PSWF [14, 15]||✓||✓||✗||✗|
|Fourier-Bessel 333The new expansion algorithm  improves the previous computational method  in terms of accuracy guarantees, computational complexity and the fact that it derives weights such that radial convolution is a diagonal operation.||✓||✓||✓||✓|
2.5. Block diagonal structure
The steerable property of the Fourier-Bessel basis implies that the 2D-covariance will be block diagonal in this basis. A full representation of an matrix requires elements, but this is reduced to nonzero entries by the block diagonal structure. This block diagonal structure follows from the form of the basis functions and that the distribution of in-plane rotations is assumed to be uniform. Indeed, suppose that For simplicity assume that has mean zero; subtracting the mean will only change the radial components, since the other components corresponding to non-vanishing angular frequencies have zero mean by merely averaging over all possible in-plane rotations. By assumption 2, the covariance function in polar coordinates satisfies
for all in . The covariance function in (3) can be expanded in a double Fourier-Bessel basis expansion as
where is a Dirac function that is equal to if and zero otherwise, and we note that the second to last equality uses the fact that . Since the coefficients in the expansion (8) are unique, it follows from (9) that when . Hence, the covariance matrix has a block diagonal structure whose blocks consist of the indices for a given value of . In the following section, we show how these properties enable a fast method to estimate the covariance matrix.
2.6. Covariance Estimation
In the Fourier-Bessel basis, (1) is written as
where , , and are coefficient vectors of in the Fourier-Bessel basis, respectively and is the vector encoding the convolution operator of §2.4, i.e., with components . The vectors are -dimensional column vectors, where is the number of basis coefficients. We use this simple structure to obtain a closed form expression for the sample covariance matrix of the . We estimate this matrix by minimizing the discrepancy between the sample covariance and the population covariance; more precisely, the estimated covariance matrix is computed by solving the least squares-problem
where is defined by
whose solution is
The least squares-solution of (11) can be determined by the following system of linear equations
It follows that
As discussed in §2.4, the covariance matrix is block diagonal in the Fourier-Bessel basis. More precisely, the only non-zero elements of the matrix are those with . Therefore, the matrices in (15) and (16) need only be calculated for this subset of indices. Since there is a total of of these indices, this reduces the computational complexity compared to computing with the full matrices.
Note that the covariance matrix estimated from (16) may not be positive semidefinite due to subtraction of the term
. Therefore, when running the method in practice it is beneficial to use an eigenvalue shrinkage method. The computational cost of eigenvalue shrinkage for a matrix with our block structure is, see for example [7, 4]. For completeness, we include this computational cost of eigenvalue shrinkage in our overall computational complexity. Informally speaking, the idea of eigenvalue shrinkage is to replace the term in (16) by , and then shrink and truncate the eigenvalues in a systematic way before Hadamard division, see  for details. The steps of the algorithm are summarized in Algorithm 1.
The complexity of Step 1 of the algorithm is . The complexity of Step 1 is . The complexity of Step 1 is , since the number of basis coefficients is . The complexity of Step 1 is , since there are at most non-zero elements with indices , The complexity of Step 1 is , where the additional term comes from the computational complexity of eigenvalue shrinkage. Thus, the total complexity of the algorithm is .
2.7. Application to image denoising
3. Synthetic data results
We compare the timings of our fast method to previous approaches , for synthetic images generated from the 3-D volume of SARS-CoV-2 (Omicron) spike complexes  (EMD-32743), from the online EM data bank . The original volume has size in each dimension, with pixel size 0.832 Å. We downsample the original volume to size , with
and 512, respectively, and show the computational times. To generate the synthetic noisy images, we first generate 10000 clean projection images of the 3-D volume from random and uniformly distributed viewing directions. We next divide the set of clean images into a number of defocus groups, where the defocus values range from 1to 4 . For all CTFs, we set the voltage as 300 and the spherical aberration as 2 . After convolving the images with their CTFs, we add colored noise with power spectral density up to a constant scale, where is the radial frequency. For both the previous method and ours, the CTFs and the noisy images are whitened before estimating the covariance. A few sample clean and noisy images are shown in Figure 5. All experiments were carried out on a machine with 72 cores. We note that our implementation does not take advantage of all cores since some of the packages we use are not fully optimized for parallel computing; improving this is a technical direction for future work.
Figure 2 shows the time required to estimate the covariance matrices as a function of the number of defocus groups. Note that the runtime for the previous approach for the largest image and defocus group sizes is infeasibly large and that our fast method exhibits a speedup of up to three orders of magnitude.
Figure 3 shows the top six principal components estimated by our method and by traditional PCA using raw images, where , SNR and defocus groups, compared to traditional PCA on images. We use the sample covariance matrix of phase-flipped images in real space for the traditional PCA. For all methods, we use to denote times the eigenvalues of the eigenimages. The eigenimages from the traditional PCA look much noisier than ours, and contain artifacts that are due to imperfect CTF correction (see e.g. the circular artifacts in top 3 eigenimages in Figure 3). The eigenimages from the traditional PCA also fail to preserve the symmetries (see e.g. 6-th eigenimage in Figure 3) that are present in our eigenimages, since they do not utilize the steerable basis and rotation-augmented images.
Figure 4 shows the quality of covariance estimation and image denoising when and using defocus groups. The quality of the covariance estimation is measured by the relative error in each angular frequency, which is defined as
where and are respectively the -th diagonal blocks of the clean and estimated covariance matrix, corresponding to all indices of the form . The clean covariance matrix was approximated by the sample covariance matrix of clean projection images. The performance of image denoising is measured by the Fourier ring correlation (FRC) between the clean and denoised images. Namely, for the -th pair of clean and denoised images and , we first compute their Fourier coefficient vectors at radial frequency by the nonuniform FFT [8, 9, 18], where and . We then compute their averaged correlation for each :
where denotes the inner product between two complex vectors. The FRC is a real-valued quantity due to a symmetry property that arises since the images are real-valued.
In addition to the speedups of Figure 2, the left panel of Figure 4 demonstrates a slight increase in the estimation quality of our proposed algorithm, compared to the previous approach. This is possibly caused by improved accuracy in the Fourier-Bessel basis approximation as well as improved accuracy by using the closed-form expression (16) compared to the approximate conjugate gradient step of previous approaches. Similarly, on the right panel, the Fourier ring correlation between the denoised and clean images shows a slight performance increase. Figure 5 shows sample denoised images for different values of the signal-to-noise ratio (SNR) where and . As a comparison, we show images denoised using the approach of this paper and the CWF method .
4. Experimental data results
We conclude by using our method on two experimental datasets obtained from the Electron Microscopy Public Image Archive , namely EMPIAR-10028 , and EMPIAR-10081 . EMPIAR-10028 is a dataset of the Plasmodium falciparum 80S ribosome bound to the anti-protozoan drug emetine whose 3-D reconstruction is available in the EM data bank as EMD-2660 . The dataset contains 105247 motion corrected and picked particle images, from 1081 defocus groups, of size with 1.34 Å pixel size. EMPIAR-10081 is a dataset of the human HCN1 hyperpolarization-activated cyclic nucleotide-gated ion channel, whose 3-D reconstruction can be found in the EM data bank as EMD-8511 . The dataset contains 55870 motion corrected and picked particle images, from 53384 defocus groups, of size with 1.3 Å pixel size.
Computational times are shown in Table 2, showing a speedup of more than two orders of magnitude for the datasets with the largest number of distinct CTFs. Note that regular CWF on EMPIAR-10081 encounters memory issues and cannot be run to completion, whereas our fast method runs seamlessly. In order to obtain a comparison, we therefore additionally downsample these images to where the original CWF can successfully run. On EMPIAR-10028, we used all images for covariance estimation, and denoised the 2014 images from -th, -th, -th, …, -th defocus groups. On EMPIAR-10081, we used all images for covariance estimation, and denoised the 502 images from -th, -th, -th, …, -th defocus groups. For both old and new methods, the covariance matrices are further refined to correct image contrast variations . Sample visualization results are shown in Figures 6 and 7.
Covariance estimation and PCA of cryo-EM images are key ingredients in many classic cryo-EM methods including multivariate statistical analysis [27, 26, 28, 29] and Kam’s method for ab-initio modeling . We propose a fast method to estimate the covariance matrix of noisy cryo-EM images, and then illustrate its application to simultaneously correct for the CTFs and denoise the images. The approach relies on recent improvements to algorithms for expanding images in the Fourier-Bessel basis , and has time complexity which is independent of the number of defocus groups. Our new approach is both significantly faster and more memory-efficient compared to the previous CWF method  and we apply our method to large experimental datasets with many distinct CTFs with speedups by factors up to more than two orders of magnitude. Our approach could potentially be extended to higher-dimensional data and to the setting where images are distorted by CTFs which are not exactly radial, using either analytical correction terms or iterative numerical steps.
-  A. S. Bandeira, B. Blum-Smith, J. Kileel, A. Perry, J. Weed, and A. S. Wein. Estimation under group actions: recovering orbits from invariants. arXiv preprint arXiv:1712.10163, 2017.
-  T. Bendory, Y. Khoo, J. Kileel, O. Mickelin, and A. Singer. Autocorrelation analysis for cryo-EM with sparsity constraints: Improved sample complexity and projection-based algorithms. arXiv preprint arXiv:2209.10531, 2022.
-  T. Bhamre, T. Zhang, and A. Singer. Orthogonal matrix retrieval in cryo-electron microscopy. In 2015 IEEE 12th International Symposium on Biomedical Imaging (ISBI), pages 1048–1052. IEEE, 2015.
-  T. Bhamre, T. Zhang, and A. Singer. Denoising and covariance estimation of single particle cryo-EM images. Journal of structural biology, 195(1):72–81, 2016.
-  T. Bhamre, T. Zhang, and A. Singer. Anisotropic twicing for single particle reconstruction using autocorrelation analysis. arXiv preprint arXiv:1704.07969, 2017.
-  NIST Digital Library of Mathematical Functions. http://dlmf.nist.gov/, Release 1.1.5 of 2022-03-15. F. W. J. Olver, A. B. Olde Daalhuis, D. W. Lozier, B. I. Schneider, R. F. Boisvert, C. W. Clark, B. R. Miller, B. V. Saunders, H. S. Cohl, and M. A. McClain, eds.
-  D. L. Donoho, M. Gavish, and I. M. Johnstone. Optimal shrinkage of eigenvalues in the spiked covariance model. Annals of statistics, 46(4):1742, 2018.
-  A. Dutt and V. Rokhlin. Fast fourier transforms for nonequispaced data. SIAM Journal on Scientific computing, 14(6):1368–1393, 1993.
-  L. Greengard and J.-Y. Lee. Accelerating the nonuniform fast Fourier transform. SIAM Review, 46(3):443–454, 2004.
-  H. Guo, Y. Gao, T. Li, T. Li, Y. Lu, L. Zheng, Y. Liu, T. Yang, F. Luo, S. Song, et al. Structures of omicron spike complexes and implications for neutralizing antibody development. Cell reports, 39(5):110770, 2022.
-  S. Huang, M. Zehni, I. Dokmanić, and Z. Zhao. Orthogonal matrix retrieval with spatial consensus for 3D unknown-view tomography. arXiv preprint arXiv:2207.02985, 2022.
-  A. Iudin, P. K. Korir, J. Salavert-Torres, G. J. Kleywegt, and A. Patwardhan. EMPIAR: a public archive for raw electron microscopy image data. Nature methods, 13(5):387–388, 2016.
-  Z. Kam. The reconstruction of structure from electron micrographs of randomly oriented particles. Journal of Theoretical Biology, 82(1):15–39, 1980.
-  B. Landa and Y. Shkolnisky. Approximation scheme for essentially bandlimited and space-concentrated functions on a disk. Applied and Computational Harmonic Analysis, 43(3):381–403, 2017.
-  B. Landa and Y. Shkolnisky. Steerable principal components for space-frequency localized images. SIAM Journal on Imaging Sciences, 10(2):508–534, 2017.
-  C. L. Lawson, A. Patwardhan, M. L. Baker, C. Hryc, E. S. Garcia, B. P. Hudson, I. Lagerstedt, S. J. Ludtke, G. Pintilie, R. Sala, et al. EMDataBank unified data resource for 3DEM. Nucleic Acids Research, 44(D1):D396–D403, 2016.
-  C.-H. Lee and R. MacKinnon. Structures of the human HCN1 hyperpolarization-activated channel. Cell, 168(1-2):111–120, 2017.
-  J.-Y. Lee and L. Greengard. The type 3 nonuniform FFT and its applications. Journal of Computational Physics, 206(1):1–5, 2005.
-  E. Levin, T. Bendory, N. Boumal, J. Kileel, and A. Singer. 3D ab initio modeling in cryo-EM by autocorrelation analysis. In 2018 IEEE 15th International Symposium on Biomedical Imaging (ISBI 2018), pages 1569–1573. IEEE, 2018.
W. Liu and J. Frank.
Estimation of variance distribution in three-dimensional reconstruction. i. theory.JOSA A, 12(12):2615–2627, 1995.
-  N. F. Marshall, O. Mickelin, and A. Singer. Fast expansion into harmonics on the disk: a steerable basis with fast radial convolutions. arXiv preprint arXiv:2207.13674, 2022.
-  P. A. Penczek, M. Kimmel, and C. M. Spahn. Identifying conformational states of macromolecules by eigen-analysis of resampled cryo-em images. Structure, 19(11):1582–1590, 2011.
-  A. Rangan, M. Spivak, J. Andén, and A. Barnett. Factorization of the translation kernel for fast rigid image alignment. Inverse Problems, 36(2):024001, 2020.
-  N. Sharon, J. Kileel, Y. Khoo, B. Landa, and A. Singer. Method of moments for 3D single particle ab initio modeling with non-uniform distribution of viewing angles. Inverse Problems, 36(4):044003, 2020.
-  Y. Shi and A. Singer. Ab-initio Contrast Estimation and Denoising of Cryo-EM Images. Computer Methods and Programs in Biomedicine, 224:107018, 2022.
M. Van Heel.
Multivariate statistical classification of noisy images (randomly oriented biological macromolecules).Ultramicroscopy, 13(1-2):165–183, 1984.
-  M. Van Heel and J. Frank. Use of multivariates statistics in analysing the images of biological macromolecules. Ultramicroscopy, 6(1):187–194, 1981.
-  M. van Heel, R. Portugal, and M. Schatz. Multivariate statistical analysis in single particle (Cryo) electron microscopy. An electronic textbook: electron microscopy in life science. 3D-EM Network of Excellence, London, 2009.
-  M. van Heel, R. V. Portugal, M. Schatz, et al. Multivariate statistical analysis of large datasets: single particle electron microscopy. Open Journal of Statistics, 6(04):701, 2016.
-  W. Wong, X.-c. Bai, A. Brown, I. S. Fernandez, E. Hanssen, M. Condron, Y. H. Tan, J. Baum, and S. H. Scheres. Cryo-EM structure of the Plasmodium falciparum 80S ribosome bound to the anti-protozoan drug emetine. Elife, 3:e03080, 2014.
-  Z. Zhao, Y. Shkolnisky, and A. Singer. Fast steerable principal component analysis. IEEE Transactions on Computational Imaging, 2(1):1–12, 2016.
-  Z. Zhao and A. Singer. Rotationally invariant image representation for viewing direction classification in cryo-EM. Journal of structural biology, 186(1):153–166, 2014.