---
name: "PDB Resolution Bias by Source Organism"
description: "Tests whether X-ray crystal structures from different source organisms have systematically different resolution distributions using Kolmogorov-Smirnov tests with permutation-based p-values"
version: "1.0.0"
author: "Claw \U0001F99E, David Austin"
tags:
  - claw4s-2026
  - structural-biology
  - protein-data-bank
  - resolution-bias
  - permutation-test
  - kolmogorov-smirnov
python_version: ">=3.8"
dependencies: []
---

# PDB Resolution Bias by Source Organism

## Motivation

The Protein Data Bank (PDB) contains structures solved at varying resolutions. If resolution
distributions differ systematically by source organism, this could confound cross-species
structural comparisons. We test whether human-derived X-ray structures have significantly
different resolution distributions than E. coli, yeast, and mouse structures using
Kolmogorov-Smirnov tests with permutation-based p-values (5,000 shuffles), complemented by Welch's t-test.

**Methodological hook:** Most comparative structural studies assume resolution is organism-
independent. We provide a rigorous permutation-based test of this assumption across four
major model organisms, with bootstrap confidence intervals for effect sizes and multi-
parameter sensitivity analysis.

## Step 1: Create workspace

```bash
mkdir -p /tmp/claw4s_auto_pdb-resolution-bias-by-organism
```

**Expected output:** Directory created (exit code 0).

## Step 2: Write analysis script

