From In-Silico to Wet-Lab: Evaluating AI Protein Design Performance
In this tutorial, we use Anthropic’s claude-protein-binder-design dataset, which incorporates 1,440 AI-designed miniprotein binders examined in opposition to 16 targets. Because the discharge contains each computational predictions and actual wet-lab outcomes from two impartial labs, we will transcend merely finding out the designs. We consider how nicely construction predictors determine profitable binders, whether or not combining predictions improves efficiency, how rankings translate into sensible testing budgets, and the way a lot disagreement comes from the assays themselves. Also, we prepare a target-aware classifier to check whether or not these alerts can reliably predict experimental success.
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("putting in:", ", ".be part of(_missing))
subprocess.run([sys.executable, "-m", "pip", "install", "-q", *_missing], examine=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("show.width", 200)
pd.set_option("show.max_columns", 100)
plt.rcParams.replace({"determine.dpi": 110, "font.dimension": 9, "axes.grid": True,
"grid.alpha": 0.25, "axes.spines.high": False, "axes.spines.proper": 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("knowledge/tables/") and f.endswith(".parquet"):
key = f[len("data/tables/"): -len(".parquet")].change("/", "_")
TABLES[key] = f
print(f"Found {len(TABLES)} Parquet tables:")
for okay in sorted(TABLES):
print(f" - {okay:38s} {TABLES[k]}")
def load_table(title: str) -> pd.DataFrame:
"""Load a subset by its viewer title, with a datasets-library fallback."""
if title in TABLES:
return pd.read_parquet(hf_hub_download(REPO, TABLES[name], repo_type="dataset"))
from datasets import load_dataset
return load_dataset(REPO, title, break up="full").to_pandas()
ds = load_table("design_summary")
print(f"ndesign_summary: {ds.form[0]:,} rows x {ds.form[1]} columns")
We begin by putting in solely what the runtime is definitely lacking, then enumerate the repository as soon as and construct a {subset to path} map as an alternative of hard-coding file areas. This issues as a result of the naming shouldn’t be uniform; the subset wetlab_summary lives at knowledge/tables/wetlab/abstract.parquet, and a guessed path would fail silently. With the map in place we pull design_summary, one row per design, 1,440 rows broad sufficient to carry each be part of we want downstream.
head(2, "SCHEMA + EVALUABLE SET")
CALLS = {"binder", "non_binder"}
examined = 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 name): {len(ev):,}")
print(f"Confirmed binders : {int(ev['y'].sum()):,} "
f"({100 * ev['y'].imply():.1f}% base fee)")
print(f"Never measured : {len(ds) - len(ev):,}")
print("nCategorical ranges:")
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)}): {', '.be part of(vals.index[:6])}"
+ (" ..." if len(vals) > 6 else ""))
print(f"nTargets ({ds['target'].nunique()}): {', '.be part of(sorted(ds['target'].distinctive()))}")
print(f"Binder size: {ds.binder_length.min()}-{ds.binder_length.max()} aa "
f"(median {ds.binder_length.median():.0f})")
head(3, "HIT-RATE LANDSCAPE")
def wilson(okay, n, z=1.96):
if n == 0:
return (np.nan, np.nan, np.nan)
p = okay / 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, hello = wilson(int(g.y.sum()), len(g))
rows.append({by: key, "n": len(g), "hits": int(g.y.sum()),
"fee": p, "lo": lo, "hello": hello})
return pd.DataFrame(rows).sort_values("fee", 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 fee by {dim} ---")
print(t.to_string(index=False,
formatters={"fee": "{:.3f}".format, "lo": "{:.3f}".format, "hello": "{:.3f}".format}))
tt = rate_table(ev, "goal")
fig, ax = plt.subplots(figsize=(9, 4.2))
ax.bar(tt.goal, tt.fee, colour="#4C72B0")
ax.errorbar(tt.goal, tt.fee,
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.imply(), ls="--", c="crimson", lw=1, label=f"pooled {ev.y.imply():.2f}")
ax.set_ylabel("experimental hit fee"); ax.set_title("Hit fee by goal (Wilson 95% CI)")
ax.tick_params(axis="x", rotation=55); ax.legend(); plt.tight_layout(); plt.present()
print("nReadvert this plot because the dominant impact dimension within the dataset: goal alternative "
"swamps generator alternative. Any mannequin comparability that doesn't stratify by "
"goal is usually measuring which targets that mannequin was pointed at.")
We outline the evaluable set by filtering on precise vendor calls slightly than on binder_final, as a result of that column is a bool and so data the 120 never-measured designs as False slightly than lacking. From there we compute hit charges by mannequin, marketing campaign, generator, and goal, wrapping every in a Wilson interval since a number of subgroups sit within the small-n regime the place the conventional approximation misbehaves. The goal plot is the one to learn first: it exhibits antigen alternative swamping each different issue we examine.
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)}): {', '.be part of(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.distinctive(y)) < 2:
return dict(auc=np.nan, lo=np.nan, hello=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 vary(n_boot):
b = rng.alternative(idx, len(idx), change=True)
if len(np.distinctive(y[b])) > 1:
boots.append(roc_auc_score(y[b], s[b]))
lo, hello = (np.percentile(boots, [2.5, 97.5]) if boots else (np.nan, np.nan))
return dict(auc=base, lo=lo, hello=hello, 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(decrease=0).fillna(0)
hi_err = (sub.hello - sub.auc).clip(decrease=0).fillna(0)
ax.bar(x + (i - 0.5) * w, sub.auc, w, label=metric, colour=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="proper")
ax.set_ylabel("AUC"); ax.set_ylim(0.35, None)
ax.set_title("In-silico rating vs wet-lab binding, by construction predictor")
ax.legend(); plt.tight_layout(); plt.present()
print("Interpretation: AUCs land nicely above likelihood however far beneath the ~0.9 you "
"would wish to belief a single filter. That hole is the whole sensible "
"purpose 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="hold")
R_ips, R_dq = pct_rank(ev, ips_cols), pct_rank(ev, dq_cols)
ev["cons_ipsae"] = R_ips.imply(axis=1)
ev["cons_dockq"] = R_dq.imply(axis=1)
ev["cons_all"] = pd.concat([R_ips, R_dq], axis=1).imply(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 title in ["cons_ipsae", "cons_dockq", "cons_all", "cons_median", "cons_min", "cons_disagree"]:
r = auc_ci(ev.y, ev[name])
print(f" {title:16s} AUC={r['auc']:.3f} [{r['lo']:.3f}, {r['hi']:.3f}] AP={r['ap']:.3f}")
corr = ev[ips_cols].corr(methodology="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(vary(len(lbl))); ax.set_xticklabels(lbl, rotation=90)
ax.set_yticks(vary(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.present()
print("nIf each off-diagonal cell have been ~1.0 there could be no ensemble achieve to "
"harvest. The average correlations are why cons_all usually edges out "
"the very best single predictor — and why disagreement itself carries sign.")
We rating all ten predictors in opposition to the wet-lab label, on each ipSAE and self-consistency DockQ, with bootstrapped confidence intervals so we will see which variations are actual. We then rank-normalize every column to percentiles and combination them, which retains the comparability scale-free throughout metrics that dwell on completely different ranges and pile up otherwise at zero. The Spearman heatmap explains why the ensemble helps in any respect; if the predictors agreed completely there could be nothing left to harvest.
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)
okay = 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, fashion 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, fashion, lw=1.8, label=lab)
ax.axhline(ev.y.imply(), ls=":", c="crimson", lw=1.4, label=f"random baseline ({ev.y.imply():.2f})")
ax.set_xlabel("designs ordered for wet-lab testing (N, best-first)")
ax.set_ylabel("hit fee amongst high N"); ax.set_title("How a lot does in-silico triage purchase you?")
ax.legend(); plt.tight_layout(); plt.present()
print("Enrichment at small budgets:")
for N in [25, 50, 100, 200]:
line = f" N={N:4d} | random {ev.y.imply():.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.imply():.2f}x)"
print(line)
head(7, "VENDOR CONCORDANCE")
each = ev[ev.adaptyv_binding.isin(CALLS) & ev.twist_binding.isin(CALLS)]
ct = pd.crosstab(each.adaptyv_binding, each.twist_binding)
print(f"Designs with calls from BOTH distributors: {len(each):,}n")
print(ct.to_string())
if len(each) > 10:
kappa = cohen_kappa_score(each.adaptyv_binding, each.twist_binding)
agree = (each.adaptyv_binding == each.twist_binding).imply()
print(f"nRaw settlement: {agree:.3f} Cohen's kappa: {kappa:.3f}")
print("Kappa nicely beneath 1.0 means a part of the 'unpredictable' variance above "
"is assay disagreement, not mannequin 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.present()
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 deal with absolute KD throughout distributors as ordinal, not interchangeable.")
We convert rating efficiency into precision@N, as a result of no lab orders 1,300 constructs and AUC quietly hides how a rating behaves on the high of the listing. The enrichment desk then tells us what triage really buys at budgets of 25, 50, 100, and 200. We comply with it with Cohen’s κ and a log-log KD comparability between distributors, which units the ceiling: label noise bounds how excessive any AUC above can truthfully climb.
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 uncooked 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"],
colour=["#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 consequence")
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("imply")
axes[1].barh(ex.index, ex["mean"], colour="#DD8452")
axes[1].set_xlabel("hit fee"); axes[1].set_title("Hit fee by Adaptyv expression class")
plt.tight_layout(); plt.present()
print("nTakeaway: if expression alone scores meaningfully above 0.5, then a part of "
"each AUC in part 4 is a solubility sign using alongside. To isolate "
"interface high quality, re-run part 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 solely : {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.break up(";"):
tok = tok.strip()
if not tok:
proceed
out.add(tok.break up(":")[-1])
return frozenset(out)
ev["epi"] = ev["epitope_residues"].apply(parse_epitope)
def mean_pairwise_jaccard(units, max_pairs=4000, seed=SEED):
units = [s for s in sets if len(s) > 0]
if len(units) < 2:
return np.nan
pairs = listing(itertools.mixtures(vary(len(units)), 2))
rng = np.random.default_rng(seed)
if len(pairs) > max_pairs:
pairs = [pairs[i] for i in rng.alternative(len(pairs), max_pairs, change=False)]
vals = []
for i, j in pairs:
u = len(units[i] | units[j])
vals.append(len(units[i] & units[j]) / u if u else 0.0)
return float(np.imply(vals))
rows = []
for tgt, g in ev.groupby("goal"):
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({"goal": 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 throughout targets: p={w.pvalue:.4f} "
f"(binders extra epitope-convergent than failures?)")
tgt = epi.sort_values("n_bind", ascending=False).goal.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()
high = freq_b.head(18).index
fig, ax = plt.subplots(figsize=(9, 3.8))
xx = np.arange(len(high))
ax.bar(xx - 0.2, (freq_b.reindex(high).fillna(0) / max(1, (sub.y == 1).sum())), 0.4,
label="binders", colour="#4C72B0")
ax.bar(xx + 0.2, (freq_n.reindex(high).fillna(0) / max(1, (sub.y == 0).sum())), 0.4,
label="non-binders", colour="#BBBBBB")
ax.set_xticks(xx); ax.set_xticklabels(high, rotation=70, ha="proper")
ax.set_ylabel("fraction of designs contacting"); ax.set_title(f"Epitope utilization on {tgt}")
ax.legend(); plt.tight_layout(); plt.present()
We check whether or not expression titer alone discriminates binders, and if it does, we all know a part of each rating from above is solubility using alongside beneath one other title. Re-running consensus on expressed-only designs isolates interface high quality from biophysics. We then parse the epitope contact lists into residue units and ask, per goal and paired throughout targets, whether or not confirmed binders converge on a shared patch greater than the failures do.
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 = {"Ok": 1, "R": 1, "H": 0.1, "D": -1, "E": -1}
def seq_features(seq):
seq = "".be part of(ch for ch in str(seq).higher() if ch in AAS)
L = max(1, len(seq))
counts = {a: seq.rely(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.imply([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, finest = 0, 0
for c in seq:
run = run + 1 if KD_HYDRO[c] > 1.5 else 0
finest = max(finest, run)
f["max_hydrophobic_run"] = finest
return f
SF = pd.DataFrame([seq_features(s) for s in ev.sequence], index=ev.index)
seq_cols = listing(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
teams = ev.goal.values
FEATURE_SETS = {
"in-silico solely": sil_cols + meta_cols,
"sequence solely": seq_cols,
"in-silico + sequence": sil_cols + meta_cols + seq_cols,
}
def cv_auc(X, y, splitter, teams=None):
aucs = []
it = splitter.break up(X, y, teams) if teams shouldn't be None else splitter.break up(X, y)
for tr, te in it:
if len(np.distinctive(y[te])) < 2:
proceed
m = HistGradientBoostingClassifier(max_depth=4, max_iter=250,
learning_rate=0.06, random_state=SEED)
m.match(X.iloc[tr], y[tr])
aucs.append(roc_auc_score(y[te], m.predict_proba(X.iloc[te])[:, 1]))
return float(np.imply(aucs)), float(np.std(aucs)), len(aucs)
print(f"{'characteristic set':24s} {'random 5-fold':>18s} {'grouped-by-target':>20s}")
print("-" * 66)
outcomes = {}
for title, cols in FEATURE_SETS.objects():
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), teams=teams)
outcomes[name] = (r_mean, g_mean)
print(f"{title:24s} {r_mean:.3f} +/- {r_sd:.3f} {g_mean:.3f} +/- {g_sd:.3f}")
hole = outcomes["in-silico + sequence"][0] - outcomes["in-silico + sequence"][1]
print(f"nRandom-CV minus grouped-CV for the complete characteristic set: {hole:+.3f}")
print("That hole is leakage: options that encode goal id (epitope dimension, "
"size priors, generator habits) let a randomly-split mannequin recuperate the "
"per-target base fee as an alternative of studying what makes a binder. Report the "
"grouped quantity; the random one is what a target-blind reviewer will catch.")
Xt = pd.get_dummies(pd.Series(teams, index=ev.index), prefix="tgt")
r_mean, _, _ = cv_auc(Xt, y, StratifiedKFold(5, shuffle=True, random_state=SEED))
print(f"nControl - goal one-hot ONLY, random CV: AUC={r_mean:.3f} "
"(pure base-rate memorisation, zero design sign).")
gkf = GroupKFold(n_splits=5)
tr, te = subsequent(iter(gkf.break up(X_all, y, teams)))
mannequin = HistGradientBoostingClassifier(max_depth=4, max_iter=250,
learning_rate=0.06, random_state=SEED).match(
X_all[FEATURE_SETS["in-silico + sequence"]].iloc[tr], y[tr])
imp = permutation_importance(mannequin, 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], colour="#55A868")
ax.set_xlabel("drop in AUC when permuted")
ax.set_title("Permutation significance (held-out goal block)")
plt.tight_layout(); plt.present()
head("", "SUMMARY")
print(f"""
Evaluable designs : {len(ev):,} base hit fee {ev.y.imply():.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 {outcomes['in-silico + sequence'][1]:.3f} <- the one to report
Same mannequin, random CV : AUC {outcomes['in-silico + sequence'][0]:.3f} (hole = {hole:+.3f} leakage)
Five issues this dataset teaches {that a} design paper normally can not:
1. Target id dominates each different issue; at all times stratify.
2. Structure-predictor confidence is actual however weak sign (AUC ~0.6-0.75),
nowhere close to a standalone go/no-go filter.
3. Ensembling throughout predictors is an affordable, dependable few-points-of-AUC win.
4. Cross-vendor label noise caps how excessive any AUC right here can truthfully go.
5. Expression failure masquerades as binding failure. Condition on it.
Extensions value making an attempt:
- load_table('insilico_cofold_predictions') for all 5 seeds/predictor, and
check whether or not seed VARIANCE beats seed-best as a confidence sign
- load_table('adaptyv_fit_curves') to refit kinetics your self and flag
designs whose reported KD rests on a poorly-conditioned match
- load_table('insilico_provenance_steps') to relate optimisation-round rely
to eventual success
- snapshot_download(..., allow_patterns='knowledge/designs/EGFR/<title>/*') for
mmCIF buildings + PAE matrices on a single design
""")
We featurize sequences by composition, cost, hydropathy, entropy, and hydrophobic run size, then match gradient boosting beneath two schemes: random folds and target-grouped folds. The hole between them is the leakage, since designs nest inside targets with very completely different base charges and a random break up lets the mannequin memorize which antigens are simple. The one-hot management makes that express, and permutation significance on a held-out goal block exhibits what survives once we take away the shortcut.
In conclusion, in-silico scoring was useful, however it didn’t inform the entire story. The goal strongly influenced the outcomes, so evaluating fashions with out accounting for it might simply give us a deceptive image. The construction predictors confirmed helpful alerts, and mixing them gave a modest enchancment, however they have been nonetheless not dependable sufficient to use on their very own. We additionally discovered that variations between experiments and protein expression might make a design appear like a binding failure even when the actual concern was poor expression. Overall, we realized that cautious analysis mattered greater than chasing spectacular particular person metrics. By grouping our cross-validation by goal, we obtained a extra reasonable view of how nicely the fashions might generalize to new targets.
Check out the FULL CODES here. Also, be happy to comply with us on Twitter and don’t neglect to be part of our 150k+ML SubReddit and Subscribe to our Newsletter. Wait! are you on telegram? now you can join us on telegram as well.
Need to associate with us for selling your GitHub Repo OR Hugging Face Page OR Product Release OR Webinar and so on.? Connect with us
The submit From In-Silico to Wet-Lab: Evaluating AI Protein Design Performance appeared first on MarkTechPost.
