От in silico до wet lab: оценка на ефективността на дизайна на протеини с ИИ

import subprocess, sys, warnings, itertools, math
warnings.filterwarnings("ignore")
import importlib.util
_needed = {"huggingface_hub": "huggingface_hub>=0.24", "pyarrow": "pyarrow",
          "pandas": "pandas", "sklearn": "scikit-learn",
          "matplotlib": "matplotlib", "scipy": "scipy"}
_missing = [pkg for mod, pkg in _needed.items() if importlib.util.find_spec(mod) is None]
if _missing:
   print("installing:", ", ".join(_missing))
   subprocess.run([sys.executable, "-m", "pip", "install", "-q", *_missing], check=False)
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy import stats
from huggingface_hub import HfApi, hf_hub_download
from sklearn.metrics import roc_auc_score, cohen_kappa_score, average_precision_score
from sklearn.model_selection import GroupKFold, StratifiedKFold
from sklearn.ensemble import HistGradientBoostingClassifier
from sklearn.inspection import permutation_importance
SEED = 0
rng_global = np.random.default_rng(SEED)
pd.set_option("display.width", 200)
pd.set_option("display.max_columns", 100)
plt.rcParams.update({"figure.dpi": 110, "font.size": 9, "axes.grid": True,
                    "grid.alpha": 0.25, "axes.spines.top": False, "axes.spines.right": False})
REPO = "Anthropic/claude-protein-binder-design"
BAR = "=" * 78
def head(n, title):
   prefix = f"{n}. " if str(n) else ""
   print(f"\n{BAR}\n  {prefix}{title}\n{BAR}")
head(1, "TABLE DISCOVERY")
api = HfApi()
repo_files = api.list_repo_files(REPO, repo_type="dataset")
TABLES = {}
for f in repo_files:
   if f.startswith("data/tables/") and f.endswith(".parquet"):
       key = f[len("data/tables/"): -len(".parquet")].replace("/", "_")
       TABLES[key] = f
print(f"Found {len(TABLES)} Parquet tables:")
for k in sorted(TABLES):
   print(f"   - {k:38s} {TABLES[k]}")
def load_table(name: str) -> pd.DataFrame:
   """Load a subset by its viewer name, with a datasets-library fallback."""
   if name in TABLES:
       return pd.read_parquet(hf_hub_download(REPO, TABLES[name], repo_type="dataset"))
   from datasets import load_dataset
   return load_dataset(REPO, name, split="full").to_pandas()
ds = load_table("design_summary")
print(f"\ndesign_summary: {ds.shape[0]:,} rows x {ds.shape[1]} columns")

head(2, "SCHEMA + EVALUABLE SET")
CALLS = {"binder", "non_binder"}
tested = ds["adaptyv_binding"].isin(CALLS) | ds["twist_binding"].isin(CALLS)
ev = ds[tested].copy()
ev["y"] = ev["binder_final"].astype(int)
print(f"All designs               : {len(ds):,}")
print(f"Evaluable (>=1 vendor call): {len(ev):,}")
print(f"Confirmed binders          : {int(ev['y'].sum()):,}  "
     f"({100 * ev['y'].mean():.1f}% base rate)")
print(f"Never measured             : {len(ds) - len(ev):,}")
print("\nCategorical levels:")
for c in ["design_model", "campaign", "generator", "sequence_design_method", "vendor_agreement"]:
   vals = ds[c].astype(str).value_counts()
   print(f"  {c:24s} ({len(vals)}): {', '.join(vals.index[:6])}"
         + (" ..." if len(vals) > 6 else ""))
print(f"\nTargets ({ds['target'].nunique()}): {', '.join(sorted(ds['target'].unique()))}")
print(f"Binder length: {ds.binder_length.min()}-{ds.binder_length.max()} aa "
     f"(median {ds.binder_length.median():.0f})")
