uchrom.recon

Single-cell reconstruction

EMber: noise-aware single-cell Hi-C 3-D genome reconstruction.

EMber (EM + annealing / sampling) is the native engine’s protocol "robust" (uchrom_nucdyn; native/nucdyn/src/robust.rs): the NucDynamics 2017 hierarchical annealing with a mixture observation model — every contact is a true proximity or a false contact (rate ε), contact restraints are re-weighted by their posterior probability of being true (EM; ε and the contact-kernel width learned from the data, plus structure-free local-support evidence) — a short schedule (2,500 MD steps per stage), a sampling stage at finite temperature (calibrated ensembles) and adaptive resolution below 25 kb.

It was pre-registered, developed on a development set, frozen (tag robust-frozen) and tested once against NucDynamics (benchmarks/screcon/PREREG.md, TEST_RESULTS.md). The defaults below are the frozen method.

from uchrom.recon.sc import reconstruct_ember cd = reconstruct_ember(“cell1.pairs”, n_models=10, device=”auto”) cd.layers[“model_0”] # one model per layer cd.results[“ember.contact_weights”] # posterior weight of every input contact

uchrom.recon.sc.ember.EMBER_PROTOCOL = 'robust'

engine protocol implementing EMber

uchrom.recon.sc.ember.reconstruct_ember(contacts: str | PathLike | DataFrame, *, chrom: None | str | Sequence[str] = None, n_models: int = 10, params: Mapping[str, Any] | None = None, device: str = 'auto', n_threads: int = 0, seed: int = 0, particle_sizes: Sequence[float] = (8000000.0, 4000000.0, 2000000.0, 400000.0, 200000.0, 100000.0), cell_id: str | None = None, key_added: str | None = 'ember', copy: bool = False, verbose: bool = False, **engine_params: Any)

Single-cell Hi-C genome structure ensemble with EMber.

Parameters:

• contacts – Path (.ncc, .pairs[.gz], GEO chr_A pos_A chr_B pos_B table) or DataFrame of contacts (see uchrom.recon.sc.nucdyn.load_contacts()).

• chrom – None (default): the whole genome. A chromosome name / range ("chr19", "chr1:1000000-5000000") or a list of them restricts the contacts to those regions (calling convention: chrom=None runs everything).

• n_models – Models (structures) in the ensemble.

• params – Engine parameters overriding the frozen defaults (see uchrom_nucdyn.default_params(); the robust_* fields are EMber’s); keyword arguments are merged into it.

• device – "auto" (GPU when available), "cpu" (f64, multi-core) or "gpu" (f32, wgpu).

• n_threads – As uchrom.recon.sc.nucdyn.reconstruct_nucdyn().

• seed – As uchrom.recon.sc.nucdyn.reconstruct_nucdyn().

• particle_sizes – As uchrom.recon.sc.nucdyn.reconstruct_nucdyn().

• cell_id – As uchrom.recon.sc.nucdyn.reconstruct_nucdyn().

• verbose – As uchrom.recon.sc.nucdyn.reconstruct_nucdyn().

• key_added – Key of cd.uns[key_added] (method, engine, protocol, parameters) and of the results <key>.stages / <key>.em (EM log) / <key>.contact_weights; None stores neither.

• copy – Reconstruction always creates a new ChromData from the contacts; accepted for the calling convention (no effect).

Returns:

One spot per particle; coords = model 0, layers['model_<k>'] = every model; units: particle radii. Same layout as uchrom.recon.sc.nucdyn.reconstruct_nucdyn() with protocol="robust" (identical coordinates for the same seed).

Return type:

ChromData

NucDynamics: single-cell Hi-C genome structure calculation.

reconstruct_nucdyn(contacts, ...) returns a ChromData ensemble computed by the native engine (uchrom_nucdyn: Rust, multi-core CPU and wgpu GPU; pip install u-chrom[nucdyn]). main (also python -m uchrom.recon.sc.nucdyn) is the file-level entry point; it falls back to the deprecated Taichi port when the native engine is not installed.

uchrom.recon.sc.nucdyn.engine_info() → Dict[str, Any]

Version of the native engine and the GPU it would use.

