Statistics, seeds, refinement

The building blocks fit is made of. Stable enough to call directly, but the round loop is the supported surface.

pystreme.statistics

pystreme.statistics

Vectorized enrichment statistics.

See DESIGNDOC.md §2.5. Phase 2. The vectorized optimizer is validated (L2) in tests/test_statistics.py against a brute-force per-threshold scipy.stats.fisher_exact loop -- that reference loop is written as a test fixture, not thrown away, since it's the oracle for this file for good.

ThresholdResult dataclass

Result of optimal_threshold: the best cut point and its significance.

log_pvalue is a natural-log p-value (not log10), computed via scipy's logsf rather than log(sf(...)) so it stays finite deep into underflow territory instead of collapsing to -inf -- the whole point of working in log-p space here (§2.5's "single most likely way to end up slower" concern is about speed, but STREME-style refinement also needs p-values sharp enough to rank candidates once sf itself would return machine zero).

table property

table: tuple[int, int, int, int]

2x2 contingency table as (primary_above, primary_below, control_above, control_below).

CentralEnrichment dataclass

Result of central_enrichment: log_pvalue is Bonferroni-adjusted (natural log) over the candidate central-window half-widths tried; half_width is the winning window (sites with |offset| <= half_width), n_central how many of the n_sites fell inside it, and expected how many would have under a uniform positional null.

optimal_threshold

optimal_threshold(scores_primary, scores_control) -> ThresholdResult

Sort-once, all-cut-points-at-once threshold optimization (§2.5).

scores_primary/scores_control are 1-D best-site scores per sequence (e.g. ScanResult.scores, one motif's row) -- -inf entries (no valid window in that sequence, §2.2) are kept as real population members that simply never clear any threshold, not dropped.

Candidate thresholds are every distinct finite score value observed in either set: any other real cut point would partition the data exactly like the nearest such value, so restricting to observed values loses nothing and keeps the search over one array instead of a continuum. For each candidate t, "above" means score >= t.

Fisher's exact test at every candidate (hypergeometric tail, with the bounded search of _fisher_logsf_min_candidates), whatever the two set sizes -- fixed-width windows make every sequence the same length, which is the condition STREME's Fisher branch actually needs (§8.2; see _log_enrichment_pvalues for why the old size-based Binomial branch is gone). One-sided for enrichment (more primary sequences above threshold than expected), matching STREME/SEA's convention.

Returns the ThresholdResult for whichever candidate threshold gives the smallest (most significant) p-value.

enrichment_at_threshold

enrichment_at_threshold(scores_primary, scores_control, threshold: float) -> ThresholdResult

Enrichment test at one fixed score cut (score >= threshold).

The hold-out half of STREME's significance story: optimal_threshold picks the cut on the training sequences (and so its p-value is optimistic -- it was selected), then this re-tests that same cut on sequences the search never saw, giving Motif.holdout_logp. Same test conventions as optimal_threshold (Fisher's exact test whatever the two set sizes, -inf scores never count as "above").

central_enrichment

central_enrichment(offsets, n_windows: int) -> CentralEnrichment | None

CentriMo-style central enrichment (§2.7) of best-site positions.

offsets are site centers relative to the sequence center, in bases (as Motif.positions reports them -- half-integers when the motif width and sequence length have different parity); n_windows is the number of possible window starts per sequence (L - w + 1), which fixes the null: a site is equally likely at any of them, so the chance it lands within h of center is (number of window starts whose center is within h) / n_windows.

For every candidate half-width h (each distinct |offset| a window center can take, out to the sequence edge) the number of sites within h is tested against that null with a one-sided Binomial test, and the best h is reported with its log p-value Bonferroni-multiplied by the number of candidates tried -- the same max-over-widths-then-correct shape CentriMo uses. Returns None when there are no sites.

hypergeom_logsf

hypergeom_logsf(k, M, n, N, n_terms: int = 512) -> np.ndarray

log P(X > k) for X ~ Hypergeometric(M, n, N), vectorized and fast.

scipy.stats.hypergeom.logsf is a Python loop that logsumexps the whole upper tail per element (this was 95% of fit's runtime, and the seed ranking hands it every candidate word at once -- millions on a real peak set). hypergeom.sf is a C ufunc but underflows to 0 for the strongly enriched words that matter most, so: log(sf) wherever sf is representable, and for the rest the tail series summed directly in log space -- pmf(k+1) (vectorized logpmf) times the geometric-like sum of successive pmf ratios, which by then (far above the mode, or sf wouldn't have underflowed) decay fast enough that n_terms terms are exact to double precision.

binom_logsf

binom_logsf(k, n, p, n_terms: int = 512) -> np.ndarray