head(3, "HIT-RATE LANDSCAPE")
def wilson(k, n, z=1.96):
   if n == 0:
       return (np.nan, np.nan, np.nan)
   p = k / n
   d = 1 + z**2 / n
   c = (p + z**2 / (2 * n)) / d
   h = z * math.sqrt(p * (1 - p) / n + z**2 / (4 * n**2)) / d
   return p, max(0.0, c - h), min(1.0, c + h)
def rate_table(df, by):
   rows = []
   for key, g in df.groupby(by, dropna=False):
       p, lo, hi = wilson(int(g.y.sum()), len(g))
       rows.append({by: key, "n": len(g), "hits": int(g.y.sum()),
                    "rate": p, "lo": lo, "hi": hi})
   return pd.DataFrame(rows).sort_values("rate", ascending=False).reset_index(drop=True)
for dim in ["design_model", "campaign", "generator", "sequence_design_method"]:
   t = rate_table(ev, dim)
   print(f"\n--- hit rate by {dim} ---")
   print(t.to_string(index=False,
                     formatters={"rate": "{:.3f}".format, "lo": "{:.3f}".format, "hi": "{:.3f}".format}))
tt = rate_table(ev, "target")
fig, ax = plt.subplots(figsize=(9, 4.2))
ax.bar(tt.target, tt.rate, color="#4C72B0")
ax.errorbar(tt.target, tt.rate,
           yerr=[(tt.rate - tt.lo).clip(lower=0), (tt.hi - tt.rate).clip(lower=0)],
           fmt="none", ecolor="0.25", capsize=3, lw=1)
ax.axhline(ev.y.mean(), ls="--", c="crimson", lw=1, label=f"pooled {ev.y.mean():.2f}")
ax.set_ylabel("experimental hit rate"); ax.set_title("Hit rate by target (Wilson 95% CI)")
ax.tick_params(axis="x", rotation=55); ax.legend(); plt.tight_layout(); plt.show()
print("\nRead this plot as the dominant effect size in the dataset: target choice "
     "swamps generator choice. Any model comparison that does not stratify by "
     "target is mostly measuring which targets that model was pointed at.")

head(4, "PER-PREDICTOR DISCRIMINATIVE POWER")
PREDICTORS = sorted({c[len("ipsae_min_"):] for c in ds.columns if c.startswith("ipsae_min_")})
print(f"Predictors ({len(PREDICTORS)}): {', '.join(PREDICTORS)}")
def auc_ci(y, s, n_boot=300, seed=SEED):
   s = np.asarray(s, dtype=float); y = np.asarray(y, dtype=int)
   m = ~np.isnan(s)
   y, s = y[m], s[m]
   if len(y) < 30 or len(np.unique(y)) < 2:
       return dict(auc=np.nan, lo=np.nan, hi=np.nan, n=len(y), ap=np.nan)
   base = roc_auc_score(y, s)
   ap = average_precision_score(y, s)
   rng = np.random.default_rng(seed)
   idx, boots = np.arange(len(y)), []
   for _ in range(n_boot):
       b = rng.choice(idx, len(idx), replace=True)
       if len(np.unique(y[b])) > 1:
           boots.append(roc_auc_score(y[b], s[b]))
   lo, hi = (np.percentile(boots, [2.5, 97.5]) if boots else (np.nan, np.nan))
   return dict(auc=base, lo=lo, hi=hi, n=len(y), ap=ap)
rows = []
for p in PREDICTORS:
   for metric in ["ipsae_min", "sc_dockq"]:
       col = f"{metric}_{p}"
       if col in ev.columns:
           r = auc_ci(ev.y, ev[col])
           rows.append({"predictor": p, "metric": metric, **r})
perf = pd.DataFrame(rows)
piv = perf.pivot(index="predictor", columns="metric", values="auc").sort_values("ipsae_min", ascending=False)
print("\nAUC vs experimental binder_final:")
print(perf.sort_values("auc", ascending=False).to_string(
   index=False, formatters={c: "{:.3f}".format for c in ["auc", "lo", "hi", "ap"]}))
