Bio Core Domains

Sequence Motifs and Protein Domains

When to Use

• Building a PWM/PPM from a set of aligned TF binding sites and scanning a promoter, enhancer, or genome window for matches on both strands
• Computing per-position information content and drawing a sequence logo for a DNA or protein motif
• Converting a PROSITE pattern (e.g. N-{P}-[ST]-{P} for N-glycosylation) into a regex and scanning a protein sequence for hits
• Parsing HMMER hmmscan --domtblout output or Pfam/InterPro accessions to identify and annotate protein domains
• Visualizing or comparing domain architecture across a protein family (e.g. SH3-SH2-kinase in Src-family kinases)

Version Compatibility

• Python ≥3.9, NumPy ≥1.24, matplotlib ≥3.7 (for logos/heatmaps)
• HMMER ≥3.3 (hmmscan, hmmsearch) with a hmmpress-indexed Pfam-A.hmm (release ≥35)
• InterProScan ≥5.60 for combined Pfam/SMART/PROSITE/CDD scans
• No external motif library required for the core code below (pure NumPy/re); logomaker ≥0.8 is a drop-in for publication-quality logos

Prerequisites

• pip install numpy matplotlib
• For real domain scans: conda install -c bioconda hmmer plus a local Pfam-A.hmm database
• Familiarity with FASTA/protein sequences (bio-sequence-io-read-sequences) and basic regex syntax

Key Concepts

• Motif: short functional pattern (<20 aa/bp), does not fold independently (e.g. NLS, phosphorylation site, TF binding site)
• Domain: structurally independent unit (50–300 aa), folds on its own, recurs across proteins (e.g. SH2, kinase domain)
• PFM → PPM → PWM: PFM is raw counts; PPM = PFM/N with pseudocounts; PWM = log2(PPM/background) — only the PWM is a valid additive scoring matrix
• Profile HMM (Pfam): represents a domain family with match/insert/delete states; scores in bits vs a null model

Goal: build a PWM from aligned binding sites, score candidate windows, and scan both strands of a longer sequence.

Approach: count bases per position with a pseudocount (PFM→PPM), convert to log-odds against a uniform background (PWM), then slide a window of the motif's length across the sequence and its reverse complement, keeping windows scoring above a fraction of the maximum possible score.

import numpy as np
import re

BASES = ['A', 'C', 'G', 'T']


def build_ppm(sequences, pseudocount=0.1):
    """Position Probability Matrix from a list of equal-length aligned sequences."""
    n_pos = len(sequences[0])
    counts = np.full((4, n_pos), pseudocount)
    for seq in sequences:
        for pos, base in enumerate(seq.upper()):
            if base in BASES:
                counts[BASES.index(base), pos] += 1
    return counts / counts.sum(axis=0)


def ppm_to_pwm(ppm, background=None):
    """Convert a PPM to a log-odds PWM (bits) against a background composition."""
    if background is None:
        background = np.array([0.25, 0.25, 0.25, 0.25])
    bg = background[:, np.newaxis]
    return np.log2((ppm + 1e-10) / bg)


def information_content(ppm):
    """Per-position information content in bits. Max = 2 bits for DNA (log2(4))."""
    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


def score_sequence(pwm, sequence):
    """Sum of log-odds scores for a sequence exactly matching the PWM width."""
    return sum(pwm[BASES.index(b), i]
               for i, b in enumerate(sequence.upper()) if b in BASES)


def scan_sequence(pwm, sequence, threshold_pct=0.6):
    """Slide the PWM across `sequence`; return (position, subsequence, score) hits above threshold."""
    motif_len = pwm.shape[1]
    max_score = np.sum(np.max(pwm, axis=0))
    threshold = threshold_pct * 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):
    """Return the reverse complement of a DNA sequence."""
    comp = str.maketrans('ATGCatgc', 'TACGtacg')
    return seq.translate(comp)[::-1]


def scan_both_strands(pwm, sequence, threshold_pct=0.6):
    """Scan forward and reverse strands; report all hits in forward-strand coordinates."""
    motif_len = pwm.shape[1]
    fwd = [(p, s, sc, '+') for p, s, sc in scan_sequence(pwm, sequence, threshold_pct)]
    rev = scan_sequence(pwm, reverse_complement(sequence), threshold_pct)
    rev_fwd = [(len(sequence) - p - motif_len, s, sc, '-') for p, s, sc in rev]
    return sorted(fwd + rev_fwd, key=lambda x: -x[2])


# Example: CRP (cAMP receptor protein) binding sites from E. coli
crp_sites = [
    "AAATGTGATCTAGATCACATTT", "GAATGTGATCTATATCACATTT",
    "CAATGTGATCGAGATCACATAA", "TAATGTGATCTTAATCACATAT",
]
ppm = build_ppm(crp_sites, pseudocount=0.5)
pwm = ppm_to_pwm(ppm)
print("Consensus:", ''.join(BASES[i] for i in np.argmax(ppm, axis=0)))
print("Total IC (bits):", information_content(ppm).sum().round(2))

