Bio Core Motif Discovery

Motif Discovery

When to Use

• Building a PWM from a set of aligned transcription-factor binding sites (ChIP-seq summits, SELEX reads, JASPAR sites)
• Scoring or genome-wide scanning a sequence for matches to a known motif
• Computing per-position information content, KDIC, or an IUPAC consensus for a motif
• Choosing a statistically defensible score threshold for calling motif hits
• Testing whether a motif is enriched in a foreground set (e.g. peaks) vs background, with multiple-testing correction
• Matching a de novo discovered motif against JASPAR/HOCOMOCO (TomTom-style)

Version Compatibility

NumPy >= 1.24, SciPy >= 1.10, statsmodels >= 0.14, pandas >= 2.0, Python >= 3.10. Concepts match MEME Suite 5.x (STREME/TomTom) and JASPAR 2024 / HOCOMOCO v12 motif databases.

Prerequisites

pip install numpy scipy pandas statsmodels matplotlib. Assumes basic familiarity with DNA sequences (bio-sequence-manipulation-motif-search) and, for genome-scale scanning, BED coordinates (bio-genome-intervals-bed-file-basics).

Matrix Representations

Information content: IC[i] = 2 + sum(p[i,k] * log2(p[i,k])). IC=0 bits -> degenerate position; IC=2 bits -> fully conserved.

KDIC: mean IC / 2, normalized to [0,1]. A simplified teaching proxy (assumes uniform background) for ranking motifs by specificity — not the true KL-divergence KDIC of Kulakovskiy et al.

Goal: turn a set of aligned binding sites into a scoring matrix and quantify how informative each position is.

Approach: count bases per column (PFM) -> add pseudocounts and normalize (PPM) -> take log-odds vs background (PWM). Never score with a PPM directly (rare bases at rare positions would swamp the log line); only PWM log-odds are additive across the true probability model.

import numpy as np

ALPHABET = "ACGT"
BASE_TO_IDX = {b: i for i, b in enumerate(ALPHABET)}

def make_ppm(seqs, pseudocount=0.1):
    """Aligned equal-length sequences -> Position Probability Matrix (L x 4)."""
    L = len(seqs[0])
    counts = np.zeros((L, 4))
    for s in seqs:
        for i, b in enumerate(s.upper()):
            if b in BASE_TO_IDX:
                counts[i, BASE_TO_IDX[b]] += 1
    return (counts + pseudocount) / (len(seqs) + 4 * pseudocount)

def make_pwm(ppm, bg=0.25):
    """PPM -> log2-odds PWM against a uniform (or per-base array) background."""
    return np.log2(ppm / bg)

def information_content(ppm):
    """Per-position IC in bits, clipped to [0, 2] for numerical safety."""
    ic = 2 + np.sum(ppm * np.log2(np.clip(ppm, 1e-9, 1)), axis=1)
    return np.clip(ic, 0, 2)

def kdic(ppm):
    """Mean IC normalized to [0, 1]; higher = more specific motif."""
    return information_content(ppm).mean() / 2

IUPAC Consensus

IUPAC_MAP = {
    frozenset("A"): "A", frozenset("C"): "C", frozenset("G"): "G", frozenset("T"): "T",
    frozenset("AG"): "R", frozenset("CT"): "Y", frozenset("GC"): "S", frozenset("AT"): "W",
    frozenset("GT"): "K", frozenset("AC"): "M", frozenset("CGT"): "B", frozenset("AGT"): "D",
    frozenset("ACT"): "H", frozenset("ACG"): "V", frozenset("ACGT"): "N",
}

def iupac_consensus(ppm, ic_threshold=1.0, freq_threshold=0.25):
    """Degenerate-code consensus: 'N' for low-IC positions, else the IUPAC code
    for the set of bases at/above freq_threshold."""
    ic = information_content(ppm)
    cons = []
    for row, ic_val in zip(ppm, ic):
        if ic_val < ic_threshold:
            cons.append("N")
        else:
            dominant = frozenset(ALPHABET[j] for j, p in enumerate(row) if p >= freq_threshold)
            cons.append(IUPAC_MAP.get(dominant, "N"))
    return "".join(cons)

Scoring, Scanning, and Threshold Selection

Goal: score a sequence against a PWM, scan both strands of a longer sequence, and pick a threshold with a known false-positive rate.

Approach: sum log-odds at each aligned position for scoring. For threshold calibration, enumerate the exact background score distribution when the motif is short (L <= 10, 4**L possibilities); otherwise estimate it by Monte Carlo sampling under the background model.

from itertools import product

