---
name: palindrome-deserts
description: >
  Measures palindromic sequence avoidance in bacterial genomes and tests whether
  the most-avoided 4-mer palindromes correspond to restriction enzyme recognition
  sites, replicating Gelfand & Koonin (1997). Fetches 5 bacterial genomes (E. coli,
  B. subtilis, P. aeruginosa, M. tuberculosis, S. coelicolor) and 2 organellar
  controls (human mitochondria, Arabidopsis chloroplast) from NCBI RefSeq using
  hardcoded accessions. Computes observed/expected ratios for all 16 palindromic
  4-mers and asserts that CTAG is the most-avoided palindrome in E. coli.
  Triggers: palindrome avoidance, restriction enzyme recognition sites, bacterial
  genomes, O/E ratio, Gelfand Koonin 1997, restriction-modification.
allowed-tools: Bash(python3 *), Bash(mkdir *), Bash(cat *), Bash(echo *)
---

# Palindrome Deserts

Replicates and extends Gelfand & Koonin (1997): bacterial genomes systematically
under-represent 4-mer palindromic sequences that overlap with restriction enzyme
recognition sites, while organellar genomes lacking restriction-modification (R-M)
systems show no such avoidance.

**Data source:** NCBI RefSeq — 5 bacterial genomes + 2 organellar controls  
**Analysis:** Observed/Expected ratio for all 16 palindromic 4-mers per genome  
**Output:** `output/results.json` with per-genome O/E tables, rankings, and comparison  
**Verification:** CTAG is most avoided in E. coli; bacterial mean O/E < organellar mean O/E

---

## Step 1: Setup Workspace

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

Expected output:
```
(no output — directories data/, scripts/, output/ created)
```

---

## Step 2: Fetch Bacterial and Organellar Genomes from NCBI

```bash
cd workspace
python3 - <<'PY'
"""Fetch 7 complete genomes (5 bacterial + 2 organellar) from NCBI EFetch."""
import urllib.request
import urllib.error
import time
import pathlib

GENOMES = [
    ("NC_000913.3", "Escherichia coli K-12 MG1655",           "bacterial"),
    ("NC_000964.3", "Bacillus subtilis 168",                   "bacterial"),
    ("NC_002516.2", "Pseudomonas aeruginosa PAO1",             "bacterial"),
    ("NC_000962.3", "Mycobacterium tuberculosis H37Rv",        "bacterial"),
    ("NC_003888.3", "Streptomyces coelicolor A3(2)",           "bacterial"),
    ("NC_012920.1", "Homo sapiens mitochondrion",              "organellar"),
    ("NC_000932.1", "Arabidopsis thaliana chloroplast",        "organellar"),
]

EFETCH_URL = (
    "https://eutils.ncbi.nlm.nih.gov/entrez/eutils/efetch.fcgi"
    "?db=nuccore&id={acc}&rettype=fasta&retmode=text"
)

def fetch_with_retry(url, retries=3):
    for attempt in range(retries):
        try:
            req = urllib.request.Request(url, headers={"User-Agent": "PalindromeDeserts/1.0"})
            with urllib.request.urlopen(req, timeout=120) as r:
                return r.read().decode("utf-8")
        except (urllib.error.URLError, urllib.error.HTTPError) as e:
            if attempt < retries - 1:
                wait = 2 ** attempt
                print(f"  Attempt {attempt+1} failed ({e}), retrying in {wait}s...")
                time.sleep(wait)
            else:
                raise RuntimeError(f"Failed to fetch {url} after {retries} attempts: {e}") from e

pathlib.Path("data").mkdir(exist_ok=True)
for acc, name, category in GENOMES:
    out_path = pathlib.Path(f"data/{acc}.fasta")
    if out_path.exists() and out_path.stat().st_size > 1000:
        print(f"  CACHED {acc} ({name})")
        continue
    url = EFETCH_URL.format(acc=acc)
    print(f"  Fetching {acc} ({name})...")
    content = fetch_with_retry(url)
    out_path.write_text(content, encoding="utf-8")
    size_kb = out_path.stat().st_size // 1024
    print(f"  Saved {acc}: {size_kb} KB")
    time.sleep(0.4)  # NCBI rate limit: max 3 requests/second without API key

print(f"\nAll {len(GENOMES)} genomes fetched.")
PY
```