log P(X > k) for X ~ Binomial(n, p): log(binom.sf) (a C ufunc) with the same log-space tail-series fallback as hypergeom_logsf where sf underflows (scipy's binom.logsf is just log(sf), so it returns -inf there).

pystreme.seeds

pystreme.seeds

Seed finding: k-mer counting and ranking.

See DESIGNDOC.md §2.1 (direct counting replaces the suffix tree for this job) and §5 (memory: 4^w count arrays, sort-based fallback above w=12). Phase 3.

Words are counted pooled with their reverse complement -- canonicalized to whichever of {word, revcomp(word)} sorts first as a packed integer -- since the scanner searches both strands (§2.4): a word and its RC are the same candidate motif, not two different ones splitting the same signal.

That holds only while the scan really is double-stranded. Under revcomp=False the two orientations are not the same candidate, and emitting the lexicographically smaller one hands refinement a seed the forward-only scanner can never match, so words are then counted in the orientation they were observed in.

Seed dataclass

top_seeds

top_seeds(primary_store: SequenceStore, control_store: SequenceStore, widths, n_per_width: int = 4, n_exact: int = 2000, revcomp: bool = True) -> dict[int, list[Seed]]

Rank candidate seed words per width by primary-vs-control enrichment.

One vectorized counting pass per width (§2.1): direct dense bincount into a 4**w-sized array when that fits comfortably (w <= 12, the design doc's own cutoff -- 4**12 is ~16.7M entries), else a sort/unique -based count over just the words actually observed (packed ids fit an int64 up to w=31, far past the default 6-15 width range, so nothing changes about how words are represented -- only how they're counted).

revcomp must match the scan the seeds are refined under: with False, words are counted in their observed orientation rather than pooled with their reverse complements (see the module docstring).

Enrichment is a single vectorized statistics.hypergeom_logsf call across every candidate word at once (the same one-sided "pool of primary+control occurrences, how many landed in primary" test statistics.py uses for score thresholds, applied here to raw word-occurrence counts instead of a scan-score cutoff).

pystreme.refine

pystreme.refine

PWM refinement: the iterative loop that dominates STREME's runtime.

See DESIGNDOC.md §2.1 (best-site lookup as batched convolution -- the reason this is fast here) and §3 (Refiner). Phase 3.

Each iteration: scan primary+control with the current PWM, find the enrichment-optimal score threshold (§2.5), then re-estimate the PWM from the primary sites clearing that threshold (their actual best-matching window, on whichever strand won). This is the same NREF-seeds x NITER-iterations structure §2.1 costs out for STREME -- a handful of seeds per width, ~20 iterations each, each iteration one batched GPU scan rather than a suffix tree walk.

RefinedMotif dataclass

refine

refine(seed_words: list[str], primary_store: SequenceStore, control_store: SequenceStore, n_iter: int = 20, revcomp: bool = True, kernel: str = 'auto') -> RefinedMotif

Refine a batch of same-width seed words; return the single best result.

Both stores must already have a fitted background of the same order (§2.3's "fit on the control set" -- see MotifDiscovery's callers of this function for how that's set up).

seed_to_pwm

seed_to_pwm(word: str, sharpness: float = 0.8) -> np.ndarray

A discrete seed word -> a (4, w) PWM: sharpness probability at the observed base, the remainder split evenly across the other three -- concentrated enough to scan meaningfully, soft enough that refinement (rather than the seed itself) does the real work.

trim_flanks

trim_flanks(pwm: ndarray, min_ic: float, min_width: int = 1) -> tuple[np.ndarray, int, int]

Drop leading/trailing columns with information below min_ic bits, never going below min_width columns. Returns (trimmed_pwm, n_left, n_right) -- the counts trimmed from each end, so a site position for the original PWM can be shifted by n_left to point at the trimmed one.

A refined PWM one or two columns wider than the real site keeps its extra columns near-flat (they were estimated from whatever base happens to flank each site), and the enrichment objective is blind to them -- see MotifDiscovery._discover_on. Trimming them is how the reported width ends up meaning something.

extend_flanks

extend_flanks(pwm: ndarray, threshold: float, store: SequenceStore, min_ic: float, max_width: int, revcomp: bool = True, kernel: str = 'auto', control_store: SequenceStore | None = None) -> np.ndarray

Grow a refined PWM outward while the column just beyond either edge is informative: re-estimate a (4, w+2) matrix from the sequences whose best site clears threshold (their windows widened by one base on each side, in motif orientation), keep an edge column if its information is at least min_ic bits, and repeat until neither edge qualifies or max_width is reached.

The mirror image of trim_flanks: a seed word that sits one base off the real site refines into a motif missing a column at one end (the fixed-width refinement can't shift its window), and the extra column at the other end is flat. Trimming removes the flat one, this puts the missing one back. Sites too close to a sequence edge to widen are simply left out of the re-estimate. The caller re-scores the result (its threshold changes with width).

threshold is only the starting point: score scale moves with width, so a threshold fitted to the pre-trim PWM is the wrong cut for the trimmed one (systematically too strict -- fewer columns, lower scores) and for each widened one in turn. Given control_store, the enrichment-optimal threshold is refitted for the current PWM at every step, so sites are selected by a cut that belongs to the matrix doing the selecting. Without it the passed threshold is used throughout, as before.

column_information

column_information(pwm: ndarray) -> np.ndarray

Per-column information content in bits, relative to a uniform base distribution: sum_b p log2(4 p), 0 for a flat column, 2 for a fixed base.