---
name: protein-druggability-dossier
description: >
  Given a PDB ID, produces a comprehensive druggability dossier by running
  fpocket for pocket geometry, benchmarking against a curated druggable/undruggable
  reference set, querying ChEMBL for bioactivity evidence, analysing binding site
  amino acid composition, computing pocket residue B-factor flexibility, and
  assessing pocket conservation across homologous structures. Outputs a composite
  druggability confidence score and a self-contained HTML report.
allowed-tools: Bash(pip *), Bash(conda *), Bash(python *), Bash(curl *), Bash(wget *), Bash(echo *), Bash(cat *), Bash(mkdir *), Bash(find *), Bash(ls *)
---

# Protein Pocket Druggability Dossier

## Parameters

Edit only this section to run on a different target.

```bash
cat > params.py << 'EOF'
# Target structure
PDB_ID = "2ITY"   # EGFR kinase domain — swap for any 4-letter PDB ID

# fpocket druggability thresholds (literature-derived)
DRUGGABILITY_SCORE_MIN = 0.5   # fpocket drug_score >= this → druggable
POCKET_VOLUME_MIN      = 300   # ų — minimum meaningful pocket volume
TOP_N_POCKETS          = 5     # how many pockets to report in detail

# Coverage analysis
TANIMOTO_THRESHOLD = 0.4       # for ChEMBL nearest-neighbour lookup

# Reference set (PDB IDs with known approved drugs — druggable ground truth)
DRUGGABLE_REFS = [
    "2ITY",  # EGFR kinase
    "3ERT",  # Estrogen receptor alpha
    "1IEP",  # BCR-ABL (imatinib)
    "1YVJ",  # JAK3 kinase
    "4DJV",  # BRAF kinase
    "2OIQ",  # Src kinase
    "1S9D",  # VEGFR2 kinase
    "3O96",  # AKT1 kinase
]

# Known undruggable / difficult targets
UNDRUGGABLE_REFS = [
    "2MZD",  # p53+CBP complex (historically difficult)
    "1TF6",  # phage repressor (no human ChEMBL data — use as geometry-only ref)
    "1A1T",  # MYC bHLH domain (canonical undruggable)
    "2LZK",  # BCL6 BTB domain (NMR — protein-protein interface)
    "1CE8",  # CPS1 (bacterial — no ChEMBL data)
]

RANDOM_SEED = 42
OUTPUT_DIR  = "druggability_output"
EOF
```

> **Important:** Run all commands from the **project root** (directory containing
> `params.py`). Use `python scripts/NN_name.py` — never `cd scripts/`.
> Each script inserts the project root into `sys.path` explicitly.

---

## Step 1 — Environment Setup

**Expected time:** ~5 minutes. fpocket is installed via conda.

```bash
# Python packages
pip install biopython==1.83 \
            requests==2.31.0 \
            pandas==2.1.4 \
            numpy==1.26.4 \
            matplotlib==3.9.0 \
            seaborn==0.13.2 \
            scipy==1.13.0 \
            chembl-webresource-client==0.10.9 \
            jinja2==3.1.4

# fpocket via conda — use conda-forge (bioconda builds are outdated)
conda install -y -c conda-forge fpocket

mkdir -p druggability_output/structures \
         druggability_output/fpocket \
         druggability_output/data \
         druggability_output/figures \
         scripts
```

**Validation:**
```bash
python -c "import Bio, requests, pandas, chembl_webresource_client; print('Python OK')"
fpocket --version
```
Both must succeed without error.

---

## Step 2 — Download PDB Structure(s)

Downloads the query structure and all reference structures (druggable + undruggable)
from the RCSB PDB. Also fetches all available structures of the query protein
for multi-structure pocket stability analysis (Step 8).

```python
# scripts/01_download_structures.py
import sys, os, requests, time
sys.path.insert(0, os.path.dirname(os.path.dirname(os.path.abspath(__file__))))
from params import PDB_ID, DRUGGABLE_REFS, UNDRUGGABLE_REFS, OUTPUT_DIR

ALL_IDS = list({PDB_ID} | set(DRUGGABLE_REFS) | set(UNDRUGGABLE_REFS))

def download_pdb(pdb_id, out_dir):
    pdb_id = pdb_id.upper()
    out_path = os.path.join(out_dir, f"{pdb_id}.pdb")
    if os.path.exists(out_path):
        print(f"  {pdb_id} already cached")
        return out_path
    url = f"https://files.rcsb.org/download/{pdb_id}.pdb"
    resp = requests.get(url, timeout=30)
    if resp.status_code != 200:
        print(f"  WARNING: Could not download {pdb_id} (HTTP {resp.status_code})")
        return None
    with open(out_path, "w") as f:
        f.write(resp.text)
    print(f"  Downloaded {pdb_id} → {out_path}")
    return out_path

print(f"Downloading {len(ALL_IDS)} structures ...")
results = {}
for pid in ALL_IDS:
    path = download_pdb(pid, f"{OUTPUT_DIR}/structures")
    results[pid] = path
    time.sleep(0.3)  # be polite to RCSB

# Fetch homologous structures of the query protein via RCSB search API
import json
query = {
    "query": {
        "type": "terminal",
        "service": "sequence",
        "parameters": {
            "evalue_cutoff": 1e-10,
            "identity_cutoff": 0.9,
            "target": "pdb_protein_sequence",
            "value": open(f"{OUTPUT_DIR}/structures/{PDB_ID.upper()}.pdb").read()
        }
    },
    "return_type": "entry",
    "request_options": {"results_slice": {"start": 0, "rows": 20}}
}
# Use simpler entity search instead — query by PDB ID cluster
url = f"https://data.rcsb.org/rest/v1/holdings/entry/{PDB_ID.upper()}"
resp = requests.get(url, timeout=15)
homologs = []
if resp.status_code == 200:
    data = resp.json()
    entity_ids = data.get("rcsb_entry_container_identifiers", {}).get(
        "polymer_entity_ids", [])
    if entity_ids:
        cluster_url = (f"https://data.rcsb.org/rest/v1/holdings/entry/"
                       f"{PDB_ID.upper()}/polymer_entities/{entity_ids[0]}"
                       f"/polymer_entity_instances")
        # Fall back to a simpler strategy: fetch related entries from RCSB
        search_url = ("https://search.rcsb.org/rcsbsearch/v2/query?"
                      f"json=" + requests.utils.quote(json.dumps({
                          "query": {"type": "terminal", "service": "full_text",
                                    "parameters": {"value": PDB_ID.upper()}},
                          "return_type": "entry",
                          "request_options": {"paginate": {"start": 0, "rows": 10}}
                      })))
        sr = requests.get(search_url, timeout=15)
        if sr.status_code == 200:
            hits = sr.json().get("result_set", [])
            homologs = [h["identifier"] for h in hits
                        if h["identifier"].upper() != PDB_ID.upper()][:5]

print(f"\nHomologous structures to compare: {homologs}")
import pickle
with open(f"{OUTPUT_DIR}/data/homologs.pkl", "wb") as f:
    pickle.dump(homologs, f)

# Download homologs
for pid in homologs:
    download_pdb(pid, f"{OUTPUT_DIR}/structures")
    time.sleep(0.3)

print(f"\nAll structures saved to {OUTPUT_DIR}/structures/")
```

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

**Validation:** `ls druggability_output/structures/*.pdb | wc -l` must be ≥ 3.

---

## Step 3 — Run fpocket on All Structures

Runs fpocket on every downloaded PDB file and parses the output into a
structured DataFrame.