def score_sequence(seq, pwm):
    """Sum of PWM log-odds along seq; non-ACGT bases and overhang are skipped."""
    return sum(pwm[i, BASE_TO_IDX[b]] for i, b in enumerate(seq.upper())
               if b in BASE_TO_IDX and i < len(pwm))

def scan_sequence(genome_seq, pwm, rc=True):
    """Scan sequence (and reverse complement if rc=True). Returns [(pos, score, strand)]."""
    L = len(pwm)
    rc_map = str.maketrans("ACGT", "TGCA")
    hits = []
    def _scan(seq, strand):
        for i in range(len(seq) - L + 1):
            hits.append((i, score_sequence(seq[i:i + L], pwm), strand))
    _scan(genome_seq, "+")
    if rc:
        _scan(genome_seq[::-1].translate(rc_map), "-")
    return hits

def score_threshold(pwm, pctile=99.99, mc_samples=50_000, seed=42):
    """Background score distribution -> threshold at the given percentile.
    Exact enumeration for L<=10, Monte Carlo otherwise."""
    L = len(pwm)
    rng = np.random.default_rng(seed)
    if L <= 10:
        scores = np.array([sum(pwm[i, b] for i, b in enumerate(bases))
                            for bases in product(range(4), repeat=L)])
    else:
        scores = np.array([score_sequence("".join(rng.choice(list(ALPHABET), size=L)), pwm)
                            for _ in range(mc_samples)])
    return np.percentile(scores, pctile)

Enrichment Testing (Foreground vs Background)

Goal: decide whether a motif is enriched in a set of peak/foreground sequences relative to background, correcting for testing many motifs.

Approach: count PWM hits above threshold in each group, run Fisher's exact test on the 2x2 contingency table, then apply Benjamini-Hochberg FDR when testing multiple motifs.

from scipy.stats import fisher_exact
from statsmodels.stats.multitest import multipletests

def count_hits(sequences, pwm, threshold):
    """Number of sequences with a PWM score >= threshold."""
    return sum(1 for s in sequences if score_sequence(s, pwm) >= threshold)

def motif_enrichment(peak_seqs, bg_seqs, pwm, threshold):
    """Fisher's exact test (one-sided, 'greater') for motif enrichment in peaks vs background.
    Returns (fold_enrichment, p_value)."""
    hits_peak = count_hits(peak_seqs, pwm, threshold)
    hits_bg = count_hits(bg_seqs, pwm, threshold)
    table = [[hits_peak, len(peak_seqs) - hits_peak],
             [hits_bg, len(bg_seqs) - hits_bg]]
    _, p_value = fisher_exact(table, alternative="greater")
    fold = (hits_peak / len(peak_seqs)) / max(hits_bg / len(bg_seqs), 1e-9)
    return fold, p_value

def correct_multiple_motifs(pvals):
    """Benjamini-Hochberg FDR correction across several motifs' p-values."""
    _, qvals, _, _ = multipletests(pvals, method="fdr_bh")
    return qvals

De novo discovery (MEME) treats binding-site position as hidden per sequence and alternates: E-step estimates P(site at each window | current PWM), M-step re-estimates the PWM from those weighted windows. It converges to a local optimum — run multiple seeds, or use STREME for large datasets. To identify a discovered motif, compare its PWM against JASPAR 2024 / HOCOMOCO v12 (e.g. via pyjaspar) using TomTom-style column correlation.

Pitfalls

• PFM -> PPM -> PWM pipeline: PFM = raw counts; PPM = frequencies (add pseudocount first); PWM = log2(PPM/background). Only PWM is suitable for scoring.
• Pseudocount choice: too small -> log(0) = -inf; too large -> drowns signal. Wasserman & Sandelin (2004) recommend a total pseudocount of sqrt(N) spread by background frequency; fixed 0.1-0.5 is common for small motifs.
• Background model: equal (0.25) is simplest; a dinucleotide Markov model is better for GC-rich genomes.
• MEME EM: converges to a local optimum. Run multiple times or use STREME for large datasets.
• Threshold selection: no universal threshold. Options: p-value from score distribution, fraction of max score (e.g. 80%), or empirical calibration from ChIP-seq. Always report the threshold and its false-positive rate.
• Coordinate systems: BED = 0-based half-open; VCF/GFF = 1-based inclusive — matters when converting scan hits back to genome coordinates.

See Also

• bio-sequence-manipulation-motif-search — regex/IUPAC pattern search in sequences
• bio-chip-seq-motif-analysis — applying motif scanning to ChIP-seq peak sets
• bio-atac-seq-motif-deviation — motif accessibility deviation scores (chromVAR-style)
• jaspar-database — fetching curated TF motifs (PFMs) for comparison