Bio.motifs

Core Principle: Bio.motifs models sequence binding motifs as count matrices, Position Weight Matrices (PWMs), and Position-Specific Scoring Matrices (PSSMs) across JASPAR, MEME, and TRANSFAC specifications.

Building a motif

The matrix progression

counts  →  PWM (normalised)  →  PSSM (log-odds vs background)

pwm = m.counts.normalize(pseudocounts=0.5)        # position weight matrix
pssm = pwm.log_odds()                              # position-specific scoring matrix
pssm = pwm.log_odds(background={"A": .3, "C": .2, "G": .2, "T": .3})

pssm.max; pssm.min                                 # best/worst possible score
pssm.mean(background); pssm.std(background)        # expected score distribution

Pseudocounts are not optional. Without them, a position where a base was never observed gets probability 0, so its log-odds is −∞, and any sequence with that base scores −∞ regardless of how well it matches elsewhere. With five instances you have not proven a base is impossible.

pseudocounts=0.5 (Jeffreys prior) or pseudocounts="sqrt" are conventional.

Scanning a sequence

seq = Seq("TACACTGCATTACAACCCTAATT")

for pos, score in pssm.search(seq, threshold=3.0):
    strand = "+" if pos >= 0 else "-"
    start = pos if pos >= 0 else len(seq) + pos
    print(f"{strand} strand, position {start}, score {score:.2f}")

scores = pssm.calculate(seq)          # score at every position, as a NumPy array

A negative position means the reverse strand, and the value is measured from the end of the sequence. This is easy to mishandle. Transcription factors bind both strands, so scanning both is correct: just be careful converting back to coordinates.

pssm.calculate(seq) returning a NumPy array is the efficient route when you want to plot a score profile or find all local maxima. See ndarray.

Choosing a threshold

dist = pssm.distribution(background=background, precision=10**4)

threshold = dist.threshold_fpr(0.001)       # 0.1% false positive rate
threshold = dist.threshold_fnr(0.1)          # 10% false negative rate
threshold = dist.threshold_balanced(1000)    # balance FPR and FNR
threshold = dist.threshold_patser()          # a common heuristic

Do not pick a threshold by eye. A PWM scan over a 3 Gb genome at an arbitrary cutoff produces millions of "hits", nearly all of them noise: motif sites are short and information-poor, so statistical significance is not biological significance. Real TF binding depends on chromatin accessibility, cofactors and cooperativity, none of which a PWM knows about.

The honest use of a PWM scan is: restrict to regions you have independent evidence for (ATAC-seq peaks, ChIP-seq peaks, promoters), then scan.

Reading motif databases

Supported formats: "jaspar", "pfm", "transfac", "meme", "minimal", "sites", "clusterbuster", "xms", "pfm-four-columns", "pfm-four-rows".

JASPAR is the usual source for curated TF motifs; MEME format is what the MEME suite emits.

Sequence logos

For local, controllable logos, use logomaker:

import logomaker

df = pd.DataFrame(m.counts).div(pd.DataFrame(m.counts).sum(axis=1), axis=0)
info = logomaker.transform_matrix(df, from_type="probability", to_type="information")

fig, ax = plt.subplots(figsize=(4, 2))
logomaker.Logo(info, ax=ax, color_scheme="classic")
ax.set(xlabel="Position", ylabel="Information (bits)")

The y-axis in bits is the informative one: total height at a position is 2 - H for DNA, so a fully constrained position is 2 bits and an unconstrained one is 0.

Information content

Total information content tells you how specific a motif is. A 6-bp motif with 12 bits occurs by chance roughly every 4,096 bp: i.e. about 750,000 times in the human genome. That number is the reality check every motif analysis needs.

RNA motifs

Biopython 1.85 added RNA support to Bio.motifs, so you can build motifs over ACGU for RNA-binding-protein sites and splice signals.

Bioinformatics examples

# build a splice donor motif from annotated introns
donors = [genome[start:start+9] for start in intron_starts]
m = motifs.create(donors)
print(m.degenerate_consensus)                     # ~ 'MAGGTRAGT'
pwm = m.counts.normalize(pseudocounts=0.5)

# score a variant's effect on the splice site
ref_score = pwm.log_odds().calculate(Seq(ref_context))
alt_score = pwm.log_odds().calculate(Seq(alt_context))
print(f"Δ score = {alt_score.max() - ref_score.max():.2f}")

# scan ATAC peaks (not the whole genome) for a TF motif
with open("MA0139.1.jaspar") as fh:
    ctcf = motifs.parse(fh, "jaspar")[0]
pssm = ctcf.counts.normalize(pseudocounts=0.5).log_odds()
thr = pssm.distribution(precision=10**4).threshold_fpr(1e-4)

hits = []
for peak in peaks.itertuples():
    seq = genome[peak.chrom].seq[peak.start:peak.end]
    for pos, score in pssm.search(seq, threshold=thr):
        hits.append({"chrom": peak.chrom, "pos": peak.start + abs(pos),
                     "score": score, "strand": "+" if pos >= 0 else "-"})
hits = pd.DataFrame(hits)

# motif score profile across a promoter
scores = pssm.calculate(promoter_seq)
fig, ax = plt.subplots(figsize=(8, 2))
ax.plot(scores, lw=0.8)
ax.axhline(thr, ls="--", c="crimson", label=f"FPR 1e-4")
ax.set(xlabel="Position in promoter", ylabel="PSSM score")

The variant-effect example is a genuinely useful pattern: the difference in best motif score between reference and alternate context is a cheap, interpretable predictor of a regulatory or splicing effect.

Common mistakes

• No pseudocounts, giving −∞ scores.
• Arbitrary thresholds instead of threshold_fpr.
• Genome-wide scans treated as binding predictions.
• Ignoring the reverse strand, or mishandling the negative position convention.
• Wrong background frequencies. Human genome GC is ~41%, not 50%; using uniform background biases GC-rich motifs.
• Confusing PWM and PSSM. PWM is probabilities, PSSM is log-odds.
• Comparing scores across motifs of different lengths. Longer motifs score higher; normalise or use p-values.
• Forgetting weblogo() needs network access: it fails silently offline in some setups.

See also

Seq · ndarray · Bio.SeqUtils · Applied - Sequence Composition Analysis · Linear Algebra with NumPy

Exercise

This exercise rebuilds the counts to PWM to PSSM progression by hand, because Biopython is not installed here and because the arithmetic is where the understanding lives. Add a pseudocount of 0.5 to every base at every position before normalising, use a uniform background of 0.25 for all four bases, and take base-2 logarithms.

Write scan(seq, instances), where instances are motif sequences of equal length over ACGT and seq is the sequence to scan. Return the consensus, the total information content, and the best-scoring window. The consensus takes the highest-probability base at each position, breaking ties in favour of the base that sorts first. Information content is the sum over all positions and bases of p * log2(p / 0.25). Score a window by summing log2(p / 0.25) for its bases, and skip any window containing a character other than A, C, G or T: a run of N must not score. Print consensus=<s>, then ic=<x> rounded to three decimals, then pos=<i> score=<x> with the 0-based start of the best window and its score rounded to three decimals, breaking ties in favour of the lower position. If no window can be scored, print no hit in place of that last line. The sequence arrives on the first line of input and the instances follow, one per line.