```python
# scripts/02_run_fpocket.py
import sys, os, subprocess, re, json
sys.path.insert(0, os.path.dirname(os.path.dirname(os.path.abspath(__file__))))
import pandas as pd
from params import OUTPUT_DIR, PDB_ID, DRUGGABLE_REFS, UNDRUGGABLE_REFS

STRUCT_DIR  = f"{OUTPUT_DIR}/structures"
FPOCKET_DIR = f"{OUTPUT_DIR}/fpocket"
os.makedirs(FPOCKET_DIR, exist_ok=True)

def run_fpocket(pdb_path, out_dir):
    """Run fpocket and return path to output directory.

    fpocket 4.0 ignores the -o flag and always writes output next to the
    input PDB file as {basename}_out/. We move it to out_dir afterward.
    """
    pdb_id   = os.path.basename(pdb_path).replace(".pdb", "")
    # fpocket 4.0 writes here regardless of -o
    auto_out = os.path.join(os.path.dirname(pdb_path), f"{pdb_id}_out")
    dest_out = os.path.join(out_dir, f"{pdb_id}_out")

    if os.path.isdir(dest_out):
        print(f"  {pdb_id}: fpocket output already exists, skipping")
        return dest_out

    cmd = ["fpocket", "-f", pdb_path]
    try:
        subprocess.run(cmd, capture_output=True, timeout=120, check=True,
                       cwd=os.path.dirname(pdb_path))
        print(f"  {pdb_id}: fpocket completed")
    except subprocess.CalledProcessError as e:
        print(f"  {pdb_id}: fpocket failed — {e.stderr.decode()[:200]}")
        return None
    except subprocess.TimeoutExpired:
        print(f"  {pdb_id}: fpocket timed out")
        return None

    # Move output to the designated fpocket directory
    if os.path.isdir(auto_out):
        import shutil
        shutil.move(auto_out, dest_out)
        print(f"  {pdb_id}: moved output → {dest_out}")
    elif os.path.isdir(dest_out):
        pass  # already there
    else:
        print(f"  {pdb_id}: WARNING — could not find fpocket output at {auto_out}")
        return None

    return dest_out

def parse_fpocket_info(result_dir, pdb_id):
    """Parse the fpocket _info.txt file into a list of pocket dicts."""
    info_file = os.path.join(result_dir, f"{pdb_id}_info.txt")
    if not os.path.exists(info_file):
        # try alternate naming
        candidates = [f for f in os.listdir(result_dir) if f.endswith("_info.txt")]
        if not candidates:
            return []
        info_file = os.path.join(result_dir, candidates[0])

    pockets = []
    current = {}
    with open(info_file) as f:
        for line in f:
            line = line.strip()
            if line.startswith("Pocket"):
                if current:
                    pockets.append(current)
                pocket_num = int(re.search(r"Pocket\s+(\d+)", line).group(1))
                current = {"pocket_id": pocket_num, "pdb_id": pdb_id.upper()}
            elif ":" in line:
                key, _, val = line.partition(":")
                key = key.strip().lower().replace(" ", "_").replace("(", "").replace(")", "").replace(".", "")
                try:
                    current[key] = float(val.strip())
                except ValueError:
                    current[key] = val.strip()
    if current:
        pockets.append(current)
    return pockets

all_pockets = []
pdb_files = [f for f in os.listdir(STRUCT_DIR) if f.endswith(".pdb")]
print(f"Running fpocket on {len(pdb_files)} structures ...")

for pdb_file in sorted(pdb_files):
    pdb_path  = os.path.join(STRUCT_DIR, pdb_file)
    pdb_id    = pdb_file.replace(".pdb", "")
    result_dir = run_fpocket(pdb_path, FPOCKET_DIR)
    if result_dir:
        pockets = parse_fpocket_info(result_dir, pdb_id)
        # Label each pocket with its reference category
        if pdb_id.upper() == PDB_ID.upper():
            cat = "query"
        elif pdb_id.upper() in [r.upper() for r in DRUGGABLE_REFS]:
            cat = "druggable_ref"
        elif pdb_id.upper() in [r.upper() for r in UNDRUGGABLE_REFS]:
            cat = "undruggable_ref"
        else:
            cat = "homolog"
        for p in pockets:
            p["category"] = cat
        all_pockets.extend(pockets)
        print(f"  {pdb_id}: {len(pockets)} pockets parsed")

df = pd.DataFrame(all_pockets)
df.to_csv(f"{OUTPUT_DIR}/data/all_pockets.csv", index=False)
print(f"\nTotal pockets across all structures: {len(df)}")
print(f"Columns: {list(df.columns)}")
```

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

**Validation:** `druggability_output/data/all_pockets.csv` must exist and have > 0 rows.
```bash
python -c "
import pandas as pd
df = pd.read_csv('druggability_output/data/all_pockets.csv')
print(f'{len(df)} pockets, {df[\"pdb_id\"].nunique()} structures')
assert len(df) > 0
print('OK')
"
```

---

## ⚠️ Checkpoint — Inspect fpocket Output Before Continuing

**Do not proceed to Step 4 until this checkpoint passes.**

fpocket output directory naming and column names vary by version. Run the
following inspection commands and verify the results match expectations.
If they do not match, fix the paths/column names in all subsequent scripts
before continuing.

```bash
# 1. Show the actual fpocket output directory structure for the query protein
ls druggability_output/fpocket/

# 2. Show contents of the query protein's fpocket output directory
ls druggability_output/fpocket/$(ls druggability_output/fpocket/ | head -1)/

# 3. Show the actual column names in all_pockets.csv
python -c "
import pandas as pd
df = pd.read_csv('druggability_output/data/all_pockets.csv')
print('Columns:', list(df.columns))
print()
# Identify the drug score column
drug_col = next((c for c in df.columns if 'drug' in c.lower()), None)
vol_col  = next((c for c in df.columns if 'volume' in c.lower() or 'vol' in c.lower()), None)
print(f'Drug score column detected: {drug_col}')
print(f'Volume column detected:     {vol_col}')
print()
# Show sample values for the query structure
query_rows = df[df['pdb_id'].str.upper() == df['pdb_id'].str.upper().iloc[0]]
print(f'Sample rows for {query_rows[\"pdb_id\"].iloc[0]}:')
print(query_rows[[c for c in [\"pocket_id\", drug_col, vol_col] if c]].head(5).to_string())
"

# 4. Check pocket PDB files exist and follow expected naming
ls druggability_output/fpocket/$(ls druggability_output/fpocket/ | grep -i $(python -c "from params import PDB_ID; print(PDB_ID.upper())" 2>/dev/null || echo "2ITY") | head -1)/pocket*.pdb 2>/dev/null | head -5 || echo "WARNING: pocket PDB files not found — check fpocket version"
```

**Expected results:**
- The fpocket output directory should be named `{PDB_ID}_out` (e.g. `2ITY_out`)
- Pocket atom files should be named `pocket{N}_atm.pdb` inside that directory
- The drug score column should contain the word `drug` (e.g. `drug_score`)
- The volume column should contain `vol` or `volume` (e.g. `volume`)

**If the directory or file naming differs**, update the following variables at
the top of scripts 05, 06, and 07 before running them:

```python
# If fpocket uses a different output directory name:
fpocket_dir = f"{OUTPUT_DIR}/fpocket/YOUR_ACTUAL_DIR_NAME"

# If fpocket uses a different pocket file naming convention:
pocket_pdb = os.path.join(fpocket_dir, f"YOUR_POCKET_FILE_PREFIX{pocket_id}_atm.pdb")
```

**If the drug score column name differs**, it will be auto-detected by the
fuzzy match (`"drug" in column name`) already in scripts 03, 07, and 08 —
but print and confirm the detected name matches a meaningful score column.

---

## Step 4 — Druggability Benchmarking