fig, ax = plt.subplots(figsize=(9, 4.2))
x = np.arange(len(piv)); w = 0.38
for i, (metric, colr) in enumerate([("ipsae_min", "#4C72B0"), ("sc_dockq", "#DD8452")]):
   sub = perf[perf.metric == metric].set_index("predictor").reindex(piv.index)
   lo_err = (sub.auc - sub.lo).clip(lower=0).fillna(0)
   hi_err = (sub.hi - sub.auc).clip(lower=0).fillna(0)
   ax.bar(x + (i - 0.5) * w, sub.auc, w, label=metric, color=colr)
   ax.errorbar(x + (i - 0.5) * w, sub.auc,
               yerr=[lo_err, hi_err], fmt="none", ecolor="0.3", capsize=2, lw=0.9)
ax.axhline(0.5, ls="--", c="crimson", lw=1)
ax.set_xticks(x); ax.set_xticklabels(piv.index, rotation=45, ha="right")
ax.set_ylabel("AUC"); ax.set_ylim(0.35, None)
ax.set_title("In-silico score vs wet-lab binding, by structure predictor")
ax.legend(); plt.tight_layout(); plt.show()
print("Interpretation: AUCs land well above chance but far below the ~0.9 you "
     "would need to trust a single filter. That gap is the entire practical "
     "reason this dataset exists.")
head(5, "CONSENSUS SCORING")
ips_cols = [f"ipsae_min_{p}" for p in PREDICTORS if f"ipsae_min_{p}" in ev.columns]
dq_cols = [f"sc_dockq_{p}" for p in PREDICTORS if f"sc_dockq_{p}" in ev.columns]
def pct_rank(df, cols):
   return df[cols].rank(pct=True, na_option="keep")
R_ips, R_dq = pct_rank(ev, ips_cols), pct_rank(ev, dq_cols)
ev["cons_ipsae"] = R_ips.mean(axis=1)
ev["cons_dockq"] = R_dq.mean(axis=1)
ev["cons_all"] = pd.concat([R_ips, R_dq], axis=1).mean(axis=1)
ev["cons_median"] = pd.concat([R_ips, R_dq], axis=1).median(axis=1)
ev["cons_min"] = pd.concat([R_ips, R_dq], axis=1).min(axis=1)
ev["cons_disagree"] = pd.concat([R_ips, R_dq], axis=1).std(axis=1)
best_single = perf.loc[perf.auc.idxmax()]
print(f"Best single column: {best_single.metric}_{best_single.predictor}  AUC={best_single.auc:.3f}")
print()
for name in ["cons_ipsae", "cons_dockq", "cons_all", "cons_median", "cons_min", "cons_disagree"]:
   r = auc_ci(ev.y, ev[name])
   print(f"  {name:16s} AUC={r['auc']:.3f}  [{r['lo']:.3f}, {r['hi']:.3f}]  AP={r['ap']:.3f}")
corr = ev[ips_cols].corr(method="spearman")
fig, ax = plt.subplots(figsize=(6.2, 5.2))
im = ax.imshow(corr.values, cmap="viridis", vmin=0, vmax=1)
lbl = [c.replace("ipsae_min_", "") for c in ips_cols]
ax.set_xticks(range(len(lbl))); ax.set_xticklabels(lbl, rotation=90)
ax.set_yticks(range(len(lbl))); ax.set_yticklabels(lbl)
ax.set_title("Spearman correlation between predictors (ipSAE)")
ax.grid(False); fig.colorbar(im, shrink=0.8); plt.tight_layout(); plt.show()
print("\nIf every off-diagonal cell were ~1.0 there would be no ensemble gain to "
     "harvest. The moderate correlations are why cons_all typically edges out "
     "the best single predictor — and why disagreement itself carries signal.")

head(6, "BUDGET CURVES (precision@N)")
def budget_curve(df, score_col, max_n=400):
   d = df[[score_col, "y"]].dropna().sort_values(score_col, ascending=False)
   hits = d.y.values.cumsum()
   n = np.arange(1, len(d) + 1)
   k = min(max_n, len(d))
   return n[:k], (hits / n)[:k]
