---
name: cpg-camouflage
description: >
  Measures how closely RNA viruses mimic their host's CpG dinucleotide suppression
  and tests whether zoonotic spillover viruses show measurable mismatch to their new
  host. Fetches 10 hardcoded NCBI RefSeq viral genomes (human-adapted, bat-associated,
  and spillover), computes observed/expected (O/E) ratios for all 16 dinucleotides,
  calculates camouflage distance to the human CpG O/E baseline (0.23), and asserts
  that human-adapted viruses are better camouflaged than bat-associated viruses.
  Triggers: CpG suppression, ZAP evasion, viral dinucleotide composition, RNA virus
  host adaptation, zoonotic spillover analysis, CpG camouflage benchmark.
allowed-tools: Bash(python3 *), Bash(mkdir *), Bash(cat *), Bash(echo *)
---

## Overview

This skill tests the **CpG camouflage hypothesis**: RNA viruses that have co-evolved
with a specific host evolve to suppress their CpG dinucleotide frequency to match host
levels, evading detection by the zinc-finger antiviral protein (ZAP). Human-adapted
viruses should show lower camouflage distance to the human CpG O/E baseline than
bat-associated viruses that have not yet adapted to human hosts.

**Panel:** 10 hardcoded NCBI RefSeq accessions — 5 human-adapted, 3 bat-associated,
2 outgroups (equine host, bat/human spillover).

**Key metric:** Camouflage distance = |virus_CpG_OE − 0.23|, where 0.23 is the
established human genome CpG O/E (Karlin & Mrazek, *Genome Research* 1997).

**Verification:** `assert human_adapted_mean_distance < bat_associated_mean_distance`
then `print("cpg_camouflage_verified")`

---

## Step 1: Create Workspace

```bash
mkdir -p workspace && cd workspace && mkdir -p data/genomes scripts output
```

Expected output:
```
(no output — directories created silently)
```

---

## Step 2: Fetch Viral Genomes from NCBI

```bash
cd workspace && cat > scripts/fetch_genomes.py <<'PY'
#!/usr/bin/env python3
"""Fetch 10 viral genomes from NCBI EFetch. Rate-limited, with retry logic."""
import urllib.request
import urllib.error
import time
import pathlib
import sys

# Fixed panel — never use "latest" or search-based queries
ACCESSIONS = {
    # Human-adapted viruses (long co-evolutionary history with Homo sapiens)
    "NC_045512.2": "SARS-CoV-2",
    "NC_001802.1": "HIV-1",
    "NC_001474.2": "Dengue-2",
    "NC_004102.1": "HCV",
    "NC_002549.1": "Ebola-Zaire",
    # Bat-associated viruses (primary reservoir: bats; not yet human-adapted)
    "NC_014470.1": "Bat-CoV-HKU9",
    "NC_009019.1": "Bat-CoV-HKU4",
    "NC_025217.1": "Bat-CoV-BM48-31",
    # Outgroups
    "NC_001608.3": "Equine-Arteritis-Virus",  # horse host
    "NC_002640.1": "Nipah-Virus",             # bat reservoir, human spillover
}

NCBI_EFETCH = (
    "https://eutils.ncbi.nlm.nih.gov/entrez/eutils/efetch.fcgi"
    "?db=nuccore&id={acc}&rettype=fasta&retmode=text"
)
MAX_RETRIES = 3
RATE_LIMIT_SLEEP = 0.35  # NCBI allows ~3 req/s without API key


def fetch_with_retry(url, retries=MAX_RETRIES):
    for attempt in range(retries):
        try:
            with urllib.request.urlopen(url, timeout=60) as r:
                return r.read().decode("utf-8")
        except urllib.error.URLError as e:
            if attempt < retries - 1:
                wait = 2 ** attempt  # exponential backoff: 1s, 2s, 4s
                print(f"  Retry {attempt+1}/{retries-1} after {wait}s: {e}", file=sys.stderr)
                time.sleep(wait)
            else:
                raise RuntimeError(f"Failed after {retries} attempts: {e}") from e


out_dir = pathlib.Path("data/genomes")
out_dir.mkdir(parents=True, exist_ok=True)

for acc, name in ACCESSIONS.items():
    url = NCBI_EFETCH.format(acc=acc)
    print(f"Fetching {acc} ({name})...")
    content = fetch_with_retry(url)
    fasta_path = out_dir / f"{acc}.fasta"
    fasta_path.write_text(content)
    size_kb = len(content) / 1024
    print(f"  Saved {acc}.fasta ({size_kb:.1f} KB)")
    time.sleep(RATE_LIMIT_SLEEP)

print(f"\nFetched {len(ACCESSIONS)} genomes to data/genomes/")
PY
python3 scripts/fetch_genomes.py
```

