# Bio Core Domains

> Build/scan PWMs (PPM, log-odds, info content, logos) for DNA motifs, convert PROSITE patterns to regex, parse HMMER domtblout/Pfam hits. Use for TF binding scans, promoter motifs, PROSITE matches, or Pfam/InterPro domain mapping.

- Skill: `pavel-kravchenko/bio-core-domains` (Agent Skill)
- Install (CLI): `npx skillmds@latest add pavel-kravchenko/bio-core-domains`
- Raw SKILL.md: https://api.skillmd.com/api/skills/pavel-kravchenko/bio-core-domains/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-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.

```python
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.

```python
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.

```python
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

| Database | Content | Best use |
|----------|---------|----------|
| Pfam | Profile HMMs, ~20k families | General domain annotation |
| InterPro | Integrates Pfam, SMART, CDD, PROSITE | First-pass comprehensive scan |
| SMART | Signaling domains | Kinases, receptors |
| CDD | NCBI-curated, integrates Pfam/SMART | Free with BLAST |
| PROSITE | Patterns and profiles | Active sites, short motifs |

```bash
# 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