Practical: Single-Particle Analysis

13. Practical: Single-Particle Analysis#

13.1. The reconstruction problem#

In electron microscopy we want to deduce the 3D structure of a structurally homogeneous ‘particle’ — for example a protein — from cryo-EM images. Each image shows the particle in a different, unknown orientation. To make the problem tractable we work in 2D: we treat the unknown structure as a 2D image \(A\) and each observation as a rotated, noisy version:

(13.1)#\[ X_i = R^{\theta_i} A + \sigma G_i \]

where \(R^{\theta_i}\) is rotation by the unknown angle \(\theta_i\) and \(G_i\) is i.i.d. Gaussian noise. If we can estimate \(\theta_i\) for each image, we can recover \(A\) by rotating each image back and averaging.

We use the letter Q as our 2D ‘particle’ because it is immediately recognisable and its asymmetry makes orientation unambiguous. In real SPA the structure is of course unknown.

Interactive elements

Click Live Code in the top toolbar to activate the kernel, then expand each Show code toggle and click ▶ to run the cell.

13.1.1. Helper functions#

The cell below defines all functions used throughout this chapter. Run it first.

Show code cell source

Hide code cell source

import io, warnings
import numpy as np
import matplotlib
matplotlib.use('agg')
import matplotlib.pyplot as plt
from scipy import ndimage
from scipy.ndimage import zoom as nd_zoom, rotate as nd_rotate, gaussian_filter
from ipywidgets import (IntSlider, FloatSlider, Dropdown, VBox, HBox,
                        Layout, Output)
from IPython.display import display, Image, clear_output
warnings.filterwarnings('ignore')

# ── image helpers ─────────────────────────────────────────────────────────────

def make_letter(text, size=64):
    """Render a letter as a 2D float array via matplotlib Agg canvas."""
    fig, ax = plt.subplots(figsize=(2, 2))
    ax.axis('off')
    fig.text(0.23, 0.26, text, fontsize=100)
    fig.canvas.draw()
    buf = np.frombuffer(fig.canvas.buffer_rgba(), dtype=np.uint8)
    w, h = fig.canvas.get_width_height()
    data = buf.reshape(h, w, 4)[:, :, :3].mean(2)
    plt.close(fig)
    data = (data - data.min()) / (data.max() - data.min() + 1e-10)
    data = (~data.astype(bool))[::-1].astype(float)
    return nd_zoom(data, size / data.shape[0])

def make_circle(size=64, radius_frac=0.35):
    """Filled circle — a featureless starting reference."""
    y, x = np.mgrid[0:size, 0:size].astype(float)
    return ((x - size/2)**2 + (y - size/2)**2 < (size*radius_frac)**2).astype(float)

def low_pass_filter(im, sigma=3.0):
    """Gaussian low-pass filter (approximates tanh LP filter)."""
    return gaussian_filter(im, sigma=sigma)

def show(im, ax=None, title='', vmin=None, vmax=None, cmap='gray'):
    if ax is None:
        _, ax = plt.subplots()
    ax.imshow(im, cmap=cmap, origin='lower',
              vmin=vmin or im.min(), vmax=vmax or im.max())
    ax.axis('off')
    if title:
        ax.set_title(title, fontsize=9)

def fig2img(fig):
    buf = io.BytesIO(); fig.savefig(buf, format='png', dpi=100, bbox_inches='tight')
    buf.seek(0); plt.close(fig)
    return Image(data=buf.read())

# ── simulation ────────────────────────────────────────────────────────────────

def simulate_image(image, angle, sigma=1, rng=None):
    """Rotate image and add Gaussian noise."""
    rng = rng or np.random.default_rng()
    return nd_rotate(image, angle, reshape=False) + sigma * rng.standard_normal(image.shape)

def simulate_images(image, sigma, N, seed=42):
    """Generate N noisy randomly-rotated images."""
    rng = np.random.default_rng(seed)
    angles = rng.uniform(0, 360, N)
    imgs = np.array([simulate_image(image, a, sigma, rng) for a in angles])
    return angles, imgs