```bash
cat << 'SCRIPT_EOF' > /tmp/claw4s_auto_pdb-resolution-bias-by-organism/pdb_resolution_bias.py
#!/usr/bin/env python3
"""
PDB Resolution Bias by Source Organism
======================================
Tests whether X-ray crystal structures from different source organisms
have systematically different resolution distributions.

Uses Kolmogorov-Smirnov tests with permutation-based p-values,
bootstrap confidence intervals for effect sizes, and sensitivity analysis.

Data: RCSB Protein Data Bank (10,000 recent X-ray entries)
Null model: Permutation test (2,000 label shuffles)
Python 3.8+ standard library only.
"""

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

# ─── Configuration ───────────────────────────────────────────────────
SEED = 42
N_PERMUTATIONS = 5000
N_BOOTSTRAP = 2000
SEARCH_API = "https://search.rcsb.org/rcsbsearch/v2/query"
GRAPHQL_API = "https://data.rcsb.org/graphql"
ENTRY_API_BASE = "https://data.rcsb.org/rest/v1/core/entry"
TARGET_N = 10000
GQL_BATCH_SIZE = 100
CACHE_FILE = "pdb_metadata_cache.json"
RESULTS_FILE = "results.json"
REPORT_FILE = "report.md"

TARGET_ORGANISMS = {
    "Homo sapiens": "Human",
    "Escherichia coli": "E. coli",
    "Saccharomyces cerevisiae": "Yeast",
    "Mus musculus": "Mouse",
}
REFERENCE = "Homo sapiens"
COMPARISONS = ["Escherichia coli", "Saccharomyces cerevisiae", "Mus musculus"]


# ─── Basic Statistics (stdlib only) ──────────────────────────────────
def mean(xs):
    return sum(xs) / len(xs) if xs else 0.0


def median(xs):
    s = sorted(xs)
    n = len(s)
    if n == 0:
        return 0.0
    if n % 2 == 1:
        return s[n // 2]
    return (s[n // 2 - 1] + s[n // 2]) / 2.0


def stdev(xs):
    if len(xs) < 2:
        return 0.0
    m = mean(xs)
    return math.sqrt(sum((x - m) ** 2 for x in xs) / (len(xs) - 1))


def percentile(xs, p):
    """Percentile using linear interpolation."""
    s = sorted(xs)
    n = len(s)
    if n == 0:
        return 0.0
    k = (n - 1) * p / 100.0
    f = math.floor(k)
    c = math.ceil(k)
    if f == c:
        return s[int(k)]
    return s[int(f)] * (c - k) + s[int(c)] * (k - f)


def iqr(xs):
    return percentile(xs, 75) - percentile(xs, 25)


# ─── Statistical Tests ──────────────────────────────────────────────
def ks_statistic(s1, s2):
    """Two-sample Kolmogorov-Smirnov statistic via sorted merge.
    O(n log n) where n = len(s1) + len(s2)."""
    n1, n2 = len(s1), len(s2)
    if n1 == 0 or n2 == 0:
        return 0.0
    tagged = sorted([(v, 0) for v in s1] + [(v, 1) for v in s2])
    e1 = e2 = max_d = 0.0
    for val, grp in tagged:
        if grp == 0:
            e1 += 1.0 / n1
        else:
            e2 += 1.0 / n2
        d = abs(e1 - e2)
        if d > max_d:
            max_d = d
    return max_d


def permutation_ks_test(s1, s2, n_perm=N_PERMUTATIONS, seed=SEED):
    """Permutation-based two-sample KS test.
    Returns (observed_ks, p_value, n_exceeding).
    P-value uses conservative estimator: (exceed + 1) / (n_perm + 1)."""
    rng = random.Random(seed)
    observed = ks_statistic(s1, s2)
    combined = list(s1) + list(s2)
    n1 = len(s1)
    exceed = 0
    for i in range(n_perm):
        rng.shuffle(combined)
        if ks_statistic(combined[:n1], combined[n1:]) >= observed:
            exceed += 1
        if (i + 1) % 500 == 0:
            sys.stdout.write(f"    {i + 1}/{n_perm} permutations done\n")
            sys.stdout.flush()
    p_value = (exceed + 1) / (n_perm + 1)
    return observed, p_value, exceed


def cohens_d(s1, s2):
    """Cohen's d effect size with pooled standard deviation."""
    m1, m2 = mean(s1), mean(s2)
    n1, n2 = len(s1), len(s2)
    if n1 < 2 or n2 < 2:
        return 0.0
    v1 = sum((x - m1) ** 2 for x in s1) / (n1 - 1)
    v2 = sum((x - m2) ** 2 for x in s2) / (n2 - 1)
    pooled = math.sqrt(((n1 - 1) * v1 + (n2 - 1) * v2) / (n1 + n2 - 2))
    return (m1 - m2) / pooled if pooled > 0 else 0.0


def welch_t_test(s1, s2):
    """Welch's t-test (unequal variances). Returns (t_statistic, dof, p_value).
    P-value approximated via the regularized incomplete beta function."""
    m1, m2 = mean(s1), mean(s2)
    n1, n2 = len(s1), len(s2)
    if n1 < 2 or n2 < 2:
        return 0.0, 0, 1.0
    v1 = sum((x - m1) ** 2 for x in s1) / (n1 - 1)
    v2 = sum((x - m2) ** 2 for x in s2) / (n2 - 1)
    se = math.sqrt(v1 / n1 + v2 / n2)
    if se == 0:
        return 0.0, 0, 1.0
    t_stat = (m1 - m2) / se
    # Welch-Satterthwaite degrees of freedom
    num = (v1 / n1 + v2 / n2) ** 2
    den = (v1 / n1) ** 2 / (n1 - 1) + (v2 / n2) ** 2 / (n2 - 1)
    dof = num / den if den > 0 else 1.0
    # P-value via beta function approximation
    p_value = _t_pvalue(abs(t_stat), dof)
    return t_stat, dof, p_value


def _t_pvalue(t, dof):
    """Two-tailed p-value for t-distribution using incomplete beta function."""
    x = dof / (dof + t * t)
    # Regularized incomplete beta function via continued fraction
    p = _betainc(dof / 2.0, 0.5, x)
    return min(p, 1.0)


def _betainc(a, b, x):
    """Regularized incomplete beta function I_x(a, b) via continued fraction."""
    if x < 0 or x > 1:
        return 0.0
    if x == 0 or x == 1:
        return x
    # Use continued fraction (Lentz's method)
    lbeta = math.lgamma(a) + math.lgamma(b) - math.lgamma(a + b)
    front = math.exp(math.log(x) * a + math.log(1 - x) * b - lbeta) / a
    # Modified Lentz's continued fraction
    f = 1.0
    c = 1.0
    d = 1.0 - (a + b) * x / (a + 1)
    if abs(d) < 1e-30:
        d = 1e-30
    d = 1.0 / d
    f = d
    for i in range(1, 201):
        m = i
        # Even step
        numerator = m * (b - m) * x / ((a + 2 * m - 1) * (a + 2 * m))
        d = 1.0 + numerator * d
        if abs(d) < 1e-30:
            d = 1e-30
        c = 1.0 + numerator / c
        if abs(c) < 1e-30:
            c = 1e-30
        d = 1.0 / d
        f *= d * c
        # Odd step
        numerator = -(a + m) * (a + b + m) * x / ((a + 2 * m) * (a + 2 * m + 1))
        d = 1.0 + numerator * d
        if abs(d) < 1e-30:
            d = 1e-30
        c = 1.0 + numerator / c
        if abs(c) < 1e-30:
            c = 1e-30
        d = 1.0 / d
        delta = d * c
        f *= delta
        if abs(delta - 1.0) < 1e-10:
            break
    return front * f


def median_diff(s1, s2):
    """Difference of medians."""
    return median(s1) - median(s2)


def bootstrap_ci(s1, s2, stat_fn, n_boot=N_BOOTSTRAP, seed=SEED, alpha=0.05):
    """Bootstrap confidence interval for a two-sample statistic.
    Returns (observed, ci_lower, ci_upper)."""
    rng = random.Random(seed)
    observed = stat_fn(s1, s2)
    boot_vals = []
    for _ in range(n_boot):
        b1 = [s1[rng.randint(0, len(s1) - 1)] for _ in range(len(s1))]
        b2 = [s2[rng.randint(0, len(s2) - 1)] for _ in range(len(s2))]
        boot_vals.append(stat_fn(b1, b2))
    boot_vals.sort()
    lo = boot_vals[max(0, int(n_boot * alpha / 2))]
    hi = boot_vals[min(n_boot - 1, int(n_boot * (1 - alpha / 2)))]
    return observed, lo, hi


# ─── Data Fetching ──────────────────────────────────────────────────
def url_fetch(url, data=None, headers=None, timeout=60, retries=3):
    """Fetch URL with retry logic and exponential backoff."""
    if headers is None:
        headers = {}
    ctx = ssl.create_default_context()
    for attempt in range(retries):
        try:
            req = urllib.request.Request(url, data=data, headers=headers)
            with urllib.request.urlopen(req, timeout=timeout, context=ctx) as resp:
                return resp.read().decode()
        except Exception as e:
            if attempt < retries - 1:
                wait = 2 ** (attempt + 1)
                print(f"  Retry {attempt + 1}/{retries} after {wait}s: {e}")
                time.sleep(wait)
            else:
                raise


def fetch_entry_ids():
    """Fetch 10,000 recent X-ray diffraction entry IDs from RCSB Search API."""
    query = {
        "query": {
            "type": "terminal",
            "service": "text",
            "parameters": {
                "attribute": "exptl.method",
                "operator": "exact_match",
                "value": "X-RAY DIFFRACTION",
            },
        },
        "return_type": "entry",
        "request_options": {
            "sort": [
                {
                    "sort_by": "rcsb_accession_info.initial_release_date",
                    "direction": "desc",
                }
            ],
            "paginate": {"start": 0, "rows": TARGET_N},
        },
    }
    resp = url_fetch(
        SEARCH_API,
        data=json.dumps(query).encode(),
        headers={"Content-Type": "application/json"},
        timeout=120,
    )
    result = json.loads(resp)
    ids = [r["identifier"] for r in result["result_set"]]
    print(f"  Retrieved {len(ids)} entry IDs from RCSB Search API")
    return ids


def fetch_metadata(entry_ids):
    """Fetch resolution and organism data via RCSB GraphQL API in batches.
    Caches results locally with SHA256 integrity verification."""
    if os.path.exists(CACHE_FILE):
        print(f"  Found cache: {CACHE_FILE}")
        with open(CACHE_FILE) as f:
            cached = json.load(f)
        verify_bytes = json.dumps(cached["entries"], sort_keys=True).encode()
        sha = hashlib.sha256(verify_bytes).hexdigest()
        if sha == cached.get("sha256"):
            print(f"  Cache integrity verified (SHA256: {sha[:16]}...)")
            print(f"  Loaded {len(cached['entries'])} cached entries")
            return cached["entries"]
        print("  Cache integrity check FAILED, re-fetching")

    all_entries = []
    n_batches = (len(entry_ids) + GQL_BATCH_SIZE - 1) // GQL_BATCH_SIZE
    failed_batches = 0

    for bi in range(n_batches):
        start = bi * GQL_BATCH_SIZE
        batch = entry_ids[start : start + GQL_BATCH_SIZE]

        if (bi + 1) % 20 == 0 or bi == 0 or bi == n_batches - 1:
            pct = 100.0 * (bi + 1) / n_batches
            print(
                f"  GraphQL batch {bi + 1}/{n_batches} ({pct:.0f}%) "
                f"- {len(all_entries)} entries fetched"
            )

        # Build GraphQL query with aliases
        fragments = []
        for i, eid in enumerate(batch):
            fragments.append(
                f'e{i}: entry(entry_id: "{eid}") {{ '
                f"rcsb_id "
                f"rcsb_entry_info {{ resolution_combined }} "
                f"polymer_entities {{ "
                f"rcsb_entity_source_organism {{ ncbi_scientific_name }} "
                f"}} }}"
            )
        gql = "{ " + " ".join(fragments) + " }"

        try:
            resp = url_fetch(
                GRAPHQL_API,
                data=json.dumps({"query": gql}).encode(),
                headers={"Content-Type": "application/json"},
                timeout=120,
            )
            result = json.loads(resp)
            data_block = result.get("data", {})

            if "errors" in result and not data_block:
                msg = result["errors"][0].get("message", "unknown")
                print(f"  GraphQL error batch {bi + 1}: {msg}")
                failed_batches += 1
                if failed_batches > 20:
                    print("  Too many GraphQL failures, stopping")
                    break
                time.sleep(2)
                continue

            for i, eid in enumerate(batch):
                entry_data = data_block.get(f"e{i}")
                if entry_data is None:
                    continue

                resolution = None
                info = entry_data.get("rcsb_entry_info")
                if info and info.get("resolution_combined"):
                    resolution = info["resolution_combined"][0]

                organisms = []
                for pe in entry_data.get("polymer_entities") or []:
                    if pe is None:
                        continue
                    for org in pe.get("rcsb_entity_source_organism") or []:
                        if org and org.get("ncbi_scientific_name"):
                            organisms.append(org["ncbi_scientific_name"])

                all_entries.append(
                    {
                        "entry_id": eid,
                        "resolution": resolution,
                        "organisms": list(set(organisms)),
                    }
                )
        except Exception as e:
            print(f"  ERROR in batch {bi + 1}: {e}")
            failed_batches += 1
            if failed_batches > 20:
                print("  Too many failures, stopping")
                break
            time.sleep(3)
            continue

        # Rate limiting
        if bi < n_batches - 1:
            time.sleep(0.25)

    # Save cache with SHA256
    cache_bytes = json.dumps(all_entries, sort_keys=True).encode()
    sha = hashlib.sha256(cache_bytes).hexdigest()
    cache_obj = {
        "entries": all_entries,
        "sha256": sha,
        "n_entries": len(all_entries),
        "fetched_at": time.strftime("%Y-%m-%dT%H:%M:%SZ"),
        "source": "RCSB PDB (search: https://search.rcsb.org, data: https://data.rcsb.org/graphql)",
    }
    with open(CACHE_FILE, "w") as f:
        json.dump(cache_obj, f, indent=2)
    print(f"  Cached {len(all_entries)} entries (SHA256: {sha[:16]}...)")
    return all_entries


# ─── Analysis ───────────────────────────────────────────────────────
def organize_data(entries):
    """Group resolution values by target organism."""
    org_res = defaultdict(list)
    no_resolution = 0
    no_target_org = 0

    for e in entries:
        if e["resolution"] is None:
            no_resolution += 1
            continue
        matched = False
        for org_name in e["organisms"]:
            if org_name in TARGET_ORGANISMS:
                org_res[org_name].append(e["resolution"])
                matched = True
        if not matched:
            no_target_org += 1

    print(f"  Entries without resolution: {no_resolution}")
    print(f"  Entries not matching target organisms: {no_target_org}")
    for org in [REFERENCE] + COMPARISONS:
        vals = org_res.get(org, [])
        if vals:
            print(
                f"  {TARGET_ORGANISMS[org]:10s}: n={len(vals):5d}, "
                f"median={median(vals):.2f} A, mean={mean(vals):.2f} A, "
                f"SD={stdev(vals):.2f} A"
            )
    return dict(org_res)


def run_comparisons(org_data):
    """KS tests with permutation p-values and bootstrap CIs."""
    ref = org_data.get(REFERENCE, [])
    if not ref:
        print("  ERROR: No data for reference organism (Homo sapiens)")
        return []

    results = []
    for comp_org in COMPARISONS:
        comp = org_data.get(comp_org, [])
        if len(comp) < 30:
            print(
                f"  SKIP {TARGET_ORGANISMS[comp_org]}: only {len(comp)} entries (need >= 30)"
            )
            continue

        label = f"{TARGET_ORGANISMS[REFERENCE]} vs {TARGET_ORGANISMS[comp_org]}"
        print(f"\n  --- {label} ---")
        print(f"  Sample sizes: n1={len(ref)}, n2={len(comp)}")

        # KS test with permutation p-value
        print(f"  Computing KS statistic with {N_PERMUTATIONS} permutations...")
        ks_stat, ks_pval, ks_exceed = permutation_ks_test(ref, comp)
        print(f"  KS = {ks_stat:.4f}, p = {ks_pval:.6f} ({ks_exceed}/{N_PERMUTATIONS} exceeded)")

        # Bootstrap CI for median difference
        print(f"  Bootstrap CI for median difference ({N_BOOTSTRAP} resamples)...")
        md_obs, md_lo, md_hi = bootstrap_ci(ref, comp, median_diff, N_BOOTSTRAP, SEED)
        print(f"  Median diff = {md_obs:.4f} A [95% CI: {md_lo:.4f}, {md_hi:.4f}]")

        # Cohen's d with bootstrap CI
        print(f"  Bootstrap CI for Cohen's d ({N_BOOTSTRAP} resamples)...")
        d_obs, d_lo, d_hi = bootstrap_ci(ref, comp, cohens_d, N_BOOTSTRAP, SEED + 100)
        print(f"  Cohen's d = {d_obs:.4f} [95% CI: {d_lo:.4f}, {d_hi:.4f}]")

        # Welch's t-test (parametric complement)
        t_stat, t_dof, t_pval = welch_t_test(ref, comp)
        print(f"  Welch's t = {t_stat:.4f}, dof = {t_dof:.1f}, p = {t_pval:.2e}")

        # Bonferroni correction
        bonf_pval = min(ks_pval * len(COMPARISONS), 1.0)
        bonf_t_pval = min(t_pval * len(COMPARISONS), 1.0)
        sig = bonf_pval < 0.05

        # Note: p-value floor = 1/(N_PERMUTATIONS+1)
        pval_floor = 1.0 / (N_PERMUTATIONS + 1)

        results.append(
            {
                "comparison": label,
                "reference": REFERENCE,
                "comparison_organism": comp_org,
                "n_reference": len(ref),
                "n_comparison": len(comp),
                "ref_median": round(median(ref), 4),
                "ref_mean": round(mean(ref), 4),
                "ref_sd": round(stdev(ref), 4),
                "comp_median": round(median(comp), 4),
                "comp_mean": round(mean(comp), 4),
                "comp_sd": round(stdev(comp), 4),
                "ks_statistic": round(ks_stat, 6),
                "permutation_p_value": round(ks_pval, 6),
                "permutations_exceeding": ks_exceed,
                "n_permutations": N_PERMUTATIONS,
                "permutation_p_value_floor": round(pval_floor, 6),
                "bonferroni_p_value": round(bonf_pval, 6),
                "welch_t_statistic": round(t_stat, 4),
                "welch_dof": round(t_dof, 2),
                "welch_p_value": round(t_pval, 10),
                "welch_bonferroni_p_value": round(bonf_t_pval, 10),
                "median_diff_angstrom": round(md_obs, 4),
                "median_diff_ci_lo": round(md_lo, 4),
                "median_diff_ci_hi": round(md_hi, 4),
                "cohens_d": round(d_obs, 4),
                "cohens_d_ci_lo": round(d_lo, 4),
                "cohens_d_ci_hi": round(d_hi, 4),
                "significant_bonferroni_005": sig,
            }
        )

    return results


def sensitivity_analysis(org_data):
    """Test robustness under different parameter choices."""
    ref = org_data.get(REFERENCE, [])
    results = {}

    # 1. Vary permutation count
    print("\n  Sensitivity: varying permutation count...")
    perm_sensitivity = {}
    for comp_org in COMPARISONS:
        comp = org_data.get(comp_org, [])
        if len(comp) < 30:
            continue
        name = TARGET_ORGANISMS[comp_org]
        perm_sensitivity[name] = {}
        for np_count in [1000, 2000, 5000]:
            _, pv, _ = permutation_ks_test(ref, comp, n_perm=np_count, seed=SEED)
            perm_sensitivity[name][str(np_count)] = round(pv, 6)
            print(f"    {name}, n_perm={np_count}: p={pv:.6f}")
    results["permutation_count_sensitivity"] = perm_sensitivity

    # 2. Vary sample fraction (subsampling stability)
    print("\n  Sensitivity: varying sample fraction...")
    rng = random.Random(SEED + 200)
    sample_sensitivity = {}
    for comp_org in COMPARISONS:
        comp = org_data.get(comp_org, [])
        if len(comp) < 30:
            continue
        name = TARGET_ORGANISMS[comp_org]
        sample_sensitivity[name] = {}
        for frac in [0.5, 0.75, 1.0]:
            n1 = max(10, int(len(ref) * frac))
            n2 = max(10, int(len(comp) * frac))
            sub_ref = rng.sample(ref, min(n1, len(ref)))
            sub_comp = rng.sample(comp, min(n2, len(comp)))
            ks = ks_statistic(sub_ref, sub_comp)
            sample_sensitivity[name][str(frac)] = round(ks, 6)
            print(f"    {name}, frac={frac}: KS={ks:.4f}")
    results["sample_fraction_sensitivity"] = sample_sensitivity

    # 3. Resolution threshold sensitivity
    print("\n  Sensitivity: varying resolution threshold...")
    threshold_sensitivity = {}
    for comp_org in COMPARISONS:
        comp_all = org_data.get(comp_org, [])
        if len(comp_all) < 30:
            continue
        name = TARGET_ORGANISMS[comp_org]
        threshold_sensitivity[name] = {}
        for thresh in [3.0, 4.0, 5.0, 10.0]:
            ref_f = [x for x in ref if x <= thresh]
            comp_f = [x for x in comp_all if x <= thresh]
            if len(ref_f) < 10 or len(comp_f) < 10:
                continue
            ks = ks_statistic(ref_f, comp_f)
            threshold_sensitivity[name][f"<={thresh}"] = {
                "ks_statistic": round(ks, 6),
                "n_ref": len(ref_f),
                "n_comp": len(comp_f),
            }
            print(
                f"    {name}, res<={thresh}A: KS={ks:.4f} "
                f"(n={len(ref_f)}+{len(comp_f)})"
            )
    results["resolution_threshold_sensitivity"] = threshold_sensitivity

    return results


def write_results(org_data, comparisons, sensitivity):
    """Write results.json and report.md."""
    # Descriptive statistics
    descriptive = {}
    for org, vals in org_data.items():
        if org not in TARGET_ORGANISMS:
            continue
        descriptive[TARGET_ORGANISMS[org]] = {
            "scientific_name": org,
            "n": len(vals),
            "mean": round(mean(vals), 4),
            "median": round(median(vals), 4),
            "sd": round(stdev(vals), 4),
            "iqr": round(iqr(vals), 4),
            "min": round(min(vals), 4),
            "max": round(max(vals), 4),
            "pct_25": round(percentile(vals, 25), 4),
            "pct_75": round(percentile(vals, 75), 4),
        }

    results_obj = {
        "analysis": "PDB Resolution Bias by Source Organism",
        "data_source": "RCSB Protein Data Bank (https://data.rcsb.org)",
        "search_api": SEARCH_API,
        "graphql_api": GRAPHQL_API,
        "entry_api_base": ENTRY_API_BASE,
        "method": "Kolmogorov-Smirnov test with permutation-based p-values",
        "n_entries_queried": TARGET_N,
        "n_total_entries_analyzed": sum(len(v) for v in org_data.values()),
        "seed": SEED,
        "n_permutations": N_PERMUTATIONS,
        "n_bootstrap": N_BOOTSTRAP,
        "reference_organism": REFERENCE,
        "descriptive_statistics": descriptive,
        "comparisons": comparisons,
        "sensitivity_analysis": sensitivity,
        "bonferroni_correction": True,
        "n_comparisons": len(COMPARISONS),
    }

    with open(RESULTS_FILE, "w") as f:
        json.dump(results_obj, f, indent=2)
    print(f"  Wrote {RESULTS_FILE}")

    # Write report.md
    lines = [
        "# PDB Resolution Bias by Source Organism",
        "",
        "## Summary",
        "",
        "Analyzed resolution distributions of X-ray crystal structures from "
        "four model organisms in the RCSB Protein Data Bank.",
        f"Reference organism: {TARGET_ORGANISMS[REFERENCE]} ({REFERENCE}).",
        "",
        "## Descriptive Statistics",
        "",
        "| Organism | N | Median (A) | Mean (A) | SD (A) | IQR (A) |",
        "|----------|---|-----------|----------|--------|---------|",
    ]
    for org in [REFERENCE] + COMPARISONS:
        if org not in org_data:
            continue
        d = descriptive[TARGET_ORGANISMS[org]]
        lines.append(
            f"| {TARGET_ORGANISMS[org]} | {d['n']} | {d['median']:.2f} | "
            f"{d['mean']:.2f} | {d['sd']:.2f} | {d['iqr']:.2f} |"
        )

    lines += ["", "## Statistical Comparisons", ""]

    for c in comparisons:
        sig_str = "**significant**" if c["significant_bonferroni_005"] else "not significant"
        floor_note = ""
        if c["permutations_exceeding"] == 0:
            floor_note = f" (at floor; true p < {c['permutation_p_value_floor']:.4f})"
        lines += [
            f"### {c['comparison']}",
            "",
            f"- **KS statistic:** {c['ks_statistic']:.4f}",
            f"- **Permutation p-value:** {c['permutation_p_value']:.6f} "
            f"({c['permutations_exceeding']}/{c['n_permutations']} exceeded){floor_note}",
            f"- **Welch's t-test:** t = {c['welch_t_statistic']:.3f}, "
            f"dof = {c['welch_dof']:.1f}, p = {c['welch_p_value']:.2e}",
            f"- **Bonferroni-corrected p (permutation):** {c['bonferroni_p_value']:.6f} ({sig_str})",
            f"- **Median difference:** {c['median_diff_angstrom']:.3f} A "
            f"[95% CI: {c['median_diff_ci_lo']:.3f}, {c['median_diff_ci_hi']:.3f}]",
            f"- **Cohen's d:** {c['cohens_d']:.3f} "
            f"[95% CI: {c['cohens_d_ci_lo']:.3f}, {c['cohens_d_ci_hi']:.3f}]",
            "",
        ]

    lines += [
        "## Sensitivity Analysis",
        "",
        "Results tested across:",
        "- Permutation counts: 1000, 2000, 5000",
        "- Sample fractions: 50%, 75%, 100%",
        "- Resolution thresholds: <=3.0, <=4.0, <=5.0, <=10.0 A",
        "",
        "See results.json for full sensitivity analysis details.",
        "",
        "## Limitations",
        "",
        "1. Resolution is self-reported by depositors; accuracy varies across labs",
        "2. Multi-organism entries (complexes) are counted in all matching organism groups",
        "3. Only four model organisms are compared; other species may show different patterns",
        "4. Results reflect a temporal snapshot of the PDB and may shift with new depositions",
        "5. Confounders (protein size, crystal packing, beamline technology) are not controlled",
        "6. Analysis restricted to X-ray diffraction; cryo-EM and NMR are excluded",
        "7. The search API returns entries sorted by release date, introducing temporal correlation",
    ]

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


# ─── Verification ───────────────────────────────────────────────────
def verify_results():
    """Machine-checkable verification with 11 assertions."""
    print("=" * 60)
    print("VERIFICATION MODE")
    print("=" * 60)
    passed = 0
    failed = 0

    def check(name, condition):
        nonlocal passed, failed
        status = "PASS" if condition else "FAIL"
        print(f"  [{status}] {name}")
        if condition:
            passed += 1
        else:
            failed += 1

    # 1. results.json exists and is valid JSON
    results = None
    try:
        with open(RESULTS_FILE) as f:
            results = json.load(f)
        check("results.json exists and is valid JSON", True)
    except Exception as e:
        check(f"results.json exists and is valid JSON ({e})", False)
        print(f"\nVERIFICATION FAILED: {failed} checks failed")
        sys.exit(1)

    # 2. Has comparisons with expected count
    comps = results.get("comparisons", [])
    check(f"results.json has >= 1 comparison (found {len(comps)})", len(comps) >= 1)

    # 3. All permutation p-values in [0, 1]
    pvals_ok = all(0 <= c["permutation_p_value"] <= 1 for c in comps)
    check("All permutation p-values in [0, 1]", pvals_ok)

    # 4. All KS statistics in [0, 1]
    ks_ok = all(0 <= c["ks_statistic"] <= 1 for c in comps)
    check("All KS statistics in [0, 1]", ks_ok)

    # 5. Cohen's d point estimate within bootstrap CI
    d_ok = all(c["cohens_d_ci_lo"] <= c["cohens_d"] <= c["cohens_d_ci_hi"] for c in comps)
    check("Cohen's d point estimates within CI bounds", d_ok)

    # 6. Median diff CI is ordered
    md_ok = all(c["median_diff_ci_lo"] <= c["median_diff_ci_hi"] for c in comps)
    check("Median difference CI is properly ordered", md_ok)

    # 7. Total entries analyzed > 1000
    total = results.get("n_total_entries_analyzed", 0)
    check(f"Total entries analyzed > 1000 (got {total})", total > 1000)

    # 8. Each organism has >= 30 entries
    desc = results.get("descriptive_statistics", {})
    for org_label in ["Human", "E. coli", "Yeast", "Mouse"]:
        if org_label in desc:
            n = desc[org_label]["n"]
            check(f"{org_label} has >= 30 entries (got {n})", n >= 30)

    # 9. Sensitivity analysis present and populated
    sens = results.get("sensitivity_analysis", {})
    check(
        "Sensitivity analysis has permutation count results",
        "permutation_count_sensitivity" in sens and len(sens["permutation_count_sensitivity"]) > 0,
    )
    check(
        "Sensitivity analysis has sample fraction results",
        "sample_fraction_sensitivity" in sens and len(sens["sample_fraction_sensitivity"]) > 0,
    )
    check(
        "Sensitivity analysis has threshold results",
        "resolution_threshold_sensitivity" in sens
        and len(sens["resolution_threshold_sensitivity"]) > 0,
    )

    # 10. report.md exists and has content
    report_ok = os.path.exists(REPORT_FILE) and os.path.getsize(REPORT_FILE) > 100
    check("report.md exists and has content", report_ok)

    # 11. Cache SHA256 integrity
    cache_ok = False
    try:
        with open(CACHE_FILE) as f:
            cache = json.load(f)
        verify_bytes = json.dumps(cache["entries"], sort_keys=True).encode()
        sha = hashlib.sha256(verify_bytes).hexdigest()
        cache_ok = sha == cache.get("sha256")
    except Exception:
        pass
    check("Cache file SHA256 integrity verified", cache_ok)

    # 12. Bonferroni correction correctly applied
    n_comp = results.get("n_comparisons", 3)
    bonf_ok = all(
        abs(c["bonferroni_p_value"] - min(c["permutation_p_value"] * n_comp, 1.0)) < 0.0001
        for c in comps
    )
    check("Bonferroni correction correctly applied", bonf_ok)

    print(f"\n  Results: {passed} passed, {failed} failed out of {passed + failed} checks")
    if failed > 0:
        print("VERIFICATION FAILED")
        sys.exit(1)
    else:
        print("VERIFICATION PASSED")


# ─── Main Entry Point ───────────────────────────────────────────────
def main():
    if "--verify" in sys.argv:
        verify_results()
        return

    t0 = time.time()

    print("[1/7] Fetching entry IDs from RCSB Search API...")
    entry_ids = fetch_entry_ids()

    print(f"\n[2/7] Fetching metadata for {len(entry_ids)} entries via GraphQL...")
    entries = fetch_metadata(entry_ids)
    print(f"  Total entries with metadata: {len(entries)}")

    print(f"\n[3/7] Organizing data by organism...")
    org_data = organize_data(entries)

    print(f"\n[4/7] Running KS tests with permutation p-values ({N_PERMUTATIONS} shuffles)...")
    comparisons = run_comparisons(org_data)

    print(f"\n[5/7] Running sensitivity analyses...")
    sensitivity = sensitivity_analysis(org_data)

    print(f"\n[6/7] Writing results...")
    write_results(org_data, comparisons, sensitivity)

    elapsed = time.time() - t0
    print(f"\n[7/7] Analysis complete in {elapsed:.1f} seconds")
    print("\nANALYSIS COMPLETE")


if __name__ == "__main__":
    main()
SCRIPT_EOF
```

