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_recon; packages/uchrom-recon/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

Phased (homolog-resolved) contacts — hickit .pairs with phase0 / phase1, Dip-C .con — give a diploid reconstruction: one particle chain per chromosome copy, the copy pair of every contact with an unphased end inferred by EM from the structure (the copy pair is a latent variable of the same mixture model):

cd = reconstruct_ember("cell.contacts.pairs.gz", ploidy="auto")   # traces chr1(pat), chr1(mat), ...
cd.results["ember.imputed_phase"]           # most probable (hap_a, hap_b) of every contact

The diploid extension is new (not part of the frozen, pre-registered test); its haploid behaviour is unchanged.

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, ploidy: Any = None, hap_names: Sequence[str] = ('pat', 'mat'), phase_prior: 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_recon.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).

• ploidy – Phased (diploid) input, see uchrom.recon.sc.nucdyn.reconstruct_nucdyn(): copies per chromosome (None = "auto": male when chrY is present; "haploid" ignores the phases), haplotype names in the trace ids, and whether hickit’s imputed phase probabilities (when present) are used as an extra prior of the copy pairs (default False: the frozen configuration, EMber’s own neighbourhood prior only; benchmarks/screcon/PREREG.md §6). Adds the results <key>.phase_probs / <key>.imputed_phase and uns[key]["diploid"] (phasing statistics).

• hap_names – Phased (diploid) input, see uchrom.recon.sc.nucdyn.reconstruct_nucdyn(): copies per chromosome (None = "auto": male when chrY is present; "haploid" ignores the phases), haplotype names in the trace ids, and whether hickit’s imputed phase probabilities (when present) are used as an extra prior of the copy pairs (default False: the frozen configuration, EMber’s own neighbourhood prior only; benchmarks/screcon/PREREG.md §6). Adds the results <key>.phase_probs / <key>.imputed_phase and uns[key]["diploid"] (phasing statistics).

• phase_prior – Phased (diploid) input, see uchrom.recon.sc.nucdyn.reconstruct_nucdyn(): copies per chromosome (None = "auto": male when chrY is present; "haploid" ignores the phases), haplotype names in the trace ids, and whether hickit’s imputed phase probabilities (when present) are used as an extra prior of the copy pairs (default False: the frozen configuration, EMber’s own neighbourhood prior only; benchmarks/screcon/PREREG.md §6). Adds the results <key>.phase_probs / <key>.imputed_phase and uns[key]["diploid"] (phasing statistics).

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_recon: Rust, multi-core CPU and wgpu GPU; pip install u-chrom[recon]). main (also python -m uchrom.recon.sc.nucdyn) is the file-level entry point.

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).

Phased (homolog-resolved) contacts add hap_a / hap_b (int8: 0, 1, -1 = unknown) and, for hickit’s imputed pairs, phase_prob ((n, 4), index 2 hap_a + hap_b): files for which is_phased_file() is true (hickit .pairs with phase0 / phase1; Dip-C .con), or DataFrames with hap1 / hap2 (or phase0 / phase1; . = unknown) and optionally phase_prob00 .. phase_prob11 (uchrom.io.read_phased_pairs()).

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). out_file: .chromdata.zarr / .cdz / .h5cd (the whole ensemble: coords = model 0, layers['model_<k>'] = every model) or .csv (model 0 only).

engine: the native engine (uchrom_recon; "auto" / "native"). The Taichi port ("taichi") was removed: it needs u-chrom 0.2 from before the monorepo (commit 0063b92). 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 those of the former Taichi port (dyns, hot, cold, random_seed, … are mapped to the native ones).

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

True when the native engine (uchrom_recon) 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, ploidy: Any = None, hap_names: Sequence[str] = ('pat', 'mat'), phase_prior: 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_recon.default_params(): temp_steps, dynamics_steps, temp_start, temp_end, time_step, contact_dist_lower, …); keyword arguments are merged into it.

• ploidy – Phased (homolog-resolved) contacts — hickit .pairs with phase0 / phase1, Dip-C .con, a DataFrame with hap1 / hap2 (see load_contacts()) — are reconstructed as a diploid genome (protocols "nuc_dynamics_2017" and "robust" / EMber): one particle chain per chromosome copy, contacts with an unphased end ambiguous between the candidate copies (2017: NucDynamics ambiguous restraint; EMber: the copy pair is inferred by EM from the structure). ploidy gives the copies per chromosome: None / "auto" (male — X and Y single copy — when a Y chromosome is present, as hickit), "diploid", "male", a dict {chrom: 1 | 2}, or "haploid" (ignore the phases: one chain per chromosome).

• hap_names – Names of haplotypes 0 and 1 in the trace ids (default Dip-C’s ("pat", "mat"): traces chr1(pat) / chr1(mat); single-copy chromosomes keep the plain name).

• phase_prior – Also use hickit’s imputed phase_prob00 .. phase_prob11 (when present) as a prior over the candidate copy pairs (default False, the frozen diploid EMber configuration: EMber’s own 2-D neighbourhood prior robust_dip_nbr_prior only, no external prior).

Returns:

One spot per particle; one trace per chromosome (diploid: per chromosome copy, spot column haplotype = 0 / 1; both copies share the bins). 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.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.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.