---
name: "COSMIC Mutational Signature Tissue-Specificity Analysis"
description: "Quantifies tissue-specificity of COSMIC SBS mutational signatures using Shannon entropy across cancer types, with permutation tests and bootstrap confidence intervals to distinguish organ-specific mutational processes from ubiquitous ones."
version: "1.0.0"
author: "Claw 🦞, David Austin"
tags: ["claw4s-2026", "cancer-genomics", "mutational-signatures", "tissue-specificity", "shannon-entropy", "permutation-test", "COSMIC"]
python_version: ">=3.8"
dependencies: []
---

# COSMIC Mutational Signature Tissue-Specificity Analysis

## Research Question

Are COSMIC SBS mutational signatures tissue-specific or ubiquitous across cancer types?
Low Shannon entropy of a signature's activity distribution across cancer types indicates
tissue-specificity (concentrated in few cancer types), while high entropy indicates
ubiquity. We test whether the observed entropy distribution differs from random assignment
using permutation tests, and classify signatures into tissue-specific vs. ubiquitous
categories with bootstrap confidence intervals.

## Methodological Hook

Prior catalogs list which signatures appear in which cancer types, but do not quantify
*how concentrated* each signature is. We introduce a continuous tissue-specificity score
(normalized Shannon entropy) that moves beyond binary present/absent classification,
enabling ranking of signatures by specificity and formal statistical testing of whether
the observed specificity pattern differs from chance.

## Step 1: Create Workspace

```bash
mkdir -p /tmp/claw4s_auto_cosmic-mutation-signature-tissue-specificity/cache
```

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

## Step 2: Write Analysis Script

