---
name: chembl-drug-readiness-audit
description: >
  Given a ChEMBL target ID and activity threshold, download curated actives,
  calculate ADMET and physicochemical properties, and produce a drug-discovery
  readiness report with scaffold analysis, activity cliff detection, chemical
  space visualisation, and ADMET bottleneck analysis. Outputs a self-contained
  HTML report and machine-readable CSVs.
allowed-tools: Bash(pip *), Bash(python *), Bash(curl *), Bash(echo *), Bash(cat *)
---

# Drug Discovery Readiness Audit via ChEMBL

## Parameters

All parameters are declared here. Change only this section to re-run on a
different target.

```python
# params.py  — written by Step 1, consumed by all later steps
TARGET_CHEMBL_ID = "CHEMBL279"        # EGFR; swap for any ChEMBL target ID
ACTIVITY_TYPE    = "IC50"             # IC50 | Ki | Kd
ACTIVITY_UNIT    = "nM"
ACTIVE_THRESHOLD = 1000               # nM — compounds with value <= this are "active"
CHEMBL_VERSION   = 34                 # pin the release for reproducibility
RANDOM_SEED      = 42
OUTPUT_DIR       = "audit_output"
```

---

## Step 1 — Environment Setup

**Expected time:** ~3 minutes on a clean environment.

```bash
pip install chembl-webresource-client==0.10.9 \
            datamol==0.12.3 \
            rdkit==2024.3.1 \
            admet-ai==1.4.0 \
            umap-learn==0.5.6 \
            pandas==2.1.4 \
            numpy==1.26.4 \
            matplotlib==3.9.0 \
            seaborn==0.13.2 \
            upsetplot==0.9.0 \
            tqdm==4.66.4 \
            jinja2==3.1.4

mkdir -p audit_output/figures audit_output/data
```

Write `params.py` to disk so all later steps share a single config:

```bash
cat > params.py << 'EOF'
TARGET_CHEMBL_ID = "CHEMBL279"
ACTIVITY_TYPE    = "IC50"
ACTIVITY_UNIT    = "nM"
ACTIVE_THRESHOLD = 1000
CHEMBL_VERSION   = 34
RANDOM_SEED      = 42
OUTPUT_DIR       = "audit_output"
EOF
```

> **Important:** All scripts must be run from the **project root** (the
> directory containing `params.py`), not from inside `scripts/`. Every
> command in this skill uses the form `python scripts/NN_name.py` and
> assumes the working directory is the project root. Each script also
> inserts the project root into `sys.path` explicitly as a safeguard —
> see the `sys.path.insert` line at the top of each script below.

**Validation:** `python -c "import datamol, chembl_webresource_client, umap; print('OK')"` must print `OK`.

---

## Step 2 — Data Download

Downloads raw bioactivity records from ChEMBL for the configured target.

```python
# scripts/01_download.py
import sys, os
sys.path.insert(0, os.path.dirname(os.path.dirname(os.path.abspath(__file__))))
import pandas as pd
from chembl_webresource_client.new_client import new_client
from params import TARGET_CHEMBL_ID, ACTIVITY_TYPE, ACTIVITY_UNIT, OUTPUT_DIR

activity_api = new_client.activity
records = activity_api.filter(
    target_chembl_id=TARGET_CHEMBL_ID,
    standard_type=ACTIVITY_TYPE,
    standard_units=ACTIVITY_UNIT,
).only([
    "molecule_chembl_id",
    "canonical_smiles",
    "standard_value",
    "standard_type",
    "standard_units",
    "assay_chembl_id",
    "document_chembl_id",
])

df = pd.DataFrame(list(records))
raw_path = f"{OUTPUT_DIR}/data/raw_activities.csv"
df.to_csv(raw_path, index=False)

print(f"Downloaded {len(df)} records.")
print(f"Columns: {list(df.columns)}")
print(f"Saved to {raw_path}")
```

```bash
python scripts/01_download.py
```

**Expected output:** A line like `Downloaded 4321 records.`  
**Validation:** `wc -l audit_output/data/raw_activities.csv` should be > 2 (header + at least one row).

---

## Step 3 — Curation