Computes a **druggability percentile** for the query protein's top pocket
against the reference set distribution. This contextualises raw fpocket scores.

```python
# scripts/03_benchmark.py
import sys, os, warnings, json
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 matplotlib.patches as mpatches
from scipy import stats
from params import (OUTPUT_DIR, PDB_ID, DRUGGABILITY_SCORE_MIN,
                    POCKET_VOLUME_MIN, TOP_N_POCKETS)

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

# Identify the drug_score column (fpocket naming can vary slightly)
score_col = next((c for c in df.columns if "drug" in c.lower()), None)
vol_col   = next((c for c in df.columns if "volume" in c.lower() or "vol" in c.lower()), None)

if score_col is None:
    raise ValueError(f"Could not find drug score column. Columns: {list(df.columns)}")

print(f"Using score column: '{score_col}', volume column: '{vol_col}'")

# Best pocket per structure (highest drug score)
best = df.groupby("pdb_id")[score_col].max().reset_index()
best.columns = ["pdb_id", "best_drug_score"]
best = best.merge(
    df.groupby("pdb_id")["category"].first().reset_index(), on="pdb_id")

# Query best score
query_score = best.loc[
    best["pdb_id"].str.upper() == PDB_ID.upper(), "best_drug_score"].values

if len(query_score) == 0:
    raise ValueError(f"Query PDB {PDB_ID} not found in results.")
query_score = float(query_score[0])

druggable_scores   = best.loc[best["category"] == "druggable_ref",
                                "best_drug_score"].values
undruggable_scores = best.loc[best["category"] == "undruggable_ref",
                                "best_drug_score"].values

# Percentile relative to druggable refs
if len(druggable_scores) > 0:
    pct = float(stats.percentileofscore(druggable_scores, query_score))
else:
    pct = float("nan")

# Separation between druggable and undruggable refs
if len(druggable_scores) > 0 and len(undruggable_scores) > 0:
    ks_stat, ks_p = stats.ks_2samp(druggable_scores, undruggable_scores)
else:
    ks_stat, ks_p = float("nan"), float("nan")

result = {
    "query_pdb":            PDB_ID.upper(),
    "query_best_drug_score": round(query_score, 4),
    "druggability_percentile_vs_refs": round(pct, 1),
    "verdict": ("DRUGGABLE" if query_score >= DRUGGABILITY_SCORE_MIN
                else "BORDERLINE" if query_score >= DRUGGABILITY_SCORE_MIN * 0.7
                else "CHALLENGING"),
    "ref_druggable_mean":   round(float(np.mean(druggable_scores)), 4) if len(druggable_scores) else None,
    "ref_undruggable_mean": round(float(np.mean(undruggable_scores)), 4) if len(undruggable_scores) else None,
    "ref_separation_ks":    round(ks_stat, 4),
    "ref_separation_p":     round(ks_p, 6),
    "score_column_used":    score_col,
}
print(json.dumps(result, indent=2))

with open(f"{OUTPUT_DIR}/data/benchmark_result.json", "w") as f:
    json.dump(result, f, indent=2)

# Figure: score distribution
fig, ax = plt.subplots(figsize=(9, 5))
if len(druggable_scores):
    ax.hist(druggable_scores, bins=15, alpha=0.65, color="#1565C0",
            label=f"Druggable refs (n={len(druggable_scores)})", density=True)
if len(undruggable_scores):
    ax.hist(undruggable_scores, bins=15, alpha=0.65, color="#E53935",
            label=f"Undruggable refs (n={len(undruggable_scores)})", density=True)
ax.axvline(query_score, color="#2E7D32", lw=2.5, ls="--",
           label=f"Query: {PDB_ID.upper()} ({query_score:.3f})")
ax.axvline(DRUGGABILITY_SCORE_MIN, color="grey", lw=1.5, ls=":",
           label=f"Threshold ({DRUGGABILITY_SCORE_MIN})")
ax.set_xlabel("fpocket Drug Score"); ax.set_ylabel("Density")
ax.set_title(f"Druggability Benchmarking — {PDB_ID.upper()} "
             f"(percentile: {pct:.0f}%)")
ax.legend(fontsize=9)
plt.tight_layout()
plt.savefig(f"{OUTPUT_DIR}/figures/benchmark.png", dpi=150)
print("Saved benchmark.png")

# Figure: all pockets of query ranked
query_pockets = df[df["pdb_id"].str.upper() == PDB_ID.upper()].copy()
query_pockets = query_pockets.sort_values(score_col, ascending=False).head(TOP_N_POCKETS)

fig, ax = plt.subplots(figsize=(8, 4))
colours = ["#1565C0" if s >= DRUGGABILITY_SCORE_MIN else
           "#FFA726" if s >= DRUGGABILITY_SCORE_MIN * 0.7 else
           "#E53935"
           for s in query_pockets[score_col]]
bars = ax.barh(range(len(query_pockets)), query_pockets[score_col],
               color=colours, alpha=0.85)
ax.set_yticks(range(len(query_pockets)))
ax.set_yticklabels([f"Pocket {int(r['pocket_id'])}" for _, r in query_pockets.iterrows()])
ax.invert_yaxis()
ax.axvline(DRUGGABILITY_SCORE_MIN, color="grey", lw=1.5, ls=":", label="Threshold")
ax.set_xlabel("Drug Score"); ax.set_title(f"Top {TOP_N_POCKETS} Pockets — {PDB_ID.upper()}")
ax.legend()
plt.tight_layout()
plt.savefig(f"{OUTPUT_DIR}/figures/pocket_ranking.png", dpi=150)
print("Saved pocket_ranking.png")

query_pockets.to_csv(f"{OUTPUT_DIR}/data/query_top_pockets.csv", index=False)
```

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

**Validation:** `druggability_output/data/benchmark_result.json` must contain `verdict`.

---

## Step 5 — ChEMBL Bioactivity Evidence

Resolves the query PDB ID to a UniProt ID via RCSB, then queries ChEMBL
for known small molecule binders and their potency distribution.

