Bio Applied Regulatory Analysis

Promoter and Regulatory Sequence Analysis

When to Use

• Searching a promoter/upstream region for a TATA box, Inr, or other core promoter motif
• Calling CpG islands (or checking a region against the classic ≥200bp / GC≥50% / O:E≥0.6 rule)
• Building a PFM/PWM from a set of aligned transcription factor binding sites and scanning a sequence for hits
• Profiling GC content, dinucleotide frequency, or motif density relative to a known/candidate TSS
• Testing whether a TFBS is enriched in one gene set (e.g. stress-responsive) vs. another (e.g. housekeeping)

Version Compatibility

• Python ≥ 3.10, NumPy ≥ 1.24, pandas ≥ 2.0, SciPy ≥ 1.11 (for scipy.stats.mannwhitneyu)
• No genome-specific databases required — all examples work on plain strings; swap in real sequences via BioPython (Bio.SeqIO) as needed

Prerequisites

• pip install numpy pandas scipy biopython
• Familiarity with Python regex (re module) and basic sequence coordinates (TSS = +1, upstream = negative)
• Related skill: bio-sequence-manipulation-motif-search for generic motif scanning; bio-chip-seq-motif-analysis for ChIP-seq-derived motifs

Key Regulatory Elements

• TSS = position +1; upstream positions are negative
• TATA box: consensus TATAAA at ~−30; present in ~10–20% of human genes (TATA-less promoters use Inr, DPE instead)
• CpG islands: near ~70% of human gene promoters; classic criteria: length ≥200 bp, GC ≥50%, CpG observed/expected (O/E) ≥0.6

TATA Box Detection

Goal: find TATA-box-like motifs in a promoter sequence and report their position relative to the TSS. Approach: regex search for the exact consensus (TATAAA) or the relaxed IUPAC pattern (TATAWAW, W = A/T).

import re


def find_tata_boxes(sequence, strict=True):
    """Find TATA box motifs in a DNA sequence.

    strict=True:  exact consensus TATAAA
    strict=False: relaxed IUPAC pattern TATA[AT]A[AT] (TATAWAW)
    Returns a list of (start_position, matched_motif) tuples.
    """
    sequence = sequence.upper()
    pattern = 'TATAAA' if strict else 'TATA[AT]A[AT]'
    return [(m.start(), m.group()) for m in re.finditer(pattern, sequence)]


# Example: TSS at position 1000, TATA box inserted at -30
tss = 1000
promoter_seq = "N" * 970 + "TATAAAG" + "N" * 1023
for pos, motif in find_tata_boxes(promoter_seq):
    print(f"Position {pos} (TSS{pos - tss:+d}): {motif}")

CpG Island Detection

Goal: call CpG islands in a promoter and merge overlapping sliding-window hits into island coordinates. Approach: slide a window across the sequence, compute GC% and CpG O/E per window, keep windows passing both thresholds, then merge adjacent/overlapping passing windows.

def cpg_island_scanner(sequence, window=200, step=10, gc_thresh=0.5, oe_thresh=0.6):
    """Sliding-window CpG island detection.

    Returns a list of (start, end, gc_content, cpg_oe_ratio) for qualifying windows.
    Use step=1 for exact boundaries (slower); step=10-50 for a fast first pass.
    """
    sequence = sequence.upper()
    islands = []
    for i in range(0, len(sequence) - window + 1, step):
        win = sequence[i:i + window]
        n_c, n_g = win.count('C'), win.count('G')
        n_cpg = win.count('CG')
        gc = (n_c + n_g) / window
        oe = (n_cpg * window) / (n_c * n_g) if n_c > 0 and n_g > 0 else 0.0
        if gc >= gc_thresh and oe >= oe_thresh:
            islands.append((i, i + window, gc, oe))
    return islands


def merge_island_windows(windows, max_gap=0):
    """Merge overlapping/adjacent qualifying windows into contiguous CpG islands."""
    if not windows:
        return []
    sorted_wins = sorted(windows, key=lambda x: x[0])
    merged = [list(sorted_wins[0])]
    for start, end, gc, oe in sorted_wins[1:]:
        if start <= merged[-1][1] + max_gap:
            merged[-1][1] = max(merged[-1][1], end)
            merged[-1][2] = max(merged[-1][2], gc)
            merged[-1][3] = max(merged[-1][3], oe)
        else:
            merged.append([start, end, gc, oe])
    return merged