# ── alignment ─────────────────────────────────────────────────────────────────

def rotate_reference(reference, delta_angle=5):
    """Rotate reference through all candidate angles in [0, 360)."""
    angles = np.arange(0, 360, delta_angle)
    stack  = np.array([nd_rotate(reference, a, reshape=False) for a in angles])
    return angles, stack

def correlation(image, reference_stack):
    """Cross-correlation score for each reference in the stack."""
    return (image[None] * reference_stack).mean(axis=(1, 2))

def align_images(images_stack, reference_stack, cand_angles):
    """Return estimated angle for each image."""
    estimated = []
    for img in images_stack:
        scores   = correlation(img, reference_stack)
        best_idx = int(scores.argmax())
        estimated.append(cand_angles[best_idx])
    return np.array(estimated)

# ── reconstruction ─────────────────────────────────────────────────────────────

def reconstruct(images_stack, estimated_angles):
    """Average back-rotated images."""
    recon = np.zeros_like(images_stack[0])
    for angle, img in zip(estimated_angles, images_stack):
        recon += nd_rotate(img, -angle, reshape=False)
    return recon / len(images_stack)

print("Helper functions loaded.")

13.2. Part 0 – The test image#

13.3. Simulating cryo-EM images#

The interactive below shows a single simulated particle image. The image formation model is:

\[ X = R^\theta A + \sigma G \]

Drag the angle slider to rotate the particle and the σ slider to add noise.

13.4. Part 1 – Choosing an initial reference#

The alignment algorithm needs a starting reference to compare images against. Three strategies are common:

Reference	Advantage	Risk
Low-pass filtered true structure	Fast convergence	Only works if structure is already known (circular reasoning)
Featureless circle	No reference bias	Slow convergence; poor angular discrimination
Random noise / another structure	Quick to obtain	Risk of converging to wrong answer (model bias)

In practice, a featureless circle or a low-pass-filtered version of a prior reconstruction is used. The interactive below lets you choose the starting reference and see how much angular discrimination it provides for a single noisy image.

Show code cell source

Hide code cell source

_sl2 = Layout(width="420px")
ref_dd  = Dropdown(options=['Q (oracle)', 'Blurred Q (σ=5)', 'Circle',
                             'P (wrong)', 'O (wrong)'],
                   value='Blurred Q (σ=5)', description="Reference:", style=_sty)
sigma_sl1 = FloatSlider(value=2.0, min=0.0, max=8.0, step=0.5,
                         description="Noise σ", style=_sty, layout=_sl2)

_REFS = {
    'Q (oracle)'       : Q,
    'Blurred Q (σ=5)'  : low_pass_filter(Q, sigma=5),
    'Circle'           : make_circle(64),
    'P (wrong)'        : make_letter('P', 64),
    'O (wrong)'        : make_letter('O', 64),
}
# Pre-compute candidate rotations at 5° for each reference
_CAND5 = np.arange(0, 360, 5)
_STACKS5 = {name: np.array([nd_rotate(ref, a, reshape=False) for a in _CAND5])
            for name, ref in _REFS.items()}