Standardises structures and applies four sequential filters. Drop-off counts are
printed at each stage and saved to `curation_log.json` for the report.

```python
# scripts/02_curate.py
import sys, os
sys.path.insert(0, os.path.dirname(os.path.dirname(os.path.abspath(__file__))))
import json, hashlib
import pandas as pd
import datamol as dm
from rdkit.Chem.FilterCatalog import FilterCatalog, FilterCatalogParams
from params import ACTIVE_THRESHOLD, OUTPUT_DIR

df = pd.read_csv(f"{OUTPUT_DIR}/data/raw_activities.csv")
log = {"raw": len(df)}

# 1. Drop rows missing SMILES or activity value
df = df.dropna(subset=["canonical_smiles", "standard_value"])
df["standard_value"] = pd.to_numeric(df["standard_value"], errors="coerce")
df = df.dropna(subset=["standard_value"])
log["after_missing_drop"] = len(df)

# 2. Keep only actives below threshold
df = df[df["standard_value"] <= ACTIVE_THRESHOLD].copy()
log["after_activity_filter"] = len(df)

# 3. Standardise with datamol; drop molecules that fail
def safe_standardize(smi):
    try:
        mol = dm.to_mol(smi, sanitize=True)
        if mol is None:
            return None
        mol = dm.standardize_mol(mol)
        mol = dm.fix_mol(mol)
        return dm.to_smiles(mol)
    except Exception:
        return None

try:
    dm.disable_logs()
except AttributeError:
    pass  # not available in all datamol versions
df["std_smiles"] = df["canonical_smiles"].apply(safe_standardize)
df = df.dropna(subset=["std_smiles"])
log["after_standardization"] = len(df)

# 4. Deduplicate on canonical SMILES, keep lowest IC50 per molecule
df = df.sort_values("standard_value")
df = df.drop_duplicates(subset="std_smiles", keep="first").reset_index(drop=True)
log["after_deduplication"] = len(df)

# 5. PAINS filter
params = FilterCatalogParams()
params.AddCatalog(FilterCatalogParams.FilterCatalogs.PAINS)
catalog = FilterCatalog(params)

def is_pains(smi):
    mol = dm.to_mol(smi)
    return mol is not None and catalog.HasMatch(mol)

df["is_pains"] = df["std_smiles"].apply(is_pains)
log["pains_flagged"] = int(df["is_pains"].sum())
df_clean = df[~df["is_pains"]].copy().reset_index(drop=True)
log["after_pains"] = len(df_clean)

# Save
df_clean.to_csv(f"{OUTPUT_DIR}/data/curated_actives.csv", index=False)
with open(f"{OUTPUT_DIR}/data/curation_log.json", "w") as f:
    json.dump(log, f, indent=2)

print(json.dumps(log, indent=2))
```

```bash
python scripts/02_curate.py
```

**Expected output:** JSON showing progressive drop-off across curation stages.  
**Validation:** `audit_output/data/curated_actives.csv` must exist and have > 0 rows.

---

## Step 4 — Property Calculation

Calculates RDKit 2D descriptors, Lipinski / Veber / QED rule-based flags, and
ADMET-AI predictions. All outputs are merged into one enriched CSV.

