API Reference

This section documents the internal modules of synth-cryo-em.

synth_cryo_em.core

apply_ctf(data, voxel_size, defoc=2.0, cs=2.7, voltage=300, amplitude_contrast=0.1, b_factor=0.0)

Apply a simple Contrast Transfer Function (CTF) to the 3D data.

Parameters:

Name	Type	Description	Default
data	ndarray	The input 3D density map.	required
voxel_size	tuple	Voxel size in Angstroms (x, y, z).	required
defoc	float	Defocus in micrometers. Defaults to 2.0.	2.0
cs	float	Spherical aberration in mm. Defaults to 2.7.	2.7
voltage	float	Acceleration voltage in kV. Defaults to 300.	300
amplitude_contrast	float	Amplitude contrast fraction. Defaults to 0.1.	0.1
b_factor	float	Envelope function B-factor. Defaults to 0.0.	0.0

Returns:

Type	Description
NDArray[float64]	numpy.ndarray: The CTF-corrupted 3D density map.

Source code in src/synth_cryo_em/core.py

def apply_ctf(
    data: npt.NDArray[np.float64],
    voxel_size: tuple[float, float, float],
    defoc: float = 2.0,
    cs: float = 2.7,
    voltage: float = 300,
    amplitude_contrast: float = 0.1,
    b_factor: float = 0.0,
) -> npt.NDArray[np.float64]:
    """
    Apply a simple Contrast Transfer Function (CTF) to the 3D data.

    Args:
        data (numpy.ndarray): The input 3D density map.
        voxel_size (tuple): Voxel size in Angstroms (x, y, z).
        defoc (float, optional): Defocus in micrometers. Defaults to 2.0.
        cs (float, optional): Spherical aberration in mm. Defaults to 2.7.
        voltage (float, optional): Acceleration voltage in kV. Defaults to 300.
        amplitude_contrast (float, optional): Amplitude contrast fraction. Defaults to 0.1.
        b_factor (float, optional): Envelope function B-factor. Defaults to 0.0.

    Returns:
        numpy.ndarray: The CTF-corrupted 3D density map.
    """
    if any(v <= 0 for v in voxel_size):
        raise ValueError(f"Voxel size dimensions must be positive, got {voxel_size}")
    if voltage <= 0:
        raise ValueError(f"Voltage must be positive, got {voltage}")
    if not (0 <= abs(amplitude_contrast) <= 1):
        raise ValueError(f"Amplitude contrast must be between -1 and 1, got {amplitude_contrast}")

    # Constants
    wl = 12.26 / np.sqrt(voltage * 1000 + 0.9784 * voltage**2)  # wavelength in Angstroms
    cs_a = cs * 1e7  # cs in Angstroms
    defoc_a = defoc * 10000  # defocus in Angstroms

    nz, ny, nx = data.shape
    # Frequencies
    kz = np.fft.fftfreq(nz, d=voxel_size[0])
    ky = np.fft.fftfreq(ny, d=voxel_size[1])
    kx = np.fft.fftfreq(nx, d=voxel_size[2])

    Kz, Ky, Kx = np.meshgrid(kz, ky, kx, indexing="ij")
    k2 = Kz**2 + Ky**2 + Kx**2

    # Phase shift
    chi = np.pi * wl * k2 * (defoc_a - 0.5 * wl**2 * k2 * cs_a)

    # CTF
    ctf = -(np.sqrt(1 - amplitude_contrast**2) * np.sin(chi) + amplitude_contrast * np.cos(chi))

    # Envelope function
    if b_factor > 0:
        envelope = np.exp(-b_factor * k2 / 4.0)
        ctf *= envelope

    # Apply in Fourier domain
    data_f = np.fft.fftn(data)
    data_f *= ctf
    return np.real(np.fft.ifftn(data_f))

compute_ccc(data1, data2)

Compute the Cross-Correlation Coefficient (CCC) between two 3D maps.

Parameters:

Name	Type	Description	Default
data1	ndarray	The first 3D density map.	required
data2	ndarray	The second 3D density map.	required

Returns:

Name	Type	Description
float	float	The Pearson cross-correlation coefficient between the two maps.

Source code in src/synth_cryo_em/core.py

