Bias and variance reduction and denoising for CTF Estimation

Ayelet Heimowitz; Joakim And\'en; Amit Singer

arXiv:1908.03454·eess.IV·January 29, 2020

Bias and variance reduction and denoising for CTF Estimation

Ayelet Heimowitz, Joakim And\'en, Amit Singer

PDF

Open Access

TL;DR

This paper introduces a novel CTF estimation method for electron microscopy images that reduces bias and variance using multi-taper spectral estimation and leverages known properties of the CTF to improve accuracy.

Contribution

The proposed method combines multi-taper spectral estimation with known CTF properties to enhance the accuracy of CTF estimation in electron microscopy images.

Findings

01

Captures zero-crossings of the CTF in low-mid frequencies

02

Reduces bias and variance in spectral estimates

03

Improves CTF estimation accuracy

Abstract

When using an electron microscope for imaging of particles embedded in vitreous ice, the objective lens will inevitably corrupt the projection images. This corruption manifests as a band-pass filter on the micrograph. In addition, it causes the phase of several frequency bands to be flipped and distorts frequency bands. As a precursor to compensating for this distortion, the corrupting point spread function, which is termed the contrast transfer function (CTF) in reciprocal space, must be estimated. In this paper, we will present a novel method for CTF estimation. Our method is based on the multi-taper method for power spectral density estimation, which aims to reduce the bias and variance of the estimator. Furthermore, we use known properties of the CTF and of the background of the power spectrum to increase the accuracy of our estimation. We will show that the resulting estimates…

Tables3

Table 1. Table 1: Description of the EMPIAR-10002, EMPIAR-10028, EMPIAR-10242 and EMPIAR-10249 datasets.

Dataset	Molecule	Pixel size (Å)	Spherical	Voltage (kV)	Microscope	Detector
			aberration
EMPIAR-10002	80S ribosome	1.77	2.0	300	Polara	Falcon I
EMPIAR-10028	80S ribosome	1.34	2.0	300	Polara	Falcon II
EMPIAR-10242	2N3R tau filaments	1.04	2.7	300	Titan Krios	Gatan K2 Summit
EMPIAR-10249	HLA dehydrogenase	0.56	2.7	200	Talos Arctica	Gatan K2 Summit

Table 2. Table 2: Comparison between parameters estimated by ASPIRE-CTF, CTFFIND4 and Gctf. We present the mean (over each dataset) of normalized differences between each two CTF estimation methods as detailed in ( 31 ). The subscripts a 𝑎 a , g 𝑔 g and c 𝑐 c indicate ASPIRE-CTF, Gctf and CTFFIND4, respectively.

Dataset	Molecule	$ϵ_{a, g} (Δ f_{1})$	$ϵ_{a, c} (Δ f_{1})$	$ϵ_{c, g} (Δ f_{1})$	$ϵ_{a, g} (Δ f_{2})$	$ϵ_{a, c} (Δ f_{2})$	$ϵ_{c, g} (Δ f_{2})$	$ϵ_{a, g} (α_{f})$	$ϵ_{a, c} (α_{f})$	$ϵ_{c, g} (α_{f})$
001	GroEL	0.0157	0.2370	0.2364	0.0155	0.1282	0.6081	1.6599	2.1429	1.3848
002	GroEL	1.3236	1.6500	0.1515	1.4876	0.9746	0.2374	3.6137	5.3979	1.9592
003	60S ribosome	0.0029	0.0049	0.0029	0.0030	0.0032	0.0017	6.1042	6.6618	2.9159
004	60S ribosome	0.0051	0.0130	0.0098	0.0057	0.0406	0.2990	4.2653	2.3667	2.6050
005	apoferritin	0.0036	0.1744	0.8809	0.0036	0.1997	1.5230	0.6597	0.8927	0.7413
006	apoferritin	0.0099	0.2094	0.6392	0.0069	0.3409	1.5217	0.2670	1.3014	1.6142
007	TMV virus	0.0224	0.2144	0.5163	0.0265	0.2674	1.1320	0.4254	0.9891	0.8985
008	TMV virus	0.1961	0.6043	0.3988	0.2163	0.3016	0.5707	3.0253	1.8514	2.8394

Table 3. Table 3: Comparison between parameters estimated by ASPIRE-CTF, CTFFIND4 and Gctf. We present the standard deviation (over each dataset) of normalized differences between each two CTF estimation methods as detailed in ( 31 ). The subscripts a 𝑎 a , g 𝑔 g and c 𝑐 c indicate ASPIRE-CTF, Gctf and CTFFIND4, respectively.

Dataset	Molecule	$ϵ_{a, g} (Δ f_{1})$	$ϵ_{a, c} (Δ f_{1})$	$ϵ_{c, g} (Δ f_{1})$	$ϵ_{a, g} (Δ f_{2})$	$ϵ_{a, c} (Δ f_{2})$	$ϵ_{c, g} (Δ f_{2})$	$ϵ_{a, g} (α_{f})$	$ϵ_{a, c} (α_{f})$	$ϵ_{c, g} (α_{f})$
001	GroEL	0.0144	0.4710	0.5276	0.0267	0.2489	1.9206	3.1585	2.0435	2.2632
002	GroEL	2.9445	2.2763	0.2655	3.3142	1.9597	0.5196	9.2965	15.0486	2.7796
003	60S ribosome	0.0022	0.0031	0.0036	0.0031	0.0047	0.0020	20.8434	22.3327	15.8681
004	60S ribosome	0.0028	0.0328	0.0387	0.0039	0.1778	1.4504	9.4786	3.5234	3.9324
005	apoferritin	0.0037	0.3265	1.8456	0.0029	0.3591	3.1834	0.9381	0.7127	0.6918
006	apoferritin	0.0126	0.2902	1.9539	0.0064	0.3330	2.7270	0.3485	1.1711	2.7480
007	TMV virus	0.0192	0.2698	0.8581	0.0172	0.3335	2.0805	0.6437	0.9296	0.9966
008	TMV virus	0.8337	0.7936	0.8543	0.8947	0.4088	0.9907	4.7943	2.6911	6.1583

Equations77

S_{x} (g) = n_{1} = - \infty \sum \infty n_{2} = - \infty \sum \infty R_{x} [m, n] e^{j 2 π (g_{1} n_{1} + g_{2} n_{2})},

S_{x} (g) = n_{1} = - \infty \sum \infty n_{2} = - \infty \sum \infty R_{x} [m, n] e^{j 2 π (g_{1} n_{1} + g_{2} n_{2})},

y = h_{ϕ} * x + e,

y = h_{ϕ} * x + e,

H_{ϕ} (g) = - sin (χ_{ϕ} (g)),

H_{ϕ} (g) = - sin (χ_{ϕ} (g)),

χ_{ϕ} (g) = \frac{1}{2 p ^{2}} π λ r^{2} Δ f_{ϕ} (α) - \frac{1}{2 p ^{4}} π λ^{3} r^{4} C_{s} + w,

χ_{ϕ} (g) = \frac{1}{2 p ^{2}} π λ r^{2} Δ f_{ϕ} (α) - \frac{1}{2 p ^{4}} π λ^{3} r^{4} C_{s} + w,

Δ f_{ϕ} (α) = Δ f_{1} + Δ f_{2} + (Δ f_{1} - Δ f_{2}) cos (2 α - 2 α_{f}),

Δ f_{ϕ} (α) = Δ f_{1} + Δ f_{2} + (Δ f_{1} - Δ f_{2}) cos (2 α - 2 α_{f}),

S_{y} (g) = ∣ H_{ϕ} (g) ∣^{2} S_{x} (g) + S_{e} (g) .

S_{y} (g) = ∣ H_{ϕ} (g) ∣^{2} S_{x} (g) + S_{e} (g) .

\hat{S}_{y}^{(p)} (g) = \frac{1}{N ^{2}} k_{1}, k_{2} = 0 \sum N - 1 y [k_{1}, k_{2}] e^{- j 2 π (g_{1} k_{1} + g_{2} k_{2})}^{2},

\hat{S}_{y}^{(p)} (g) = \frac{1}{N ^{2}} k_{1}, k_{2} = 0 \sum N - 1 y [k_{1}, k_{2}] e^{- j 2 π (g_{1} k_{1} + g_{2} k_{2})}^{2},

M_{N} = {- \frac{1}{2}, - \frac{1}{2} + \frac{2}{N}, \dots, \frac{1}{2} - \frac{2}{N}}^{2} .

M_{N} = {- \frac{1}{2}, - \frac{1}{2} + \frac{2}{N}, \dots, \frac{1}{2} - \frac{2}{N}}^{2} .

MSE (\hat{S}_{y}^{(p)} (g)) = E [∣ \hat{S}_{y}^{(p)} (g) - S_{y} (g) ∣^{2}] .

MSE (\hat{S}_{y}^{(p)} (g)) = E [∣ \hat{S}_{y}^{(p)} (g) - S_{y} (g) ∣^{2}] .

Bias (\hat{S}_{y}^{(p)} (g)) = E [\hat{S}_{y}^{(p)} (g)] - S_{y} (g)

Bias (\hat{S}_{y}^{(p)} (g)) = E [\hat{S}_{y}^{(p)} (g)] - S_{y} (g)

Var (\hat{S}_{y}^{(p)} (g)) = E [\hat{S}_{y}^{(p)} (g) - E [\hat{S}_{y}^{(p)} (g)]^{2}]

Var (\hat{S}_{y}^{(p)} (g)) = E [\hat{S}_{y}^{(p)} (g) - E [\hat{S}_{y}^{(p)} (g)]^{2}]

MSE (\hat{S}_{y}^{(p)} (g)) = Bias^{2} (\hat{S}_{y}^{(p)} (g)) + Var (\hat{S}_{y}^{(p)} (g)) .

