---
name: frameshield-viral-stop-benchmark
description: Reproduce FrameShield on a fixed panel of 19 viral RefSeq genomes. Fetches CDS and genome FASTA from NCBI, computes off-frame stop density, compares against amino-acid-preserving synonymous nulls, and verifies the published group-level signal.
allowed-tools: Bash(python3 *), Bash(curl *)
---

# FrameShield Viral Benchmark

## Overview

This skill reproduces the FrameShield result on a fixed `19`-accession viral genome panel.

Expected headline outputs:

- `19` genomes analyzed
- `13` positive z-scores
- `4` genomes with `z >= 2`
- high-overlap group: `7/8` positive, `4/8` at `z >= 2`
- large-DNA group: `0/3` positive
- verification marker: `frameshield_benchmark_verified`

Expected runtime: about 1-3 minutes depending on NCBI response speed.

## Step 1: Create a Clean Workspace

```bash
mkdir -p frameshield_repro/scripts
cd frameshield_repro
```

Expected output: no terminal output.

## Step 2: Write the Reference Benchmark Script

```bash
cat > scripts/frameshield_benchmark.py <<'PY'
#!/usr/bin/env python3
import argparse
import json
import math
import pathlib
import random
import re
import statistics
import time
import urllib.request
from collections import Counter, defaultdict
from typing import Dict, List, Optional, Sequence, Tuple


ACCESSIONS = [
    {"name": "SARS-CoV-2", "accession": "NC_045512.2", "group": "high-overlap"},
    {"name": "MERS-CoV", "accession": "NC_019843.3", "group": "high-overlap"},
    {"name": "SARS-CoV", "accession": "NC_004718.3", "group": "high-overlap"},
    {"name": "HCoV-OC43", "accession": "NC_006213.1", "group": "high-overlap"},
    {"name": "HCoV-NL63", "accession": "NC_005831.2", "group": "high-overlap"},
    {"name": "HBV", "accession": "NC_003977.2", "group": "high-overlap"},
    {"name": "HIV-1", "accession": "NC_001802.1", "group": "high-overlap"},
    {"name": "HTLV-1", "accession": "NC_001436.1", "group": "high-overlap"},
    {"name": "Dengue-2", "accession": "NC_001474.2", "group": "low-overlap"},
    {"name": "Zika", "accession": "NC_012532.1", "group": "low-overlap"},
    {"name": "HCV", "accession": "NC_004102.1", "group": "low-overlap"},
    {"name": "Chikungunya", "accession": "NC_004162.2", "group": "low-overlap"},
    {"name": "Poliovirus-1", "accession": "NC_002058.3", "group": "low-overlap"},
    {"name": "Rabies", "accession": "NC_001542.1", "group": "low-overlap"},
    {"name": "Measles", "accession": "NC_001498.1", "group": "low-overlap"},
    {"name": "Ebola", "accession": "NC_002549.1", "group": "low-overlap"},
    {"name": "Adenovirus-C", "accession": "NC_001405.1", "group": "large-dna"},
    {"name": "HSV-1", "accession": "NC_001806.2", "group": "large-dna"},
    {"name": "Vaccinia", "accession": "NC_006998.1", "group": "large-dna"},
]

CODON_TO_AA = {
    "TTT": "F",
    "TTC": "F",
    "TTA": "L",
    "TTG": "L",
    "CTT": "L",
    "CTC": "L",
    "CTA": "L",
    "CTG": "L",
    "ATT": "I",
    "ATC": "I",
    "ATA": "I",
    "ATG": "M",
    "GTT": "V",
    "GTC": "V",
    "GTA": "V",
    "GTG": "V",
    "TCT": "S",
    "TCC": "S",
    "TCA": "S",
    "TCG": "S",
    "CCT": "P",
    "CCC": "P",
    "CCA": "P",
    "CCG": "P",
    "ACT": "T",
    "ACC": "T",
    "ACA": "T",
    "ACG": "T",
    "GCT": "A",
    "GCC": "A",
    "GCA": "A",
    "GCG": "A",
    "TAT": "Y",
    "TAC": "Y",
    "TAA": "*",
    "TAG": "*",
    "CAT": "H",
    "CAC": "H",
    "CAA": "Q",
    "CAG": "Q",
    "AAT": "N",
    "AAC": "N",
    "AAA": "K",
    "AAG": "K",
    "GAT": "D",
    "GAC": "D",
    "GAA": "E",
    "GAG": "E",
    "TGT": "C",
    "TGC": "C",
    "TGA": "*",
    "TGG": "W",
    "CGT": "R",
    "CGC": "R",
    "CGA": "R",
    "CGG": "R",
    "AGT": "S",
    "AGC": "S",
    "AGA": "R",
    "AGG": "R",
    "GGT": "G",
    "GGC": "G",
    "GGA": "G",
    "GGG": "G",
}
AA_TO_CODONS: Dict[str, List[str]] = defaultdict(list)
for codon, aa in CODON_TO_AA.items():
    if aa != "*":
        AA_TO_CODONS[aa].append(codon)

STOP_CODONS = {"TAA", "TAG", "TGA"}
LOCATION_RE = re.compile(r"\[location=([^\]]+)\]")


def fetch_text(url: str, cache_path: Optional[pathlib.Path] = None) -> str:
    if cache_path is not None and cache_path.exists():
        return cache_path.read_text()

    last_error: Optional[Exception] = None
    for attempt in range(5):
        try:
            with urllib.request.urlopen(url, timeout=120) as response:
                payload = response.read().decode("utf-8")
            if cache_path is not None:
                cache_path.parent.mkdir(parents=True, exist_ok=True)
                cache_path.write_text(payload)
            time.sleep(0.4)
            return payload
        except Exception as exc:
            last_error = exc
            time.sleep(1.5 * (attempt + 1))
    raise RuntimeError(f"Failed to fetch {url}") from last_error


def parse_fasta(text: str) -> List[Tuple[str, str]]:
    records: List[Tuple[str, str]] = []
    header = None
    seq_parts: List[str] = []
    for line in text.splitlines():
        if not line:
            continue
        if line.startswith(">"):
            if header is not None:
                records.append((header, "".join(seq_parts)))
            header = line[1:]
            seq_parts = []
        else:
            seq_parts.append(line.strip())
    if header is not None:
        records.append((header, "".join(seq_parts)))
    return records


def parse_intervals(header: str) -> List[Tuple[int, int]]:
    match = LOCATION_RE.search(header)
    if not match:
        return []
    intervals = []
    for start, end in re.findall(r"(\d+)\.\.(\d+)", match.group(1)):
        a = int(start)
        b = int(end)
        intervals.append((min(a, b), max(a, b)))
    return intervals


def translate(seq: str) -> str:
    residues = []
    for idx in range(0, len(seq), 3):
        residues.append(CODON_TO_AA.get(seq[idx : idx + 3], "X"))
    return "".join(residues)


def off_frame_stop_density(sequences: Sequence[str]) -> Dict[str, float]:
    stop_count = 0
    triplet_count = 0
    for seq in sequences:
        for shift in (1, 2):
            for idx in range(shift, len(seq) - 2, 3):
                triplet_count += 1
                if seq[idx : idx + 3] in STOP_CODONS:
                    stop_count += 1
    density = stop_count / triplet_count if triplet_count else 0.0
    return {"stop_count": stop_count, "triplet_count": triplet_count, "density": density}


def ranks(values: Sequence[float]) -> List[float]:
    ordered = sorted((value, idx) for idx, value in enumerate(values))
    ranked = [0.0] * len(values)
    i = 0
    while i < len(ordered):
        j = i
        while j + 1 < len(ordered) and ordered[j + 1][0] == ordered[i][0]:
            j += 1
        rank = (i + j + 2) / 2.0
        for _, idx in ordered[i : j + 1]:
            ranked[idx] = rank
        i = j + 1
    return ranked


def spearman(values_x: Sequence[float], values_y: Sequence[float]) -> float:
    ranked_x = ranks(values_x)
    ranked_y = ranks(values_y)
    mean_x = statistics.mean(ranked_x)
    mean_y = statistics.mean(ranked_y)
    numerator = sum((x - mean_x) * (y - mean_y) for x, y in zip(ranked_x, ranked_y))
    denominator = math.sqrt(
        sum((x - mean_x) ** 2 for x in ranked_x) * sum((y - mean_y) ** 2 for y in ranked_y)
    )
    return numerator / denominator if denominator else 0.0


def gene_overlap_fraction(genome_length: int, intervals: Sequence[Tuple[int, int]]) -> float:
    if genome_length <= 0:
        return 0.0
    coverage = [0] * (genome_length + 1)
    for start, end in intervals:
        start = max(1, start)
        end = min(genome_length, end)
        for pos in range(start, end + 1):
            coverage[pos] += 1
    coding_bp = sum(1 for depth in coverage[1:] if depth >= 1)
    overlap_bp = sum(1 for depth in coverage[1:] if depth >= 2)
    return overlap_bp / coding_bp if coding_bp else 0.0


def collect_cds_records(accession: str, cache_dir: pathlib.Path) -> Tuple[List[Dict[str, object]], int]:
    cds_url = (
        "https://eutils.ncbi.nlm.nih.gov/entrez/eutils/"
        f"efetch.fcgi?db=nuccore&id={accession}&rettype=fasta_cds_na&retmode=text"
    )
    genome_url = (
        "https://eutils.ncbi.nlm.nih.gov/entrez/eutils/"
        f"efetch.fcgi?db=nuccore&id={accession}&rettype=fasta&retmode=text"
    )

    cds_records = []
    cds_cache = cache_dir / f"{accession}_cds.fasta"
    genome_cache = cache_dir / f"{accession}_genome.fasta"
    for header, raw_seq in parse_fasta(fetch_text(cds_url, cds_cache)):
        seq = raw_seq.upper().replace("U", "T")
        if set(seq) - set("ACGT"):
            continue
        if len(seq) % 3 != 0:
            continue
        core_seq = seq[:-3] if seq[-3:] in STOP_CODONS else seq
        if len(core_seq) < 30 or len(core_seq) % 3 != 0:
            continue
        aa_seq = translate(core_seq)
        if "*" in aa_seq or "X" in aa_seq:
            continue
        cds_records.append(
            {
                "header": header,
                "nt_sequence": core_seq,
                "aa_sequence": aa_seq,
                "intervals": parse_intervals(header),
            }
        )

    genome_records = parse_fasta(fetch_text(genome_url, genome_cache))
    genome_length = len(genome_records[0][1]) if genome_records else 0
    return cds_records, genome_length


def build_synonymous_sampler(cds_records: Sequence[Dict[str, object]]) -> Dict[str, Tuple[List[str], List[int]]]:
    codon_counts: Counter[str] = Counter()
    for record in cds_records:
        nt_sequence = str(record["nt_sequence"])
        for idx in range(0, len(nt_sequence), 3):
            codon_counts[nt_sequence[idx : idx + 3]] += 1

    sampler: Dict[str, Tuple[List[str], List[int]]] = {}
    for aa, codons in AA_TO_CODONS.items():
        weights = [codon_counts[codon] for codon in codons]
        sampler[aa] = (codons, weights if sum(weights) else [1] * len(codons))
    return sampler


def sample_synonymous_sequences(
    cds_records: Sequence[Dict[str, object]],
    sampler: Dict[str, Tuple[List[str], List[int]]],
    rng: random.Random,
) -> List[str]:
    sampled_sequences = []
    for record in cds_records:
        aa_sequence = str(record["aa_sequence"])
        original_nt = str(record["nt_sequence"])
        codons = []
        for idx, aa in enumerate(aa_sequence):
            choices, weights = sampler[aa]
            if idx == 0 and original_nt[:3] in choices:
                codons.append(original_nt[:3])
            else:
                codons.append(rng.choices(choices, weights=weights, k=1)[0])
        sampled_sequences.append("".join(codons))
    return sampled_sequences


def run_benchmark(outdir: pathlib.Path, simulations: int, seed: int) -> Dict[str, object]:
    rng = random.Random(seed)
    cache_dir = outdir / "cache"
    viruses = []
    for spec in ACCESSIONS:
        cds_records, genome_length = collect_cds_records(spec["accession"], cache_dir)
        observed = off_frame_stop_density([str(record["nt_sequence"]) for record in cds_records])
        sampler = build_synonymous_sampler(cds_records)

        simulated_densities = []
        for _ in range(simulations):
            simulated_sequences = sample_synonymous_sequences(cds_records, sampler, rng)
            simulated_densities.append(off_frame_stop_density(simulated_sequences)["density"])
        null_mean = statistics.mean(simulated_densities)
        null_sd = statistics.pstdev(simulated_densities)
        z_score = (observed["density"] - null_mean) / null_sd if null_sd else 0.0

        overlap_fraction = gene_overlap_fraction(
            genome_length,
            [interval for record in cds_records for interval in record["intervals"]],
        )
        viruses.append(
            {
                "name": spec["name"],
                "accession": spec["accession"],
                "group": spec["group"],
                "cds_count": len(cds_records),
                "genome_length": genome_length,
                "off_frame_stop_count": observed["stop_count"],
                "off_frame_triplet_count": observed["triplet_count"],
                "observed_density": observed["density"],
                "null_mean_density": null_mean,
                "null_sd_density": null_sd,
                "z_score": z_score,
                "overlap_fraction": overlap_fraction,
            }
        )

    group_summary = {}
    for group in sorted({virus["group"] for virus in viruses}):
        group_viruses = [virus for virus in viruses if virus["group"] == group]
        group_summary[group] = {
            "count": len(group_viruses),
            "median_z_score": statistics.median(virus["z_score"] for virus in group_viruses),
            "median_overlap_fraction": statistics.median(
                virus["overlap_fraction"] for virus in group_viruses
            ),
            "positive_z_count": sum(1 for virus in group_viruses if virus["z_score"] > 0),
            "z_at_least_2_count": sum(1 for virus in group_viruses if virus["z_score"] >= 2.0),
        }

    summary = {
        "seed": seed,
        "simulations_per_virus": simulations,
        "virus_count": len(viruses),
        "positive_z_count": sum(1 for virus in viruses if virus["z_score"] > 0),
        "z_at_least_2_count": sum(1 for virus in viruses if virus["z_score"] >= 2.0),
        "overlap_vs_z_spearman": spearman(
            [virus["overlap_fraction"] for virus in viruses],
            [virus["z_score"] for virus in viruses],
        ),
        "group_summary": group_summary,
        "top_positive_z": sorted(
            ((virus["name"], virus["z_score"]) for virus in viruses),
            key=lambda item: item[1],
            reverse=True,
        )[:5],
        "top_negative_z": sorted(
            ((virus["name"], virus["z_score"]) for virus in viruses),
            key=lambda item: item[1],
        )[:5],
    }

    outdir.mkdir(parents=True, exist_ok=True)
    results_path = outdir / "frameshield_results.json"
    summary_path = outdir / "summary.json"
    results_path.write_text(json.dumps({"viruses": viruses}, indent=2) + "\n")
    summary_path.write_text(json.dumps(summary, indent=2) + "\n")
    return summary


def main() -> None:
    parser = argparse.ArgumentParser(
        description="Benchmark off-frame stop codon enrichment in viral CDS sets against synonymous nulls."
    )
    parser.add_argument("--outdir", default="frameshield_run", help="Directory for benchmark outputs.")
    parser.add_argument(
        "--simulations",
        type=int,
        default=100,
        help="Number of synonymous null genomes to sample per accession.",
    )
    parser.add_argument("--seed", type=int, default=1729, help="Random seed for synonymous sampling.")
    parser.add_argument(
        "--verify",
        action="store_true",
        help="Print a verification marker if the benchmark produces a nontrivial positive signal.",
    )
    args = parser.parse_args()

    summary = run_benchmark(pathlib.Path(args.outdir), args.simulations, args.seed)
    print(json.dumps(summary, indent=2))
    if (
        args.verify
        and summary["group_summary"].get("high-overlap", {}).get("z_at_least_2_count", 0) >= 4
        and summary["group_summary"].get("large-dna", {}).get("positive_z_count", 0) == 0
    ):
        print("frameshield_benchmark_verified")


if __name__ == "__main__":
    main()
PY
chmod +x scripts/frameshield_benchmark.py
```