Expected output:
```
  Fetching NC_000913.3 (Escherichia coli K-12 MG1655)...
  Saved NC_000913.3: 4814 KB
  Fetching NC_000964.3 (Bacillus subtilis 168)...
  Saved NC_000964.3: 4161 KB
  Fetching NC_002516.2 (Pseudomonas aeruginosa PAO1)...
  Saved NC_002516.2: 6264 KB
  Fetching NC_000962.3 (Mycobacterium tuberculosis H37Rv)...
  Saved NC_000962.3: 4411 KB
  Fetching NC_003888.3 (Streptomyces coelicolor A3(2))...
  Saved NC_003888.3: 8977 KB
  Fetching NC_012920.1 (Homo sapiens mitochondrion)...
  Saved NC_012920.1: 16 KB
  Fetching NC_000932.1 (Arabidopsis thaliana chloroplast)...
  Saved NC_000932.1: 153 KB

All 7 genomes fetched.
```

---

## Step 3: Write Palindrome Analysis Script

```bash
cd workspace
cat > scripts/palindrome_analysis.py <<'PY'
#!/usr/bin/env python3
"""
Compute observed/expected ratios for all 16 palindromic 4-mers across
7 genomes (5 bacterial + 2 organellar). Replicates Gelfand & Koonin (1997).
"""
import json
import pathlib
import collections


# ---------------------------------------------------------------------------
# Constants
# ---------------------------------------------------------------------------

GENOMES = [
    ("NC_000913.3", "Escherichia coli K-12 MG1655",           "bacterial"),
    ("NC_000964.3", "Bacillus subtilis 168",                   "bacterial"),
    ("NC_002516.2", "Pseudomonas aeruginosa PAO1",             "bacterial"),
    ("NC_000962.3", "Mycobacterium tuberculosis H37Rv",        "bacterial"),
    ("NC_003888.3", "Streptomyces coelicolor A3(2)",           "bacterial"),
    ("NC_012920.1", "Homo sapiens mitochondrion",              "organellar"),
    ("NC_000932.1", "Arabidopsis thaliana chloroplast",        "organellar"),
]

# All 16 palindromic 4-mers (sequence == reverse complement of itself)
PALINDROMES_4MER = [
    "AATT", "ATAT", "ACGT", "AGCT",
    "TATA", "TTAA", "TCGA", "TGCA",
    "CATG", "CTAG", "CCGG", "CGCG",
    "GATC", "GTAC", "GCGC", "GGCC",
]

# Known restriction enzyme recognition sites mapping (4-mer level)
# Source: REBASE / standard molecular biology references
RESTRICTION_ENZYMES = {
    "AATT": "EcoRI core (GAATTC contains AATT); also MfeI",
    "ATAT": "EcoRV-related (GATATC contains ATAT)",
    "ACGT": "AciI/ScrFI recognition",
    "AGCT": "HindIII core (AAGCTT contains AGCT); AluI (AGCT)",
    "TATA": "Complement of TATA; near TATA-box motifs",
    "TTAA": "MseI recognition (TTAA); AseI (ATTAAT contains TTAA)",
    "TCGA": "TaqI recognition (TCGA); SalI core",
    "TGCA": "PstI core (CTGCAG contains TGCA)",
    "CATG": "NlaIII recognition (CATG); SphI core (GCATGC contains CATG)",
    "CTAG": "XbaI core (TCTAGA contains CTAG); PvuI-related; SpeI core (ACTAGT contains CTAG)",
    "CCGG": "HpaII/MspI recognition (CCGG); EagI-related",
    "CGCG": "ThaI recognition (CGCG); BssHII core",
    "GATC": "MboI/Sau3AI recognition (GATC); BamHI core (GGATCC contains GATC); Dam methylase target",
    "GTAC": "KpnI core (GGTACC contains GTAC); Acc65I",
    "GCGC": "BstUI recognition (CGCG complement); related to CpG methylation sites",
    "GGCC": "HaeIII recognition (GGCC); NotI core (GCGGCCGC contains GGCC)",
}


# ---------------------------------------------------------------------------
# Helper functions
# ---------------------------------------------------------------------------

def reverse_complement(seq):
    """Return the reverse complement of a DNA sequence."""
    comp = str.maketrans("ACGTacgt", "TGCAtgca")
    return seq.translate(comp)[::-1]


def verify_palindromes():
    """Assert that all listed 4-mers are truly self-complementary."""
    for pal in PALINDROMES_4MER:
        rc = reverse_complement(pal)
        assert pal == rc, f"{pal} is not palindromic: RC={rc}"
    return True


def parse_fasta(fasta_text):
    """Parse a FASTA file and return the concatenated sequence (uppercase, ACGT only)."""
    lines = fasta_text.strip().splitlines()
    seq_parts = []
    for line in lines:
        if line.startswith(">"):
            continue
        cleaned = line.strip().upper()
        # Keep only canonical bases
        for ch in cleaned:
            if ch in "ACGT":
                seq_parts.append(ch)
    return "".join(seq_parts)


def nucleotide_frequencies(seq):
    """Compute mononucleotide frequencies from a sequence."""
    counts = collections.Counter(seq)
    total = sum(counts[b] for b in "ACGT")
    if total == 0:
        return {"A": 0.25, "C": 0.25, "G": 0.25, "T": 0.25}
    return {b: counts.get(b, 0) / total for b in "ACGT"}


def count_4mers(seq):
    """Count all 4-mer occurrences in seq using a sliding window."""
    counts = collections.Counter()
    n = len(seq)
    for i in range(n - 3):
        counts[seq[i:i+4]] += 1
    return counts


def expected_count(kmer, freqs, genome_length):
    """
    Compute expected count of a k-mer under the mononucleotide independence model.
    expected = freq(pos1) * freq(pos2) * freq(pos3) * freq(pos4) * (L - 3)
    """
    prob = 1.0
    for base in kmer:
        prob *= freqs.get(base, 0.0)
    return prob * (genome_length - 3)


def analyze_genome(acc, name, category):
    """Load and analyze one genome. Returns dict with O/E table and ranking."""
    fasta_path = pathlib.Path(f"data/{acc}.fasta")
    assert fasta_path.exists(), f"FASTA not found: {fasta_path}"

    fasta_text = fasta_path.read_text(encoding="utf-8")
    seq = parse_fasta(fasta_text)
    L = len(seq)

    freqs = nucleotide_frequencies(seq)
    gc_content = freqs["G"] + freqs["C"]

    kmer_counts = count_4mers(seq)

    palindrome_results = {}
    for pal in PALINDROMES_4MER:
        obs = kmer_counts.get(pal, 0)
        exp = expected_count(pal, freqs, L)
        oe = obs / exp if exp > 0 else float("nan")
        palindrome_results[pal] = {
            "observed": obs,
            "expected": round(exp, 2),
            "oe_ratio": round(oe, 4),
            "restriction_enzyme_note": RESTRICTION_ENZYMES.get(pal, ""),
        }

    # Rank by O/E ratio (ascending = most avoided first)
    ranking = sorted(
        PALINDROMES_4MER,
        key=lambda p: palindrome_results[p]["oe_ratio"]
    )

    return {
        "accession": acc,
        "name": name,
        "category": category,
        "genome_length_bp": L,
        "gc_content": round(gc_content, 4),
        "nucleotide_freqs": {b: round(freqs[b], 6) for b in "ACGT"},
        "palindrome_oe": palindrome_results,
        "ranking_most_to_least_avoided": ranking,
        "most_avoided": ranking[0],
        "top3_avoided": ranking[:3],
        "top3_mean_oe": round(
            sum(palindrome_results[p]["oe_ratio"] for p in ranking[:3]) / 3, 4
        ),
    }


# ---------------------------------------------------------------------------
# Main
# ---------------------------------------------------------------------------

def main():
    verify_palindromes()
    print(f"Verified: all {len(PALINDROMES_4MER)} listed 4-mers are palindromic")

    results = {}
    for acc, name, category in GENOMES:
        print(f"\nAnalyzing {acc} ({name})...")
        genome_result = analyze_genome(acc, name, category)
        results[acc] = genome_result
        L = genome_result["genome_length_bp"]
        gc = genome_result["gc_content"] * 100
        most_av = genome_result["most_avoided"]
        most_oe = genome_result["palindrome_oe"][most_av]["oe_ratio"]
        top3_oe = genome_result["top3_mean_oe"]
        print(f"  Length: {L:,} bp  GC: {gc:.1f}%")
        print(f"  Most avoided: {most_av} (O/E={most_oe:.4f})")
        print(f"  Top-3 mean O/E: {top3_oe:.4f}")
        print(f"  Full ranking (most→least avoided):")
        for i, pal in enumerate(genome_result["ranking_most_to_least_avoided"], 1):
            oe = genome_result["palindrome_oe"][pal]["oe_ratio"]
            re_note = genome_result["palindrome_oe"][pal]["restriction_enzyme_note"]
            flag = " ***" if oe < 0.5 else ""
            print(f"    {i:2d}. {pal}  O/E={oe:.4f}  {re_note[:50]}{flag}")

    # Comparison: bacterial vs organellar mean O/E for top-3 avoided palindromes
    bacterial_top3_oe = [
        results[acc]["top3_mean_oe"]
        for acc, _, cat in GENOMES if cat == "bacterial"
    ]
    organellar_top3_oe = [
        results[acc]["top3_mean_oe"]
        for acc, _, cat in GENOMES if cat == "organellar"
    ]
    mean_bacterial = sum(bacterial_top3_oe) / len(bacterial_top3_oe)
    mean_organellar = sum(organellar_top3_oe) / len(organellar_top3_oe)

    print(f"\n=== Bacterial vs Organellar Avoidance Comparison ===")
    print(f"Bacterial genomes top-3 mean O/E:   {[round(x, 4) for x in bacterial_top3_oe]}")
    print(f"Organellar genomes top-3 mean O/E:  {[round(x, 4) for x in organellar_top3_oe]}")
    print(f"Mean bacterial top-3 O/E:   {mean_bacterial:.4f}")
    print(f"Mean organellar top-3 O/E:  {mean_organellar:.4f}")
    print(f"Bacterial < Organellar: {mean_bacterial < mean_organellar}")

    # E. coli specific summary
    ecoli = results["NC_000913.3"]
    print(f"\n=== E. coli K-12 Most Avoided Palindromes ===")
    for i, pal in enumerate(ecoli["ranking_most_to_least_avoided"][:5], 1):
        oe = ecoli["palindrome_oe"][pal]["oe_ratio"]
        re_note = ecoli["palindrome_oe"][pal]["restriction_enzyme_note"]
        print(f"  {i}. {pal}  O/E={oe:.4f}  ({re_note[:60]})")

    # Write output JSON
    output = {
        "genomes": results,
        "comparison": {
            "bacterial_accessions": [acc for acc, _, cat in GENOMES if cat == "bacterial"],
            "organellar_accessions": [acc for acc, _, cat in GENOMES if cat == "organellar"],
            "bacterial_top3_mean_oe_per_genome": {
                acc: results[acc]["top3_mean_oe"]
                for acc, _, cat in GENOMES if cat == "bacterial"
            },
            "organellar_top3_mean_oe_per_genome": {
                acc: results[acc]["top3_mean_oe"]
                for acc, _, cat in GENOMES if cat == "organellar"
            },
            "mean_bacterial_top3_oe": round(mean_bacterial, 4),
            "mean_organellar_top3_oe": round(mean_organellar, 4),
            "bacterial_avoidance_stronger": mean_bacterial < mean_organellar,
        },
        "palindrome_list": PALINDROMES_4MER,
        "restriction_enzyme_notes": RESTRICTION_ENZYMES,
        "limitations": [
            "Only 4-mer palindromes analyzed; 6-mer palindromes (the standard in Gelfand & Koonin 1997) would be more informative but require larger genomes to achieve statistical power.",
            "Expected counts use a mononucleotide independence model; dinucleotide or higher-order models would be more accurate.",
            "Small panel of 5 bacterial genomes; a larger phylogenetically diverse panel would strengthen the generalizability.",
            "O/E ratios are asymmetric (bounded below at 0, unbounded above); log transformation would be more appropriate for statistical comparison.",
            "The organellar genomes (mitochondria, chloroplast) have distinct replication mechanisms and evolutionary pressures beyond just the absence of R-M systems.",
            "Observed counts on + strand only; some analyses use both strands (results differ by factor close to 2 for non-palindromes, but identical for palindromes due to symmetry).",
        ],
    }

    out_path = pathlib.Path("output/results.json")
    out_path.write_text(json.dumps(output, indent=2), encoding="utf-8")
    print(f"\nResults written to {out_path}")
    print(f"Output file size: {out_path.stat().st_size:,} bytes")


if __name__ == "__main__":
    main()
PY
echo "Script written to scripts/palindrome_analysis.py"
```

