# Bio Core Sequence Motifs > 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. - Skill: `pavel-kravchenko/bio-core-sequence-motifs` (Agent Skill) - Install (CLI): `npx skillmds@latest add pavel-kravchenko/bio-core-sequence-motifs` - Raw SKILL.md: https://api.skillmd.com/api/skills/pavel-kravchenko/bio-core-sequence-motifs/raw - Safety review: pending - Works with: Claude Code, Claude.ai, OpenAI Codex - Category: Coding & Dev Tools - Author: pavel-kravchenko (https://skillmd.com/u/pavel-kravchenko) - Updated: 2026-09-17 - Page: https://skillmd.com/skills/pavel-kravchenko/bio-core-sequence-motifs --- # 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). ```python 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", "ATGACTCG"] pfm = build_pfm(binding_sites) ppm = pfm_to_ppm(pfm) pwm = ppm_to_pwm(ppm) consensus = "".join(BASES[i] for i in np.argmax(ppm, axis=0)) print(f"Consensus: {consensus}") print(f"PWM score range: {pwm.min(axis=0).sum():.2f} to {pwm.max(axis=0).sum():.2f}") ``` ## 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. ```python 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 # 60% of max is a common default 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: """Return the reverse complement of a DNA sequence.""" comp = str.maketrans("ATGCatgc", "TACGtacg") return seq.translate(comp)[::-1] def scan_both_strands(pwm: np.ndarray, sequence: str, threshold: float | None = None) -> list[tuple]: """Scan forward and reverse strand; convert rev-strand hit positions to forward coords.""" motif_len = pwm.shape[1] fwd = [(pos, seq, score, "+") for pos, seq, score in scan_sequence(pwm, sequence, threshold)] rev = scan_sequence(pwm, reverse_complement(sequence), threshold) rev_fwd = [(len(sequence) - pos - motif_len, seq, score, "-") for pos, seq, score in rev] return sorted(fwd + rev_fwd, key=lambda x: -x[2]) target = "GCCTAGATGACTCAGGTTTCCCGTGACTCAATGCAATGACTTACCC" for pos, subseq, score, strand in scan_both_strands(pwm, target): print(f"{pos:4d} {strand:>2} {subseq} {score:7.2f}") ``` ## 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. ```python 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 # log2(4) - H 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]) # ascending, so tallest letter drawn last y_bottom = 0.0 for idx in order: height = ppm[idx, pos] * ic[pos] if height > 0.01: base = BASES[idx] ax.bar(pos + 1, height, bottom=y_bottom, width=0.8, color=BASE_COLORS[base], edgecolor="none") y_bottom += height ax.set_xlim(0.4, n_pos + 0.6) ax.set_ylim(0, 2.1) ax.set_xlabel("Position") ax.set_ylabel("Information content (bits)") ax.set_title(title) ax.set_xticks(range(1, n_pos + 1)) plt.tight_layout() plt.show() ic = information_content(ppm) print(f"IC per position: {ic.round(2)}; total: {ic.sum():.2f} bits (max {2 * len(ic)})") # CTCF motifs: ~15-16 total bits; TATA box: ~10-12 bits ``` ## 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. ```python 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) # Examples: # N-glycosylation: N-{P}-[ST]-{P} -> N[^P][ST][^P] # Kinase C site: [ST]-x-[RK] -> [ST].[RK] # Zinc finger C2H2: C-x(2,4)-C-x(3)-[LIVMFYWC]-x(8)-H-x(3,5)-H assert prosite_to_regex("N-{P}-[ST]-{P}") == "N[^P][ST][^P]" ``` ## 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.