Expected output:
```
Fetching NC_045512.2 (SARS-CoV-2)...
  Saved NC_045512.2.fasta (30.0 KB)
Fetching NC_001802.1 (HIV-1)...
  Saved NC_001802.1.fasta (9.3 KB)
Fetching NC_001474.2 (Dengue-2)...
  Saved NC_001474.2.fasta (10.8 KB)
Fetching NC_004102.1 (HCV)...
  Saved NC_004102.1.fasta (9.8 KB)
Fetching NC_002549.1 (Ebola-Zaire)...
  Saved NC_002549.1.fasta (19.1 KB)
Fetching NC_014470.1 (Bat-CoV-HKU9)...
  Saved NC_014470.1.fasta (31.2 KB)
Fetching NC_009019.1 (Bat-CoV-HKU4)...
  Saved NC_009019.1.fasta (30.4 KB)
Fetching NC_025217.1 (Bat-CoV-BM48-31)...
  Saved NC_025217.1.fasta (30.1 KB)
Fetching NC_001608.3 (Equine-Arteritis-Virus)...
  Saved NC_001608.3.fasta (12.9 KB)
Fetching NC_002640.1 (Nipah-Virus)...
  Saved NC_002640.1.fasta (18.2 KB)

Fetched 10 genomes to data/genomes/
```

---

## Step 3: Write Dinucleotide Analysis Script

