Bio Analytical Validation

Version Compatibility

Reference examples tested with: numpy 1.26+, scipy 1.12+, statsmodels 0.14+

Before using code patterns, verify installed versions match. If versions differ:

• Python: pip show <package> then help(module.function) to check signatures

If code throws ImportError, AttributeError, or TypeError, introspect the installed package and adapt the example to match the actual API rather than retrying.

Analytical Validation and Detection Limits

"What is the real limit of detection of my ctDNA assay, and can I trust the number I am about to report?" -> Quantify the Poisson sampling ceiling, the error-suppression floor, and the LoB/LoD/LoQ that together define a defensible sensitivity claim.

• Python: scipy.stats.poisson for detection-probability math, scipy.stats.norm for CLSI LoB/LoD, statsmodels Probit/Logit for a dilution-series LoD95 fit.

The Single Most Important Modern Insight -- LoD Is Set by Genome Equivalents Sampled and Error Suppression, NOT by Sequencing Depth or the Caller

A ctDNA assay is a molecule-counting experiment at the Poisson edge. The mutant signal is a fixed, tiny number of physical template molecules in the tube, and the limit of detection is governed by two ceilings: how many genome equivalents were sampled (Poisson), and how low the background error floor was driven (error suppression). 1 ng of human DNA is ~330 haploid genome equivalents; the expected mutant-molecule count is lambda = input_GE x VAF. A 0.1% variant on 1,000 GE (~3.0 ng) has lambda = 1, so e^-1 ~= 37% of the time the mutant template was never in the tube and a perfect sequencer detects nothing. Past the point where every input molecule has been read once (sampling saturation, visible as a deduplication plateau in UMI families), additional read depth re-sequences the same physical molecules and adds zero information. Reporting an LoD as a bare VAF -- with no input mass, no unique-molecule (consensus) depth, no replicate detection rate -- is reporting an undefined quantity.

The second ceiling is the per-base background error rate, which sets the VAF floor independently: a 0.1% variant cannot be distinguished from noise if the assay manufactures that base at 0.1%. Error suppression is a ladder (raw NGS ~1e-3 -> single-strand UMI consensus ~1e-4/1e-5 -> duplex <1e-7), and single-strand consensus does NOT remove template-resident damage (C->T deamination, G->T 8-oxoG) because every PCR copy of that strand inherits the lesion -- only duplex strand-concordance catches it. The achieved LoD is the worse of the two ceilings: error dominates above ~0.1% VAF for tumor-naive single-locus calling, sampling dominates below it. The escape hatch is integration -- a bespoke panel summing mutant molecules across 16-50 loci against summed background reaches single-ppm even though each locus alone is ~1e-3 to 1e-4 (per-locus vs panel-integrated LoD).

Methods Landscape

Decision Tree by Scenario

Genome-Equivalent and Poisson Detection Calculator

Goal: Convert an input mass and target VAF into an expected mutant-molecule count and a detection probability, so a sensitivity claim is anchored to molecules rather than to a VAF alone.

Approach: Convert ng to haploid genome equivalents (~330/ng), set lambda = input_GE x VAF, and read the detection probability as a Poisson tail P(X >= k) = 1 - cdf(k-1, lambda); invert for the minimum GE that puts lambda at the >=3 sampling-detection threshold.

import numpy as np
from scipy.stats import poisson

GE_PER_NG = 330  # haploid ~3.3 pg -> strict 1 ng / 3.3 pg = 303; 330 is the common diploid-6.6 pg/rounding convention

def genome_equivalents(input_ng):
    return input_ng * GE_PER_NG

def detection_probability(input_ng, vaf, min_mutant_molecules=1):
    '''P(at least min_mutant_molecules present) under Poisson(lambda = GE * VAF).'''
    lam = genome_equivalents(input_ng) * vaf
    return float(poisson.sf(min_mutant_molecules - 1, lam))

def ge_for_sampling_detection(vaf, target_lambda=3.0):
    '''GE needed so lambda >= 3 -> ~95% chance the mutant molecule is present at all.'''
    return target_lambda / vaf

# A 0.1% variant on 3.0 ng (~990 GE) has lambda ~= 1 -> ~63% detected, ~37% missed by sampling alone.
detection_probability(3.0, 0.001)          # ~0.63
ge_for_sampling_detection(1e-4)            # 30000 GE (~91 ng) for a single 0.01% variant

LoB and LoD95 from a Dilution Series (CLSI EP17 style)

Goal: Estimate the VAF at which the assay detects 95% of the time, from a contrived dilution series, and anchor it to the blank-derived false-positive floor.

Approach: Compute LoB from blank replicates (mean + 1.645*SD, one-sided 95th pct), then fit a probit GLM of binary detection on log10(VAF) (CLSI EP17 fits on log concentration) across the dilution series and invert it for the 95% detection point. The series must bracket the 0.95 crossing — all-detected upper levels cause near-complete separation and an unstable slope.

import numpy as np
import statsmodels.api as sm
from scipy.stats import norm

def limit_of_blank(blank_signals):
    '''LoB = mean + 1.645*SD; one-sided 95th percentile of analyte-free blanks.'''
    blank_signals = np.asarray(blank_signals, dtype=float)
    return blank_signals.mean() + 1.645 * blank_signals.std(ddof=1)

def lod95_probit(vaf_levels, detected):
    '''Probit fit of detection (0/1) on log10(VAF); returns the VAF where P(detect) = 0.95.'''
    log_vaf = np.log10(np.asarray(vaf_levels, dtype=float))
    y = np.asarray(detected, dtype=float)
    X = sm.add_constant(log_vaf)
    fit = sm.GLM(y, X, family=sm.families.Binomial(link=sm.families.links.Probit())).fit()
    intercept, slope = fit.params
    return 10 ** ((norm.ppf(0.95) - intercept) / slope)

Per-Locus to Panel-Integrated LoD

Goal: Combine independent per-locus detection probabilities into the panel-level detection probability that a bespoke MRD assay actually achieves, and find the integrated LoD.

Approach: Treat each tracked locus as an independent Poisson sampler at the same tumor VAF; a panel positive call requires at least k loci detected, so the panel detection probability is the binomial-tail over the per-locus probabilities -- this is why summing 16-50 loci reaches ppm.

import numpy as np
from scipy.stats import poisson, binom

def panel_detection_probability(input_ng, vaf, n_loci, min_loci_positive=2):
    '''P(>= min_loci_positive of n_loci detected); >=2-of-N is the Signatera-style positivity rule.'''
    per_locus = float(poisson.sf(0, input_ng * 330 * vaf))
    return float(binom.sf(min_loci_positive - 1, n_loci, per_locus))

def panel_integrated_lod95(input_ng, n_loci, min_loci_positive=2, grid=None):
    '''Lowest VAF on a log grid where the >=k-of-N panel call hits 95%.'''
    grid = np.logspace(-6, -2, 400) if grid is None else np.asarray(grid)
    probs = [panel_detection_probability(input_ng, v, n_loci, min_loci_positive) for v in grid]
    hits = grid[np.asarray(probs) >= 0.95]
    return float(hits.min()) if hits.size else float('nan')