def compute_ccc(data1: npt.NDArray[Any], data2: npt.NDArray[Any]) -> float:
    """
    Compute the Cross-Correlation Coefficient (CCC) between two 3D maps.

    Args:
        data1 (numpy.ndarray): The first 3D density map.
        data2 (numpy.ndarray): The second 3D density map.

    Returns:
        float: The Pearson cross-correlation coefficient between the two maps.
    """
    if data1.shape != data2.shape:
        raise ValueError(f"Maps have different shapes: {data1.shape} vs {data2.shape}")

    # Flatten and convert to float64 for precision
    d1 = data1.astype(np.float64).ravel()
    d2 = data2.astype(np.float64).ravel()

    d1 = d1 - np.mean(d1)
    d2 = d2 - np.mean(d2)

    num = np.sum(d1 * d2)
    den = np.sqrt(np.sum(d1**2) * np.sum(d2**2))

    if den == 0:
        return 0.0
    return float(np.clip(num / den, -1.0, 1.0))

generate_density_map(input_path, resolution, grid_spacing=None, use_bfactors=False, margin=None)

Generate a density map from an atomic model file (PDB, mmCIF, BCIF) using gemmi.

Parameters:

Name	Type	Description	Default
input_path	str	Path to the input structure file.	required
resolution	float	Target resolution in Angstroms.	required
grid_spacing	float	Grid spacing (pixel size) in Angstroms. Defaults to resolution / 3.0.	None
use_bfactors	bool	If True, use atomic B-factors for local resolution. Defaults to False.	False
margin	float	Margin around the structure in Angstroms. Defaults to resolution * 2.0.	None

Returns:

Name	Type	Description
tuple	tuple[FloatGrid, NDArray[float64]]	A tuple containing: - numpy.ndarray: The 3D density grid. - numpy.ndarray: The origin coordinate of the map.

Source code in src/synth_cryo_em/core.py

def generate_density_map(
    input_path: str,
    resolution: float,
    grid_spacing: float | None = None,
    use_bfactors: bool = False,
    margin: float | None = None,
) -> tuple[gemmi.FloatGrid, npt.NDArray[np.float64]]:
    """
    Generate a density map from an atomic model file (PDB, mmCIF, BCIF) using gemmi.

    Args:
        input_path (str): Path to the input structure file.
        resolution (float): Target resolution in Angstroms.
        grid_spacing (float, optional): Grid spacing (pixel size) in Angstroms. Defaults to resolution / 3.0.
        use_bfactors (bool, optional): If True, use atomic B-factors for local resolution. Defaults to False.
        margin (float, optional): Margin around the structure in Angstroms. Defaults to resolution * 2.0.

    Returns:
        tuple: A tuple containing:
            - numpy.ndarray: The 3D density grid.
            - numpy.ndarray: The origin coordinate of the map.
    """
    if resolution <= 0:
        raise ValueError(f"Resolution must be positive, got {resolution}")
    if grid_spacing is not None and grid_spacing <= 0:
        raise ValueError(f"Grid spacing must be positive, got {grid_spacing}")
    if margin is not None and margin < 0:
        raise ValueError(f"Margin must be non-negative, got {margin}")

    st = gemmi.read_structure(input_path)
    # If grid_spacing is not provided, use a rule of thumb (resolution / 3 or 4)
    if grid_spacing is None:
        grid_spacing = resolution / 3.0

    # Get all atomic positions
    positions_list: list[list[float]] = []
    for model in st:
        for chain in model:
            for residue in chain:
                for atom in residue:
                    positions_list.append(atom.pos.tolist())

    if not positions_list:
        raise ValueError("No atoms found in structure")

    positions = np.array(positions_list)
    if margin is None:
        margin = resolution * 2.0

    min_pos = positions.min(axis=0) - margin
    max_pos = positions.max(axis=0) + margin
    size = max_pos - min_pos

    st_shifted = st.clone()
    # Ensure spacegroup is P1 for the synthetic box to avoid incorrect symmetry applications
    st_shifted.spacegroup_hm = "P1"
    cell = gemmi.UnitCell(size[0], size[1], size[2], 90, 90, 90)
    st_shifted.cell = cell

    for model in st_shifted:
        for chain in model:
            for residue in chain:
                for atom in residue:
                    atom.pos.x -= min_pos[0]
                    atom.pos.y -= min_pos[1]
                    atom.pos.z -= min_pos[2]

    # Map the atoms to the grid using DensityCalculatorE
    calc = gemmi.DensityCalculatorE()
    calc.d_min = resolution

    # Calculate sampling rate to match grid_spacing
    # spacing = d_min / (2 * rate)  => rate = d_min / (2 * spacing)
    calc.rate = resolution / (2.0 * grid_spacing)

    # Initialize grid
    calc.set_grid_cell_and_spacegroup(st_shifted)

    # If use_bfactors is True, we use the atomic B-factors.
    # Gemmi's DensityCalculatorE uses atomic B-factors by default
    # when calling put_model_density_on_grid.
    # However, we can add a constant "base" blur to represent resolution.
    # resolution (d) relates to B-factor roughly by B = 8 * pi^2 * (d/2)^2 = 2 * pi^2 * d^2
    # But gemmi also uses d_min as a cutoff.

    if not use_bfactors:
        # If not using B-factors, we set them all to 0 and use a global blur
        # equivalent to the target resolution.
        for model in st_shifted:
            for chain in model:
                for residue in chain:
                    for atom in residue:
                        atom.b_iso = 0.0
        # Set a global blur to match the target resolution
        # A common heuristic is B = 8 * res^2 for synthetic maps
        calc.blur = 8.0 * resolution**2

    calc.initialize_grid()
    if len(st_shifted) > 0:
        calc.put_model_density_on_grid(st_shifted[0])

    return calc.grid, min_pos