```python
# scripts/03_properties.py
import sys, os, warnings
sys.path.insert(0, os.path.dirname(os.path.dirname(os.path.abspath(__file__))))
warnings.filterwarnings("ignore")
import pandas as pd
import numpy as np
import datamol as dm
from rdkit.Chem import Descriptors, QED
from admet_ai import ADMETModel
from params import OUTPUT_DIR

df = pd.read_csv(f"{OUTPUT_DIR}/data/curated_actives.csv")

# --- 2D descriptors ---
def calc_descriptors(smi):
    mol = dm.to_mol(smi)
    if mol is None:
        return {}
    return {
        "MW":       Descriptors.MolWt(mol),
        "LogP":     Descriptors.MolLogP(mol),
        "TPSA":     Descriptors.TPSA(mol),
        "HBD":      Descriptors.NumHDonors(mol),
        "HBA":      Descriptors.NumHAcceptors(mol),
        "RotBonds": Descriptors.NumRotatableBonds(mol),
        "Fsp3":     Descriptors.FractionCSP3(mol),
        "QED":      QED.qed(mol),
        "RingCount": Descriptors.RingCount(mol),
    }

desc = df["std_smiles"].apply(calc_descriptors).apply(pd.Series)
df = pd.concat([df, desc], axis=1)

# --- Rule-based filters ---
df["Ro5_pass"] = (
    (df["MW"]   <= 500) &
    (df["LogP"] <=   5) &
    (df["HBD"]  <=   5) &
    (df["HBA"]  <=  10)
)
df["Veber_pass"] = (df["TPSA"] <= 140) & (df["RotBonds"] <= 10)
df["LeadLike_pass"] = (
    (df["MW"]   <= 350) &
    (df["LogP"] <=   3) &
    (df["RingCount"] <= 4) &
    (df["RotBonds"]  <= 7)
)

# --- ADMET-AI ---
model = ADMETModel()
admet_preds = model.predict(smiles=df["std_smiles"].tolist())
admet_df = pd.DataFrame(admet_preds).reset_index(drop=True)

# Binary flags for key endpoints (model-specific threshold)
ADMET_ENDPOINTS = {
    "HIA_Hou":            ("HIA_pass",    0.5, ">="),   # Human intestinal absorption
    "hERG":               ("hERG_safe",   0.5, "<"),    # Low hERG risk
    "BBB_Martini":        ("BBB_pass",    0.5, ">="),   # BBB penetration
    "CYP3A4_Substrate":   ("CYP3A4_flag", 0.5, "<"),    # Not a CYP3A4 substrate
    "Caco2_Wang":         ("Caco2_pass",  0.5, ">="),   # Permeability
}
for col, (flag, thr, op) in ADMET_ENDPOINTS.items():
    if col in admet_df.columns:
        if op == ">=":
            admet_df[flag] = (admet_df[col] >= thr)
        else:
            admet_df[flag] = (admet_df[col] < thr)

df = pd.concat([df, admet_df], axis=1)
out = f"{OUTPUT_DIR}/data/enriched_actives.csv"
df.to_csv(out, index=False)
print(f"Saved enriched dataset: {len(df)} compounds → {out}")
print(df[["MW","LogP","QED","Ro5_pass","Veber_pass"]].describe().round(2))
```

```bash
python scripts/03_properties.py
```

**Validation:** `audit_output/data/enriched_actives.csv` must contain columns `MW`, `LogP`, `QED`, `Ro5_pass`, `hERG_safe`.

---

## Step 5 — Chemical Space Visualisation (UMAP)

Projects Morgan fingerprints to 2D via UMAP, coloured by pIC50 and by Ro5
pass/fail status.

```python
# scripts/04_umap.py
import sys, os, warnings
sys.path.insert(0, os.path.dirname(os.path.dirname(os.path.abspath(__file__))))
warnings.filterwarnings("ignore")
import matplotlib
matplotlib.use("Agg")
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import umap
from rdkit.Chem import AllChem
import datamol as dm
from params import OUTPUT_DIR, RANDOM_SEED

df = pd.read_csv(f"{OUTPUT_DIR}/data/enriched_actives.csv")

def morgan_fp(smi, radius=2, nbits=2048):
    mol = dm.to_mol(smi)
    if mol is None:
        return np.zeros(nbits)
    return np.array(AllChem.GetMorganFingerprintAsBitVect(mol, radius, nbits))

print("Computing fingerprints...")
fps = np.vstack(df["std_smiles"].apply(morgan_fp).values)

print("Running UMAP...")
reducer = umap.UMAP(n_components=2, random_state=RANDOM_SEED, n_neighbors=15, min_dist=0.1)
coords = reducer.fit_transform(fps)
df["UMAP_1"], df["UMAP_2"] = coords[:, 0], coords[:, 1]
df["pIC50"] = -np.log10(df["standard_value"] * 1e-9)

fig, axes = plt.subplots(1, 2, figsize=(14, 5))

sc = axes[0].scatter(df["UMAP_1"], df["UMAP_2"], c=df["pIC50"],
                     cmap="plasma", s=12, alpha=0.7)
plt.colorbar(sc, ax=axes[0], label="pIC50")
axes[0].set_title("Chemical Space — coloured by pIC50")
axes[0].set_xlabel("UMAP 1"); axes[0].set_ylabel("UMAP 2")

colours = df["Ro5_pass"].map({True: "#2196F3", False: "#F44336"})
axes[1].scatter(df["UMAP_1"], df["UMAP_2"], c=colours, s=12, alpha=0.7)
axes[1].set_title("Chemical Space — Ro5 pass (blue) / fail (red)")
axes[1].set_xlabel("UMAP 1"); axes[1].set_ylabel("UMAP 2")

plt.tight_layout()
plt.savefig(f"{OUTPUT_DIR}/figures/umap.png", dpi=150)
print(f"Saved {OUTPUT_DIR}/figures/umap.png")

df[["molecule_chembl_id","std_smiles","pIC50","UMAP_1","UMAP_2"]].to_csv(
    f"{OUTPUT_DIR}/data/umap_coords.csv", index=False)
```