out_ref = Output()
def _update_ref(_=None):
    rng = np.random.default_rng(3)
    true_ang = 130
    img = simulate_image(Q, true_ang, sigma_sl1.value, rng)
    ref_name = ref_dd.value
    ref = _REFS[ref_name]; stack = _STACKS5[ref_name]
    scores = correlation(img, stack)
    best_ang = _CAND5[scores.argmax()]

    fig, axes = plt.subplots(1, 4, figsize=(14, 3.5))
    show(ref, axes[0], title=f"Reference\n({ref_name})")
    v = max(abs(img.min()), abs(img.max()))
    show(img, axes[1], vmin=-v, vmax=v,
         title=f"Noisy image (σ={sigma_sl1.value:.1f})\ntrue angle: {true_ang}°")
    axes[2].plot(_CAND5, scores, color='steelblue', linewidth=1.4)
    axes[2].axvline(true_ang, color='red',    linewidth=1.4, ls='--', label=f'True {true_ang}°')
    axes[2].axvline(best_ang, color='orange', linewidth=1.4, ls=':',  label=f'Best {best_ang}°')
    axes[2].set_xlabel("Candidate angle (°)"); axes[2].set_ylabel("Correlation")
    axes[2].set_title("Alignment score curve"); axes[2].legend(fontsize=8); axes[2].grid(alpha=0.3)
    aligned = nd_rotate(img, -best_ang, reshape=False)
    show(aligned, axes[3], title=f"Image rotated back {best_ang}°")
    with out_ref:
        clear_output(wait=True); display(fig2img(fig))

for w in [ref_dd, sigma_sl1]:
    w.observe(_update_ref, names='value')
display(VBox([ref_dd, sigma_sl1, out_ref]))
_update_ref()

13.5. Part 2 – Alignment algorithm#

For a set of \(N\) images \(\{X_i\}\), we align each to the reference by searching over discrete candidate angles \(\Delta\theta\) apart. For each image we rotate the reference to all candidate angles and pick the best-matching one:

\[ \hat\theta_i = \arg\max_{\theta_j} \bigl(X_i \cdot R^{\theta_j} A_\text{ref}\bigr) \]

After estimating all angles, the scatter plot of true vs estimated angles reveals how well alignment worked. A tight diagonal indicates good alignment; spread indicates confusion.

Task

Part 2, step 1: Implement the alignment loop in the code cell below. For each image, compute the correlation against every rotated reference and record the best angle. Run the cell to see the scatter plot of true vs estimated angles.

13.6. Part 3 – Reconstruction#

Given estimated angles, rotating each image back and averaging gives the reconstruction:

\[ A \leftarrow \frac{1}{N} \sum_{i=1}^N (R^{\hat\theta_i})^{-1} X_i \]

Task

Part 3: Complete the reconstruction: for each image, rotate it back by \(-\hat\theta_i\) and sum all rotated images. Divide by \(N\) to get the average.

13.7. Interactive reconstruction#

The widget below runs the complete pipeline — from noisy images to reconstruction — with adjustable parameters. It also shows the scatter of estimated vs true angles and the correlation between reconstruction and true image as a quality metric.

Show code cell source

Hide code cell source

_sl3 = Layout(width="440px")
n_sl_r    = IntSlider(  value=100,  min=10,  max=500, step=10,
                         description="# images",       style=_sty, layout=_sl3)
sigma_sl_r= FloatSlider(value=2.0,  min=0.0, max=10.0, step=0.5,
                         description="Noise σ",        style=_sty, layout=_sl3)
blur_sl_r = FloatSlider(value=5.0,  min=0.0, max=15.0, step=1.0,
                         description="LP blur σ",      style=_sty, layout=_sl3)
da_sl_r   = IntSlider(  value=5,    min=1,   max=20,  step=1,
                         description="Δangle (°)",     style=_sty, layout=_sl3)

