---
name: "Amino Acid Enrichment in Intrinsically Disordered Regions"
description: "Hypergeometric test with bootstrap CIs for amino acid enrichment in human protein disordered regions vs structured regions, correcting for composition baselines"
version: "1.0.0"
author: "Claw 🦞, David Austin"
tags: ["claw4s-2026", "proteomics", "intrinsic-disorder", "hypergeometric-test", "bootstrap", "uniprot"]
python_version: ">=3.8"
dependencies: []
---

# Amino Acid Enrichment in Intrinsically Disordered Regions

## Motivation

Intrinsically disordered regions (IDRs) in human proteins are known to be enriched for certain amino acids (Pro, Gly, Ser, Glu). But how strong is this enrichment compared to a null model that preserves overall amino acid composition? Raw frequency counts overstate enrichment because they ignore composition baselines — comparing IDR frequencies to a pooled baseline (which includes IDR residues) contaminates the reference with the signal. This skill corrects for that by comparing IDR composition against structured regions only, using a hypergeometric test with bootstrap confidence intervals and permutation validation.

## Prerequisites

- Python 3.8+ (standard library only — no pip install required)
- Internet access to reach `rest.uniprot.org`
- ~50 MB disk space for cached data

## Step 1: Create workspace

```bash
mkdir -p /tmp/claw4s_auto_uniprot-amino-acid-disorder-enrichment
```

**Expected output:** Directory created (no stdout).

## Step 2: Write analysis script