```bash
python scripts/04_umap.py
```

**Validation:** `audit_output/figures/umap.png` must exist and be > 10 KB.

---

## Step 6 — Scaffold Analysis

Bemis-Murcko scaffold decomposition: top-10 scaffolds by frequency, scaffold
diversity index, and concentration of activity in top scaffold.

```python
# scripts/05_scaffolds.py
import sys, os, warnings
sys.path.insert(0, os.path.dirname(os.path.dirname(os.path.abspath(__file__))))
warnings.filterwarnings("ignore")
import matplotlib
matplotlib.use("Agg")
import pandas as pd
import matplotlib.pyplot as plt
import datamol as dm
from rdkit.Chem.Scaffolds.MurckoScaffold import GetScaffoldForMol
from rdkit.Chem import MolToSmiles
from params import OUTPUT_DIR

df = pd.read_csv(f"{OUTPUT_DIR}/data/enriched_actives.csv")

def get_scaffold(smi):
    mol = dm.to_mol(smi)
    if mol is None:
        return None
    try:
        scaffold = GetScaffoldForMol(mol)
        return MolToSmiles(scaffold) if scaffold else None
    except Exception:
        return None

df["scaffold"] = df["std_smiles"].apply(get_scaffold)
df = df.dropna(subset=["scaffold"])

scaffold_counts = df["scaffold"].value_counts()
n_unique   = len(scaffold_counts)
n_total    = len(df)
diversity  = n_unique / n_total
top1_pct   = scaffold_counts.iloc[0] / n_total * 100
top10_pct  = scaffold_counts.iloc[:10].sum() / n_total * 100

print(f"Total compounds:    {n_total}")
print(f"Unique scaffolds:   {n_unique}")
print(f"Scaffold diversity: {diversity:.3f}  (1.0 = all unique)")
print(f"Top scaffold share: {top1_pct:.1f}%")
print(f"Top-10 scaffolds:   {top10_pct:.1f}% of all actives")

top10 = scaffold_counts.head(10).reset_index()
top10.columns = ["scaffold_smiles", "count"]
top10["pct"] = (top10["count"] / n_total * 100).round(1)
top10.to_csv(f"{OUTPUT_DIR}/data/top_scaffolds.csv", index=False)

fig, ax = plt.subplots(figsize=(8, 4))
ax.barh(range(10), top10["count"], color="#1976D2")
ax.set_yticks(range(10))
ax.set_yticklabels([f"Scaffold {i+1}" for i in range(10)], fontsize=9)
ax.invert_yaxis()
ax.set_xlabel("Number of actives")
ax.set_title("Top 10 Bemis-Murcko Scaffolds")
plt.tight_layout()
plt.savefig(f"{OUTPUT_DIR}/figures/scaffolds.png", dpi=150)
print(f"Saved {OUTPUT_DIR}/figures/scaffolds.png")

import json
with open(f"{OUTPUT_DIR}/data/scaffold_stats.json", "w") as f:
    json.dump({"n_total": n_total, "n_unique_scaffolds": n_unique,
               "diversity_index": round(diversity, 4),
               "top1_pct": round(top1_pct, 2),
               "top10_pct": round(top10_pct, 2)}, f, indent=2)
```

```bash
python scripts/05_scaffolds.py
```