```python
# scripts/04_chembl_evidence.py
import sys, os, warnings, json, time
sys.path.insert(0, os.path.dirname(os.path.dirname(os.path.abspath(__file__))))
warnings.filterwarnings("ignore")
import requests
import pandas as pd
import numpy as np
from chembl_webresource_client.new_client import new_client
from params import PDB_ID, OUTPUT_DIR

# Step A: PDB → UniProt via RCSB API
print(f"Resolving {PDB_ID.upper()} → UniProt ...")
url = f"https://data.rcsb.org/rest/v1/core/entry/{PDB_ID.upper()}"
resp = requests.get(url, timeout=20)
uniprot_ids = []
if resp.status_code == 200:
    entry = resp.json()
    polymer_ids = (entry.get("rcsb_entry_container_identifiers", {})
                        .get("polymer_entity_ids", []))
    for eid in polymer_ids:
        eu = f"https://data.rcsb.org/rest/v1/core/polymer_entity/{PDB_ID.upper()}/{eid}"
        er = requests.get(eu, timeout=15)
        if er.status_code == 200:
            ed = er.json()
            refs = (ed.get("rcsb_polymer_entity_container_identifiers", {})
                      .get("uniprot_ids", []))
            uniprot_ids.extend(refs)
        time.sleep(0.2)

uniprot_ids = list(set(uniprot_ids))
print(f"UniProt IDs found: {uniprot_ids}")

# Step B: UniProt → ChEMBL target
chembl_targets = []
target_api = new_client.target
for uid in uniprot_ids:
    hits = target_api.filter(target_components__accession=uid).only(
        ["target_chembl_id", "pref_name", "target_type"])
    chembl_targets.extend(list(hits))
    time.sleep(0.2)

print(f"ChEMBL targets: {[t['target_chembl_id'] for t in chembl_targets]}")

# Step C: Fetch bioactivity data
activity_api = new_client.activity
all_acts = []
for tgt in chembl_targets:
    tid = tgt["target_chembl_id"]
    acts = activity_api.filter(
        target_chembl_id=tid,
        standard_type__in=["IC50", "Ki", "Kd"],
        standard_units="nM",
    ).only(["molecule_chembl_id", "standard_value", "standard_type",
            "canonical_smiles", "assay_chembl_id"])
    batch = list(acts)
    for a in batch:
        a["chembl_target_id"] = tid
        a["target_name"] = tgt.get("pref_name", "")
    all_acts.extend(batch)
    time.sleep(0.3)

result_meta = {
    "query_pdb":      PDB_ID.upper(),
    "uniprot_ids":    uniprot_ids,
    "chembl_targets": [t["target_chembl_id"] for t in chembl_targets],
    "total_bioactivity_records": len(all_acts),
}

if all_acts:
    df = pd.DataFrame(all_acts)
    df["standard_value"] = pd.to_numeric(df["standard_value"], errors="coerce")
    df = df.dropna(subset=["standard_value"])
    df["pIC50"] = -np.log10(df["standard_value"] * 1e-9)
    df.to_csv(f"{OUTPUT_DIR}/data/chembl_bioactivity.csv", index=False)

    result_meta.update({
        "n_unique_compounds":   int(df["molecule_chembl_id"].nunique()),
        "best_potency_nM":      round(float(df["standard_value"].min()), 2),
        "median_potency_nM":    round(float(df["standard_value"].median()), 2),
        "n_sub_100nM":          int((df["standard_value"] <= 100).sum()),
        "n_sub_1uM":            int((df["standard_value"] <= 1000).sum()),
        "activity_types":       df["standard_type"].value_counts().to_dict(),
        "bioactivity_verdict":  (
            "WELL_CHARACTERISED" if df["molecule_chembl_id"].nunique() > 50 else
            "SOME_EVIDENCE"      if df["molecule_chembl_id"].nunique() > 5  else
            "SPARSE_EVIDENCE"
        ),
    })
else:
    pd.DataFrame().to_csv(f"{OUTPUT_DIR}/data/chembl_bioactivity.csv", index=False)
    result_meta["bioactivity_verdict"] = "NO_EVIDENCE"

print(json.dumps(result_meta, indent=2))
with open(f"{OUTPUT_DIR}/data/chembl_evidence.json", "w") as f:
    json.dump(result_meta, f, indent=2)
```

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

**Validation:** `druggability_output/data/chembl_evidence.json` must contain `bioactivity_verdict`.

---

## Step 6 — Binding Site Amino Acid Composition

Parses the PDB file and fpocket pocket residue lists to characterise the
amino acid composition of each top pocket. Computes hydrophobicity fraction,
aromaticity, and charge — all predictors of druggability.

```python
# scripts/05_pocket_composition.py
import sys, os, warnings, json
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 Bio import PDB
from params import PDB_ID, OUTPUT_DIR, TOP_N_POCKETS

# Amino acid property lookup
AA_PROPS = {
    # (hydrophobic, aromatic, positively_charged, negatively_charged, polar)
    "ALA": (1,0,0,0,0), "VAL": (1,0,0,0,0), "ILE": (1,0,0,0,0),
    "LEU": (1,0,0,0,0), "MET": (1,0,0,0,0), "PHE": (1,1,0,0,0),
    "TRP": (1,1,0,0,0), "PRO": (1,0,0,0,0),
    "TYR": (0,1,0,0,1), "CYS": (0,0,0,0,1), "SER": (0,0,0,0,1),
    "THR": (0,0,0,0,1), "ASN": (0,0,0,0,1), "GLN": (0,0,0,0,1),
    "HIS": (0,1,1,0,0),
    "LYS": (0,0,1,0,0), "ARG": (0,0,1,0,0),
    "ASP": (0,0,0,1,0), "GLU": (0,0,0,1,0),
    "GLY": (0,0,0,0,0),
}

def composition_stats(residues):
    """Given list of 3-letter AA codes, return property fractions.
    Returns zero-filled dict if residues is empty (e.g. NMR structures
    where fpocket pocket PDB files contain no ATOM records).
    """
    zero = {
        "frac_hydrophobic": 0.0, "frac_aromatic": 0.0,
        "frac_positive":    0.0, "frac_negative":  0.0,
        "frac_polar":       0.0, "charge_balance":  0.0,
        "n_residues":       0,
    }
    if not residues:
        return zero
    n = len(residues)
    hydro = sum(AA_PROPS.get(r, (0,)*5)[0] for r in residues) / n
    arom  = sum(AA_PROPS.get(r, (0,)*5)[1] for r in residues) / n
    pos   = sum(AA_PROPS.get(r, (0,)*5)[2] for r in residues) / n
    neg   = sum(AA_PROPS.get(r, (0,)*5)[3] for r in residues) / n
    polar = sum(AA_PROPS.get(r, (0,)*5)[4] for r in residues) / n
    return {
        "frac_hydrophobic": round(hydro, 3),
        "frac_aromatic":    round(arom, 3),
        "frac_positive":    round(pos, 3),
        "frac_negative":    round(neg, 3),
        "frac_polar":       round(polar, 3),
        "charge_balance":   round(pos - neg, 3),
        "n_residues":       n,
    }

# fpocket 4.0 places pocket atom PDB files in a pockets/ subdirectory.
# Try both locations for compatibility with older fpocket versions.
fpocket_dir = f"{OUTPUT_DIR}/fpocket/{PDB_ID.upper()}_out"

def find_pocket_pdb(fpocket_dir, pocket_num):
    """Return path to pocket PDB file, trying both fpocket 3.x and 4.x locations."""
    candidates = [
        os.path.join(fpocket_dir, "pockets", f"pocket{pocket_num}_atm.pdb"),  # fpocket 4.x
        os.path.join(fpocket_dir, f"pocket{pocket_num}_atm.pdb"),              # fpocket 3.x
        os.path.join(fpocket_dir, f"pocket{pocket_num}.pdb"),                  # older still
    ]
    for p in candidates:
        if os.path.exists(p):
            return p
    return None

pocket_compositions = []
top_pockets = pd.read_csv(f"{OUTPUT_DIR}/data/query_top_pockets.csv")

for _, row in top_pockets.iterrows():
    pocket_num = int(row["pocket_id"])
    pocket_pdb = find_pocket_pdb(fpocket_dir, pocket_num)

    residues = []
    if pocket_pdb:
        with open(pocket_pdb) as f:
            for line in f:
                if line.startswith("ATOM") or line.startswith("HETATM"):
                    res_name = line[17:20].strip()
                    if res_name in AA_PROPS:
                        residues.append(res_name)
    else:
        print(f"  Pocket {pocket_num} PDB not found in {fpocket_dir}")

    comp = composition_stats(list(set(residues)))
    comp["pocket_id"] = pocket_num
    comp["residue_list"] = ",".join(sorted(set(residues)))

    # Druggability prediction from composition
    # High hydrophobicity + aromaticity → druggable
    score = (comp.get("frac_hydrophobic", 0) * 0.5 +
             comp.get("frac_aromatic", 0) * 0.4 -
             abs(comp.get("charge_balance", 0)) * 0.3)
    comp["composition_druggability_score"] = round(score, 3)
    pocket_compositions.append(comp)

comp_df = pd.DataFrame(pocket_compositions)
comp_df.to_csv(f"{OUTPUT_DIR}/data/pocket_composition.csv", index=False)
print(comp_df[["pocket_id","frac_hydrophobic","frac_aromatic",
               "charge_balance","composition_druggability_score"]].to_string())

# Figure: composition radar for top 3 pockets
props = ["frac_hydrophobic","frac_aromatic","frac_positive","frac_negative","frac_polar"]
labels = ["Hydrophobic","Aromatic","Positive","Negative","Polar"]
angles = np.linspace(0, 2*np.pi, len(props), endpoint=False).tolist()
angles += angles[:1]

fig, ax = plt.subplots(figsize=(7, 7), subplot_kw=dict(polar=True))
colours = ["#1565C0","#E53935","#2E7D32"]
for i, (_, row) in enumerate(comp_df.head(3).iterrows()):
    vals = [row.get(p, 0) for p in props]
    vals += vals[:1]
    ax.plot(angles, vals, color=colours[i], lw=2,
            label=f"Pocket {int(row['pocket_id'])}")
    ax.fill(angles, vals, color=colours[i], alpha=0.1)
ax.set_xticks(angles[:-1])
ax.set_xticklabels(labels, fontsize=10)
ax.set_title(f"Binding Site Composition — {PDB_ID.upper()}", pad=20)
ax.legend(loc="upper right", bbox_to_anchor=(1.3, 1.1))
plt.tight_layout()
plt.savefig(f"{OUTPUT_DIR}/figures/pocket_composition.png", dpi=150,
            bbox_inches="tight")
print("Saved pocket_composition.png")
```

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

