Enrichment and discovery¶
Testing motifs for enrichment, co-occurrence and positional enrichment in a set of sequences, and finding new motifs de novo.
enrich_motifs(
motifs: Motif | Iterable[Motif],
sequences: Sequence[str] | Mapping[str, str],
bkg_sequences: Sequence[str]
| Mapping[str, str]
| None = None,
pvalue: float = 0.0001,
qvalue: float = 0.1,
test: str = "seqs",
rc: bool = True,
shuffle_k: int = 2,
seed: int | None = None,
pseudocount: int = 1,
nthreads: int = 1,
) -> pl.DataFrame
Test motifs for enrichment in target sequences.
Target and background sequences of any alphabet are scanned with
scan_sequences(), keeping hits with
P-values below pvalue, and each motif gets a one-sided Fisher's exact
test with pseudocount added to every cell. With test="seqs" it
tests the number of sequences with a hit, and with test="sites" the
number of hits out of all the positions the motif could occupy.
Q-values are Benjamini-Hochberg adjusted across all motifs.
Without bkg_sequences, the background is the targets shuffled with
shuffle_sequences() using the
Euler method with k=shuffle_k and seed.
Returns a DataFrame of the motifs with q-values up to qvalue, in
order of P-value, with the columns motif, motif_i,
consensus, target_seq_n, target_seq_hits,
target_site_hits, bkg_seq_n, bkg_seq_hits, bkg_site_hits,
enrichment, log2_enrichment, pvalue, log10_pvalue,
qvalue and log10_qvalue. motif_i is 0-based, enrichment
is the ratio of hit rates, with a floor of half a hit on the background
rate, and the log columns stay finite when the P-values are too small
for a double.
Examples:
The G-box and W-box planted in the example targets are found enriched over the controls, while the TATA-box is not:
>>> import universalmotif as um
>>> motifs = [um.example_motifs()[i] for i in (0, 3, 5)]
>>> targets = um.example_sequences("targets")
>>> controls = um.example_sequences("controls")
>>> res = um.enrich_motifs(motifs, targets, controls, nthreads=1)
>>> res.select("motif", "target_seq_hits", "bkg_seq_hits").rows()
[('G-box', 21, 3), ('W-box', 30, 9)]
>>> res["log2_enrichment"].round(1).to_list()
[2.8, 1.7]
motif_coocc(
motifs: Motif | Iterable[Motif],
sequences: Sequence[str]
| Mapping[str, str]
| None = None,
hits: DataFrame | None = None,
n_sequences: int | None = None,
pvalue: float = 0.0001,
rc: bool = True,
no_overlaps: bool = False,
no_overlaps_by: str = "pvalue",
max_distance: int | None = None,
min_coocc: int = 1,
pseudocount: int = 0,
self_pairs: bool = False,
nthreads: int = 1,
) -> pl.DataFrame
Find motif pairs that share sequences more often than expected.
Every pair of motifs gets a one-sided Fisher's exact test on the 2x2
table of sequences with hits from one motif, the other, both or neither,
with pseudocount added to every cell. Pairs sharing fewer than
min_coocc sequences are not tested. Q-values are Benjamini-Hochberg
adjusted across the tested pairs.
Give either sequences or a hits table plus n_sequences.
Sequences are scanned with
scan_sequences() using pvalue and
rc, and no_overlaps removes overlapping hits of the same motif. A
hits table, which works for any alphabet, needs a 0-based
sequence_i column and either a 0-based motif_i column or a
motif column of names.
With max_distance, two descriptive columns that do not affect the
test are added: both_clustered, the number of shared sequences whose
nearest pair of hits by start position is within
max_distance, and median_distance, the median of those nearest
distances. self_pairs also pairs each motif with itself.
nthreads applies to the scan and to testing the pairs.
Returns a DataFrame with the columns motif_a, motif_b,
a_only, b_only, both, neither, odds_ratio,
pvalue and qvalue, sorted by q-value then P-value, with untested
pairs and their null P-values last.
Examples:
>>> import universalmotif as um
>>> motifs = [um.example_motifs()[i] for i in (3, 5)]
>>> targets = um.example_sequences("targets")
>>> res = um.motif_coocc(motifs, targets, nthreads=1)
>>> res.select("motif_a", "motif_b", "a_only", "b_only", "both", "neither").row(0)
('G-box', 'W-box', 13, 22, 8, 57)
>>> round(res["pvalue"][0], 3)
0.256
motif_peaks(
hits: DataFrame,
seq_length: int,
mode: str = "central",
seq_strand: Mapping[object, str] | None = None,
qvalue: float = 0.1,
min_window: int = 10,
max_window: int | None = None,
window_step: int = 10,
position_step: int = 5,
nthreads: int = 1,
) -> pl.DataFrame
Test whether motif hits cluster along sequences of equal length.
A CentriMo-style test (Bailey and Machanick 2012). Only each motif's
highest-scoring hit per sequence is used, at its centre (start + end) /
2. If centres were uniform along the sequences, the number falling in a
window of width w would be binomial with probability w /
seq_length. Windows from min_window to max_window in steps of
window_step are tested, centred on the middle of the sequences with
mode="central" or every position_step along them with
mode="local", and max_window defaults to half the sequence length.
The best window's P-value is Bonferroni-corrected for the number of windows
tested, and q-values are Benjamini-Hochberg adjusted across motifs.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
hits
|
DataFrame
|
A hit table such as from
|
required |
seq_strand
|
Mapping[object, str] | None
|
A mapping from the values in the sequence column to |
None
|
Returns:
| Type | Description |
|---|---|
DataFrame
|
A DataFrame of the motifs with q-values up to |
Examples:
W-box sites planted towards the middle of the example targets give a central peak:
>>> import universalmotif as um
>>> wbox = um.example_motifs()[5]
>>> hits = um.scan_sequences(wbox, um.example_sequences("targets"), nthreads=1)
>>> peaks = um.motif_peaks(hits, seq_length=1000, nthreads=1)
>>> peaks.select("motif", "nhits", "best_window", "best_centre", "hits_in").row(0)
('W-box', 30, 490, 499.0, 27)
>>> f"{peaks['pvalue'][0]:.1e}"
'1.3e-04'
motif_finder(
sequences: Sequence[str] | Mapping[str, str],
bkg_sequences: Sequence[str]
| Mapping[str, str]
| None = None,
min_width: int = 6,
max_width: int = 15,
nmotifs: int = 10,
hit_pvalue: float = 0.0001,
stop_pvalue: float = 0.001,
qvalue: float = 0.05,
dedup_overlap: float = 0.5,
rc: bool = True,
shuffle_k: int = 2,
seed: int | None = None,
pseudocount: int = 1,
nthreads: int = 1,
) -> pl.DataFrame
Discover motifs enriched in DNA or RNA sequences.
A minimal STREME-like search (Bailey 2021). At each width from
min_width to max_width, which is at most 30, k-mers that a
one-sided Fisher's exact test finds in more positive than background
sequences seed PPMs built from their Hamming-1 neighbours, which are
refined over two rescans. A motif whose hits, the sites scoring with a
P-value below hit_pvalue, are enriched with a P-value of at most
stop_pvalue is kept and its sites masked, up to nmotifs per width.
Motifs from different widths sharing more than dedup_overlap of their
sites are deduplicated, Benjamini-Hochberg q-values are computed,
low-information flanks are trimmed, and motifs with q-values above
qvalue are dropped.
Each motif is the best of many candidate words, so its P-value, from
Fisher's exact test on the sequences with hits, is multiplied by the
number of words tested in the search that found it, given in
candidates, and the q-values count every word tested at every width.
This keeps random sequences from giving motifs, though the refinement of
each motif from its word is not accounted for, so the P-values are still
somewhat optimistic. log10_pvalue and log10_qvalue stay finite
when the P-values are too small for a double.
Without bkg_sequences, the background is the positives shuffled
with shuffle_sequences() using the
Euler method with k=shuffle_k and seed. Parallelism is across
widths.
Returns a DataFrame in order of P-value, with the motifs in an object
column motif alongside name, consensus, rank, width,
nsites, seqs_pos, seqs_neg, sites_pos, sites_neg,
n_pos, n_neg, candidates, pvalue, log10_pvalue,
qvalue and log10_qvalue. The motifs are PPMs
named motif_1, motif_2 and so on, with their q-values as
E-values. seqs_pos and seqs_neg count sequences with a hit, and
sites_pos and sites_neg count hits.
Examples:
A site planted in each of 20 random sequences is found again, here as its reverse complement:
>>> import universalmotif as um
>>> site = um.create_motif("TTGACCAT", nsites=50)
>>> seqs = um.create_sequences(seqnum=20, seqlen=100, seed=1)
>>> seqs = um.implant_motifs(site, seqs, seed=1)
>>> res = um.motif_finder(seqs, min_width=6, max_width=10, seed=1)
>>> res.select("consensus", "seqs_pos", "seqs_neg").rows()
[('TGGTCAA', 20, 0)]