raw = cpg_island_scanner(promoter_seq, window=200, step=10)
for start, end, gc, oe in merge_island_windows(raw, max_gap=50):
    print(f"Island {start}-{end} ({end - start} bp)  GC={gc:.2f}  CpG O/E={oe:.2f}")

PWM/PFM Construction and Scanning

Goal: build a transcription-factor motif model from aligned binding sites and scan a sequence for matches (JASPAR/TRANSFAC-style TFBS scan). Approach: count bases per position (PFM), convert to log2-odds vs. background (PWM), then slide the motif across the target sequence and sum log-odds scores.

import numpy as np


def build_pfm(sites):
    """Build a position frequency matrix from aligned binding sites (same length)."""
    length = len(sites[0])
    pfm = {b: [0] * length for b in 'ACGT'}
    for site in sites:
        for i, base in enumerate(site.upper()):
            pfm[base][i] += 1
    return pfm


def pfm_to_pwm(pfm, pseudocount=0.5, bg=None):
    """Convert a PFM to a log2-odds PWM (adds a pseudocount to avoid -inf scores)."""
    bg = bg or {'A': 0.25, 'C': 0.25, 'G': 0.25, 'T': 0.25}
    n = sum(pfm[b][0] for b in 'ACGT')
    length = len(pfm['A'])
    pwm = {b: [0.0] * length for b in 'ACGT'}
    for b in 'ACGT':
        for i in range(length):
            freq = (pfm[b][i] + pseudocount) / (n + 4 * pseudocount)
            pwm[b][i] = np.log2(freq / bg[b])
    return pwm


def score_sequence(subseq, pwm):
    """Score one sequence window against a PWM: sum of per-position log-odds."""
    return sum(pwm[b][i] for i, b in enumerate(subseq.upper()) if b in pwm)


def scan_with_pwm(sequence, pwm, threshold=0.0):
    """Slide a PWM across a sequence; return (pos, subseq, score) above threshold."""
    sequence = sequence.upper()
    motif_len = len(pwm['A'])
    hits = []
    for i in range(len(sequence) - motif_len + 1):
        sub = sequence[i:i + motif_len]
        score = score_sequence(sub, pwm)
        if score >= threshold:
            hits.append((i, sub, score))
    return hits


tata_sites = ['TATAAAG', 'TATAAAT', 'TATAAAA', 'TATATAG', 'TATAAAC',
              'TATAAAT', 'TATAAAG', 'TATAAAT', 'TATAAAA', 'TATAAAG']
pwm = pfm_to_pwm(build_pfm(tata_sites))
for pos, sub, score in sorted(scan_with_pwm(promoter_seq, pwm, threshold=5.0),
                               key=lambda x: -x[2])[:5]:
    print(f"Position {pos} (TSS{pos - tss:+d}): {sub}  score={score:.2f}")

TSS Prediction Signals

TF Motif Databases

Pitfalls

• CpG underrepresentation: vertebrate bulk genome has CpG O/E ~0.2 due to methylation-driven mutation; islands (O/E ≥0.6) are real anomalies, not noise
• Coordinate systems: BED = 0-based half-open; VCF/GFF = 1-based inclusive — mixing causes off-by-one errors when reporting TSS-relative positions
• PWM threshold selection: score ≥80% of the max possible PWM score is a common heuristic; too low a threshold floods results with false positives
• Pseudocount matters: without pseudocounts, a single zero count in the PFM gives a −∞ PWM score for any sequence with that base at that position
• Multiple testing: scanning a genome-wide promoter set for TFBS enrichment requires FDR correction (e.g. Benjamini-Hochberg), not raw p-values
• Window step size: step=1 in cpg_island_scanner gives exact boundaries but is O(n); use a larger step for a fast first pass, then refine only around merged candidate islands

See Also

• bio-sequence-manipulation-motif-search — generic sequence motif searching
• bio-chip-seq-motif-analysis — motif discovery/enrichment from ChIP-seq peaks
• bio-chip-seq-peak-annotation — annotating peaks relative to TSS/promoters
• jaspar-database — fetching real PFMs/PWMs from JASPAR