MSE (\hat{S}_{y}^{(p)} (g)) = Bias^{2} (\hat{S}_{y}^{(p)} (g)) + Var (\hat{S}_{y}^{(p)} (g)) .

\hat{S}_{y}^{(b)} (g) = \frac{1}{B} b = 0 \sum B - 1 \hat{S}_{y_{b}}^{(p)} (g) .

\hat{S}_{y}^{(b)} (g) = \frac{1}{B} b = 0 \sum B - 1 \hat{S}_{y_{b}}^{(p)} (g) .

\hat{S}_{y}^{(w)} (g) = \frac{1}{B} b = 1 \sum B \hat{S}_{y_{b} \cdot w}^{(p)} (g)

\hat{S}_{y}^{(w)} (g) = \frac{1}{B} b = 1 \sum B \hat{S}_{y_{b} \cdot w}^{(p)} (g)

\hat{S}_{y}^{(mt)} (g) = \frac{1}{L B} b = 0 \sum B - 1 ℓ = 0 \sum L - 1 \hat{S}_{y_{b} \cdot w_{ℓ}}^{(p)} (g),

\hat{S}_{y}^{(mt)} (g) = \frac{1}{L B} b = 0 \sum B - 1 ℓ = 0 \sum L - 1 \hat{S}_{y_{b} \cdot w_{ℓ}}^{(p)} (g),

S_{y} (g) - S_{e} (g) = ∣ H_{ϕ} (g) ∣^{2} S_{x} (g) .

S_{y} (g) - S_{e} (g) = ∣ H_{ϕ} (g) ∣^{2} S_{x} (g) .

\begin{array}[]{ll}\displaystyle\operatorname*{minimize}_{\hat{S}_{\mathbf{e}}}&\underset{r=0,\frac{1}{K},\dots,\frac{m}{K}}{\sum}{\hat{S}^{(\mathrm{mt})}_{\mathbf{y}}\left(r\right)}-\hat{S}_{\mathbf{e}}\left(r\right)\\ \text{subject to}&\hat{S}_{\mathbf{e}}(r)\leq{\hat{S}^{(\mathrm{mt})}_{\mathbf{y}}\left(r\right)},\quad r=0,\dots,\frac{m}{K}\\ &\hat{S}_{\mathbf{e}}\left(r+1\right)+\hat{S}_{\mathbf{e}}\left(r-1\right)\geq 2\hat{S}_{\mathbf{e}}\left(r\right),r=1,\dots,\frac{m}{K}\\ &\hat{S}_{\mathbf{e}}\left(r\right)\geq 0,\quad r=0,\dots,\frac{m}{K},\end{array}

\begin{array}[]{ll}\displaystyle\operatorname*{minimize}_{\hat{S}_{\mathbf{e}}}&\underset{r=0,\frac{1}{K},\dots,\frac{m}{K}}{\sum}{\hat{S}^{(\mathrm{mt})}_{\mathbf{y}}\left(r\right)}-\hat{S}_{\mathbf{e}}\left(r\right)\\ \text{subject to}&\hat{S}_{\mathbf{e}}(r)\leq{\hat{S}^{(\mathrm{mt})}_{\mathbf{y}}\left(r\right)},\quad r=0,\dots,\frac{m}{K}\\ &\hat{S}_{\mathbf{e}}\left(r+1\right)+\hat{S}_{\mathbf{e}}\left(r-1\right)\geq 2\hat{S}_{\mathbf{e}}\left(r\right),r=1,\dots,\frac{m}{K}\\ &\hat{S}_{\mathbf{e}}\left(r\right)\geq 0,\quad r=0,\dots,\frac{m}{K},\end{array}

x (r, α) = k = - \infty \sum \infty q = 0 \sum \infty a_{k, q} f_{k, q} (r) e^{j k α},

x (r, α) = k = - \infty \sum \infty q = 0 \sum \infty a_{k, q} f_{k, q} (r) e^{j k α},

H_{ϕ} (g) = n = 0 n even \sum P \frac{( - 1 ) ^{\frac{n}{2} + 1}}{n !} sin (χ_{ϕ}^{0} (r)) C_{n, ϕ} (g) + n = 1 n odd \sum P \frac{( - 1 ) ^{\frac{n + 1}{2}}}{n !} cos (χ_{ϕ}^{0} (r)) C_{n, ϕ} (g) + R_{P} (g),

H_{ϕ} (g) = n = 0 n even \sum P \frac{( - 1 ) ^{\frac{n}{2} + 1}}{n !} sin (χ_{ϕ}^{0} (r)) C_{n, ϕ} (g) + n = 1 n odd \sum P \frac{( - 1 ) ^{\frac{n + 1}{2}}}{n !} cos (χ_{ϕ}^{0} (r)) C_{n, ϕ} (g) + R_{P} (g),

C_{n, ϕ} (g) = (\frac{1}{2} π λ (Δ f_{1} - Δ f_{2}) cos (2 (α - α_{f})) \frac{r ^{2}}{p ^{2}})^{n},

C_{n, ϕ} (g) = (\frac{1}{2} π λ (Δ f_{1} - Δ f_{2}) cos (2 (α - α_{f})) \frac{r ^{2}}{p ^{2}})^{n},

χ_{ϕ}^{0} (r) = \frac{1}{2} π λ r^{2} (Δ f_{1} + Δ f_{2}) - \frac{1}{2} π λ^{3} r^{4} C_{s} + w

χ_{ϕ}^{0} (r) = \frac{1}{2} π λ r^{2} (Δ f_{1} + Δ f_{2}) - \frac{1}{2} π λ^{3} r^{4} C_{s} + w

(\frac{Δ f _{1} - Δ f _{2}}{Δ f _{1} + Δ f _{2}})^{(P + 1)}

(\frac{Δ f _{1} - Δ f _{2}}{Δ f _{1} + Δ f _{2}})^{(P + 1)}

H_{ϕ} (g) \approx - sin (χ_{ϕ}^{0} (r)) - cos (χ_{ϕ}^{0} (r)) C_{1, ϕ} (g)

H_{ϕ} (g) \approx - sin (χ_{ϕ}^{0} (r)) - cos (χ_{ϕ}^{0} (r)) C_{1, ϕ} (g)

H_{ϕ} (g) \approx - sin (χ_{ϕ}^{0} (r)) - \frac{1}{4 p ^{2}} cos (χ_{ϕ}^{0} (r)) π λ (Δ f_{1} - Δ f_{2}) e^{- j 2 α_{f}} r^{2} e^{j 2 α} - \frac{1}{4 p ^{2}} cos (χ_{ϕ}^{0} (r)) π λ (Δ f_{1} - Δ f_{2}) e^{j 2 α_{f}} r^{2} e^{- j 2 α} .

H_{ϕ} (g) \approx - sin (χ_{ϕ}^{0} (r)) - \frac{1}{4 p ^{2}} cos (χ_{ϕ}^{0} (r)) π λ (Δ f_{1} - Δ f_{2}) e^{- j 2 α_{f}} r^{2} e^{j 2 α} - \frac{1}{4 p ^{2}} cos (χ_{ϕ}^{0} (r)) π λ (Δ f_{1} - Δ f_{2}) e^{j 2 α_{f}} r^{2} e^{- j 2 α} .

a_{k, q} = \frac{1}{K ^{2}} g \in M_{K} \sum \hat{S}_{z}^{(lp)} (g) f_{k, q} (r) e^{j k α},

a_{k, q} = \frac{1}{K ^{2}} g \in M_{K} \sum \hat{S}_{z}^{(lp)} (g) f_{k, q} (r) e^{j k α},

P_{cc} (ϕ) = \frac{\sum _{g \in R} ∣ H _{ϕ} ( g ) ∣ S ^ _{z}^{1/2} ( g )}{( \sum _{g \in R} ∣ H _{ϕ} ( g ) ∣ ^{2} \sum _{g \in R} S ^ _{z} ( g ) ) ^{1/2}},

P_{cc} (ϕ) = \frac{\sum _{g \in R} ∣ H _{ϕ} ( g ) ∣ S ^ _{z}^{1/2} ( g )}{( \sum _{g \in R} ∣ H _{ϕ} ( g ) ∣ ^{2} \sum _{g \in R} S ^ _{z} ( g ) ) ^{1/2}},

R = {m_{1}, m_{1} + \frac{2}{N}, \dots, \frac{3}{8}} \times {0},

R = {m_{1}, m_{1} + \frac{2}{N}, \dots, \frac{3}{8}} \times {0},

M_{1, 1} = \sum_{g \in M_{K}} g_{1}^{2} \hat{S}_{z}^{1/2} (g), M_{1, 2} = M_{2, 1} = \sum_{g \in M_{K}} g_{1} g_{2} \hat{S}_{z}^{1/2} (g), M_{2, 2} = \sum_{g \in M_{K}} g_{2}^{2} \hat{S}_{z}^{1/2} (g) .

M_{1, 1} = \sum_{g \in M_{K}} g_{1}^{2} \hat{S}_{z}^{1/2} (g), M_{1, 2} = M_{2, 1} = \sum_{g \in M_{K}} g_{1} g_{2} \hat{S}_{z}^{1/2} (g), M_{2, 2} = \sum_{g \in M_{K}} g_{2}^{2} \hat{S}_{z}^{1/2} (g) .

\frac{1}{2} (Δ f_{1} + Δ f_{2})

\frac{1}{2} (Δ f_{1} + Δ f_{2})

\frac{Δ f _{1}}{Δ f _{2}}

Δ f_{1, ⋆}

Peer Reviews

No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.

Videos

No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.

Taxonomy