Expected output:
```
Script written to scripts/palindrome_analysis.py
```

---

## Step 4: Run Palindrome Analysis

```bash
cd workspace
python3 scripts/palindrome_analysis.py
```

Expected output:
```
Verified: all 16 listed 4-mers are palindromic

Analyzing NC_000913.3 (Escherichia coli K-12 MG1655)...
  Length: 4,641,652 bp  GC: 50.8%
  Most avoided: CTAG (O/E=0.0488)
  Top-3 mean O/E: 0.4115
  Full ranking (most->least avoided):
     1. CTAG  O/E=0.0488  XbaI core (TCTAGA contains CTAG); PvuI-related; Sp ***
     2. TATA  O/E=0.5338  Complement of TATA; near TATA-box motifs
     ...

Analyzing NC_000964.3 (Bacillus subtilis 168)...
  Most avoided: CTAG (O/E=0.1953)
  ...

=== Bacterial vs Organellar Avoidance Comparison ===
Bacterial genomes top-3 mean O/E:   [0.4115, 0.4176, 0.2731, 0.3146, 0.1039]
Organellar genomes top-3 mean O/E:  [0.355, 0.5855]
Mean bacterial top-3 O/E:   0.3041
Mean organellar top-3 O/E:  0.4703
Bacterial < Organellar: True

=== E. coli K-12 Most Avoided Palindromes ===
  1. CTAG  O/E=0.0488  (XbaI core (TCTAGA contains CTAG); PvuI-related; ...)
  2. TATA  O/E=0.5338  (Complement of TATA; near TATA-box motifs)
  ...

Results written to output/results.json
Output file size: 31,747 bytes
```