out_recon = Output()
def _update_recon(_=None):
    sigma = sigma_sl_r.value; N = n_sl_r.value
    blur  = blur_sl_r.value;  da = da_sl_r.value
    true_angs, imgs = simulate_images(Q, sigma, N, seed=42)
    ref_init = low_pass_filter(Q, sigma=blur) if blur > 0 else Q.copy()
    c_angs, r_stack = rotate_reference(ref_init, delta_angle=da)
    est_angs = align_images(imgs, r_stack, c_angs)
    recon    = reconstruct(imgs, est_angs)
    # quality: correlation with true image
    qc = float((recon * Q).mean() /
               (np.sqrt((recon**2).mean() * (Q**2).mean()) + 1e-10))
    err = np.abs(((est_angs - true_angs + 180) % 360) - 180)

    fig, axes = plt.subplots(1, 4, figsize=(16, 3.8))
    show(ref_init, axes[0], title="Initial reference")
    show(recon,    axes[1], title=f"Reconstruction (CC={qc:.2f})")
    axes[2].scatter(true_angs, est_angs, s=1, c='steelblue', alpha=0.5)
    axes[2].plot([0,360],[0,360],'r--',lw=1)
    axes[2].set_xlabel("True angle (°)"); axes[2].set_ylabel("Estimated (°)")
    axes[2].set_title(f"Angle accuracy (med. err={np.median(err):.0f}°)")
    axes[3].hist(err, bins=36, color='steelblue', ec='k', lw=0.4)
    axes[3].set_xlabel("Angle error (°)"); axes[3].set_ylabel("Count")
    axes[3].set_title("Error histogram")
    plt.tight_layout()
    with out_recon:
        clear_output(wait=True); display(fig2img(fig))

for w in [n_sl_r, sigma_sl_r, blur_sl_r, da_sl_r]:
    w.observe(_update_recon, names='value')
display(VBox([n_sl_r, sigma_sl_r, blur_sl_r, da_sl_r, out_recon]))
_update_recon()

13.8. Bonus 1 – Adding translation#

So far we assumed particles are perfectly centred. In reality, particle positions within the extracted box vary. The model extends to:

\[ X_i = T^{t_i}(R^{\theta_i} A) + \sigma G_i \]

where \(T^{t_i}\) is a translation by vector \(\mathbf{t}_i = (dx_i, dy_i)\). We can handle integer-pixel translations by rolling the image array. The search grid is now \(N_\theta \times N_{dx} \times N_{dy}\), so runtime increases by a factor of \((2 t_\text{max}+1)^2\).

Show code cell source

Hide code cell source

def build_rotation_translation_stack(reference, delta_angle=10, max_t=8, t_step=4):
    """
    Build a reference stack covering all combinations of rotation and translation.
    Returns (stack, angle_list, dx_list, dy_list).
    """
    angles = np.arange(0, 360, delta_angle)
    ts     = np.arange(-max_t, max_t+1, t_step)
    stack, ang_list, dx_list, dy_list = [], [], [], []
    for a in angles:
        rot = nd_rotate(reference, a, reshape=False)
        for dx in ts:
            for dy in ts:
                stack.append(translate_image(rot, dx, dy))
                ang_list.append(a); dx_list.append(dx); dy_list.append(dy)
    return np.array(stack), np.array(ang_list), np.array(dx_list), np.array(dy_list)

def align_with_translation(images, ref_stack, ang_list, dx_list, dy_list):
    """Vectorised correlation alignment over rotation × translation grid."""
    n, h, w = images.shape
    R = ref_stack.shape[0]
    A = images.reshape(n, h*w).astype(float)
    B = ref_stack.reshape(R, h*w).astype(float)
    cc = (A @ B.T) / (h * w)
    best = cc.argmax(1)
    return ang_list[best], dx_list[best], dy_list[best]

def reconstruct_from_rotation_translation(images, est_angles, est_dx, est_dy):
    recon = np.zeros_like(images[0])
    for img, ang, dx, dy in zip(images, est_angles, est_dx, est_dy):
        undone = translate_image(img, -dx, -dy)
        recon += nd_rotate(undone, -ang, reshape=False)
    return recon / len(images)

# Run
N_t  = 80; sigma_t = 1.0; max_t = 10
imgs_t, true_ang_t, true_dx, true_dy = simulate_images_with_translation(
    Q, sigma_t, N_t, max_t=max_t, seed=42)
ref_t = low_pass_filter(Q, sigma=5)
stack_t, ang_l, dx_l, dy_l = build_rotation_translation_stack(ref_t, delta_angle=10, max_t=max_t, t_step=4)
est_ang_t, est_dx_t, est_dy_t = align_with_translation(imgs_t, stack_t, ang_l, dx_l, dy_l)
recon_t = reconstruct_from_rotation_translation(imgs_t, est_ang_t, est_dx_t, est_dy_t)