TopicsAdvanced Electron Microscopy Techniques and Applications · Electron and X-Ray Spectroscopy Techniques · Advanced Fluorescence Microscopy Techniques

Full text

Reducing Bias and Variance for CTF Estimation in Single Particle Cryo-EM

Ayelet Heimowitz

[email protected]

Joakim Andén

[email protected]

Amit Singer

[email protected]

The Program in Applied and Computational Mathematics, Princeton University, Princeton, NJ

Center for Computational Biology, Flatiron Institute, New York, NY

Department of Mathematics, Princeton University, Princeton, NJ

Abstract

When using an electron microscope for imaging of particles embedded in vitreous ice, the recorded image, or micrograph, is a significantly degraded version of the tomographic projection of the sample. Apart from noise, the image is affected by the optical configuration of the microscope. This transformation is typically modeled as a convolution with a point spread function. The Fourier transform of this function, known as the contrast transfer function (CTF), is oscillatory, attenuating and amplifying different frequency bands, and sometimes flipping their signs. High-resolution reconstruction requires this CTF to be accounted for, but as its form depends on experimental parameters, it must first be estimated from the micrograph. We present a new method for CTF estimation based on multitaper methods, which reduces bias and variance in the estimate. We also use known properties of the CTF and the background noise power spectrum to further reduce the variance through background subtraction and steerable basis projection. We show that the resulting power spectrum estimates better capture the zero-crossings of the CTF and yield accurate CTF estimates on several experimental micrographs.

keywords:

contrast transfer function, cryo-electron microscopy, linear programming, multitaper estimator, spectral estimation, steerable basis expansion

††journal: Ultramicroscopy

1 Introduction

In recent years, single particle cryo-electron microscopy (cryo-EM) has emerged as a leading tool for resolving the 3D structure of macromolecules from multiple 2D projections of a specimen [Cheng et al., 2017]. In this technique, multiple copies of a particle are embedded in vitreous ice and imaged in an electron microscope. This yields a set of micrographs, each containing several 2D particle projections.

The micrograph does not contain clean particle projections but is contaminated by several factors, including noise, ice aggregates and carbon film projection. The noise stems from an inherent limitation on the number of imaging electrons that can be applied to the specimen. The interference from carbon film and ice aggregates are due to the particular sample preparation techniques used.

The 2D projections in the micrograph are also distorted by convolution with a point spread function. This point spread function is due to the electron microscope configuration. It attenuates certain frequencies and flips the sign of certain frequency bands. A 3D density map reconstructed from distorted projections yields an unreliable representation of the particle [Frank, 1996]. It is therefore important to estimate the point spread function and account for it during reconstruction.

To estimate these parameters, it is convenient to consider the Fourier transform of the point spread function, known as the contrast transfer function (CTF). This is due to two factors. First, the CTF has a simple expression in the polar coordinates of the spatial frequency. Second, its effect is directly visible in the frequency domain where the CTF acts as a pointwise multiplication rather than a convolution [Erickson & Klug, 1971].

CTF estimation is one of the first steps in the single particle cryo-EM pipeline. Indeed, accounting for the CTF is needed in a variety of tasks, such as particle picking [Heimowitz et al., 2018], denoising [Bhamre et al., 2016], class averaging [Scheres, 2012], ab initio reconstruction [Punjani et al., 2017], refinement [Scheres, 2012, Punjani et al., 2017, Tang et al., 2007, Grant et al., 2018] and heterogeneity analysis [Scheres, 2012].

The CTF is typically modeled as a sine function whose argument depends on the spatial frequency and several parameters of the objective lens of the microscope [Frank, 1996]. The parameters we focus on in this paper are the defocus and astigmatism of the objective lens as these are unknown and must be estimated from the data. Additionally, the CTF is multiplied by a damping envelope, which suppresses the information in high frequencies [Sorzano et al., 2007].

When estimating the CTF parameters, it is common to first estimate the power spectrum of the micrograph. The observed micrograph image is typically modeled as a CTF-dependent term plus a noise term unaffected by CTF. The first term corresponds to a noiseless micrograph, that is, the tomographic projection of the sample filtered by the CTF Zhu et al. [1997]. Modeling the unfiltered and filtered micrographs as 2D random fields, we find that their power spectra are closely related: the latter equals the former multiplied by the squared CTF. This multiplication induces concentric rings, known as Thon rings [Thon, 1971], in the power spectrum of the filtered micrograph (see Fig. 2). Estimating the CTF therefore reduces to fitting the parameters of the CTF to the estimated power spectrum.

The vast majority of CTF estimation methods use a variant of the periodogram when estimating the power spectrum of the micrograph. This is due to its speed and simplicity. Unfortunately, the periodogram produces a biased and inconsistent estimate of the micrograph power spectrum.

Beyond these issues with the spectral estimators, fitting the CTF model to the estimated power spectrum is complicated by factors such as the high levels of noise present in the micrograph, coincidence loss at the detector and more. The expected power spectrum equals the power spectrum of the clean, filtered micrograph plus a background term caused by the aforementioned complications. The background masks the true oscillations of the power spectrum of the particle projection. It is therefore important to estimate and remove the background from the estimated power spectrum [Zhu et al., 1997, Rohou & Grigorieff, 2015, Mindell & Grigorieff, 2003, Zhang, 2016].

Assuming that the micrograph power spectrum and the background were both estimated perfectly, the background-subtracted power spectrum equals the power spectrum of the filtered, clean micrograph. One way to estimate the CTF parameters is then to maximize the correlation of the background-subtracted power spectrum estimate and a squared CTF (or some monotonic function thereof) [Rohou & Grigorieff, 2015, Mindell & Grigorieff, 2003, Zhang, 2016, Tani et al., 1996]. Optimizing the correlation then provides an estimate of the defocus and astigmatism. Another approach identifies a single ring in the estimated power spectrum and uses it to derive a closed-form solution of the CTF parameters Yan et al. [2017]. In order to formulate this solution, all prior knowledge regarding the spherical aberration must be ignored.

In this paper, we present ASPIRE-CTF, which is a new method for CTF estimation, available as part of the ASPIRE package.111https://github.com/ComputationalCryoEM/ASPIRE-Python We first estimate the power spectrum using a multitaper estimator [Babadi & Brown, 2014], further reducing the variance by averaging estimates from multiple regions of the micrograph. Using this estimated power spectrum, we estimate the background noise spectrum using linear programming (LP). Instead of using an approximate background model, our scheme ensures that the background-subtracted power spectrum estimate is non-negative and convex. We also show that the CTF is contained in the span of a small number of steerable basis functions. Thus, we further reduce the variance in our power spectrum estimate by projecting onto this span.

Given the power spectrum estimate, we provide two solutions for estimating the CTF parameters. Our first solution is similar to [Rohou & Grigorieff, 2015, Mindell & Grigorieff, 2003, Zhang, 2016, Tani et al., 1996], where CTF parameters are estimated by maximizing the correlation of the square root of the power spectrum estimate with the absolute value of simulated CTFs. The second solution uses the spatial frequencies of several zero-crossings. Since we expect these zero-crossings to coincide with those of the squared CTF, we use them to define an overdetermined system of equations over the CTF parameters that we then solve. We note that, while our first solution is more robust, our second solution is faster to compute.

Our method is experimentally verified in Section 3. This is done via a comparison of the defocus estimates with that of Rohou & Grigorieff [2015], Zhang [2016] on several datasets from the CTF challenge [Marabini et al., 2015]. We show that our power spectrum estimation method is usually in agreement with one of the state-of-the-art methods Rohou & Grigorieff [2015], Zhang [2016].

The main contribution of this paper, appearing in Sections 2.2.2-2.3, is our method for estimating the power spectrum of a micrograph. We reduce the variability of the power spectrum estimate, and are therefore the first to obtain an estimate where several zero-crossing rings of the CTF are easily recovered without additional assumptions.

We present the pipeline of our method in Fig. 1. For each step of our suggested framework, we refer the reader to the appropriate section of the paper.

Notation

Given a 2D stationary random field $\mathbf{x}$ defined over $\mathbb{Z}^{2}$ , we denote its autocovariance function by $R_{\mathbf{x}}$ . The Fourier transform of $R_{\mathbf{x}}$ is known as the power spectrum of $\mathbf{x}$ and is given by

[TABLE]

for $\mathbf{g}=(g_{1},g_{2})\in[-1/2,1/2]^{2}$ and $\mathrm{j}=\sqrt{-1}$ . We denote magnitude of the spatial frequency vector $\mathbf{g}$ by ${r}$ and its counterclockwise angle with the positive x-axis by ${\alpha}$ .

2 Materials and Methods

2.1 Problem formulation

In the sample preparation stage of the single particle cryo-EM pipeline, many copies of a particle are embedded in vitreous ice. The imaging process uses an electron microscope to obtain a micrograph containing 2D projections of each instance. Under the weak-phase object approximation, we may describe this process by the linear model [Frank, 1996, Thon, 1971, Mindell & Grigorieff, 2003]

[TABLE]

where the clean tomographic projection $\mathbf{x}$ and the additive noise $\mathbf{e}$ are modeled as 2D stationary random fields [Frank, 1996]. Since convolution preserves stationarity, the observed micrograph $\mathbf{y}$ is also a stationary random field. In this model, the clean projection $\mathbf{x}$ is convolved with the point spread function of the microscope $h_{\phi}$ which depends on a parameter vector $\phi$ . We will at times denote this clean, but filtered, micrograph by $\mathbf{z}=h_{\phi}*\mathbf{x}$ .

The CTF is the Fourier transform of the point spread function $H_{\phi}$ and may be modeled by [Rohou & Grigorieff, 2015]

[TABLE]

