Sequences, controls, scanning¶

`pystreme.sequence_store`¶

pystreme.sequence_store.SequenceStore ¶

Fixed-width (N, L) sequence universe with an order-k background model.

Construct via from_bed or from_fasta. See DESIGNDOC.md §2.2 for the extraction and edge-handling policy this follows.

from_bed `classmethod` ¶

from_bed(peaks, genome, size: int, max_n_frac: float = 0.2, device: str = 'cpu') -> 'SequenceStore'

Fixed-width extraction: center ± size/2, HOMER -size semantics.

peaks is a BED path (or an already-parsed list of Interval). genome is a FASTA path (or an open pyfaidx.Fasta) with a .fai index alongside it. Sequences whose window would run off a chromosome end, or reference a chromosome not in the genome, are dropped with a warning; so are sequences whose N-fraction exceeds max_n_frac. See DESIGNDOC.md §2.2.

from_fasta `classmethod` ¶

from_fasta(path, device: str = 'cpu', max_n_frac: float = 0.2) -> 'SequenceStore'

Sequences from a FASTA of already-equal-length records.

Unlike from_bed there is no window to re-center: every record is used as-is (and must be the same length as the first record), only the max_n_frac filter applies.

from_sequences `classmethod` ¶

from_sequences(seqs: list[str], device: str = 'cpu') -> 'SequenceStore'

Build directly from a list of equal-length ACGTN strings, no FASTA/BED.

For in-memory sequences that didn't come from a genome+peaks pair -- e.g. control.py's dinucleotide-shuffle and GC-matched pools, or tests. intervals get a synthetic "control"/seq_i name since there are no real genomic coordinates to record.

subset ¶

subset(idx) -> 'SequenceStore'

A new store over a subset of rows, keeping the fitted background (indexed the same way) if one exists rather than dropping it -- useful for e.g. MotifDiscovery.enrichment's idx= and control.Explicit's index-array form, where re-fitting on the subset would silently change what "background" means.

to_strings ¶

to_strings() -> list[str]

Decode every row back to an ACGTN string (masked/erased positions read as N).

to_fasta ¶

to_fasta(path) -> None

Write the store as FASTA, one record per row named after its interval (chrom:start-end plus the BED name) -- the round trip needed to hand the same sequences to an external tool (STREME/FIMO for the L1/L6 comparisons).

erase ¶

erase(rows, starts, width: int) -> int

Mask out the windows [starts[i], starts[i]+width) of rows[i], in place: the bases become N (code 4, mask False), so nothing built on this store afterwards -- k-mer counts, scan windows, background fits -- can see them. This is STREME's "erasing" of a found motif's sites before the next round, so the next motif isn't the same one again.

Invalidates the fitted background (bg_cumsum is set back to None; call fit_background again) and the one-hot cache. Returns the number of windows erased.

onehot ¶

onehot(dtype: dtype = torch.float16) -> torch.Tensor

(N, 4, L) one-hot encoding, cached per dtype. N positions are all-zero.

background_model ¶

background_model(order: int = 2, pseudocount: float = 1.0) -> 'BackgroundModel'

Estimate an order-order Markov background from this store's sequences and return it, without scoring anything.

Estimation and application are separate because the round loop applies one model to four or five stores per round (training and hold-out primary/control, plus the full primary set it reports from). Re-estimating per store made every round count the same control set that many times over -- and since the graded fallback needs a table per order 0..k, that redundancy multiplied. Fit once, apply_background many.

apply_background ¶

apply_background(model: 'BackgroundModel') -> None

Score this store's own codes/mask under an already-estimated model and set bg_cumsum (plus bg_order/bg_table/ bg_tables/bg_marginals/bg_model).

A position with fewer than model.order valid bases of preceding context -- near the start of a sequence, or after an N or an erased base -- backs off to the highest order its context does support, and to the order-0 marginal only when the base immediately before it is invalid. Any window touching an N is masked out at scan time regardless of what value ends up here, so N positions are given a value of 0 rather than something meaningful -- but the valid positions after an N are scored for real, which is why their context cannot be silently read as T.

fit_background ¶

fit_background(order: int = 2, source: 'SequenceStore | None' = None, pseudocount: float = 1.0)

Fit an order-order Markov background and set bg_cumsum.

source defaults to self (fit on the sequences being scored); DESIGNDOC.md §2.3 recommends fitting on the control set once one exists (Phase 2), by passing that store's SequenceStore here instead -- the tables are estimated from source's sequences, then applied to self's own codes/mask to produce self.bg_cumsum. source need not have the same number of sequences or the same length as self.

This is background_model followed by apply_background, kept as the one-liner for the common single-store case. Call the two halves directly when one model is applied to several stores -- see background_model for why that matters.

background_window_sum ¶

background_window_sum(start: int, width: int) -> torch.Tensor

Sum of background log-likelihood over [start, start+width), all sequences.

Requires fit_background to have been called first.

background_window_sums ¶

background_window_sums(width: int) -> torch.Tensor

Background log-likelihood sum for every window of width, all starts at once.

Shape (N, L-width+1). This is the vectorized form the scanner uses: one subtraction of two cumsum slices covers every window position at once, regardless of background order (DESIGNDOC.md §2.3).

pystreme.sequence_store.Interval `dataclass` ¶

pystreme.sequence_store.Exclusion `dataclass` ¶

`pystreme.control`¶