fig, axes = plt.subplots(1, 3, figsize=(10, 3.5))
show(Q,       axes[0], title="True A")
show(imgs_t[0], axes[1], vmin=imgs_t[0].min(), vmax=imgs_t[0].max(),
     title=f"Sample image (σ={sigma_t})")
show(recon_t, axes[2], title=f"Reconstruction with translation\n(N={N_t})")
plt.tight_layout()
display(fig2img(fig))

13.9. Bonus 2 – Bayesian maximum-likelihood alignment#

Hard assignment to the single best-fitting angle ignores alignment uncertainty. The maximum-likelihood approach instead computes a probability for each angle:

\[ P(X_i \mid \theta_j, A) \propto \exp\!\left[-\frac{\|X_i - R^{\theta_j}A\|^2}{2\sigma^2}\right] \]

The soft reconstruction then uses all angles weighted by their posterior:

\[ w_{ij} = \frac{P(X_i \mid \theta_j, A)}{\sum_{j'} P(X_i \mid \theta_{j'}, A)}, \qquad A \leftarrow \sum_i \sum_j w_{ij} (R^{\theta_j})^{-1} X_i \]

At high SNR the weights collapse onto the single best angle (recovering hard assignment); at low SNR they spread across many angles, effectively averaging out alignment noise.

13.11. Bonus 4 – 2D Classification#

Real datasets contain multiple particle types (different proteins, conformations, or viewing directions). 2D classification assigns each image to one of \(K\) classes simultaneously with the alignment:

Maintain \(K\) references \(\{A_1, \ldots, A_K\}\)
For each image, find the best (class, angle) pair by correlation
Reconstruct each class from the images assigned to it
Update references and iterate

The interactive below mixes \(K\) letters (Q, P, R) and attempts to recover them. The classification accuracy depends on the noise level and the number of images per class.

Show code cell source

Hide code cell source

K_sl     = IntSlider(  value=2,   min=2,  max=3,   step=1,
                        description="# classes K",  style=_sty, layout=_sl3)
nper_sl  = IntSlider(  value=40,  min=10, max=100, step=5,
                        description="Images/class",  style=_sty, layout=_sl3)
sigma_sl3= FloatSlider(value=1.5, min=0.0, max=6.0, step=0.5,
                        description="Noise σ",       style=_sty, layout=_sl3)

_CLS_LETTERS = ['Q', 'P', 'R']
_CLS_IMGS    = {lt: make_letter(lt, 64) for lt in _CLS_LETTERS}
_CLS_CAND    = np.arange(0, 360, 10)

def _precompute_cls_stack(letter):
    return np.array([nd_rotate(_CLS_IMGS[letter], a, reshape=False) for a in _CLS_CAND])

_CLS_STACKS = {lt: _precompute_cls_stack(lt) for lt in _CLS_LETTERS}