where $\mathbf{g}$ is the spatial frequency. Its phase is given by

[TABLE]

where $\lambda$ is the electron wavelength, $C_{s}$ is the spherical aberration, $w$ is the amplitude contrast, and $p$ is the pixel size. We also have the astigmatic defocus depth

[TABLE]

where ${\alpha}$ is the polar angle of $\mathbf{g}$ and $\Delta f_{1}$ , $\Delta f_{2}$ , and $\alpha_{f}$ are the major and minor defocus depths and the defocus angle, respectively. These together form the defocus vector $\phi=(\Delta f_{1},\Delta f_{2},\alpha_{f})$ , which parametrizes the CTF. The values $\Delta f_{1}$ and $\Delta f_{2}$ determine the amount of defocus along two perpendicular axes, while $\alpha_{f}$ specifies the counterclockwise angle between the major defocus axis and the positive x-axis. The difference $\Delta f_{1}-\Delta f_{2}$ measures the amount of astigmatism in the CTF. A visualization of the effect of astigmatism is provided in Fig. 2.

The model (3) allows us to discern several properties of the CTF. First, $H_{\phi}$ is real and oscillates between positive and negative values. As a result, it has several zero crossings. Second, the CTF is radially symmetric when $\Delta f_{1}=\Delta f_{2}$ (the non-astigmatic case). Third, with no spherical aberration (i.e., $C_{s}=0$ ) the level sets of the CTF consist of ellipses centered at the origin. The spherical aberration $C_{s}$ thus accounts for small deviations from the elliptical shape.

While the parameters $\lambda$ , $C_{s}$ , and $w$ are typically known from the microscope configuration, the defocus parameters $\phi$ vary widely between experiments. We must therefore estimate them to obtain an accurate model of the CTF.

To estimate $\phi$ , we turn to the power spectrum of the micrograph. The power spectra $S_{\mathbf{x}}$ , $S_{\mathbf{y}}$ , and $S_{\mathbf{e}}$ of $\mathbf{x}$ , $\mathbf{y}$ , and $\mathbf{e}$ , respectively, are related by

[TABLE]

This follows from (2) and the fact that convolving a stationary random field with $h_{\phi}$ multiplies its power spectrum by the square Fourier transform magnitude $|H_{\phi}|^{2}$ .

Equation (6) suggests that estimates of the power spectra $S_{\mathbf{y}}$ , $S_{\mathbf{x}}$ , and $S_{\mathbf{e}}$ can be useful in resolving the CTF. We note that $S_{\mathbf{x}}$ and $S_{\mathbf{e}}$ are slowly decaying while $|H_{\phi}|^{2}$ oscillates rapidly in comparison. The background subtracted power spectrum is therefore approximately proportional to $|H_{\phi}|^{2}$ . It follows that in order to estimate the defocus parameters $\phi$ , we may estimate $S_{\mathbf{y}}-S_{\mathbf{e}}$ and maximize its correlation with $|H_{\phi}|^{2}$ . This approach is used in [Rohou & Grigorieff, 2015, Mindell & Grigorieff, 2003, Zhang, 2016].

Another approach is to estimate $\phi$ from zero-crossings of $S_{\mathbf{y}}-S_{\mathbf{e}}$ [Yan et al., 2017, Tani et al., 1996]. Specifically, for spatial frequencies where $H_{\phi}(\mathbf{g})=0$ , we have $S_{\mathbf{y}}(\mathbf{g})-S_{\mathbf{e}}(\mathbf{g})=0$ . Identifying these zero-crossings from estimates of $S_{\mathbf{y}}-S_{\mathbf{e}}$ thus constrains the zeros of $H_{\phi}$ and lets us estimate its defocus parameters $\phi$ .

For both approaches, the first step is to estimate the background-subtracted power spectrum $S_{\mathbf{y}}-S_{\mathbf{e}}$ . In the following, we propose an estimation method and show how the resulting estimate may be used to estimate $\phi$ by either maximizing correlation or matching zero-crossings.

As mentioned above, the CTF is also multiplied by an exponentially decreasing envelope function [Sorzano et al., 2007], which effectively acts as a low-pass filter on $h_{\phi}*\mathbf{x}$ . In this paper, rather than include the envelope function in our analysis, we ignore high frequencies as they are strongly attenuated by the envelope. We also reduce the effect of the envelope function by estimating the CTF using the square root of our power spectrum estimate as in [Rohou & Grigorieff, 2015, Mindell & Grigorieff, 2003]. In this way, the effect of the envelope function on the two methods discussed is smaller.

2.2 Power spectrum estimation

In this section we present several methods for estimating the power spectrum of the micrograph. We first present the periodogram estimator and then show different methods for reducing its bias and variance.

2.2.1 Periodogram estimator

In an experimental setting, we only have access to an $N\times N$ sample of $\mathbf{y}$ , given by the values $\mathbf{y}[k_{1},k_{2}]$ for $(k_{1},k_{2})\in\{0,1,\ldots,N-1\}^{2}$ . Given these values, a common power spectrum estimator is provided by the periodogram [Oppenheim & Schafer, 1989]

[TABLE]

for $\mathbf{g}\in[-1/2,1/2]^{2}$ . While $\hat{S}^{(\mathrm{p})}_{\mathbf{y}}(\mathbf{g})$ may be calculated for any $\mathbf{g}$ , it is typically calculated on the $N\times N$ grid

[TABLE]

This enables the use of fast Fourier transforms (FFTs) for computing the periodogram with $O(N^{2}\log N)$ computational complexity. Due to this and its ease of implementation, the periodogram is a popular spectral estimator in cryo-EM.

Since our goal is to estimate $S_{\mathbf{y}}$ , let us consider how well it is estimated by the periodogram. The mean square error (MSE) of $\hat{S}^{(\mathrm{p})}_{\mathbf{y}}$ at $\mathbf{g}$ is given by

[TABLE]

To analyze the source of error, it is useful to define the bias and variance of the periodogram. The bias is defined as

[TABLE]

and measures the deviation of the expectation from the true value, while the variance

[TABLE]

measures the average deviation of the periodogram from its expectation. Both contribute to the MSE through the identity

[TABLE]

A low MSE therefore requires low bias and low variance.

The periodogram, however, fails on both counts. First, while the periodogram is asymptotically unbiased [Percival & Walden, 1993, Thomson, 1982], its bias remains large for small samples. Second, the periodogram is an inconsistent estimator, that is, its variance does not decrease with an increase in sample size. Therefore, a periodogram that extends over the entire micrograph will have variance approximately equal to that of a periodogram that extends over some section of the micrograph.

In the following sections, we will therefore consider different approaches to reducing both the bias and variance of the periodogram.

2.2.2 Bartlett’s method

We first consider an approach for reducing variance called Bartlett’s method [Oppenheim & Schafer, 1989]. In this approach, the periodogram estimate is computed for several non-overlapping regions of the image. These estimates are then averaged, reducing the variance by a factor approximately equal to the number of regions used. It may therefore be tempting to drastically reduce the size of these regions. However, in experimental data, averaging over regions that are too small will increase the bias. Among other things, this would prevent us from properly estimating the low spatial frequencies.

We thus divide our image into $B$ non-overlapping blocks $\mathbf{y}_{0},\ldots,\mathbf{y}_{B-1}$ of size $K\times K$ . The averaged periodogram is

[TABLE]

If each block $\mathbf{y}_{b}$ is independent of the others, we have $\operatorname{Var}(\hat{S}^{(\mathrm{b})}_{\mathbf{y}}(\mathbf{g}))=B^{-1}\operatorname{Var}(\hat{S}^{(\mathrm{p})}_{\mathbf{y}}(\mathbf{g}))$ . Note that, since the block size is now $K\times K$ , we sample $\mathbf{g}$ on $M_{K}$ .

2.2.3 Welch’s method

The expected value of the periodogram estimator is known to be a convolution between the true power spectrum of the micrograph and a 2D Fejér kernel [Percival & Walden, 1993]. As the Fejér kernel has high sidelobes, this convolution leads to frequency leakage and therefore a high bias.

One method of lowering the bias of the periodogram estimation is tapering [Percival & Walden, 1993]. This multiplies the data $\mathbf{y}$ by a data taper $\mathbf{w}$ prior to computing the periodogram, resulting in a modified periodogram. While many options for data tapers exist, such as the Hann window [Vulovic et al., 2012], Babadi and Brown [Babadi & Brown, 2014] suggest the use of the zeroth-order discrete prolate spheroidal sequence (DPSS) [Slepian, 1978]. The expected value of this modified periodogram is a convolution between the true power spectrum of the micrograph and a kernel with smaller sidelobes than those of the Fejér kernel [Percival & Walden, 1993]. This reduces the frequency leakage, and, therefore, the bias of the estimator.

While the taper may be applied to the entire micrograph, it is also possible to apply it to each block in Bartlett’s method (13). The resulting approach is known as Welch’s method [Welch, 1967]. Welch also showed that further variance reduction is possible using overlapping (typically half-overlapping) blocks [Fernández et al., 1997, Huang et al., 2003, Frank, 1996, Zhu et al., 1997]. This yields the modified averaged periodogram,

[TABLE]

where $\mathbf{y}_{b}\cdot\mathbf{w}$ is the pointwise multiplication of $\mathbf{y}_{b}$ and $\mathbf{w}$ .

2.2.4 Multitaper estimators

