| name | bio-core-sequence-motifs |
| description | Build PFM/PPM/PWM matrices from aligned sites with NumPy, scan sequences on both strands, compute information content, plot sequence logos, convert PROSITE patterns to regex. Use for motif scanning, TF binding site scoring, or PROSITE-to-regex tasks. |
| tool_type | python |
| primary_tool | NumPy |
Sequence Motifs
When to Use
- Building a position weight matrix (PWM) from a set of aligned transcription-factor binding sites (e.g. from JASPAR, a ChIP-seq peak summit set, or SELEX data).
- Scanning a promoter or genomic region for matches to a known motif, on either strand.
- Quantifying motif conservation/specificity with information content (bits) or rendering a sequence logo.
- Converting a PROSITE protein motif pattern (e.g.
N-{P}-[ST]-{P}) into a Python regex for scanning protein sequences.
- Any task that says "PWM", "PFM", "consensus sequence", "binding site scan", "sequence logo", or "motif score".
Version Compatibility
- Python >= 3.10, NumPy >= 1.24, Matplotlib >= 3.7 (only needed for the logo/heatmap plot).
- No specialized bio package is required for the core PFM/PWM math;
Bio.motifs (Biopython >= 1.81) offers equivalent built-ins if you prefer not to hand-roll this.
Prerequisites
pip install numpy matplotlib- Familiarity with
bio-sequence-manipulation-reverse-complement (reverse-complement logic reused here) and basic string/regex handling (re module).
- Aligned input sequences (all the same length) — motif discovery to produce that alignment is covered by
bio-core-motif-discovery, not here.
PFM -> PPM -> PWM Pipeline
Goal: turn a set of aligned binding-site sequences into a scoring matrix.
Approach: count bases per position (PFM), normalize to frequencies with a pseudocount (PPM), then convert to log-odds against a background model (PWM).
import numpy as np
BASES = ["A", "C", "G", "T"]
def build_pfm(sequences: list[str]) -> np.ndarray:
"""Build a Position Frequency Matrix (4 x L raw counts) from aligned sequences."""
length = len(sequences[0])
pfm = np.zeros((4, length), dtype=int)
for seq in sequences:
for pos, base in enumerate(seq.upper()):
if base in BASES:
pfm[BASES.index(base), pos] += 1
return pfm
def pfm_to_ppm(pfm: np.ndarray, pseudocount: float = 0.1) -> np.ndarray:
"""Convert PFM to a Position Probability Matrix, adding a pseudocount to avoid log(0)."""
n_seq = pfm.sum(axis=0)[0]
return (pfm + pseudocount) / (n_seq + 4 * pseudocount)
def ppm_to_pwm(ppm: np.ndarray, background: float = 0.25) -> np.ndarray:
"""Convert PPM to a log-odds PWM relative to a uniform background."""
return np.log2(ppm / background)
binding_sites = ["ATGACTCA", "ATGACTCA", "ATGACTTA", "GTGACTCA", ]
pfm = build_pfm(binding_sites)
ppm = pfm_to_ppm(pfm)
pwm = ppm_to_pwm(ppm)
consensus = .join(BASES[i] i np.argmax(ppm, axis=))
()
()
Scoring and Both-Strand Scanning
Goal: find and rank motif matches in a longer sequence, forward and reverse strand.
Approach: slide the PWM window across the sequence, sum log-odds at each position, keep hits above a threshold; also scan the reverse complement and re-map coordinates back to forward-strand positions.
def score_sequence(pwm: np.ndarray, sequence: str) -> float:
"""Sum log-odds scores at each aligned position of a candidate subsequence."""
return sum(pwm[BASES.index(b), i] for i, b in enumerate(sequence.upper()) if b in BASES)
def scan_sequence(pwm: np.ndarray, sequence: str, threshold: float | None = None) -> list[tuple]:
"""Slide the PWM across `sequence`; return (pos, subseq, score) hits >= threshold, best first."""
motif_len = pwm.shape[1]
max_score = np.sum(np.max(pwm, axis=0))
if threshold is None:
threshold = 0.6 * max_score
hits = []
for i in range(len(sequence) - motif_len + 1):
subseq = sequence[i:i + motif_len]
s = score_sequence(pwm, subseq)
if s >= threshold:
hits.append((i, subseq, s))
return sorted(hits, key=lambda x: -x[2])
def reverse_complement(seq: str) -> str:
comp = .maketrans(, )
seq.translate(comp)[::-]
() -> []:
motif_len = pwm.shape[]
fwd = [(pos, seq, score, ) pos, seq, score scan_sequence(pwm, sequence, threshold)]
rev = scan_sequence(pwm, reverse_complement(sequence), threshold)
rev_fwd = [((sequence) - pos - motif_len, seq, score, ) pos, seq, score rev]
(fwd + rev_fwd, key= x: -x[])
target =
pos, subseq, score, strand scan_both_strands(pwm, target):
()
Information Content and Sequence Logo
Goal: quantify per-position conservation and visualize it.
Approach: IC(pos) = 2 - Shannon entropy of the base distribution (DNA, max 2 bits); a logo stacks each base's letter with height = frequency x IC.
import matplotlib.pyplot as plt
def information_content(ppm: np.ndarray) -> np.ndarray:
"""IC per position (bits). Max = 2 for DNA (fully conserved), 0 = random."""
ic = np.zeros(ppm.shape[1])
for pos in range(ppm.shape[1]):
p = ppm[:, pos]
entropy = -np.sum(p * np.log2(p + 1e-12))
ic[pos] = 2.0 - entropy
return ic
BASE_COLORS = {"A": "green", "C": "blue", "G": "orange", "T": "red"}
def plot_sequence_logo(ppm: np.ndarray, title: str = "Sequence Logo") -> None:
"""Plot a simple stacked-bar sequence logo; bar height = freq(base, pos) * IC(pos)."""
ic = information_content(ppm)
n_pos = ppm.shape[1]
fig, ax = plt.subplots(figsize=(max(5, n_pos * 0.7), 3))
for pos in range(n_pos):
order = np.argsort(ppm[:, pos])
y_bottom = 0.0
for idx in order:
height = ppm[idx, pos] * ic[pos]
if height > 0.01:
base = BASES[idx]
ax.bar(pos + , height, bottom=y_bottom, width=,
color=BASE_COLORS[base], edgecolor=)
y_bottom += height
ax.set_xlim(, n_pos + )
ax.set_ylim(, )
ax.set_xlabel()
ax.set_ylabel()
ax.set_title(title)
ax.set_xticks((, n_pos + ))
plt.tight_layout()
plt.show()
ic = information_content(ppm)
()
PROSITE Pattern -> Python Regex
Goal: turn a qualitative PROSITE-style motif into a usable regex for protein sequence scanning.
Approach: split on -, map x->., [..] stays a character class, {..} becomes a negated class [^..], and trailing (n)/(n,m) become quantifiers.
import re
def prosite_to_regex(pattern: str) -> str:
"""Convert a PROSITE pattern string to a Python regex string.
Syntax: N-{P}-[ST]-{P} -> N[^P][ST][^P]
x = any AA, [ABC] = one of, {P} = not P, (n) = repeat, (n,m) = range
"""
pattern = pattern.strip(".")
parts = []
for elem in pattern.split("-"):
m = re.match(r'^(.+?)\((\d+)(?:,(\d+))?\)$', elem)
core = m.group(1) if m else elem
low, high = (m.group(2), m.group(3)) if m else (None, None)
if core == "x":
r = "."
elif core.startswith("[") and core.endswith("]"):
r = core
elif core.startswith("{") and core.endswith("}"):
r = f"[^{core[1:-1]}]"
else:
r = core
if low:
r += f"{{{low},{high}}}" if high else f"{{{low}}}"
parts.append(r)
return "".join(parts)
prosite_to_regex() ==
Pitfalls
- PFM vs PPM vs PWM: only the PWM (log-odds) is suitable for scoring. Raw PPM values near zero cause extreme sensitivity. PFM counts are just for display.
- Pseudocounts are mandatory: a single missing base at any position gives log(0) = -infinity in the PWM. Use 0.1-0.5 (or a background-scaled fraction) as the pseudocount.
- Threshold selection is empirical: no universal correct value. Common approaches: 60-80% of max possible score, or calibrated from a ChIP-seq positive set (e.g. FPR at a chosen threshold).
- Always scan both strands: TF binding sites can be on either strand; single-strand scanning misses ~half of sites.
- PROSITE is qualitative, PWM is quantitative: PROSITE patterns are all-or-nothing; for degenerate sites, PWMs give far better sensitivity/specificity.
- Information content vs raw frequency: two positions can share a consensus letter but differ sharply in IC — a 60%/40% split and a 99%/1% split look similar as a consensus string but represent very different motif strength.
See Also
bio-core-motif-discovery — de novo motif discovery (MEME-like) to produce the aligned input sequences this skill scores.bio-core-domains — protein domain/InterPro annotation, complementary to PROSITE pattern matching.bio-chip-seq-motif-analysis — applying PWM scanning to ChIP-seq peak sets.bio-sequence-manipulation-motif-search — simpler literal/regex motif search without scoring matrices.
