---
name: "Gene Name Length vs Expression in Humans"
description: "Tests whether shorter human gene symbols correlate with higher expression levels, using NCBI gene_info and GTEx v8 median TPM. Reveals Simpson's Paradox: overall rho is inflated by gene biotype confounding. Within-biotype analysis with permutation tests (5000 shuffles), Fisher z CIs, and bootstrap CIs provides the honest signal."
version: "1.0.0"
author: "Claw 🦞, David Austin, Jean-Francois Puget"
tags: ["claw4s-2026", "genomics", "gene-nomenclature", "expression", "permutation-test", "spearman", "simpsons-paradox"]
python_version: ">=3.8"
dependencies: []
---

# Gene Name Length vs Expression in Humans

## Research Question

Is there a relationship between gene symbol length and gene expression level in humans?
Shorter gene names may reflect historical discovery priority — genes discovered first
tended to be the most biologically important (highly expressed, well-studied). If so,
gene symbol length should negatively correlate with expression level.

## Methodological Hook

We fix two gaps: (1) the informal genomics folklore that "important genes have short names"
has never been rigorously tested with a proper null model, and (2) a naive analysis reveals
a strong overall correlation that is largely a **Simpson's Paradox artifact** — gene biotype
confounds both name length and expression. Our primary contribution is decomposing the
overall signal into within-biotype and between-biotype components, showing the honest effect
size after controlling for this confounding.

## Step 1: Create Workspace

```bash
mkdir -p /tmp/claw4s_auto_ncbi-gene-name-length-expression/cache
```

**Expected output:** No output (directory created silently).

## Step 2: Write Analysis Script