As discussed in Section 2.2.2, one way to lower the variance in the periodogram is to average several estimates. For this reason, Thomson [Thomson, 1982] suggested combining the estimates obtained from multiple tapers. Each taper yields a different estimate of the power spectrum, and averaging them significantly reduces the variance. A large number of tapers, however, results in significant smoothing of the power spectrum estimate, so the variance reduction needs to be balanced with an increase in bias for non-smooth power spectra. Thomson found that higher-order DPSSs were well-suited to this task and called the resulting power spectrum estimator the multitaper estimator. These estimators have recently demonstrated their usefulness for noise power spectrum estimation in cryo-EM [Andén & Singer, 2017, Andén & Romero, 2019]. For details of the DPSS data tapers we refer the reader to Appendix A.

Combining all the above methods for variance and bias reduction, we arrive at the multitaper estimator

[TABLE]

where $\mathbf{w}_{\ell}$ is the $\ell$ th out $L$ DPSSs for grids of size $K\times K$ .

Figs. 3-4 present a comparison between $\hat{S}^{(\mathrm{p})}_{\mathbf{y}}$ , $\hat{S}^{(\mathrm{b})}_{\mathbf{y}}$ , $\hat{S}^{(\mathrm{w})}_{\mathbf{y}}$ , and $\hat{S}^{(\mathrm{mt})}_{\mathbf{y}}$ . The CTF oscillations are best resolved by the multitaper estimator $\hat{S}^{(\mathrm{mt})}_{\mathbf{y}}$ .

2.3 Background subtraction

In this section we present a method for removing the background spectrum and further reducing the variability of the estimator $\hat{S}^{(\mathrm{mt})}_{\mathbf{y}}$ . We do this by first estimating the radial profile of the background spectrum and removing this estimation from $\hat{S}^{(\mathrm{mt})}_{\mathbf{y}}$ . Further, we show that the set of squared CTFs is contained in the span of a steerable basis and project our estimate onto that basis.

2.3.1 Estimating the background

We saw in (6) that the micrograph power spectrum can be expressed as the sum of two spectra: the clean, filtered power spectrum $|H_{\phi}(\mathbf{g})|^{2}\,S_{\mathbf{x}}(\mathbf{g})$ and the background $S_{\mathbf{e}}$ . The background-subtracted power spectrum is therefore

[TABLE]

An estimate of $S_{\mathbf{y}}-S_{\mathbf{e}}$ is used by many methods to estimate the CTF parameters $\phi$ [Rohou & Grigorieff, 2015, Zhu et al., 1997]. Their success therefore depends on accurate estimation of the background $S_{\mathbf{e}}$ .

The background is influenced by many factors and accurately modeling these factors is an open challenge. Many background estimation methods instead treat the background as a radially symmetric and slowly varying function [Frank, 1996].

The background estimation problem can be formulated as a curve fitting problem [Frank, 1996]. We note that the background should coincide with $\hat{S}^{(\mathrm{mt})}_{\mathbf{y}}$ at the zero-crossings of $H_{\phi}$ . Furthermore, the background should be strictly smaller than $\hat{S}^{(\mathrm{mt})}_{\mathbf{y}}$ at spatial frequencies where $H_{\phi}$ does not have a zero-crossing (since $S_{\mathbf{x}}$ is strictly positive). We therefore estimate the background by minimizing the difference between $\hat{S}^{(\mathrm{mt})}_{\mathbf{y}}$ and $S_{\mathbf{e}}$ .

While the radial profile of the background is monotonically decreasing in most settings, this is not the case when a Gatan K2 direct detector is used in counting mode with a high dose rate. Rather, the background will be monotonically decreasing in the lower frequencies and monotonically increasing in higher frequencies [Li et al., 2013]. Since a monotonically decreasing function, as well as a function that is at first monotonically decreasing and later monotonically increasing, must be convex, we model the background as the non-negative, convex function that is closest to, and no larger than, $\hat{S}^{(\mathrm{mt})}_{\mathbf{y}}$ .

We propose estimating the background $S_{\mathbf{e}}$ through LP. Specifically, we minimize the $\ell_{1}$ norm of the background-subtracted power spectrum estimate subject to several linear constraints. The first constraint ensures that the background-subtracted power spectrum estimate is non-negative, while the remaining constraints ensure that $\hat{S}_{\mathbf{e}}$ is a non-negative and convex.

Since we assume the background is radially symmetric, we consider its radial profile. To this end, we calculate the radial average of $\hat{S}^{(\mathrm{mt})}_{\mathbf{y}}(\mathbf{g})$ , which we denote, by a slight abuse of notation, $\hat{S}^{(\mathrm{mt})}_{\mathbf{y}}(r)$ . The radial averaging is performed by projecting $\hat{S}^{(\mathrm{mt})}_{\mathbf{y}}(\mathbf{g})$ on the circularly symmetric (i.e., purely radial) elements of a steerable basis (see Section 2.3.2).

The resulting linear program, whose result we denote by $\hat{S}_{\mathbf{e}}^{(\mathrm{lp})}$ , is then

[TABLE]

where $\hat{S}_{\mathbf{e}}=\begin{bmatrix}\hat{S}(0),\ldots,\hat{S}(m/K))\end{bmatrix}^{T}$ , and $0<m/K\leq 0.5$ is the spatial frequency above which $\hat{S}^{(\mathrm{mt})}_{\mathbf{y}}\left(r\right)$ is typically dominated by noise. As its default, ASPIRE-CTF sets $m/K=3/8$ .

We present the result of our linear program in Fig. 5.

Expanding the 1D background spectrum to a 2D function, we again abuse notation slightly and set $\hat{S}_{\mathbf{e}}^{(\mathrm{lp})}(\mathbf{g})=\hat{S}_{\mathbf{e}}^{(\mathrm{lp})}(r)$ for all $\mathbf{g}\in[-1/2,1/2]^{2}$ . We denote the background-subtracted power spectrum estimate $\hat{S}^{(\mathrm{mt})}_{\mathbf{y}}(\mathbf{g})-\hat{S}_{\mathbf{e}}^{(\mathrm{lp})}(\mathbf{g})$ by $\hat{S}^{(\mathrm{lp})}_{\mathbf{z}}(\mathbf{g})$ where $\mathbf{z}=h_{\phi}*\mathbf{x}$ .

A different LP that can be used for background estimation was suggested in [Huang et al., 2003]. However, contrary to our non-parametric approach which assumes convexity alone, [Huang et al., 2003] suggests a LP based on parametric estimation.

2.3.2 Expansion over a steerable basis

In this section, we show that any function of the form (3)-(5) is contained in a low-dimensional subspace spanned by a set of steerable basis functions, such as a Fourier–Bessel basis [Zhao & Singer, 2013, Zhao et al., 2016] or prolate spheroidal wave functions (PSWFs) [Landa & Shkolnisky, 2017, 2018]. We will use this property to further reduce the variability of the power spectrum estimator by projecting the background-subtracted power spectrum estimate $\hat{S}^{(\mathrm{lp})}_{\mathbf{z}}(\mathbf{g})$ onto this subspace.

A steerable basis consists of functions $f_{k,q}({r})\,\mathrm{e}^{\mathrm{j}k{\alpha}}$ , where $k\in\mathbb{Z}$ and $q=0,\ldots,p_{k}-1$ for some $p_{k}\geq 0$ . The radial part $f_{k,q}(r)$ depends on the specific choice of basis (e.g., in a Fourier–Bessel basis, it is a scaled Bessel function of order $q$ ) and does not enter explicitly into our analysis. We shall therefore leave it unspecified. A given function in polar coordinates may be decomposed in the basis as

[TABLE]

where $a_{k,q}\in\mathbb{C}$ is the coefficient corresponding to angular frequency $k$ and radial frequency $q$ .

To determine the steerable basis expansion of the CTF (3), we consider its Taylor expansion around $\Delta f_{1}-\Delta f_{2}=0$ ,

[TABLE]

where $R_{P}(\mathbf{g})$ is the remainder term,

[TABLE]

and

[TABLE]

is the non-astigmatic phase function.

The remainder term is bounded by a function of

[TABLE]

and is therefore small when astigmatism is small, which is the case for experimental cryo-EM data. We therefore conclude that

[TABLE]

is a good approximation of the CTF.

Since $\cos(\alpha)=\frac{1}{2}(\mathrm{e}^{-\mathrm{j}\alpha}+\mathrm{e}^{-\mathrm{j}\alpha})$ , we rewrite (20) as

[TABLE]

Comparing (21) and (18), we see that only terms corresponding to $k=-2$ , [math], and $2$ are present.

Concretely, we compute the coefficients $a_{k,q}$ of the expansion of $\sqrt{\hat{S}^{(\mathrm{lp})}_{\mathbf{z}}}$ over the steerable basis functions with radial frequencies to $k=0,\pm 2$ . The coefficients are computed through an inner product on a $K\times K$ grid:

[TABLE]

where $r$ and ${\alpha}$ are the polar coordinates of $\mathbf{g}$ . Evaluating (18) for these $a_{k,q}$ and squaring the result then gives a new power spectrum estimate, which we denote as $\hat{S}_{\mathbf{z}}$ .

In Fig. 6, we present the results of our power spectrum estimation method on a micrograph from the EMPIAR-10028 dataset Wong et al. [2014]. This includes multitaper estimate $\hat{S}^{(\mathrm{mt})}_{\mathbf{y}}$ as well as the background-subtracted estimate $\hat{S}^{(\mathrm{lp})}_{\mathbf{z}}$ and the projection onto the steerable basis with $k=0,\pm 2$ . The result is smooth enough that many of the zero-crossings of the power spectrum can be easily resolved.

2.4 CTF Parameter Estimation

In the previous sections we have introduced our method for estimating the background-subtracted power spectrum. In this section we will discuss two methods that use this estimate to recover the defocus and astigmatism of the micrograph.

2.4.1 CTF estimation through correlation