fig, ax = plt.subplots(figsize=(8, 4.4))
best_col = f"{best_single.metric}_{best_single.predictor}"
for col, lab, style in [(best_col, f"best single ({best_col})", "-"),
                       ("cons_all", "consensus (rank-avg, all)", "-"),
                       ("cons_min", "consensus (unanimity/min)", "--")]:
   n, prec = budget_curve(ev, col)
   ax.plot(n, prec, style, lw=1.8, label=lab)
ax.axhline(ev.y.mean(), ls=":", c="crimson", lw=1.4, label=f"random baseline ({ev.y.mean():.2f})")
ax.set_xlabel("designs ordered for wet-lab testing (N, best-first)")
ax.set_ylabel("hit rate among top N"); ax.set_title("How much does in-silico triage buy you?")
ax.legend(); plt.tight_layout(); plt.show()
print("Enrichment at small budgets:")
for N in [25, 50, 100, 200]:
   line = f"  N={N:4d} | random {ev.y.mean():.3f}"
   for col, lab in [(best_col, "best-single"), ("cons_all", "consensus")]:
       n, prec = budget_curve(ev, col, max_n=N)
       line += f" | {lab} {prec[-1]:.3f} ({prec[-1] / ev.y.mean():.2f}x)"
   print(line)
head(7, "VENDOR CONCORDANCE")
both = ev[ev.adaptyv_binding.isin(CALLS) & ev.twist_binding.isin(CALLS)]
ct = pd.crosstab(both.adaptyv_binding, both.twist_binding)
print(f"Designs with calls from BOTH vendors: {len(both):,}\n")
print(ct.to_string())
if len(both) > 10:
   kappa = cohen_kappa_score(both.adaptyv_binding, both.twist_binding)
   agree = (both.adaptyv_binding == both.twist_binding).mean()
   print(f"\nRaw agreement: {agree:.3f}   Cohen's kappa: {kappa:.3f}")
   print("Kappa well under 1.0 means part of the 'unpredictable' variance above "
         "is assay disagreement, not model failure.")
kd = ev[["adaptyv_kd_nM", "twist_kd_nM"]].dropna()
kd = kd[(kd > 0).all(axis=1)]
if len(kd) > 10:
   rho, pv = stats.spearmanr(kd.adaptyv_kd_nM, kd.twist_kd_nM)
   fig, ax = plt.subplots(figsize=(4.8, 4.6))
   ax.scatter(kd.adaptyv_kd_nM, kd.twist_kd_nM, s=16, alpha=0.6, c="#4C72B0", edgecolor="none")
   lims = [min(kd.min()) * 0.5, max(kd.max()) * 2]
   ax.plot(lims, lims, "k--", lw=1)
   ax.set_xscale("log"); ax.set_yscale("log")
   ax.set_xlabel("Adaptyv KD (nM)"); ax.set_ylabel("Twist KD (nM)")
   ax.set_title(f"Cross-vendor KD, n={len(kd)}, Spearman rho={rho:.2f}")
   plt.tight_layout(); plt.show()
   med_ratio = np.median(kd.twist_kd_nM / kd.adaptyv_kd_nM)
   print(f"Median KD ratio (Twist/Adaptyv): {med_ratio:.2f}x  -> systematic format offset, "
         "so treat absolute KD across vendors as ordinal, not interchangeable.")