**Validation:** `druggability_output/data/pocket_composition.csv` must exist.

---

## Step 7 — B-Factor Flexibility Analysis

Parses ATOM records from the PDB to compute mean B-factors of pocket-lining
residues. High B-factor = flexible/disordered = harder to drug.

```python
# scripts/06_bfactor.py
import sys, os, warnings, json
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 params import PDB_ID, OUTPUT_DIR

pdb_path = f"{OUTPUT_DIR}/structures/{PDB_ID.upper()}.pdb"
comp_df  = pd.read_csv(f"{OUTPUT_DIR}/data/pocket_composition.csv")

# Parse all ATOM B-factors from PDB
atom_records = []
with open(pdb_path) as f:
    for line in f:
        if line.startswith("ATOM"):
            try:
                res_seq  = int(line[22:26].strip())
                res_name = line[17:20].strip()
                b_factor = float(line[60:66].strip())
                atom_records.append({"res_seq": res_seq,
                                     "res_name": res_name,
                                     "b_factor": b_factor})
            except (ValueError, IndexError):
                continue

atom_df = pd.DataFrame(atom_records)

# Detect NMR structures: B-factors are uniformly 0 in NMR PDB files
is_nmr = len(atom_df) > 0 and atom_df["b_factor"].std() < 0.01
if is_nmr:
    print("WARNING: B-factors are uniformly zero — this is likely an NMR structure.")
    print("B-factor flexibility analysis is not meaningful for NMR structures.")
    print("Flexibility will be recorded as NMR_STRUCTURE for all pockets.")

# Mean B-factor per residue
res_bfactor    = atom_df.groupby("res_seq")["b_factor"].mean()
global_mean_b  = float(atom_df["b_factor"].mean()) if len(atom_df) else 0.0
global_std_b   = float(atom_df["b_factor"].std())  if len(atom_df) else 0.0
print(f"Global mean B-factor: {global_mean_b:.2f} ± {global_std_b:.2f}"
      + (" [NMR — not meaningful]" if is_nmr else ""))

# fpocket 4.0 pocket PDB helper (same as script 05)
fpocket_dir = f"{OUTPUT_DIR}/fpocket/{PDB_ID.upper()}_out"

def find_pocket_pdb(fpocket_dir, pocket_num):
    candidates = [
        os.path.join(fpocket_dir, "pockets", f"pocket{pocket_num}_atm.pdb"),
        os.path.join(fpocket_dir, f"pocket{pocket_num}_atm.pdb"),
        os.path.join(fpocket_dir, f"pocket{pocket_num}.pdb"),
    ]
    for p in candidates:
        if os.path.exists(p):
            return p
    return None

# For each pocket, compute mean B-factor of lining residues
bfactor_results = []
for _, row in comp_df.iterrows():
    pocket_id  = int(row["pocket_id"])
    pocket_pdb = find_pocket_pdb(fpocket_dir, pocket_id)

    res_seqs = set()
    if pocket_pdb:
        with open(pocket_pdb) as f:
            for line in f:
                if line.startswith("ATOM") or line.startswith("HETATM"):
                    try:
                        res_seqs.add(int(line[22:26].strip()))
                    except ValueError:
                        pass

    if is_nmr:
        mean_b       = 0.0
        normalized_b = 0.0
        flexibility  = "NMR_STRUCTURE"
    elif res_seqs:
        pocket_b = res_bfactor[res_bfactor.index.isin(res_seqs)]
        mean_b   = float(pocket_b.mean()) if len(pocket_b) else float("nan")
        normalized_b = (mean_b - global_mean_b) / (global_std_b + 1e-9)
        flexibility  = ("RIGID"    if normalized_b < -0.5 else
                        "MODERATE" if normalized_b <  0.5 else
                        "FLEXIBLE")
    else:
        mean_b       = float("nan")
        normalized_b = 0.0
        flexibility  = "UNKNOWN"

    bfactor_results.append({
        "pocket_id":           pocket_id,
        "mean_b_factor":       round(mean_b, 2) if not np.isnan(mean_b) else None,
        "normalized_b_factor": round(normalized_b, 3),
        "flexibility_class":   flexibility,
    })
    print(f"  Pocket {pocket_id}: mean B = {mean_b:.2f} ({flexibility})"
          if not np.isnan(mean_b) else f"  Pocket {pocket_id}: ({flexibility})")

bf_df = pd.DataFrame(bfactor_results)
bf_df.to_csv(f"{OUTPUT_DIR}/data/bfactor_analysis.csv", index=False)

with open(f"{OUTPUT_DIR}/data/bfactor_summary.json", "w") as f:
    json.dump({
        "global_mean_b": round(global_mean_b, 2),
        "global_std_b":  round(global_std_b, 2),
        "is_nmr":        is_nmr,
        "pockets":       bfactor_results,
    }, f, indent=2)

# Figure
fig, ax = plt.subplots(figsize=(7, 4))
colours_map = {"RIGID": "#2E7D32", "MODERATE": "#F9A825", "FLEXIBLE": "#E53935"}
colours = [colours_map[r["flexibility_class"]] for r in bfactor_results]
ax.bar([f"Pocket {r['pocket_id']}" for r in bfactor_results],
       [r["mean_b_factor"] for r in bfactor_results],
       color=colours, alpha=0.85)
ax.axhline(global_mean_b, color="black", ls="--", lw=1.5,
           label=f"Global mean ({global_mean_b:.1f})")
ax.set_ylabel("Mean B-factor (Ų)")
ax.set_title(f"Pocket Flexibility — {PDB_ID.upper()}")
ax.legend()
patches = [plt.Rectangle((0,0),1,1, color=v, alpha=0.85, label=k)
           for k, v in colours_map.items()]
ax.legend(handles=patches + [plt.Line2D([0],[0], color="black",
          ls="--", label=f"Global mean ({global_mean_b:.1f})")], fontsize=9)
plt.tight_layout()
plt.savefig(f"{OUTPUT_DIR}/figures/bfactor.png", dpi=150)
print("Saved bfactor.png")
```

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

**Validation:** `druggability_output/data/bfactor_analysis.csv` must exist.

---

## Step 8 — Multi-Structure Pocket Stability

Runs pocket analysis across all homologous structures to assess whether the
top-ranked pocket is consistently present (stable site) or variable (induced-fit
or artifact).