---

## Step 5: Smoke Tests

```bash
cd workspace
python3 - <<'PY'
"""Smoke tests for palindrome-deserts analysis."""
import json
import pathlib

print("Running smoke tests...")

results_path = pathlib.Path("output/results.json")
assert results_path.exists(), "output/results.json not found — run Step 4 first"
data = json.loads(results_path.read_text())

# Test 1: All 7 genomes fetched
genomes = data["genomes"]
assert len(genomes) == 7, f"Expected 7 genomes, got {len(genomes)}"
print(f"  PASS: 7 genomes present in output")

# Test 2: Exactly 16 palindromes analyzed per genome
for acc, gdata in genomes.items():
    n_pals = len(gdata["palindrome_oe"])
    assert n_pals == 16, f"Expected 16 palindromes for {acc}, got {n_pals}"
print(f"  PASS: Exactly 16 palindromes analyzed per genome")

# Test 3: All O/E ratios > 0 and < 5 (biologically plausible)
for acc, gdata in genomes.items():
    for pal, info in gdata["palindrome_oe"].items():
        oe = info["oe_ratio"]
        assert 0 < oe < 5, (
            f"O/E ratio out of plausible range for {pal} in {acc}: {oe}. "
            "Expected 0 < O/E < 5."
        )
print(f"  PASS: All O/E ratios in biologically plausible range (0, 5)")

# Test 4: Bacterial genomes show at least some O/E < 0.5 (strong avoidance)
bacterial_accs = data["comparison"]["bacterial_accessions"]
for acc in bacterial_accs:
    oe_vals = [info["oe_ratio"] for info in genomes[acc]["palindrome_oe"].values()]
    min_oe = min(oe_vals)
    assert min_oe < 0.5, (
        f"Bacterial genome {acc} has minimum O/E={min_oe:.4f}. "
        "Expected at least one palindrome with O/E < 0.5 (strong avoidance)."
    )
print(f"  PASS: All bacterial genomes show at least one O/E < 0.5 (strong avoidance)")

# Test 5: Organellar genomes do NOT show CTAG as most avoided palindrome
# (CTAG avoidance is the hallmark of bacterial R-M systems; organelles lack R-M)
# Also: organellar CTAG O/E should be > 0.5 (not strongly avoided)
organellar_accs = data["comparison"]["organellar_accessions"]
for acc in organellar_accs:
    most_avoided = genomes[acc]["most_avoided"]
    ctag_oe = genomes[acc]["palindrome_oe"]["CTAG"]["oe_ratio"]
    assert most_avoided != "CTAG", (
        f"Organellar genome {acc} shows CTAG as most avoided (O/E={ctag_oe:.4f}). "
        "Organelles lack R-M systems and should not show bacterial-style CTAG avoidance."
    )
    assert ctag_oe > 0.5, (
        f"Organellar genome {acc} has CTAG O/E={ctag_oe:.4f} < 0.5. "
        "Organelles lack restriction-modification systems; CTAG should not be strongly avoided."
    )
print(f"  PASS: No organellar genome shows CTAG as most avoided palindrome (bacterial-specific pattern absent)")

# Test 6: Palindrome list verified (all 16 are self-complementary)
COMP = str.maketrans("ACGT", "TGCA")
palindromes = data["palindrome_list"]
assert len(palindromes) == 16, f"Expected 16 palindromes, got {len(palindromes)}"
for pal in palindromes:
    rc = pal.translate(COMP)[::-1]
    assert pal == rc, f"{pal} is not palindromic (RC={rc})"
print(f"  PASS: All 16 palindromes in list are verified self-complementary")

print("\nAll smoke tests passed.")
print("smoke_tests_passed")
PY
```

