MotifDiscovery and Motif¶
The public entry point and the result object. See the usage guide for how the pieces fit together.
pystreme.discovery.MotifDiscovery ¶
Fixed-width, GPU-batched de novo motif discovery.
from_bed + scan are the Phase 1 scanner; enrichment the Phase 2
SEA-equivalent; discover the Phase 3 single-motif finder; fit the
Phase 4 round loop -- the motifs = disc.fit(idx) call from §3's API
sketch. The constructor arguments are the per-object defaults fit
and discover fall back to when not given their own.
from_bed
classmethod
¶
from_bed(peaks, genome, size: int, bg_order: int = 2, max_n_frac: float = 0.2, device: str = 'cpu', **kwargs) -> 'MotifDiscovery'
Build a store and immediately fit its background (§2.3, order 2 by
default). scan() requires a fitted background, and there's no
obviously-right implicit default for a self-contained discovery
object the way there might be for a bare scanner call -- see
pystreme.scanner.scan's docstring for why that stays explicit
there. Remaining keyword arguments (widths, control,
holdout_frac, revcomp, seed) are the constructor's defaults.
from_fasta
classmethod
¶
from_fasta(path, bg_order: int = 2, max_n_frac: float = 0.2, device: str = 'cpu', **kwargs) -> 'MotifDiscovery'
Equal-length sequences already in a FASTA file (e.g. exported by
another tool), no genome/BED needed. Same background fitting and
constructor defaults as from_bed.
from_sequences
classmethod
¶
from_sequences(seqs: list[str], bg_order: int = 2, device: str = 'cpu', **kwargs) -> 'MotifDiscovery'
Equal-length sequences in memory (a list of ACGTN strings).
scan ¶
scan(pwm, idx=None, revcomp: bool | None = None, kernel: str = 'auto')
Standalone PWM scanner (§"Phase 1" deliverable). See pystreme.scanner.scan.
pwm is a single (4, w) probability matrix, or a list of them for a
batched multi-motif scan. A single pwm returns squeezed (n_seqs,)
tensors instead of scanner.scan's (1, n_seqs); a list returns the
full ScanResult with the motif axis intact.
enrichment ¶
enrichment(pwm, control=None, idx=None, revcomp: bool | None = None, kernel: str = 'auto', n_per_seq: int = 1, bg_order: int | None = None, seed: int | None = None) -> ThresholdResult
SEA-equivalent enrichment test for one known motif (§"Phase 2"
deliverable): scan pwm against the primary sequences and against a
matched control set, then find the score threshold that most
significantly enriches primary over control (§2.5).
control is "dinuc_shuffle" (the constructor default), a
control.DinucShuffle/control.GCMatched/control.Explicit
instance, or anything control.Explicit accepts directly -- a
SequenceStore, list[str], or an index array into self.store
(see control.resolve_control). idx restricts which of
self.store's sequences count as "primary" (the control set is then
built to match that subset, not the whole store). bg_order
defaults to whatever from_bed already fit (self.store.bg_order).
discover ¶
discover(widths=None, control=None, idx=None, n_per_width: int = 4, n_iter: int = 20, revcomp: bool | None = None, kernel: str = 'auto', n_per_seq: int = 1, bg_order: int | None = None, seed: int | None = None, min_flank_ic: float = 0.3, width_tolerance: float = 2.0) -> Motif
Single-motif de novo discovery (§"Phase 3" deliverable): find seed
words per width (§2.1), refine each width's best seeds (§2.1's
batched-convolution refinement loop), and return the single best
motif across all widths (see _discover_on for how widths are
compared: min_flank_ic trims flat flank columns, width_tolerance
breaks near-ties by information content).
This is one round of fit without the hold-out split, erasing, or
stopping rules -- one motif, one significance number computed on
the same data it was found in (Motif.holdout_logp stays None;
see Motif's docstring). Positions and central enrichment are
filled in the same way fit does.
fit ¶
fit(primary_idx=None, n_motifs: int = 3, pvalue_thresh: float = 0.05, *, widths=None, control=None, holdout_frac: float | None = None, seed: int | None = None, n_per_width: int = 4, n_iter: int = 20, revcomp: bool | None = None, kernel: str = 'auto', n_per_seq: int = 1, bg_order: int | None = None, patience: int = 3, keep_insignificant: bool = False, min_holdout: int = 10, min_flank_ic: float = 0.3, width_tolerance: float = 2.0, verbose: bool = False) -> list[Motif]
The round loop (§"Phase 4" deliverable): motifs = disc.fit(idx).
- Resolve
primary_idx(all ofself.storeby default) and build the control set (control, seeded byseed); fit the background on the control (§2.3). - Split primary and control into training / hold-out subsets
(
holdout_frac, seeded). Hold-out is skipped, with a warning, when either side would have fewer thanmin_holdoutsequences. - Each round: seeds -> refinement -> best motif on the training
set (
discover's machinery); re-test its threshold on the hold-out set (Motif.holdout_logp); report sites/positions/ central enrichment on the full primary set; then erase the motif's sites from every working set and go again. - Stop after
n_motifssignificant motifs (holdout_logp <= log(pvalue_thresh), or the training p-value when there is no hold-out), or afterpatienceconsecutive insignificant ones, or when no seed words remain.
Insignificant motifs are dropped unless keep_insignificant (they
are still erased before continuing, as in STREME). Motifs come back
in discovery order, Motif.round numbered from 1. min_flank_ic
and width_tolerance control how the per-round width is settled
(see _discover_on).
Hold-out caveat (§5): the hold-out p-value is valid for one call.
Calling fit repeatedly over overlapping subsets and picking the
best result reintroduces the selection the split was meant to
remove -- vary seed per call and keep a never-touched set for a
final check if that's the workflow.
to_bed ¶
to_bed(motifs, path, passing_only: bool = True) -> int
Write the motifs' sites as BED6 (genomic coordinates, motif|peak
name, log-odds score, strand). See pystreme.results.sites_bed.
to_meme ¶
to_meme(motifs, path, background: ndarray | None = None)
Write motifs to MEME minimal format. See pystreme.meme_io.write_meme.
motifs is a list of Motif (as fit returns; named
<round>-<consensus>, with nsites and the hold-out p-value -- or
the training one when there is none -- in the E= slot), or the
{name: (4, w) array} / list-of-MemeMotif forms write_meme
accepts directly.
pystreme.discovery.Motif
dataclass
¶
A discovered motif: PWM, consensus, significance, sites, positions.
Two significance numbers, deliberately kept apart (§5's hold-out
warning): log_pvalue is the enrichment found on the training
sequences -- the same sequences the motif was searched for in, so it's
an optimistic, selected number, fine for ranking but not for reporting.
holdout_logp re-tests the training threshold on the held-out
sequences the search never saw; that's the number fit uses to decide
whether a motif counts, and the one to quote. It's None when there
was no hold-out set (discover() directly, or fit(holdout_frac=0),
or too few sequences to split) -- then log_pvalue is all there is.
sites covers every primary sequence (best score/position/strand,
with -inf where no valid window exists); n_sites counts those at or
above threshold. positions are the centers of those passing sites
relative to the sequence center, in bases (§2.7: fixed-width,
peak-centered windows make this free), and central /central_logp
is the CentriMo-style central-enrichment test over them.