head(8, "EXPRESSION CONFOUND")
if "twist_expression_mg_per_mL" in ev.columns:
   g = ev.dropna(subset=["twist_expression_mg_per_mL"])
   a = g.loc[g.y == 1, "twist_expression_mg_per_mL"]
   b = g.loc[g.y == 0, "twist_expression_mg_per_mL"]
   if len(a) > 5 and len(b) > 5:
       u, pv = stats.mannwhitneyu(a, b)
       print(f"Titer (mg/mL)  binders median {a.median():.2f} (n={len(a)})  |  "
             f"non-binders {b.median():.2f} (n={len(b)})   Mann-Whitney p={pv:.2e}")
   r = auc_ci(g.y, g.twist_expression_mg_per_mL)
   print(f"AUC of raw expression titer alone as a 'binder' predictor: {r['auc']:.3f}")
   fig, axes = plt.subplots(1, 2, figsize=(9, 3.6))
   axes[0].hist([b, a], bins=25, label=["non-binder", "binder"],
                color=["#BBBBBB", "#4C72B0"], density=True)
   axes[0].set_xlabel("Twist titer (mg/mL)"); axes[0].set_ylabel("density"); axes[0].legend()
   axes[0].set_title("Expression by outcome")
   if "adaptyv_expression" in ev.columns:
       ex = ev.groupby(ev.adaptyv_expression.astype(str)).y.agg(["mean", "size"])
       ex = ex[ex["size"] >= 10].sort_values("mean")
       axes[1].barh(ex.index, ex["mean"], color="#DD8452")
       axes[1].set_xlabel("hit rate"); axes[1].set_title("Hit rate by Adaptyv expression class")
   plt.tight_layout(); plt.show()
print("\nTakeaway: if expression alone scores meaningfully above 0.5, then part of "
     "every AUC in section 4 is a solubility signal riding along. To isolate "
     "interface quality, re-run section 4 restricted to designs that expressed.")
expressed = ev[ev.adaptyv_expression.astype(str).isin(["medium", "high"])] if "adaptyv_expression" in ev.columns else ev
if len(expressed) > 100:
   r_all = auc_ci(ev.y, ev.cons_all)
   r_exp = auc_ci(expressed.y, expressed.cons_all)
   print(f"  consensus AUC, all evaluable   : {r_all['auc']:.3f} (n={r_all['n']})")
   print(f"  consensus AUC, expressed only  : {r_exp['auc']:.3f} (n={r_exp['n']})")
head(9, "EPITOPE CONVERGENCE")
def parse_epitope(s):
   if not isinstance(s, str) or not s.strip():
       return frozenset()
   out = set()
   for tok in s.split(";"):
       tok = tok.strip()
       if not tok:
           continue
       out.add(tok.split(":")[-1])
   return frozenset(out)
ev["epi"] = ev["epitope_residues"].apply(parse_epitope)
def mean_pairwise_jaccard(sets, max_pairs=4000, seed=SEED):
   sets = [s for s in sets if len(s) > 0]
   if len(sets) < 2:
       return np.nan
   pairs = list(itertools.combinations(range(len(sets)), 2))
   rng = np.random.default_rng(seed)
   if len(pairs) > max_pairs:
       pairs = [pairs[i] for i in rng.choice(len(pairs), max_pairs, replace=False)]
   vals = []
   for i, j in pairs:
       u = len(sets[i] | sets[j])
       vals.append(len(sets[i] & sets[j]) / u if u else 0.0)
   return float(np.mean(vals))
rows = []
for tgt, g in ev.groupby("target"):
   B = g.loc[g.y == 1, "epi"].tolist()
   N = g.loc[g.y == 0, "epi"].tolist()
   if len(B) >= 3 and len(N) >= 3:
       rows.append({"target": tgt, "n_bind": len(B), "n_non": len(N),
                    "J_binders": mean_pairwise_jaccard(B),
                    "J_nonbinders": mean_pairwise_jaccard(N)})