Expected output:
```
Running smoke tests...
  PASS: 7 genomes present in output
  PASS: Exactly 16 palindromes analyzed per genome
  PASS: All O/E ratios in biologically plausible range (0, 5)
  PASS: All bacterial genomes show at least one O/E < 0.5 (strong avoidance)
  PASS: No organellar genome shows CTAG as most avoided palindrome (bacterial-specific pattern absent)
  PASS: All 16 palindromes in list are verified self-complementary

All smoke tests passed.
smoke_tests_passed
```

---

## Step 6: Verify Results

```bash
cd workspace
python3 - <<'PY'
"""Final verification for palindrome-deserts skill."""
import json
import pathlib

results_path = pathlib.Path("output/results.json")
assert results_path.exists(), "output/results.json not found — run Steps 4 and 5 first"
assert results_path.stat().st_size > 0, "output/results.json is empty"

data = json.loads(results_path.read_text())
genomes = data["genomes"]
comparison = data["comparison"]

# Assertion 1: Most-avoided palindrome in E. coli is CTAG
# (Gelfand & Koonin 1997, PMID 9171096)
ecoli = genomes["NC_000913.3"]
ecoli_most_avoided = ecoli["most_avoided"]
ecoli_ctag_oe = ecoli["palindrome_oe"]["CTAG"]["oe_ratio"]
assert ecoli_most_avoided == "CTAG", (
    f"Expected CTAG to be most avoided in E. coli K-12, but got {ecoli_most_avoided} "
    f"(CTAG O/E={ecoli_ctag_oe:.4f}). "
    "This replicates the central finding of Gelfand & Koonin (1997) PMID:9171096."
)

# Assertion 2: Mean bacterial O/E for top-3 avoided < mean organellar O/E
mean_bact = comparison["mean_bacterial_top3_oe"]
mean_org  = comparison["mean_organellar_top3_oe"]
assert mean_bact < mean_org, (
    f"Expected mean bacterial top-3 O/E ({mean_bact:.4f}) < "
    f"mean organellar top-3 O/E ({mean_org:.4f}). "
    "Bacteria with R-M systems should show stronger palindrome avoidance than organelles."
)

# Print summary
print(f"E. coli K-12 most avoided palindrome: {ecoli_most_avoided} (O/E={ecoli_ctag_oe:.4f})")
print(f"  -> Matches Gelfand & Koonin (1997) CTAG avoidance result: CONFIRMED")
print(f"Mean bacterial top-3 O/E:   {mean_bact:.4f}")
print(f"Mean organellar top-3 O/E:  {mean_org:.4f}")
print(f"  -> Bacterial avoidance stronger than organellar: CONFIRMED")

# Print per-genome summary
print("\nPer-genome top-3 most avoided palindromes:")
for acc, _, cat in [
    ("NC_000913.3", "E. coli K-12",        "bacterial"),
    ("NC_000964.3", "B. subtilis 168",      "bacterial"),
    ("NC_002516.2", "P. aeruginosa PAO1",   "bacterial"),
    ("NC_000962.3", "M. tuberculosis H37Rv","bacterial"),
    ("NC_003888.3", "S. coelicolor A3(2)",  "bacterial"),
    ("NC_012920.1", "Human mitochondria",   "organellar"),
    ("NC_000932.1", "A. thaliana chloroplast","organellar"),
]:
    g = genomes[acc]
    top3 = g["top3_avoided"]
    oes = [g["palindrome_oe"][p]["oe_ratio"] for p in top3]
    gc_pct = g["gc_content"] * 100
    print(f"  [{cat[:4]}] {acc} GC={gc_pct:.0f}%: {top3[0]}({oes[0]:.3f}) {top3[1]}({oes[1]:.3f}) {top3[2]}({oes[2]:.3f})")

print("palindrome_deserts_verified")
PY
```