```python
# scripts/07_pocket_stability.py
import sys, os, warnings, json, pickle
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 params import PDB_ID, OUTPUT_DIR, DRUGGABILITY_SCORE_MIN

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

with open(f"{OUTPUT_DIR}/data/homologs.pkl", "rb") as f:
    homologs = pickle.load(f)

score_col = next((c for c in df.columns if "drug" in c.lower()), None)
vol_col   = next((c for c in df.columns if "volume" in c.lower() or "vol" in c.lower()), None)

# Structures to compare: query + homologs
structs_to_compare = [PDB_ID.upper()] + [h.upper() for h in homologs]
compare_df = df[df["pdb_id"].str.upper().isin(structs_to_compare)]

stability = []
for pid in structs_to_compare:
    sub = compare_df[compare_df["pdb_id"].str.upper() == pid]
    if len(sub) == 0:
        continue
    best_score = float(sub[score_col].max())
    n_druggable = int((sub[score_col] >= DRUGGABILITY_SCORE_MIN).sum())
    stability.append({
        "pdb_id":           pid,
        "is_query":         pid == PDB_ID.upper(),
        "n_pockets":        len(sub),
        "best_drug_score":  round(best_score, 4),
        "n_druggable_pockets": n_druggable,
        "has_druggable_pocket": best_score >= DRUGGABILITY_SCORE_MIN,
    })

stab_df = pd.DataFrame(stability)
stab_df.to_csv(f"{OUTPUT_DIR}/data/pocket_stability.csv", index=False)

n_total    = len(stab_df)
n_drugg    = int(stab_df["has_druggable_pocket"].sum())
stability_pct = round(n_drugg / n_total * 100, 1) if n_total > 0 else 0

summary = {
    "structures_analysed": n_total,
    "structures_with_druggable_pocket": n_drugg,
    "pocket_stability_pct": stability_pct,
    "stability_verdict": ("STABLE"   if stability_pct >= 80 else
                          "VARIABLE" if stability_pct >= 50 else
                          "UNSTABLE"),
}
print(json.dumps(summary, indent=2))
with open(f"{OUTPUT_DIR}/data/stability_summary.json", "w") as f:
    json.dump(summary, f, indent=2)

# Figure
fig, ax = plt.subplots(figsize=(max(6, n_total * 1.2), 4))
colours = ["#1565C0" if r["is_query"] else
           "#2E7D32" if r["has_druggable_pocket"] else "#E53935"
           for _, r in stab_df.iterrows()]
ax.bar(stab_df["pdb_id"], stab_df["best_drug_score"], color=colours, alpha=0.85)
ax.axhline(DRUGGABILITY_SCORE_MIN, color="grey", ls="--", lw=1.5,
           label=f"Threshold ({DRUGGABILITY_SCORE_MIN})")
ax.set_ylabel("Best Pocket Drug Score")
ax.set_title(f"Pocket Stability Across Homologous Structures ({stability_pct}% druggable)")
ax.tick_params(axis="x", rotation=30)
ax.legend()
plt.tight_layout()
plt.savefig(f"{OUTPUT_DIR}/figures/pocket_stability.png", dpi=150)
print("Saved pocket_stability.png")
```

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

**Validation:** `druggability_output/data/stability_summary.json` must contain `stability_verdict`.

---

## Step 9 — Composite Druggability Score

Combines all five evidence streams into a weighted composite score
and final druggability classification.

```python
# scripts/08_composite_score.py
import sys, os, json
sys.path.insert(0, os.path.dirname(os.path.dirname(os.path.abspath(__file__))))
import pandas as pd
import numpy as np
from params import PDB_ID, OUTPUT_DIR, DRUGGABILITY_SCORE_MIN

bench  = json.load(open(f"{OUTPUT_DIR}/data/benchmark_result.json"))
chembl = json.load(open(f"{OUTPUT_DIR}/data/chembl_evidence.json"))
stab   = json.load(open(f"{OUTPUT_DIR}/data/stability_summary.json"))
bf     = json.load(open(f"{OUTPUT_DIR}/data/bfactor_summary.json"))
comp   = pd.read_csv(f"{OUTPUT_DIR}/data/pocket_composition.csv")

# --- Component scores (all normalised 0-1) ---

# 1. fpocket geometry score
geom_raw = bench["query_best_drug_score"]
geom_score = min(geom_raw / 1.0, 1.0)   # fpocket drug_score range ~ 0-1

# 2. Benchmarking percentile
pct_score  = bench["druggability_percentile_vs_refs"] / 100.0

# 3. Bioactivity evidence
bact_map   = {"WELL_CHARACTERISED": 1.0, "SOME_EVIDENCE": 0.5,
              "SPARSE_EVIDENCE": 0.2, "NO_EVIDENCE": 0.0}
bact_score = bact_map.get(chembl.get("bioactivity_verdict", "NO_EVIDENCE"), 0.0)

# 4. Pocket stability
stab_score = stab["pocket_stability_pct"] / 100.0

# 5. Binding site composition (top pocket)
comp_score_raw = (comp["composition_druggability_score"].iloc[0]
                  if len(comp) else 0.0)
comp_score = max(0.0, min(comp_score_raw, 1.0))

# 6. Flexibility (inverse — rigid = good; NMR structures get neutral 0.5)
pockets_bf = bf.get("pockets", [])
is_nmr     = bf.get("is_nmr", False)
if is_nmr or not pockets_bf:
    flex_score = 0.5   # neutral — B-factors not meaningful for NMR
else:
    top_b_norm = pockets_bf[0].get("normalized_b_factor", 0.0) or 0.0
    flex_score = max(0.0, min(1.0 - (top_b_norm + 1) / 2, 1.0))

# Weighted composite
WEIGHTS = {
    "fpocket_geometry":    0.25,
    "benchmark_percentile": 0.20,
    "bioactivity_evidence": 0.25,
    "pocket_stability":    0.15,
    "pocket_composition":  0.10,
    "flexibility":         0.05,
}
SCORES = {
    "fpocket_geometry":    geom_score,
    "benchmark_percentile": pct_score,
    "bioactivity_evidence": bact_score,
    "pocket_stability":    stab_score,
    "pocket_composition":  comp_score,
    "flexibility":         flex_score,
}

composite = sum(WEIGHTS[k] * SCORES[k] for k in WEIGHTS)
final_verdict = ("HIGH CONFIDENCE DRUGGABLE"  if composite >= 0.70 else
                 "LIKELY DRUGGABLE"            if composite >= 0.50 else
                 "BORDERLINE"                  if composite >= 0.35 else
                 "CHALLENGING")

result = {
    "query_pdb":       PDB_ID.upper(),
    "composite_score": round(composite, 4),
    "final_verdict":   final_verdict,
    "component_scores": {k: round(SCORES[k], 4) for k in SCORES},
    "weights":          WEIGHTS,
}
print(json.dumps(result, indent=2))
with open(f"{OUTPUT_DIR}/data/composite_score.json", "w") as f:
    json.dump(result, f, indent=2)
```

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

**Validation:** `druggability_output/data/composite_score.json` must contain `final_verdict`.

---

## Step 10 — Compile Report