out_cls4 = Output()
def _update_cls4(_=None):
    K    = K_sl.value; n_per = nper_sl.value; sig = sigma_sl3.value
    letters = _CLS_LETTERS[:K]
    rng  = np.random.default_rng(42)

    # Generate mixed dataset
    all_imgs, true_labels = [], []
    for k, lt in enumerate(letters):
        idx = rng.integers(0, len(_CLS_CAND), n_per)
        noisy = _CLS_STACKS[lt][idx] + rng.standard_normal((n_per, 64, 64)) * sig
        all_imgs.append(noisy); true_labels.extend([k]*n_per)
    all_imgs   = np.concatenate(all_imgs)
    true_labels= np.array(true_labels)
    shuf = rng.permutation(len(all_imgs))
    all_imgs = all_imgs[shuf]; true_labels = true_labels[shuf]

    # Build combined reference stack
    ref_stack_all = np.concatenate([_CLS_STACKS[lt] for lt in letters])
    n_cand_per = len(_CLS_CAND)
    n_all = all_imgs.shape[0]
    A = all_imgs.reshape(n_all, 64*64).astype(float)
    B = ref_stack_all.reshape(len(ref_stack_all), 64*64).astype(float)
    cc = (A @ B.T) / (64*64)
    best_global  = cc.argmax(1)
    pred_class   = best_global // n_cand_per
    best_ang_idx = best_global %  n_cand_per

    acc = (pred_class == true_labels).mean() * 100
    class_recons = []
    for k in range(K):
        mask = (pred_class == k)
        recon_k = np.zeros((64, 64))
        count = 0
        for i in np.where(mask)[0]:
            recon_k += nd_rotate(all_imgs[i], -_CLS_CAND[best_ang_idx[i]], reshape=False)
            count   += 1
        class_recons.append(recon_k / max(count, 1))

    naive_avg = all_imgs.mean(0)
    fig, axes = plt.subplots(1, K+3, figsize=(3.5*(K+3), 3.5))
    show(all_imgs[0], axes[0], vmin=all_imgs[0].min(), vmax=all_imgs[0].max(),
         title="Example image\n(mixed)")
    show(naive_avg, axes[1], title="Naive average\n(blurry)")
    for k, (recon_k, lt) in enumerate(zip(class_recons, letters)):
        n_in_cls = (pred_class==k).sum()
        axes[k+2].imshow(recon_k, cmap='gray', origin='lower')
        axes[k+2].axis('off')
        axes[k+2].set_title(f"Class {k+1}: {lt}\n({n_in_cls} imgs assigned)", fontsize=8)
    for k, lt in enumerate(letters):
        axes[K+2].imshow(_CLS_IMGS[lt], cmap='gray', origin='lower',
                         extent=[k, k+0.9, 0, 0.9], aspect='auto')
    axes[K+2].set_xlim(-0.1, K); axes[K+2].set_ylim(-0.05, 1)
    axes[K+2].axis('off'); axes[K+2].set_title(f"True classes\nAcc={acc:.0f}%", fontsize=8)
    plt.suptitle(f"2D Classification: K={K} classes, σ={sig}, {n_per} imgs/class", fontsize=9)
    plt.tight_layout()
    with out_cls4:
        clear_output(wait=True); display(fig2img(fig))

for w in [K_sl, nper_sl, sigma_sl3]:
    w.observe(_update_cls4, names='value')
display(VBox([K_sl, nper_sl, sigma_sl3, out_cls4]))
_update_cls4()

13.12. Summary#

Topic	Key concept
Observation model	\(X_i = R^{\theta_i}A + \sigma G_i\)
Reference selection	Circular, blurred prior, or oracle; avoid reference bias
Alignment	Maximise correlation \(\langle X_i, R^{\theta_j}A\rangle\) over discrete candidates
Reconstruction	Average back-rotated images
Model bias	Wrong reference → biased reconstruction; correct with iterative refinement
ML alignment	Soft weights via \(\exp(-\|X_i - R^\theta A\|^2/2\sigma^2)\); better at low SNR
Iterative refinement	Use reconstruction as new reference; converges toward true structure
Classification	\(K\)-class joint alignment+assignment; separates heterogeneous datasets

Practical: Single-Particle Analysis

Contents

13. Practical: Single-Particle Analysis#

13.1. The reconstruction problem#

13.1.1. Helper functions#

13.2. Part 0 – The test image#

13.3. Simulating cryo-EM images#

13.4. Part 1 – Choosing an initial reference#

13.5. Part 2 – Alignment algorithm#

13.6. Part 3 – Reconstruction#

13.7. Interactive reconstruction#

13.8. Bonus 1 – Adding translation#

13.9. Bonus 2 – Bayesian maximum-likelihood alignment#

13.10. Bonus 3 – Iterative refinement#

13.11. Bonus 4 – 2D Classification#

13.12. Summary#