Bio Motif Search

Version Compatibility

Reference examples tested with: BioPython 1.83+

Before using code patterns, verify installed versions match. If versions differ:

• Python: pip show <package> then help(module.function) to check signatures

If code throws ImportError, AttributeError, or TypeError, introspect the installed package and adapt the example to match the actual API rather than retrying.

Motif Search

"Search for a sequence motif or binding-site pattern" -> Scan sequences for a fixed motif, a degenerate IUPAC consensus, or a probabilistic PWM, on one or both strands, and locate transcription-factor binding sites, regulatory elements, or custom patterns.

• Python: Bio.SeqUtils.nt_search (IUPAC + overlaps), re (regex/lookahead), Bio.motifs (PWM/PSSM scoring + matrix file parsing)

The Governing Principle

Two silent failures dominate motif searching; both return a plausible-but-wrong answer with no error:

1. Overlapping matches are dropped. str.count, str.find, and re.findall consume the string left to right, so a motif that overlaps its own next occurrence is undercounted. Target AAGCGCGCGAA, motif GCGC: str.count returns 1, the true answer is 2 (starts at positions 2 and 4). Use a zero-width lookahead re.finditer(r'(?=(GCGC))', target) or Bio.SeqUtils.nt_search, both of which report overlaps.
2. A PSSM score is a likelihood in bits, not a probability. pssm.calculate returns log2-odds versus background. A "high-looking" threshold chosen by eye is arbitrary and non-reproducible; derive the threshold from the score distribution at a chosen false-positive rate. And a PSSM scans only the strand it is given, so scoring just the forward strand silently misses roughly half of real sites on double-stranded DNA.

Which Approach for Which Question

IUPAC Degenerate Motifs

A degenerate motif expands each ambiguity code to a regex character class:

B, D, H, V each exclude A, C, G, T respectively (the code preceding the one they drop is a mnemonic).

IUPAC_DNA = {'N': '[ACGT]', 'R': '[AG]', 'Y': '[CT]', 'W': '[AT]', 'S': '[GC]',
             'K': '[GT]', 'M': '[AC]', 'B': '[CGT]', 'D': '[AGT]', 'H': '[ACT]', 'V': '[ACG]'}

def iupac_to_regex(pattern):
    return ''.join(IUPAC_DNA.get(base, base) for base in pattern)

# 'GATNNTC' -> 'GAT[ACGT][ACGT]TC'

Bio.SeqUtils.nt_search expands IUPAC ambiguity in the query motif automatically and reports overlapping hits, so it is the shortest correct path for a degenerate consensus.

Overlapping Matches (the count trap)

Goal: Report every start position of a motif, including self-overlapping occurrences.

Approach: Use a zero-width lookahead so the regex engine never consumes the matched text; recover the match string from the inner capture group. For IUPAC motifs, prefer nt_search, which both expands ambiguity and reports overlaps.

Reference (BioPython 1.83+):

import re
from Bio.SeqUtils import nt_search

target = 'AAGCGCGCGAA'

starts = [match.start(1) for match in re.finditer(r'(?=(GCGC))', target)]  # [2, 4]
hits = [(match.start(1), match.group(1)) for match in re.finditer(r'(?=([AG]CG[CT]))', target)]

result = nt_search(target, 'GCGC')  # ['GCGC', 2, 4]
pattern, positions = result[0], result[1:]

nt_search returns a heterogeneous list: result[0] is the (expanded) pattern string and result[1:] are the 0-based start positions. When there are no hits it returns just [pattern] (length 1), so test len(result) > 1 before indexing rather than truthiness.

Bio.motifs PWM / PSSM Pipeline

Goal: Build a probabilistic model from a set of aligned binding sites and score a target sequence for graded matches.

Approach: Create a motif from instances or a matrix file, set pseudocounts and background, read the recomputed PSSM, then scan. The count matrix m.counts['A', 0] is indexed [base, position].

Reference (BioPython 1.83+):

from Bio import motifs
from Bio.Seq import Seq

m = motifs.create([Seq('TACAA'), Seq('TACGA'), Seq('TACTA'), Seq('TGCAA')])  # alphabet defaults to ACGT

m.pseudocounts = 0.5        # set BEFORE reading m.pssm (see trap below)
m.background = None         # None gives uniform 0.25; or pass a dict of base frequencies

pwm = m.pwm                 # normalized frequencies (property)
pssm = m.pssm               # log2-odds vs background (property; RECOMPUTED on each access)

m.consensus                 # most frequent base per column
m.degenerate_consensus      # IUPAC-degenerate consensus

m.counts.normalize(pseudocounts=0.5) returns a position weight matrix and pwm.log_odds() returns a PSSM; these are equivalent to reading m.pwm / m.pssm after setting m.pseudocounts.

The Pseudocount / -inf Trap (silent)

A column where some base has count 0 gives that base frequency 0 and a log-odds of negative infinity; any target carrying that base at that position then scores -inf and is unmatchable. This is common with short motifs or few instances. Setting m.pseudocounts (a flat 0.5, or sqrt(N) with N the number of instances; scalar or per-base dict) makes every cell finite by shrinking toward background.

Critically, m.pssm is recomputed from m.pseudocounts and m.background on every access. Set both BEFORE reading m.pssm (or pwm.log_odds()), or the matrix is silently wrong.