Goal: match short functional motifs (PROSITE patterns) inside a protein sequence.

Approach: convert the dash-separated PROSITE syntax (x, [..], {..}, (n), (n,m)) into an equivalent Python regex, then use re.finditer to report 1-based match positions.

def prosite_to_regex(pattern):
    """Convert a PROSITE pattern string (e.g. 'N-{P}-[ST]-{P}') to a Python regex."""
    pattern = pattern.strip('.')
    regex_parts = []
    for elem in pattern.split('-'):
        m = re.match(r'^(.+?)\((\d+)(?:,(\d+))?\)$', elem)
        core, low, high = (m.group(1), m.group(2), m.group(3)) if m else (elem, None, None)

        if core == 'x':
            r = '.'
        elif core.startswith('['):
            r = core
        elif core.startswith('{'):
            r = f'[^{core[1:-1]}]'
        else:
            r = core

        if low is not None:
            r += f'{{{low},{high}}}' if high else f'{{{low}}}'
        regex_parts.append(r)
    return ''.join(regex_parts)


def scan_prosite(sequence, prosite_pattern, pattern_name="pattern"):
    """Scan a protein sequence for a PROSITE pattern; returns 1-based (start, end, match) tuples."""
    regex = prosite_to_regex(prosite_pattern)
    return [(m.start() + 1, m.end(), m.group()) for m in re.finditer(regex, sequence)]


# Examples:
# N-glycosylation:      'N-{P}-[ST]-{P}'  -> 'N[^P][ST][^P]'
# PKC phosphorylation:  '[ST]-x-[RK]'     -> '[ST].[RK]'
# Zinc finger C2H2:     'C-x(2,4)-C-x(3)-[LIVMFYWC]-x(8)-H-x(3,5)-H'
protein = "MKVLLFAANISTHRGALVNVTPKSCNLTKVDYKNQTLLGSNLSECVFAIDNATASEKFLNYTRAR"
print(scan_prosite(protein, "N-{P}-[ST]-{P}", "N-glycosylation"))

Goal: turn a real hmmscan --domtblout file into annotated protein domains.

Approach: parse the fixed-width whitespace fields (domain name/acc, sequence/domain E-values and scores, alignment coordinates), filter on domain-level E-value, then map hits onto the protein for an architecture plot.

def parse_domtblout(text):
    """Parse HMMER --domtblout text into a list of domain-hit dicts, sorted by start position."""
    hits = []
    for line in text.strip().split("\n"):
        if line.startswith("#") or not line.strip():
            continue
        fields = line.split()
        if len(fields) >= 23:
            hits.append({
                "domain_name": fields[0],
                "domain_acc": fields[1],
                "query_name": fields[3],
                "e_value": float(fields[6]),        # sequence-level
                "score": float(fields[7]),
                "dom_e_value": float(fields[12]),    # domain-level (use this to filter)
                "dom_score": float(fields[13]),
                "ali_from": int(fields[17]),
                "ali_to": int(fields[18]),
            })
    return sorted(hits, key=lambda h: h["ali_from"])


def domains_above_threshold(hits, dom_evalue_max=1e-5):
    """Keep only domain hits passing a domain-level E-value cutoff (not the sequence-level one)."""
    return [h for h in hits if h["dom_e_value"] <= dom_evalue_max]

Domain Database Overview

# HMMER search against a local, hmmpress-indexed Pfam-A
hmmscan --domtblout results.domtbl Pfam-A.hmm protein.fasta

# InterProScan (Pfam + SMART + PROSITE + CDD in one pass)
interproscan.sh -i protein.fasta -f tsv -o results.tsv

Pitfalls

• Profile HMM E-values: hmmscan reports two E-values — sequence-level and domain-level; use the domain E-value (dom_e_value above) when annotating individual domain hits, not the sequence-level one
• Gathering threshold (GA): Pfam uses family-specific GA bit-score thresholds, not a universal cutoff; respect per-family thresholds rather than a fixed E-value everywhere
• Domain boundaries are approximate: Pfam HMM boundaries reflect consensus, not exact structural extents — check the alignment for precise boundaries
• PROSITE false positives: short degenerate patterns (e.g. N-x-[ST]) match frequently by chance; cross-validate with structural or conservation evidence
• Scan both strands: TF binding sites occur on either strand of dsDNA — always scan the reverse complement too, as scan_both_strands does
• PWM threshold choice: 50–70% of the maximum score is a reasonable starting point; tune against known sites in your organism's genome
• Pseudocounts matter: without pseudocounts, one missing base at a position gives −∞ log-odds; use ≥0.1 pseudocount per position (0.5 is a common default for small site sets)

See Also

• bio-sequence-manipulation-motif-search — simpler single-motif regex/consensus searches
• bio-chip-seq-motif-analysis — de novo motif discovery from ChIP-seq peaks (MEME-style)
• jaspar-database — fetching real PFMs/PWMs for known TFs instead of building them from scratch
• bio-genome-annotation-functional-annotation — genome-scale Pfam/InterPro annotation pipelines