**Validation:** `audit_output/data/scaffold_stats.json` must contain `diversity_index`.

---

## Step 7 — Activity Cliff Detection

Pairwise Tanimoto scan: finds compound pairs with high structural similarity
(≥ 0.85) but large activity difference (ΔpIC50 ≥ 2, i.e. 100× in IC50).

> **Note:** Cliff detection is capped at `MAX_CLIFF_COMPOUNDS = 500` for
> computational tractability — the O(N²) pairwise loop over the full dataset
> is intractable on large ChEMBL targets (e.g. ~31M pairs for EGFR).
> Analysis on a 500-compound random sample is standard practice and sufficient
> to characterise cliff density. Expected runtime: < 2 minutes.

```python
# scripts/06_cliffs.py
import sys, os, warnings
sys.path.insert(0, os.path.dirname(os.path.dirname(os.path.abspath(__file__))))
warnings.filterwarnings("ignore")
import matplotlib
matplotlib.use("Agg")
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from rdkit.Chem import AllChem, DataStructs
import datamol as dm
from params import OUTPUT_DIR, RANDOM_SEED

df = pd.read_csv(f"{OUTPUT_DIR}/data/enriched_actives.csv")
df["pIC50"] = -np.log10(df["standard_value"] * 1e-9)

def morgan_fp(smi):
    mol = dm.to_mol(smi)
    if mol is None:
        return None
    return AllChem.GetMorganFingerprintAsBitVect(mol, 2, 2048)

fps = df["std_smiles"].apply(morgan_fp).tolist()
valid = [(i, fp) for i, fp in enumerate(fps) if fp is not None]

# Cap both cliff detection and scatter plot at MAX_CLIFF_COMPOUNDS
MAX_CLIFF_COMPOUNDS = 500
rng = np.random.default_rng(RANDOM_SEED)
if len(valid) > MAX_CLIFF_COMPOUNDS:
    indices = rng.choice(len(valid), MAX_CLIFF_COMPOUNDS, replace=False)
    valid = [valid[k] for k in sorted(indices)]
print(f"Running cliff detection on {len(valid)} compounds (capped at {MAX_CLIFF_COMPOUNDS})")

SIM_THRESHOLD   = 0.85
DELTA_THRESHOLD = 2.0

cliffs = []
for i in range(len(valid)):
    for j in range(i + 1, len(valid)):
        idx_i, fp_i = valid[i]
        idx_j, fp_j = valid[j]
        sim   = DataStructs.TanimotoSimilarity(fp_i, fp_j)
        delta = abs(df.iloc[idx_i]["pIC50"] - df.iloc[idx_j]["pIC50"])
        if sim >= SIM_THRESHOLD and delta >= DELTA_THRESHOLD:
            cliffs.append({
                "chembl_id_A":  df.iloc[idx_i]["molecule_chembl_id"],
                "smiles_A":     df.iloc[idx_i]["std_smiles"],
                "pIC50_A":      round(df.iloc[idx_i]["pIC50"], 2),
                "chembl_id_B":  df.iloc[idx_j]["molecule_chembl_id"],
                "smiles_B":     df.iloc[idx_j]["std_smiles"],
                "pIC50_B":      round(df.iloc[idx_j]["pIC50"], 2),
                "tanimoto":     round(sim, 3),
                "delta_pIC50":  round(delta, 2),
            })

cliff_df = pd.DataFrame(cliffs)
if len(cliff_df):
    cliff_df = cliff_df.sort_values("delta_pIC50", ascending=False)
cliff_df.to_csv(f"{OUTPUT_DIR}/data/activity_cliffs.csv", index=False)
print(f"Activity cliffs detected: {len(cliff_df)}")
if len(cliff_df):
    print(cliff_df[["chembl_id_A","chembl_id_B","tanimoto","delta_pIC50"]].head(5).to_string())

# Scatter: all pairs within the same capped set
all_pairs = []
for i in range(len(valid)):
    for j in range(i + 1, len(valid)):
        idx_i, fp_i = valid[i]; idx_j, fp_j = valid[j]
        all_pairs.append((
            DataStructs.TanimotoSimilarity(fp_i, fp_j),
            abs(df.iloc[idx_i]["pIC50"] - df.iloc[idx_j]["pIC50"])
        ))
sims, deltas = zip(*all_pairs)
fig, ax = plt.subplots(figsize=(7, 5))
ax.scatter(sims, deltas, s=4, alpha=0.3, color="#90A4AE")
ax.axvline(SIM_THRESHOLD,   color="red",  ls="--", lw=1, label=f"Sim ≥ {SIM_THRESHOLD}")
ax.axhline(DELTA_THRESHOLD, color="blue", ls="--", lw=1, label=f"ΔpIC50 ≥ {DELTA_THRESHOLD}")
ax.set_xlabel("Tanimoto Similarity"); ax.set_ylabel("ΔpIC50")
ax.set_title(f"Activity Cliff Map (sample of {len(valid)} compounds)")
ax.legend()
plt.tight_layout()
plt.savefig(f"{OUTPUT_DIR}/figures/activity_cliffs.png", dpi=150)
print(f"Saved {OUTPUT_DIR}/figures/activity_cliffs.png")
```