```bash
cd workspace && cat > scripts/analyze_cpg.py <<'PY'
#!/usr/bin/env python3
"""Compute dinucleotide O/E ratios and CpG camouflage distance for each viral genome."""
import json
import pathlib
import statistics

# ---------------------------------------------------------------------------
# Host reference
# Human genome CpG O/E = 0.23 (Karlin & Mrazek, Genome Research 1997;
# consistent with Takata et al. J Virol 2017 who confirm ~0.20-0.25 in human
# transcriptome; Fros et al. 2020 PLOS Pathog similarly use 0.25 as human
# baseline for ZAP-evasion analyses).
# ---------------------------------------------------------------------------
HOST_CPG_OE = 0.23  # human genome representative CpG O/E ratio

# Group membership for comparison test
HUMAN_ADAPTED = {"NC_045512.2", "NC_001802.1", "NC_001474.2", "NC_004102.1", "NC_002549.1"}
BAT_ASSOCIATED = {"NC_014470.1", "NC_009019.1", "NC_025217.1"}

ACCESSION_NAMES = {
    "NC_045512.2": "SARS-CoV-2",
    "NC_001802.1": "HIV-1",
    "NC_001474.2": "Dengue-2",
    "NC_004102.1": "HCV",
    "NC_002549.1": "Ebola-Zaire",
    "NC_014470.1": "Bat-CoV-HKU9",
    "NC_009019.1": "Bat-CoV-HKU4",
    "NC_025217.1": "Bat-CoV-BM48-31",
    "NC_001608.3": "Equine-Arteritis-Virus",
    "NC_002640.1": "Nipah-Virus",
}

NUCLEOTIDES = list("ACGT")
DINUCLEOTIDES = [a + b for a in NUCLEOTIDES for b in NUCLEOTIDES]  # 16 pairs


def parse_fasta_sequence(fasta_text):
    """Return the concatenated nucleotide sequence from a FASTA string (uppercase, ACGT only)."""
    lines = fasta_text.strip().splitlines()
    seq_lines = [l for l in lines if not l.startswith(">")]
    seq = "".join(seq_lines).upper()
    # Keep only unambiguous ACGT characters
    seq = "".join(c for c in seq if c in "ACGT")
    return seq


def compute_oe_ratios(seq):
    """Compute O/E for all 16 dinucleotides.

    O/E(XY) = count(XY) / (count(X) * count(Y) / total)
    Returns a dict mapping dinucleotide → O/E ratio.
    """
    n = len(seq)
    if n < 2:
        raise ValueError(f"Sequence too short: {n} nt")

    # Mononucleotide counts
    mono = {nt: seq.count(nt) for nt in NUCLEOTIDES}

    # Dinucleotide counts (overlapping — standard method for genome composition)
    di_counts = {}
    for di in DINUCLEOTIDES:
        count = 0
        for i in range(n - 1):
            if seq[i] == di[0] and seq[i + 1] == di[1]:
                count += 1
        di_counts[di] = count

    total_di = n - 1  # number of dinucleotide positions

    oe = {}
    for di in DINUCLEOTIDES:
        x, y = di[0], di[1]
        expected = (mono[x] * mono[y]) / n  # expected count given total_di ~ n
        if expected == 0:
            oe[di] = 0.0
        else:
            oe[di] = di_counts[di] / expected
    return oe


def main():
    genome_dir = pathlib.Path("data/genomes")
    results = {}

    for acc, name in ACCESSION_NAMES.items():
        fasta_path = genome_dir / f"{acc}.fasta"
        if not fasta_path.exists():
            raise FileNotFoundError(f"Missing genome file: {fasta_path}")

        fasta_text = fasta_path.read_text()
        seq = parse_fasta_sequence(fasta_text)
        oe = compute_oe_ratios(seq)
        cpg_oe = oe["CG"]
        camouflage_distance = abs(cpg_oe - HOST_CPG_OE)

        if acc in HUMAN_ADAPTED:
            group = "human_adapted"
        elif acc in BAT_ASSOCIATED:
            group = "bat_associated"
        else:
            group = "outgroup"

        results[acc] = {
            "name": name,
            "group": group,
            "genome_length": len(seq),
            "cpg_oe": round(cpg_oe, 4),
            "camouflage_distance": round(camouflage_distance, 4),
            "host_cpg_oe": HOST_CPG_OE,
            "all_dinucleotide_oe": {k: round(v, 4) for k, v in oe.items()},
        }
        print(f"{acc:15s} {name:30s} group={group:15s} CpG_OE={cpg_oe:.4f}  dist={camouflage_distance:.4f}")

    # Ranking: lower distance = better camouflaged for humans
    ranking = sorted(results.keys(), key=lambda a: results[a]["camouflage_distance"])

    # Group comparison
    ha_distances = [results[a]["camouflage_distance"] for a in results if results[a]["group"] == "human_adapted"]
    ba_distances = [results[a]["camouflage_distance"] for a in results if results[a]["group"] == "bat_associated"]

    ha_mean = statistics.mean(ha_distances)
    ba_mean = statistics.mean(ba_distances)

    summary = {
        "host_cpg_oe": HOST_CPG_OE,
        "host_cpg_oe_source": "Karlin & Mrazek, Genome Research 1997",
        "human_adapted_mean_distance": round(ha_mean, 4),
        "bat_associated_mean_distance": round(ba_mean, 4),
        "ranking_best_to_worst_camouflage": ranking,
        "viruses": results,
    }

    output_dir = pathlib.Path("output")
    output_dir.mkdir(parents=True, exist_ok=True)
    (output_dir / "results.json").write_text(json.dumps(summary, indent=2))

    print("\n--- Group Summary ---")
    print(f"Human-adapted mean camouflage distance:  {ha_mean:.4f}")
    print(f"Bat-associated mean camouflage distance: {ba_mean:.4f}")
    print(f"\nRanking (best to worst camouflage for humans):")
    for rank, acc in enumerate(ranking, 1):
        v = results[acc]
        print(f"  {rank:2d}. {acc:15s} {v['name']:30s} dist={v['camouflage_distance']:.4f}")

    print("\nResults written to output/results.json")


if __name__ == "__main__":
    main()
PY
python3 scripts/analyze_cpg.py
```

Expected output:
```
NC_045512.2     SARS-CoV-2                     group=human_adapted   CpG_OE=0.4077  dist=0.1777
NC_001802.1     HIV-1                          group=human_adapted   CpG_OE=0.2061  dist=0.0239
NC_001474.2     Dengue-2                       group=human_adapted   CpG_OE=0.4114  dist=0.1814
NC_004102.1     HCV                            group=human_adapted   CpG_OE=0.7307  dist=0.5007
NC_002549.1     Ebola-Zaire                    group=human_adapted   CpG_OE=0.6049  dist=0.3749
NC_014470.1     Bat-CoV-HKU9                   group=bat_associated  CpG_OE=0.5110  dist=0.2810
NC_009019.1     Bat-CoV-HKU4                   group=bat_associated  CpG_OE=0.5115  dist=0.2815
NC_025217.1     Bat-CoV-BM48-31                group=bat_associated  CpG_OE=0.5891  dist=0.3591
NC_001608.3     Equine-Arteritis-Virus         group=outgroup        CpG_OE=0.5300  dist=0.3000
NC_002640.1     Nipah-Virus                    group=outgroup        CpG_OE=0.3842  dist=0.1542

--- Group Summary ---
Human-adapted mean camouflage distance:  0.2517
Bat-associated mean camouflage distance: 0.3072

Ranking (best to worst camouflage for humans):
   1. NC_001802.1     HIV-1                          dist=0.0239
   2. NC_002640.1     Nipah-Virus                    dist=0.1542
   3. NC_045512.2     SARS-CoV-2                     dist=0.1777
   4. NC_001474.2     Dengue-2                       dist=0.1814
   5. NC_014470.1     Bat-CoV-HKU9                   dist=0.2810
   6. NC_009019.1     Bat-CoV-HKU4                   dist=0.2815
   7. NC_001608.3     Equine-Arteritis-Virus         dist=0.3000
   8. NC_025217.1     Bat-CoV-BM48-31                dist=0.3591
   9. NC_002549.1     Ebola-Zaire                    dist=0.3749
  10. NC_004102.1     HCV                            dist=0.5007

Results written to output/results.json
```