save_mrc(data, output_path, origin=(0, 0, 0), spacing=(1, 1, 1))

Save numpy array to an MRC file.

Parameters:

Name	Type	Description	Default
data	ndarray	The 3D density map to save.	required
output_path	str	Path to the output MRC file.	required
origin	tuple	The origin coordinates (x, y, z). Defaults to (0, 0, 0).	(0, 0, 0)
spacing	tuple	The voxel size in Angstroms (x, y, z). Defaults to (1, 1, 1).	(1, 1, 1)

Source code in src/synth_cryo_em/core.py

def save_mrc(
    data: npt.NDArray[Any],
    output_path: str,
    origin: npt.ArrayLike = (0, 0, 0),
    spacing: npt.ArrayLike = (1, 1, 1),
) -> None:
    """
    Save numpy array to an MRC file.

    Args:
        data (numpy.ndarray): The 3D density map to save.
        output_path (str): Path to the output MRC file.
        origin (tuple, optional): The origin coordinates (x, y, z). Defaults to (0, 0, 0).
        spacing (tuple, optional): The voxel size in Angstroms (x, y, z). Defaults to (1, 1, 1).
    """
    # Convert origin and spacing to 1D arrays for easier indexing if they aren't already
    origin_arr = np.asarray(origin)
    spacing_arr = np.asarray(spacing)

    with mrcfile.new(output_path, overwrite=True) as mrc:
        mrc.set_data(data.astype(np.float32))
        mrc.voxel_size = (float(spacing_arr[0]), float(spacing_arr[1]), float(spacing_arr[2]))
        # mrcfile uses x, y, z for origin
        mrc.header.origin.x = float(origin_arr[0])
        mrc.header.origin.y = float(origin_arr[1])
        mrc.header.origin.z = float(origin_arr[2])
        mrc.update_header_from_data()

synth_cryo_em.cli

main(input_path, output_path, resolution, spacing, snr, defocus, voltage, cs, bfactor, bfactors, apply_physics)

Generate a synthetic Cryo-EM map from an atomic model (PDB, mmCIF, or BCIF).

This function acts as the CLI entrypoint for map generation. It parses arguments and coordinates the voxelization, CTF physics, and noise simulation.

Source code in src/synth_cryo_em/cli.py