Expected output:
```
E. coli K-12 most avoided palindrome: CTAG (O/E=0.0488)
  -> Matches Gelfand & Koonin (1997) CTAG avoidance result: CONFIRMED
Mean bacterial top-3 O/E:   0.3041
Mean organellar top-3 O/E:  0.4703
  -> Bacterial avoidance stronger than organellar: CONFIRMED

Per-genome top-3 most avoided palindromes:
  [bact] NC_000913.3 GC=51%: CTAG(0.049) TATA(0.534) GGCC(0.652)
  [bact] NC_000964.3 GC=44%: CTAG(0.195) GTAC(0.490) ACGT(0.568)
  [bact] NC_002516.2 GC=67%: CTAG(0.099) TTAA(0.239) TATA(0.482)
  [bact] NC_000962.3 GC=66%: CTAG(0.281) TTAA(0.330) TATA(0.333)
  [bact] NC_003888.3 GC=72%: CTAG(0.064) TATA(0.116) TTAA(0.131)
  [orga] NC_012920.1 GC=44%: CGCG(0.216) ACGT(0.405) GATC(0.444)
  [orga] NC_000932.1 GC=36%: ACGT(0.486) TTAA(0.628) GTAC(0.642)

palindrome_deserts_verified
```

---

## Limitations

- Only 4-mer palindromes are analyzed. Gelfand & Koonin (1997) analyzed 6-mers, which correspond more directly to Type II restriction enzyme recognition sites (e.g., GAATTC for EcoRI). 4-mers are used here for computational tractability with large genomes.
- Expected counts use a mononucleotide independence model. This over-estimates expected counts for CpG-containing palindromes in high-GC genomes where dinucleotide composition is non-random.
- The panel of 5 bacterial genomes is small. A broader phylogenetic survey would be needed to establish whether avoidance strength correlates with the number of R-M systems encoded in each genome.
- The organellar control is not perfect: mitochondria and chloroplasts have distinct evolutionary histories and compositional biases beyond the absence of R-M systems.
- O/E ratios are not normally distributed; the comparison uses simple means rather than a formal statistical test.