uchrom.recon.sc.nucdyn.load_contacts(contacts, genome_ranges=None) → Dict[str, ndarray]

Contacts as arrays chrom_a, pos_a, chrom_b, pos_b, ambiguity, active.

contacts is a path (.ncc as written by NucProcess; .pairs / .pairs.gz; a whitespace table chr_A pos_A chr_B pos_B such as the GEO files of Stevens et al. 2017) or a DataFrame with columns chr1, pos1, chr2, pos2 (as uchrom.io.read_pairs() returns) or chrom1, pos1, chrom2, pos2; optional ambiguity / active. genome_ranges keeps only contacts with both ends in the given ranges ("chr19", "chr1:1000000-5000000", or a list).

uchrom.recon.sc.nucdyn.main(in_file, out_file, arch='gpu', device_memory_fraction=0.9, cell_id=None, engine='auto', n_models=None, device=None, n_threads=0, seed=None, protocol='nuc_dynamics_2017', size_steps=None, genome_ranges=None, **kwargs)

Calculate a single-cell genome structure from contacts and write it.

in_file: contacts (.pairs[.gz], .ncc, GEO contact table; the Taichi port also reads .cool). out_file: .chromdata.zarr / .cdz / .h5cd (the whole ensemble: coords = model 0, layers['model_<k>'] = every model) or .csv (model 0 only).

engine: "auto" uses the native engine (uchrom_nucdyn) when it is installed, otherwise the deprecated Taichi port; "native" / "taichi" force one. device (auto / cpu / gpu) defaults from arch (cpu -> cpu; gpu / cuda / metal -> gpu). size_steps: particle sizes in Mb (default 8 4 2 0.4 0.2 0.1). Other keyword arguments: native engine parameters, or the Taichi port’s (dyns, hot, cold, random_seed, … are mapped to the native ones).

uchrom.recon.sc.nucdyn.native_available() → bool

True when the native engine (uchrom_nucdyn) can be imported.

uchrom.recon.sc.nucdyn.reconstruct_nucdyn(contacts: str | PathLike | DataFrame, *, n_models: int = 10, device: str = 'auto', n_threads: int = 0, seed: int = 0, protocol: str = 'nuc_dynamics_2017', particle_sizes: Sequence[float] = (8000000.0, 4000000.0, 2000000.0, 400000.0, 200000.0, 100000.0), genome_ranges=None, cell_id: str | None = None, key_added: str | None = 'nucdyn', params: Mapping[str, Any] | None = None, verbose: bool = False, _method: str = 'nucdyn', **engine_params: Any)

Single-cell Hi-C genome structure ensemble with NucDynamics.

Runs the hierarchical simulated-annealing protocol of NucDynamics (Stevens et al. 2017, Nature 544:59) on the native engine and returns the ensemble as a new ChromData.

Parameters:

• contacts – Path (.ncc, .pairs[.gz], GEO chr_A pos_A chr_B pos_B table) or DataFrame of contacts, see load_contacts().

• n_models – Number of structures (models) in the ensemble. All models are computed at once (in parallel on the CPU, in one batch on the GPU).

• device – "auto" (GPU when available, else CPU), "cpu" (f64, multi-core) or "gpu" (f32, wgpu: Metal / Vulkan / DX12).

• n_threads – CPU worker threads (0 = all cores).

• seed – Random seed (start coordinates and velocities); on the CPU a seed gives bitwise-identical results.

• protocol – "nuc_dynamics_2017" (default): the code that produced the published 2017 structures (tjs23/nuc_dynamics master). "release_1.3": the later bead-size-scaled version with ambiguity resolution and model selection (2 x n_models above 1 Mb, the n_models closest to the mean kept). "robust" (public name EMber, alias "ember"; see uchrom.recon.sc.ember.reconstruct_ember()): the 2017 protocol with a noise-aware observation model (contact restraints re-weighted by their posterior probability of being true contacts; false-contact rate and contact-kernel width learned by EM; robust_* engine parameters). Adds the results <key_added>.em (EM log) and <key_added>.contact_weights (final weight of every input contact, NaN = filtered out).

• particle_sizes – Hierarchical particle sizes in bp (coarse to fine).