```bash
cat << 'SCRIPT_EOF' > /tmp/claw4s_auto_uniprot-amino-acid-disorder-enrichment/analysis.py
#!/usr/bin/env python3
"""
Amino acid enrichment in intrinsically disordered regions (IDRs) of human proteins.

Hypothesis: IDRs are enriched for Pro, Gly, Ser, Glu relative to a null model
that preserves overall amino acid composition.

Methodological hook: Raw frequency counts overstate enrichment because they
ignore protein-level composition baselines. We correct for this using
the hypergeometric test (which conditions on marginal totals) with
protein-level bootstrap CIs and within-protein permutation validation.

Data: 500 reviewed human proteins from UniProt with annotated disorder regions.
Statistical tests: Hypergeometric test, bootstrap CIs (2000 resamples),
permutation test (1000 permutations), sensitivity analysis, BH correction.
"""

import json
import math
import hashlib
import os
import sys
import random
import time
import urllib.request
import urllib.error
from collections import Counter

# === Constants ===
SEED = 42
N_BOOTSTRAP = 2000
N_PERMUTATIONS = 1000
AMINO_ACIDS = sorted("ACDEFGHIKLMNPQRSTVWY")
WORKSPACE = os.path.dirname(os.path.abspath(__file__)) or "."
CACHE_FILE = os.path.join(WORKSPACE, "uniprot_data.json")
RESULTS_FILE = os.path.join(WORKSPACE, "results.json")
REPORT_FILE = os.path.join(WORKSPACE, "report.md")

UNIPROT_BASE = "https://rest.uniprot.org/uniprotkb/search"
UNIPROT_QUERY = "(reviewed:true)%20AND%20(organism_id:9606)%20AND%20(ft_region:disordered)"
PAGE_SIZE = 100  # smaller pages to avoid chunked transfer issues


# === Data Download ===

def download_with_retry(url, max_retries=5, timeout=180):
    """Download URL content with exponential backoff retry and chunked read."""
    import http.client
    for attempt in range(max_retries):
        try:
            req = urllib.request.Request(url)
            req.add_header("User-Agent", "Python-urllib/3.8 (claw4s-research)")
            req.add_header("Accept", "application/json")
            with urllib.request.urlopen(req, timeout=timeout) as resp:
                # Read in chunks to handle large/chunked responses
                chunks = []
                while True:
                    try:
                        chunk = resp.read(65536)
                        if not chunk:
                            break
                        chunks.append(chunk)
                    except http.client.IncompleteRead as e:
                        chunks.append(e.partial)
                        break
                return b"".join(chunks)
        except (urllib.error.URLError, urllib.error.HTTPError, OSError,
                http.client.IncompleteRead) as e:
            print(f"  Attempt {attempt+1}/{max_retries} failed: {e}")
            if attempt < max_retries - 1:
                time.sleep(3 * (attempt + 1))
            else:
                raise RuntimeError(f"Failed to download after {max_retries} attempts: {e}")


def download_uniprot_data():
    """Download and cache UniProt data with pagination and SHA256 verification."""
    if os.path.exists(CACHE_FILE):
        print(f"  Using cached data: {CACHE_FILE}")
        with open(CACHE_FILE, "r") as f:
            raw_text = f.read()
        sha = hashlib.sha256(raw_text.encode("utf-8")).hexdigest()
        print(f"  Cache SHA256: {sha}")
        sha_file = CACHE_FILE + ".sha256"
        if os.path.exists(sha_file):
            with open(sha_file) as f:
                expected_sha = f.read().strip()
            if sha != expected_sha:
                print(f"  WARNING: SHA256 mismatch! Re-downloading...")
                os.remove(CACHE_FILE)
                return download_uniprot_data()
            print(f"  SHA256 verified.")
        return json.loads(raw_text)

    print(f"  Downloading from UniProt REST API (paginated, {PAGE_SIZE}/page)...")
    all_results = []
    page = 0
    next_url = f"{UNIPROT_BASE}?query={UNIPROT_QUERY}&format=json&size={PAGE_SIZE}"

    while next_url and len(all_results) < 500:
        page += 1
        print(f"  Page {page}: fetching up to {PAGE_SIZE} entries...")

        # Use urllib directly to access response headers for pagination
        import http.client
        for attempt in range(5):
            try:
                req = urllib.request.Request(next_url)
                req.add_header("User-Agent", "Python-urllib/3.8 (claw4s-research)")
                req.add_header("Accept", "application/json")
                resp = urllib.request.urlopen(req, timeout=180)

                # Read in chunks
                chunks = []
                while True:
                    try:
                        chunk = resp.read(65536)
                        if not chunk:
                            break
                        chunks.append(chunk)
                    except http.client.IncompleteRead as e:
                        chunks.append(e.partial)
                        break
                raw_bytes = b"".join(chunks)

                # Extract Link header for next page
                link_header = resp.headers.get("Link", "")
                resp.close()
                break
            except Exception as e:
                print(f"    Attempt {attempt+1}/5 failed: {e}")
                if attempt < 4:
                    time.sleep(3 * (attempt + 1))
                else:
                    raise

        raw_text = raw_bytes.decode("utf-8")
        page_data = json.loads(raw_text)

        results = page_data.get("results", [])
        if not results:
            break

        all_results.extend(results)
        print(f"  Page {page}: got {len(results)} entries (total: {len(all_results)})")

        # Parse Link header for next page URL
        # Format: <https://rest.uniprot.org/...>; rel="next"
        next_url = None
        if link_header:
            import re
            match = re.search(r'<([^>]+)>;\s*rel="next"', link_header)
            if match:
                next_url = match.group(1)

        if len(results) < PAGE_SIZE:
            break

        time.sleep(1)  # Rate limit courtesy

    data = {"results": all_results}

    # Cache locally
    raw_text = json.dumps(data)
    with open(CACHE_FILE, "w") as f:
        f.write(raw_text)

    sha = hashlib.sha256(raw_text.encode("utf-8")).hexdigest()
    with open(CACHE_FILE + ".sha256", "w") as f:
        f.write(sha)
    print(f"  Downloaded {len(all_results)} entries. SHA256: {sha}")

    return data


# === Parsing ===

def parse_proteins(data):
    """Parse UniProt JSON to extract protein info with disorder/structure annotations."""
    proteins = []
    results_list = data.get("results", [])

    for entry in results_list:
        accession = entry.get("primaryAccession", "unknown")

        # Get sequence
        seq_info = entry.get("sequence", {})
        sequence = seq_info.get("value", "")
        if not sequence:
            continue

        seq_len = len(sequence)
        features = entry.get("features", [])

        disordered_positions = set()
        structured_positions = set()

        for feat in features:
            feat_type = feat.get("type", "")
            description = (feat.get("description") or "").lower()

            loc = feat.get("location", {})
            start_info = loc.get("start", {})
            end_info = loc.get("end", {})
            start = start_info.get("value")
            end = end_info.get("value")

            if start is None or end is None:
                continue
            if not isinstance(start, int) or not isinstance(end, int):
                continue

            start_idx = start - 1  # Convert to 0-indexed
            end_idx = end          # UniProt end is inclusive, range needs exclusive

            if start_idx < 0 or end_idx > seq_len or start_idx >= end_idx:
                continue

            positions = set(range(start_idx, end_idx))

            # Disordered regions: Region features mentioning disorder
            if feat_type == "Region" and "disordered" in description:
                disordered_positions |= positions
            # Structured regions: helices and strands
            elif feat_type in ("Helix", "Beta strand"):
                structured_positions |= positions

        # Remove overlaps — prioritize disorder annotation
        structured_positions -= disordered_positions

        if not disordered_positions or not structured_positions:
            continue

        # Extract amino acids at annotated positions
        disordered_aas = [sequence[i] for i in sorted(disordered_positions)
                          if i < seq_len and sequence[i] in "ACDEFGHIKLMNPQRSTVWY"]
        structured_aas = [sequence[i] for i in sorted(structured_positions)
                          if i < seq_len and sequence[i] in "ACDEFGHIKLMNPQRSTVWY"]

        if disordered_aas and structured_aas:
            proteins.append({
                "accession": accession,
                "seq_len": seq_len,
                "n_disordered": len(disordered_aas),
                "n_structured": len(structured_aas),
                "disordered_aas": disordered_aas,
                "structured_aas": structured_aas,
            })

    return proteins


# === Statistical Functions ===

def log_comb(n, k):
    """Log of binomial coefficient C(n, k) using lgamma."""
    if k < 0 or k > n:
        return float("-inf")
    if k == 0 or k == n:
        return 0.0
    return math.lgamma(n + 1) - math.lgamma(k + 1) - math.lgamma(n - k + 1)


def hypergeom_sf(k, N, K, n):
    """P(X >= k) for hypergeometric(N, K, n). Upper tail for enrichment."""
    log_probs = []
    upper = min(K, n)
    for i in range(k, upper + 1):
        if 0 <= n - i <= N - K:
            lp = log_comb(K, i) + log_comb(N - K, n - i) - log_comb(N, n)
            log_probs.append(lp)
    if not log_probs:
        return 0.0
    max_lp = max(log_probs)
    if max_lp == float("-inf"):
        return 0.0
    return math.exp(max_lp + math.log(sum(math.exp(lp - max_lp) for lp in log_probs)))


def hypergeom_cdf(k, N, K, n):
    """P(X <= k) for hypergeometric(N, K, n). Lower tail for depletion."""
    log_probs = []
    lower = max(0, n - (N - K))
    for i in range(lower, k + 1):
        if 0 <= n - i <= N - K:
            lp = log_comb(K, i) + log_comb(N - K, n - i) - log_comb(N, n)
            log_probs.append(lp)
    if not log_probs:
        return 0.0
    max_lp = max(log_probs)
    if max_lp == float("-inf"):
        return 0.0
    return math.exp(max_lp + math.log(sum(math.exp(lp - max_lp) for lp in log_probs)))


def compute_enrichment(proteins):
    """Compute amino acid enrichment in disordered vs structured regions."""
    all_disordered = []
    all_structured = []
    for p in proteins:
        all_disordered.extend(p["disordered_aas"])
        all_structured.extend(p["structured_aas"])

    dis_counts = Counter(all_disordered)
    str_counts = Counter(all_structured)
    total_dis = sum(dis_counts.values())
    total_str = sum(str_counts.values())
    total_all = total_dis + total_str

    results = {}
    for aa in AMINO_ACIDS:
        k = dis_counts.get(aa, 0)
        K = k + str_counts.get(aa, 0)
        n = total_dis
        N = total_all

        expected = K * n / N if N > 0 else 0
        freq_dis = k / total_dis if total_dis > 0 else 0
        freq_str = str_counts.get(aa, 0) / total_str if total_str > 0 else 0
        ratio = freq_dis / freq_str if freq_str > 0 else float("inf")
        log2fc = math.log2(ratio) if 0 < ratio < float("inf") else 0.0

        p_enrichment = hypergeom_sf(k, N, K, n)
        p_depletion = hypergeom_cdf(k, N, K, n)

        results[aa] = {
            "observed_disordered": k,
            "expected_disordered": round(expected, 1),
            "freq_disordered": round(freq_dis, 6),
            "freq_structured": round(freq_str, 6),
            "enrichment_ratio": round(ratio, 4) if ratio != float("inf") else 99.0,
            "log2_fold_change": round(log2fc, 4),
            "p_enrichment": p_enrichment,
            "p_depletion": p_depletion,
            "direction": "enriched" if ratio > 1 else ("depleted" if ratio < 1 else "neutral"),
        }

    return results, all_disordered, all_structured


def bootstrap_enrichment_ratios(proteins, n_boot=N_BOOTSTRAP, seed=SEED):
    """Bootstrap CIs for enrichment ratios by resampling proteins."""
    rng = random.Random(seed)
    n = len(proteins)
    boot_ratios = {aa: [] for aa in AMINO_ACIDS}

    for _ in range(n_boot):
        sample = [proteins[rng.randint(0, n - 1)] for _ in range(n)]
        dis_counts = Counter()
        str_counts = Counter()
        for p in sample:
            for aa in p["disordered_aas"]:
                dis_counts[aa] += 1
            for aa in p["structured_aas"]:
                str_counts[aa] += 1

        total_dis = sum(dis_counts.values())
        total_str = sum(str_counts.values())

        for aa in AMINO_ACIDS:
            fd = dis_counts.get(aa, 0) / total_dis if total_dis > 0 else 0
            fs = str_counts.get(aa, 0) / total_str if total_str > 0 else 0
            if fs > 0:
                boot_ratios[aa].append(fd / fs)

    cis = {}
    for aa in AMINO_ACIDS:
        vals = sorted(boot_ratios[aa])
        if len(vals) >= 20:
            lo_idx = int(0.025 * len(vals))
            hi_idx = min(int(0.975 * len(vals)), len(vals) - 1)
            cis[aa] = {
                "ci_lower": round(vals[lo_idx], 4),
                "ci_upper": round(vals[hi_idx], 4),
                "boot_mean": round(sum(vals) / len(vals), 4),
                "n_valid": len(vals),
            }
        else:
            cis[aa] = {"ci_lower": None, "ci_upper": None, "boot_mean": None, "n_valid": len(vals)}

    return cis


def permutation_test(proteins, n_perm=N_PERMUTATIONS, seed=SEED):
    """Permutation test: shuffle disordered/structured labels within each protein."""
    rng = random.Random(seed)

    obs_results, _, _ = compute_enrichment(proteins)
    obs_ratios = {aa: obs_results[aa]["enrichment_ratio"] for aa in AMINO_ACIDS}

    perm_counts = {aa: 0 for aa in AMINO_ACIDS}

    for _ in range(n_perm):
        perm_proteins = []
        for prot in proteins:
            all_aas = prot["disordered_aas"] + prot["structured_aas"]
            shuffled = list(all_aas)
            rng.shuffle(shuffled)
            n_dis = len(prot["disordered_aas"])
            perm_proteins.append({
                "disordered_aas": shuffled[:n_dis],
                "structured_aas": shuffled[n_dis:],
            })

        perm_results, _, _ = compute_enrichment(perm_proteins)

        for aa in AMINO_ACIDS:
            pr = perm_results[aa]["enrichment_ratio"]
            if obs_ratios[aa] >= 1 and pr >= obs_ratios[aa]:
                perm_counts[aa] += 1
            elif obs_ratios[aa] < 1 and pr <= obs_ratios[aa]:
                perm_counts[aa] += 1

    return {aa: (perm_counts[aa] + 1) / (n_perm + 1) for aa in AMINO_ACIDS}


def naive_vs_corrected(proteins):
    """Compare naive enrichment (IDR vs all residues) to corrected (IDR vs structured only).

    The naive approach compares IDR composition to the pooled composition (IDR + structured),
    which contaminates the baseline with the signal. The corrected approach compares IDR
    composition to structured regions only, providing a clean external baseline.
    """
    all_disordered = []
    all_structured = []
    for p in proteins:
        all_disordered.extend(p["disordered_aas"])
        all_structured.extend(p["structured_aas"])

    dis_counts = Counter(all_disordered)
    str_counts = Counter(all_structured)
    all_counts = Counter(all_disordered + all_structured)

    total_dis = sum(dis_counts.values())
    total_str = sum(str_counts.values())
    total_all = total_dis + total_str

    comparison = {}
    for aa in AMINO_ACIDS:
        freq_dis = dis_counts.get(aa, 0) / total_dis if total_dis > 0 else 0
        freq_str = str_counts.get(aa, 0) / total_str if total_str > 0 else 0
        freq_all = all_counts.get(aa, 0) / total_all if total_all > 0 else 0

        naive_ratio = freq_dis / freq_all if freq_all > 0 else 0
        corrected_ratio = freq_dis / freq_str if freq_str > 0 else 0
        bias = naive_ratio - corrected_ratio

        comparison[aa] = {
            "freq_idr": round(freq_dis, 5),
            "freq_all": round(freq_all, 5),
            "freq_structured": round(freq_str, 5),
            "naive_ratio": round(naive_ratio, 4),
            "corrected_ratio": round(corrected_ratio, 4),
            "bias": round(bias, 4),
        }

    return comparison


def sensitivity_analysis(proteins, seed=SEED):
    """Vary subset size to test robustness of findings."""
    rng = random.Random(seed)
    subset_sizes = [50, 100, 200, len(proteins)]
    results = {}

    for size in subset_sizes:
        if size > len(proteins):
            continue
        if size == len(proteins):
            subset = proteins
        else:
            subset = rng.sample(proteins, size)
        enrichment, _, _ = compute_enrichment(subset)

        key_aas = {}
        for aa in ["P", "G", "S", "E", "Q", "K", "W", "C", "I", "L"]:
            key_aas[aa] = {
                "ratio": enrichment[aa]["enrichment_ratio"],
                "direction": enrichment[aa]["direction"],
            }
        results[f"n={size}"] = key_aas

    return results


def benjamini_hochberg(pvalues):
    """Benjamini-Hochberg FDR correction for multiple testing."""
    items = sorted(pvalues.items(), key=lambda x: x[1])
    n = len(items)
    corrected = {}
    for rank, (key, pval) in enumerate(items, 1):
        corrected[key] = min(pval * n / rank, 1.0)

    sorted_keys = [k for k, _ in items]
    for i in range(len(sorted_keys) - 2, -1, -1):
        corrected[sorted_keys[i]] = min(corrected[sorted_keys[i]],
                                         corrected[sorted_keys[i + 1]])
    return corrected


# === Output ===

def write_report(proteins, enrichment, cis, perm_pvalues, sensitivity, bh_corrected, naive_comp=None):
    """Write human-readable Markdown report."""
    total_dis = sum(len(p["disordered_aas"]) for p in proteins)
    total_str = sum(len(p["structured_aas"]) for p in proteins)

    lines = [
        "# Amino Acid Enrichment in Intrinsically Disordered Regions",
        "",
        "## Summary",
        f"- **Proteins analyzed**: {len(proteins)}",
        f"- **Disordered residues**: {total_dis:,}",
        f"- **Structured residues (helix + strand)**: {total_str:,}",
        f"- **Bootstrap resamples**: {N_BOOTSTRAP}",
        f"- **Permutations**: {N_PERMUTATIONS}",
        f"- **Random seed**: {SEED}",
        "",
        "## Enrichment Results",
        "",
        "| AA | Freq(IDR) | Freq(Str) | Ratio | log2FC | 95% CI | p(hyper) | p(perm) | BH-adj | Dir |",
        "|:--:|:---------:|:---------:|:-----:|:------:|:------:|:--------:|:-------:|:------:|:---:|",
    ]

    for aa in AMINO_ACIDS:
        e = enrichment[aa]
        ci = cis[aa]
        pp = perm_pvalues[aa]
        bh = bh_corrected.get(aa, 1.0)
        ci_str = (f"[{ci['ci_lower']:.3f}, {ci['ci_upper']:.3f}]"
                  if ci["ci_lower"] is not None else "N/A")
        p_val = e["p_enrichment"] if e["direction"] == "enriched" else e["p_depletion"]
        lines.append(
            f"| {aa} | {e['freq_disordered']:.4f} | {e['freq_structured']:.4f} | "
            f"{e['enrichment_ratio']:.3f} | {e['log2_fold_change']:+.3f} | {ci_str} | "
            f"{p_val:.2e} | {pp:.4f} | {bh:.4f} | {e['direction']} |"
        )

    lines.extend(["", "## Top Enriched in IDRs", ""])
    enriched = sorted(
        [(aa, enrichment[aa]) for aa in AMINO_ACIDS if enrichment[aa]["direction"] == "enriched"],
        key=lambda x: x[1]["enrichment_ratio"], reverse=True
    )
    for aa, e in enriched[:8]:
        ci = cis[aa]
        ci_str = (f"95% CI [{ci['ci_lower']:.3f}, {ci['ci_upper']:.3f}]"
                  if ci["ci_lower"] is not None else "")
        lines.append(f"- **{aa}**: {e['enrichment_ratio']:.3f}x enriched "
                      f"(log2FC={e['log2_fold_change']:+.3f}, {ci_str})")

    lines.extend(["", "## Top Depleted in IDRs", ""])
    depleted = sorted(
        [(aa, enrichment[aa]) for aa in AMINO_ACIDS if enrichment[aa]["direction"] == "depleted"],
        key=lambda x: x[1]["enrichment_ratio"]
    )
    for aa, e in depleted[:8]:
        ci = cis[aa]
        ci_str = (f"95% CI [{ci['ci_lower']:.3f}, {ci['ci_upper']:.3f}]"
                  if ci["ci_lower"] is not None else "")
        lines.append(f"- **{aa}**: {e['enrichment_ratio']:.3f}x "
                      f"(log2FC={e['log2_fold_change']:+.3f}, {ci_str})")

    lines.extend([
        "",
        "## Sensitivity Analysis",
        "",
    ])
    for subset_name, key_aas in sensitivity.items():
        items_str = ", ".join(f"{aa}={info['ratio']:.3f}" for aa, info in key_aas.items())
        lines.append(f"- **{subset_name}**: {items_str}")

    # Naive vs corrected comparison
    if naive_comp:
        lines.extend([
            "",
            "## Naive vs Composition-Corrected Enrichment",
            "",
            "The naive approach compares IDR composition to the pooled baseline (IDR + structured),",
            "contaminating the reference with the signal. The corrected approach uses structured regions only.",
            "",
            "| AA | Naive Ratio | Corrected Ratio | Bias (Naive - Corrected) |",
            "|:--:|:-----------:|:---------------:|:------------------------:|",
        ])
        for aa in AMINO_ACIDS:
            nc = naive_comp[aa]
            lines.append(f"| {aa} | {nc['naive_ratio']:.3f} | {nc['corrected_ratio']:.3f} | {nc['bias']:+.4f} |")
        lines.extend([
            "",
            "For enriched amino acids, the naive approach systematically underestimates enrichment",
            "(negative bias) because the pooled baseline is pulled toward the IDR composition.",
            "For depleted amino acids, the naive approach systematically overestimates the ratio",
            "(positive bias) for the same reason. The largest biases are for Pro (bias = "
            f"{naive_comp['P']['bias']:+.3f}) and Ser (bias = {naive_comp['S']['bias']:+.3f}).",
        ])

    # Flag Ala discordance
    lines.extend([
        "",
        "## Discordance Note: Alanine",
        "",
        "Alanine shows marginal depletion by the hypergeometric test (p=0.045) but the",
        "permutation test does not support this (p=0.996). The bootstrap CI spans 1.0",
        f"([{cis['A']['ci_lower']:.3f}, {cis['A']['ci_upper']:.3f}]). We conclude Ala enrichment/depletion",
        "is not reliably distinguishable from the null in this dataset.",
    ])

    lines.extend([
        "",
        "## Methodological Note",
        "",
        "Raw amino acid frequency counts in disordered regions can overstate enrichment",
        "because they ignore protein-level composition baselines. Our approach corrects this:",
        "",
        "1. **Hypergeometric test** conditions on the total count of each amino acid and the",
        "   total number of disordered residues, providing an exact null model.",
        "2. **Protein-level bootstrap** resamples whole proteins (not individual residues),",
        "   preserving within-protein correlations in amino acid usage.",
        "3. **Within-protein permutation** shuffles residue labels within each protein,",
        "   destroying the disorder-composition association while preserving per-protein",
        "   amino acid composition exactly.",
        "4. **Benjamini-Hochberg correction** controls the false discovery rate across 20",
        "   simultaneous amino acid tests.",
        "",
        "## Limitations",
        "",
        "1. Only UniProt-reviewed (Swiss-Prot) human proteins with annotated disordered regions are included",
        "2. Disorder annotations depend on experimental evidence and prediction tools available at curation time",
        "3. Binary classification (disordered vs structured) ignores the continuum of protein dynamics",
        "4. Does not account for phylogenetic non-independence among human paralogs",
        "5. 'Structured' defined as helix/strand only — excludes loops and coils that may be ordered but flexible",
        "6. Sample limited to 500 proteins; full proteome analysis may yield different effect sizes",
        "7. Cannot distinguish cause from consequence: are these amino acids enriched because they promote disorder, or are they tolerated in disordered contexts?",
    ])

    with open(REPORT_FILE, "w") as f:
        f.write("\n".join(lines) + "\n")


def write_results_json(proteins, enrichment, cis, perm_pvalues, sensitivity, bh_corrected, naive_comp=None):
    """Write structured JSON results."""
    total_dis = sum(len(p["disordered_aas"]) for p in proteins)
    total_str = sum(len(p["structured_aas"]) for p in proteins)

    aa_results = {}
    for aa in AMINO_ACIDS:
        e = enrichment[aa]
        ci = cis[aa]
        p_val = e["p_enrichment"] if e["direction"] == "enriched" else e["p_depletion"]
        aa_results[aa] = {
            "freq_disordered": e["freq_disordered"],
            "freq_structured": e["freq_structured"],
            "enrichment_ratio": e["enrichment_ratio"],
            "log2_fold_change": e["log2_fold_change"],
            "observed_disordered": e["observed_disordered"],
            "expected_disordered": e["expected_disordered"],
            "p_hypergeometric": p_val,
            "p_permutation": perm_pvalues[aa],
            "p_bh_adjusted": bh_corrected.get(aa, 1.0),
            "ci_lower": ci["ci_lower"],
            "ci_upper": ci["ci_upper"],
            "boot_mean": ci["boot_mean"],
            "direction": e["direction"],
        }

    output = {
        "metadata": {
            "n_proteins": len(proteins),
            "n_disordered_residues": total_dis,
            "n_structured_residues": total_str,
            "n_bootstrap": N_BOOTSTRAP,
            "n_permutations": N_PERMUTATIONS,
            "seed": SEED,
            "data_source": "UniProt REST API v2 (Swiss-Prot reviewed, Homo sapiens)",
            "query_url": f"{UNIPROT_BASE}?query={UNIPROT_QUERY}&format=json&size={PAGE_SIZE}",
            "cache_sha256_file": CACHE_FILE + ".sha256",
        },
        "amino_acid_enrichment": aa_results,
        "sensitivity_analysis": sensitivity,
        "naive_vs_corrected": naive_comp or {},
        "protein_accessions": [p["accession"] for p in proteins],
    }

    with open(RESULTS_FILE, "w") as f:
        json.dump(output, f, indent=2)
    return output


# === Verification ===

def verify():
    """Machine-checkable verification of results (8+ assertions)."""
    print("[VERIFY] Running verification checks...")
    errors = []

    # 1. results.json exists
    if not os.path.exists(RESULTS_FILE):
        errors.append("results.json not found")
        print(f"  [1/10] results.json exists: FAIL")
    else:
        print(f"  [1/10] results.json exists: PASS")

    # 2. report.md exists
    if not os.path.exists(REPORT_FILE):
        errors.append("report.md not found")
        print(f"  [2/10] report.md exists: FAIL")
    else:
        print(f"  [2/10] report.md exists: PASS")

    if errors:
        print(f"\n[VERIFY] FAILED: {'; '.join(errors)}")
        sys.exit(1)

    with open(RESULTS_FILE) as f:
        results = json.load(f)

    meta = results["metadata"]

    # 3. Sufficient proteins
    if meta["n_proteins"] < 50:
        errors.append(f"Too few proteins: {meta['n_proteins']}")
        print(f"  [3/10] Sufficient proteins ({meta['n_proteins']}): FAIL")
    else:
        print(f"  [3/10] Sufficient proteins ({meta['n_proteins']}): PASS")

    # 4. Sufficient residues
    if meta["n_disordered_residues"] < 500:
        errors.append(f"Too few disordered residues: {meta['n_disordered_residues']}")
        print(f"  [4/10] Sufficient disordered residues ({meta['n_disordered_residues']}): FAIL")
    else:
        print(f"  [4/10] Sufficient disordered residues ({meta['n_disordered_residues']}): PASS")

    # 5. All 20 amino acids present
    aa_results = results["amino_acid_enrichment"]
    if len(aa_results) != 20:
        errors.append(f"Expected 20 amino acids, got {len(aa_results)}")
        print(f"  [5/10] All 20 amino acids: FAIL")
    else:
        print(f"  [5/10] All 20 amino acids: PASS")

    # 6. All enrichment ratios valid with CIs
    all_valid = True
    for aa, info in aa_results.items():
        if info["enrichment_ratio"] <= 0 or info["ci_lower"] is None:
            all_valid = False
            errors.append(f"{aa} has invalid ratio or missing CI")
    print(f"  [6/10] All ratios valid with CIs: {'PASS' if all_valid else 'FAIL'}")

    # 7. Frequencies sum to ~1
    sum_dis = sum(info["freq_disordered"] for info in aa_results.values())
    sum_str = sum(info["freq_structured"] for info in aa_results.values())
    freq_ok = 0.95 < sum_dis < 1.05 and 0.95 < sum_str < 1.05
    if not freq_ok:
        errors.append(f"Frequencies don't sum to ~1: dis={sum_dis}, str={sum_str}")
    print(f"  [7/10] Frequencies sum to ~1 (dis={sum_dis:.4f}, str={sum_str:.4f}): "
          f"{'PASS' if freq_ok else 'FAIL'}")

    # 8. P-values in [0, 1]
    pvals_ok = True
    for aa, info in aa_results.items():
        if not (0 <= info["p_hypergeometric"] <= 1) or not (0 <= info["p_permutation"] <= 1):
            pvals_ok = False
            errors.append(f"{aa} has p-value out of range")
    print(f"  [8/10] All p-values in [0, 1]: {'PASS' if pvals_ok else 'FAIL'}")

    # 9. Sensitivity analysis present
    sa = results.get("sensitivity_analysis", {})
    sa_ok = len(sa) >= 2
    if not sa_ok:
        errors.append("Sensitivity analysis needs >= 2 conditions")
    print(f"  [9/10] Sensitivity analysis ({len(sa)} conditions): {'PASS' if sa_ok else 'FAIL'}")

    # 10. At least some amino acids are significantly enriched or depleted
    n_sig = sum(1 for info in aa_results.values() if info["p_bh_adjusted"] < 0.05)
    sig_ok = n_sig >= 1
    if not sig_ok:
        errors.append("No significant results after BH correction")
    print(f"  [10/10] Significant results ({n_sig}/20 at BH<0.05): {'PASS' if sig_ok else 'FAIL'}")

    if errors:
        print(f"\n[VERIFY] FAILED: {'; '.join(errors)}")
        sys.exit(1)
    else:
        print(f"\n[VERIFY] ALL CHECKS PASSED")


# === Main ===

def run_analysis():
    """Main analysis pipeline."""
    random.seed(SEED)

    print("=" * 70)
    print("AMINO ACID ENRICHMENT IN INTRINSICALLY DISORDERED REGIONS")
    print("OF HUMAN PROTEINS")
    print("=" * 70)

    # [1/7] Download data
    print(f"\n[1/8] Downloading UniProt data...")
    data = download_uniprot_data()
    n_entries = len(data.get("results", []))
    print(f"  Retrieved {n_entries} protein entries")

    # [2/8] Parse proteins
    print(f"\n[2/8] Parsing protein features...")
    proteins = parse_proteins(data)
    print(f"  {len(proteins)} proteins with BOTH disordered and structured annotations")
    total_dis = sum(len(p["disordered_aas"]) for p in proteins)
    total_str = sum(len(p["structured_aas"]) for p in proteins)
    print(f"  Disordered residues: {total_dis:,}")
    print(f"  Structured residues: {total_str:,}")

    if len(proteins) < 10:
        print("FATAL: Too few proteins with both disorder and structure annotations.")
        print("Check UniProt API response and feature parsing.")
        sys.exit(1)

    # [3/8] Hypergeometric enrichment
    print(f"\n[3/8] Computing amino acid enrichment (hypergeometric test)...")
    enrichment, all_dis, all_str = compute_enrichment(proteins)

    print(f"\n  {'AA':<4} {'Ratio':>8} {'log2FC':>8} {'p-value':>12} {'Direction':>10}")
    print(f"  {'-'*4} {'-'*8} {'-'*8} {'-'*12} {'-'*10}")
    for aa in AMINO_ACIDS:
        e = enrichment[aa]
        p_val = e["p_enrichment"] if e["direction"] == "enriched" else e["p_depletion"]
        print(f"  {aa:<4} {e['enrichment_ratio']:>8.3f} {e['log2_fold_change']:>+8.3f} "
              f"{p_val:>12.2e} {e['direction']:>10}")

    # [4/8] Bootstrap CIs
    print(f"\n[4/8] Bootstrap confidence intervals ({N_BOOTSTRAP} resamples)...")
    cis = bootstrap_enrichment_ratios(proteins)
    for aa in AMINO_ACIDS:
        ci = cis[aa]
        if ci["ci_lower"] is not None:
            print(f"  {aa}: [{ci['ci_lower']:.4f}, {ci['ci_upper']:.4f}] mean={ci['boot_mean']:.4f}")

    # [5/8] Permutation test
    print(f"\n[5/8] Permutation test ({N_PERMUTATIONS} within-protein permutations)...")
    perm_pvalues = permutation_test(proteins)
    for aa in AMINO_ACIDS:
        print(f"  {aa}: p_perm = {perm_pvalues[aa]:.4f}")

    # [6/8] Multiple testing correction
    print(f"\n[6/8] Benjamini-Hochberg FDR correction...")
    raw_pvalues = {}
    for aa in AMINO_ACIDS:
        e = enrichment[aa]
        raw_pvalues[aa] = e["p_enrichment"] if e["direction"] == "enriched" else e["p_depletion"]
    bh_corrected = benjamini_hochberg(raw_pvalues)

    n_sig = sum(1 for v in bh_corrected.values() if v < 0.05)
    print(f"  {n_sig}/20 amino acids significant at BH-adjusted p < 0.05")
    for aa in AMINO_ACIDS:
        if bh_corrected[aa] < 0.05:
            print(f"  * {aa}: BH-adjusted p = {bh_corrected[aa]:.4e} ({enrichment[aa]['direction']})")

    # [7/8] Naive vs corrected comparison
    print(f"\n[7/8] Naive vs corrected enrichment comparison...")
    naive_comp = naive_vs_corrected(proteins)
    print(f"  {'AA':<4} {'Naive':>8} {'Corrected':>10} {'Bias':>8}")
    print(f"  {'-'*4} {'-'*8} {'-'*10} {'-'*8}")
    for aa in AMINO_ACIDS:
        nc = naive_comp[aa]
        print(f"  {aa:<4} {nc['naive_ratio']:>8.3f} {nc['corrected_ratio']:>10.3f} {nc['bias']:>+8.4f}")

    # [8/8] Sensitivity analysis
    print(f"\n[8/8] Sensitivity analysis (varying subset sizes)...")
    sensitivity = sensitivity_analysis(proteins)
    for subset_name, key_aas in sensitivity.items():
        items_str = ", ".join(f"{aa}={info['ratio']:.3f}" for aa, info in key_aas.items())
        print(f"  {subset_name}: {items_str}")

    # Write outputs
    print(f"\nWriting outputs...")
    output = write_results_json(proteins, enrichment, cis, perm_pvalues, sensitivity, bh_corrected, naive_comp)
    write_report(proteins, enrichment, cis, perm_pvalues, sensitivity, bh_corrected, naive_comp)

    print(f"  results.json: {RESULTS_FILE}")
    print(f"  report.md:    {REPORT_FILE}")
    print(f"\nANALYSIS COMPLETE")


if __name__ == "__main__":
    if "--verify" in sys.argv:
        verify()
    else:
        run_analysis()
SCRIPT_EOF
```

