Bio.motifs

In one line: sequence motifs as count matrices, PWMs and PSSMs — construct them, score sequences with them, and read the standard motif database formats.

Building a motif

from Bio import motifs
from Bio.Seq import Seq

instances = [Seq("TACAA"), Seq("TACGC"), Seq("TACAC"), Seq("TACCC"), Seq("AACCC")]
m = motifs.create(instances)

m.counts                # position frequency matrix (PFM)
m.consensus             # 'TACAC'
m.degenerate_consensus  # 'TACRC'  — IUPAC ambiguity codes
m.length
m.counts["A"]           # counts of A at each position

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

with open("jaspar.pfm") as fh:
    ms = motifs.parse(fh, "jaspar")

for m in ms:
    print(m.matrix_id, m.name, m.length)

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.

m.format("jaspar")           # write back out
m.format("transfac")

Sequence logos

m.weblogo("logo.png")        # calls the WebLogo web service — needs network

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

from math import log2

def position_ic(pwm, pos, background=0.25):
    return sum(p * log2(p / background)
               for base in "ACGT"
               if (p := pwm[base][pos]) > 0)

ic = [position_ic(pwm, i) for i in range(m.length)]
total_ic = sum(ic)

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