---

## Step 4: Run Smoke Tests

```bash
cd workspace && python3 - <<'PY'
#!/usr/bin/env python3
"""Smoke tests: validate genome files, O/E plausibility, and output structure."""
import json
import pathlib
import sys

EXPECTED_ACCESSIONS = [
    "NC_045512.2", "NC_001802.1", "NC_001474.2", "NC_004102.1", "NC_002549.1",
    "NC_014470.1", "NC_009019.1", "NC_025217.1",
    "NC_001608.3", "NC_002640.1",
]
DINUCLEOTIDES = [a + b for a in "ACGT" for b in "ACGT"]

errors = []

# ---- Test 1: All 10 genome files exist and have non-zero size ----
genome_dir = pathlib.Path("data/genomes")
for acc in EXPECTED_ACCESSIONS:
    fpath = genome_dir / f"{acc}.fasta"
    if not fpath.exists():
        errors.append(f"MISSING genome file: {fpath}")
    elif fpath.stat().st_size == 0:
        errors.append(f"EMPTY genome file: {fpath}")
    else:
        print(f"  [OK] {acc}.fasta ({fpath.stat().st_size:,} bytes)")

print(f"\nTest 1 (genome files): {'PASS' if not errors else 'FAIL'}")

# ---- Test 2: Load results JSON and verify structure ----
results_path = pathlib.Path("output/results.json")
if not results_path.exists():
    errors.append("MISSING output/results.json")
    print("Test 2 FAIL — output file missing, cannot continue")
    sys.exit(1)

data = json.loads(results_path.read_text())

required_top_keys = [
    "host_cpg_oe", "host_cpg_oe_source",
    "human_adapted_mean_distance", "bat_associated_mean_distance",
    "ranking_best_to_worst_camouflage", "viruses"
]
for key in required_top_keys:
    if key not in data:
        errors.append(f"MISSING top-level key in results.json: {key}")

print(f"\nTest 2 (output JSON keys): {'PASS' if not errors else 'FAIL'}")

# ---- Test 3: Verify all 10 viruses are in results ----
viruses = data["viruses"]
for acc in EXPECTED_ACCESSIONS:
    if acc not in viruses:
        errors.append(f"MISSING virus in results: {acc}")
    else:
        print(f"  [OK] {acc} present in results")

print(f"\nTest 3 (all 10 viruses in results): {'PASS' if not errors else 'FAIL'}")

# ---- Test 4: All O/E ratios are biologically plausible (0 < OE < 2) ----
for acc, v in viruses.items():
    for di, oe_val in v["all_dinucleotide_oe"].items():
        if not (0.0 <= oe_val <= 2.0):
            errors.append(
                f"O/E out of range [0,2] for {acc} dinucleotide {di}: {oe_val}"
            )

print(f"\nTest 4 (all O/E in [0, 2]): {'PASS' if not errors else 'FAIL'}")

# ---- Test 5: CpG O/E < 1.0 for all viruses (universal CpG suppression) ----
for acc, v in viruses.items():
    cpg_oe = v["cpg_oe"]
    if cpg_oe >= 1.0:
        errors.append(
            f"CpG O/E >= 1.0 for {acc} ({v['name']}): {cpg_oe} — unexpected (CpG suppression not present)"
        )
    else:
        print(f"  [OK] {acc} ({v['name']}) CpG_OE={cpg_oe:.4f} < 1.0")

print(f"\nTest 5 (CpG O/E < 1.0 for all): {'PASS' if not errors else 'FAIL'}")

# ---- Test 6: All camouflage distances are non-negative ----
for acc, v in viruses.items():
    dist = v["camouflage_distance"]
    if dist < 0:
        errors.append(f"Negative camouflage_distance for {acc}: {dist}")

print(f"\nTest 6 (camouflage distances non-negative): {'PASS' if not errors else 'FAIL'}")

# ---- Test 7: Ranking has all 10 entries ----
ranking = data["ranking_best_to_worst_camouflage"]
if len(ranking) != 10:
    errors.append(f"Ranking length {len(ranking)} != 10")

print(f"\nTest 7 (ranking has 10 entries): {'PASS' if not errors else 'FAIL'}")

# ---- Final report ----
if errors:
    print(f"\n{'='*60}")
    print(f"SMOKE TESTS FAILED — {len(errors)} error(s):")
    for e in errors:
        print(f"  ERROR: {e}")
    sys.exit(1)
else:
    print(f"\n{'='*60}")
    print("All 7 smoke tests passed.")
    print("smoke_tests_passed")
PY
```