```bash
python scripts/06_cliffs.py
```

**Validation:** `audit_output/data/activity_cliffs.csv` must exist (may be empty if no cliffs found — that itself is a finding).

---

## Step 8 — ADMET Bottleneck Analysis

Computes pass rates for each filter independently and as intersections.
Produces an UpSet plot and a summary table.

```python
# scripts/07_admet_bottleneck.py
import sys, os, warnings
sys.path.insert(0, os.path.dirname(os.path.dirname(os.path.abspath(__file__))))
warnings.filterwarnings("ignore")
import matplotlib
matplotlib.use("Agg")
import pandas as pd
import matplotlib.pyplot as plt
import json
from upsetplot import UpSet, from_indicators
from params import OUTPUT_DIR

df = pd.read_csv(f"{OUTPUT_DIR}/data/enriched_actives.csv")

FILTER_COLS = {
    "Ro5":       "Ro5_pass",
    "Veber":     "Veber_pass",
    "Lead-like": "LeadLike_pass",
    "hERG safe": "hERG_safe",
    "HIA":       "HIA_pass",
}

available = {k: v for k, v in FILTER_COLS.items() if v in df.columns}
filter_df = df[[v for v in available.values()]].copy()
filter_df.columns = list(available.keys())
filter_df = filter_df.fillna(False).astype(bool)

pass_rates = {k: round(filter_df[k].mean() * 100, 1) for k in filter_df.columns}
print("Individual pass rates (%):")
for k, v in pass_rates.items():
    print(f"  {k}: {v}%")

all_pass = filter_df.all(axis=1)
print(f"\nPass ALL filters: {all_pass.sum()} / {len(df)} ({all_pass.mean()*100:.1f}%)")

# UpSet plot
upset_data = from_indicators(filter_df.columns.tolist(), data=filter_df)
fig = plt.figure(figsize=(10, 5))
UpSet(upset_data, subset_size="count", show_counts=True).plot(fig)
plt.suptitle("ADMET Filter Intersections", y=1.02)
plt.savefig(f"{OUTPUT_DIR}/figures/admet_upset.png", dpi=150, bbox_inches="tight")
print(f"Saved {OUTPUT_DIR}/figures/admet_upset.png")

with open(f"{OUTPUT_DIR}/data/admet_bottleneck.json", "w") as f:
    json.dump({"pass_rates_pct": pass_rates,
               "pass_all_count": int(all_pass.sum()),
               "pass_all_pct":   round(all_pass.mean()*100, 1),
               "total":          len(df)}, f, indent=2)
```

```bash
python scripts/07_admet_bottleneck.py
```

**Validation:** `audit_output/data/admet_bottleneck.json` must contain `pass_all_pct`.

---

## Step 9 — Compile Report

Generates a self-contained HTML report from all outputs. The agent
fills in the narrative template with computed statistics.