@click.command()
@click.argument("input_path", type=click.Path(exists=True))
@click.argument("output_path", type=click.Path())
@click.option("--resolution", "-r", default=4.0, help="Resolution in Angstroms")
@click.option("--spacing", "-s", default=None, type=float, help="Grid spacing in Angstroms")
@click.option("--snr", default=None, type=float, help="Signal-to-noise ratio")
@click.option("--defocus", default=2.0, help="Defocus in micrometers")
@click.option("--voltage", default=300.0, help="Acceleration voltage in kV")
@click.option("--cs", default=2.7, help="Spherical aberration in mm")
@click.option("--bfactor", default=0.0, help="Envelope B-factor")
@click.option(
    "--bfactors/--no-bfactors", default=False, help="Use atomic B-factors for local resolution"
)
@click.option("--apply-physics/--no-physics", default=False, help="Apply CTF effects")
def main(
    input_path: str,
    output_path: str,
    resolution: float,
    spacing: float | None,
    snr: float | None,
    defocus: float,
    voltage: float,
    cs: float,
    bfactor: float,
    bfactors: bool,
    apply_physics: bool,
) -> None:
    """
    Generate a synthetic Cryo-EM map from an atomic model (PDB, mmCIF, or BCIF).

    This function acts as the CLI entrypoint for map generation. It parses
    arguments and coordinates the voxelization, CTF physics, and noise simulation.
    """
    try:
        click.echo(f"Generating map for {input_path} at {resolution}A resolution...")

        grid, origin = generate_density_map(
            input_path, resolution, grid_spacing=spacing, use_bfactors=bfactors
        )

        data = np.array(grid, copy=True)

        # Voxel size is from the unit cell
        uc = grid.unit_cell
        vox_size = (uc.a / grid.nu, uc.b / grid.nv, uc.c / grid.nw)

        if apply_physics:
            click.echo(
                f"Applying CTF (defocus={defocus}um, voltage={voltage}kV, B-factor={bfactor})..."
            )
            data = apply_ctf(
                data, vox_size, defoc=defocus, cs=cs, voltage=voltage, b_factor=bfactor
            )

        if snr is not None:
            click.echo(f"Adding Gaussian noise (SNR={snr})...")
            data = add_gaussian_noise(data, snr)

        save_mrc(data, output_path, origin=origin, spacing=vox_size)
        click.echo(f"Saved synthetic map to {output_path}")
    except ValueError as e:
        click.echo(f"Error: {e}", err=True)
        raise click.Abort() from e

synth_cryo_em.validate

compute_and_report_fsc(freqs, fsc, ccc, output=None)

Report FSC and CCC results to the console and optionally to a file.

Source code in src/synth_cryo_em/validate.py

def compute_and_report_fsc(
    freqs: np.ndarray, fsc: np.ndarray, ccc: float, output: str | None = None
) -> None:
    """
    Report FSC and CCC results to the console and optionally to a file.
    """
    click.echo(f"Overall Cross-Correlation Coefficient (CCC): {ccc:.4f}")
    click.echo("\nFSC Summary:")
    click.echo(f"{'Resolution (A)':<15} | {'FSC':<10}")
    click.echo("-" * 30)
    # Ensure step is at least 1
    step = max(1, len(freqs) // 10)
    for i in range(0, len(freqs), step):
        res = 1.0 / freqs[i] if freqs[i] > 0 else float("inf")
        click.echo(f"{res:<15.2f} | {fsc[i]:<10.4f}")

    # Find 0.5 and 0.143 crossings
    for val in [0.5, 0.143]:
        cross_idx = np.where(fsc < val)[0]
        if len(cross_idx) > 0:
            res = 1.0 / freqs[cross_idx[0]]
            click.echo(f"\nFSC={val} crossing at {res:.2f} Angstroms")

    if output:
        np.savetxt(output, np.column_stack((freqs, fsc)), delimiter=",", header="frequency,fsc")
        click.echo(f"\nFSC data saved to {output}")

main(map1_path, map2_path, output)

Compare two Cryo-EM maps using Fourier Shell Correlation (FSC) and CCC.

This function acts as the CLI entrypoint for map validation. It loads two MRC maps, verifies their dimensions, and computes correlation metrics.

Source code in src/synth_cryo_em/validate.py

@click.command()
@click.argument("map1_path", type=click.Path(exists=True))
@click.argument("map2_path", type=click.Path(exists=True))
@click.option("--output", "-o", help="Path to save FSC data (CSV)")
def main(map1_path: str, map2_path: str, output: str | None) -> None:
    """
    Compare two Cryo-EM maps using Fourier Shell Correlation (FSC) and CCC.

    This function acts as the CLI entrypoint for map validation. It loads two
    MRC maps, verifies their dimensions, and computes correlation metrics.
    """
    try:
        with mrcfile.open(map1_path) as m1, mrcfile.open(map2_path) as m2:
            data1 = m1.data
            data2 = m2.data
            voxel_size = (m1.voxel_size.z, m1.voxel_size.y, m1.voxel_size.x)

        ccc = compute_ccc(data1, data2)
        freqs, fsc = compute_fsc(data1, data2, voxel_size)

        compute_and_report_fsc(freqs, fsc, ccc, output)
    except ValueError as e:
        click.echo(f"Error: {e}", err=True)
        raise click.Abort() from e