**Expected output:** File `pdb_resolution_bias.py` written to workspace (exit code 0).

## Step 3: Run analysis

```bash
cd /tmp/claw4s_auto_pdb-resolution-bias-by-organism && python3 pdb_resolution_bias.py
```

**Expected output:**
- Sectioned progress `[1/7]` through `[7/7]`
- Descriptive statistics for 4 organisms
- KS statistics with permutation p-values for 3 comparisons
- Sensitivity analysis results
- Final line: `ANALYSIS COMPLETE`
- Files created: `results.json`, `report.md`, `pdb_metadata_cache.json`

**Expected runtime:** 5-30 minutes depending on RCSB API response times and permutation computation (5,000 permutations per comparison).

## Step 4: Verify results

```bash
cd /tmp/claw4s_auto_pdb-resolution-bias-by-organism && python3 pdb_resolution_bias.py --verify
```

**Expected output:**
- 12+ verification checks, each showing `[PASS]` or `[FAIL]`
- Final line: `VERIFICATION PASSED`
- Exit code 0

## Success Criteria

1. All 12+ verification checks pass
2. `results.json` contains valid statistical results for >= 3 organism comparisons
3. `report.md` contains descriptive statistics, comparison results, sensitivity analysis, and limitations
4. Cache file has valid SHA256 integrity hash
5. All p-values computed via permutation test (5,000 shuffles) with Welch's t-test as parametric complement
6. Bootstrap 95% CIs computed for both median difference and Cohen's d
7. Sensitivity analysis covers permutation count, sample fraction, and resolution threshold

## Failure Conditions

1. Any network error fetching from RCSB API (after retries exhausted)
2. Fewer than 1,000 entries successfully fetched
3. Any target organism has fewer than 30 entries
4. Verification reports any `[FAIL]` check
5. Script crashes or produces non-zero exit code

## Data Sources

- **RCSB Search API:** https://search.rcsb.org/rcsbsearch/v2/query
- **RCSB GraphQL API:** https://data.rcsb.org/graphql
- **RCSB REST API:** https://data.rcsb.org/rest/v1/core/entry/{entry_id}
- **Authority:** The RCSB PDB is the US member of the Worldwide Protein Data Bank (wwPDB),
  the authoritative archive for 3D macromolecular structure data.