bio-core-sequence-motifs

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  # 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:
    
    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  # 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 + , 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.