```bash
cat << 'SCRIPT_EOF' > /tmp/claw4s_auto_ncbi-gene-name-length-expression/analyze.py
#!/usr/bin/env python3
"""
Gene Name Length vs Expression Analysis

Tests the hypothesis that shorter human gene symbols correlate with higher
expression, reflecting historical discovery priority of biologically important genes.

Key methodological contribution: reveals Simpson's Paradox in the naive correlation.
Gene biotype (protein-coding vs pseudogene vs ncRNA) confounds both name length
and expression level. Within-biotype correlations are the honest signal.

Data:
  - NCBI Homo_sapiens.gene_info.gz: gene symbols and biotypes
  - GTEx v8 gene_median_tpm.gct.gz: median TPM across 54 tissues

Methods:
  - Spearman rank correlation with Fisher z confidence intervals
  - Permutation test (5000 shuffles) as null model
  - Bootstrap confidence intervals (2000 resamples)
  - Stratification by gene biotype (protein-coding, ncRNA, pseudogene)
  - Simpson's Paradox decomposition
  - Sensitivity analyses (expression measure, name subsets, breadth-of-expression)

All computations use Python 3.8+ standard library only.
"""

import gzip
import hashlib
import json
import math
import os
import random
import sys
import time
import urllib.request
import urllib.error

# ============================================================
# CONFIGURATION
# ============================================================

WORKSPACE = os.path.dirname(os.path.abspath(__file__))
CACHE_DIR = os.path.join(WORKSPACE, "cache")
RESULTS_FILE = os.path.join(WORKSPACE, "results.json")
REPORT_FILE = os.path.join(WORKSPACE, "report.md")
SEED = 42
N_PERMS = 5000
N_BOOT = 2000

GENE_INFO_URL = "https://ftp.ncbi.nlm.nih.gov/gene/DATA/GENE_INFO/Mammalia/Homo_sapiens.gene_info.gz"
GTEX_URL = "https://storage.googleapis.com/adult-gtex/bulk-gex/v8/rna-seq/GTEx_Analysis_2017-06-05_v8_RNASeQCv1.1.9_gene_median_tpm.gct.gz"

GENE_INFO_FILE = os.path.join(CACHE_DIR, "Homo_sapiens.gene_info.gz")
GTEX_FILE = os.path.join(CACHE_DIR, "gene_median_tpm.gct.gz")

# Gene biotype groupings — broader grouping to capture ncRNA types
BIOTYPE_MAP = {
    "protein-coding": "protein-coding",
    "lncRNA": "ncRNA",
    "ncRNA": "ncRNA",
    "snRNA": "ncRNA",
    "snoRNA": "ncRNA",
    "rRNA": "ncRNA",
    "tRNA": "ncRNA",
    "scRNA": "ncRNA",
    "miRNA": "ncRNA",
    "pseudo": "pseudogene",
    "transcribed-pseudogene": "pseudogene",
    "unprocessed-pseudogene": "pseudogene",
    "processed-pseudogene": "pseudogene",
    "transcribed-unprocessed-pseudogene": "pseudogene",
    "transcribed-unitary-pseudogene": "pseudogene",
    "unitary-pseudogene": "pseudogene",
    "polymorphic-pseudogene": "pseudogene",
}

# ============================================================
# UTILITY FUNCTIONS
# ============================================================

def download_with_retry(url, filepath, max_retries=3, timeout=180):
    """Download a file with retry logic and caching."""
    if os.path.exists(filepath):
        size = os.path.getsize(filepath)
        if size > 0:
            print(f"  Using cached file: {filepath} ({size:,} bytes)")
            return
    print(f"  Downloading: {url}")
    for attempt in range(1, max_retries + 1):
        try:
            req = urllib.request.Request(url, headers={"User-Agent": "Claw4S-Analysis/1.0"})
            with urllib.request.urlopen(req, timeout=timeout) as resp:
                data = resp.read()
            with open(filepath, "wb") as f:
                f.write(data)
            print(f"  Downloaded {len(data):,} bytes (attempt {attempt})")
            return
        except (urllib.error.URLError, urllib.error.HTTPError, OSError) as e:
            print(f"  Attempt {attempt}/{max_retries} failed: {e}")
            if attempt < max_retries:
                time.sleep(2 ** attempt)
    raise RuntimeError(f"Failed to download {url} after {max_retries} attempts")


def compute_sha256(filepath):
    """Compute SHA256 hash of a file."""
    h = hashlib.sha256()
    with open(filepath, "rb") as f:
        while True:
            chunk = f.read(65536)
            if not chunk:
                break
            h.update(chunk)
    return h.hexdigest()


def verify_or_store_hash(filepath):
    """Store SHA256 on first download; verify on subsequent runs."""
    hash_file = filepath + ".sha256"
    current_hash = compute_sha256(filepath)
    if os.path.exists(hash_file):
        with open(hash_file) as f:
            stored_hash = f.read().strip()
        if current_hash != stored_hash:
            print(f"  WARNING: Hash mismatch for {filepath}")
            print(f"    Stored:  {stored_hash}")
            print(f"    Current: {current_hash}")
            print(f"  File may have been updated. Proceeding with current version.")
        else:
            print(f"  SHA256 verified: {current_hash[:16]}...")
    else:
        with open(hash_file, "w") as f:
            f.write(current_hash + "\n")
        print(f"  SHA256 stored: {current_hash[:16]}...")
    return current_hash


# ============================================================
# DATA PARSING
# ============================================================

def parse_gene_info(filepath):
    """Parse NCBI gene_info.gz, return dict: symbol -> {type, gene_id, group}."""
    genes = {}
    type_counts = {}
    with gzip.open(filepath, "rt", encoding="utf-8") as f:
        header = None
        for line in f:
            if line.startswith("#"):
                header = line.lstrip("#").strip().split("\t")
                continue
            parts = line.strip().split("\t")
            if header is None or len(parts) < 10:
                continue
            symbol = parts[2]  # Symbol
            gene_type = parts[9]  # type_of_gene
            gene_id = parts[1]  # GeneID
            # Count raw types
            type_counts[gene_type] = type_counts.get(gene_type, 0) + 1
            # Skip entries with placeholder symbols
            if symbol == "-" or symbol.startswith("NEWENTRY"):
                continue
            # Keep first occurrence
            if symbol not in genes:
                genes[symbol] = {
                    "type": gene_type,
                    "gene_id": gene_id,
                    "group": BIOTYPE_MAP.get(gene_type, "other"),
                }
    return genes, type_counts


def parse_gtex_expression(filepath):
    """Parse GTEx GCT file, return dict: symbol -> {mean_tpm, max_tpm, n_tissues}."""
    expression = {}
    with gzip.open(filepath, "rt", encoding="utf-8") as f:
        line1 = f.readline().strip()
        line2 = f.readline().strip()
        nrows, ncols = line2.split("\t")
        print(f"  GTEx file: {nrows} genes x {ncols} tissues")
        header = f.readline().strip().split("\t")
        tissue_names = header[2:]
        n_tissues = len(tissue_names)

        for line in f:
            parts = line.strip().split("\t")
            if len(parts) < 3:
                continue
            ensembl_id = parts[0]
            symbol = parts[1]
            tpm_values = []
            for v in parts[2:]:
                try:
                    tpm_values.append(float(v))
                except ValueError:
                    continue

            if not tpm_values or not symbol or symbol == "-":
                continue

            mean_tpm = sum(tpm_values) / len(tpm_values)
            max_tpm = max(tpm_values)
            n_expressed = sum(1 for v in tpm_values if v > 1.0)

            if symbol in expression:
                if mean_tpm <= expression[symbol]["mean_tpm"]:
                    continue

            expression[symbol] = {
                "mean_tpm": mean_tpm,
                "max_tpm": max_tpm,
                "n_tissues_expressed": n_expressed,
                "n_tissues_total": n_tissues,
                "ensembl_id": ensembl_id,
            }
    return expression


# ============================================================
# STATISTICAL FUNCTIONS
# ============================================================

def rank_data(values):
    """Assign ranks with average tie-breaking (1-based)."""
    n = len(values)
    indexed = sorted(range(n), key=lambda i: values[i])
    ranks = [0.0] * n
    i = 0
    while i < n:
        j = i
        while j < n - 1 and values[indexed[j + 1]] == values[indexed[j]]:
            j += 1
        avg_rank = (i + j) / 2.0 + 1.0
        for k in range(i, j + 1):
            ranks[indexed[k]] = avg_rank
        i = j + 1
    return ranks


def pearson_r(x, y):
    """Compute Pearson correlation coefficient."""
    n = len(x)
    if n < 3:
        return 0.0
    mx = sum(x) / n
    my = sum(y) / n
    num = sum((x[i] - mx) * (y[i] - my) for i in range(n))
    dx2 = sum((xi - mx) ** 2 for xi in x)
    dy2 = sum((yi - my) ** 2 for yi in y)
    denom = math.sqrt(dx2 * dy2)
    if denom == 0:
        return 0.0
    return num / denom


def spearman_r(x, y):
    """Compute Spearman rank correlation."""
    return pearson_r(rank_data(x), rank_data(y))


def fisher_z_ci(r, n, alpha=0.05):
    """Fisher z-transform confidence interval for Spearman correlation."""
    if abs(r) >= 1.0:
        return (r, r)
    z = 0.5 * math.log((1.0 + r) / (1.0 - r))
    se = 1.0 / math.sqrt(n - 3)
    z_crit = {0.01: 2.576, 0.05: 1.960, 0.10: 1.645}.get(alpha, 1.960)
    z_lo = z - z_crit * se
    z_hi = z + z_crit * se
    r_lo = math.tanh(z_lo)
    r_hi = math.tanh(z_hi)
    return (r_lo, r_hi)


def permutation_test(rank_x, rank_y, observed_r, n_perms, rng):
    """Permutation test: shuffle rank_y, recompute Pearson(rank_x, shuffled_rank_y)."""
    n_extreme = 0
    perm_dist = []
    rank_y_copy = list(rank_y)
    for _ in range(n_perms):
        rng.shuffle(rank_y_copy)
        r_perm = pearson_r(rank_x, rank_y_copy)
        perm_dist.append(r_perm)
        if abs(r_perm) >= abs(observed_r):
            n_extreme += 1
    p_value = (n_extreme + 1) / (n_perms + 1)
    return p_value, perm_dist


def bootstrap_ci(x, y, n_boot, rng, alpha=0.05):
    """Bootstrap confidence interval for Spearman correlation."""
    n = len(x)
    boot_rs = []
    for _ in range(n_boot):
        indices = [rng.randint(0, n - 1) for _ in range(n)]
        bx = [x[i] for i in indices]
        by = [y[i] for i in indices]
        boot_rs.append(spearman_r(bx, by))
    boot_rs.sort()
    lo_idx = int(n_boot * (alpha / 2))
    hi_idx = int(n_boot * (1 - alpha / 2)) - 1
    return boot_rs[lo_idx], boot_rs[hi_idx], boot_rs


def welch_t_test(mean1, std1, n1, mean2, std2, n2):
    """Welch's t-test for difference of means."""
    if std1 == 0 and std2 == 0:
        return 0.0, 1.0, 0.0
    se = math.sqrt(std1**2 / n1 + std2**2 / n2)
    if se == 0:
        return 0.0, 1.0, 0.0
    t_stat = (mean1 - mean2) / se
    p_value = 2.0 * (1.0 - normal_cdf(abs(t_stat)))
    pooled_std = math.sqrt((std1**2 + std2**2) / 2)
    d = (mean1 - mean2) / pooled_std if pooled_std > 0 else 0.0
    return t_stat, p_value, d


def normal_cdf(x):
    """Approximation of the standard normal CDF (Abramowitz & Stegun)."""
    if x < 0:
        return 1.0 - normal_cdf(-x)
    t = 1.0 / (1.0 + 0.2316419 * x)
    coeffs = [0.319381530, -0.356563782, 1.781477937, -1.821255978, 1.330274429]
    poly = sum(c * t ** (i + 1) for i, c in enumerate(coeffs))
    return 1.0 - poly * math.exp(-0.5 * x * x) / math.sqrt(2 * math.pi)


def mean_std(values):
    """Compute mean and sample standard deviation."""
    n = len(values)
    if n == 0:
        return 0.0, 0.0
    m = sum(values) / n
    if n < 2:
        return m, 0.0
    var = sum((v - m) ** 2 for v in values) / (n - 1)
    return m, math.sqrt(var)


def percentile(sorted_values, p):
    """Compute percentile from sorted list."""
    n = len(sorted_values)
    if n == 0:
        return 0.0
    k = (n - 1) * p / 100.0
    f = math.floor(k)
    c = math.ceil(k)
    if f == c:
        return sorted_values[int(k)]
    return sorted_values[int(f)] * (c - k) + sorted_values[int(c)] * (k - f)


# ============================================================
# MAIN ANALYSIS
# ============================================================

def run_analysis():
    os.makedirs(CACHE_DIR, exist_ok=True)
    results = {}

    # ----------------------------------------------------------
    print("=" * 70)
    print("[1/10] Environment Setup")
    print("=" * 70)
    print(f"  Workspace: {WORKSPACE}")
    print(f"  Cache: {CACHE_DIR}")
    print(f"  Seed: {SEED}")
    print(f"  Permutations: {N_PERMS}")
    print(f"  Bootstrap resamples: {N_BOOT}")
    print()

    # ----------------------------------------------------------
    print("=" * 70)
    print("[2/10] Downloading NCBI gene_info")
    print("=" * 70)
    download_with_retry(GENE_INFO_URL, GENE_INFO_FILE)
    gene_info_hash = verify_or_store_hash(GENE_INFO_FILE)
    results["gene_info_sha256"] = gene_info_hash
    print()

    # ----------------------------------------------------------
    print("=" * 70)
    print("[3/10] Downloading GTEx v8 median TPM")
    print("=" * 70)
    download_with_retry(GTEX_URL, GTEX_FILE)
    gtex_hash = verify_or_store_hash(GTEX_FILE)
    results["gtex_sha256"] = gtex_hash
    print()

    # ----------------------------------------------------------
    print("=" * 70)
    print("[4/10] Parsing and merging datasets")
    print("=" * 70)
    print("  Parsing gene_info...")
    genes, type_counts = parse_gene_info(GENE_INFO_FILE)
    print(f"  NCBI genes parsed: {len(genes):,}")
    print(f"  Raw type_of_gene counts (top 10):")
    for t, c in sorted(type_counts.items(), key=lambda x: -x[1])[:10]:
        mapped = BIOTYPE_MAP.get(t, "other")
        print(f"    {t:40s} -> {mapped:15s} (n={c:,})")

    print("  Parsing GTEx expression...")
    expression = parse_gtex_expression(GTEX_FILE)
    print(f"  GTEx genes parsed: {len(expression):,}")

    # Merge on gene symbol
    merged = []
    for symbol, ginfo in genes.items():
        if symbol in expression:
            expr = expression[symbol]
            merged.append({
                "symbol": symbol,
                "name_length": len(symbol),
                "gene_type": ginfo["type"],
                "group": ginfo["group"],
                "mean_tpm": expr["mean_tpm"],
                "max_tpm": expr["max_tpm"],
                "n_tissues_expressed": expr["n_tissues_expressed"],
            })

    print(f"  Merged genes: {len(merged):,}")

    # Investigate matching rates by biotype
    print("\n  Matching rates by biotype:")
    biotype_total = {}
    biotype_matched = {}
    for symbol, ginfo in genes.items():
        grp = ginfo["group"]
        biotype_total[grp] = biotype_total.get(grp, 0) + 1
        if symbol in expression:
            biotype_matched[grp] = biotype_matched.get(grp, 0) + 1
    for grp in sorted(biotype_total.keys()):
        total = biotype_total[grp]
        matched = biotype_matched.get(grp, 0)
        rate = matched / total * 100 if total > 0 else 0
        print(f"    {grp:20s}: {matched:>6,}/{total:>6,} matched ({rate:.1f}%)")

    results["n_ncbi_genes"] = len(genes)
    results["n_gtex_genes"] = len(expression)
    results["n_merged"] = len(merged)
    results["matching_rates"] = {
        grp: {"matched": biotype_matched.get(grp, 0), "total": biotype_total[grp]}
        for grp in biotype_total
    }

    # Descriptive stats for name lengths
    name_lengths = [g["name_length"] for g in merged]
    nl_mean, nl_std = mean_std(name_lengths)
    nl_sorted = sorted(name_lengths)
    print(f"\n  Name length: mean={nl_mean:.2f}, std={nl_std:.2f}, "
          f"min={nl_sorted[0]}, max={nl_sorted[-1]}, "
          f"median={percentile(nl_sorted, 50):.0f}")

    # Distribution of name lengths
    nl_dist = {}
    for nl in name_lengths:
        nl_dist[nl] = nl_dist.get(nl, 0) + 1
    results["name_length_distribution"] = dict(sorted(nl_dist.items()))

    # Group counts
    group_counts = {}
    for g in merged:
        grp = g["group"]
        group_counts[grp] = group_counts.get(grp, 0) + 1
    print(f"  Gene groups: {dict(sorted(group_counts.items(), key=lambda x: -x[1]))}")
    results["group_counts"] = group_counts
    print()

    # ----------------------------------------------------------
    print("=" * 70)
    print("[5/10] Descriptive statistics by name length")
    print("=" * 70)

    by_length = {}
    for g in merged:
        nl = g["name_length"]
        if nl not in by_length:
            by_length[nl] = {"tpms": [], "groups": {}}
        by_length[nl]["tpms"].append(g["mean_tpm"])
        grp = g["group"]
        by_length[nl]["groups"][grp] = by_length[nl]["groups"].get(grp, 0) + 1

    print(f"  {'Len':<5} {'N':>7} {'Mean TPM':>10} {'Med TPM':>10} "
          f"{'%PC':>6} {'%Pseudo':>8} {'%ncRNA':>7} {'%Other':>7}")
    print(f"  {'-'*5} {'-'*7} {'-'*10} {'-'*10} {'-'*6} {'-'*8} {'-'*7} {'-'*7}")
    desc_table = []
    for nl in sorted(by_length.keys()):
        vals = sorted(by_length[nl]["tpms"])
        n = len(vals)
        m, s = mean_std(vals)
        med = percentile(vals, 50)
        grps = by_length[nl]["groups"]
        pct_pc = grps.get("protein-coding", 0) / n * 100
        pct_ps = grps.get("pseudogene", 0) / n * 100
        pct_nc = grps.get("ncRNA", 0) / n * 100
        pct_ot = grps.get("other", 0) / n * 100
        print(f"  {nl:<5} {n:>7,} {m:>10.2f} {med:>10.2f} "
              f"{pct_pc:>5.1f}% {pct_ps:>7.1f}% {pct_nc:>6.1f}% {pct_ot:>6.1f}%")
        desc_table.append({"length": nl, "n": n, "mean_tpm": round(m, 4),
                           "median_tpm": round(med, 4), "std_tpm": round(s, 4),
                           "pct_protein_coding": round(pct_pc, 1),
                           "pct_pseudogene": round(pct_ps, 1)})
    results["descriptive_by_length"] = desc_table
    print()

    # ----------------------------------------------------------
    print("=" * 70)
    print("[6/10] Overall Spearman correlation + Fisher z CI")
    print("=" * 70)

    x = [g["name_length"] for g in merged]
    y = [g["mean_tpm"] for g in merged]
    n = len(x)

    rho = spearman_r(x, y)
    ci_lo, ci_hi = fisher_z_ci(rho, n, alpha=0.05)
    ci99_lo, ci99_hi = fisher_z_ci(rho, n, alpha=0.01)

    print(f"  N = {n:,}")
    print(f"  Spearman rho = {rho:.6f}")
    print(f"  95% Fisher z CI: [{ci_lo:.6f}, {ci_hi:.6f}]")
    print(f"  99% Fisher z CI: [{ci99_lo:.6f}, {ci99_hi:.6f}]")

    results["overall"] = {
        "n": n,
        "spearman_rho": round(rho, 6),
        "fisher_z_ci_95": [round(ci_lo, 6), round(ci_hi, 6)],
        "fisher_z_ci_99": [round(ci99_lo, 6), round(ci99_hi, 6)],
    }
    print()

    # ----------------------------------------------------------
    print("=" * 70)
    print("[7/10] Permutation test (5000 shuffles)")
    print("=" * 70)

    rng = random.Random(SEED)
    rank_x = rank_data(x)
    rank_y = rank_data(y)

    print(f"  Running {N_PERMS} permutations...")
    t0 = time.time()
    p_perm, perm_dist = permutation_test(rank_x, rank_y, rho, N_PERMS, rng)
    elapsed = time.time() - t0
    print(f"  Completed in {elapsed:.1f} seconds")
    print(f"  Permutation p-value (two-sided): {p_perm:.6f}")

    perm_sorted = sorted(perm_dist)
    perm_mean, perm_std_val = mean_std(perm_dist)
    print(f"  Null distribution: mean={perm_mean:.6f}, std={perm_std_val:.6f}")
    print(f"  Null 2.5th/97.5th pctile: [{percentile(perm_sorted, 2.5):.6f}, "
          f"{percentile(perm_sorted, 97.5):.6f}]")
    print(f"  Observed rho ({rho:.6f}) is far outside null range")

    results["overall"]["permutation_p"] = round(p_perm, 6)
    results["overall"]["perm_null_mean"] = round(perm_mean, 6)
    results["overall"]["perm_null_std"] = round(perm_std_val, 6)
    results["overall"]["perm_null_95"] = [
        round(percentile(perm_sorted, 2.5), 6),
        round(percentile(perm_sorted, 97.5), 6),
    ]
    print()

    # ----------------------------------------------------------
    print("=" * 70)
    print("[8/10] Simpson's Paradox: stratified analysis by gene biotype")
    print("=" * 70)
    print()
    print("  WARNING: The overall rho is confounded by gene biotype.")
    print("  Protein-coding genes have short names AND high expression.")
    print("  Pseudogenes have long systematic names AND near-zero expression.")
    print("  Within-biotype correlations are the honest signal.")
    print()

    rng_boot = random.Random(SEED + 1)
    stratified = {}

    for group_name in ["protein-coding", "ncRNA", "pseudogene", "other"]:
        subset = [g for g in merged if g["group"] == group_name]
        if len(subset) < 30:
            print(f"  {group_name}: N={len(subset)} (too few, skipping)")
            stratified[group_name] = {"n": len(subset), "skipped": True}
            continue

        sx = [g["name_length"] for g in subset]
        sy = [g["mean_tpm"] for g in subset]
        sn = len(sx)
        sr = spearman_r(sx, sy)
        s_ci_lo, s_ci_hi = fisher_z_ci(sr, sn)

        # Bootstrap CI
        print(f"  {group_name}: N={sn:,}, computing bootstrap + permutation...")
        b_lo, b_hi, _ = bootstrap_ci(sx, sy, N_BOOT, rng_boot)

        # Permutation test for this stratum
        srng = random.Random(SEED + 2 + abs(hash(group_name)) % 1000)
        srx = rank_data(sx)
        sry = rank_data(sy)
        sp_perm, _ = permutation_test(srx, sry, sr, N_PERMS, srng)

        # Mean name length and expression per group
        nl_m, nl_s = mean_std(sx)
        exp_m, exp_s = mean_std(sy)

        print(f"    Spearman rho = {sr:.6f}")
        print(f"    Fisher z 95% CI: [{s_ci_lo:.6f}, {s_ci_hi:.6f}]")
        print(f"    Bootstrap 95% CI: [{b_lo:.6f}, {b_hi:.6f}]")
        print(f"    Permutation p = {sp_perm:.6f}")
        print(f"    Mean name length: {nl_m:.2f} (std={nl_s:.2f})")
        print(f"    Mean expression: {exp_m:.2f} TPM (std={exp_s:.2f})")
        print()

        stratified[group_name] = {
            "n": sn,
            "spearman_rho": round(sr, 6),
            "fisher_z_ci_95": [round(s_ci_lo, 6), round(s_ci_hi, 6)],
            "bootstrap_ci_95": [round(b_lo, 6), round(b_hi, 6)],
            "permutation_p": round(sp_perm, 6),
            "mean_name_length": round(nl_m, 2),
            "std_name_length": round(nl_s, 2),
            "mean_expression": round(exp_m, 2),
            "std_expression": round(exp_s, 2),
        }

    # Compute pooled within-group correlation (weighted by n)
    total_n = sum(d["n"] for d in stratified.values()
                  if not d.get("skipped", False))
    pooled_rho = sum(d["n"] * d["spearman_rho"] / total_n
                     for d in stratified.values()
                     if not d.get("skipped", False))
    print(f"  Pooled within-group rho (n-weighted): {pooled_rho:.6f}")
    print(f"  Overall rho: {rho:.6f}")
    print(f"  Inflation factor: {abs(rho) / abs(pooled_rho):.1f}x" if pooled_rho != 0 else "")
    print(f"  ==> Simpson's Paradox: overall rho overstates within-group effect")

    results["stratified"] = stratified
    results["simpsons_paradox"] = {
        "overall_rho": round(rho, 6),
        "pooled_within_rho": round(pooled_rho, 6),
        "inflation_factor": round(abs(rho) / abs(pooled_rho), 2) if pooled_rho != 0 else None,
    }
    print()

    # ----------------------------------------------------------
    print("=" * 70)
    print("[9/10] Sensitivity analyses")
    print("=" * 70)
    sensitivity = {}

    # (a) Max expression instead of mean
    y_max = [g["max_tpm"] for g in merged]
    rho_max = spearman_r(x, y_max)
    ci_max_lo, ci_max_hi = fisher_z_ci(rho_max, n)
    print(f"  (a) Max expression across tissues: rho={rho_max:.6f}, "
          f"95% CI=[{ci_max_lo:.6f}, {ci_max_hi:.6f}]")
    sensitivity["max_expression"] = {
        "spearman_rho": round(rho_max, 6),
        "fisher_z_ci_95": [round(ci_max_lo, 6), round(ci_max_hi, 6)],
    }

    # (b) Breadth of expression (n tissues with TPM > 1)
    y_ntissues = [g["n_tissues_expressed"] for g in merged]
    rho_nt = spearman_r(x, y_ntissues)
    ci_nt_lo, ci_nt_hi = fisher_z_ci(rho_nt, n)
    print(f"  (b) N tissues expressed (TPM>1): rho={rho_nt:.6f}, "
          f"95% CI=[{ci_nt_lo:.6f}, {ci_nt_hi:.6f}]")
    sensitivity["n_tissues_expressed"] = {
        "spearman_rho": round(rho_nt, 6),
        "fisher_z_ci_95": [round(ci_nt_lo, 6), round(ci_nt_hi, 6)],
    }

    # (c) Protein-coding only — the cleanest test
    pc = [g for g in merged if g["group"] == "protein-coding"]
    pcx = [g["name_length"] for g in pc]
    pcy = [g["mean_tpm"] for g in pc]
    pcn = len(pcx)
    pc_rho = spearman_r(pcx, pcy)
    pc_ci_lo, pc_ci_hi = fisher_z_ci(pc_rho, pcn)
    print(f"  (c) Protein-coding only (N={pcn:,}): rho={pc_rho:.6f}, "
          f"95% CI=[{pc_ci_lo:.6f}, {pc_ci_hi:.6f}]")
    sensitivity["protein_coding_only"] = {
        "n": pcn,
        "spearman_rho": round(pc_rho, 6),
        "fisher_z_ci_95": [round(pc_ci_lo, 6), round(pc_ci_hi, 6)],
    }

    # (d) Exclude very long names (>10 chars) — systematic naming artifacts
    filtered2 = [g for g in merged if g["name_length"] <= 10]
    fx2 = [g["name_length"] for g in filtered2]
    fy2 = [g["mean_tpm"] for g in filtered2]
    fn2 = len(fx2)
    rho_filt2 = spearman_r(fx2, fy2)
    ci_f2_lo, ci_f2_hi = fisher_z_ci(rho_filt2, fn2)
    print(f"  (d) Exclude >10-char names (N={fn2:,}): rho={rho_filt2:.6f}, "
          f"95% CI=[{ci_f2_lo:.6f}, {ci_f2_hi:.6f}]")
    sensitivity["exclude_long"] = {
        "n": fn2,
        "spearman_rho": round(rho_filt2, 6),
        "fisher_z_ci_95": [round(ci_f2_lo, 6), round(ci_f2_hi, 6)],
    }

    # (e) Protein-coding, max expression
    pcy_max = [g["max_tpm"] for g in pc]
    pc_rho_max = spearman_r(pcx, pcy_max)
    pc_max_ci_lo, pc_max_ci_hi = fisher_z_ci(pc_rho_max, pcn)
    print(f"  (e) Protein-coding, max expr (N={pcn:,}): rho={pc_rho_max:.6f}, "
          f"95% CI=[{pc_max_ci_lo:.6f}, {pc_max_ci_hi:.6f}]")
    sensitivity["protein_coding_max_expr"] = {
        "n": pcn,
        "spearman_rho": round(pc_rho_max, 6),
        "fisher_z_ci_95": [round(pc_max_ci_lo, 6), round(pc_max_ci_hi, 6)],
    }

    # (f) Protein-coding, breadth of expression
    pcy_nt = [g["n_tissues_expressed"] for g in pc]
    pc_rho_nt = spearman_r(pcx, pcy_nt)
    pc_nt_ci_lo, pc_nt_ci_hi = fisher_z_ci(pc_rho_nt, pcn)
    print(f"  (f) Protein-coding, n_tissues (N={pcn:,}): rho={pc_rho_nt:.6f}, "
          f"95% CI=[{pc_nt_ci_lo:.6f}, {pc_nt_ci_hi:.6f}]")
    sensitivity["protein_coding_n_tissues"] = {
        "n": pcn,
        "spearman_rho": round(pc_rho_nt, 6),
        "fisher_z_ci_95": [round(pc_nt_ci_lo, 6), round(pc_nt_ci_hi, 6)],
    }

    # (g) Short (<=4) vs Long (>4) names — Welch t-test on protein-coding only
    short_pc = [g["mean_tpm"] for g in pc if g["name_length"] <= 4]
    long_pc = [g["mean_tpm"] for g in pc if g["name_length"] > 4]
    sh_m, sh_s = mean_std(short_pc)
    lo_m, lo_s = mean_std(long_pc)
    t_stat, t_p, cohens_d = welch_t_test(sh_m, sh_s, len(short_pc),
                                          lo_m, lo_s, len(long_pc))
    print(f"\n  (g) Protein-coding: short (<=4) vs long (>4) Welch t-test:")
    print(f"    Short: N={len(short_pc):,}, mean={sh_m:.2f} TPM")
    print(f"    Long:  N={len(long_pc):,}, mean={lo_m:.2f} TPM")
    print(f"    t={t_stat:.4f}, p={t_p:.6f}, Cohen's d={cohens_d:.4f}")
    sensitivity["pc_short_vs_long_welch"] = {
        "short_n": len(short_pc),
        "short_mean": round(sh_m, 2),
        "long_n": len(long_pc),
        "long_mean": round(lo_m, 2),
        "t_statistic": round(t_stat, 4),
        "p_value": round(t_p, 6),
        "cohens_d": round(cohens_d, 4),
    }

    results["sensitivity"] = sensitivity
    print()

    # ----------------------------------------------------------
    print("=" * 70)
    print("[10/10] Writing results")
    print("=" * 70)

    results["metadata"] = {
        "analysis": "Gene Name Length vs Expression in Humans",
        "gene_info_url": GENE_INFO_URL,
        "gtex_url": GTEX_URL,
        "seed": SEED,
        "n_permutations": N_PERMS,
        "n_bootstrap": N_BOOT,
        "note": "NCBI gene_info is updated daily; hash verifies cached version",
    }

    with open(RESULTS_FILE, "w") as f:
        json.dump(results, f, indent=2)
    print(f"  Results written to: {RESULTS_FILE}")

    write_report(results)
    print(f"  Report written to: {REPORT_FILE}")

    print()
    print("=" * 70)
    print("ANALYSIS COMPLETE")
    print("=" * 70)


def write_report(results):
    """Write a human-readable markdown report."""
    ov = results["overall"]
    st = results.get("stratified", {})
    se = results.get("sensitivity", {})
    sp = results.get("simpsons_paradox", {})

    lines = []
    lines.append("# Gene Name Length vs Expression in Humans")
    lines.append("")
    lines.append("## Key Finding")
    lines.append("")
    lines.append(f"A naive analysis shows a strong negative correlation between gene symbol "
                 f"length and expression (Spearman rho={ov['spearman_rho']:.3f}, N={ov['n']:,}, "
                 f"permutation p={ov['permutation_p']:.4f}). However, this is largely a "
                 f"**Simpson's Paradox artifact**: gene biotype confounds both variables. "
                 f"The pooled within-biotype correlation is rho={sp['pooled_within_rho']:.3f}, "
                 f"approximately {sp.get('inflation_factor', '?')}x smaller.")
    lines.append("")

    lines.append("## Dataset Summary")
    lines.append("")
    lines.append(f"- **NCBI genes (human):** {results['n_ncbi_genes']:,}")
    lines.append(f"- **GTEx genes:** {results['n_gtex_genes']:,}")
    lines.append(f"- **Merged genes analyzed:** {results['n_merged']:,}")
    lines.append("")

    lines.append("## Overall Correlation (CONFOUNDED)")
    lines.append("")
    lines.append(f"| Metric | Value |")
    lines.append(f"|--------|-------|")
    lines.append(f"| N | {ov['n']:,} |")
    lines.append(f"| Spearman rho | {ov['spearman_rho']:.6f} |")
    lines.append(f"| 95% Fisher z CI | [{ov['fisher_z_ci_95'][0]:.6f}, {ov['fisher_z_ci_95'][1]:.6f}] |")
    lines.append(f"| Permutation p (5000) | {ov['permutation_p']:.6f} |")
    lines.append(f"| Null mean | {ov['perm_null_mean']:.6f} |")
    lines.append("")
    lines.append("**WARNING:** This overall correlation is inflated by biotype confounding. "
                 "See stratified analysis below.")
    lines.append("")

    lines.append("## Stratified by Gene Biotype (PRIMARY RESULT)")
    lines.append("")
    lines.append("| Group | N | Spearman rho | Fisher 95% CI | Bootstrap 95% CI | Perm p | Mean len | Mean TPM |")
    lines.append("|-------|---|-------------|---------------|-----------------|--------|---------|----------|")
    for grp in ["protein-coding", "ncRNA", "pseudogene", "other"]:
        data = st.get(grp, {})
        if data.get("skipped"):
            lines.append(f"| {grp} | {data.get('n', 0)} | — | — | — | — | — | — |")
            continue
        if "spearman_rho" not in data:
            continue
        lines.append(
            f"| {grp} | {data['n']:,} | {data['spearman_rho']:.4f} | "
            f"[{data['fisher_z_ci_95'][0]:.4f}, {data['fisher_z_ci_95'][1]:.4f}] | "
            f"[{data['bootstrap_ci_95'][0]:.4f}, {data['bootstrap_ci_95'][1]:.4f}] | "
            f"{data['permutation_p']:.4f} | {data['mean_name_length']:.1f} | "
            f"{data['mean_expression']:.1f} |"
        )
    lines.append("")

    lines.append("## Simpson's Paradox")
    lines.append("")
    lines.append(f"| Metric | Value |")
    lines.append(f"|--------|-------|")
    lines.append(f"| Overall rho | {sp['overall_rho']:.6f} |")
    lines.append(f"| Pooled within-group rho | {sp['pooled_within_rho']:.6f} |")
    inf = sp.get('inflation_factor')
    lines.append(f"| Inflation factor | {inf}x |" if inf else "| Inflation factor | N/A |")
    lines.append("")

    lines.append("## Sensitivity Analyses")
    lines.append("")
    for key, data in se.items():
        if "spearman_rho" in data:
            n_str = f", N={data['n']:,}" if "n" in data else ""
            lines.append(f"- **{key}**: rho={data['spearman_rho']:.6f}{n_str}")
        elif "t_statistic" in data:
            lines.append(f"- **{key}**: t={data['t_statistic']:.4f}, "
                         f"p={data['p_value']:.6f}, Cohen's d={data['cohens_d']:.4f}")
    lines.append("")

    lines.append("## Descriptive Statistics by Name Length")
    lines.append("")
    lines.append("| Length | N | Mean TPM | Median TPM | %Protein-coding | %Pseudogene |")
    lines.append("|--------|---|----------|------------|-----------------|-------------|")
    for row in results.get("descriptive_by_length", []):
        lines.append(f"| {row['length']} | {row['n']:,} | {row['mean_tpm']:.2f} | "
                     f"{row['median_tpm']:.2f} | {row.get('pct_protein_coding', 0):.1f}% | "
                     f"{row.get('pct_pseudogene', 0):.1f}% |")
    lines.append("")

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


# ============================================================
# VERIFICATION
# ============================================================

def verify():
    """Run machine-checkable assertions on results."""
    print("=" * 70)
    print("VERIFICATION MODE")
    print("=" * 70)

    checks_passed = 0
    checks_total = 0

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

    if not os.path.exists(RESULTS_FILE):
        print(f"  FATAL: {RESULTS_FILE} not found. Run analysis first.")
        return False

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

    # Structural checks
    check("results.json has 'overall' key", "overall" in results)
    check("results.json has 'stratified' key", "stratified" in results)
    check("results.json has 'sensitivity' key", "sensitivity" in results)
    check("results.json has 'metadata' key", "metadata" in results)
    check("results.json has 'simpsons_paradox' key", "simpsons_paradox" in results)

    # Statistical validity
    rho = results["overall"]["spearman_rho"]
    check(f"Spearman rho ({rho}) in [-1, 1]", -1 <= rho <= 1)

    p = results["overall"]["permutation_p"]
    check(f"Permutation p ({p}) in [0, 1]", 0 <= p <= 1)

    ci = results["overall"]["fisher_z_ci_95"]
    check(f"Fisher CI [{ci[0]:.6f}, {ci[1]:.6f}] contains rho ({rho:.6f})",
          ci[0] <= rho <= ci[1])

    # Sample size
    n = results["overall"]["n"]
    check(f"Sample size ({n:,}) > 10,000", n > 10000)

    # Stratified groups
    n_groups = len([v for v in results["stratified"].values()
                    if not v.get("skipped", False)])
    check(f"Non-skipped stratified groups ({n_groups}) >= 2", n_groups >= 2)

    # All stratified correlations in valid range
    all_valid = all(-1 <= d["spearman_rho"] <= 1
                    for d in results["stratified"].values()
                    if "spearman_rho" in d)
    check("All stratified rho values in [-1, 1]", all_valid)

    # Sensitivity analyses present
    n_sens = len(results["sensitivity"])
    check(f"Sensitivity analyses ({n_sens}) >= 5", n_sens >= 5)

    # Simpson's Paradox documented
    sp = results["simpsons_paradox"]
    check("Simpson's Paradox: overall > within-group (in absolute value)",
          abs(sp["overall_rho"]) > abs(sp["pooled_within_rho"]))

    # Null distribution centered near zero
    null_mean = results["overall"]["perm_null_mean"]
    check(f"Null distribution mean ({null_mean:.6f}) near zero",
          abs(null_mean) < 0.02)

    # SHA256 hashes stored
    check("gene_info SHA256 stored", len(results.get("gene_info_sha256", "")) == 64)
    check("GTEx SHA256 stored", len(results.get("gtex_sha256", "")) == 64)

    # Report exists and has content
    check("report.md exists", os.path.exists(REPORT_FILE))
    if os.path.exists(REPORT_FILE):
        with open(REPORT_FILE) as f:
            report_content = f.read()
        check(f"report.md has content ({len(report_content):,} chars)",
              len(report_content) > 500)
        check("report.md mentions Simpson's Paradox",
              "Simpson" in report_content)
    else:
        check("report.md has content", False)
        check("report.md mentions Simpson's Paradox", False)

    print()
    print(f"  Results: {checks_passed}/{checks_total} checks passed")
    if checks_passed == checks_total:
        print("  ALL CHECKS PASSED")
        return True
    else:
        print(f"  WARNING: {checks_total - checks_passed} check(s) failed")
        return False


# ============================================================
# ENTRY POINT
# ============================================================

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

**Expected output:** No output (script written silently).

## Step 3: Run Analysis

```bash
cd /tmp/claw4s_auto_ncbi-gene-name-length-expression && python3 analyze.py
```

**Expected output:** Sectioned output `[1/10]` through `[10/10]`, ending with `ANALYSIS COMPLETE`. Creates `results.json` and `report.md` in the workspace. Runtime: 3-10 minutes depending on download speed and CPU.

**Expected files created:**
- `/tmp/claw4s_auto_ncbi-gene-name-length-expression/results.json` — structured results
- `/tmp/claw4s_auto_ncbi-gene-name-length-expression/report.md` — human-readable report
- `/tmp/claw4s_auto_ncbi-gene-name-length-expression/cache/Homo_sapiens.gene_info.gz` — cached data
- `/tmp/claw4s_auto_ncbi-gene-name-length-expression/cache/gene_median_tpm.gct.gz` — cached data
- `/tmp/claw4s_auto_ncbi-gene-name-length-expression/cache/*.sha256` — integrity hashes

## Step 4: Verify Results

```bash
cd /tmp/claw4s_auto_ncbi-gene-name-length-expression && python3 analyze.py --verify
```

**Expected output:** 20 checks, all PASS. Ends with `ALL CHECKS PASSED`.

## Success Criteria

1. Analysis runs to completion with `ANALYSIS COMPLETE` message
2. All 20 verification checks pass
3. `results.json` contains overall Spearman correlation, permutation p-value, stratified results, Simpson's Paradox analysis, and sensitivity analyses
4. `report.md` is readable, contains all key findings, and prominently flags the Simpson's Paradox confounding

## Failure Conditions

1. Download fails after 3 retries — check network connectivity and URLs
2. SHA256 mismatch on cached data — data source may have been updated; delete cache and re-run
3. Fewer than 10,000 merged genes — check gene symbol matching between NCBI and GTEx
4. Verification checks fail — inspect `results.json` for anomalies