• genome_ranges – Restrict the contacts to these ranges (e.g. "chr19").

• cell_id – Written to spots['cell_id'].

• key_added – Key of the provenance record in cd.results (per-stage log) and of cd.uns[key_added]; None stores neither.

• params – Further engine parameters (see uchrom_nucdyn.default_params(): temp_steps, dynamics_steps, temp_start, temp_end, time_step, contact_dist_lower, …); keyword arguments are merged into it.

Returns:

One spot per particle; one trace per chromosome. coords holds model 0 and layers['model_<k>'] holds every model k (model 0 included), so the whole ensemble travels with the object. Units are particle radii. Bins are the particles’ genomic intervals: in the 2017 protocol a particle at sequence position p collects the contacts in (p - size, p] (bins.start = p - size, bins.end = p); in release_1.3 [p, p + size). uns[key_added] records the engine, protocol, parameters and per-stage log.

Return type:

ChromData

Bulk reconstruction (MDS)

uchrom.recon.bulk.mds.apply_distance_decay_prior(contact_mat, weight=0.05)

Apply distance decay prior to smooth contact frequencies. Expected values computed from nonzero contacts only (matching miniMDS).

uchrom.recon.bulk.mds.apply_transform(coords, rotation=None, translation=None, scale=1.0)

Apply rotation, translation and scaling.

uchrom.recon.bulk.mds.compute_radius_of_gyration(coords)

Rg = sqrt(mean(||x - centroid||^2)).

uchrom.recon.bulk.mds.compute_stress(coords, dist_mat, weights=None)

Compute MDS stress, skipping missing-data pairs (dist_mat == 0).

uchrom.recon.bulk.mds.contact_to_distance(contact_mat, alpha=4.0)

Convert contact frequencies to distances: d = c^(-1/alpha). Zero contacts are treated as missing data (distance = 0).

uchrom.recon.bulk.mds.contacts_to_matrix(bin1, bin2, counts, n, dtype=torch.float64)

Build symmetric contact matrix from sparse contact data.

uchrom.recon.bulk.mds.fill_missing_distances(dist_mat, contact_mat=None)

Fill zero (missing) distances using genomic distance prior. After contact_to_distance, zeros represent missing data, not zero distance.

uchrom.recon.bulk.mds.inter_mds(input_path, resolution_inter=1000000, resolution_intra=100000, chroms=None, alpha=4.0, weight=0.05, n_iter=1000, device='auto', output_dir=None, verbose=True)

Whole-genome 3D reconstruction with inter-chromosomal contacts.

Parameters:

• input_path – Path to .hic or .mcool file

• resolution_inter – Resolution for inter-chromosomal scaffold (default 1Mb)

• resolution_intra – Resolution for intra-chromosomal structures (default 100kb)

• chroms – List of chromosomes (default: autosomes + X)

• alpha – Contact-to-distance exponent

• weight – Distance decay prior weight

• n_iter – MDS iterations

• device – ‘auto’, ‘cpu’, ‘cuda’, ‘mps’

• output_dir – Output directory (None = don’t save)

• verbose – Print progress

Returns:

DataFrame with chrom, start, end, x, y, z for all bins

Return type:

genome_df

uchrom.recon.bulk.mds.normalize_distances(dist_mat)

Normalize distance matrix to have unit mean. Includes zeros in mean calculation to match miniMDS behavior (miniMDS divides by np.mean(distMat) which includes zeros).

uchrom.recon.bulk.mds.partitioned_mds(contact_mat, tad_regions=None, device='auto', res_ratio=10, alpha=4.0, alpha2=2.5, weight=0.05, n_iter=1000, verbose=False, n_workers=1)

Partitioned MDS for high-resolution Hi-C data.

uchrom.recon.bulk.mds.procrustes_alignment(source, target, scale=True)

Align source to target using SVD-based Procrustes analysis.

uchrom.recon.bulk.mds.run_mds(contact_mat, alpha=4.0, device='auto', weight=0.05, **kwargs)

Full MDS pipeline: contact matrix -> 3D coordinates.

Zero-contact bins (rows/columns with no observed contacts) are removed before MDS, matching miniMDS behavior. Returns coordinates only for non-zero bins.

