File size: 3,692 Bytes
31d1d63 | 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 | import json
from collections import defaultdict
from pathlib import Path
import numpy as np
from scipy.stats import spearmanr
RS = Path(__file__).resolve().parents[2] / "data" / "p53" / "rmsd_screen"
COVER = 0.9
N_SPLITS = 400
rng = np.random.default_rng(0)
def load_matrix():
manifest = json.loads((RS / "receptors.json").read_text())
recs = {(r["state"], r["frame"]): r for r in manifest}
per = defaultdict(dict)
for f in sorted((RS / "scores").glob("*.json")):
d = json.loads(f.read_text())
for s in d["scores"]:
if s.get("affinity") is None:
continue
r = recs.get((s.get("state"), s.get("frame")))
if r:
per[d["ligand"]][r["frame"]] = float(s["affinity"])
ligs_all = sorted(per)
all_frames = sorted({r["frame"] for r in manifest})
cov = {i: sum(i in per[l] for l in ligs_all) / len(ligs_all) for i in all_frames}
frames = [i for i in all_frames if cov[i] >= COVER]
lig_cov = {l: sum(i in per[l] for i in frames) / len(frames) for l in ligs_all}
ligs = [l for l in ligs_all if lig_cov[l] >= COVER]
rec_by_frame = {r["frame"]: r for r in manifest}
vol = np.array([rec_by_frame[i]["site_volume"] for i in frames], dtype=float)
S = np.full((len(ligs), len(frames)), np.nan)
for li, l in enumerate(ligs):
for ri, i in enumerate(frames):
if i in per[l]:
S[li, ri] = per[l][i]
row_mean = np.nanmean(S, axis=1, keepdims=True)
S = np.where(np.isnan(S), row_mean, S)
return S, ligs, frames, vol, rec_by_frame, manifest
def volume_trend(Z, vol, nbins=10):
"""Return the pooled volume-decile trend for each receptor."""
edges = np.quantile(vol, np.linspace(0, 1, nbins + 1))
edges[-1] += 1e-9
bin_id = np.clip(np.searchsorted(edges, vol, side="right") - 1, 0, nbins - 1)
g = np.array(
[Z[:, bin_id == b].mean() if (bin_id == b).any() else 0.0 for b in range(nbins)]
)
return g[bin_id]
def gamma_matrix():
"""Exactly the manuscript's gamma: S -> z-score by ligand -> remove f(V_r) -> remove beta_r."""
S, ligs, frames, vol, rec_by_frame, manifest = load_matrix()
Z = (S - S.mean(axis=1, keepdims=True)) / S.std(axis=1, keepdims=True)
G = volume_trend(Z, vol)
R = Z - G[None, :]
gamma = R - R.mean(axis=0)[None, :]
strat = np.array([rec_by_frame[i]["state"] for i in frames])
return gamma, np.array(ligs), np.array(frames), strat
def pref_corr(mat):
"""Ligand-preference profile = ligand x ligand correlation over the given receptors."""
return np.corrcoef(mat)
def split_half_within(gamma, strat, n_splits=N_SPLITS, shuffle_strata=False):
"""Split receptors into halves INSIDE each stratum, then compare ligand-preference
structure between the halves. Between-stratum contrasts cannot drive this."""
n_lig = gamma.shape[0]
iu = np.triu_indices(n_lig, k=1)
out = []
for _ in range(n_splits):
s = strat.copy()
if shuffle_strata:
s = rng.permutation(s)
A, B = [], []
for lab in np.unique(s):
idx = np.where(s == lab)[0]
perm = rng.permutation(idx)
A.extend(perm[: len(idx) // 2])
B.extend(perm[len(idx) // 2 :])
A, B = np.array(A), np.array(B)
if len(A) < 3 or len(B) < 3:
continue
gc = gamma.copy()
for lab in np.unique(s):
m = s == lab
gc[:, m] -= gc[:, m].mean(axis=1, keepdims=True)
r_, _ = spearmanr(pref_corr(gc[:, A])[iu], pref_corr(gc[:, B])[iu])
if np.isfinite(r_):
out.append(r_)
return np.array(out)
|