epi = pd.DataFrame(rows)
if len(epi):
   epi["delta"] = epi.J_binders - epi.J_nonbinders
   print(epi.sort_values("delta", ascending=False).to_string(
       index=False, formatters={c: "{:.3f}".format for c in ["J_binders", "J_nonbinders", "delta"]}))
   w = stats.wilcoxon(epi.J_binders, epi.J_nonbinders) if len(epi) >= 6 else None
   if w:
       print(f"\nPaired Wilcoxon across targets: p={w.pvalue:.4f}  "
             f"(binders more epitope-convergent than failures?)")
   tgt = epi.sort_values("n_bind", ascending=False).target.iloc[0]
   sub = ev[ev.target == tgt]
   freq_b = pd.Series([r for s in sub[sub.y == 1].epi for r in s]).value_counts()
   freq_n = pd.Series([r for s in sub[sub.y == 0].epi for r in s]).value_counts()
   top = freq_b.head(18).index
   fig, ax = plt.subplots(figsize=(9, 3.8))
   xx = np.arange(len(top))
   ax.bar(xx - 0.2, (freq_b.reindex(top).fillna(0) / max(1, (sub.y == 1).sum())), 0.4,
          label="binders", color="#4C72B0")
   ax.bar(xx + 0.2, (freq_n.reindex(top).fillna(0) / max(1, (sub.y == 0).sum())), 0.4,
          label="non-binders", color="#BBBBBB")
   ax.set_xticks(xx); ax.set_xticklabels(top, rotation=70, ha="right")
   ax.set_ylabel("fraction of designs contacting"); ax.set_title(f"Epitope usage on {tgt}")
   ax.legend(); plt.tight_layout(); plt.show()

head(10, "MODELLING WITH HONEST CROSS-VALIDATION")
AAS = "ACDEFGHIKLMNPQRSTVWY"
KD_HYDRO = dict(zip(AAS, [1.8, 2.5, -3.5, -3.5, 2.8, -0.4, -3.2, 4.5, -3.9, 3.8,
                         1.9, -3.5, -1.6, -3.5, -4.5, -0.8, -0.7, 4.2, -0.9, -1.3]))
CHARGE = {"K": 1, "R": 1, "H": 0.1, "D": -1, "E": -1}
def seq_features(seq):
   seq = "".join(ch for ch in str(seq).upper() if ch in AAS)
   L = max(1, len(seq))
   counts = {a: seq.count(a) / L for a in AAS}
   f = {f"aa_{a}": counts[a] for a in AAS}
   f["length"] = len(seq)
   f["net_charge"] = sum(CHARGE.get(c, 0) for c in seq)
   f["charge_density"] = f["net_charge"] / L
   f["gravy"] = float(np.mean([KD_HYDRO[c] for c in seq])) if seq else 0.0
   f["aromatic"] = sum(counts[a] for a in "FWY")
   f["helix_prone"] = sum(counts[a] for a in "AELM")
   f["beta_prone"] = sum(counts[a] for a in "VIYFT")
   f["gly_pro"] = counts["G"] + counts["P"]
   p = np.array([counts[a] for a in AAS]); p = p[p > 0]
   f["entropy"] = float(-(p * np.log2(p)).sum())
   run, best = 0, 0
   for c in seq:
       run = run + 1 if KD_HYDRO[c] > 1.5 else 0
       best = max(best, run)
   f["max_hydrophobic_run"] = best
   return f
SF = pd.DataFrame([seq_features(s) for s in ev.sequence], index=ev.index)
seq_cols = list(SF.columns)
sil_cols = [c for c in ev.columns if c.startswith(("ipsae_min_", "sc_dockq_"))] + \
          ["cons_all", "cons_min", "cons_disagree"]
meta_cols = [c for c in ["rank", "n_optimization_rounds", "epitope_n_residues"] if c in ev.columns]
X_all = pd.concat([ev[sil_cols + meta_cols], SF], axis=1)
y = ev.y.values
groups = ev.target.values
FEATURE_SETS = {
   "in-silico only": sil_cols + meta_cols,
   "sequence only": seq_cols,
   "in-silico + sequence": sil_cols + meta_cols + seq_cols,
}
def cv_auc(X, y, splitter, groups=None):
   aucs = []
   it = splitter.split(X, y, groups) if groups is not None else splitter.split(X, y)
   for tr, te in it:
       if len(np.unique(y[te])) < 2:
           continue
       m = HistGradientBoostingClassifier(max_depth=4, max_iter=250,
                                          learning_rate=0.06, random_state=SEED)
       m.fit(X.iloc[tr], y[tr])
       aucs.append(roc_auc_score(y[te], m.predict_proba(X.iloc[te])[:, 1]))
   return float(np.mean(aucs)), float(np.std(aucs)), len(aucs)