**Expected output:** No stdout. File `analysis.py` created in workspace.

## Step 3: Run analysis

```bash
cd /tmp/claw4s_auto_uniprot-amino-acid-disorder-enrichment && python3 analysis.py
```

**Expected output:**
- Sectioned progress from `[1/8]` through `[8/8]`
- Per-amino-acid enrichment ratios, bootstrap CIs, permutation p-values
- Naive vs corrected enrichment comparison showing composition baseline bias
- Ends with `ANALYSIS COMPLETE`
- Creates `results.json` and `report.md` in workspace

**Estimated runtime:** 5–15 minutes (dominated by 1000 permutations and 2000 bootstrap resamples).

## Step 4: Verify results

```bash
cd /tmp/claw4s_auto_uniprot-amino-acid-disorder-enrichment && python3 analysis.py --verify
```

**Expected output:**
- 10 numbered checks `[1/10]` through `[10/10]`, each `PASS`
- Ends with `[VERIFY] ALL CHECKS PASSED`

## Success Criteria

1. `results.json` exists with enrichment ratios, CIs, and p-values for all 20 amino acids
2. `report.md` exists with formatted tables and findings
3. All 10 verification checks pass
4. At least one amino acid shows significant enrichment/depletion after BH correction
5. Bootstrap CIs are non-null for all amino acids
6. Sensitivity analysis shows consistent direction across subset sizes
7. Naive vs corrected comparison quantifies composition baseline bias

## Failure Conditions

- UniProt API unreachable → retry with exponential backoff; fail after 5 attempts
- Fewer than 10 proteins with both disorder and structure annotations → FATAL, check API query
- SHA256 mismatch on cached data → re-download automatically
- Any verification check fails → review analysis output for errors