```python
# scripts/09_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 PDB_ID, OUTPUT_DIR

def img_b64(path):
    p = pathlib.Path(path)
    if not p.exists():
        return ""
    return base64.b64encode(p.read_bytes()).decode()

bench   = json.load(open(f"{OUTPUT_DIR}/data/benchmark_result.json"))
chembl  = json.load(open(f"{OUTPUT_DIR}/data/chembl_evidence.json"))
stab    = json.load(open(f"{OUTPUT_DIR}/data/stability_summary.json"))
bf      = json.load(open(f"{OUTPUT_DIR}/data/bfactor_summary.json"))
comp    = json.load(open(f"{OUTPUT_DIR}/data/composite_score.json"))
comp_df = pd.read_csv(f"{OUTPUT_DIR}/data/pocket_composition.csv")
topk_df = pd.read_csv(f"{OUTPUT_DIR}/data/query_top_pockets.csv")

score_col = next((c for c in topk_df.columns if "drug" in c.lower()), topk_df.columns[-1])

VERDICT_COLOUR = {
    "HIGH CONFIDENCE DRUGGABLE": "#1B5E20",
    "LIKELY DRUGGABLE":          "#2E7D32",
    "BORDERLINE":                "#F57F17",
    "CHALLENGING":               "#B71C1C",
}
v_colour = VERDICT_COLOUR.get(comp["final_verdict"], "#333")

TEMPLATE = """<!DOCTYPE html><html><head><meta charset="utf-8">
<title>Druggability Dossier — {{ pdb_id }}</title>
<style>
  body{font-family:Georgia,serif;max-width:980px;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;font-size:0.93em}
  th{background:#e8eaf6}
  img{max-width:100%;border:1px solid #eee;border-radius:4px;margin:8px 0}
  .grid{display:grid;grid-template-columns:repeat(3,1fr);gap:14px;margin:20px 0}
  .card{background:#f5f5f5;border-radius:6px;padding:14px;text-align:center}
  .stat{font-size:1.8em;font-weight:bold;color:#1565c0}
  .label{font-size:0.82em;color:#555}
  .verdict{padding:16px;border-radius:8px;color:white;text-align:center;
           font-size:1.5em;font-weight:bold;margin:20px 0;
           background:{{ v_colour }}}
  .bar-wrap{background:#eee;border-radius:4px;height:20px;margin:4px 0}
  .bar{height:20px;border-radius:4px;background:#1565C0;min-width:4px}
</style></head><body>
<h1>Druggability Dossier: {{ pdb_id }}</h1>

<div class="verdict">{{ comp.final_verdict }} — Composite Score: {{ "%.3f"|format(comp.composite_score) }}</div>

<div class="grid">
  <div class="card"><div class="stat">{{ "%.3f"|format(bench.query_best_drug_score) }}</div>
    <div class="label">fpocket Drug Score</div></div>
  <div class="card"><div class="stat">{{ bench.druggability_percentile_vs_refs }}%</div>
    <div class="label">Percentile vs. known druggable targets</div></div>
  <div class="card"><div class="stat">{{ chembl.bioactivity_verdict.replace("_"," ") }}</div>
    <div class="label">ChEMBL Bioactivity Evidence</div></div>
</div>

<h2>1. Pocket Geometry Ranking</h2>
<img src="data:image/png;base64,{{ pocket_ranking_b64 }}" alt="Pocket ranking">
<table>
  <tr><th>Pocket</th><th>Drug Score</th></tr>
  {% for _, row in topk.iterrows() %}
  <tr><td>Pocket {{ row.pocket_id | int }}</td>
      <td>{{ "%.4f"|format(row[score_col]) }}</td></tr>
  {% endfor %}
</table>

<h2>2. Druggability Benchmarking</h2>
<img src="data:image/png;base64,{{ benchmark_b64 }}" alt="Benchmarking">
<p>The query protein's best pocket scores at the
<strong>{{ bench.druggability_percentile_vs_refs }}th percentile</strong>
of known druggable reference targets
(fpocket drug score: {{ "%.3f"|format(bench.query_best_drug_score) }}).
Reference set KS separation statistic: {{ bench.ref_separation_ks }}
(p = {{ bench.ref_separation_p }}).</p>

<h2>3. ChEMBL Bioactivity Evidence</h2>
{% if chembl.total_bioactivity_records > 0 %}
<table>
  <tr><th>Metric</th><th>Value</th></tr>
  <tr><td>UniProt IDs</td><td>{{ chembl.uniprot_ids | join(", ") }}</td></tr>
  <tr><td>ChEMBL targets</td><td>{{ chembl.chembl_targets | join(", ") }}</td></tr>
  <tr><td>Total bioactivity records</td><td>{{ chembl.total_bioactivity_records }}</td></tr>
  <tr><td>Unique compounds</td><td>{{ chembl.n_unique_compounds }}</td></tr>
  <tr><td>Best potency</td><td>{{ chembl.best_potency_nM }} nM</td></tr>
  <tr><td>Median potency</td><td>{{ chembl.median_potency_nM }} nM</td></tr>
  <tr><td>Compounds ≤ 100 nM</td><td>{{ chembl.n_sub_100nM }}</td></tr>
  <tr><td>Compounds ≤ 1 µM</td><td>{{ chembl.n_sub_1uM }}</td></tr>
</table>
{% else %}
<p>No bioactivity records found in ChEMBL for this target.</p>
{% endif %}

<h2>4. Binding Site Composition</h2>
<img src="data:image/png;base64,{{ composition_b64 }}" alt="Pocket composition">
<table>
  <tr><th>Pocket</th><th>Hydrophobic</th><th>Aromatic</th><th>Positive</th>
      <th>Negative</th><th>Polar</th><th>Composition Score</th></tr>
  {% for _, row in comp_df.iterrows() %}
  <tr>
    <td>Pocket {{ row.pocket_id | int }}</td>
    <td>{{ "%.2f"|format(row.frac_hydrophobic) }}</td>
    <td>{{ "%.2f"|format(row.frac_aromatic) }}</td>
    <td>{{ "%.2f"|format(row.frac_positive) }}</td>
    <td>{{ "%.2f"|format(row.frac_negative) }}</td>
    <td>{{ "%.2f"|format(row.frac_polar) }}</td>
    <td>{{ "%.3f"|format(row.composition_druggability_score) }}</td>
  </tr>
  {% endfor %}
</table>

<h2>5. Pocket Flexibility (B-factors)</h2>
<img src="data:image/png;base64,{{ bfactor_b64 }}" alt="B-factors">
<p>Global structure mean B-factor: {{ bf.global_mean_b }} ± {{ bf.global_std_b }} Ų</p>
<table>
  <tr><th>Pocket</th><th>Mean B-factor (Ų)</th><th>Normalised</th><th>Classification</th></tr>
  {% for p in bf.pockets %}
  <tr>
    <td>Pocket {{ p.pocket_id }}</td>
    <td>{{ p.mean_b_factor }}</td>
    <td>{{ p.normalized_b_factor }}</td>
    <td>{{ p.flexibility_class }}</td>
  </tr>
  {% endfor %}
</table>

<h2>6. Pocket Stability Across Homologous Structures</h2>
<img src="data:image/png;base64,{{ stability_b64 }}" alt="Stability">
<p>A druggable pocket was detected in
<strong>{{ stab.structures_with_druggable_pocket }} of {{ stab.structures_analysed }}</strong>
homologous structures ({{ stab.pocket_stability_pct }}% — <strong>{{ stab.stability_verdict }}</strong>).</p>

<h2>7. Composite Druggability Score</h2>
<table>
  <tr><th>Evidence Stream</th><th>Weight</th><th>Score</th><th>Contribution</th></tr>
  {% for key, score in comp.component_scores.items() %}
  <tr>
    <td>{{ key.replace("_"," ").title() }}</td>
    <td>{{ "%.0f"|format(comp.weights[key]*100) }}%</td>
    <td>{{ "%.3f"|format(score) }}</td>
    <td>
      <div class="bar-wrap">
        <div class="bar" style="width:{{ (score*100)|int }}%"></div>
      </div>
    </td>
  </tr>
  {% endfor %}
  <tr style="font-weight:bold">
    <td>COMPOSITE</td><td>100%</td>
    <td>{{ "%.3f"|format(comp.composite_score) }}</td>
    <td><div class="bar-wrap"><div class="bar"
        style="width:{{ (comp.composite_score*100)|int }}%;background:{{ v_colour }}">
        </div></div></td>
  </tr>
</table>
</body></html>
"""

html = Template(TEMPLATE).render(
    pdb_id          = PDB_ID.upper(),
    bench           = bench,
    chembl          = chembl,
    stab            = stab,
    bf              = bf,
    comp            = comp,
    comp_df         = comp_df,
    topk            = topk_df,
    score_col       = score_col,
    v_colour        = v_colour,
    pocket_ranking_b64 = img_b64(f"{OUTPUT_DIR}/figures/pocket_ranking.png"),
    benchmark_b64      = img_b64(f"{OUTPUT_DIR}/figures/benchmark.png"),
    composition_b64    = img_b64(f"{OUTPUT_DIR}/figures/pocket_composition.png"),
    bfactor_b64        = img_b64(f"{OUTPUT_DIR}/figures/bfactor.png"),
    stability_b64      = img_b64(f"{OUTPUT_DIR}/figures/pocket_stability.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/09_report.py
```

