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).

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.

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.

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.

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.