Expected output:
```
  [OK] NC_045512.2.fasta (...) bytes)
  ...
Test 1 (genome files): PASS
Test 2 (output JSON keys): PASS
  [OK] NC_045512.2 present in results
  ...
Test 3 (all 10 viruses in results): PASS
Test 4 (all O/E in [0, 2]): PASS
  [OK] NC_045512.2 (SARS-CoV-2) CpG_OE=0.xxxx < 1.0
  ...
Test 5 (CpG O/E < 1.0 for all): PASS
Test 6 (camouflage distances non-negative): PASS
Test 7 (ranking has 10 entries): PASS
============================================================
All 7 smoke tests passed.
smoke_tests_passed
```

---

## Step 5: Verify Results

```bash
cd workspace && python3 - <<'PY'
#!/usr/bin/env python3
"""Final verification: assert core scientific hypothesis and print marker."""
import json
import pathlib

results_path = pathlib.Path("output/results.json")
assert results_path.exists(), "output/results.json not found — run analysis first"

data = json.loads(results_path.read_text())

ha_mean = data["human_adapted_mean_distance"]
ba_mean = data["bat_associated_mean_distance"]

print(f"Human-adapted mean camouflage distance:  {ha_mean:.4f}")
print(f"Bat-associated mean camouflage distance: {ba_mean:.4f}")

# Core hypothesis assertion
assert ha_mean < ba_mean, (
    f"Hypothesis FAILED: human-adapted mean ({ha_mean:.4f}) >= "
    f"bat-associated mean ({ba_mean:.4f}). "
    "Expected human-adapted viruses to be better camouflaged for the human host."
)

# Sanity: group means are plausible distances
assert 0.0 <= ha_mean <= 1.0, f"Human-adapted mean out of plausible range: {ha_mean}"
assert 0.0 <= ba_mean <= 1.0, f"Bat-associated mean out of plausible range: {ba_mean}"

# Sanity: all 10 viruses present
assert len(data["viruses"]) == 10, f"Expected 10 viruses, got {len(data['viruses'])}"

# Sanity: ranking has 10 entries in correct order
ranking = data["ranking_best_to_worst_camouflage"]
assert len(ranking) == 10, f"Ranking should have 10 entries, got {len(ranking)}"
distances = [data["viruses"][acc]["camouflage_distance"] for acc in ranking]
assert distances == sorted(distances), "Ranking is not sorted by camouflage distance"

print("\nAll assertions passed.")
print("cpg_camouflage_verified")
PY
```

Expected output:
```
Human-adapted mean camouflage distance:  0.2517
Bat-associated mean camouflage distance: 0.3072

All assertions passed.
cpg_camouflage_verified
```

---

## Notes / Limitations

- **CpG is one of 16 dinucleotides.** Full dinucleotide O/E profiles are computed and saved for all 16 pairs. CpG is the focus because it is the ZAP-detected signal; other suppressed dinucleotides (TpA, CpA) are available in `output/results.json` for follow-up analysis.
- **Host CpG O/E is approximated from literature.** The value 0.23 is derived from Karlin & Mrazek (*Genome Research* 1997) and corroborated by Takata et al. (*J Virol* 2017) and Fros et al. (*PLOS Pathogens* 2020). It is not re-derived here from a human reference genome assembly (which would require a >3 GB download outside the scope of this skill).
- **Small panel (n=10).** The 5 vs. 3 comparison (human-adapted vs. bat-associated) has limited statistical power. The direction of the result is the primary finding.
- **Phylogenetic non-independence.** The three bat coronaviruses share common ancestry; they are not statistically independent observations.
- **Single-sequence treatment.** CpG composition is computed across the full genome sequence. Some viruses encode CpG-rich accessory genes embedded in otherwise suppressed genomes; per-gene analysis would add resolution.
- **Ebola and Nipah** are biosafety-critical pathogens analyzed here solely at the sequence-composition level. No infectious material is implied or produced.