**Expected final output:** `druggability_output/report.html` — open in any browser.

---

## Expected Outputs

```
druggability_output/
├── structures/
│   ├── 2ITY.pdb            # query structure
│   ├── 3ERT.pdb            # druggable refs
│   └── ...
├── fpocket/
│   ├── 2ITY_out/           # fpocket results per structure
│   └── ...
├── data/
│   ├── all_pockets.csv             # all pocket descriptors
│   ├── query_top_pockets.csv       # top N pockets for query
│   ├── benchmark_result.json       # percentile vs. reference set
│   ├── chembl_evidence.json        # bioactivity from ChEMBL
│   ├── pocket_composition.csv      # AA composition per pocket
│   ├── bfactor_analysis.csv        # flexibility per pocket
│   ├── bfactor_summary.json
│   ├── pocket_stability.csv        # scores across homologs
│   ├── stability_summary.json
│   ├── composite_score.json        # final weighted verdict
│   └── homologs.pkl
├── figures/
│   ├── benchmark.png
│   ├── pocket_ranking.png
│   ├── pocket_composition.png
│   ├── bfactor.png
│   └── pocket_stability.png
└── report.html                     ← primary deliverable
```

---

## Step 11 — Batch Reports for All Structures

Generates a druggability dossier for every structure that has an fpocket
output directory, not just the query. Produces a cross-target summary table.

```python
# scripts/10_batch_reports.py
import sys, os, json, importlib, shutil
sys.path.insert(0, os.path.dirname(os.path.dirname(os.path.abspath(__file__))))
import params as P

FPOCKET_DIR  = f"{P.OUTPUT_DIR}/fpocket"
REPORTS_DIR  = f"{P.OUTPUT_DIR}/reports"
os.makedirs(REPORTS_DIR, exist_ok=True)

# Find every PDB ID that has an fpocket output directory
all_pdb_ids = sorted([
    d.replace("_out", "")
    for d in os.listdir(FPOCKET_DIR)
    if d.endswith("_out") and os.path.isdir(os.path.join(FPOCKET_DIR, d))
])
print(f"Found {len(all_pdb_ids)} structures with fpocket output: {all_pdb_ids}")

summary_rows = []
failed = []

for pdb_id in all_pdb_ids:
    print(f"\n{'='*60}")
    print(f"Processing {pdb_id} ...")
    out_dir = os.path.join(REPORTS_DIR, pdb_id)
    os.makedirs(f"{out_dir}/data", exist_ok=True)
    os.makedirs(f"{out_dir}/figures", exist_ok=True)

    # Temporarily override params for this PDB ID
    P.PDB_ID     = pdb_id
    P.OUTPUT_DIR = out_dir

    # Symlink structures/ and fpocket/ so scripts can find them
    for sub in ["structures", "fpocket"]:
        link = os.path.join(out_dir, sub)
        target = os.path.abspath(f"{P.OUTPUT_DIR.replace(out_dir, '').lstrip('/')}"
                                  f"/../{sub}").replace("//", "/")
        src = os.path.abspath(f"druggability_output/{sub}")
        if not os.path.exists(link):
            os.symlink(src, link)

    # Run each analysis script with overridden params
    import importlib.util

    def run_script(script_name):
        path = os.path.join("scripts", script_name)
        if not os.path.exists(path):
            print(f"  SKIP: {path} not found")
            return False
        spec = importlib.util.spec_from_file_location("script", path)
        mod  = importlib.util.module_from_spec(spec)
        try:
            spec.loader.exec_module(mod)
            return True
        except Exception as e:
            print(f"  ERROR in {script_name}: {e}")
            return False

    steps = [
        ("03_benchmark.py",         "benchmark_result.json"),
        ("04_chembl_evidence.py",   "chembl_evidence.json"),
        ("05_pocket_composition.py","pocket_composition.csv"),
        ("06_bfactor.py",           "bfactor_analysis.csv"),
        ("07_pocket_stability.py",  "stability_summary.json"),
        ("08_composite_score.py",   "composite_score.json"),
        ("09_report.py",            "report.html"),
    ]

    ok = True
    for script, output_file in steps:
        out_file_path = (os.path.join(out_dir, "data", output_file)
                         if not output_file.endswith(".html")
                         else os.path.join(out_dir, output_file))
        if os.path.exists(out_file_path):
            print(f"  {script}: already done, skipping")
            continue
        success = run_script(script)
        if not success:
            ok = False
            break

    # Read composite score for summary
    comp_path = os.path.join(out_dir, "data", "composite_score.json")
    if os.path.exists(comp_path):
        comp = json.load(open(comp_path))
        summary_rows.append({
            "pdb_id":          pdb_id,
            "composite_score": comp["composite_score"],
            "verdict":         comp["final_verdict"],
        })
    else:
        failed.append(pdb_id)

# Print summary table
print(f"\n{'='*70}")
print(f"{'PDB ID':<10} {'Score':>7}  {'Verdict'}")
print(f"{'-'*70}")
for r in sorted(summary_rows, key=lambda x: -x["composite_score"]):
    print(f"{r['pdb_id']:<10} {r['composite_score']:>7.4f}  {r['verdict']}")
if failed:
    print(f"\nFailed: {failed}")

import pandas as pd
pd.DataFrame(summary_rows).sort_values(
    "composite_score", ascending=False
).to_csv(f"{P.OUTPUT_DIR}/batch_summary.csv", index=False)
print(f"\nSummary saved to {REPORTS_DIR}/../batch_summary.csv")
```

```bash
python scripts/10_batch_reports.py
```

**Expected output:** A ranked summary table of all structures with composite
scores and verdicts, plus individual `druggability_output/reports/{PDB_ID}/report.html`
for each target.

**Validation:**
```bash
ls druggability_output/reports/
cat druggability_output/batch_summary.csv
```

---

## Expected Outputs

```
druggability_output/
├── structures/
│   ├── 2ITY.pdb            # query + all reference structures
│   └── ...
├── fpocket/
│   ├── 2ITY_out/           # fpocket results per structure
│   └── ...
├── data/                   # analysis outputs for query PDB
├── figures/                # figures for query PDB
├── report.html             # query PDB report
├── batch_summary.csv       # all-targets ranked table
└── reports/
    ├── 2ITY/report.html    # per-target dossiers
    ├── 1IEP/report.html
    └── ...
```

---

## Reproducing on a Different Target

1. Edit `params.py`: change `PDB_ID` to any 4-letter RCSB identifier
2. Re-run steps 2–11 in order
3. To run on all reference structures at once, run only Step 11 after Step 3