The power spectrum $S_{\mathbf{x}}$ is a slowly-varying function of radial frequency. It follows that the oscillations in $|H_{\phi}(\mathbf{g})|^{2}\,S_{\mathbf{x}}(\mathbf{g})$ are due to those of $|H_{\phi}(\mathbf{g})|^{2}$ . As a consequence, the square root of the background-subtracted power spectrum $|H_{\phi}(\mathbf{g})|\,S_{\mathbf{x}}^{1/2}(\mathbf{g})$ is proportional to the absolute value of the CTF (we note that using the square root instead of the actual power spectrum estimate reduces the influence of large values). It follows that the correlation test is a useful tool in estimating CTF parameters. Indeed, many CTF estimation methods therefore estimate the defocus parameters $\phi$ through correlation with a simulated CTF magnitude [Rohou & Grigorieff, 2015, Mindell & Grigorieff, 2003, Zhang, 2016].

The Pearson correlation of $|H_{\phi}(\mathbf{g})|$ and $\hat{S}_{\mathbf{z}}^{1/2}(\mathbf{g})$ is

[TABLE]

where $R$ is the set of frequencies over which correlation is computed, and will be defined below. To optimize $P_{cc}(\phi)$ , we first need an initial guess for the parameters $\phi$ . For this, we follow Zhang [2016] and first consider non-astigmatic CTFs where $\Delta f_{1}=\Delta f_{2}$ , which renders the value of $\alpha_{f}$ irrelevant. We thus calculate $P_{cc}(\phi)$ for $\phi=(\Delta f,\Delta f,0)$ with $\Delta f$ on a 1D grid from $\Delta f_{\mathrm{min}}$ to $\Delta f_{\mathrm{max}}$ with a step of $\Delta f_{\mathrm{step}}$ . Since $H_{\phi}$ is considered (at this stage) to be radially symmetric, we define the set of frequencies over which correlation is computed as

[TABLE]

where $m_{1}$ is the first maximum of the radial profile of $\hat{S}^{(\mathrm{lp})}_{\mathbf{z}}$ . That is, we only consider frequencies $\mathbf{g}$ along a 1D radial profile and, furthermore, ignore the very low and very high frequencies (since the very low frequencies may dominate the cross-correlation result and the very high frequencies are strongly effected by the envelope function). The $\Delta f$ which maximizes $P_{cc}(\phi)$ on this grid is denoted $\Delta f_{\star}$ .

To estimate the astigmatism of the CTF, we compute the principal directions of the second-order moments of $\hat{S}_{\mathbf{z}}^{1/2}$ . Specifically, we form the $2\times 2$ matrix $\mathbf{M}$ given by

[TABLE]

The eigenvalues $\mu_{1}$ and $\mu_{2}$ of $\mathbf{M}$ estimate the size of the major and minor axes in $\hat{S}_{\mathbf{z}}^{1/2}$ . We therefore expect the ratio $\mu_{1}/\mu_{2}$ to approximate $\Delta f_{1}/\Delta f_{2}$ . Combining this with the estimated mean defocus $\Delta f_{\star}$ , we get

[TABLE]

which has the solution

[TABLE]

In order to improve our estimation of the defocus parameters, we run gradient descent on $P_{cc}(\phi)$ . As we no longer approximate the image as non-astigmatic, we define the set of frequencies over which correlation is computed as $R=M_{K}$ . We now initialize our gradient descent at $\phi=(\Delta f_{1,\star},\Delta f_{2,\star},\alpha_{f})$ , where $\alpha_{f}$ is set as detailed in [Zhang, 2016] to an arbitrarily selected $0\leq a<\pi/6$ (e.g. $a=\pi/12$ ) and $(a+\pi/6)$ , $(a-\pi/6)$ , $(a+\pi/3)$ , $(a-\pi/3)$ or $(a-\pi/2)$ . One run of gradient descent is performed for each value of $\alpha_{f}$ and the result with the highest value of $P_{cc}(\phi)$ is kept. The resulting $\phi$ is our defocus estimate for the micrograph.

We note that, as is done in [Rohou & Grigorieff, 2015, Mindell & Grigorieff, 2003, Zhang, 2016], we discard information in the lower and higher frequencies of $\hat{S}_{\mathbf{z}}$ . These frequencies can be determined by the user. As default values we use those suggested in [Rohou & Grigorieff, 2015].

2.4.2 CTF estimation through zero-crossings

We have seen in the previous sections that the true background-subtracted power spectrum is $|H_{\phi}(\mathbf{g})|^{2}S_{\mathbf{x}}(\mathbf{g})$ . Under the assumption that $S_{\mathbf{x}}(\mathbf{g})$ is slowly-varying, it follows that at any frequency where the background-subtracted power spectrum reaches a minimal value of zero, the CTF must reach a zero-crossing.

Furthermore, even if the aforementioned assumption did not hold true, we could still infer the frequencies where the CTF reaches a zero-crossing. This is due to the fact that the zero-crossings of the CTF are known to create concentric, nearly elliptical rings, centered around the origin (see Section 2.1). Therefore, this can be used a cue to differentiate between any minima of $|H_{\phi}(\mathbf{g})|^{2}\,S_{\mathbf{x}}(\mathbf{g})$ that stem from the zero-crossings of $H_{\phi}$ and any minima that stem from $S_{\mathbf{x}}$ .

As we show in Fig. 6(c), our estimation of the background-subtracted power spectrum enables easy detection of several elliptical rings where $\hat{S}_{\mathbf{z}}$ reaches its minima. To do this, we define any pixel with a value smaller than that of at least six out of its eight neighbors as a zero-crossing. Once the minima of $\hat{S}_{\mathbf{z}}$ are found, we discard any frequency that is not on a closed ring. Furthermore, we verify that the spatial frequencies of pixels residing on closed rings representing the minima of $\hat{S}_{\mathbf{z}}$ form ellipses approximately centered at the origin. We are then left with frequencies of several zero-crossings of the CTF.

Since we have $H_{\phi}(\mathbf{g})=-\sin(\chi_{\phi}(\mathbf{g}))$ , we reach a zero-crossing of the CTF when $\chi_{\phi}(\mathbf{g})$ is an integer multiple of $\pi$ . Formally, the set of spatial frequencies on the $\ell$ th ring of zero-crossings, denoted by $G_{\ell}$ , satisfies

[TABLE]

Empirically, we are typically able to identify at least three rings, that is, three different values of $\ell$ .

Combining (29) for all $\mathbf{g}$ in $\hat{G}_{\ell}$ and combining these for different values of $\ell$ , we obtain an overdetermined system of equations. Solving it yields an estimate for the defocus vector $\phi$ . To solve the system, we use the trust-region-dogleg method implemented in MATLAB (a variant of [Powell, 1970]).

We note that the estimated positions of CTF zeroes are extremely sensitive to the method of background subtraction. This method should therefore be used to obtain an initial estimate of the defocus. Refinement of the solution can be done as in Section 2.4.1 using gradient-based optimization.

We note that typically both methods suggested in this section achieve similar results. However, while the zero-crossings-based method has lower computational complexity, the correlation-based method is more robust to noise. Therefore, for micrographs with very low SNR we recommend using the correlation-based method, while for cleaner micrographs we suggest using the zero-crossings method.

3 Results

We present experimental results for the ASPIRE-CTF framework presented in this paper. We apply our framework to datasets that are publicly available from the EMPIAR database [Iudin et al., 2016] or the CTF challenge [Marabini et al., 2015]. Unless otherwise stated, in the experiments below we use $L=4$ tapers and project the power spectrum onto the steerable PSWF basis.

3.1 Estimating CTF from movie frames

The CTF may be estimated either from motion-corrected micrographs, from several frames averaged in real space or, alternatively, directly from the frames. This is done by estimating $\hat{S}^{(\mathrm{mt})}$ individually from each frame and averaging the estimates (see Fig. 1). A benefit of estimating the CTF directly from the frames is that this practice enables us to correct for motion and estimate the CTF concurrently, thus speeding up the pipeline. Furthermore, any errors added by motion-correction will have no effect on the CTF estimation [Bartesaghi et al., 2014].

In this section, we compare the CTF estimates produced from motion-corrected micrographs to the estimation produced directly from the frames. We do this over several publicly available datasets, namely, EMPIAR-10002 [Bai et al., 2013], EMPIAR-10028 [Wong et al., 2014], EMPIAR-10242 [Zhang et al., 2019], and EMPIAR-10249 [Herzik et al., 2019]. A summary of these datasets appears in Table 1.

While the EMPIAR-10028 and EMPIAR-10242 datasets contain both movies and motion-corrected micrographs, EMPIAR-10002 and EMPIAR-10249 contain movies alone. We therefore use MotionCor2 [Zheng et al., 2017] to produce the motion-corrected micrograph for these two datasets. We present in Fig. 7 a comparison of the astigmatism ( $\Delta f_{1}-\Delta f_{2}$ ), average defocus ( $\Delta f_{1}/2+\Delta f_{2}/2$ ), and astigmatism angle ( $\alpha_{f}$ ) as estimated from a motion-corrected micrograph with an estimate produced from the raw movie frames. We note that, as expected, the parameters estimated from each of these methods are nearly identical.

3.2 CTF Challenge

The CTF challenge [Marabini et al., 2015] consists of nearly $200$ micrographs of GroEL, 60S ribosome, apoferritin and TMV virus. These micrographs are taken from eight experimental datasets and one synthetic dataset, each referred to by a number ranging from $001$ to $009$ . In the following, we restrict our attention to the experimental datasets, that is, sets $001$ through $008$ .

