Skip to content

iscc.data.estimate_dna

Method-of-moments / 1-D MLE estimate_dna() for the DNA assay.

The DNA analogue of the scRNA estimate(): fits the technical hyper-parameters of the DNA batch model (DNABatchHyperParams in batch.py) from a real DNA-seq dataset reduced to counts-level summary statistics — per-locus coverage, alt/ref counts, and a per-locus called copy number — so realistic technical magnitudes are learned, not guessed. The fitted DNABatchHyperParams drive bulkDNA / scDNA / run_dna_batches directly (same field names, no redefinition).

This is not ABC: every DNA assay parameter maps directly to a marginal summary statistic, so estimation is closed-form / 1-D MLE — fast and simulation-free. ABC stays reserved for the biology.

Conditional on called copy number (the one wrinkle). Our depth model is copy-number-scaled (coverage ∝ CN × efficiency). Real data does not hand you true CN, so estimation is a two-pass: (1) run a CN caller (ASCAT/Sequenza for bulk; Ginkgo/HMMcopy for single-cell) to get per-locus CN, then (2) estimate the residual technical noise given that CN. We are explicit that the technical params are fit conditional on the caller's CN — exactly the CN regime re-injected when simulating. This is the DNA analogue of the scRNA estimator's library-size de-biasing of the dispersion: separate the technical layer from the part that is really biology.

What is fit (and from what), per modality:

mu_depth          mean per-locus coverage                              (bulk + sc)
kappa (DM)        DM concentration: coverage overdispersion AFTER       (bulk + sc)
                  dividing out CN x efficiency (index of dispersion)
nb_dispersion(NB) NB per-locus phi from the same CN-conditioned residual(bulk + sc; depth_model=nb)
gc_curve_sigma    curvature (log-space sd) of the coverage-vs-GC fit    (bulk + sc; needs GC)
capture_sigma     LogNormal sd of per-target/amplicon coverage after    (WES/panel only)
                  removing GC+CN (depth-noise de-biased)
error_rate        alt fraction at non-variant / hom-ref loci            (bulk + sc; needs mask)
depth_batch_sigma sd of log mean-depth across batches                   (>=2 batches)
ado_rate          fraction of known-het loci collapsing to ~0/~1 BAF    (single-cell only)
beta_binom_conc   BAF overdispersion at het loci beyond Binomial        (single-cell only)
ffpe_ct_rate      excess alt at C>T-eligible loci vs the error floor     (bulk + sc; needs ct mask)
doublet_rate      NOT identifiable from counts -> prior-only (flagged)  (carried from preset)
breadth, depth_model  chosen, not estimated (you know the assay)

The _PRIOR_ONLY escape that the scRNA estimate() uses (ambient/doublet) carries the genuinely unidentifiable fields honestly instead of pretending to fit them. ado_rate / beta_binom_conc are the per-cell amplification/dropout layer — fittable only from single-cell data; they are flagged not-fit (and carried from the preset) on bulk data, and kappa means different things per modality (large from pooled bulk, small from lumpy MDA/MALBAC single-cell).

numpy / scipy only (no JAX / R); stateless and config-driven.

Classes:

Name Description
DNAEstimate

Fitted DNA technical parameters, ready to drive bulkDNA / scDNA / run_dna_batches.

Functions:

Name Description
estimate_dna

Fit DNABatchHyperParams from counts-level DNA summary statistics (MoM/MLE).

estimate_dna_from_assay

Fit from a run bulkDNA / scDNA instance (extracts the reduced statistics).

DNAEstimate dataclass

DNAEstimate(hypers: DNABatchHyperParams, modality: str, breadth: str, depth_model: str, n_loci: int, n_cells: int, n_batches: int, fitted: list, diagnostics: dict = dict())

Fitted DNA technical parameters, ready to drive bulkDNA / scDNA / run_dna_batches.

hypers is a DNABatchHyperParams with the fitted magnitudes; dna_kwargs() returns the hyper fields (minus breadth / depth_model, which are explicit assay constructor args) as a kwargs dict, so the fit splats straight back into the assay without redefinition. fitted lists the fields actually learned from this data; everything else is carried from the (modality, breadth) preset and flagged honestly. diagnostics keeps realized intermediate quantities (per-sample kappa, realized GC/capture spread, het-BAF variance) for reporting.

Methods:

Name Description
dna_kwargs

Hyper-parameter overrides to splat into bulkDNA / scDNA (minus breadth/model).

dna_kwargs

dna_kwargs()

Hyper-parameter overrides to splat into bulkDNA / scDNA (minus breadth/model).

breadth and depth_model are explicit assay constructor args, so drive the assay as scDNA(breadth=est.breadth, depth_model=est.depth_model, **est.dna_kwargs()).

estimate_dna

estimate_dna(coverage, alt_counts, called_cn, *, modality='bulk', breadth='wgs', depth_model='dm', gc=None, mappability=None, ct_sites=None, variant_mask=None, het_mask=None, batch_depths=None)

Fit DNABatchHyperParams from counts-level DNA summary statistics (MoM/MLE).

Parameters:

Name Type Description Default
coverage array

Per-locus coverage. 1D (n_loci,) for bulk; 2D (n_cells, n_loci) for single-cell.

required
alt_counts array

Alt read counts, same shape as coverage.

required
called_cn array

Caller per-locus total copy number (the two-pass conditioning), broadcastable to coverage's shape. Estimation is explicitly conditional on this CN.

required
modality (bulk, sc)

Bulk vs single-cell. ado_rate / beta_binom_conc are fit only for "sc"; on bulk they are flagged not-fit and carried from the preset (and vice-versa is moot — they don't exist for bulk). kappa means different regimes per modality.

"bulk"
breadth (wgs, wes, panel)

Capture breadth (chosen, not estimated). capture_sigma is fit only for WES/panel.

"wgs"
depth_model (dm, nb)

Depth emission (chosen). Selects kappa (DM) vs nb_dispersion (NB) as the fit target.

"dm"
gc array | None

Reference per-locus GC content / mappability (n_loci,). GC enables gc_curve_sigma and (with CN) sharpens the dispersion / capture fits; mappability is divided out when given.

None
mappability array | None

Reference per-locus GC content / mappability (n_loci,). GC enables gc_curve_sigma and (with CN) sharpens the dispersion / capture fits; mappability is divided out when given.

None
ct_sites array | None

Boolean per-locus C>T-eligible mask -> ffpe_ct_rate.

None
variant_mask array | None

Boolean per-locus mask, True at known variant sites. Hom-ref loci (~variant_mask) give the error_rate floor.

None
het_mask array | None

Boolean per-locus mask of known-het loci (single-cell ado_rate / beta_binom_conc).

None
batch_depths sequence | None

Per-batch mean depths (>=2) -> depth_batch_sigma.

None

Returns:

Type Description
DNAEstimate

estimate_dna_from_assay

estimate_dna_from_assay(assay, *, batch_depths=None, het_atol=0.15)

Fit from a run bulkDNA / scDNA instance (extracts the reduced statistics).

Pulls coverage / alt / true CN (the perfect-caller stand-in: in synthetic recovery the CN is exact, exactly the two-pass conditioning) / GC / mappability / C-site / variant + het masks from the assay's surfaced ground truth, then defers to :func:estimate_dna. Modality and breadth are read off the assay. het_mask is loci whose true alt fraction is ~0.5.