Returns:

np.ndarray of shape (n_nonzero, 3) nonzero_mask: np.ndarray boolean mask of shape (n_total,)

indicating which bins were kept

Return type:

coords

uchrom.recon.bulk.mds.torch_mds(dist_mat, device='auto', n_iter=1000, lr=0.01, tol=1e-06, init='cmds', verbose=False, method='smacof')

Run iterative MDS.

Parameters:

method – ‘smacof’ (default, fast) or ‘adam’ (gradient descent).

uchrom.recon.bulk.mds.torch_mds.cmds_init(dist_mat)

Classical MDS initialization via eigendecomposition.

uchrom.recon.bulk.mds.torch_mds.compute_stress(coords, dist_mat, weights=None)

Compute MDS stress, skipping missing-data pairs (dist_mat == 0).

uchrom.recon.bulk.mds.torch_mds.get_device(device='auto')

Get appropriate torch device.

uchrom.recon.bulk.mds.torch_mds.get_dtype(device)

MPS only supports float32.

uchrom.recon.bulk.mds.torch_mds.run_mds(contact_mat, alpha=4.0, device='auto', weight=0.05, **kwargs)

Full MDS pipeline: contact matrix -> 3D coordinates.

Zero-contact bins (rows/columns with no observed contacts) are removed before MDS, matching miniMDS behavior. Returns coordinates only for non-zero bins.

Returns:

np.ndarray of shape (n_nonzero, 3) nonzero_mask: np.ndarray boolean mask of shape (n_total,)

indicating which bins were kept

Return type:

coords

uchrom.recon.bulk.mds.torch_mds.smacof(dist_mat, device='auto', n_iter=1000, tol=1e-06, init='cmds', verbose=False)

Run SMACOF (Scaling by MAjorizing a Complicated Function) MDS.

Unlike the Adam-based approach, SMACOF uses a majorization algorithm that does not require autograd, resulting in much lower per-iteration overhead on CPU.

uchrom.recon.bulk.mds.torch_mds.torch_mds(dist_mat, device='auto', n_iter=1000, lr=0.01, tol=1e-06, init='cmds', verbose=False, method='smacof')

Run iterative MDS.

Parameters:

method – ‘smacof’ (default, fast) or ‘adam’ (gradient descent).

uchrom.recon.bulk.mds.inter.inter_mds(input_path, resolution_inter=1000000, resolution_intra=100000, chroms=None, alpha=4.0, weight=0.05, n_iter=1000, device='auto', output_dir=None, verbose=True)

Whole-genome 3D reconstruction with inter-chromosomal contacts.

Parameters:

• input_path – Path to .hic or .mcool file

• resolution_inter – Resolution for inter-chromosomal scaffold (default 1Mb)

• resolution_intra – Resolution for intra-chromosomal structures (default 100kb)

• chroms – List of chromosomes (default: autosomes + X)

• alpha – Contact-to-distance exponent

• weight – Distance decay prior weight

• n_iter – MDS iterations

• device – ‘auto’, ‘cpu’, ‘cuda’, ‘mps’

• output_dir – Output directory (None = don’t save)

• verbose – Print progress

Returns:

DataFrame with chrom, start, end, x, y, z for all bins

Return type:

genome_df

uchrom.recon.bulk.mds.transforms.align_substructure_to_scaffold(high_res_coords, low_res_coords, scaffold_coords, res_ratio=10)

Align high-res substructure to global scaffold via Procrustes.

uchrom.recon.bulk.mds.transforms.apply_transform(coords, rotation=None, translation=None, scale=1.0)

Apply rotation, translation and scaling.

uchrom.recon.bulk.mds.transforms.center_coords(coords)

Center coordinates to zero mean.

uchrom.recon.bulk.mds.transforms.compute_radius_of_gyration(coords)

Rg = sqrt(mean(||x - centroid||^2)).

uchrom.recon.bulk.mds.transforms.downsample_coords(coords, res_ratio, method='mean')

Downsample coordinates by resolution ratio.

uchrom.recon.bulk.mds.transforms.procrustes_alignment(source, target, scale=True)

Align source to target using SVD-based Procrustes analysis.