The advantage of the CTF challenge is that each dataset is acquired using a different combination of microscope and camera, allowing for a qualitative comparison of CTF estimation methods for a variety of experimental setups. Notably, datasets $003$ and $004$ use a Gatan K2 direct detector in counting mode with a high electron dose, causing $S_{\mathbf{e}}$ to increase at high frequencies [Li et al., 2013]. Additionally, dataset $008$ has an especially low signal-to-noise ratio (SNR), rendering CTF estimation difficult. A summary of these datasets is presented in [Marabini et al., 2015].

The estimate of each micrograph’s power spectrum is computed as detailed in Section 2.2.4. Specifically, we divide the micrograph into half-overlapping blocks of size $K\times K$ , where $K=512$ , and use $L=4$ tapers in the estimation. We then estimate the background spectrum as detailed in Section 2.3.1 and expand the background-subtracted power spectrum over the PSWF basis in order to reduce variability in the power spectrum estimate (Section 2.3.2). We use the correlation-based method (Section 2.4.1) to estimate defocus parameters, and denote the resulting vector of defocus parameters as $\phi_{a}^{(512)}$ .

In some cases, a power spectrum of size $512\times 512$ may not capture the oscillations of the power spectrum with sufficient accuracy [Rohou & Grigorieff, 2015]. We therefore compute a second estimate of the power spectrum using half-overlapping blocks of size $1024\times 1024$ . As this reduces the number of blocks, the variance of the estimator will grow. We therefore use $L=16$ data tapers in this case. We use this estimate of the power spectrum to estimate a vector of defocus parameters which we denote by $\phi_{a}^{(1024)}$ .

We compare our results to the estimates produced by CTFFIND4 (version 4.1.13) and Gctf (version 1.06). We denote the vector of estimated defocus parameters produced by CTFFIND4 when using block of size $512\times 512$ and $1024\times 1024$ as $\phi_{c}^{(512)}$ and $\phi_{c}^{(1024)}$ , respectively. We further denote the vector of estimated defocus parameters produced by Gctf when using block of size $512\times 512$ and $1024\times 1024$ as $\phi_{g}^{(512)}$ and $\phi_{g}^{(1024)}$ , respectively. For each estimation method we select the vector of estimated defocus parameters that leads to highest correlation with the estimated power spectrum, that is

[TABLE]

where $j\in\{a,c,g\}$ , $s_{t}$ is the number of micrographs in the $t$ th dataset and $P_{cc}^{m}$ is the correlation for the $m$ th micrograph in the dataset, computed as in (23). We note that the correlation is computed with the power spectrum estimate suggested in [Mindell & Grigorieff, 2003]. That is, we compute the Pearson correlation coefficient between the background subtracted power spectrum computed as in [Mindell & Grigorieff, 2003] and using blocks of size $512\times 512$ with $H_{\phi_{j}^{(512)}}$ , and between the background subtracted power spectrum computed using blocks of size $1024\times 1024$ with $H_{\phi_{j}^{(1024)}}$ . In this manner, we choose the block size that best captures the oscillations of each dataset.

In order to compare the consistency of these $3$ methods, we present the differences between $\phi_{a}$ , $\phi_{c}$ and $\phi_{g}$ in Tables 2-3. That is, for each micrograph $m$ we compute

[TABLE]

where $j$ and $k$ are two estimation methods (ASPIRE-CTF, CTFFIND4 or Gctf). We report the mean and variance of $\epsilon$ . We note that either $\epsilon_{a,g}$ or $\epsilon_{a,c}$ are often smaller than $\epsilon_{c,g}$ , thus showing the ASPIRE-CTF estimate to be in the consensus of the three estimation vectors.

Fig. 8 contains a visual comparison between the power spectrum computed by our suggested framework and the power spectra computed by Gctf and CTFFIND4. We present the comparison over a micrograph from the eighth set of the CTF challenge as this set is known to be difficult. We note that the oscillations of the ASPIRE-CTF power spectrum are highly noticeable. In comparison, the variability of the power spectra computed by Gctf and CTFFIND4 make visual detection of oscillations challenging.

3.3 Runtime

We compute runtime of ASPIRE-CTF and CTFFIND4 over dataset $001$ of the CTF challenge. For both methods, we partition the micrograph into blocks of size $512\times 512$ . When running ASPIRE-CTF we employ $L=4$ data tapers. Furthermore, we use the exhaustive search option for CTFFIND4, and perform an exhaustive 1D search in ASPIRE-CTF. While the CTF estimation results are comparable, there is a significant speedup when using ASPIRE-CTF. Runtime for ASPIRE-CTF is 22.5 seconds on average per micrograph, while runtime for CTFFIND4 is 541 seconds.

These experiments are run on a $2.6$ GHz Intel Core i7 CPU with four cores and $16$ GB of memory. We do not compare to the runtime of Gctf as it must be run on a GPU.

3.4 Consistency in low SNR

To test consistency of results with changing SNR, we turn to the EMPIAR-10249 dataset Herzik et al. [2019]. This dataset consists of movies with $44$ frames per movie. Usually, all these frames, except for a few frames at the beginning and a few at the end, are motion-corrected and summed to create a micrograph. This is due to the fact that a micrograph created from as many motion-corrected frames as possible will have the best SNR.

We disregard the first frame and use MotionCor2 [Zheng et al., 2017] to create $9$ motion-corrected micrographs. These consist of summing $5$ , $8$ , $13$ , $18$ , $23$ , $28$ , $33$ , $38$ , and $43$ motion-corrected frames, respectively. This gives us a sequence of micrographs with increasing SNR.

We estimated the CTF parameters independently from each micrograph in the manner detailed in Section 3.2. Fig. 9 shows the astigmatism $|\Delta f_{1}-\Delta f_{2}|$ vs. mean defocus $(\Delta f_{1}+\Delta f_{2})/2$ of the CTF estimation for each method and over each micrograph. We see that while the average defocus values remain similar for all three methods, Gctf incurs a larger error in the astigmatism when 23 frames are used. On the other hand, our method and CTFFIND4 achieve consistent estimates regardless of the amount of frames averaged.

4 Conclusion

In this paper we have presented a novel approach for power spectrum estimation of cryo-EM experimental data. Our approach uses the multitaper estimator, which often leads to reduced mean square error over Bartlett’s and Welch’s methods. Additionally, we presented a method for error reduction that is driven directly by the mathematical model of the contrast transfer function. We did this by projecting the power spectrum estimate onto a steerable basis and discarding any basis function where the CTF must be negligible. We showed that the combination of these two contributions leads to greatly reduced variability in our estimator.

We presented experimental results on twelve datasets, and showed that our method is well suited to both motion-corrected micrographs and raw movies data.

Acknowledgments

This work was partially supported by the Simons Foundation Math+X Investigator Award and the Moore Foundation Data-Driven Discovery Investigator Award. The authors thank B. Landa and I. Sason for help optimizing the PSWF code. The authors are also indebted to B. Landa, Y. Shkolnisky and A. Rohou for helpful comments and discussions. The Flatiron Institute is a division of the Simons Foundation.

Appendix A

Zeroth-order discrete prolate spheroidal sequence (DPSS) [Slepian, 1978] were proposed as data tapers in (Babadi & Brown, 2014; Section 2.2.3). Here we describe their generation.

The zeroth-order discrete prolate spheroidal sequence is a sequence of $d$ 1-D vectors, determined as the $d$ leading eigenvectors of the matrix $\mathbf{\mathcal{L}}\in\mathbb{R}^{K\times K}$ , where

[TABLE]

and $R=\frac{2d}{N}$ . We denote the resulting data tapers by $\mathbf{t}_{1},\dots,\mathbf{t}_{N}$ .

As the blocks $\mathbf{y}_{b}$ are 2D, that is $\mathbf{y}_{b}\in\mathbb{R}^{K\times K}$ , the data tapers we use are defined as

[TABLE]

Lastly, we note that $d$ is selected such that

[TABLE]

References

Andén & Romero [2019]

Andén, J., & Romero, J. L. (2019).

Multitaper estimation on arbitrary domains.

Submitted to SIAM J. Imag. Sci., arXiv:1812.03225.

Andén & Singer [2017]

Andén, J., & Singer, A. (2017).

Factor analysis for spectral estimation.

In Proc. SampTA (pp. 169–173).

doi:10.1109/SAMPTA.2017.8024447.

Babadi & Brown [2014]

Babadi, B., & Brown, E. N. (2014).

A review of multitaper spectral analysis.

IEEE Trans. Biomed. Eng., 61, 1555–1564.

Bai et al. [2013]

Bai, X.-C., Fernandez, I. S., McMullan, G., & Scheres, S. H. W. (2013).

Ribosome structures to near-atomic resolution from thirty thousand cryo-EM particles.

eLife, 2, e00461.

Bartesaghi et al. [2014]

Bartesaghi, A., Matthies, D., Banerjee, S., Merk, A., & Subramaniam, S. (2014).

Structure of $\beta$ -galactosidase at 3.2-Å resolution obtained by cryo-electron microscopy.

Proceedings of the National Academy of Sciences, 111, 11709–11714. doi:10.1073/pnas.1402809111.

Bhamre et al. [2016]

Bhamre, T., Zhang, T., & Singer, A. (2016).

Denoising and covariance estimation of single particle cryo-EM images.

J. Struct. Biol., 195, 72–81.

Cheng et al. [2017]

Cheng, Y., Glaeser, R. M., & Nogales, E. (2017).

How cryo-EM became so hot.

Cell, 171, 1229–1231.

Erickson & Klug [1971]

Erickson, H. P., & Klug, A. (1971).

Measurement and compensation of defocusing and aberrations by Fourier processing of electron micrographs.

Phil. Trans. R. Soc. Lond. B, 261, 105–118.

Fernández et al. [1997]

Fernández, J. J., Sanjurjo, J. R., & Carazo, J.-M. (1997).