print(f"{'feature set':24s} {'random 5-fold':>18s} {'grouped-by-target':>20s}")
print("-" * 66)
results = {}
for name, cols in FEATURE_SETS.items():
   X = X_all[cols]
   r_mean, r_sd, _ = cv_auc(X, y, StratifiedKFold(5, shuffle=True, random_state=SEED))
   g_mean, g_sd, nf = cv_auc(X, y, GroupKFold(n_splits=5), groups=groups)
   results[name] = (r_mean, g_mean)
   print(f"{name:24s} {r_mean:.3f} +/- {r_sd:.3f}   {g_mean:.3f} +/- {g_sd:.3f}")
gap = results["in-silico + sequence"][0] - results["in-silico + sequence"][1]
print(f"\nRandom-CV minus grouped-CV for the full feature set: {gap:+.3f}")
print("That gap is leakage: features that encode target identity (epitope size, "
     "length priors, generator habits) let a randomly-split model recover the "
     "per-target base rate instead of learning what makes a binder. Report the "
     "grouped number; the random one is what a target-blind reviewer will catch.")
Xt = pd.get_dummies(pd.Series(groups, index=ev.index), prefix="tgt")
r_mean, _, _ = cv_auc(Xt, y, StratifiedKFold(5, shuffle=True, random_state=SEED))
print(f"\nControl - target one-hot ONLY, random CV: AUC={r_mean:.3f} "
     "(pure base-rate memorisation, zero design signal).")
gkf = GroupKFold(n_splits=5)
tr, te = next(iter(gkf.split(X_all, y, groups)))
model = HistGradientBoostingClassifier(max_depth=4, max_iter=250,
                                      learning_rate=0.06, random_state=SEED).fit(
   X_all[FEATURE_SETS["in-silico + sequence"]].iloc[tr], y[tr])
imp = permutation_importance(model, X_all[FEATURE_SETS["in-silico + sequence"]].iloc[te],
                            y[te], n_repeats=12, random_state=SEED, scoring="roc_auc")
order = np.argsort(imp.importances_mean)[-18:]
names = np.array(FEATURE_SETS["in-silico + sequence"])[order]
fig, ax = plt.subplots(figsize=(7, 5))
ax.barh(names, imp.importances_mean[order],
       xerr=imp.importances_std[order], color="#55A868")
ax.set_xlabel("drop in AUC when permuted")
ax.set_title("Permutation importance (held-out target block)")
plt.tight_layout(); plt.show()
head("", "SUMMARY")
print(f"""
Evaluable designs        : {len(ev):,}   base hit rate {ev.y.mean():.3f}
Best single in-silico    : {best_col}  AUC {best_single.auc:.3f}
Rank-average consensus   : AUC {auc_ci(ev.y, ev.cons_all)['auc']:.3f}
Honest ML (grouped CV)   : AUC {results['in-silico + sequence'][1]:.3f}   <- the one to report
Same model, random CV    : AUC {results['in-silico + sequence'][0]:.3f}   (gap = {gap:+.3f} leakage)
Five things this dataset teaches that a design paper usually cannot:
  1. Target identity dominates every other factor; always stratify.
  2. Structure-predictor confidence is real but weak signal (AUC ~0.6-0.75),
     nowhere near a standalone go/no-go filter.
  3. Ensembling across predictors is a cheap, reliable few-points-of-AUC win.
  4. Cross-vendor label noise caps how high any AUC here can honestly go.
  5. Expression failure masquerades as binding failure. Condition on it.
Extensions worth trying:
  - load_table('insilico_cofold_predictions') for all 5 seeds/predictor, and
    test whether seed VARIANCE beats seed-best as a confidence signal
  - load_table('adaptyv_fit_curves') to refit kinetics yourself and flag
    designs whose reported KD rests on a poorly-conditioned fit
  - load_table('insilico_provenance_steps') to relate optimisation-round count
    to eventual success
  - snapshot_download(..., allow_patterns='data/designs/EGFR/<name>/*') for
    mmCIF structures + PAE matrices on a single design
""")