```python
# scripts/08_report.py
import sys, os
sys.path.insert(0, os.path.dirname(os.path.dirname(os.path.abspath(__file__))))
import json, base64, pathlib
import pandas as pd
from jinja2 import Template
from params import TARGET_CHEMBL_ID, ACTIVITY_TYPE, ACTIVE_THRESHOLD, OUTPUT_DIR

def img_b64(path):
    data = pathlib.Path(path).read_bytes()
    return base64.b64encode(data).decode()

cur  = json.load(open(f"{OUTPUT_DIR}/data/curation_log.json"))
scaf = json.load(open(f"{OUTPUT_DIR}/data/scaffold_stats.json"))
admt = json.load(open(f"{OUTPUT_DIR}/data/admet_bottleneck.json"))
try:
    cliff_df = pd.read_csv(f"{OUTPUT_DIR}/data/activity_cliffs.csv")
except (pd.errors.EmptyDataError, FileNotFoundError):
    cliff_df = pd.DataFrame(columns=["chembl_id_A","chembl_id_B","tanimoto","delta_pIC50"])

TEMPLATE = """
<!DOCTYPE html><html><head><meta charset="utf-8">
<title>Drug Readiness Audit — {{ target }}</title>
<style>
  body{font-family:Georgia,serif;max-width:900px;margin:40px auto;line-height:1.6;color:#222}
  h1{color:#1a237e} h2{color:#283593;border-bottom:1px solid #ccc;padding-bottom:4px}
  table{border-collapse:collapse;width:100%} th,td{border:1px solid #ddd;padding:8px}
  th{background:#e8eaf6} img{max-width:100%;border:1px solid #eee;border-radius:4px}
  .stat{font-size:1.8em;font-weight:bold;color:#1565c0}
  .label{font-size:0.85em;color:#555}
  .grid{display:grid;grid-template-columns:1fr 1fr 1fr;gap:16px;margin:20px 0}
  .card{background:#f5f5f5;border-radius:6px;padding:16px;text-align:center}
</style></head><body>
<h1>Drug Discovery Readiness Audit</h1>
<p><strong>Target:</strong> {{ target }} &nbsp;|&nbsp;
   <strong>Activity:</strong> {{ act_type }} ≤ {{ threshold }} nM</p>

<div class="grid">
  <div class="card"><div class="stat">{{ cur.raw }}</div><div class="label">Raw records</div></div>
  <div class="card"><div class="stat">{{ cur.after_pains }}</div><div class="label">Curated actives</div></div>
  <div class="card"><div class="stat">{{ admt.pass_all_pct }}%</div><div class="label">Pass all ADMET filters</div></div>
</div>

<h2>1. Curation Funnel</h2>
<table>
  <tr><th>Stage</th><th>Count</th><th>Dropped</th></tr>
  <tr><td>Raw download</td><td>{{ cur.raw }}</td><td>—</td></tr>
  <tr><td>After missing-value drop</td><td>{{ cur.after_missing_drop }}</td>
      <td>{{ cur.raw - cur.after_missing_drop }}</td></tr>
  <tr><td>After activity threshold</td><td>{{ cur.after_activity_filter }}</td>
      <td>{{ cur.after_missing_drop - cur.after_activity_filter }}</td></tr>
  <tr><td>After standardisation</td><td>{{ cur.after_standardization }}</td>
      <td>{{ cur.after_activity_filter - cur.after_standardization }}</td></tr>
  <tr><td>After deduplication</td><td>{{ cur.after_deduplication }}</td>
      <td>{{ cur.after_standardization - cur.after_deduplication }}</td></tr>
  <tr><td><strong>After PAINS filter</strong></td><td><strong>{{ cur.after_pains }}</strong></td>
      <td>{{ cur.pains_flagged }} PAINS removed</td></tr>
</table>

<h2>2. Chemical Space</h2>
<img src="data:image/png;base64,{{ umap_b64 }}" alt="UMAP">

<h2>3. Scaffold Analysis</h2>
<p>The {{ cur.after_pains }} curated actives map to
<strong>{{ scaf.n_unique_scaffolds }} unique Bemis-Murcko scaffolds</strong>
(diversity index: {{ scaf.diversity_index }}).
The most frequent scaffold accounts for <strong>{{ scaf.top1_pct }}%</strong>
of all actives; the top 10 scaffolds together cover <strong>{{ scaf.top10_pct }}%</strong>.</p>
<img src="data:image/png;base64,{{ scaf_b64 }}" alt="Scaffolds">

<h2>4. Activity Cliffs</h2>
<p>Detected <strong>{{ n_cliffs }} activity cliff pairs</strong>
(Tanimoto ≥ 0.85, ΔpIC50 ≥ 2.0).</p>
<img src="data:image/png;base64,{{ cliff_b64 }}" alt="Activity Cliffs">
{% if top_cliffs %}
<table>
  <tr><th>Compound A</th><th>Compound B</th><th>Tanimoto</th><th>ΔpIC50</th></tr>
  {% for r in top_cliffs %}
  <tr><td>{{ r.chembl_id_A }}</td><td>{{ r.chembl_id_B }}</td>
      <td>{{ r.tanimoto }}</td><td>{{ r.delta_pIC50 }}</td></tr>
  {% endfor %}
</table>
{% endif %}

<h2>5. ADMET Bottleneck</h2>
<img src="data:image/png;base64,{{ admt_b64 }}" alt="ADMET UpSet">
<table>
  <tr><th>Filter</th><th>Pass rate (%)</th></tr>
  {% for k,v in admt.pass_rates_pct.items() %}
  <tr><td>{{ k }}</td><td>{{ v }}%</td></tr>
  {% endfor %}
  <tr><td><strong>All filters</strong></td><td><strong>{{ admt.pass_all_pct }}%</strong></td></tr>
</table>

<h2>6. Conclusions</h2>
<p>Of {{ cur.raw }} raw ChEMBL records for {{ target }},
{{ cur.after_pains }} compounds survived full curation
({{ "%.1f"|format(cur.after_pains / cur.raw * 100) }}% retention).
The active chemical matter shows a scaffold diversity index of {{ scaf.diversity_index }},
indicating {% if scaf.diversity_index > 0.5 %}broad{% else %}narrow{% endif %} structural coverage.
Only {{ admt.pass_all_pct }}% of actives pass all five ADMET filters simultaneously,
with {{ admt.pass_rates_pct | dictsort | first | last }}% being the most restrictive individual gate.
{{ n_cliffs }} activity cliffs were identified, highlighting regions of high SAR sensitivity
that warrant focused medicinal chemistry attention.</p>
</body></html>
"""

html = Template(TEMPLATE).render(
    target     = TARGET_CHEMBL_ID,
    act_type   = ACTIVITY_TYPE,
    threshold  = ACTIVE_THRESHOLD,
    cur        = cur,
    scaf       = scaf,
    admt       = admt,
    n_cliffs   = len(cliff_df),
    top_cliffs = cliff_df.head(5).to_dict("records") if len(cliff_df) else [],
    umap_b64   = img_b64(f"{OUTPUT_DIR}/figures/umap.png"),
    scaf_b64   = img_b64(f"{OUTPUT_DIR}/figures/scaffolds.png"),
    cliff_b64  = img_b64(f"{OUTPUT_DIR}/figures/activity_cliffs.png"),
    admt_b64   = img_b64(f"{OUTPUT_DIR}/figures/admet_upset.png"),
)

out = f"{OUTPUT_DIR}/report.html"
pathlib.Path(out).write_text(html, encoding="utf-8")
print(f"Report saved to {out}")
```

```bash
python scripts/08_report.py
```

**Expected final output:** `audit_output/report.html` — a self-contained HTML file with all figures inlined as base64. Open in any browser.

---

## Expected Outputs

```
audit_output/
├── data/
│   ├── raw_activities.csv        # ChEMBL download
│   ├── curated_actives.csv       # After curation
│   ├── enriched_actives.csv      # + descriptors + ADMET
│   ├── umap_coords.csv
│   ├── top_scaffolds.csv
│   ├── activity_cliffs.csv
│   ├── curation_log.json
│   ├── scaffold_stats.json
│   └── admet_bottleneck.json
├── figures/
│   ├── umap.png
│   ├── scaffolds.png
│   ├── activity_cliffs.png
│   └── admet_upset.png
└── report.html                   ← primary deliverable
```

---

## Reproducing on a Different Target

1. Edit `params.py`: change `TARGET_CHEMBL_ID` (e.g. `"CHEMBL3227"` for BRAF)
2. Re-run steps 2–9 in order
3. All outputs regenerate automatically in the same `audit_output/` directory