A spectral estimation approach to contrast transfer function detection in electron microscopy.

Ultramicroscopy, 68, 267–295.

Frank [1996]

Frank, J. (1996).

Three-Dimensional Electron Microscopy of Macromolecular Assemblies.

Academic Press.

Grant et al. [2018]

Grant, T., Rohou, A., & Grigorieff, N. (2018).

cisTEM, user-friendly software for single-particle image processing.

eLife, 7, e35383.

Heimowitz et al. [2018]

Heimowitz, A., Andén, J., & Singer, A. (2018).

APPLE picker: Automatic particle picking, a low-effort cryo-EM framework.

J. Struct. Biol., 204, 215–227.

Herzik et al. [2019]

Herzik, M. A., Wu, M., & Lander, G. C. (2019).

High-resolution structure determination of sub-100 kDa complexes using conventional cryo-EM.

Nat. Commun., 10, 1032. doi:10.1038/s41467-019-08991-8.

Huang et al. [2003]

Huang, Z., Baldwin, P. R., Mullapudi, S., & Penczek, P. A. (2003).

Automated determination of parameters describing power spectra of micrograph images in electron microscopy.

J. Struct. Biol., 144, 79–94.

Iudin et al. [2016]

Iudin, A., Korir, A. K., Salavert-Torres, J., Kleywegt, G. J., & Patwardhan, A. (2016).

EMPIAR: A public archive for raw electron microscopy image data.

Nat. Methods, 13.

Landa & Shkolnisky [2017]

Landa, B., & Shkolnisky, Y. (2017).

Approximation scheme for essentially bandlimited and space-concentrated functions on a disk.

Appl. Comput. Harmon. Anal., 43, 381–403.

Landa & Shkolnisky [2018]

Landa, B., & Shkolnisky, Y. (2018).

Steerable principal components for space-frequency localized images.

SIAM J. Imaging Sci., 10, 508–534.

Li et al. [2013]

Li, X., Zheng, S. Q., Egami, K., Agard, D. A., & Cheng, Y. (2013).

Influence of electron dose rate on electron counting images recorded with the K2 camera.

J. Struct. Biol., 184, 251–260.

Marabini et al. [2015]

Marabini, R., Carragher, B., Cheni, S., Chen, J., Cheng, A., Downing, K. H. et al. (2015).

CTF Challenge: Result summary.

J. Struct. Biol., 190, 348–359.

Mindell & Grigorieff [2003]

Mindell, J. A., & Grigorieff, N. (2003).

Accurate determination of local defocus and specimen tilt in electron microscopy.

J. Struct. Biol., 142, 334–347.

Oppenheim & Schafer [1989]

Oppenheim, A. V., & Schafer, R. W. (1989).

Discrete-Time Signal Processing.

(1st ed.).

Prentice Hall.

Percival & Walden [1993]

Percival, D. B., & Walden, A. T. (1993).

Spectral Analysis for Physical Applications.

Cambridge University Press.

Powell [1970]

Powell, M. J. D. (1970).

Numerical Methods for Nonlinear Algebraic Equations.

Punjani et al. [2017]

Punjani, A., Rubinstein, J. L., Fleet, D. J., & Brubaker, M. A. (2017).

cryoSPARC: algorithms for rapid unsupervised cryo-EM structure determination.

Nat. Methods, 14, 290.

Rohou & Grigorieff [2015]

Rohou, A., & Grigorieff, N. (2015).

CTFFIND4: Fast and accurate defocus estimation from electron micrographs.

J. Struct. Biol., 192, 216–221.

Scheres [2015]

Scheres, S. H. (2015).

Semi-automated selection of cryo-EM particles in RELION-1.3.

J. Struct. Biol., 189, 114–122. doi:https://doi.org/10.1016/j.jsb.2014.11.010.

Scheres [2012]

Scheres, S. H. W. (2012).

RELION: Implementation of a Bayesian approach to cryo-EM structure determination.

J. Struct. Biol., 180, 519–530.

Slepian [1978]

Slepian, D. (1978).

Prolate spheroidal wave functions, Fourier analysis, and uncertainty—V: The discrete case.

Bell Syst. Tech. J., 57, 1371–1430.

Sorzano et al. [2007]

Sorzano, C., Jonic, S., Núñez-Ramírez, R., Boisset, N., & Carazo, J. M. (2007).

Fast, robust, and accurate determination of transmission electron microscopy contrast transfer function.

J. Struct. Biol., 160, 249–262.

Tang et al. [2007]

Tang, G., Peng, L., Baldwin, P. R., Mann, D. S., Jiang, W., Rees, I. et al. (2007).

EMAN2: An extensible image processing suite for electron microscopy.

J. Struct. Biol., 157, 38–46.

Tani et al. [1996]

Tani, K., Sasab, H., & Toyoshima, C. (1996).

A set of computer programs for determining defocus and astigmatism in electron images.

Ultramicroscopy, 65, 31–44.

Thomson [1982]

Thomson, D. J. (1982).

Spectrum estimation and harmonic analysis.

Proc. IEEE, 70, 1055–1096.

Thon [1971]

Thon, F. (1971).

Phase contrast electron microscopy.

In U. Valdré (Ed.), Electron Microscopy in Material Science.

Academic Press.

Vulovic et al. [2012]

Vulovic, M., Franken, E. M., Ravelli, R. B. G., van Vliet, L. J., & Rieger, B. (2012).

Precise and unbiased estimation of astigmatism and defocus in transmission electron microscopy.

Ultramicroscopy, 116, 115–134.

Welch [1967]

Welch, P. D. (1967).

The use of fast Fourier transform for the estimation of power spectra: A method based on time averaging over short, modified periodograms.

IEEE Trans. Audio Electroacoust., 15, 70–73.

Wong et al. [2014]

Wong, W., Bai, X.-C., Brown, A., Fernandez, I. S., Hanssen, E., Condron, M. et al. (2014).

Cryo-EM structure of the Plasmodium falciparum 80S ribosome bound to the anti-protozoan drug emetine.

eLife, 3.

Yan et al. [2017]

Yan, R., Li, K., & Jiang, W. (2017).

Real-time detection and single-pass minimization of TEM objective lens astigmatism.

J. Struct. Biol., 197, 210–219.

Zhang [2016]

Zhang, K. (2016).

Gctf: Real-time CTF determination and correction.

J. Struct. Biol., 193, 1–12.

Zhang et al. [2019]

Zhang, W., Falcon, B., Murzin, A. G., Fan, J., Crowther, R. A., Goedert, M. et al. (2019).

Heparin-induced tau filaments are polymorphic and differ from those in Alzheimer’s and Pick’s diseases.

eLife, 8. doi:10.7554/eLife.43584.

Zhao et al. [2016]

Zhao, Z., Shkolnisky, Y., & Singer, A. (2016).

Fast steerable principal component analysis.

IEEE Trans. Comput. Imaging, 2, 1–12.

Zhao & Singer [2013]

Zhao, Z., & Singer, A. (2013).

Fourier–Bessel rotational invariant eigenimages.

J. Opt. Soc. Am. A, 30, 871–877. doi:10.1364/JOSAA.30.000871.

Zheng et al. [2017]

Zheng, S. Q., Palovcak, E., Armache, J.-P., Verba, K. A., Cheng, Y., & Agard, D. A. (2017).

MotionCor2: anisotropic correction of beam-induced motion for improved cryo-electron microscopy.

Nat. Methods, 14, 331–331.

Zhu et al. [1997]

Zhu, J., Penczek, P. A., Schröder, R., & Frank, J. (1997).

Three-dimensional reconstruction with contrast transfer function correction from energy-filtered cryoelectron micrographs: Procedure and application to the 70S Escherichia coli ribosome.

J. Struct. Biol., 118, 197–219.

Bibliography43

The reference list from the paper itself. Each links out to its DOI / PubMed record.

1Andén & Romero [2019] Andén, J., & Romero, J. L. (2019). Multitaper estimation on arbitrary domains. Submitted to SIAM J. Imag. Sci., ar Xiv:1812.03225 .
2Andén & Singer [2017] Andén, J., & Singer, A. (2017). Factor analysis for spectral estimation. In Proc. Samp TA (pp. 169–173). doi: 10.1109/SAMPTA.2017.8024447 . · doi ↗
3Babadi & Brown [2014] Babadi, B., & Brown, E. N. (2014). A review of multitaper spectral analysis. IEEE Trans. Biomed. Eng. , 61 , 1555–1564.
4Bai et al. [2013] Bai, X.-C., Fernandez, I. S., Mc Mullan, G., & Scheres, S. H. W. (2013). Ribosome structures to near-atomic resolution from thirty thousand cryo-EM particles. e Life , 2 , e 00461.
5Bartesaghi et al. [2014] Bartesaghi, A., Matthies, D., Banerjee, S., Merk, A., & Subramaniam, S. (2014). Structure of β 𝛽 \beta -galactosidase at 3.2-Å resolution obtained by cryo-electron microscopy. Proceedings of the National Academy of Sciences , 111 , 11709–11714. doi: 10.1073/pnas.1402809111 . · doi ↗
6Bhamre et al. [2016] Bhamre, T., Zhang, T., & Singer, A. (2016). Denoising and covariance estimation of single particle cryo-EM images. J. Struct. Biol. , 195 , 72–81.
7Cheng et al. [2017] Cheng, Y., Glaeser, R. M., & Nogales, E. (2017). How cryo-EM became so hot. Cell , 171 , 1229–1231.
8Erickson & Klug [1971] Erickson, H. P., & Klug, A. (1971). Measurement and compensation of defocusing and aberrations by Fourier processing of electron micrographs. Phil. Trans. R. Soc. Lond. B , 261 , 105–118.