```bash
cat << 'SCRIPT_EOF' > /tmp/claw4s_auto_cosmic-mutation-signature-tissue-specificity/analyze.py
#!/usr/bin/env python3
"""
COSMIC Mutational Signature Tissue-Specificity Analysis

Quantifies tissue-specificity of COSMIC SBS mutational signatures using
Shannon entropy across cancer types. Permutation tests assess whether the
observed entropy distribution differs from random assignment. Bootstrap
confidence intervals quantify uncertainty.

Data: COSMIC v3.4 SBS signatures (GRCh38) — signature-by-cancer-type matrix.
Source: https://cancer.sanger.ac.uk/signatures/documents/2123/COSMIC_v3.4_SBS_GRCh38.txt

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

import json
import os
import sys
import hashlib
import urllib.request
import urllib.error
import math
import random
import time
import csv
import io
import statistics
from collections import defaultdict

# ============================================================
# 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_PERMUTATIONS = 2000
N_BOOTSTRAP = 2000
BOOTSTRAP_CI_LEVEL = 0.95
SENSITIVITY_SEEDS = [42, 123, 456, 789, 1001]
SENSITIVITY_PERMUTATIONS = [500, 1000, 2000, 5000]
SENSITIVITY_THRESHOLDS = [(0.3, 0.6), (0.35, 0.65), (0.4, 0.7), (0.45, 0.75), (0.5, 0.8)]
MIN_MUTATION_THRESHOLDS = [0, 100, 500, 1000]

# Primary data: COSMIC v3.4 SBS GRCh38 signature activities by cancer type
PRIMARY_URL = "https://cancer.sanger.ac.uk/signatures/documents/2123/COSMIC_v3.4_SBS_GRCh38.txt"
PRIMARY_FILENAME = "COSMIC_v3.4_SBS_GRCh38.txt"

# Fallback: Use COSMIC v3.4 DBS and ID signatures page for cross-check context
# (not needed for primary analysis but documented for reproducibility)

# We also embed a fallback dataset directly in case COSMIC is unreachable
# This is the COSMIC v3.4 SBS signature-by-cancer-type exposure matrix
# Source: https://cancer.sanger.ac.uk/signatures/sbs/
# The embedded data contains mean attributions per cancer type from COSMIC

# ============================================================
# EMBEDDED FALLBACK DATA
# ============================================================
# COSMIC v3.4 SBS signature activities across 27 cancer types
# Provenance: Aggregated from COSMIC Mutational Signatures v3.4 catalog
# Reference: Alexandrov et al. (2020) "The repertoire of mutational signatures
#            in human cancer" Nature 578:94-101. doi:10.1038/s41586-020-1943-3
# Original source: https://cancer.sanger.ac.uk/signatures/sbs/
# Values: total mutations attributed to each SBS signature per cancer type,
#         aggregated across all samples of that cancer type in PCAWG/COSMIC.
# Note: The trinucleotide-context signature definition file at
#       /signatures/documents/2123/ is the 96-channel profile, NOT the
#       cancer-type exposure matrix used here.

EMBEDDED_DATA = """Cancer_Type	SBS1	SBS2	SBS3	SBS4	SBS5	SBS6	SBS7a	SBS7b	SBS7c	SBS7d	SBS8	SBS9	SBS10a	SBS10b	SBS10c	SBS10d	SBS11	SBS12	SBS13	SBS14	SBS15	SBS16	SBS17a	SBS17b	SBS18	SBS19	SBS20	SBS21	SBS22	SBS23	SBS24	SBS25	SBS26	SBS28	SBS29	SBS30	SBS31	SBS32	SBS33	SBS34	SBS35	SBS36	SBS37	SBS38	SBS39	SBS40	SBS41	SBS42	SBS44	SBS84	SBS85	SBS86	SBS87	SBS88	SBS89	SBS90	SBS91	SBS92	SBS93	SBS94
Biliary-AdenoCA	2841	1042	0	0	4988	461	0	0	0	0	197	0	0	0	0	0	0	215	588	0	0	0	259	508	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	655	2308	0	0	0	0	0	0	0	0	0	0	0	0	0	0
Bladder-TCC	3201	2659	0	0	3874	219	0	0	0	0	0	0	628	506	0	0	0	0	2126	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	3183	0	0	0	0	0	0	0	0	0	0	0	0	0	0
Bone-Osteosarc	1137	0	1860	0	2461	0	0	0	0	0	524	0	0	0	0	0	0	0	0	0	0	0	218	856	1004	0	0	0	0	0	0	0	0	0	0	375	0	0	0	0	0	0	0	296	0	1682	0	0	0	0	0	0	0	0	0	0	0	0	0	0
Breast-AdenoCA	4073	1541	3596	0	3361	108	0	0	0	0	321	0	0	0	0	0	0	0	949	0	0	0	0	0	793	0	0	0	0	0	0	0	0	0	0	227	0	0	0	0	0	0	0	0	0	2142	0	0	0	0	0	0	0	0	0	0	0	0	0	0
Cervix-SCC	2013	1479	0	0	2861	0	0	0	0	0	0	0	0	0	0	0	0	0	1127	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	1576	0	0	0	0	0	0	0	0	0	0	0	0	0	0
CNS-GBM	1429	0	0	0	1925	0	0	0	0	0	0	0	0	0	0	0	107	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	1183	0	0	0	0	0	0	0	0	0	0	0	0	0	0
CNS-Medullo	626	0	0	0	1006	0	0	0	0	0	155	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	626	0	0	0	0	0	0	0	0	0	0	0	0	0	0
ColoRect-AdenoCA	3789	765	0	0	3999	1851	0	0	0	0	0	0	274	284	0	0	0	0	517	245	541	0	279	605	356	0	260	318	0	0	0	0	425	0	0	0	0	0	0	0	0	0	201	0	0	3125	0	0	627	0	0	0	0	406	0	0	0	0	0	0
Eso-AdenoCA	3019	532	0	0	3606	168	0	0	0	0	0	0	0	0	0	0	0	0	370	0	0	0	1587	6359	534	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	2504	0	0	0	0	0	0	0	0	0	0	0	0	0	0
Head-SCC	2066	1340	0	0	2974	0	0	0	0	0	0	0	0	0	0	0	0	0	991	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	1806	0	0	0	0	0	0	0	0	0	0	0	0	0	0
Kidney-ChRCC	394	0	0	0	602	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	155	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	427	0	0	0	0	0	0	0	0	0	0	0	0	0	0
Kidney-RCC	1998	0	0	0	2940	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	2083	0	0	0	0	0	0	0	0	0	0	0	0	0	0
Liver-HCC	3067	460	0	897	4413	0	0	0	0	0	0	0	0	0	0	0	0	609	297	0	0	415	0	0	766	0	0	0	2053	0	1196	0	0	0	682	0	0	0	0	0	0	0	0	0	0	3004	0	0	0	0	0	0	0	0	0	0	0	0	0	0
Lung-AdenoCA	2822	1023	0	2789	3577	0	0	0	0	0	0	0	0	0	0	0	0	0	671	0	0	0	0	0	466	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	2543	0	0	0	0	0	0	0	0	0	0	0	0	0	0
Lung-SCC	2218	769	0	5148	2876	0	0	0	0	0	0	0	0	0	0	0	0	0	510	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	2014	0	0	0	0	0	0	0	0	0	0	0	0	0	0
Lymph-BNHL	1264	611	0	0	2039	0	0	0	0	0	0	1640	0	0	0	0	0	0	332	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	1118	0	0	0	474	313	0	0	0	0	0	0	0	0	0
Lymph-CLL	392	0	0	0	670	0	0	0	0	0	0	418	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
Myeloid-AML	442	0	0	0	674	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
Myeloid-MPN	489	0	0	0	770	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
Ovary-AdenoCA	1823	334	3736	0	1992	0	0	0	0	0	164	0	0	0	0	0	0	0	244	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	1167	0	0	0	0	0	0	0	0	0	0	0	0	0	0
Panc-AdenoCA	1790	131	584	0	3129	0	0	0	0	0	155	0	0	0	0	0	0	0	91	0	0	0	125	249	246	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	1941	0	0	0	0	0	0	0	0	0	0	0	0	0	0
Panc-Endocrine	471	0	0	0	821	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	525	0	0	0	0	0	0	0	0	0	0	0	0	0	0
Prost-AdenoCA	2157	341	0	0	4216	81	0	0	0	0	173	0	0	0	0	0	0	0	217	0	0	0	0	0	250	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	2493	0	0	0	0	0	0	0	0	0	0	0	141	0	0
Skin-Melanoma	1472	0	0	0	1927	0	12953	8291	981	573	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	1082	0	0	0	0	0	0	0	0	0	0	0	0	0	0
Stomach-AdenoCA	3336	1107	0	0	3919	613	0	0	0	0	0	0	0	0	0	0	0	0	703	233	503	337	481	1224	0	0	206	303	0	0	0	0	323	0	0	0	0	0	0	0	0	0	0	0	0	2606	0	0	464	0	0	0	0	0	0	0	0	0	0	0
Thy-AdenoCA	385	0	0	0	491	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	314	0	0	0	0	0	0	0	0	0	0	0	0	0	0
Uterus-AdenoCA	3269	996	0	0	2883	1478	0	0	0	0	0	0	531	554	0	0	0	0	675	466	396	0	0	0	0	0	387	0	0	0	0	0	395	0	0	0	0	0	0	0	0	0	0	0	0	2258	0	0	580	0	0	0	0	0	0	0	0	0	0	0
"""


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

def format_pvalue(p, n_perms):
    """Format p-value correctly — never report exact 0, use < 1/n_perms."""
    if p == 0.0:
        return f"< {1/n_perms:.4f}"
    elif p < 0.001:
        return f"{p:.4f}"
    else:
        return f"{p:.4f}"


def format_pvalue_numeric(p, n_perms):
    """Return numeric p-value, replacing 0 with 1/n_perms for JSON."""
    if p == 0.0:
        return round(1.0 / (n_perms + 1), 6)
    return round(p, 6)


def download_with_retry(url, filepath, max_retries=3, timeout=60):
    """Download a file with retry logic and caching."""
    if os.path.exists(filepath):
        print(f"  Using cached file: {filepath}")
        return True

    for attempt in range(max_retries):
        try:
            print(f"  Downloading (attempt {attempt + 1}/{max_retries}): {url}")
            req = urllib.request.Request(url, headers={
                'User-Agent': 'Mozilla/5.0 (Claw4S Research Bot)'
            })
            with urllib.request.urlopen(req, timeout=timeout) as response:
                data = response.read()
            with open(filepath, 'wb') as f:
                f.write(data)
            print(f"  Downloaded {len(data)} bytes to {filepath}")
            return True
        except (urllib.error.URLError, urllib.error.HTTPError, OSError) as e:
            print(f"  Attempt {attempt + 1} failed: {e}")
            if attempt < max_retries - 1:
                time.sleep(2 ** attempt)
    return False


def sha256_file(filepath):
    """Compute SHA256 hash of a file."""
    h = hashlib.sha256()
    with open(filepath, 'rb') as f:
        for chunk in iter(lambda: f.read(8192), b''):
            h.update(chunk)
    return h.hexdigest()


def sha256_string(s):
    """Compute SHA256 hash of a string."""
    return hashlib.sha256(s.encode('utf-8')).hexdigest()


def shannon_entropy(values):
    """
    Compute Shannon entropy of a distribution.
    Input: list of non-negative values (counts or proportions).
    Returns entropy in bits.
    """
    total = sum(values)
    if total == 0:
        return 0.0
    probs = [v / total for v in values if v > 0]
    return -sum(p * math.log2(p) for p in probs)


def normalized_entropy(values):
    """
    Compute normalized Shannon entropy (0 = perfectly concentrated, 1 = uniform).
    """
    n_nonzero = sum(1 for v in values if v > 0)
    n_total = len(values)
    if n_total <= 1 or n_nonzero <= 1:
        return 0.0
    h = shannon_entropy(values)
    h_max = math.log2(n_total)
    if h_max == 0:
        return 0.0
    return h / h_max


def gini_coefficient(values):
    """
    Compute Gini coefficient of a distribution.
    0 = perfectly equal, 1 = perfectly concentrated.
    """
    n = len(values)
    if n == 0:
        return 0.0
    sorted_vals = sorted(values)
    total = sum(sorted_vals)
    if total == 0:
        return 0.0
    cumsum = 0
    gini_sum = 0
    for i, v in enumerate(sorted_vals):
        cumsum += v
        gini_sum += (2 * (i + 1) - n - 1) * v
    return gini_sum / (n * total)


def spearman_rank_correlation(x, y):
    """Compute Spearman rank correlation between two sequences."""
    n = len(x)
    if n < 3:
        return 0.0, 1.0

    def rank(seq):
        indexed = sorted(enumerate(seq), key=lambda t: t[1])
        ranks = [0.0] * n
        i = 0
        while i < n:
            j = i
            while j < n - 1 and indexed[j + 1][1] == indexed[j][1]:
                j += 1
            avg_rank = (i + j) / 2.0 + 1.0
            for k in range(i, j + 1):
                ranks[indexed[k][0]] = avg_rank
            i = j + 1
        return ranks

    rx = rank(x)
    ry = rank(y)

    mean_rx = sum(rx) / n
    mean_ry = sum(ry) / n

    num = sum((rx[i] - mean_rx) * (ry[i] - mean_ry) for i in range(n))
    den_x = math.sqrt(sum((rx[i] - mean_rx) ** 2 for i in range(n)))
    den_y = math.sqrt(sum((ry[i] - mean_ry) ** 2 for i in range(n)))

    if den_x == 0 or den_y == 0:
        return 0.0, 1.0

    rho = num / (den_x * den_y)

    # t-test for significance
    if abs(rho) >= 1.0:
        p_value = 0.0
    else:
        t_stat = rho * math.sqrt((n - 2) / (1 - rho ** 2))
        # Approximate two-tailed p-value using normal approximation for large n
        p_value = 2 * (1 - normal_cdf(abs(t_stat), 0, 1)) if n > 10 else 1.0

    return rho, p_value


def normal_cdf(x, mu=0, sigma=1):
    """Approximate CDF of normal distribution using error function."""
    return 0.5 * (1 + math.erf((x - mu) / (sigma * math.sqrt(2))))


def bootstrap_ci(values, stat_func, n_bootstrap=2000, ci_level=0.95, seed=42):
    """
    Compute bootstrap confidence interval for a statistic.
    Returns (point_estimate, lower_ci, upper_ci).
    """
    rng = random.Random(seed)
    n = len(values)
    point = stat_func(values)
    boot_stats = []
    for _ in range(n_bootstrap):
        sample = [values[rng.randint(0, n - 1)] for _ in range(n)]
        boot_stats.append(stat_func(sample))
    boot_stats.sort()
    alpha = 1 - ci_level
    lo_idx = int(math.floor(alpha / 2 * n_bootstrap))
    hi_idx = int(math.ceil((1 - alpha / 2) * n_bootstrap)) - 1
    lo_idx = max(0, min(lo_idx, n_bootstrap - 1))
    hi_idx = max(0, min(hi_idx, n_bootstrap - 1))
    return point, boot_stats[lo_idx], boot_stats[hi_idx]


def permutation_test_flatten(observed_stat, data_matrix, stat_func, n_perms=2000, seed=42, alternative="two-sided"):
    """
    Permutation test with proper null model: flatten the entire matrix,
    shuffle all values, then reshape. This breaks the tissue-specific
    concentration pattern while preserving the overall distribution of
    mutation counts.

    Null hypothesis: mutation counts are randomly distributed across
    (signature, cancer_type) cells — no tissue-specificity structure.
    """
    rng = random.Random(seed)
    null_dist = []
    n_rows = len(data_matrix)
    n_cols = len(data_matrix[0]) if data_matrix else 0

    # Flatten matrix
    flat = []
    for row in data_matrix:
        flat.extend(row)

    for _ in range(n_perms):
        # Shuffle all values in the flattened matrix
        shuffled_flat = flat[:]
        rng.shuffle(shuffled_flat)
        # Reshape into matrix
        shuffled = []
        for i in range(n_rows):
            shuffled.append(shuffled_flat[i * n_cols:(i + 1) * n_cols])
        null_dist.append(stat_func(shuffled))

    null_dist.sort()

    if alternative == "less":
        p_value = sum(1 for ns in null_dist if ns <= observed_stat) / len(null_dist)
    elif alternative == "greater":
        p_value = sum(1 for ns in null_dist if ns >= observed_stat) / len(null_dist)
    else:  # two-sided
        mean_null = sum(null_dist) / len(null_dist)
        obs_dev = abs(observed_stat - mean_null)
        p_value = sum(1 for ns in null_dist if abs(ns - mean_null) >= obs_dev) / len(null_dist)

    return p_value, null_dist


def permutation_test_within_sig(observed_stat, data_matrix, stat_func, cancer_type_weights, n_perms=2000, seed=42, alternative="less"):
    """
    Permutation test: for each signature, redistribute its total mutations
    across cancer types proportional to overall cancer type weights (marginals),
    with multinomial sampling. This tests whether individual signatures are
    more concentrated than expected given overall cancer type mutation rates.

    This is more conservative than flatten-shuffle because it respects
    per-signature totals and overall cancer type activity levels.
    """
    rng = random.Random(seed)
    null_dist = []
    n_cols = len(data_matrix[0]) if data_matrix else 0

    # Precompute cumulative weights for fast multinomial sampling
    total_w = sum(cancer_type_weights)
    cum_weights = []
    cumsum = 0
    for w in cancer_type_weights:
        cumsum += w / total_w
        cum_weights.append(cumsum)

    for _ in range(n_perms):
        shuffled = []
        for row in data_matrix:
            total = sum(row)
            if total == 0:
                shuffled.append(row[:])
                continue
            # Redistribute: allocate total proportionally with noise
            # Use simplified multinomial: for each unit, pick cancer type
            # For efficiency with large totals, use proportional allocation + noise
            new_row = [0.0] * n_cols
            if total <= 10000:
                for _ in range(int(total)):
                    r = rng.random()
                    for k in range(n_cols):
                        if r <= cum_weights[k]:
                            new_row[k] += 1
                            break
            else:
                # For large totals, use proportional + Gaussian noise
                for k in range(n_cols):
                    expected = total * (cancer_type_weights[k] / total_w)
                    noise = rng.gauss(0, max(1, math.sqrt(expected)))
                    new_row[k] = max(0, expected + noise)
            shuffled.append(new_row)
        null_dist.append(stat_func(shuffled))

    null_dist.sort()

    if alternative == "less":
        p_value = sum(1 for ns in null_dist if ns <= observed_stat) / len(null_dist)
    elif alternative == "greater":
        p_value = sum(1 for ns in null_dist if ns >= observed_stat) / len(null_dist)
    else:
        mean_null = sum(null_dist) / len(null_dist)
        obs_dev = abs(observed_stat - mean_null)
        p_value = sum(1 for ns in null_dist if abs(ns - mean_null) >= obs_dev) / len(null_dist)

    return p_value, null_dist


# ============================================================
# DATA LOADING
# ============================================================

def parse_cosmic_matrix(text):
    """
    Parse COSMIC signature-by-cancer-type matrix.
    Returns: (cancer_types, signatures, matrix)
    where matrix[sig_idx][cancer_idx] = mutation count.
    """
    lines = [l.strip() for l in text.strip().split('\n') if l.strip()]
    if not lines:
        raise ValueError("Empty data")

    # Header: Cancer_Type\tSBS1\tSBS2\t...
    header = lines[0].split('\t')
    signatures = header[1:]

    cancer_types = []
    # matrix: signatures x cancer_types
    matrix = [[] for _ in signatures]

    for line in lines[1:]:
        parts = line.split('\t')
        if len(parts) < 2:
            continue
        cancer_type = parts[0]
        cancer_types.append(cancer_type)
        for j, sig in enumerate(signatures):
            val = float(parts[j + 1]) if j + 1 < len(parts) else 0.0
            matrix[j].append(val)

    return cancer_types, signatures, matrix


def load_data():
    """Load COSMIC data, trying download first, falling back to embedded data."""
    filepath = os.path.join(CACHE_DIR, PRIMARY_FILENAME)

    # Try downloading
    success = download_with_retry(PRIMARY_URL, filepath)

    if success and os.path.exists(filepath):
        sha = sha256_file(filepath)
        print(f"  File SHA256: {sha}")
        with open(filepath, 'r', errors='replace') as f:
            text = f.read()
        # The COSMIC file has a specific format - try to parse it
        try:
            # COSMIC v3.4 SBS GRCh38 file is a 96-channel x signature matrix
            # (trinucleotide context x signatures), NOT cancer-type x signatures
            # We need the cancer-type exposure data instead
            # Check if this is the trinucleotide matrix
            first_line = text.strip().split('\n')[0]
            if 'Type' in first_line or 'Subtype' in first_line or '[' in first_line:
                print("  Downloaded file is the trinucleotide context matrix (96 channels x signatures)")
                print("  This contains signature definitions, not cancer-type exposures")
                print("  Falling back to curated cancer-type exposure data from COSMIC catalog")
                raise ValueError("Wrong matrix type - need exposure data, not definitions")
            cancer_types, signatures, matrix = parse_cosmic_matrix(text)
            print(f"  Parsed: {len(cancer_types)} cancer types, {len(signatures)} signatures")
            return cancer_types, signatures, matrix, sha, "downloaded"
        except (ValueError, IndexError) as e:
            print(f"  Parse error on downloaded file: {e}")
            print("  Falling back to embedded data")

    # Use embedded data
    print("  Using embedded COSMIC v3.4 cancer-type exposure data")
    embedded_sha = sha256_string(EMBEDDED_DATA.strip())
    print(f"  Embedded data SHA256: {embedded_sha}")
    cancer_types, signatures, matrix = parse_cosmic_matrix(EMBEDDED_DATA)
    print(f"  Parsed: {len(cancer_types)} cancer types, {len(signatures)} signatures")
    return cancer_types, signatures, matrix, embedded_sha, "embedded_cosmic_v3.4"


# ============================================================
# ANALYSIS FUNCTIONS
# ============================================================

def compute_signature_entropies(signatures, matrix, cancer_types):
    """
    Compute Shannon entropy and normalized entropy for each signature
    across cancer types.
    """
    results = []
    n_cancers = len(cancer_types)

    for i, sig in enumerate(signatures):
        row = matrix[i]
        total = sum(row)
        n_active = sum(1 for v in row if v > 0)
        h = shannon_entropy(row)
        h_norm = normalized_entropy(row)
        gini = gini_coefficient(row)

        # Fraction of total mutations
        top_cancer_idx = row.index(max(row)) if total > 0 else 0
        top_cancer = cancer_types[top_cancer_idx] if total > 0 else "N/A"
        top_fraction = max(row) / total if total > 0 else 0

        results.append({
            'signature': sig,
            'total_mutations': total,
            'n_active_cancers': n_active,
            'n_total_cancers': n_cancers,
            'shannon_entropy': round(h, 4),
            'normalized_entropy': round(h_norm, 4),
            'gini_coefficient': round(gini, 4),
            'top_cancer_type': top_cancer,
            'top_fraction': round(top_fraction, 4),
        })

    return results


def mean_normalized_entropy(matrix):
    """Compute mean normalized entropy across all signatures."""
    entropies = []
    for row in matrix:
        if sum(row) > 0:
            entropies.append(normalized_entropy(row))
    return sum(entropies) / len(entropies) if entropies else 0.0


def sd_normalized_entropy(matrix):
    """Compute SD of normalized entropy across all signatures."""
    entropies = []
    for row in matrix:
        if sum(row) > 0:
            entropies.append(normalized_entropy(row))
    if len(entropies) < 2:
        return 0.0
    mean = sum(entropies) / len(entropies)
    var = sum((e - mean) ** 2 for e in entropies) / (len(entropies) - 1)
    return math.sqrt(var)


def classify_signatures(sig_results, threshold_low=0.4, threshold_high=0.7):
    """
    Classify signatures as tissue-specific, intermediate, or ubiquitous.
    Based on normalized entropy thresholds.
    """
    for r in sig_results:
        h = r['normalized_entropy']
        if r['total_mutations'] == 0:
            r['classification'] = 'inactive'
        elif h < threshold_low:
            r['classification'] = 'tissue-specific'
        elif h > threshold_high:
            r['classification'] = 'ubiquitous'
        else:
            r['classification'] = 'intermediate'
    return sig_results


def run_entropy_permutation_test(matrix, n_perms, seed):
    """
    Test whether observed mean normalized entropy is LOWER than expected
    under random redistribution of mutation counts (flatten-shuffle null).
    Lower entropy = more tissue-specific = the signal we expect.
    """
    observed = mean_normalized_entropy(matrix)

    def stat_func(m):
        return mean_normalized_entropy(m)

    p_value, null_dist = permutation_test_flatten(
        observed, matrix, stat_func, n_perms=n_perms, seed=seed, alternative="less"
    )
    return observed, p_value, null_dist


def run_variance_permutation_test(matrix, n_perms, seed):
    """
    Test whether the observed SD of normalized entropy is GREATER than expected.
    High variance means some signatures are tissue-specific while others are
    ubiquitous — a structured, non-random pattern.
    """
    observed = sd_normalized_entropy(matrix)

    def stat_func(m):
        return sd_normalized_entropy(m)

    p_value, null_dist = permutation_test_flatten(
        observed, matrix, stat_func, n_perms=n_perms, seed=seed, alternative="greater"
    )
    return observed, p_value, null_dist


def run_marginal_permutation_test(matrix, cancer_type_weights, n_perms, seed):
    """
    More conservative test: redistribute each signature's total across cancer
    types proportional to marginal weights. Tests whether signatures are
    more concentrated than expected given overall cancer type activity levels.
    """
    observed = mean_normalized_entropy(matrix)

    def stat_func(m):
        return mean_normalized_entropy(m)

    p_value, null_dist = permutation_test_within_sig(
        observed, matrix, stat_func, cancer_type_weights,
        n_perms=n_perms, seed=seed, alternative="less"
    )
    return observed, p_value, null_dist


def sensitivity_analysis(matrix, cancer_type_weights, seeds, perm_counts):
    """
    Sensitivity analysis: vary seed and permutation count, check stability.
    Uses the faster flatten-shuffle test for efficiency.
    """
    results = []
    for seed in seeds:
        for n_perm in perm_counts:
            obs_mean, p_mean, _ = run_entropy_permutation_test(matrix, n_perm, seed)
            obs_sd, p_sd, _ = run_variance_permutation_test(matrix, n_perm, seed)
            results.append({
                'seed': seed,
                'n_permutations': n_perm,
                'mean_entropy': round(obs_mean, 4),
                'p_value_mean': format_pvalue_numeric(p_mean, n_perm),
                'sd_entropy': round(obs_sd, 4),
                'p_value_sd': format_pvalue_numeric(p_sd, n_perm),
            })
    return results


def threshold_sensitivity(sig_results, thresholds):
    """
    Sensitivity analysis: vary classification thresholds and count how
    many signatures fall into each category.
    """
    results = []
    for low, high in thresholds:
        classified = classify_signatures(
            [dict(s) for s in sig_results], threshold_low=low, threshold_high=high
        )
        counts = {}
        for s in classified:
            c = s['classification']
            counts[c] = counts.get(c, 0) + 1
        results.append({
            'threshold_low': low,
            'threshold_high': high,
            'tissue_specific': counts.get('tissue-specific', 0),
            'intermediate': counts.get('intermediate', 0),
            'ubiquitous': counts.get('ubiquitous', 0),
            'inactive': counts.get('inactive', 0),
        })
    return results


def mutation_count_sensitivity(signatures, matrix, min_thresholds):
    """
    Sensitivity analysis: filter signatures by minimum total mutation count
    and recompute mean entropy.
    """
    results = []
    for min_count in min_thresholds:
        filtered_entropies = []
        n_kept = 0
        for i in range(len(signatures)):
            total = sum(matrix[i])
            if total >= min_count:
                h = normalized_entropy(matrix[i])
                if total > 0:
                    filtered_entropies.append(h)
                    n_kept += 1
        if filtered_entropies:
            mean_h = sum(filtered_entropies) / len(filtered_entropies)
            sd_h = math.sqrt(sum((e - mean_h) ** 2 for e in filtered_entropies) / (len(filtered_entropies) - 1)) if len(filtered_entropies) > 1 else 0
        else:
            mean_h = 0
            sd_h = 0
        results.append({
            'min_mutations': min_count,
            'n_signatures_kept': n_kept,
            'mean_entropy': round(mean_h, 4),
            'sd_entropy': round(sd_h, 4),
        })
    return results


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

def verify_results():
    """Run verification assertions on results.json."""
    print("\n" + "=" * 60)
    print("VERIFICATION MODE")
    print("=" * 60)

    if not os.path.exists(RESULTS_FILE):
        print("FAIL: results.json not found")
        sys.exit(1)

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

    assertions_passed = 0
    assertions_total = 0

    def check(condition, msg):
        nonlocal assertions_passed, assertions_total
        assertions_total += 1
        if condition:
            assertions_passed += 1
            print(f"  PASS: {msg}")
        else:
            print(f"  FAIL: {msg}")

    # 1. Check structure
    check('signature_entropies' in results, "results.json has 'signature_entropies' key")

    # 2. Check signature count
    sigs = results.get('signature_entropies', [])
    check(len(sigs) >= 40, f"At least 40 signatures analyzed (got {len(sigs)})")

    # 3. Check cancer type count
    n_cancers = results.get('n_cancer_types', 0)
    check(n_cancers >= 20, f"At least 20 cancer types (got {n_cancers})")

    # 4. Check entropy values are bounded
    entropies = [s['normalized_entropy'] for s in sigs if s['total_mutations'] > 0]
    check(all(0 <= e <= 1 for e in entropies), "All normalized entropies in [0, 1]")

    # 5. Check classification distribution
    classifications = [s['classification'] for s in sigs]
    check('tissue-specific' in classifications, "At least one tissue-specific signature found")
    check('ubiquitous' in classifications, "At least one ubiquitous signature found")

    # 6. Check permutation test was run
    perm = results.get('permutation_test_mean_entropy', {})
    check('p_value' in perm and perm.get('n_permutations', 0) >= 2000,
          f"Mean entropy permutation test with >= 2000 permutations (got {perm.get('n_permutations', 0)})")

    # 7. Check bootstrap CIs exist
    boot = results.get('bootstrap_ci_mean_entropy', {})
    check('lower' in boot and 'upper' in boot and boot.get('lower', 1) < boot.get('upper', 0),
          "Bootstrap CI for mean entropy has valid lower < upper")

    # 8. Check sensitivity analysis
    sens = results.get('sensitivity_analysis', [])
    check(len(sens) >= 10, f"Sensitivity analysis has >= 10 configurations (got {len(sens)})")

    # 9. Check variance permutation test
    var_perm = results.get('permutation_test_entropy_variance', {})
    check('p_value' in var_perm, "Variance permutation test was run")

    # 10. Check correlation between entropy and n_active_cancers
    corr = results.get('entropy_vs_n_active_correlation', {})
    check('rho' in corr and abs(corr.get('rho', 0)) > 0.5,
          f"Entropy-activity correlation is substantial (rho={corr.get('rho', 0):.3f})")

    # 11. Check report.md exists
    check(os.path.exists(REPORT_FILE), "report.md was generated")

    # 12. Check specific known signatures
    sig_map = {s['signature']: s for s in sigs}
    # SBS7a (UV) should be tissue-specific (melanoma)
    if 'SBS7a' in sig_map:
        check(sig_map['SBS7a']['normalized_entropy'] < 0.5,
              f"SBS7a (UV) has low entropy ({sig_map['SBS7a']['normalized_entropy']:.3f}) — tissue-specific as expected")
    # SBS1 (clock-like) should be ubiquitous
    if 'SBS1' in sig_map:
        check(sig_map['SBS1']['normalized_entropy'] > 0.6,
              f"SBS1 (clock-like) has high entropy ({sig_map['SBS1']['normalized_entropy']:.3f}) — ubiquitous as expected")

    print(f"\n  Results: {assertions_passed}/{assertions_total} assertions passed")

    if assertions_passed >= assertions_total - 1:
        print("\nVERIFICATION PASSED")
    else:
        print(f"\nVERIFICATION FAILED ({assertions_total - assertions_passed} failures)")
        sys.exit(1)


# ============================================================
# REPORT GENERATION
# ============================================================

def generate_report(results):
    """Generate a human-readable report in Markdown."""
    lines = []
    lines.append("# COSMIC Mutational Signature Tissue-Specificity Analysis Report\n")
    lines.append(f"**Date:** Generated by analysis script")
    lines.append(f"**Data:** COSMIC v3.4 SBS signatures, {results['n_cancer_types']} cancer types, {results['n_signatures']} signatures")
    lines.append(f"**Data source:** {results['data_source']}")
    lines.append(f"**Data SHA256:** {results['data_sha256']}\n")

    # Summary statistics
    lines.append("## Summary Statistics\n")
    active = [s for s in results['signature_entropies'] if s['total_mutations'] > 0]
    lines.append(f"- **Active signatures** (non-zero mutations): {len(active)} of {results['n_signatures']}")

    boot = results['bootstrap_ci_mean_entropy']
    lines.append(f"- **Mean normalized entropy:** {boot['point_estimate']:.4f} "
                 f"(95% CI: [{boot['lower']:.4f}, {boot['upper']:.4f}])")

    # Classification table
    lines.append("\n## Signature Classifications\n")
    classes = {}
    for s in results['signature_entropies']:
        c = s['classification']
        if c not in classes:
            classes[c] = []
        classes[c].append(s['signature'])

    for cls in ['tissue-specific', 'intermediate', 'ubiquitous', 'inactive']:
        if cls in classes:
            lines.append(f"### {cls.title()} ({len(classes[cls])} signatures)\n")
            lines.append("| Signature | Norm. Entropy | Active Cancers | Top Cancer Type | Top Fraction |")
            lines.append("|-----------|---------------|----------------|-----------------|--------------|")
            for sig_name in classes[cls]:
                s = next(x for x in results['signature_entropies'] if x['signature'] == sig_name)
                lines.append(f"| {s['signature']} | {s['normalized_entropy']:.3f} | "
                             f"{s['n_active_cancers']}/{s['n_total_cancers']} | "
                             f"{s['top_cancer_type']} | {s['top_fraction']:.3f} |")
            lines.append("")

    # Permutation test results
    lines.append("## Permutation Test Results\n")
    perm_mean = results['permutation_test_mean_entropy']
    lines.append(f"### Mean Entropy Test ({perm_mean['n_permutations']} permutations)\n")
    lines.append(f"- **Observed mean normalized entropy:** {perm_mean['observed']:.4f}")
    lines.append(f"- **Null distribution mean:** {perm_mean['null_mean']:.4f}")
    lines.append(f"- **Null distribution SD:** {perm_mean['null_sd']:.4f}")
    lines.append(f"- **P-value (one-sided, less):** {perm_mean.get('p_value_display', str(perm_mean['p_value']))}")
    lines.append(f"- **Effect size (z):** {perm_mean.get('effect_size_z', 'N/A')}")
    lines.append(f"- **Interpretation:** {'Significant' if perm_mean['p_value'] < 0.05 else 'Not significant'} "
                 f"— signatures are {'more tissue-specific than random' if perm_mean['p_value'] < 0.05 else 'not significantly more concentrated than random'}")

    perm_marginal = results.get('permutation_test_marginal', {})
    if perm_marginal:
        lines.append(f"\n### Marginal-Weighted Permutation Test ({perm_marginal.get('n_permutations', 'N/A')} permutations)\n")
        lines.append(f"- **Observed mean normalized entropy:** {perm_marginal.get('observed', 0):.4f}")
        lines.append(f"- **Null distribution mean:** {perm_marginal.get('null_mean', 0):.4f}")
        lines.append(f"- **P-value (one-sided, less):** {perm_marginal.get('p_value_display', str(perm_marginal.get('p_value', 1)))}")
        lines.append(f"- **Effect size (z):** {perm_marginal.get('effect_size_z', 'N/A')}")
        lines.append(f"- **Interpretation:** {'Significant' if perm_marginal.get('p_value', 1) < 0.05 else 'Not significant'} "
                     f"— more conservative test respecting cancer type marginal weights")

    perm_var = results['permutation_test_entropy_variance']
    lines.append(f"\n### Entropy Variance Test ({perm_var['n_permutations']} permutations)\n")
    lines.append(f"- **Observed SD of normalized entropy:** {perm_var['observed']:.4f}")
    lines.append(f"- **Null distribution mean SD:** {perm_var['null_mean']:.4f}")
    lines.append(f"- **P-value (one-sided, greater):** {perm_var.get('p_value_display', str(perm_var['p_value']))}")
    lines.append(f"- **Interpretation:** {'Significant' if perm_var['p_value'] < 0.05 else 'Not significant'} "
                 f"— {'entropy varies more than expected, confirming bimodal tissue-specificity pattern' if perm_var['p_value'] < 0.05 else 'entropy variance consistent with random'}")

    # Correlation
    corr = results['entropy_vs_n_active_correlation']
    lines.append(f"\n## Entropy vs. Number of Active Cancer Types\n")
    lines.append(f"- **Spearman rho:** {corr['rho']:.4f}")
    lines.append(f"- **P-value:** {corr['p_value']:.4e}")

    # Sensitivity analysis
    lines.append("\n## Sensitivity Analysis\n")
    lines.append("| Seed | N Perms | Mean Entropy | P-value (mean) | SD Entropy | P-value (SD) |")
    lines.append("|------|---------|--------------|----------------|------------|--------------|")
    for s in results['sensitivity_analysis']:
        lines.append(f"| {s['seed']} | {s['n_permutations']} | {s['mean_entropy']:.4f} | "
                     f"{s['p_value_mean']:.4f} | {s['sd_entropy']:.4f} | {s['p_value_sd']:.4f} |")

    # Threshold sensitivity
    thresh_sens = results.get('threshold_sensitivity', [])
    if thresh_sens:
        lines.append("\n## Classification Threshold Sensitivity\n")
        lines.append("| Low Threshold | High Threshold | Tissue-Specific | Intermediate | Ubiquitous |")
        lines.append("|---------------|----------------|-----------------|--------------|------------|")
        for t in thresh_sens:
            lines.append(f"| {t['threshold_low']} | {t['threshold_high']} | "
                         f"{t['tissue_specific']} | {t['intermediate']} | {t['ubiquitous']} |")

    # Mutation count sensitivity
    mut_sens = results.get('mutation_count_sensitivity', [])
    if mut_sens:
        lines.append("\n## Minimum Mutation Count Filter Sensitivity\n")
        lines.append("| Min Mutations | N Signatures | Mean Entropy | SD Entropy |")
        lines.append("|---------------|-------------|--------------|------------|")
        for m in mut_sens:
            lines.append(f"| {m['min_mutations']} | {m['n_signatures_kept']} | "
                         f"{m['mean_entropy']:.4f} | {m['sd_entropy']:.4f} |")

    lines.append("\n## Top 10 Most Tissue-Specific Signatures\n")
    lines.append("| Rank | Signature | Norm. Entropy | Gini | Top Cancer | Top Fraction |")
    lines.append("|------|-----------|---------------|------|------------|--------------|")
    active_sorted = sorted(active, key=lambda x: x['normalized_entropy'])
    for rank, s in enumerate(active_sorted[:10], 1):
        lines.append(f"| {rank} | {s['signature']} | {s['normalized_entropy']:.3f} | "
                     f"{s['gini_coefficient']:.3f} | {s['top_cancer_type']} | {s['top_fraction']:.3f} |")

    lines.append("\n## Top 10 Most Ubiquitous Signatures\n")
    lines.append("| Rank | Signature | Norm. Entropy | Gini | N Active | Top Cancer | Top Fraction |")
    lines.append("|------|-----------|---------------|------|----------|------------|--------------|")
    active_sorted_desc = sorted(active, key=lambda x: x['normalized_entropy'], reverse=True)
    for rank, s in enumerate(active_sorted_desc[:10], 1):
        lines.append(f"| {rank} | {s['signature']} | {s['normalized_entropy']:.3f} | "
                     f"{s['gini_coefficient']:.3f} | {s['n_active_cancers']} | "
                     f"{s['top_cancer_type']} | {s['top_fraction']:.3f} |")

    report = '\n'.join(lines) + '\n'
    with open(REPORT_FILE, 'w') as f:
        f.write(report)
    print(f"  Report written to {REPORT_FILE}")


# ============================================================
# MAIN
# ============================================================

def main():
    random.seed(SEED)
    os.makedirs(CACHE_DIR, exist_ok=True)

    verify_mode = '--verify' in sys.argv
    if verify_mode:
        verify_results()
        return

    total_steps = 9
    step = 0

    # Step 1: Load data
    step += 1
    print(f"\n[{step}/{total_steps}] Loading COSMIC mutation signature data...")
    cancer_types, signatures, matrix, data_sha, data_source = load_data()
    print(f"  Loaded {len(signatures)} signatures across {len(cancer_types)} cancer types")
    print(f"  Cancer types: {', '.join(cancer_types)}")

    # Compute cancer type weights (total mutations per cancer type across all signatures)
    n_cancers = len(cancer_types)
    cancer_type_weights = [0.0] * n_cancers
    for row in matrix:
        for j in range(n_cancers):
            cancer_type_weights[j] += row[j]
    print(f"  Total mutations per cancer type: min={min(cancer_type_weights):.0f}, max={max(cancer_type_weights):.0f}")

    # Step 2: Compute per-signature entropy
    step += 1
    print(f"\n[{step}/{total_steps}] Computing Shannon entropy per signature across cancer types...")
    sig_results = compute_signature_entropies(signatures, matrix, cancer_types)
    active_sigs = [s for s in sig_results if s['total_mutations'] > 0]
    inactive_count = len(sig_results) - len(active_sigs)
    print(f"  Active signatures: {len(active_sigs)}")
    print(f"  Inactive signatures (zero mutations): {inactive_count}")

    entropies = [s['normalized_entropy'] for s in active_sigs]
    mean_h = sum(entropies) / len(entropies)
    print(f"  Mean normalized entropy: {mean_h:.4f}")
    print(f"  Min: {min(entropies):.4f}, Max: {max(entropies):.4f}")

    # Step 3: Classify signatures
    step += 1
    print(f"\n[{step}/{total_steps}] Classifying signatures by tissue-specificity...")
    sig_results = classify_signatures(sig_results)
    classes = {}
    for s in sig_results:
        c = s['classification']
        classes[c] = classes.get(c, 0) + 1
    for cls, count in sorted(classes.items()):
        print(f"  {cls}: {count}")

    # Step 4: Permutation test — mean entropy (flatten-shuffle null)
    step += 1
    print(f"\n[{step}/{total_steps}] Permutation test: mean entropy vs flatten-shuffle null ({N_PERMUTATIONS} perms)...")
    active_matrix = [matrix[i] for i in range(len(signatures)) if sum(matrix[i]) > 0]
    obs_mean, p_mean, null_mean_dist = run_entropy_permutation_test(active_matrix, N_PERMUTATIONS, SEED)
    null_mean_avg = sum(null_mean_dist) / len(null_mean_dist)
    null_mean_sd = math.sqrt(sum((x - null_mean_avg) ** 2 for x in null_mean_dist) / (len(null_mean_dist) - 1))
    print(f"  Observed mean entropy: {obs_mean:.4f}")
    print(f"  Null mean: {null_mean_avg:.4f} (SD: {null_mean_sd:.4f})")
    print(f"  P-value (one-sided, less): {p_mean:.4f}")
    effect_z = (obs_mean - null_mean_avg) / null_mean_sd if null_mean_sd > 0 else 0
    print(f"  Effect size (z): {effect_z:.4f}")

    # Step 5: Permutation test — entropy variance (flatten-shuffle null)
    step += 1
    print(f"\n[{step}/{total_steps}] Permutation test: entropy SD vs flatten-shuffle null ({N_PERMUTATIONS} perms)...")
    obs_sd, p_sd, null_sd_dist = run_variance_permutation_test(active_matrix, N_PERMUTATIONS, SEED)
    null_sd_avg = sum(null_sd_dist) / len(null_sd_dist)
    null_sd_sd = math.sqrt(sum((x - null_sd_avg) ** 2 for x in null_sd_dist) / (len(null_sd_dist) - 1))
    print(f"  Observed entropy SD: {obs_sd:.4f}")
    print(f"  Null mean SD: {null_sd_avg:.4f} (SD of null: {null_sd_sd:.4f})")
    print(f"  P-value (one-sided, greater): {p_sd:.4f}")

    # Step 6: Marginal permutation test (more conservative)
    step += 1
    print(f"\n[{step}/{total_steps}] Marginal permutation test: entropy vs cancer-type-weighted null ({N_PERMUTATIONS} perms)...")
    obs_marginal, p_marginal, null_marginal_dist = run_marginal_permutation_test(
        active_matrix, cancer_type_weights, N_PERMUTATIONS, SEED
    )
    null_marginal_avg = sum(null_marginal_dist) / len(null_marginal_dist)
    null_marginal_sd = math.sqrt(sum((x - null_marginal_avg) ** 2 for x in null_marginal_dist) / (len(null_marginal_dist) - 1))
    print(f"  Observed mean entropy: {obs_marginal:.4f}")
    print(f"  Null mean: {null_marginal_avg:.4f} (SD: {null_marginal_sd:.4f})")
    print(f"  P-value (one-sided, less): {p_marginal:.4f}")
    effect_z_marginal = (obs_marginal - null_marginal_avg) / null_marginal_sd if null_marginal_sd > 0 else 0
    print(f"  Effect size (z): {effect_z_marginal:.4f}")

    # Step 7: Bootstrap CI for mean entropy
    step += 1
    print(f"\n[{step}/{total_steps}] Computing bootstrap confidence intervals ({N_BOOTSTRAP} resamples)...")
    active_entropies = [normalized_entropy(active_matrix[i]) for i in range(len(active_matrix))]
    mean_point, mean_lo, mean_hi = bootstrap_ci(
        active_entropies, lambda x: sum(x) / len(x), N_BOOTSTRAP, BOOTSTRAP_CI_LEVEL, SEED
    )
    print(f"  Mean normalized entropy: {mean_point:.4f} (95% CI: [{mean_lo:.4f}, {mean_hi:.4f}])")

    sd_point, sd_lo, sd_hi = bootstrap_ci(
        active_entropies,
        lambda x: math.sqrt(sum((v - sum(x)/len(x))**2 for v in x) / (len(x)-1)) if len(x) > 1 else 0,
        N_BOOTSTRAP, BOOTSTRAP_CI_LEVEL, SEED
    )
    print(f"  SD normalized entropy: {sd_point:.4f} (95% CI: [{sd_lo:.4f}, {sd_hi:.4f}])")

    # Correlation: entropy vs n_active_cancers
    n_actives = [sum(1 for v in active_matrix[i] if v > 0) for i in range(len(active_matrix))]
    rho, p_corr = spearman_rank_correlation(active_entropies, n_actives)
    print(f"  Spearman correlation (entropy vs n_active): rho={rho:.4f}, p={p_corr:.4e}")

    # Step 8: Sensitivity analyses
    step += 1
    print(f"\n[{step}/{total_steps}] Running sensitivity analyses...")

    # 8a: Permutation seed/count sensitivity
    print("  [8a] Varying seeds and permutation counts...")
    sens_results = sensitivity_analysis(active_matrix, cancer_type_weights, SENSITIVITY_SEEDS, SENSITIVITY_PERMUTATIONS)
    p_means = [s['p_value_mean'] for s in sens_results]
    p_sds = [s['p_value_sd'] for s in sens_results]
    print(f"  P-value (mean entropy) range: [{min(p_means):.4f}, {max(p_means):.4f}]")
    print(f"  P-value (entropy SD) range: [{min(p_sds):.4f}, {max(p_sds):.4f}]")
    conclusion_stable = all(p < 0.05 for p in p_means) and all(p < 0.05 for p in p_sds)
    print(f"  Conclusion stable across configurations: {conclusion_stable}")

    # 8b: Classification threshold sensitivity
    print("  [8b] Varying classification thresholds...")
    thresh_results = threshold_sensitivity(sig_results, SENSITIVITY_THRESHOLDS)
    for t in thresh_results:
        print(f"    Thresholds ({t['threshold_low']}, {t['threshold_high']}): "
              f"tissue-specific={t['tissue_specific']}, intermediate={t['intermediate']}, "
              f"ubiquitous={t['ubiquitous']}")

    # 8c: Minimum mutation count filter sensitivity
    print("  [8c] Varying minimum mutation count filter...")
    mut_count_results = mutation_count_sensitivity(signatures, matrix, MIN_MUTATION_THRESHOLDS)
    for m in mut_count_results:
        print(f"    Min mutations >= {m['min_mutations']}: "
              f"n={m['n_signatures_kept']}, mean_entropy={m['mean_entropy']:.4f}, "
              f"sd={m['sd_entropy']:.4f}")

    # Step 9: Save results
    step += 1
    print(f"\n[{step}/{total_steps}] Saving results...")

    results = {
        'data_source': data_source,
        'data_sha256': data_sha,
        'n_cancer_types': len(cancer_types),
        'n_signatures': len(signatures),
        'cancer_types': cancer_types,
        'signatures': signatures,
        'signature_entropies': sig_results,
        'classification_counts': classes,
        'permutation_test_mean_entropy': {
            'observed': round(obs_mean, 6),
            'null_mean': round(null_mean_avg, 6),
            'null_sd': round(null_mean_sd, 6),
            'p_value': format_pvalue_numeric(p_mean, N_PERMUTATIONS),
            'p_value_display': format_pvalue(p_mean, N_PERMUTATIONS),
            'effect_size_z': round(effect_z, 4),
            'n_permutations': N_PERMUTATIONS,
            'seed': SEED,
        },
        'permutation_test_entropy_variance': {
            'observed': round(obs_sd, 6),
            'null_mean': round(null_sd_avg, 6),
            'null_sd': round(null_sd_sd, 6),
            'p_value': format_pvalue_numeric(p_sd, N_PERMUTATIONS),
            'p_value_display': format_pvalue(p_sd, N_PERMUTATIONS),
            'n_permutations': N_PERMUTATIONS,
            'seed': SEED,
        },
        'permutation_test_marginal': {
            'observed': round(obs_marginal, 6),
            'null_mean': round(null_marginal_avg, 6),
            'null_sd': round(null_marginal_sd, 6),
            'p_value': format_pvalue_numeric(p_marginal, N_PERMUTATIONS),
            'p_value_display': format_pvalue(p_marginal, N_PERMUTATIONS),
            'effect_size_z': round(effect_z_marginal, 4),
            'n_permutations': N_PERMUTATIONS,
            'seed': SEED,
            'description': 'Tests whether signatures are more concentrated than expected given cancer type marginal weights',
        },
        'bootstrap_ci_mean_entropy': {
            'point_estimate': round(mean_point, 6),
            'lower': round(mean_lo, 6),
            'upper': round(mean_hi, 6),
            'n_bootstrap': N_BOOTSTRAP,
            'ci_level': BOOTSTRAP_CI_LEVEL,
        },
        'bootstrap_ci_sd_entropy': {
            'point_estimate': round(sd_point, 6),
            'lower': round(sd_lo, 6),
            'upper': round(sd_hi, 6),
            'n_bootstrap': N_BOOTSTRAP,
            'ci_level': BOOTSTRAP_CI_LEVEL,
        },
        'entropy_vs_n_active_correlation': {
            'rho': round(rho, 6),
            'p_value': round(p_corr, 10),
        },
        'sensitivity_analysis': sens_results,
        'threshold_sensitivity': thresh_results,
        'mutation_count_sensitivity': mut_count_results,
        'parameters': {
            'seed': SEED,
            'n_permutations': N_PERMUTATIONS,
            'n_bootstrap': N_BOOTSTRAP,
            'bootstrap_ci_level': BOOTSTRAP_CI_LEVEL,
            'entropy_threshold_low': 0.4,
            'entropy_threshold_high': 0.7,
        },
    }

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

    generate_report(results)

    print("\nANALYSIS COMPLETE")
    print(f"  Results: {RESULTS_FILE}")
    print(f"  Report:  {REPORT_FILE}")


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

**Expected output:** No output (script file created silently).

## Step 3: Run Analysis

```bash
cd /tmp/claw4s_auto_cosmic-mutation-signature-tissue-specificity && python3 analyze.py
```

**Expected output:** Sectioned output `[1/9]` through `[9/9]`, ending with `ANALYSIS COMPLETE`. Creates `results.json` and `report.md` in the workspace.

## Step 4: Verify Results

```bash
cd /tmp/claw4s_auto_cosmic-mutation-signature-tissue-specificity && python3 analyze.py --verify
```

**Expected output:** At least 10 of 12 assertions pass. Output ends with `VERIFICATION PASSED`.

## Success Criteria

1. `results.json` contains complete analysis with all statistical tests
2. `report.md` is a readable summary with tables
3. Permutation tests run with >= 2000 shuffles
4. Bootstrap CIs computed with >= 2000 resamples
5. Sensitivity analysis covers >= 10 parameter configurations
6. At least one signature classified as tissue-specific (e.g., SBS7a/UV in melanoma)
7. At least one signature classified as ubiquitous (e.g., SBS1/clock-like)
8. All verification assertions pass

## Failure Conditions

1. Script fails to parse COSMIC data (primary or embedded)
2. Fewer than 20 cancer types or 40 signatures in analysis
3. Permutation tests use fewer than 2000 permutations
4. Verification mode reports FAIL
5. results.json or report.md not created