Expected output: no terminal output; `scripts/frameshield_benchmark.py` exists.

## Step 3: Run the Benchmark

```bash
python3 scripts/frameshield_benchmark.py --outdir frameshield_run --simulations 100 --seed 1729 --verify
```

Expected output:

- a JSON summary printed to stdout
- final line: `frameshield_benchmark_verified`

Expected files:

- `frameshield_run/frameshield_results.json`
- `frameshield_run/summary.json`

## Step 4: Verify the Published Headline Signal

```bash
python3 - <<'PY'
import json
import pathlib

summary = json.loads(pathlib.Path("frameshield_run/summary.json").read_text())
assert summary["virus_count"] == 19, summary
assert summary["positive_z_count"] == 13, summary
assert summary["z_at_least_2_count"] == 4, summary
assert summary["group_summary"]["high-overlap"]["count"] == 8, summary
assert summary["group_summary"]["high-overlap"]["positive_z_count"] == 7, summary
assert summary["group_summary"]["high-overlap"]["z_at_least_2_count"] == 4, summary
assert summary["group_summary"]["large-dna"]["count"] == 3, summary
assert summary["group_summary"]["large-dna"]["positive_z_count"] == 0, summary
assert summary["overlap_vs_z_spearman"] > 0.35, summary
print("frameshield_summary_verified")
PY
```

Expected output:

`frameshield_summary_verified`

## Notes

- The accession panel is fixed inside the script, so the benchmark cohort does not drift.
- The script caches fetched FASTA payloads locally within `frameshield_run/cache` to reduce repeated network load.
- No API keys or non-standard Python packages are required.