pystreme.control ¶

Control sequence construction: one interface, three modes.

See DESIGNDOC.md §2.6. Phase 2. All three modes expose the same generate(store, n_per_seq=1) -> SequenceStore method, so MotifDiscovery.enrichment (and later, the Phase 4 round loop) can treat control="dinuc_shuffle" / control="gc_matched" / control=<anything else> uniformly via resolve_control.

DinucShuffle ¶

Per-sequence Euler-path (Altschul-Erikson) k-mer-preserving shuffle.

Preserves the frequency of every kmer-length word in each sequence exactly (default kmer=2, i.e. dinucleotide frequencies) while destroying higher-order structure, and never moves an N -- each maximal run of non-N bases is shuffled independently, in place, so separator positions (§2.2's masked positions) land exactly where they started.

Built once per fit/enrichment call, not per refinement iteration, but linear in sequences x n_per_seq (~1 ms per 200 bp sequence in pure Python). engine="numba" (the default whenever numba imports) runs the same algorithm compiled, ~50x faster; engine="python" is the reference implementation and the test oracle.

Algorithm, per sequence run: build the de Bruijn-style multigraph whose nodes are (kmer-1)-mers and whose edges are the observed kmer-mers (for kmer=2 this is exactly Altschul & Erikson's original mononucleotide graph / dinucleotide edges). A uniformly random walk from every node to the run's final node, loop-erased as it goes (Wilson's algorithm), fixes one "last edge" per node; that's exactly the "spanning in-tree rooted at the last node" condition that guarantees a random ordering of the remaining edges never gets stuck before every edge is used. The run is guaranteed to have at least one valid Eulerian path -- the original sequence is one -- so this always terminates.

GCMatched ¶

Sample a background pool to match the primary set's GC histogram.

pool is a SequenceStore of candidate background regions (build one with SequenceStore.from_bed/from_fasta/from_sequences -- typically random genomic windows of the same width, blacklist-filtered and excluding the peaks themselves; that exclusion is the caller's responsibility, not checked here).

GC fraction is computed over unmasked (non-N) positions only. Both the primary set and the pool are binned into n_bins equal-width bins over [0, 1]; for each primary sequence, one pool sequence is drawn from the matching bin, without replacement while the bin still has unused pool sequences, falling back to sampling with replacement (and warning once) if a bin runs out.

The GC-rich-motif caveat from §2.6 applies: GC-matching removes exactly the compositional signal that makes a GC-box findable. Expect a GC-rich motif to weaken or drop out relative to dinucleotide shuffle -- that's the control doing its job, not a bug. Run both and compare.

Explicit ¶

User-supplied control sequences: a SequenceStore, list[str], or an index array selecting rows of a store (e.g. a held-out subset of the same primary population).

Index arrays are resolved against universe when generate is given one, else against store. The public API passes the whole MotifDiscovery.store as the universe, which is what its docstrings promise: enrichment(idx=[0, 1], control=[3]) means row 3 of the original store, not row 3 of the two-row primary subset (which would raise, or silently select some other row).

resolve_control ¶

resolve_control(control, store: SequenceStore, n_per_seq: int = 1, seed: int | None = None, universe: SequenceStore | None = None) -> SequenceStore

Dispatch the control= argument from the public API sketch (§3) to a concrete control SequenceStore, matching store's sequence count (times n_per_seq).

Accepts the string shortcuts ("dinuc_shuffle", "gc_matched"), an already-constructed DinucShuffle/GCMatched/Explicit instance, or anything Explicit itself accepts (a SequenceStore, list[str], or an index array) passed directly.

seed makes the randomized controls reproducible (L7): it seeds the "dinuc_shuffle" shortcut, a DinucShuffle instance constructed without its own seed, and GCMatched's sampling. An instance's own explicit seed always wins over this one.

universe is the store an Explicit index array indexes into (see Explicit); generated controls are always built to match store, i.e. the primary rows actually being tested. They differ exactly when the caller restricted the primary set with idx.

`pystreme.scanner`¶

pystreme.scanner ¶

PWM scanning: the batched-convolution scanner.

See DESIGNDOC.md §2.1 (best-site lookup as convolution), §2.3 (background subtraction via cumsum, split from the motif term), §2.4 (reverse complement as a second kernel), §3.1 (two scan kernels, benchmarked not assumed). Phase 1.

score(window at i) = sum_j log P_motif(x_{i+j} | j) -- the PWM kernel term - sum_j log P_bg(x_{i+j} | context) -- store.bg_cumsum

The two terms are independent (the background one doesn't depend on the motif at all), so they're computed by different mechanisms and just subtracted: a convolution/gather for the first, a cumsum difference for the second.

ScanResult `dataclass` ¶

scan ¶

scan(pwms: list[ndarray], store: SequenceStore, idx: Tensor | ndarray | None = None, revcomp: bool = True, kernel: str = 'auto') -> ScanResult

Batched PWM scan: best score/position/strand per (motif, sequence).

pwms is a list of (4, w) probability matrices (as produced by meme_io.read_meme's MemeMotif.pwm); motifs may have different widths. store must already have a background fit (store.fit_background(...)) -- there's no implicit default here, since silently picking one on your behalf is exactly the kind of surprise that makes a scan hard to trust. idx restricts the scan to a subset of sequences (row indices into store); kernel is "conv", "gather", or "auto" (device-based default, see _default_kernel).