From a8dcca549f8c498d4140e7e3720658bef5126e41 Mon Sep 17 00:00:00 2001 From: Hossam Mahmoud Date: Sat, 26 Sep 2026 05:57:20 +0300 Subject: [PATCH 1/6] Complex analyses, annotations and analysis-window selection - contact_lifetime, water_bridges, porcupine ported from the validated legacy pipeline, generalised (partners from persisted chain identity; porcupine uses its own aligned copy and an SVD instead of a 3N x 3N covariance). - core/annotations.py: display names, biological numbering offsets, domain transfer by global alignment (BLOSUM62, gap -10/-0.5) with uncertain-edge flags, motif location (exact, then local alignment, partial matches). - annotation analysis: domains/motifs/docking site in biological numbering and their interface involvement; says so when no annotations are given. - analysis_window: Chodera equilibration detection, drift and half-vs-half tests; no silent choice of the primary averaging window. - statistics.autocorr: statistical inefficiency, corrected SEM, detection. Co-Authored-By: Claude Opus 5.5 --- moldynx/analysis/analysis_window.py | 103 ++++++++++++++++++ moldynx/analysis/annotation.py | 120 +++++++++++++++++++++ moldynx/analysis/contact_lifetime.py | 115 ++++++++++++++++++++ moldynx/analysis/porcupine.py | 81 ++++++++++++++ moldynx/analysis/water_bridges.py | 76 +++++++++++++ moldynx/core/annotations.py | 154 +++++++++++++++++++++++++++ moldynx/statistics/__init__.py | 3 + moldynx/statistics/autocorr.py | 73 +++++++++++++ tests/test_interface.py | 8 +- tests/test_regions_and_window.py | 102 ++++++++++++++++++ 10 files changed, 833 insertions(+), 2 deletions(-) create mode 100644 moldynx/analysis/analysis_window.py create mode 100644 moldynx/analysis/annotation.py create mode 100644 moldynx/analysis/contact_lifetime.py create mode 100644 moldynx/analysis/porcupine.py create mode 100644 moldynx/analysis/water_bridges.py create mode 100644 moldynx/core/annotations.py create mode 100644 moldynx/statistics/autocorr.py create mode 100644 tests/test_regions_and_window.py diff --git a/moldynx/analysis/analysis_window.py b/moldynx/analysis/analysis_window.py new file mode 100644 index 0000000..46a43d1 --- /dev/null +++ b/moldynx/analysis/analysis_window.py @@ -0,0 +1,103 @@ +""" +Is there a stationary window for binding-related averages? (Chodera equilibration detection) + +For each binding-relevant observable (minimum distance, contacts, interface residues, +buried area from ``interface``; RMSD columns from ``rmsd``): the start t₀ that +maximises the number of effectively independent samples, then -- inside [t₀, end] -- +the linear drift and a half-vs-half Welch test. The conservative t₀ is the latest +over the binding observables. + +If no window is stationary the result says so; the primary window for averages +(e.g. MM-GBSA) is then a **user decision** (``binding_energy.primary_window`` in the +config), never chosen silently. The proposal offered is the full run plus the +final 20 %. +""" + +from __future__ import annotations + +import json + +import numpy as np +import pandas as pd + +from moldynx.core.base import BaseAnalysis +from moldynx.statistics import detect_equilibration, drift + +BINDING = {"min_interface_dist_nm": "Minimum inter-partner distance (nm)", + "n_inter_contacts": "Residue–residue contacts", + "n_interface_residues": "Interface residues", + "buried_area_nm2": "Buried area (nm²)"} + + +class AnalysisWindow(BaseAnalysis): + name = "analysis_window" + label = "Stationarity and analysis-window selection" + category = "quality" + required_files = {"trajectory", "topology"} + supported_systems = {"*"} + order = 160 + default_params = {"alpha": 0.01} + outputs = ["results/window_selection.json", "results/window_selection.csv"] + + def run(self, ctx) -> dict: + p = self.params(ctx) + series: dict[str, tuple[np.ndarray, np.ndarray, bool]] = {} + f = ctx.csv_path("interface_timeseries.csv") + if f.exists(): + df = pd.read_csv(f) + for col, lab in BINDING.items(): + if col in df and df[col].notna().sum() > 10: + ok = df[col].notna() + series[lab] = (df.time_ns[ok].to_numpy(), df[col][ok].to_numpy(float), True) + f = ctx.csv_path("rmsd.csv") + if f.exists(): + df = pd.read_csv(f) + tcol = "time_ns" if "time_ns" in df else df.columns[0] + for col in df.columns: + if col != tcol and df[col].dtype.kind == "f" and "rmsd" in col.lower(): + series[col] = (df[tcol].to_numpy(), df[col].to_numpy(float), False) + if not series: + return {"status": "skipped", "reason": "no interface or RMSD time series to test"} + + rows = [] + for name, (t, x, binding) in series.items(): + eq = detect_equilibration(x, step=max(1, len(x) // 100)) + i0 = eq["t0"] + d = drift(t[i0:], x[i0:]) + span = float(t[-1] - t[i0]) or 1.0 + rows.append({"observable": name, "binding_relevant": binding, + "t0_ns": float(t[i0]), "stat_ineff": eq["g"], "n_eff": eq["n_eff"], + "mean_after_t0": float(x[i0:].mean()), + "slope_per_10ns": None if d["slope"] is None else d["slope"] * 10, + "p_slope": d["p_slope"], "half_vs_half_p": d["half_p"], + "relative_drift": None if d["slope"] is None else + d["slope"] * span / (abs(x[i0:].mean()) + 1e-12), + "stationary": bool(d["p_slope"] is not None and d["p_slope"] > p["alpha"] + and d["half_p"] > p["alpha"])}) + tab = pd.DataFrame(rows) + tab.to_csv(ctx.csv_path("window_selection.csv"), index=False) + t_all = next(iter(series.values()))[0] + t_end = float(t_all[-1]) + bind = tab[tab.binding_relevant] + t0 = float(bind.t0_ns.max()) if len(bind) else float(tab.t0_ns.max()) + step = max(t_end * 0.05, 1e-9) + t0 = float(np.ceil(t0 / step) * step) + stationary = bool(len(bind) and bind.stationary.all()) + chosen = (ctx.config.binding_energy or {}).get("primary_window", "ask") + result = { + "conservative_t0_ns": t0, "end_ns": t_end, + "binding_observables_stationary": stationary, + "non_stationary": list(tab[~tab.stationary].observable), + "proposal": {"full": [float(t_all[0]), t_end], + "final_20pct": [round(t_end * 0.8, 6), t_end], + "after_t0": [t0, t_end]}, + "primary_window": chosen, + "decision_required": chosen in (None, "ask"), + "statement": ("binding-related observables are stationary after t0" + if stationary else + "no stationary window: report time-resolved values; the primary " + "window for averages is a user decision"), + } + ctx.csv_path("window_selection.json").write_text(json.dumps( + {"summary": result, "observables": rows}, indent=2, default=float), encoding="utf-8") + return result diff --git a/moldynx/analysis/annotation.py b/moldynx/analysis/annotation.py new file mode 100644 index 0000000..bf1831e --- /dev/null +++ b/moldynx/analysis/annotation.py @@ -0,0 +1,120 @@ +""" +Regions of the system from user annotations: domains (transferred by homology), +motifs, docking-site residues -- and how much of each region is at the interface. + +Runs only when the configuration has an ``annotations:`` block; otherwise it says so. +All numbers are reported in biological numbering with the trajectory offset stated. +""" + +from __future__ import annotations + +import json +from pathlib import Path + +import pandas as pd + +from moldynx.core import annotations as ann +from moldynx.core.base import BaseAnalysis + + +def local_to_bio(chain: ann.ChainAnnotation, position: int) -> int: + """1-based position in the chain's sequence -> biological residue number.""" + return chain.bio(chain.resid_first + int(position) - 1) + + +class Annotation(BaseAnalysis): + name = "annotation" + label = "Domains, motifs and docking site (user annotations)" + category = "annotation" + required_files = {"trajectory", "topology"} + supported_systems = {"*"} + order = 150 # after interface (40) and rmsf + outputs = ["results/annotation.json", "results/region_summary.csv"] + + def run(self, ctx) -> dict: + a = ctx.config.annotations or {} + chains = ctx.core_meta.get("chains", []) + ca = ann.chain_annotations(chains, a) + doc: dict = {"chains": {s: {"display": c.display, "role": c.role, + "trajectory_resids": [c.resid_first, c.resid_last], + "biological_range": list(c.bio_range), + "numbering_offset": c.offset} for s, c in ca.items()}, + "domains": {}, "motifs": [], "docking_site": {}, + "unresolved_metadata": a.get("unresolved_metadata", {})} + regions = [] # (segid, kind, name, bio_start, bio_end) + for seg, spec in (a.get("domains") or {}).items(): + if seg not in ca: + doc["domains"][seg] = {"error": f"chain {seg} not in this system"} + continue + ref = spec.get("reference_sequence") or ( + ann.read_fasta_sequence(Path(spec["reference_fasta"])) + if spec.get("reference_fasta") else "") + if not ref: + doc["domains"][seg] = {"error": "no reference_sequence / reference_fasta"} + continue + try: + mapped = ann.transfer_regions(ca[seg].sequence, ref, spec.get("regions", {})) + except ImportError: + doc["domains"][seg] = {"error": "biopython is required for domain transfer"} + continue + for m in mapped: # sequence position -> biological + for k in ("start", "end"): + if m[k]: + m[k] = local_to_bio(ca[seg], m[k]) + doc["domains"][seg] = {"reference": spec.get("reference_name", "reference"), + "regions": mapped} + regions += [(seg, "domain", m["region"], m["start"], m["end"]) for m in mapped + if m["start"] and m["end"]] + for m in a.get("motifs") or []: + segs = [m["chain"]] if m.get("chain") else list(ca) + for seg in segs: + if seg not in ca: + continue + try: + loc = ann.locate_motif(ca[seg].sequence, m["sequence"]) + except ImportError: + loc = {"match": "exact-only search (biopython missing)"} + for k in ("start", "end"): + if loc.get(k): + loc[k] = local_to_bio(ca[seg], loc[k]) + doc["motifs"].append({"name": m["name"], "chain": seg, **loc}) + if loc.get("start"): + regions.append((seg, "motif", m["name"], loc["start"], loc["end"])) + for seg, residues in (a.get("docking_site") or {}).items(): + doc["docking_site"][seg] = sorted(int(r) for r in residues) + + # interface involvement per region (from the interface analysis, if it ran) + rows = [] + res_csv = ctx.csv_path("interface_residues.csv") + iface = pd.read_csv(res_csv) if res_csv.exists() else None + for seg, kind, name, b0, b1 in regions: + row = {"chain": seg, "display": ca[seg].display, "kind": kind, "region": name, + "bio_start": b0, "bio_end": b1, "n_residues": b1 - b0 + 1} + if iface is not None: + sub = iface[(iface.partner == seg) & + iface.resid.between(ca[seg].md(b0), ca[seg].md(b1))] + row.update({"interface_residues": int((sub.interface_occupancy > 0).sum()), + "core_residues": int((sub.interface_occupancy >= 0.5).sum()), + "occupancy_sum": float(sub.interface_occupancy.sum())}) + rows.append(row) + for seg, residues in doc["docking_site"].items(): + if iface is not None and seg in ca: + md_res = [ca[seg].md(r) for r in residues] + sub = iface[(iface.partner == seg) & iface.resid.isin(md_res)] + rows.append({"chain": seg, "display": ca[seg].display, "kind": "docking site", + "region": "docking site", "n_residues": len(residues), + "interface_residues": int((sub.interface_occupancy > 0).sum()), + "core_residues": int((sub.interface_occupancy >= 0.5).sum()), + "occupancy_sum": float(sub.interface_occupancy.sum())}) + pd.DataFrame(rows).to_csv(ctx.csv_path("region_summary.csv"), index=False) + ctx.csv_path("annotation.json").write_text(json.dumps(doc, indent=2, default=str), + encoding="utf-8") + return {"annotated": bool(a), "chains": doc["chains"], + "n_domains": sum(len(d.get("regions", [])) for d in doc["domains"].values()), + "uncertain_domain_edges": [f"{s}:{m['region']}" for s, d in doc["domains"].items() + for m in d.get("regions", []) if m["uncertain"]], + "motifs": [{k: m.get(k) for k in ("name", "chain", "match", "start", "end")} + for m in doc["motifs"]], + "regions": rows, + "note": None if a else "no annotations in the configuration: display names " + "default to segids and no regions are defined"} diff --git a/moldynx/analysis/contact_lifetime.py b/moldynx/analysis/contact_lifetime.py new file mode 100644 index 0000000..0eae23c --- /dev/null +++ b/moldynx/analysis/contact_lifetime.py @@ -0,0 +1,115 @@ +""" +Inter-partner contact lifetimes and native-contact survival (complexes). + +One streaming pass tracks every inter-partner residue–residue contact (heavy atoms +< ``cutoff``): how long each contact survives once formed, how often it re-forms, +and what fraction of the frame-0 contacts are still present over time. +""" + +from __future__ import annotations + +from collections import defaultdict + +import numpy as np +import pandas as pd +from MDAnalysis.lib.distances import capped_distance + +from moldynx import plotting +from moldynx.core.base import BaseAnalysis +from moldynx.core.system import COMPLEX_SYSTEMS +from moldynx.plotting import PALETTE + + +class ContactLifetime(BaseAnalysis): + name = "contact_lifetime" + label = "Contact lifetimes and native-contact survival" + category = "interactions" + required_files = {"trajectory", "topology"} + supported_systems = COMPLEX_SYSTEMS + order = 45 + default_params = {"cutoff": 4.5, "long_lived_ns": 1.0} + outputs = ["results/contact_lifetime.csv", "results/native_contact_survival.csv", + "figures/contact_lifetime_distribution.png", "figures/native_contact_survival.png"] + + def run(self, ctx) -> dict: + from moldynx.analysis.interface import InterfaceAnalysis + p = self.params(ctx) + plotting.set_style() + u = ctx.core_universe() + pa, pb = InterfaceAnalysis()._partners(ctx, u) + if pa is None: + return {"status": "skipped", "reason": "fewer than two partners"} + (nameA, A), (nameB, B) = pa, pb + hA, hB = A.select_atoms("not name H*"), B.select_atoms("not name H*") + locA = np.searchsorted(A.residues.resindices, hA.resindices) + locB = np.searchsorted(B.residues.resindices, hB.resindices) + rA, rB = A.residues.resids, B.residues.resids + + n = len(u.trajectory) + times, survival = np.empty(n), np.zeros(n) + active, lifetimes = {}, defaultdict(list) + formed, occ = defaultdict(int), defaultdict(int) + native = None + for i, ts in enumerate(ctx.iter_frames(u, desc="[contact_lifetime]")): + times[i] = ts.time / 1000.0 + pairs = capped_distance(hA.positions, hB.positions, max_cutoff=p["cutoff"], + box=ts.dimensions, return_distances=False) + cur = set(zip(locA[pairs[:, 0]].tolist(), locB[pairs[:, 1]].tolist())) \ + if len(pairs) else set() + if native is None: + native = set(cur) + survival[i] = len(native & cur) / len(native) if native else 0.0 + for q in cur: + if q in active: + active[q] += 1 + else: + active[q] = 1 + formed[q] += 1 + occ[q] += 1 + for q in [q for q in active if q not in cur]: + lifetimes[q].append(active.pop(q)) + for q, s in active.items(): + lifetimes[q].append(s) + + dt = float(times[1] - times[0]) if n > 1 else 0.0 + rows, all_lt = [], [] + for (a, b), lts in lifetimes.items(): + ns = np.asarray(lts) * dt + all_lt.extend(ns.tolist()) + rows.append({"a_resid": int(rA[a]), "b_resid": int(rB[b]), "occupancy": occ[(a, b)] / n, + "n_events": formed[(a, b)], "mean_lifetime_ns": float(ns.mean()), + "max_lifetime_ns": float(ns.max())}) + df = pd.DataFrame(rows, columns=["a_resid", "b_resid", "occupancy", "n_events", + "mean_lifetime_ns", "max_lifetime_ns"]) + df = df.sort_values("max_lifetime_ns", ascending=False) + ctx.write_csv(df, "contact_lifetime.csv") + ctx.write_csv(pd.DataFrame({"time_ns": times, "survival_fraction": survival}), + "native_contact_survival.csv") + all_lt = np.asarray(all_lt) + + fig, ax = plotting.new_axes() + if all_lt.size: + ax.hist(all_lt, bins=40, color=PALETTE["primary"], alpha=0.85) + ax.axvline(np.median(all_lt), ls="--", color=PALETTE["accent"], + label=f"median = {np.median(all_lt):.2f} ns") + ax.legend() + ax.set_yscale("log") + ax.set_xlabel("Contact lifetime (ns)") + ax.set_ylabel("Formation events") + ax.set_title("Inter-partner contact lifetimes") + plotting.save_figure(fig, ctx.fig_path("contact_lifetime_distribution"), dpi=ctx.config.dpi) + fig, ax = plotting.new_axes() + ax.plot(times, survival * 100, color=PALETTE["secondary"], lw=2) + ax.set_ylim(0, 105) + ax.set_xlabel("Time (ns)") + ax.set_ylabel("Surviving initial contacts (%)") + ax.set_title(f"Native contact survival ({nameA}–{nameB})") + plotting.save_figure(fig, ctx.fig_path("native_contact_survival"), dpi=ctx.config.dpi) + + return {"n_native_contacts": len(native or ()), + "final_survival_fraction": float(survival[-1]) if n else None, + "n_distinct_contacts": int(len(df)), + "median_lifetime_ns": float(np.median(all_lt)) if all_lt.size else None, + "max_lifetime_ns": float(all_lt.max()) if all_lt.size else None, + "n_long_lived_contacts": int((df.max_lifetime_ns >= p["long_lived_ns"]).sum()), + "figure": "native_contact_survival"} diff --git a/moldynx/analysis/porcupine.py b/moldynx/analysis/porcupine.py new file mode 100644 index 0000000..e702baf --- /dev/null +++ b/moldynx/analysis/porcupine.py @@ -0,0 +1,81 @@ +""" +Porcupine plot of the dominant motion (PC1 of the Cα covariance) + a mode PDB. + +Uses its own copy of the solute trajectory (alignment in memory would otherwise +rotate the shared trajectory, and its box, under the other analyses). PC1 comes +from an SVD of the aligned, centred coordinates -- no 3N × 3N covariance matrix, +so large proteins stay tractable. Arrows are coloured by chain. +""" + +from __future__ import annotations + +import numpy as np +import pandas as pd + +from moldynx import plotting +from moldynx.core.base import BaseAnalysis +from moldynx.plotting import PALETTE + +ARROW_SCALE = 2.5 + + +class Porcupine(BaseAnalysis): + name = "porcupine" + label = "Porcupine plot of PC1" + category = "dynamics" + required_files = {"trajectory", "topology"} + supported_systems = {"*"} + order = 70 + default_params = {"n_models": 20, "mode_scale": 2.0} + outputs = ["results/porcupine_pc1.csv", "results/porcupine_pc1.pdb", + "figures/porcupine_pc1.png"] + + def run(self, ctx) -> dict: + import MDAnalysis as mda + from MDAnalysis.analysis import align + p = self.params(ctx) + plotting.set_style() + ctx.core_universe() + u = mda.Universe(str(ctx.config.data_dir / "core.pdb"), str(ctx.config.data_dir / "core.xtc")) + ca = u.select_atoms("protein and name CA") + if ca.n_atoms < 3: + return {"status": "skipped", "reason": "fewer than three protein Cα atoms"} + align.AlignTraj(u, u, select="protein and name CA", ref_frame=0, in_memory=True).run() + X = np.array([ca.positions.copy() for _ in u.trajectory]) + mean = X.mean(0) + Xc = (X - mean).reshape(len(X), -1) + _, s, vt = np.linalg.svd(Xc, full_matrices=False) + var = s ** 2 / max(len(X) - 1, 1) + pc1 = vt[0].reshape(-1, 3) + amp = float(np.sqrt(var[0])) + disp = pc1 * amp + mag = np.linalg.norm(disp, axis=1) + + labels = np.array(["protein"] * ca.n_atoms, dtype=object) + for rec in ctx.core_meta.get("chains", []): + m = (ca.indices >= rec["core_start"]) & (ca.indices < rec["core_stop"]) + labels[m] = rec.get("segid") or f"chain {rec.get('index')}" + ctx.write_csv(pd.DataFrame({"resid": ca.resids, "chain": labels, "dx_A": disp[:, 0], + "dy_A": disp[:, 1], "dz_A": disp[:, 2], "magnitude_A": mag}), + "porcupine_pc1.csv") + fig = plotting.style.plt.figure(figsize=(7.4, 6.6)) + ax = fig.add_subplot(111, projection="3d") + ax.plot(*mean.T, color=PALETTE["muted"], lw=0.6, alpha=0.5) + colours = [PALETTE[c] for c in ("primary", "green", "secondary", "purple")] + for k, lab in enumerate(dict.fromkeys(labels)): + m = labels == lab + d = disp[m] * ARROW_SCALE + ax.quiver(*mean[m].T, *d.T, color=colours[k % 4], linewidth=1.0, label=lab) + ax.set_title(f"PC1 ({var[0] / var.sum() * 100:.1f}% of variance)") + ax.legend(loc="upper left") + plotting.save_figure(fig, ctx.fig_path("porcupine_pc1"), dpi=ctx.config.dpi) + with mda.Writer(str(ctx.csv_path("porcupine_pc1.pdb")), n_atoms=ca.n_atoms, + multiframe=True) as w: + for ph in np.linspace(0, 2 * np.pi, p["n_models"], endpoint=False): + ca.positions = (mean + p["mode_scale"] * amp * np.sin(ph) * pc1).astype(np.float32) + w.write(ca) + per_chain = {lab: float(mag[labels == lab].mean()) for lab in dict.fromkeys(labels)} + return {"pc1_variance_fraction": float(var[0] / var.sum()), "pc1_rms_amplitude_A": amp, + "max_displacement_resid": int(ca.resids[int(np.argmax(mag))]), + "max_displacement_A": float(mag.max()), "mean_displacement_by_chain_A": per_chain, + "figure": "porcupine_pc1"} diff --git a/moldynx/analysis/water_bridges.py b/moldynx/analysis/water_bridges.py new file mode 100644 index 0000000..087539e --- /dev/null +++ b/moldynx/analysis/water_bridges.py @@ -0,0 +1,76 @@ +""" +Water-mediated bridges between two partners (partner A – water – partner B). + +MDAnalysis ``WaterBridgeAnalysis`` (order 1: one bridging water) on the full +solvated trajectory (it needs the waters and the charges/bonds of the run input), +strided because the search is expensive. Partners are addressed by atom-index +ranges from the chain records of the full topology. +""" + +from __future__ import annotations + +import numpy as np +import pandas as pd + +from moldynx import plotting +from moldynx import statistics as st +from moldynx.core.base import BaseAnalysis +from moldynx.core.system import COMPLEX_SYSTEMS +from moldynx.plotting import PALETTE + + +class WaterBridges(BaseAnalysis): + name = "water_bridges" + label = "Water-mediated bridges between partners" + category = "interactions" + required_files = {"trajectory", "topology"} + supported_systems = COMPLEX_SYSTEMS + order = 60 + default_params = {"stride": 40, "order": 1, "persistent": 0.5} + outputs = ["results/water_bridges.csv", "tables/water_bridge_summary.csv", + "figures/water_bridges.png"] + + def run(self, ctx) -> dict: + p = self.params(ctx) + chains = sorted(getattr(ctx.system, "chains", []), key=lambda c: c.n_atoms, reverse=True)[:2] + water = ctx.system.components.get("water") if hasattr(ctx.system, "components") else None + if len(chains) < 2 or water is None or not water.n_residues: + return {"status": "skipped", + "reason": "needs two protein chains and explicit water in the run input"} + from MDAnalysis.analysis.hydrogenbonds import WaterBridgeAnalysis + plotting.set_style() + chains = sorted(chains, key=lambda c: c.atom_start) + u = ctx.full_universe() + sel = [f"index {c.atom_start}:{c.atom_stop - 1}" for c in chains] + wsel = "resname " + " ".join(water.resnames) + wb = WaterBridgeAnalysis(u, sel[0], sel[1], water_selection=wsel, order=p["order"], + update_selection=True) + wb.run(step=p["stride"], verbose=False) + cbt = wb.count_by_time() + t = np.array([x[0] for x in cbt], float) / 1000.0 + c = np.array([x[1] for x in cbt], float) + ctx.write_csv(pd.DataFrame({"time_ns": t, "n_water_bridges": c}), "water_bridges.csv") + rows = [] + try: + for entry in wb.count_by_type(): + rows.append({"bridge": " | ".join(str(x) for x in entry[:-1]), + "frequency": float(entry[-1])}) + except Exception: # older MDAnalysis: summary by type unavailable + pass + tab = pd.DataFrame(rows, columns=["bridge", "frequency"]).sort_values( + "frequency", ascending=False) + tab.to_csv(ctx.table_path("water_bridge_summary.csv"), index=False) + fig, ax = plotting.new_axes() + ax.plot(t, c, color=PALETTE["primary"], alpha=0.4, lw=1) + if len(c) > 2: + ax.plot(t, st.moving_average(c, max(2, len(c) // 10)), color=PALETTE["primary"], lw=2) + ax.set_xlabel("Time (ns)") + ax.set_ylabel("Water bridges") + ax.set_title(f"Water-mediated bridges (order {p['order']}, every {p['stride']} frames)") + plotting.save_figure(fig, ctx.fig_path("water_bridges"), dpi=ctx.config.dpi) + return {"partners": [c.segid for c in chains], "stride": p["stride"], + "frames": int(len(c)), "mean_bridges": float(c.mean()) if len(c) else 0.0, + "max_bridges": float(c.max()) if len(c) else 0.0, + "n_bridge_types": int(len(tab)), + "n_persistent": int((tab.frequency >= p["persistent"]).sum()), + "figure": "water_bridges"} diff --git a/moldynx/core/annotations.py b/moldynx/core/annotations.py new file mode 100644 index 0000000..8182189 --- /dev/null +++ b/moldynx/core/annotations.py @@ -0,0 +1,154 @@ +""" +User-supplied annotations: display names, biological numbering, domains by homology, +motifs, docking-site residues. Nothing here is invented -- every region comes from +the ``annotations:`` block of the run configuration. + +Schema (YAML):: + + annotations: + chains: + - {segid: seg_0_PROA, display: "α-zein Q946V6", role: ligand} + - {segid: seg_1_PROB, display: "ZmBiP2", role: receptor, numbering_offset: 213} + domains: # transferred by alignment, never by copying numbers + seg_1_PROB: + reference_name: "UniProt P11021 (human BiP)" + reference_sequence: "MKLSLVAAMLLLLSAARA..." # or reference_fasta: path + regions: {NBD: [26, 405], SBDbeta: [418, 507]} # reference numbering + motifs: + - {name: "Motif 1", sequence: "CSQAPIASLLPPYLSPAVSSVC", chain: seg_0_PROA} + docking_site: {seg_1_PROB: [405, 434, 435, 438]} # biological numbering + unresolved_metadata: {force_field_variant: null, salt_concentration_M: null} + +Biological numbering = trajectory resid − ``numbering_offset``; the default offset +makes each chain start at 1. +""" + +from __future__ import annotations + +from dataclasses import dataclass, field +from pathlib import Path + + +@dataclass +class ChainAnnotation: + segid: str + display: str + role: str | None + resid_first: int + resid_last: int + offset: int + sequence: str + + def bio(self, resid: int) -> int: + return int(resid) - self.offset + + def md(self, bio: int) -> int: + return int(bio) + self.offset + + @property + def bio_range(self) -> tuple[int, int]: + return self.bio(self.resid_first), self.bio(self.resid_last) + + +def chain_annotations(chains: list[dict], annotations: dict) -> dict[str, ChainAnnotation]: + """Merge persisted chain records with the user's chain annotations.""" + user = {c.get("segid"): c for c in (annotations or {}).get("chains", []) if c.get("segid")} + out = {} + for rec in chains: + seg = rec.get("segid") or f"chain{rec.get('index')}" + u = user.get(seg, {}) + out[seg] = ChainAnnotation( + segid=seg, display=u.get("display") or seg, role=u.get("role"), + resid_first=int(rec["resid_first"]), resid_last=int(rec["resid_last"]), + offset=int(u.get("numbering_offset", int(rec["resid_first"]) - 1)), + sequence=rec.get("sequence", "")) + return out + + +# --------------------------------------------------------------------------- # +# sequence alignment +# --------------------------------------------------------------------------- # +def _aligner(mode: str = "global"): + from Bio import Align + from Bio.Align import substitution_matrices + a = Align.PairwiseAligner() + a.mode = mode + a.substitution_matrix = substitution_matrices.load("BLOSUM62") + a.open_gap_score, a.extend_gap_score = -10.0, -0.5 + return a + + +def residue_map(query: str, reference: str) -> tuple[dict[int, int], float]: + """1-based reference position -> 1-based query position, and % identity over aligned pairs.""" + aln = _aligner("global").align(query, reference)[0] + ref_to_q: dict[int, int] = {} + same = pairs = 0 + for (qs, qe), (rs, re_) in zip(*aln.aligned): + for k in range(qe - qs): + ref_to_q[rs + k + 1] = qs + k + 1 + pairs += 1 + same += query[qs + k] == reference[rs + k] + return ref_to_q, (100.0 * same / pairs if pairs else 0.0) + + +def transfer_regions(query: str, reference: str, regions: dict[str, list[int]], + search: int = 15) -> list[dict]: + """ + Map region edges from reference to query numbering through a global alignment + (BLOSUM62, gap −10/−0.5). An edge that falls in an alignment gap is moved to + the nearest aligned reference position and flagged ``uncertain`` with its shift. + """ + ref_to_q, ident = residue_map(query, reference) + out = [] + for name, (start, end) in regions.items(): + edges, notes = [], [] + for pos, direction in ((int(start), 1), (int(end), -1)): + if pos in ref_to_q: + edges.append(ref_to_q[pos]) + continue + found = None + for d in range(1, search + 1): + for cand in (pos + direction * d, pos - direction * d): + if cand in ref_to_q: + found = (cand, d) + break + if found: + break + if found: + edges.append(ref_to_q[found[0]]) + notes.append(f"reference {pos} is in an alignment gap; used {found[0]} (±{found[1]})") + else: + edges.append(None) + notes.append(f"reference {pos} has no aligned residue within ±{search}") + out.append({"region": name, "reference": [int(start), int(end)], + "start": edges[0], "end": edges[1], "uncertain": bool(notes), + "notes": notes, "identity_pct": round(ident, 1)}) + return out + + +def locate_motif(sequence: str, motif: str) -> dict: + """Exact match first; else the best local alignment, with partial-match bookkeeping.""" + i = sequence.find(motif) + if i >= 0: + return {"start": i + 1, "end": i + len(motif), "match": "exact", "identity_pct": 100.0, + "covered": len(motif), "length": len(motif)} + try: + aln = _aligner("local").align(sequence, motif)[0] + except Exception: + return {"match": "none", "length": len(motif)} + blocks = list(zip(*aln.aligned)) + if not blocks: + return {"match": "none", "length": len(motif)} + s0, s1 = blocks[0][0][0], blocks[-1][0][1] + m0, m1 = blocks[0][1][0], blocks[-1][1][1] + same = sum(sequence[qs + k] == motif[ms + k] + for (qs, qe), (ms, _me) in blocks for k in range(qe - qs)) + covered = sum(qe - qs for (qs, qe), _ in blocks) + return {"start": s0 + 1, "end": s1, "match": "partial" if covered < len(motif) else "similar", + "identity_pct": round(100.0 * same / max(covered, 1), 1), + "motif_positions_covered": [m0 + 1, m1], "covered": covered, "length": len(motif)} + + +def read_fasta_sequence(path: str | Path) -> str: + lines = Path(path).read_text(encoding="utf-8").splitlines() + return "".join(x.strip() for x in lines if x and not x.startswith(">")) diff --git a/moldynx/statistics/__init__.py b/moldynx/statistics/__init__.py index 8996f91..4bb082d 100644 --- a/moldynx/statistics/__init__.py +++ b/moldynx/statistics/__init__.py @@ -6,3 +6,6 @@ plateau_detection, ) from moldynx.statistics.correlation import correlation_matrices # noqa: F401 +from moldynx.statistics.autocorr import ( # noqa: F401 + statistical_inefficiency, describe_correlated, detect_equilibration, drift, +) diff --git a/moldynx/statistics/autocorr.py b/moldynx/statistics/autocorr.py new file mode 100644 index 0000000..ecf7a22 --- /dev/null +++ b/moldynx/statistics/autocorr.py @@ -0,0 +1,73 @@ +""" +Autocorrelation-aware statistics for MD time series. + +* :func:`statistical_inefficiency` -- g = 1 + 2 Σ (1 − k/N) C(k), summed until the + autocorrelation first drops to zero (Chodera et al., J. Chem. Theory Comput. + 2007, 3, 26; the estimator used by pymbar's ``timeseries``). +* :func:`describe_correlated` -- mean, SD, naive SEM and the corrected SEM with the + number of effectively independent samples N/g. +* :func:`detect_equilibration` -- the start t₀ that maximises N_eff = (N − t₀)/g + (Chodera, J. Chem. Theory Comput. 2016, 12, 1799). +* :func:`drift` -- linear trend with p-value, and a half-vs-half Welch test. +""" + +from __future__ import annotations + +import numpy as np + + +def statistical_inefficiency(x) -> float: + x = np.asarray(x, float) + n = len(x) + if n < 4: + return 1.0 + d = x - x.mean() + var = float(np.mean(d * d)) + if var == 0: + return 1.0 + g = 1.0 + for k in range(1, n - 1): + c = float(np.sum(d[:n - k] * d[k:]) / ((n - k) * var)) + if c <= 0: + break + g += 2.0 * c * (1.0 - k / n) + return max(1.0, g) + + +def describe_correlated(x) -> dict: + x = np.asarray(x, float) + n = len(x) + if n == 0: + return {"n": 0} + sd = float(x.std(ddof=1)) if n > 1 else 0.0 + g = statistical_inefficiency(x) + return {"n": n, "mean": float(x.mean()), "sd": sd, + "sem_naive": sd / np.sqrt(n) if n else float("nan"), + "stat_ineff": g, "n_eff": n / g, "sem": sd / np.sqrt(n / g) if n else float("nan")} + + +def detect_equilibration(x, step: int = 1, max_fraction: float = 0.8) -> dict: + """Return ``{"t0": index, "g": g, "n_eff": N_eff}`` maximising N_eff = (N − t0)/g.""" + x = np.asarray(x, float) + best = {"t0": 0, "g": 1.0, "n_eff": 0.0} + for t0 in range(0, max(1, int(len(x) * max_fraction)), max(1, step)): + g = statistical_inefficiency(x[t0:]) + ne = (len(x) - t0) / g + if ne > best["n_eff"]: + best = {"t0": t0, "g": g, "n_eff": ne} + return best + + +def drift(t, x) -> dict: + """Linear slope (per unit of t) with p-value, and a half-vs-half Welch t-test.""" + from scipy import stats + t, x = np.asarray(t, float), np.asarray(x, float) + ok = np.isfinite(t) & np.isfinite(x) + t, x = t[ok], x[ok] + if len(x) < 6 or np.ptp(t) == 0: + return {"slope": None, "p_slope": None, "half_p": None} + r = stats.linregress(t, x) + h = len(x) // 2 + w = stats.ttest_ind(x[:h], x[h:], equal_var=False) + return {"slope": float(r.slope), "p_slope": float(r.pvalue), "half_p": float(w.pvalue), + "first_half_mean": float(x[:h].mean()), "second_half_mean": float(x[h:].mean())} diff --git a/tests/test_interface.py b/tests/test_interface.py index 48fa69c..fd44e76 100644 --- a/tests/test_interface.py +++ b/tests/test_interface.py @@ -36,9 +36,13 @@ def ctx(tmp_path): src.atoms[187:].translate([0.5 * k, 0, 0]) if k else None src.trajectory.ts.time = k * 100.0 w.write(src.atoms) + from moldynx.core.system import _ONE_LETTER + seq = "".join(_ONE_LETTER.get(r, "X") for r in src.residues.resnames) (cfg.data_dir / "core_meta.json").write_text(json.dumps({"chains": [ - {"segid": "seg_0_PROA", "index": 0, "core_start": 0, "core_stop": 187}, - {"segid": "seg_1_PROB", "index": 1, "core_start": 187, "core_stop": 850}]})) + {"segid": "seg_0_PROA", "index": 0, "core_start": 0, "core_stop": 187, + "resid_first": 1, "resid_last": 187, "sequence": seq[:187]}, + {"segid": "seg_1_PROB", "index": 1, "core_start": 187, "core_stop": 850, + "resid_first": 188, "resid_last": 850, "sequence": seq[187:]}]})) c = AnalysisContext.__new__(AnalysisContext) c.config, c.system, c._core = cfg, _System(), mda.Universe(str(core_pdb), str(core_xtc)) return c diff --git a/tests/test_regions_and_window.py b/tests/test_regions_and_window.py new file mode 100644 index 0000000..f8988b6 --- /dev/null +++ b/tests/test_regions_and_window.py @@ -0,0 +1,102 @@ +"""Phase 5b/6: contact lifetimes, porcupine, annotations, analysis window, autocorrelation stats.""" + +from __future__ import annotations + +import json + +import numpy as np +import pandas as pd +import pytest + +from test_interface import ctx # noqa: F401 (shared real-geometry complex fixture) +from moldynx.analysis.analysis_window import AnalysisWindow +from moldynx.analysis.annotation import Annotation +from moldynx.analysis.contact_lifetime import ContactLifetime +from moldynx.analysis.interface import InterfaceAnalysis +from moldynx.analysis.porcupine import Porcupine +from moldynx.core import annotations as ann +from moldynx.statistics import describe_correlated, detect_equilibration, statistical_inefficiency + + +# --------------------------------------------------------------------------- # +def test_statistical_inefficiency(): + rng = np.random.default_rng(1) + assert statistical_inefficiency(rng.normal(size=5000)) < 1.2 # white noise + x = np.zeros(20000) + for i in range(1, len(x)): # AR(1), phi = 0.9 + x[i] = 0.9 * x[i - 1] + rng.normal() + g = statistical_inefficiency(x) + assert 14 < g < 25 # theory: 19 + d = describe_correlated(x) + assert d["sem"] == pytest.approx(d["sem_naive"] * np.sqrt(g), rel=1e-6) + + +def test_detect_equilibration_skips_transient(): + rng = np.random.default_rng(2) + x = np.r_[np.linspace(10, 0, 200), rng.normal(0, 0.3, 800)] + assert 150 <= detect_equilibration(x, step=10)["t0"] <= 260 + + +# --------------------------------------------------------------------------- # +def test_region_transfer_and_motifs(): + pytest.importorskip("Bio") + ref = "MKVLAAGIVGLLLAQPAVSAQEKEVGTVIGIDLGTTYSCVGVFKNGRVEIIANDQGNRITPSYVAFTDGERLIGD" + query = ref[4:40] + "GG" + ref[40:] # N-terminal truncation + an insertion + mapped = ann.transfer_regions(query, ref, {"core": [26, 60]})[0] + assert mapped["start"] == 26 - 4 and not mapped["uncertain"] + assert mapped["end"] == 60 - 4 + 2 # shifted by the insertion + assert ann.locate_motif(query, ref[30:40])["match"] == "exact" + # a motif hanging off the N-terminus (only its last part is in the construct) + part = ann.locate_motif(query, "XXXXXX" + query[:6]) + assert part["match"] in ("partial", "similar") and part["start"] == 1 + + +def test_numbering_offset(): + c = ann.chain_annotations([{"segid": "B", "resid_first": 188, "resid_last": 850, + "sequence": "A"}], + {"chains": [{"segid": "B", "display": "Partner", "role": "receptor"}]}) + assert c["B"].bio(188) == 1 and c["B"].bio_range == (1, 663) and c["B"].md(205) == 392 + c = ann.chain_annotations([{"segid": "B", "resid_first": 214, "resid_last": 876, + "sequence": "A"}], {"chains": [{"segid": "B", + "numbering_offset": 213}]}) + assert c["B"].bio(418) == 205 + + +# --------------------------------------------------------------------------- # +def test_lifetime_porcupine_annotation_window(ctx): # noqa: F811 + InterfaceAnalysis().run(ctx) + life = ContactLifetime().run(ctx) + assert life["n_native_contacts"] > 0 and 0 <= life["final_survival_fraction"] <= 1 + por = Porcupine().run(ctx) + assert 0 < por["pc1_variance_fraction"] <= 1 + assert set(por["mean_displacement_by_chain_A"]) == {"seg_0_PROA", "seg_1_PROB"} + + seq_b = json.loads((ctx.config.data_dir / "core_meta.json").read_text())["chains"][1]["sequence"] + ctx.config.annotations = { + "chains": [{"segid": "seg_1_PROB", "display": "Partner B", "role": "receptor"}], + "domains": {"seg_1_PROB": {"reference_name": "self", "reference_sequence": seq_b, + "regions": {"first half": [1, 300]}}}, + "motifs": [{"name": "m1", "sequence": seq_b[99:111], "chain": "seg_1_PROB"}], + "docking_site": {"seg_1_PROB": [100, 101, 102]}} + out = Annotation().run(ctx) + assert out["chains"]["seg_1_PROB"]["biological_range"] == [1, 663] + assert out["motifs"][0]["start"] == 100 and out["motifs"][0]["match"] == "exact" + reg = pd.read_csv(ctx.csv_path("region_summary.csv")) + assert set(reg.kind) == {"domain", "motif", "docking site"} + + assert AnalysisWindow().run(ctx)["status"] == "skipped" # 5 frames: too short to judge + t = np.linspace(0, 100, 101) + rng = np.random.default_rng(3) + pd.DataFrame({"time_ns": t, + "n_inter_contacts": 60 + 0.4 * t + rng.normal(0, 2, 101), # still growing + "min_interface_dist_nm": 0.27 + rng.normal(0, 0.01, 101)}).to_csv( + ctx.csv_path("interface_timeseries.csv"), index=False) + win = AnalysisWindow().run(ctx) + assert win["binding_observables_stationary"] is False + assert "Residue–residue contacts" in win["non_stationary"] + assert win["decision_required"] is True and win["proposal"]["final_20pct"] == [80.0, 100.0] + + +def test_annotation_without_config_says_so(ctx): # noqa: F811 + out = Annotation().run(ctx) + assert out["annotated"] is False and "no annotations" in out["note"] From addddbf8cb4b91b2e9041cf6345f59792ea8be11 Mon Sep 17 00:00:00 2001 From: Hossam Mahmoud Date: Sat, 26 Sep 2026 06:08:03 +0300 Subject: [PATCH 2/6] MM-GBSA/MM-PBSA with gmx_MMPBSA: prepare, probe-gate, run, analyse Built on gmx_MMPBSA (https://github.com/Valdes-Tresanca-MS/gmx_MMPBSA; Valdes-Tresanca et al. 2021) and AmberTools MMPBSA.py; both cited in outputs. - binding/prepare.py: protein-only complex.tpr (convert-tpr on Protein), complex.top with solute molecule types only, complex.ndx from persisted chain identity (0 = receptor, 1 = ligand), no Amber forcefields line, ionic strength derived from ions and box volume, explicit print_res from residues ever within 6 A of the partner (fallback: all), run script with a probe gate (bonded deltas must be zero; decomposition coverage checked). - binding/analyse.py: autocorrelation-corrected means per window, drift, GB vs PB, per-residue hotspots in biological numbering, closure check, entropy validity gate (IE < 3.6, C2 < 6.0 kcal/mol). - analysis/mmpbsa.py: analyses existing outputs or prepares the package; never reports numbers it did not compute. New 'moldynx binding-energy'. - Reproduces a manual analysis of real gmx_MMPBSA outputs exactly. Co-Authored-By: Claude Opus 5.5 --- moldynx/analysis/mmpbsa.py | 263 +++++++++++++++++++++++++---------- moldynx/binding/__init__.py | 25 ++++ moldynx/binding/analyse.py | 158 +++++++++++++++++++++ moldynx/binding/prepare.py | 262 ++++++++++++++++++++++++++++++++++ moldynx/cli/main.py | 44 ++++++ tests/test_binding_energy.py | 65 +++++++++ 6 files changed, 741 insertions(+), 76 deletions(-) create mode 100644 moldynx/binding/__init__.py create mode 100644 moldynx/binding/analyse.py create mode 100644 moldynx/binding/prepare.py create mode 100644 tests/test_binding_energy.py diff --git a/moldynx/analysis/mmpbsa.py b/moldynx/analysis/mmpbsa.py index 86ef560..177f95a 100644 --- a/moldynx/analysis/mmpbsa.py +++ b/moldynx/analysis/mmpbsa.py @@ -1,98 +1,209 @@ """ -MM/PBSA workflow generation + structural energy-decomposition proxy. +MM-GBSA / MM-PBSA binding free energy with gmx_MMPBSA +(https://github.com/Valdes-Tresanca-MS/gmx_MMPBSA; see :mod:`moldynx.binding`). -MM-GBSA/MM-PBSA computes a *binding* free energy between two partners. For a -single ligand-free solute it is undefined, so this module GENERATES a ready-to-run -``gmx_MMPBSA`` workflow (input file + commands) for when two groups are defined -(a bound ligand, or two domains/chains), auto-detecting trajectory/topology/index. -It also provides a force-field-independent per-residue **structural contribution** -proxy from persistent contact occupancy (unitless, not kcal/mol). +* If gmx_MMPBSA outputs exist in ``/binding_energy/{gb,pb}`` they are + analysed: autocorrelation-corrected means per window, drift, GB vs PB, hotspots in + biological numbering, closure of the per-residue decomposition, entropy validity. +* Otherwise a complete, ready-to-run package is prepared (protein-only system, + derived ionic strength, explicit decomposition residues, probe-gated run script) + and the module reports ``prepared`` with the command -- never invented numbers. + ``moldynx binding-energy --execute`` runs it (Linux/WSL with gmx_MMPBSA). + +Results are approximate end-point estimates, not experimental affinities. """ from __future__ import annotations +import json +import re +from pathlib import Path + import numpy as np import pandas as pd -from moldynx.core.base import BaseAnalysis from moldynx import plotting +from moldynx.binding import CITATIONS, GMX_MMPBSA_URL +from moldynx.core import annotations as ann +from moldynx.core.base import BaseAnalysis +from moldynx.core.system import COMPLEX_SYSTEMS from moldynx.plotting import PALETTE -MMPBSA_IN = """\ -&general - sys_name = "moldynx_system", - startframe = 1, endframe = 999999, interval = 10, verbose = 2, - forcefields = "oldff/leaprc.ff99SB, leaprc.gaff", -/ -&gb - igb = 5, saltcon = 0.150, -/ -&decomp - idecomp = 1, dec_verbose = 0, print_res = "within 6", -/ -""" + +def _windows(ctx, t_end: float) -> dict: + spec = (ctx.config.binding_energy or {}).get("primary_window") + win = {f"0-{t_end:g} ns": (0.0, t_end)} + if isinstance(spec, (list, tuple)) and len(spec) == 2: + win[f"{spec[0]:g}-{spec[1]:g} ns"] = (float(spec[0]), float(spec[1])) + elif isinstance(spec, str) and re.match(r"^\s*[\d.]+\s*-\s*[\d.]+", spec): + a, b = (float(x) for x in re.findall(r"[\d.]+", spec)[:2]) + win[f"{a:g}-{b:g} ns"] = (a, b) + else: + win[f"{0.8 * t_end:g}-{t_end:g} ns"] = (0.8 * t_end, t_end) + return win class MMPBSA(BaseAnalysis): name = "mmpbsa" - label = "MM/PBSA workflow + structural proxy" + label = "MM-GBSA / MM-PBSA binding free energy (gmx_MMPBSA)" category = "interactions" required_files = {"trajectory", "topology"} - supported_systems = {"*"} - order = 160 # after contact_map for the structural proxy - outputs = ["tables/mmpbsa.in", "tables/mmpbsa_run.sh", "report/mmpbsa_workflow.md"] - - def _run_script(self, ctx) -> str: - traj = (ctx.fileset.trajectory or ctx.config.trajectory) - top = (ctx.fileset.topology or ctx.config.topology) - gtop = ctx.fileset.gmx_top - return f"""#!/usr/bin/env bash -# gmx_MMPBSA workflow (run where GROMACS + AmberTools + gmx_MMPBSA are installed). -# MM/PBSA needs TWO groups -- edit the make_ndx selections to define them. -set -euo pipefail -TPR="{top.name if top else 'topol.tpr'}" -XTC="{traj.name if traj else 'traj.xtc'}" -TOP="{gtop.name if gtop else 'topol.top'}" -gmx make_ndx -f "$TPR" -o index.ndx <<'EOF' -r 1-100 -name 20 GroupA -r 101-9999 -name 21 GroupB -q -EOF -mpirun -np 4 gmx_MMPBSA -O -i mmpbsa.in -cs "$TPR" -ct "$XTC" -ci index.ndx \\ - -cg 20 21 -cp "$TOP" -o FINAL_RESULTS_MMPBSA.dat -eo FINAL_RESULTS_MMPBSA.csv \\ - -do FINAL_DECOMP_MMPBSA.dat -deo FINAL_DECOMP_MMPBSA.csv -""" + supported_systems = COMPLEX_SYSTEMS + order = 170 # after interface (residue set) and analysis_window + default_params = {"target_spacing_ns": 1.0, "ever_cutoff_column": "within_6A_ever"} + outputs = ["binding_energy/", "results/binding_energy_summary.json", + "figures/binding_energy_timeseries.png", "figures/binding_energy_hotspots.png"] def run(self, ctx) -> dict: - plotting.set_style() - ctx.table_path("mmpbsa.in").write_text(MMPBSA_IN, encoding="utf-8") - ctx.table_path("mmpbsa_run.sh").write_text(self._run_script(ctx), encoding="utf-8") - has_ligand = ctx.system.flags.get("has_ligand", False) - meaningful = has_ligand or ctx.system.n_protein_chains >= 2 - (ctx.config.report_dir / "mmpbsa_workflow.md").write_text( - f"# MM/PBSA workflow\n\nBinding-ΔG is " - f"{'applicable (define receptor/ligand groups)' if meaningful else 'undefined for this single ligand-free solute'}. " - f"Generated `tables/mmpbsa.in` and `tables/mmpbsa_run.sh` " - f"(edit the group selections). Requires GROMACS + AmberTools + gmx_MMPBSA.\n", - encoding="utf-8") + be_dir = ctx.config.output_dir / "binding_energy" + if (be_dir / "gb" / "FINAL_RESULTS_MMGBSA.csv").exists(): + return self.analyse_outputs(ctx, be_dir) + return self.prepare(ctx, be_dir) + + # ------------------------------------------------------------------ # + def _partners(self, ctx): + chains = sorted(ctx.core_meta.get("chains", []), key=lambda c: c["core_stop"] - c["core_start"], + reverse=True)[:2] + if len(chains) < 2: + return None + roles = {c.get("segid"): c.get("role") for c in (ctx.config.annotations or {}).get("chains", [])} + rec = next((c for c in chains if roles.get(c["segid"]) == "receptor"), chains[0]) + lig = next(c for c in chains if c is not rec) + order = sorted(chains, key=lambda c: c["core_start"]) + return rec, lig, {c["segid"]: "ABCDEFGH"[k] for k, c in enumerate(order)} - summary = {"workflow_generated": True, "binding_dg_meaningful": bool(meaningful)} - occ_csv = ctx.csv_path("contact_occupancy.csv") - if occ_csv.exists(): - occ = pd.read_csv(occ_csv, index_col=0) - resids = occ.columns.astype(int).to_numpy() - contribution = occ.to_numpy().sum(axis=1) - ctx.write_csv(pd.DataFrame({"resid": resids, - "structural_contribution": contribution}), - "residue_structural_contribution.csv") - fig, ax = plotting.new_axes(figsize=(8.2, 4.4)) - ax.bar(resids, contribution, color=PALETTE["accent"], width=1.0) - ax.set_xlabel("Residue"); ax.set_ylabel("Structural contribution") - ax.set_title("Per-residue structural contribution (MM/PBSA proxy)") - plotting.save_figure(fig, ctx.fig_path("residue_structural_contribution"), - dpi=ctx.config.dpi) - summary["figure"] = "residue_structural_contribution" - summary["top_residues"] = resids[np.argsort(contribution)[::-1][:5]].tolist() + def prepare(self, ctx, be_dir: Path) -> dict: + from moldynx.binding import prepare as prep + from moldynx.io import gromacs + fs = ctx.fileset + missing = [k for k in ("gmx_top", "topology") if getattr(fs, k, None) is None] + if missing or Path(fs.topology).suffix.lower() != ".tpr": + return {"status": "skipped", + "reason": f"binding energy needs the run input (.tpr) and topol.top + toppar/ " + f"(missing: {', '.join(missing) or 'a .tpr run input'})"} + partners = self._partners(ctx) + if partners is None: + return {"status": "skipped", "reason": "fewer than two protein chains"} + rec, lig, letters = partners + u = ctx.core_universe() + n = len(u.trajectory) + dt = float(u.trajectory[1].time - u.trajectory[0].time) / 1000.0 if n > 1 else 1.0 + interval = max(1, int(round(self.params(ctx)["target_spacing_ns"] / dt))) if dt else 1 + u.trajectory[0] + volume = abs(float(np.linalg.det(u.trajectory.ts.triclinic_dimensions))) / 1000.0 + ionic = prep.derive_ionic_strength(getattr(ctx.system, "ion_counts", {}) or {}, volume) + if not ionic["ions"]: + ionic = {**ionic, "ionic_strength_M": 0.0, + "note": "no ions found in the run input: ionic strength 0 (set " + "binding_energy.ionic_strength to override)"} + override = (ctx.config.binding_energy or {}).get("ionic_strength") + if isinstance(override, (int, float)): + ionic = {**ionic, "ionic_strength_M": float(override), + "note": f"user-specified ({override} M); derived value kept for reference"} + temp = 300.0 + if fs.log is not None: + ref = gromacs.parse_log(fs.log).ref_t + temp = ref[0] if ref else temp + local = {"R": [], "L": []} + res_csv = ctx.csv_path("interface_residues.csv") + if res_csv.exists(): + df = pd.read_csv(res_csv) + col = next((c for c in df.columns if c.startswith("within_") and c.endswith("_ever")), None) + if col: + for side, c in (("R", rec), ("L", lig)): + sel = df[(df.partner == c["segid"]) & df[col].astype(bool)] + local[side] = sel.partner_residue.astype(int).tolist() + solute_xtc = ctx.config.data_dir / "core.xtc" + if u.atoms.n_atoms != sum(c["core_stop"] - c["core_start"] for c in (rec, lig)): + return {"status": "skipped", + "reason": "the solute trajectory holds more than the two partners; " + "binding energy for mixed solutes is not automated yet"} + top_text = Path(fs.gmx_top).read_text(encoding="utf-8", errors="replace") + g = lambda c: {"name": c.get("segid"), "letter": letters[c["segid"]], # noqa: E731 + "atom_start1": c["core_start"] + 1, "atom_stop1": c["core_stop"]} + info = prep.write_package( + be_dir, tpr_src=Path(fs.topology), top_text=top_text, + toppar_src=Path(fs.toppar) if fs.toppar else None, solute_xtc=solute_xtc, + receptor=g(rec), ligand=g(lig), temperature=temp, ionic=ionic, + frames={"start": 1, "end": n, "interval": interval, "frame_dt_ns": dt}, + print_res_local=local, meta={"citations": CITATIONS}) + gmx = gromacs.find_gmx() + tpr_ok, why = (prep.make_complex_tpr(gmx, Path(fs.topology), be_dir) if gmx + else (False, "GROMACS not found: run make_ndx/convert-tpr as described")) + return {"status": "prepared", "directory": str(be_dir), "complex_tpr": tpr_ok, + "complex_tpr_note": why, "receptor": rec["segid"], "ligand": lig["segid"], + "ionic_strength_M": ionic["ionic_strength_M"], "temperature_K": temp, + "frames": info["frames"], "print_res": info["print_res"], + "how_to_run": f"on Linux/WSL with gmx_MMPBSA ({GMX_MMPBSA_URL}): " + f"cd {be_dir} && bash run_mmpbsa.sh (or: moldynx binding-energy " + f"--run {ctx.config.output_dir} --execute)"} + + # ------------------------------------------------------------------ # + def analyse_outputs(self, ctx, be_dir: Path) -> dict: + from moldynx.binding.analyse import analyse + meta_p = be_dir / "prepare_meta.json" + meta = json.loads(meta_p.read_text(encoding="utf-8")) if meta_p.exists() else {} + dt = meta.get("frames", {}).get("frame_dt_ns") or (ctx.config.binding_energy or {}).get( + "frame_dt_ns", 0.1) + chains = ann.chain_annotations(ctx.core_meta.get("chains", []), ctx.config.annotations or {}) + rec = chains.get(meta.get("receptor", {}).get("name")) + lig = chains.get(meta.get("ligand", {}).get("name")) + mapper = lambda c: (lambda k: (c.display, c.bio(c.resid_first + k - 1))) if c else None # noqa: E731 + gb = be_dir / "gb" + pb_csv = be_dir / "pb" / "FINAL_RESULTS_MMPBSA.csv" + dec = gb / "FINAL_DECOMP_MMGBSA.csv" + from moldynx.binding.analyse import read_delta + t_end = float(read_delta(gb / "FINAL_RESULTS_MMGBSA.csv", dt).time_ns.max()) + dat = gb / "FINAL_RESULTS_MMGBSA.dat" + res = analyse(gb / "FINAL_RESULTS_MMGBSA.csv", pb_csv if pb_csv.exists() else None, + frame_dt_ns=dt, windows=_windows(ctx, t_end), + decomp_csv=dec if dec.exists() else None, + gb_dat=dat if dat.exists() else None, + receptor=mapper(rec), ligand=mapper(lig)) + res["components"].to_csv(ctx.csv_path("binding_energy_components.csv"), index=False) + frames = res["gb"][["frame", "time_ns", "TOTAL"]].rename(columns={"TOTAL": "dG_GB"}) + if res["pb"] is not None: + frames["dG_PB"] = res["pb"].TOTAL.values + frames.to_csv(ctx.csv_path("binding_energy_per_frame.csv"), index=False) + for w, df in res.get("decomposition", {}).items(): + df.to_csv(ctx.csv_path(f"binding_energy_residues_{w.replace(' ', '')}.csv"), index=False) + self._figures(ctx, res) + summary = {k: res[k] for k in ("checks", "windows", "headline", "trend_per_ns", "gb_vs_pb")} + summary.update({"entropy": res.get("entropy"), "closure": res.get("closure"), + "residues_in_decomposition": res.get("residues_in_decomposition"), + "hotspots": {w: df.head(20)[["partner", "residue", "resname", "TOTAL"]] + .to_dict("records") + for w, df in res.get("decomposition", {}).items()}, + "method": {"tool": GMX_MMPBSA_URL, "citations": CITATIONS, + "prepared": meta}, + "caveat": "approximate end-point estimates, not experimental affinities"}) + ctx.csv_path("binding_energy_summary.json").write_text( + json.dumps(summary, indent=2, default=float), encoding="utf-8") + summary["status"] = "analysed" + summary["figure"] = "binding_energy_timeseries" return summary + + def _figures(self, ctx, res) -> None: + plotting.set_style() + plt = plotting.style.plt + gb, pb = res["gb"], res["pb"] + fig, ax = plotting.new_axes(figsize=(8, 4.4)) + ax.plot(gb.time_ns, gb.TOTAL, color=PALETTE["primary"], lw=1.2, label="MM-GBSA") + if pb is not None: + ax.plot(pb.time_ns, pb.TOTAL, color=PALETTE["secondary"], lw=1.2, label="MM-PBSA") + ax.set_xlabel("Time (ns)") + ax.set_ylabel("ΔG_bind (kcal/mol)") + ax.set_title("Binding free energy over time (end-point estimate)") + ax.legend() + plotting.save_figure(fig, ctx.fig_path("binding_energy_timeseries"), dpi=ctx.config.dpi) + dec = res.get("decomposition", {}) + if dec: + w = list(dec)[-1] + top = dec[w].head(20).iloc[::-1] + fig, ax = plt.subplots(figsize=(7, 6)) + lab = [f"{p or s} {r}{n}" for p, s, r, n in zip(top.partner, top.side, top.resname, + top.residue)] + ax.barh(lab, top.TOTAL, color=PALETTE["primary"]) + ax.set_xlabel("Per-residue ΔG (kcal/mol)") + ax.set_title(f"Top 20 residues, {w}") + plotting.save_figure(fig, ctx.fig_path("binding_energy_hotspots"), dpi=ctx.config.dpi) diff --git a/moldynx/binding/__init__.py b/moldynx/binding/__init__.py new file mode 100644 index 0000000..1ae0ab2 --- /dev/null +++ b/moldynx/binding/__init__.py @@ -0,0 +1,25 @@ +""" +Binding free energy by MM-GBSA / MM-PBSA with gmx_MMPBSA. + +MolDynX prepares, gates, runs and analyses calculations performed with +**gmx_MMPBSA** (https://github.com/Valdes-Tresanca-MS/gmx_MMPBSA; Valdés-Tresanca +et al., J. Chem. Theory Comput. 2021, 17, 6281–6291), which drives AmberTools' +MMPBSA.py (Miller et al., J. Chem. Theory Comput. 2012, 8, 3314–3321). Cite both +when you use these results. + +* :mod:`moldynx.binding.prepare` -- protein-only complex system, input files, + derived ionic strength, explicit decomposition residue set, run script with a + probe gate. +* :mod:`moldynx.binding.analyse` -- parse the outputs, autocorrelation-corrected + statistics per window, drift, GB vs PB, hotspots, closure and entropy gates. +""" + +GMX_MMPBSA_URL = "https://github.com/Valdes-Tresanca-MS/gmx_MMPBSA" +CITATIONS = [ + "Valdés-Tresanca, M. S.; Soler, M. A.; Moreno, E. et al. gmx_MMPBSA: A New Tool to " + "Perform End-State Free Energy Calculations with GROMACS. J. Chem. Theory Comput. 2021, " + "17, 6281–6291. https://doi.org/10.1021/acs.jctc.1c00645", + "Miller, B. R.; McGee, T. D.; Swails, J. M. et al. MMPBSA.py: An Efficient Program for " + "End-State Free Energy Calculations. J. Chem. Theory Comput. 2012, 8, 3314–3321. " + "https://doi.org/10.1021/ct300418h", +] diff --git a/moldynx/binding/analyse.py b/moldynx/binding/analyse.py new file mode 100644 index 0000000..cc65c8f --- /dev/null +++ b/moldynx/binding/analyse.py @@ -0,0 +1,158 @@ +""" +Analysis of gmx_MMPBSA outputs (``FINAL_RESULTS_*.csv``, ``FINAL_DECOMP_MMGBSA.csv``, ``*.dat``). + +All means carry autocorrelation-corrected standard errors (statistical inefficiency, +N_eff = N/g). Values are approximate end-point estimates, **not** experimental +affinities; entropy terms are reported only if σ(interaction energy) passes the +validity limits gmx_MMPBSA itself warns about (IE < 3.6, C2 < 6.0 kcal/mol). +""" + +from __future__ import annotations + +import io +import re +from pathlib import Path + +import numpy as np +import pandas as pd + +from moldynx.statistics import describe_correlated, drift + +BONDED = ["BOND", "ANGLE", "DIHED", "UB", "IMP", "CMAP", "1-4 VDW", "1-4 EEL"] +COMPONENTS = [("ΔE_vdW", "VDWAALS", "VDWAALS"), ("ΔE_elec", "EEL", "EEL"), + ("ΔG_polar", "EGB", "EPB"), ("ΔG_nonpolar", "ESURF", "ENPOLAR"), + ("ΔG_gas", "GGAS", "GGAS"), ("ΔG_solv", "GSOLV", "GSOLV"), + ("ΔG_bind", "TOTAL", "TOTAL")] +ENTROPY_LIMITS = {"IE": 3.6, "C2": 6.0} # kcal/mol, as gmx_MMPBSA warns + + +def read_delta(csv: str | Path, frame_dt_ns: float) -> pd.DataFrame: + """The ``Delta Energy Terms`` block, with time from the original frame numbers.""" + lines = Path(csv).read_text(encoding="utf-8").splitlines() + i = next(k for k, line in enumerate(lines) if line.strip().startswith("Delta Energy Terms")) + rows = [] + for line in lines[i + 1:]: + if line.startswith("Frame") or (line[:1].isdigit()): + rows.append(line) + elif rows: + break + df = pd.read_csv(io.StringIO("\n".join(rows))).rename(columns={"Frame #": "frame"}) + df["time_ns"] = (df.frame - 1) * frame_dt_ns + return df + + +def read_decomp(csv: str | Path, frame_dt_ns: float, section: str = "DELTAS", + kind: str = "Total Decomposition Contribution") -> pd.DataFrame: + """Per-residue decomposition rows (``L::THR:1`` / ``R::GLU:205`` labels split out).""" + lines = Path(csv).read_text(encoding="utf-8").splitlines() + i = next(k for k, line in enumerate(lines) if line.strip().startswith(section)) + j = next(k for k in range(i, len(lines)) if lines[k].strip().startswith(kind)) + rows = [lines[j + 1]] + for line in lines[j + 2:]: + if not line or not line[:1].isdigit(): + break + rows.append(line) + df = pd.read_csv(io.StringIO("\n".join(rows))).rename(columns={"Frame #": "frame", + "Residue": "label"}) + m = df.label.str.extract(r"^([LR])::([A-Z0-9]+):(-?\d+)$") + df["side"], df["resname"], df["resnum"] = m[0], m[1], m[2].astype(int) + df["time_ns"] = (df.frame - 1) * frame_dt_ns + return df + + +def read_entropy(dat: str | Path) -> dict: + """σ(interaction energy) and −TΔS for IE and C2 from a ``FINAL_RESULTS_*.dat``.""" + text = Path(dat).read_text(encoding="utf-8", errors="replace") + out = {} + for kind in ("IE", "C2"): + m = re.search(rf"^\S+\s+{kind}\s+([-\d.]+)\s+([-\d.]+)\s+([-\d.]+)\s+([-\d.]+)", text, re.M) + if m: + sigma, value = float(m.group(1)), float(m.group(2)) + out[kind] = {"sigma_int_kcal": sigma, "minus_TdS_kcal": value, + "sem_kcal": float(m.group(4)), "limit_kcal": ENTROPY_LIMITS[kind], + "valid": sigma < ENTROPY_LIMITS[kind]} + return out + + +def _window_mask(t: np.ndarray, window) -> np.ndarray: + lo, hi = window + return (t >= lo - 1e-9) & (t <= hi + 1e-9) + + +def analyse(gb_csv, pb_csv=None, frame_dt_ns: float = 0.1, windows: dict | None = None, + decomp_csv=None, gb_dat=None, receptor=None, ligand=None) -> dict: + """ + Parameters + ---------- + windows : {"0-100 ns": (0, 100), ...}; defaults to the full run and its final 20 %. + receptor / ligand : callables mapping a gmx_MMPBSA residue number (numbered from 1 + within each partner) to a (display, biological number) tuple; identity if None. + """ + gb = read_delta(gb_csv, frame_dt_ns) + pb = read_delta(pb_csv, frame_dt_ns) if pb_csv else None + t = gb.time_ns.to_numpy() + if windows is None: + end = float(t[-1]) + windows = {f"{t[0]:g}-{end:g} ns": (float(t[0]), end), + f"{0.8 * end:g}-{end:g} ns": (0.8 * end, end)} + checks = {"frames": int(len(gb)), + "max_abs_bonded_delta": float(max(np.abs(gb[c]).max() for c in BONDED if c in gb)), + "same_frames_gb_pb": bool(pb is None or (gb.frame.values == pb.frame.values).all())} + if pb is not None: + checks["max_abs_ggas_gb_minus_pb"] = float(np.abs(gb.GGAS - pb.GGAS).max()) + checks["single_trajectory_consistent"] = checks["max_abs_bonded_delta"] < 1e-6 + + rows = [] + for method, df in (("GB", gb), ("PB", pb)): + if df is None: + continue + for wname, w in windows.items(): + mask = _window_mask(df.time_ns.to_numpy(), w) + for label, cg, cp in COMPONENTS: + col = cg if method == "GB" else cp + if col not in df: + continue + d = describe_correlated(df[col].to_numpy()[mask]) + rows.append({"method": method, "window": wname, "term": label, **d}) + comp = pd.DataFrame(rows) + headline = comp[comp.term == "ΔG_bind"][["method", "window", "mean", "sem", "sd", "n", + "n_eff", "stat_ineff"]].reset_index(drop=True) + trend = {m: drift(df.time_ns.to_numpy(), df.TOTAL.to_numpy()) + for m, df in (("GB", gb), ("PB", pb)) if df is not None} + gbpb = None + if pb is not None: + r = float(np.corrcoef(gb.TOTAL, pb.TOTAL)[0, 1]) + gbpb = {"pearson_r": r, "mean_offset_pb_minus_gb": float((pb.TOTAL - gb.TOTAL).mean())} + + out = {"checks": checks, "windows": {k: list(v) for k, v in windows.items()}, + "headline": headline.to_dict("records"), "components": comp, + "trend_per_ns": {m: {"slope": d["slope"], "p": d["p_slope"], "half_p": d["half_p"]} + for m, d in trend.items()}, + "gb_vs_pb": gbpb, "gb": gb, "pb": pb} + if gb_dat: + out["entropy"] = read_entropy(gb_dat) + if decomp_csv: + out.update(_decomposition(decomp_csv, frame_dt_ns, windows, gb, receptor, ligand)) + return out + + +def _decomposition(decomp_csv, frame_dt_ns, windows, gb, receptor, ligand) -> dict: + dec = read_decomp(decomp_csv, frame_dt_ns) + ident = lambda n: (None, n) # noqa: E731 + per_window, closure = {}, {} + for wname, w in windows.items(): + sub = dec[_window_mask(dec.time_ns.to_numpy(), w)] + res = sub.groupby(["side", "resname", "resnum"]).TOTAL.mean().reset_index() + mapped = [(receptor or ident)(n) if s == "R" else (ligand or ident)(n) + for s, n in zip(res.side, res.resnum)] + res["partner"] = [m[0] for m in mapped] + res["residue"] = [m[1] for m in mapped] + res = res.sort_values("TOTAL").reset_index(drop=True) + per_window[wname] = res + total = gb.TOTAL[_window_mask(gb.time_ns.to_numpy(), w)].mean() + s = float(res.TOTAL.sum()) + closure[wname] = {"sum_residues": s, "total": float(total), + "unattributed": float(total - s), + "closes": bool(abs(total - s) <= 1.0)} + return {"decomposition": per_window, "closure": closure, + "residues_in_decomposition": int(dec.drop_duplicates("label").shape[0])} diff --git a/moldynx/binding/prepare.py b/moldynx/binding/prepare.py new file mode 100644 index 0000000..16a298c --- /dev/null +++ b/moldynx/binding/prepare.py @@ -0,0 +1,262 @@ +""" +Prepare a gmx_MMPBSA calculation (https://github.com/Valdes-Tresanca-MS/gmx_MMPBSA). + +* **Protein-only system.** Implicit solvent needs no explicit water, and a solvated + system of 1–2 M atoms does not fit a typical WSL/laptop memory. ``complex.tpr`` is + cut from the canonical run input with ``gmx convert-tpr`` (Protein group), + ``complex.top`` keeps only the solute molecule types, ``complex.ndx`` has exactly + two groups (0 = receptor, 1 = ligand) built from the persisted chain identity, and + the trajectory is the validated, PBC-treated solute trajectory. No Amber + ``forcefields=`` line: parameters come from the GROMACS topology. +* **Ionic strength derived** from the ions in the run input and the box volume + (never an unrecorded 0.150 M). +* **Explicit decomposition residues**: every residue that came within 6 Å of the + partner in *any* frame (``print_res = "within 6"`` is evaluated on the first frame + only and silently drops residues that bind later). +* **Run script with a probe gate**: ~11 frames of GB first; production runs only if + every bonded Δ term is exactly zero (single-trajectory consistency) and the + decomposition covers the requested residues (else it falls back to all residues). +""" + +from __future__ import annotations + +import json +import os +import re +import shutil +from pathlib import Path + +AVOGADRO = 6.02214076e23 +CHARGE = {"POT": 1, "K": 1, "K+": 1, "NA": 1, "SOD": 1, "NA+": 1, "LI": 1, "RB": 1, "CS": 1, + "CLA": -1, "CL": -1, "CL-": -1, "BR": -1, "IOD": -1, "I": -1, "F": -1, + "MG": 2, "MG2": 2, "CAL": 2, "CA": 2, "CA2": 2, "ZN": 2, "ZN2": 2, "MN": 2, + "FE2": 2, "FE3": 3, "CU": 2, "CO": 2, "NI": 2, "CD": 2} +SOLVENT = {"SOL", "WAT", "HOH", "TIP3", "TIP4", "TIP5", "TIP3P", "SPC", "SPCE", "T3P", "T4P", + "OPC", "TP3"} + + +def derive_ionic_strength(ion_counts: dict[str, int], volume_nm3: float) -> dict: + """I = ½ Σ cᵢ zᵢ² with cᵢ = nᵢ / (N_A V). Also the 1:1 salt concentration (minority ion).""" + litres = volume_nm3 * 1e-24 + terms, pos, neg = [], 0, 0 + for name, n in ion_counts.items(): + z = CHARGE.get(name.upper()) + if z is None: + continue + c = n / (AVOGADRO * litres) + terms.append((name, n, z, c)) + pos += n if z > 0 else 0 + neg += n if z < 0 else 0 + ionic = 0.5 * sum(c * z * z for _, _, z, c in terms) + salt = min(pos, neg) / (AVOGADRO * litres) if terms else 0.0 + return {"ionic_strength_M": round(ionic, 4), "salt_1to1_M": round(salt, 4), + "volume_nm3": volume_nm3, "ions": {n: k for n, k, _, _ in terms}, + "note": "counts from the run input, volume of the first solute-trajectory frame; " + "the box includes the solute, so this slightly underestimates the " + "bulk concentration"} + + +def complex_topology(top_text: str, keep: list[str]) -> str: + """Replace the ``[ molecules ]`` section with the kept molecule types (count 1 each, in order).""" + m = re.search(r"^\s*\[\s*molecules\s*\]\s*$", top_text, re.M | re.I) + if not m: + raise ValueError("topology has no [ molecules ] section") + head = top_text[:m.end()] + body = "\n; Compound\t#mols\n" + "".join(f"{k}\t1\n" for k in keep) + return head + body + + +def solute_molecules(top_text: str) -> list[str]: + """Molecule types of ``[ molecules ]`` that are not solvent or ions, in order.""" + m = re.search(r"^\s*\[\s*molecules\s*\]\s*$(.*)", top_text, re.M | re.I | re.S) + out = [] + for line in (m.group(1) if m else "").splitlines(): + line = line.split(";")[0].strip() + if not line or line.startswith("["): + continue + name, *count = line.split() + if name.upper() in SOLVENT or name.upper() in CHARGE: + continue + out.extend([name] * int(count[0] if count else 1)) + return out + + +def ranges(numbers: list[int]) -> str: + """[1,2,3,7,9,10] -> '1-3,7,9-10'.""" + nums = sorted(set(int(n) for n in numbers)) + out, start = [], None + for i, n in enumerate(nums): + if start is None: + start = n + if i + 1 == len(nums) or nums[i + 1] != n + 1: + out.append(f"{start}" if start == n else f"{start}-{n}") + start = None + return ",".join(out) + + +def input_files(temperature: float, ionic: float, start: int, end: int, interval: int, + print_res: str, sys_name: str = "complex") -> tuple[str, str]: + """(mmgbsa.in, mmpbsa.in) -- GB and PB run separately because PBRadii is global.""" + gb = f"""MM-GBSA generated by MolDynX Tools (gmx_MMPBSA) +&general + sys_name = "{sys_name}_GB", + startframe = {start}, + endframe = {end}, + interval = {interval}, + temperature = {temperature}, + PBRadii = 3, + interaction_entropy = 1, + ie_segment = 25, + c2_entropy = 1, + verbose = 2, +/ +&gb + igb = 5, + saltcon = {ionic:.4f}, +/ +&decomp + idecomp = 2, + dec_verbose = 3, + print_res = "{print_res}", +/ +""" + pb = f"""MM-PBSA generated by MolDynX Tools (gmx_MMPBSA) +&general + sys_name = "{sys_name}_PB", + startframe = {start}, + endframe = {end}, + interval = {interval}, + temperature = {temperature}, + PBRadii = 7, + verbose = 2, +/ +&pb + istrng = {ionic:.4f}, + fillratio = 4.0, + inp = 1, +/ +""" + return gb, pb + + +RUN_SCRIPT = r"""#!/usr/bin/env bash +# gmx_MMPBSA (https://github.com/Valdes-Tresanca-MS/gmx_MMPBSA) -- generated by MolDynX Tools. +# probe gate (GB, ~11 frames) -> GB with decomposition -> PB. Run from this folder. +set -uo pipefail +ENV_NAME="${MOLDYNX_MMPBSA_ENV:-gmxMMPBSA}" +source "$(conda info --base 2>/dev/null)/etc/profile.d/conda.sh" 2>/dev/null && conda activate "$ENV_NAME" +export OMP_NUM_THREADS=1 OMPI_MCA_btl_vader_single_copy_mechanism=none +CORES=$(( $(nproc) / __SMT__ )); [ "$CORES" -lt 1 ] && CORES=1 +FREE_GB=$(awk '/MemAvailable/{printf "%d", $2/1048576}' /proc/meminfo) +PB_RANKS=$(( FREE_GB * 10 / 35 )); [ "$PB_RANKS" -gt "$CORES" ] && PB_RANKS=$CORES; [ "$PB_RANKS" -lt 1 ] && PB_RANKS=1 +LOG=run.log; stamp(){ echo "[$(date '+%F %T')] $*" | tee -a "$LOG"; } +COMMON="-cs complex.tpr -ci complex.ndx -cg 0 1 -ct complex.xtc -cp complex.top" + +stamp "probe: GB on 11 frames" +mkdir -p probe && sed -E 's/^( *endframe *= *).*/\1__PROBE_END__,/' mmgbsa.in > probe/probe.in +( cd probe && ln -sf ../complex.* ../toppar . && gmx_MMPBSA -O -nogui -i probe.in $COMMON \ + -o R.dat -eo R.csv -do D.dat -deo D.csv > run.log 2>&1 ); rc=$? +python3 check_probe.py probe/R.csv probe/D.csv requested_residues.json | tee -a "$LOG"; gate=$? +if [ $rc -ne 0 ] || [ $gate -eq 3 ]; then stamp "PROBE FAILED (exit $rc, gate $gate) -- not running production"; exit 3; fi +if [ $gate -eq 2 ]; then stamp "decomposition did not cover the requested residues -> print_res = all" + sed -i -E 's/^( *print_res *= *).*/\1"all",/' mmgbsa.in; fi + +stamp "GB: $CORES MPI ranks" +mkdir -p gb && ( cd gb && ln -sf ../complex.* ../toppar . && cp ../mmgbsa.in . && \ + mpirun -np $CORES --bind-to none gmx_MMPBSA -O -nogui -i mmgbsa.in $COMMON \ + -o FINAL_RESULTS_MMGBSA.dat -eo FINAL_RESULTS_MMGBSA.csv \ + -do FINAL_DECOMP_MMGBSA.dat -deo FINAL_DECOMP_MMGBSA.csv > run.log 2>&1 ); stamp "GB exit $?" +stamp "PB: $PB_RANKS MPI ranks (~3.3 GB per rank)" +mkdir -p pb && ( cd pb && ln -sf ../complex.* ../toppar . && cp ../mmpbsa.in . && \ + mpirun -np $PB_RANKS --bind-to none gmx_MMPBSA -O -nogui -i mmpbsa.in $COMMON \ + -o FINAL_RESULTS_MMPBSA.dat -eo FINAL_RESULTS_MMPBSA.csv > run.log 2>&1 ); stamp "PB exit $?" +""" + +CHECK_PROBE = r'''"""Probe gate: exit 3 = inconsistent (bonded delta != 0 or no total), 2 = decomposition incomplete.""" +import csv, io, json, sys +res_csv, dec_csv, req_json = sys.argv[1:4] +lines = open(res_csv, encoding="utf-8").read().splitlines() +i = next(k for k, l in enumerate(lines) if l.startswith("Delta Energy Terms")) +rows = list(csv.DictReader(io.StringIO("\n".join(l for l in lines[i + 1:] if l)))) +bonded = ["BOND", "ANGLE", "DIHED", "UB", "IMP", "CMAP", "1-4 VDW", "1-4 EEL"] +bad = [f"{b}={r[b]}" for r in rows for b in bonded if b in r and abs(float(r[b])) > 1e-6] +if bad or not rows or "TOTAL" not in rows[0]: + print("probe: bonded delta terms are not zero:", bad[:5]); sys.exit(3) +print(f"probe: {len(rows)} frames, bonded deltas all zero, dTOTAL {rows[0]['TOTAL']}") +want = json.load(open(req_json)) +seen = set() +for l in open(dec_csv, encoding="utf-8"): + p = l.split(",") + if len(p) > 2 and p[1][:3] in ("R::", "L::"): + seen.add(p[1].split("::")[0] + ":" + p[1].rsplit(":", 1)[1]) +missing = [w for w in want if w not in seen] +print(f"probe: decomposition covers {len(want) - len(missing)} of {len(want)} requested residues") +sys.exit(2 if missing else 0) +''' + + +def write_package(out_dir: Path, *, tpr_src: Path, top_text: str, toppar_src: Path | None, + solute_xtc: Path, receptor: dict, ligand: dict, temperature: float, + ionic: dict, frames: dict, print_res_local: dict, meta: dict) -> dict: + """ + Write everything except ``complex.tpr`` (made with GROMACS by :func:`make_complex_tpr`). + + receptor / ligand : {"name", "letter", "atom_start1", "atom_stop1"} -- 1-based, inclusive + atom ranges in the solute trajectory; ``letter`` is the chain letter of the partner + in topology order (gmx_MMPBSA ``print_res`` syntax). + print_res_local : {"R": [local residue numbers], "L": [...]} numbered from 1 per partner. + """ + out_dir.mkdir(parents=True, exist_ok=True) + keep = solute_molecules(top_text) + (out_dir / "complex.top").write_text(complex_topology(top_text, keep), encoding="utf-8") + if toppar_src and toppar_src.is_dir(): + shutil.copytree(toppar_src, out_dir / "toppar", dirs_exist_ok=True) + with open(out_dir / "complex.ndx", "w", newline="\n") as fh: + for g in (receptor, ligand): + ids = list(range(g["atom_start1"], g["atom_stop1"] + 1)) + fh.write(f"[ {re.sub(r'[^A-Za-z0-9_]', '_', g['name'])} ]\n") + for k in range(0, len(ids), 15): + fh.write(" ".join(f"{i:6d}" for i in ids[k:k + 15]) + "\n") + shutil.copyfile(solute_xtc, out_dir / "complex.xtc") + parts = [] + for side, g in (("R", receptor), ("L", ligand)): + if print_res_local.get(side): + parts.append(f"{g['letter']}/{ranges(print_res_local[side])}") + # no interface residue set available -> decompose everything ("within 6" would be + # evaluated on the first frame only and miss residues that bind later) + print_res = " ".join(parts) if parts else "all" + gb, pb = input_files(temperature, ionic["ionic_strength_M"], frames["start"], frames["end"], + frames["interval"], print_res) + (out_dir / "mmgbsa.in").write_text(gb, encoding="utf-8") + (out_dir / "mmpbsa.in").write_text(pb, encoding="utf-8") + probe_end = frames["start"] + 10 * frames["interval"] + smt = 2 if (os.cpu_count() or 2) > 1 else 1 + (out_dir / "run_mmpbsa.sh").write_text( + RUN_SCRIPT.replace("__PROBE_END__", str(probe_end)).replace("__SMT__", str(smt)), + encoding="utf-8", newline="\n") + (out_dir / "check_probe.py").write_text(CHECK_PROBE, encoding="utf-8", newline="\n") + requested = [f"{s}:{n}" for s, ns in print_res_local.items() for n in ns] + (out_dir / "requested_residues.json").write_text(json.dumps(requested), encoding="utf-8") + info = {**meta, "receptor": receptor, "ligand": ligand, "solute_molecules": keep, + "temperature_K": temperature, "ionic": ionic, "frames": frames, + "print_res": print_res, "tpr_source": str(tpr_src), + "gmx_mmpbsa": "https://github.com/Valdes-Tresanca-MS/gmx_MMPBSA"} + (out_dir / "prepare_meta.json").write_text(json.dumps(info, indent=2, default=str), + encoding="utf-8") + return info + + +def make_complex_tpr(gmx, tpr_src: Path, out_dir: Path) -> tuple[bool, str]: + """``gmx make_ndx`` (default groups) + ``gmx convert-tpr`` on the Protein group.""" + from moldynx.io.gromacs import run_gmx + ndx = out_dir / "default.ndx" + r = run_gmx(gmx, ["make_ndx", "-f", str(tpr_src), "-o", str(ndx)], stdin="q\n") + if r.returncode != 0 or not ndx.exists(): + return False, "gmx make_ndx failed: " + (r.stderr or "")[-300:] + names = re.findall(r"^\[\s*(.+?)\s*\]", ndx.read_text(), re.M) + if "Protein" not in names: + return False, "no 'Protein' group in the default index" + r = run_gmx(gmx, ["convert-tpr", "-s", str(tpr_src), "-n", str(ndx), + "-o", str(out_dir / "complex.tpr")], stdin=f"{names.index('Protein')}\n") + ok = r.returncode == 0 and (out_dir / "complex.tpr").exists() + return ok, "" if ok else "gmx convert-tpr failed: " + (r.stderr or "")[-300:] diff --git a/moldynx/cli/main.py b/moldynx/cli/main.py index 5f0423f..5824d90 100644 --- a/moldynx/cli/main.py +++ b/moldynx/cli/main.py @@ -70,6 +70,12 @@ def build_parser() -> argparse.ArgumentParser: it.add_argument("--allow-ambiguous", dest="allow_ambiguous", action="store_true") it.add_argument("--include-dir", dest="include_dir", action="append") + be = sub.add_parser("binding-energy", + help="MM-GBSA/PBSA with gmx_MMPBSA: prepare, optionally run, analyse.") + _add_common(be) + be.add_argument("--execute", action="store_true", + help="Run the prepared gmx_MMPBSA script now (hours; Linux or WSL).") + d = sub.add_parser("detect", help="Detect and print the system composition.") _add_common(d) @@ -111,6 +117,43 @@ def _yaml_keys(path) -> set: return set((yaml.safe_load(Path(path).read_text()) or {}).keys()) +def _cmd_binding_energy(args) -> int: + """Prepare (or analyse) the gmx_MMPBSA calculation; --execute runs it in between.""" + import json + import subprocess + from moldynx.analysis.mmpbsa import MMPBSA + from moldynx.core.config import RunConfig + from moldynx.core.context import AnalysisContext + from moldynx.core.pipeline import build_plan + from moldynx.io.gromacs import windows_to_wsl + cfg = RunConfig.from_args(args, yaml_path=args.config) + if cfg.input_dir is None: + print("error: --input (or 'input_dir' in --config) is required.") + return 2 + if not args.output and "output_dir" not in _yaml_keys(args.config): + cfg.output_dir = default_output_dir(cfg.input_dir) + fs, val, system, _sel, _skip = build_plan(cfg) + if not val.ok: + print(val.report()) + return 2 + ctx = AnalysisContext(cfg, system, fs) + out = MMPBSA().run(ctx) + print(json.dumps({k: v for k, v in out.items() if k not in ("headline",)}, indent=1, + default=str)[:3000]) + if args.execute and out.get("status") == "prepared": + d = out["directory"] + cmd = (["wsl.exe", "--", "bash", "-c", f"cd '{windows_to_wsl(d)}' && bash run_mmpbsa.sh"] + if sys.platform.startswith("win") else ["bash", "run_mmpbsa.sh"]) + print(f"[binding-energy] running gmx_MMPBSA in {d} (this takes hours) ...") + rc = subprocess.call(cmd, cwd=None if sys.platform.startswith("win") else d) + if rc != 0: + print(f"[binding-energy] run script exited with {rc}; see {d}/run.log") + return rc + out = MMPBSA().run(ctx) + print(json.dumps(out.get("headline"), indent=1, default=float)) + return 0 + + def _cmd_intake(args) -> int: from moldynx.io.intake import run_intake, write_intake if not args.input: @@ -191,6 +234,7 @@ def main(argv=None) -> int: print(f"MolDynX Tools (moldynx) {__version__}") return 0 return {"analyze": _cmd_analyze, "detect": _cmd_detect, "intake": _cmd_intake, + "binding-energy": _cmd_binding_energy, "list-analyses": _cmd_list}[args.command](args) diff --git a/tests/test_binding_energy.py b/tests/test_binding_energy.py new file mode 100644 index 0000000..54b0c84 --- /dev/null +++ b/tests/test_binding_energy.py @@ -0,0 +1,65 @@ +"""MM-GBSA/PBSA: analysis of real gmx_MMPBSA outputs (reference values from a manual analysis).""" + +from __future__ import annotations + +import pytest + +from conftest import FIX +from moldynx.binding import analyse as be +from moldynx.binding import prepare as prep + +MM = FIX / "gmx_mmpbsa" + + +@pytest.fixture(scope="module") +def result(): + return be.analyse(MM / "gb" / "FINAL_RESULTS_MMGBSA.csv", MM / "pb" / "FINAL_RESULTS_MMPBSA.csv", + frame_dt_ns=0.1, windows={"0-100": (0, 100), "80-100": (80, 100)}, + decomp_csv=MM / "gb" / "FINAL_DECOMP_MMGBSA_frames1-3.csv", + gb_dat=MM / "gb" / "FINAL_RESULTS_MMGBSA.dat") + + +def _hl(res, method, window): + return next(h for h in res["headline"] if h["method"] == method and h["window"] == window) + + +@pytest.mark.parametrize("method, window, mean, sem", [ + ("GB", "0-100", -33.4, 2.1), ("GB", "80-100", -37.8, 2.4), + ("PB", "0-100", -40.7, 2.5), ("PB", "80-100", -47.1, 3.3), +]) +def test_reproduces_manual_analysis(result, method, window, mean, sem): + h = _hl(result, method, window) + assert h["mean"] == pytest.approx(mean, abs=0.05) + assert h["sem"] == pytest.approx(sem, abs=0.05) # autocorrelation-corrected + + +def test_checks_gb_pb_and_entropy_gate(result): + c = result["checks"] + assert c["frames"] == 101 and c["single_trajectory_consistent"] and c["same_frames_gb_pb"] + assert c["max_abs_ggas_gb_minus_pb"] < 0.02 + assert result["gb_vs_pb"]["pearson_r"] == pytest.approx(0.73, abs=0.01) + assert result["gb_vs_pb"]["mean_offset_pb_minus_gb"] == pytest.approx(-7.4, abs=0.05) + ie, c2 = result["entropy"]["IE"], result["entropy"]["C2"] + assert ie["sigma_int_kcal"] == pytest.approx(86.72) and ie["valid"] is False + assert c2["valid"] is False and c2["limit_kcal"] == 6.0 + + +def test_decomposition_parsing(result): + assert result["residues_in_decomposition"] > 50 + first = result["decomposition"]["0-100"] + assert set(first.side) == {"R", "L"} and first.TOTAL.iloc[0] < 0 + + +def test_prepare_helpers(): + ions = prep.derive_ionic_strength({"POT": 1498, "CLA": 1480}, 15605.2) + assert ions["ionic_strength_M"] == pytest.approx(0.1594, abs=1e-3) + assert ions["salt_1to1_M"] == pytest.approx(0.1575, abs=1e-3) + top = ('#include "toppar/forcefield.itp"\n[ system ]\nx\n\n[ molecules ]\n; Compound #mols\n' + "PROA 1\nPROB 1\nTIP3 518060\nPOT 1498\nCLA 1480\n") + assert prep.solute_molecules(top) == ["PROA", "PROB"] + out = prep.complex_topology(top, ["PROA", "PROB"]) + assert "TIP3" not in out and out.rstrip().endswith("PROB\t1") and "#include" in out + assert prep.ranges([5, 1, 2, 3, 9, 10]) == "1-3,5,9-10" + gb, pb = prep.input_files(303.15, 0.1594, 1, 1001, 10, "A/1-5 B/10-12") + assert 'print_res = "A/1-5 B/10-12"' in gb and "forcefields" not in gb + assert "PBRadii = 7" in pb and "istrng = 0.1594" in pb From 2ab2c8080ba2bb2f8661725354075397e9751960 Mon Sep 17 00:00:00 2001 From: Hossam Mahmoud Date: Sat, 26 Sep 2026 06:11:57 +0300 Subject: [PATCH 3/6] Data-driven documents and proper self-contained HTML - report/documents.py: PBC_VALIDATION, EQUILIBRATION and BINDING_ENERGY written from the results files (Observations / Interpretation / Limitations), wording rules enforced (no stability claim from an RMSD plateau, binding energies are end-point estimates, entropy only if valid, gmx_MMPBSA credited and cited). - Markdown rendered with python-markdown to self-contained HTML (contents, embedded images); replaces the renderer that only deleted '**'. - generator: links the companion documents, reports PBC proof, interface, preparation, stationarity and binding-energy results. Co-Authored-By: Claude Opus 5.5 --- moldynx/report/documents.py | 316 ++++++++++++++++++++++++++++++++++++ moldynx/report/generator.py | 90 +++++----- pyproject.toml | 5 +- tests/test_documents.py | 67 ++++++++ 4 files changed, 432 insertions(+), 46 deletions(-) create mode 100644 moldynx/report/documents.py create mode 100644 tests/test_documents.py diff --git a/moldynx/report/documents.py b/moldynx/report/documents.py new file mode 100644 index 0000000..8ac66f8 --- /dev/null +++ b/moldynx/report/documents.py @@ -0,0 +1,316 @@ +""" +Data-driven documents: PBC_VALIDATION, EQUILIBRATION, BINDING_ENERGY (+ HTML). + +Every number comes from a results file of the same run; no system names or +conclusions live in code. Each document separates Observations, Interpretation +and Limitations, and the wording rules are fixed here: no "stable" from an RMSD +plateau, binding energies are end-point estimates (not affinities), entropy is +quoted only when valid, unknown metadata is stated as not recorded. +""" + +from __future__ import annotations + +import base64 +import html as _html +import json +import re +from pathlib import Path + +CSS = ("body{font-family:Arial,Helvetica,sans-serif;max-width:980px;margin:2rem auto;padding:0 1rem;" + "color:#1d1d1f;line-height:1.55}h1{border-bottom:3px solid #1f6f8b;padding-bottom:.3rem}" + "h2{color:#1f6f8b;margin-top:2rem}table{border-collapse:collapse;margin:1rem 0;font-size:.92em}" + "th,td{border:1px solid #ccc;padding:5px 9px;text-align:left;vertical-align:top}" + "th{background:#f0f4f6}code{background:#f4f4f4;padding:0 3px;border-radius:3px}" + "img{max-width:100%;border:1px solid #eee}blockquote{border-left:4px solid #e07a5f;" + "margin:1rem 0;padding:.4rem 1rem;background:#fdf3ef}nav#toc{background:#f7f9fa;" + "padding:.6rem 1.2rem;border:1px solid #e3e8eb}nav#toc a{text-decoration:none}") + + +# --------------------------------------------------------------------------- # +# Markdown -> self-contained HTML +# --------------------------------------------------------------------------- # +def to_html(md_text: str, base_dir: Path, title: str) -> str: + """Render Markdown (tables, fenced code, sub) with a contents list; embed local images.""" + try: + import markdown + body = markdown.markdown(md_text, extensions=["tables", "fenced_code", "attr_list", + "sane_lists", "md_in_html"]) + except ImportError: # plain fallback keeps the document readable + body = "
" + _html.escape(md_text) + "
" + toc, n = [], 0 + + def anchor(m): + nonlocal n + n += 1 + level, text = m.group(1), m.group(2) + slug = f"s{n}" + if level == "2": + toc.append(f'
  • {re.sub("<[^>]+>", "", text)}
    • ') + return f'{text}' + body = re.sub(r"(.*?)", anchor, body) + + def embed(m): + src = m.group(1) + p = (base_dir / src).resolve() + if p.exists() and p.suffix.lower() in (".png", ".jpg", ".jpeg", ".svg"): + mime = "image/svg+xml" if p.suffix.lower() == ".svg" else f"image/{p.suffix[1:].lower()}" + return f'src="data:{mime};base64,{base64.b64encode(p.read_bytes()).decode()}"' + return m.group(0) + body = re.sub(r'src="([^"]+)"', embed, body) + nav = f'' if toc else "" + return (f"{_html.escape(title)}" + f"{nav}{body}") + + +def write_doc(path: Path, md_text: str, formats=("md", "html")) -> list[Path]: + path.parent.mkdir(parents=True, exist_ok=True) + out = [] + if "md" in formats: + path.write_text(md_text, encoding="utf-8") + out.append(path) + if "html" in formats: + title = (re.search(r"^# (.+)$", md_text, re.M) or [None, path.stem])[1] + h = path.with_suffix(".html") + h.write_text(to_html(md_text, path.parent, title), encoding="utf-8") + out.append(h) + return out + + +def _load(p: Path) -> dict | None: + try: + return json.loads(p.read_text(encoding="utf-8")) + except (OSError, ValueError): + return None + + +def _f(x, nd=2, unit=""): + try: + return f"{float(x):.{nd}f}{unit}" + except (TypeError, ValueError): + return "—" + + +def _pm(d: dict, key_mean="mean", key_err="sem", nd=1): + return f"{_f(d.get(key_mean), nd)} ± {_f(d.get(key_err), nd)}" + + +# --------------------------------------------------------------------------- # +# PBC +# --------------------------------------------------------------------------- # +def pbc_document(s: dict, fig_rel: str = "../figures/pbc_validation.png") -> str: + c = s["checks"] + ok = all(c.values()) + L = ["# Periodic-boundary validation", "", + f"Treatment: **{s['mode']}** · molecules made whole: {'yes' if s['made_whole'] else '**no**'}" + f" · frames: {s['n_frames']}", "", + "## Observations", "", + "| Molecule | Atoms | Frames split in the raw data | Whole-box translations undone | " + "Extent mean / max (nm) |", "|---|---|---|---|---|"] + for u in s["units"]: + L.append(f"| {u['label']} | {u['n_atoms']:,} | {u['frames_split_in_raw']} | " + f"{u['whole_box_translations_undone']} | {_f(u['extent_whole_nm_mean'])} / " + f"{_f(u['extent_whole_nm_max'])} |") + L.append("") + for p in s.get("pairs", []): + r = p["min_dist_raw_pbc_nm"] + L.append(f"- Minimum heavy-atom distance {p['pair'][0]}–{p['pair'][1]} (PBC-aware, raw): " + f"{_f(r[0], 3)}–{_f(r[2], 3)} nm (mean {_f(r[1], 3)}); frames without a " + f"heavy-atom contact: {p['frames_without_heavy_atom_contact']}.") + if s.get("half_box_min_nm"): + L.append(f"- Half of the smallest box edge (minimum-image limit): {_f(s['half_box_min_nm'])} nm.") + L += ["", "## Proof that processing changed only wholeness / continuity", "", + "| Check | Result |", "|---|---|"] + names = {"only_whole_box_translations": "every atom moved only by whole box vectors", + "molecules_whole": "no bond longer than 2.5 Å after processing", + "multiple_translations_only_in_split_frames": + "two translation vectors only in frames where the molecule was split", + "interchain_distance_preserved": "partner distances equal their PBC-aware raw values"} + L += [f"| {names.get(k, k)} | {'pass' if v else '**FAIL**'} |" for k, v in c.items()] + L += ["", f"Largest deviation from a pure box translation: " + f"{_f(s.get('max_dev_from_box_translation_A'), 5)} Å.", "", + f"![PBC evidence]({fig_rel})", "", "## Interpretation", "", + ("All checks pass: the processed solute trajectory differs from the raw one only by " + "whole-box translations that make molecules whole and continuous; distance-based " + "analyses are unaffected by periodic images." if ok else + "At least one check failed. Analyses that depend on distances between molecules " + "must be read with the failed check in mind (see the table above)."), "", + "## Limitations", "", + "- The diagnosis and proof cover the analysed frames only (frame slice of this run).", + "- A centre-of-mass distance is not a proximity measure for extended molecules; " + "contact is judged from minimum heavy-atom distances."] + return "\n".join(L) + "\n" + + +# --------------------------------------------------------------------------- # +# Equilibration +# --------------------------------------------------------------------------- # +LABEL = {"em": "Energy minimisation", "nvt": "NVT", "npt": "NPT", "production": "Production"} + + +def equilibration_document(s: dict, fig_rel: str = "../figures/equilibration_overview.png") -> str: + st = s["stages"] + L = ["# Preparation and equilibration", "", + "Reconstructed from the stage logs, energy files, run inputs and job scripts found at " + "intake. Missing inputs are listed, never invented.", "", + "## Stages", "", "| Stage | Integrator / dt | Length | Thermostat | Barostat | Outcome |", + "|---|---|---|---|---|---|"] + for k in ("em", "nvt", "npt", "production"): + r = st.get(k) + if not r: + continue + em = r.get("minimization") + if em: + outcome = (f"{em.get('outcome')} in {em.get('steps')} steps; Fmax " + f"{_f(em.get('fmax'), 1)} kJ mol⁻¹ nm⁻¹" + + (" — **requested tolerance not reached**" + if em.get("reached_emtol") is False else "")) + length, thermo, baro = f"{r.get('nsteps')} steps max", "—", "—" + else: + sim = r.get("simulated_ps") + length = (f"{sim / 1000:g} ns" if sim and sim >= 1000 else f"{sim:g} ps") if sim else "—" + ref = r.get("ref_t_K") or [] + thermo = f"{r.get('tcoupl')} {ref[0]:g} K" if ref else str(r.get("tcoupl")) + baro = str(r.get("pcoupl")) + (f", τp {r.get('tau_p_ps')} ps" + if r.get("pcoupl") not in (None, "No") else "") + c = r.get("counts", {}) + outcome = (f"{r.get('sessions')} session(s); LINCS warnings " + f"{c.get('lincs_warnings', 0)}, fatal errors {c.get('fatal_errors', 0)}") + L.append(f"| {LABEL[k]} | {r.get('integrator')} / {r.get('dt_ps') or '—'} ps | {length} | " + f"{thermo} | {baro} | {outcome} |") + L += ["", "## Observations", ""] + loc = st.get("em", {}).get("fmax_location") + if loc: + L.append(f"- Largest residual force after minimisation: atom {loc['atom_number']} " + f"({loc['atom']} of {loc['resname']}, chain {loc.get('chain', loc.get('segid'))}, " + f"chain residue {loc.get('chain_residue', loc.get('resid'))}).") + nvt = st.get("nvt", {}) + if nvt.get("temperature_settled_ps") is not None: + a = nvt.get("temperature_after_settling", {}) + L.append(f"- NVT temperature within ±2 K of target from {_f(nvt['temperature_settled_ps'], 0)} ps; " + f"{_f(a.get('mean'), 2)} ± {_f(a.get('sd'), 2)} K afterwards.") + npt = st.get("npt", {}) + e = npt.get("energy", {}).get("terms", {}) + if e: + parts = [] + for term, unit, nd in (("Temperature", "K", 3), ("Pressure", "bar", 2), + ("Density", "kg m⁻³", 2)): + if term in e: + parts.append(f"{term.lower()} {_f(e[term]['tail_mean'], nd)} ± " + f"{_f(e[term]['tail_sd'], nd)} {unit}") + tail = npt["energy"].get("tail_ps") + L.append(f"- NPT, last {_f(tail, 0)} ps: " + "; ".join(parts) + ".") + if npt.get("box_x_nm"): + b = npt["box_x_nm"] + v = npt.get("volume_nm3") or [None, None, None] + L.append(f"- Box {_f(b[0], 3)} → {_f(b[1], 3)} nm (volume change {_f(v[2], 1)} %); density " + f"within 0.1 % of its final value after {_f(npt.get('density_reached_99_9pct_ps'), 0)} ps.") + pr = s.get("position_restraints", {}) + if pr: + L += ["", "| Run input | Restrained atoms | By force constant (kJ mol⁻¹ nm⁻²) | Per molecule type |", + "|---|---|---|---|"] + for k, r in pr.items(): + if r.get("available"): + fc = ", ".join(f"{n} @ {f}" for f, n in r["by_force_constant"].items()) or "none" + L.append(f"| {LABEL.get(k, k)} | {r['n_restrained']:,} | {fc} | " + f"{', '.join(str(x) for x in r['by_molecule_block']) or '—'} |") + else: + L.append(f"| {LABEL.get(k, k)} | not audited | {r.get('error', '')} | — |") + prot = s.get("protonation", {}).get("non_default", []) + if prot: + L += ["", "Non-default protonation states: " + "; ".join( + f"{p['resname']} {p.get('chain_residue', p['resid'])} ({p.get('chain', '')}, {p['state']})" + for p in prot) + "."] + tl = [r for r in s.get("timeline", []) if r.get("gap_after_previous_s") is not None] + if tl: + L.append("Gaps between consecutive stages: " + ", ".join( + f"{LABEL[r['stage']]} started {_f(r['gap_after_previous_s'], 0)} s after the previous " + "stage ended" for r in tl) + ".") + L += ["", f"![Equilibration overview]({fig_rel})", "", "## Interpretation", ""] + L += [f"- {v[0].upper() + v[1:]}." for v in s.get("verdicts", [])] + L += ["", "## What could not be recovered", ""] + L += [f"- {m}" for m in s.get("missing", [])] or ["- Nothing: every expected stage file was found."] + L += ["", "## Limitations", "", + "- A Berendsen barostat relaxes the box but does not sample the NPT ensemble; " + "only production statistics should be interpreted thermodynamically.", + "- A residual drift at the end of equilibration means production started from a system " + "that was still relaxing; production convergence is assessed separately."] + return "\n".join(L) + "\n" + + +# --------------------------------------------------------------------------- # +# Binding energy +# --------------------------------------------------------------------------- # +def binding_energy_document(s: dict, window_sel: dict | None = None) -> str: + L = ["# Binding free energy (MM-GBSA / MM-PBSA)", "", + "> **End-point estimates, not experimental affinities.** Values are mean ± " + "autocorrelation-corrected standard error (kcal/mol) with the number of effectively " + "independent frames.", "", + f"Computed with [gmx_MMPBSA]({s['method']['tool']}) (AmberTools MMPBSA.py), prepared " + "and analysed by MolDynX Tools.", "", "## Results", "", + "| Method | Window | ΔGbind | N (effective) |", "|---|---|---|---|"] + for h in s["headline"]: + L.append(f"| {h['method']} | {h['window']} | {_pm(h)} | {h['n']} ({_f(h['n_eff'], 0)}) |") + tr = s.get("trend_per_ns", {}) + drifting = [m for m, d in tr.items() if d.get("p") is not None and d["p"] < 0.01] + if drifting: + L += ["", "> **Not stationary:** ΔGbind drifts over the run (" + + "; ".join(f"{m} {_f(tr[m]['slope'] * 10, 2)} kcal/mol per 10 ns, p = " + f"{tr[m]['p']:.1g}" for m in drifting) + + "). Each window describes a different state; no single value is an equilibrium " + "average."] + if window_sel and window_sel.get("summary", {}).get("decision_required"): + L += ["", "> The averaging window was **not** chosen by the user " + "(`binding_energy.primary_window`); the full run and its final 20 % are shown."] + g = s.get("gb_vs_pb") + if g: + L += ["", f"GB and PB agree frame by frame with Pearson r = {_f(g['pearson_r'], 2)}; " + f"PB − GB = {_f(g['mean_offset_pb_minus_gb'], 1)} kcal/mol on average."] + ent = s.get("entropy") or {} + if ent: + L += ["", "## Entropy", "", "| Method | σ(interaction energy) | Validity limit | Status |", + "|---|---|---|---|"] + for k, e in ent.items(): + L.append(f"| {k} | {_f(e['sigma_int_kcal'], 1)} | < {e['limit_kcal']} | " + f"{'valid' if e['valid'] else '**invalid — not reported**'} |") + clo = s.get("closure") or {} + if clo: + L += ["", "## Per-residue decomposition", "", + f"{s.get('residues_in_decomposition')} residues decomposed. Closure (Σ residues vs " + "ΔGbind, GB):", "", "| Window | Σ residues | Total | Unattributed |", + "|---|---|---|---|"] + L += [f"| {w} | {_f(c['sum_residues'], 2)} | {_f(c['total'], 2)} | " + f"{_f(c['unattributed'], 2)}{'' if c['closes'] else ' **(> 1 kcal/mol)**'} |" + for w, c in clo.items()] + for w, rows in (s.get("hotspots") or {}).items(): + L += ["", f"Top residues, {w}:", "", "| Partner | Residue | ΔG (kcal/mol) |", + "|---|---|---|"] + L += [f"| {r.get('partner') or '—'} | {r['resname']} {r['residue']} | {_f(r['TOTAL'], 2)} |" + for r in rows[:15]] + c = s["checks"] + L += ["", "## Validation", "", + f"- Frames: {c['frames']}; GB and PB on identical frames: {c['same_frames_gb_pb']}.", + f"- Largest bonded Δ term: {_f(c['max_abs_bonded_delta'], 6)} kcal/mol " + f"({'single-trajectory consistent' if c['single_trajectory_consistent'] else '**inconsistent**'}).", + "", "## Limitations", "", + "- Implicit solvent, single-trajectory protocol, no conformational entropy unless its " + "validity limit is met.", + "- Ionic strength and settings are recorded in `prepare_meta.json`.", + "", "## Please cite", ""] + [f"- {c}" for c in s["method"]["citations"]] + return "\n".join(L) + "\n" + + +def write_all(ctx, formats=("md", "html")) -> list[Path]: + """Write every document whose results exist; returns written paths.""" + rd, out = ctx.config.report_dir, [] + s = _load(ctx.csv_path("pbc_summary.json")) + if s: + out += write_doc(rd / "PBC_VALIDATION.md", pbc_document(s), formats) + s = _load(ctx.csv_path("equilibration_summary.json")) + if s: + out += write_doc(rd / "EQUILIBRATION.md", equilibration_document(s), formats) + s = _load(ctx.csv_path("binding_energy_summary.json")) + if s: + out += write_doc(rd / "BINDING_ENERGY.md", + binding_energy_document(s, _load(ctx.csv_path("window_selection.json"))), + formats) + return out diff --git a/moldynx/report/generator.py b/moldynx/report/generator.py index f26d842..91aee59 100644 --- a/moldynx/report/generator.py +++ b/moldynx/report/generator.py @@ -8,7 +8,6 @@ from __future__ import annotations -import base64 from datetime import date from pathlib import Path @@ -39,43 +38,6 @@ def _md(blocks) -> str: return "\n".join(out) -def _html(blocks) -> str: - css = ("body{font-family:Arial,Helvetica,sans-serif;max-width:960px;margin:2rem auto;" - "padding:0 1rem;color:#222;line-height:1.5}h1{border-bottom:3px solid #1f6f8b}" - "h2{color:#1f6f8b;margin-top:2rem}table{border-collapse:collapse;margin:1rem 0}" - "th,td{border:1px solid #ccc;padding:6px 10px;text-align:left}" - "th{background:#f0f4f6}img{max-width:100%;border:1px solid #eee}" - "figure{margin:1.5rem 0}figcaption{color:#666;font-size:.9em}") - out = [f"" - f""] - for kind, payload in blocks: - if kind == "h1": - out.append(f"

    {payload}

    ") - elif kind == "h2": - out.append(f"

    {payload}

    ") - elif kind == "p": - out.append(f"

    {_inline_html(payload)}

    ") - elif kind == "table": - headers, rows = payload - out.append("" + "".join(f"" for h in headers) + "") - for r in rows: - out.append("" + "".join(f"" for c in r) + "") - out.append("
    {h}
    {c}
    ") - elif kind == "img": - path, caption = payload - p = Path(path) - if p.exists(): - b64 = base64.b64encode(p.read_bytes()).decode() - out.append(f"
    " - f"
    {caption}
    ") - out.append("") - return "\n".join(out) - - -def _inline_html(text: str) -> str: - return text.replace("**", "") # markdown bold markers -> plain (kept simple) - - # --------------------------------------------------------------------------- # def _fmt(x, nd=3): try: @@ -138,22 +100,33 @@ def generate_report(ctx, results: dict, manifest, formats=("md", "html")) -> lis "and per-analysis runtimes are recorded in `manifest.json` / " "`manifest.yaml` for exact reproduction.")) + # -- companion documents (written from the results files) ------------- # + from moldynx.report import documents + docs = documents.write_all(ctx, formats=[f for f in formats if f in ("md", "html")]) + doc_names = sorted({p.stem for p in docs}) + if doc_names: + blocks.insert(2, ("h2", "Companion documents")) + blocks.insert(3, ("p", " · ".join(f"[{n}]({n}.md)" for n in doc_names))) + # -- render ----------------------------------------------------------- # ctx.config.report_dir.mkdir(parents=True, exist_ok=True) - written = [] + written = list(docs) + md_text = _md(blocks) if "md" in formats: p = ctx.config.report_dir / "report.md" - p.write_text(_md(blocks), encoding="utf-8") + p.write_text(md_text, encoding="utf-8") written.append(p) + html_text = documents.to_html(md_text.replace(".md)", ".html)"), ctx.config.report_dir, + "MD analysis report") if "html" in formats: p = ctx.config.report_dir / "report.html" - p.write_text(_html(blocks), encoding="utf-8") + p.write_text(html_text, encoding="utf-8") written.append(p) if "pdf" in formats: try: from weasyprint import HTML pdf = ctx.config.report_dir / "report.pdf" - HTML(string=_html(blocks)).write_pdf(str(pdf)) + HTML(string=html_text).write_pdf(str(pdf)) written.append(pdf) except Exception: pass # weasyprint not installed; MD/HTML still produced @@ -172,6 +145,19 @@ def _key_result_rows(results: dict) -> list[list]: rows.append(["Most flexible residue", g(results, "rmsf", "most_flexible_resid")]) if "sasa" in results: rows.append(["Total SASA (mean)", f"{_fmt(g(results,'sasa','total_sasa','mean'),1)} nm²"]) + if "pbc_validation" in results: + rows.append(["PBC proof (all checks)", + "pass" if g(results, "pbc_validation", "all_checks_pass") else "FAIL"]) + if "interface" in results and g(results, "interface", "partners"): + rows.append(["Interface: min. distance (mean)", + f"{_fmt(g(results,'interface','min_interface_dist_nm','mean'))} nm"]) + rows.append(["Interface: core residues", + " / ".join(str(x) for x in (g(results, "interface", + "n_core_interface_residues") or []))]) + rows.append(["Interface: persistent contacts", g(results, "interface", "n_persistent_contacts")]) + if "equilibration" in results: + for v in (g(results, "equilibration", "verdicts") or [])[:3]: + rows.append(["Preparation", v]) return [r for r in rows if r[1] not in (None, "n/a")] @@ -187,9 +173,23 @@ def _interpretation(results: dict) -> list[str]: lines = [] conv = _nested(results, ("rmsd", "convergence", "converged")) if conv is not None: - lines.append(f"**Stability.** Backbone RMSD " - f"{'reached a plateau' if conv else 'had not fully plateaued'} " - f"over the analysed window.") + # an RMSD plateau alone is not evidence of stability or equilibrium + lines.append(f"**RMSD time course.** The backbone RMSD " + f"{'levels off' if conv else 'had not levelled off'} over the analysed " + f"window; this alone does not establish stability — see the convergence, " + f"interface and stationarity results.") + win = _nested(results, ("analysis_window", "binding_observables_stationary")) + if win is not None: + lines.append("**Stationarity.** " + ( + "Binding-related observables are stationary after " + f"{_fmt(_nested(results, ('analysis_window', 'conservative_t0_ns')), 1)} ns." + if win else "Binding-related observables are not stationary over the run: " + "averages describe a changing state (see the window selection).")) + be = _nested(results, ("mmpbsa", "headline")) + if be: + lines.append("**Binding energy.** " + "; ".join( + f"{h['method']} {h['window']}: {_fmt(h['mean'], 1)} ± {_fmt(h['sem'], 1)} kcal/mol" + for h in be) + " — end-point estimates, not experimental affinities.") rg = _nested(results, ("rog",)) if rg: init, fin = rg.get("rg_initial_nm"), rg.get("rg_final_nm") diff --git a/pyproject.toml b/pyproject.toml index 21714c0..24c9e2f 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -29,15 +29,18 @@ dependencies = [ "networkx>=3.0", "tqdm>=4.65", "pyyaml>=6.0", + "markdown>=3.4", # self-contained HTML documents ] [project.optional-dependencies] energy = ["panedr>=0.8"] # read .edr thermodynamics +annotation = ["biopython>=1.80"] # domain transfer / motif alignment fingerprints = ["prolif>=2.0", "rdkit>=2023.3"] # interaction fingerprints pdf = ["weasyprint>=60"] # PDF report export ui = ["rich>=13.0"] # nicer CLI output dev = ["pytest>=7.0", "pytest-cov", "ruff"] -all = ["panedr>=0.8", "prolif>=2.0", "rdkit>=2023.3", "weasyprint>=60", "rich>=13.0"] +all = ["panedr>=0.8", "biopython>=1.80", "prolif>=2.0", "rdkit>=2023.3", "weasyprint>=60", + "rich>=13.0"] [project.urls] Homepage = "https://github.com/SamDozer/molecular-dynamics-forge" diff --git a/tests/test_documents.py b/tests/test_documents.py new file mode 100644 index 0000000..fec91d2 --- /dev/null +++ b/tests/test_documents.py @@ -0,0 +1,67 @@ +"""Generated documents: content from results, honest wording, self-contained HTML.""" + +from __future__ import annotations + +import json +import re + +import numpy as np +import pytest + +from conftest import FIX +from moldynx.report import documents as docs +from test_equilibration import ctx as equil_ctx # noqa: F401 +from test_pbc import _chain, _run, _universe + + +def _html_ok(html: str, md: str): + text = re.sub(r"<[^>]+>", " ", re.sub(r"", "", html, flags=re.S)) + assert "**" not in text and "`" not in text # markdown fully rendered + assert html.count("") == len(re.findall(r"(?m)^\|[-| ]+\|$", md)) + assert not re.search(r'src="https?://', html) # self-contained + assert html.count("= {"EQUILIBRATION.md", "EQUILIBRATION.html"} + html = (equil_ctx.config.report_dir / "EQUILIBRATION.html").read_text(encoding="utf-8") + assert "data:image/png;base64," in html # figure embedded + _html_ok(html, md) From 2e8a369897c68df449f5ea743d991086d32de9ed Mon Sep 17 00:00:00 2001 From: Hossam Mahmoud Date: Sat, 26 Sep 2026 06:16:16 +0300 Subject: [PATCH 4/6] Datasets: assemble, verify, package - moldynx.dataset: numbered dataset layout from a run (report + companion documents re-linked, results, solute trajectory, validation, workflow, binding energy), raw-file manifest (files + SHA-256 with --include-raw, otherwise a stub), generated README. - verify_dataset.py shipped inside each dataset (layout, links, manifest, PBC proof, Rg smoke test); package() zips, extracts and re-verifies. - CLI: moldynx dataset / verify / package. Co-Authored-By: Claude Opus 5.5 --- moldynx/cli/main.py | 28 +++++ moldynx/dataset.py | 214 +++++++++++++++++++++++++++++++++++++ moldynx/verify_template.py | 97 +++++++++++++++++ tests/test_dataset.py | 55 ++++++++++ 4 files changed, 394 insertions(+) create mode 100644 moldynx/dataset.py create mode 100644 moldynx/verify_template.py create mode 100644 tests/test_dataset.py diff --git a/moldynx/cli/main.py b/moldynx/cli/main.py index 5824d90..d12ce72 100644 --- a/moldynx/cli/main.py +++ b/moldynx/cli/main.py @@ -76,6 +76,19 @@ def build_parser() -> argparse.ArgumentParser: be.add_argument("--execute", action="store_true", help="Run the prepared gmx_MMPBSA script now (hours; Linux or WSL).") + ds = sub.add_parser("dataset", help="Assemble a shareable, self-verifying dataset from a run.") + ds.add_argument("--run", required=True, help="Output folder of 'moldynx analyze'.") + ds.add_argument("--out", required=True, help="Dataset folder to create.") + ds.add_argument("--include-raw", action="store_true", + help="Copy and SHA-256 the canonical raw files (can be many GB).") + ds.add_argument("--title", help="Dataset title for the README.") + vf = sub.add_parser("verify", help="Run a dataset's verify_dataset.py.") + vf.add_argument("dataset") + vf.add_argument("--full", action="store_true", help="Also hash the raw files.") + pk = sub.add_parser("package", help="Zip a dataset, then extract and verify the zip.") + pk.add_argument("dataset") + pk.add_argument("--zip", help="Zip path (default: .zip).") + d = sub.add_parser("detect", help="Detect and print the system composition.") _add_common(d) @@ -233,6 +246,21 @@ def main(argv=None) -> int: return 1 print(f"MolDynX Tools (moldynx) {__version__}") return 0 + if args.command == "dataset": + from moldynx.dataset import assemble + print(f"[dataset] {assemble(args.run, args.out, include_raw=args.include_raw, title=args.title)}") + return 0 + if args.command == "verify": + from moldynx.dataset import verify + return verify(args.dataset, full=args.full) + if args.command == "package": + from pathlib import Path + from moldynx.dataset import package + r = package(args.dataset, args.zip or str(Path(args.dataset).resolve()) + ".zip") + print(f"[package] {r['zip']}: {r['files']} files, {r['bytes'] / 1e6:.1f} MB, " + f"CRC {'ok' if r['crc_ok'] else 'FAILED'}, verification exit {r['verify_exit']}") + print(r["verify_output"]) + return 0 if r["crc_ok"] and r["verify_exit"] == 0 else 1 return {"analyze": _cmd_analyze, "detect": _cmd_detect, "intake": _cmd_intake, "binding-energy": _cmd_binding_energy, "list-analyses": _cmd_list}[args.command](args) diff --git a/moldynx/dataset.py b/moldynx/dataset.py new file mode 100644 index 0000000..1bfbc87 --- /dev/null +++ b/moldynx/dataset.py @@ -0,0 +1,214 @@ +""" +Assemble a shareable, self-verifying dataset from a MolDynX run, verify it, and zip it. + +Layout:: + + / + README.md / .html generated overview + 0_raw_production/ canonical raw files (only with include_raw) + MANIFEST.sha256.tsv + 1_report/ report + companion documents (.md/.html) + figures/ + 2_tables/ 3_results/ summary tables, one CSV/JSON per analysis + 4_trajectory/ solute trajectory (core.pdb/.xtc) + core_meta.json + 5_validation/ intake report, PBC proof, window selection, verify_dataset.py + 6_workflow/ run manifest, configuration + 7_binding_energy/ gmx_MMPBSA inputs, run script and raw outputs (if prepared/run) + +``verify_dataset.py`` (standard library; MDAnalysis optional) checks the layout, that +every Markdown link/image resolves, the SHA-256 manifest, the PBC proof record and -- +if MDAnalysis is available -- recomputes Rg from the shipped trajectory. +""" + +from __future__ import annotations + +import hashlib +import json +import shutil +import tempfile +import zipfile +from pathlib import Path + +LAYOUT = ["1_report", "1_report/figures", "2_tables", "3_results", "4_trajectory", + "5_validation", "6_workflow"] +STUB = "RAW_DATA_NOT_INCLUDED.md" + + +def sha256(path: Path) -> str: + h = hashlib.sha256() + with open(path, "rb") as fh: + while b := fh.read(1 << 24): + h.update(b) + return h.hexdigest() + + +def _copytree(src: Path, dst: Path) -> None: + if src.is_dir(): + shutil.copytree(src, dst, dirs_exist_ok=True, + ignore=shutil.ignore_patterns("*_offsets.npz", "__pycache__", "*.lock")) + + +def _relink(md: str) -> str: + """Run layout (report/ next to figures/) -> dataset layout (1_report/figures).""" + return md.replace("](../figures/", "](figures/") + + +def assemble(run_dir: str | Path, out_dir: str | Path, include_raw: bool = False, + title: str | None = None) -> Path: + from moldynx.report.documents import to_html + run, out = Path(run_dir), Path(out_dir) + if out.exists() and any(out.iterdir()): + raise FileExistsError(f"{out} exists and is not empty") + for d in LAYOUT: + (out / d).mkdir(parents=True, exist_ok=True) + rep = run / "report" + for f in sorted(rep.glob("*.md")): + text = _relink(f.read_text(encoding="utf-8")) + (out / "1_report" / f.name).write_text(text, encoding="utf-8") + _copytree(run / "figures", out / "1_report" / "figures") + for f in sorted((out / "1_report").glob("*.md")): + text = f.read_text(encoding="utf-8").replace(".md)", ".html)") + f.with_suffix(".html").write_text(to_html(text, f.parent, f.stem), encoding="utf-8") + _copytree(run / "tables", out / "2_tables") + for f in sorted((run / "results").glob("*")): + if f.is_file(): + dst = "5_validation" if f.name.startswith(("pbc_", "window_selection")) else "3_results" + shutil.copy2(f, out / dst / f.name) + elif f.is_dir(): + _copytree(f, out / "3_results" / f.name) + for f in ("core.pdb", "core.xtc", "core_meta.json"): + if (run / "data" / f).exists(): + shutil.copy2(run / "data" / f, out / "4_trajectory" / f) + _copytree(run / "intake", out / "5_validation") + for f in ("manifest.json", "manifest.yaml"): + if (run / f).exists(): + shutil.copy2(run / f, out / "6_workflow" / f) + if (run / "binding_energy").is_dir(): + _copytree(run / "binding_energy", out / "7_binding_energy") + + manifest = _load_json(run / "manifest.json") or {} + intake = _load_json(run / "intake" / "intake_manifest.json") or {} + raw_rows = _raw_files(intake) + (out / "0_raw_production").mkdir(exist_ok=True) + lines = ["file\tbytes\tsha256\trole\toriginal_path"] + for role, p in raw_rows: + digest = sha256(p) if include_raw else "not-computed" + if include_raw: + shutil.copy2(p, out / "0_raw_production" / p.name) + lines.append(f"{p.name}\t{p.stat().st_size}\t{digest}\t{role}\t{p}") + (out / "0_raw_production" / "MANIFEST.sha256.tsv").write_text("\n".join(lines) + "\n", + encoding="utf-8") + if not include_raw: + (out / "0_raw_production" / STUB).write_text( + "# Raw simulation files are not included\n\nThe manifest lists the canonical files " + "this dataset was built from (sizes and original paths). Copy them into this folder " + "and run `python 5_validation/verify_dataset.py --full` to hash and check them.\n", + encoding="utf-8") + shutil.copy2(Path(__file__).parent / "verify_template.py", + out / "5_validation" / "verify_dataset.py") + readme = dataset_readme(title or run.name, manifest, intake, out) + (out / "README.md").write_text(readme, encoding="utf-8") + (out / "README.html").write_text(to_html(readme.replace(".md)", ".html)"), out, "README"), + encoding="utf-8") + return out + + +def _load_json(p: Path) -> dict | None: + try: + return json.loads(p.read_text(encoding="utf-8")) + except (OSError, ValueError): + return None + + +def _raw_files(intake: dict) -> list[tuple[str, Path]]: + fsd = intake.get("fileset", {}) + rows = [] + for role in ("trajectory", "topology", "energy", "log", "structure", "gmx_top", "index"): + if fsd.get(role): + rows.append((role, Path(fsd[role]))) + for stage, kinds in (fsd.get("stages") or {}).items(): + if stage in ("em", "nvt", "npt"): + for kind in ("log", "edr", "tpr"): + for p in kinds.get(kind, [])[:1]: + rows.append((f"{stage}_{kind}", Path(p))) + return [(r, p) for r, p in rows if p.exists()] + + +def dataset_readme(title: str, manifest: dict, intake: dict, out: Path) -> str: + sysd = manifest.get("system", {}) + an = manifest.get("analyses", {}) + ev = intake.get("evidence", {}) + L = [f"# {title}", "", + f"Molecular dynamics analysis dataset produced by **MolDynX Tools** " + f"v{manifest.get('moldynx_version', '?')} (formerly mdforge). Every number in the " + "documents below is computed from files in this folder.", "", + "## System at a glance", "", "| Property | Value |", "|---|---|", + f"| System type | {sysd.get('system_type', '—')} |", + f"| Atoms (run input) | {sysd.get('n_atoms', '—')} |", + f"| Protein chains | {sysd.get('n_protein_chains', '—')} " + f"{sysd.get('protein_chain_lengths', '')} |", + f"| Ions | {', '.join(f'{k} {v}' for k, v in (sysd.get('ion_counts') or {}).items()) or '—'} |", + f"| Production run | {ev.get('trajectory_choice', '—')} |", ""] + pl = ev.get("production_log") or {} + if pl: + L += [f"Production log: {pl.get('sessions')} mdrun session(s), " + f"{(pl.get('simulated_ps') or 0) / 1000:g} ns simulated.", ""] + L += ["## Key results", "", "| Analysis | Status | Headline |", "|---|---|---|"] + for name, rec in sorted(an.items()): + s = rec.get("summary") or {} + head = "" + if name == "pbc_validation": + head = "all proof checks pass" if s.get("all_checks_pass") else "see PBC_VALIDATION" + elif name == "mmpbsa" and s.get("headline"): + head = "; ".join(f"{h['method']} {h['window']} {h['mean']:.1f} ± {h['sem']:.1f}" + for h in s["headline"]) + elif name == "interface" and s.get("partners"): + head = f"{s.get('n_persistent_contacts')} persistent contacts" + elif name == "analysis_window": + head = s.get("statement", "") + status = s.get("status", rec.get("status")) + L.append(f"| {name} | {status} | {head} |") + docs = sorted(p.name for p in (out / "1_report").glob("*.md")) + L += ["", "## Documents", ""] + [f"- [{d[:-3]}](1_report/{d})" for d in docs] + L += ["", "## Folder map", "", "| Folder | Contents |", "|---|---|", + "| `0_raw_production/` | manifest of the canonical raw files (and the files, if included) |", + "| `1_report/` | report and companion documents, figures |", + "| `2_tables/`, `3_results/` | summary tables; one CSV/JSON per analysis |", + "| `4_trajectory/` | processed solute trajectory and its chain/PBC record |", + "| `5_validation/` | intake report, PBC proof, window selection, `verify_dataset.py` |", + "| `6_workflow/` | run manifest: versions, parameters, input fingerprints |", + "| `7_binding_energy/` | gmx_MMPBSA inputs, run script, outputs (if prepared) |", + "", "## Verifying this dataset", "", + "```bash", "python 5_validation/verify_dataset.py # --full also hashes raw files", + "```", "", "## Caveats", "", + "- Binding energies are end-point estimates, not experimental affinities.", + "- An RMSD plateau alone is not evidence of stability; see the stationarity analysis.", + "- Metadata not recorded in the simulation files (e.g. force-field variant) is stated " + "as not recorded, never guessed."] + return "\n".join(L) + "\n" + + +def package(dataset_dir: str | Path, zip_path: str | Path) -> dict: + """Deterministic zip (one top folder), then extract and run the shipped verifier.""" + import subprocess + import sys + src, zp = Path(dataset_dir), Path(zip_path) + files = sorted(p for p in src.rglob("*") if p.is_file() and "__pycache__" not in p.parts) + with zipfile.ZipFile(zp, "w", zipfile.ZIP_DEFLATED, compresslevel=9) as z: + for p in files: + z.write(p, f"{src.name}/{p.relative_to(src).as_posix()}") + with zipfile.ZipFile(zp) as z: + bad = z.testzip() + with tempfile.TemporaryDirectory() as tmp: + z.extractall(tmp) + r = subprocess.run([sys.executable, str(Path(tmp) / src.name / "5_validation" / + "verify_dataset.py")], + capture_output=True, text=True) + return {"zip": str(zp), "files": len(files), "bytes": zp.stat().st_size, + "crc_ok": bad is None, "verify_exit": r.returncode, "verify_output": r.stdout[-2000:]} + + +def verify(dataset_dir: str | Path, full: bool = False) -> int: + import subprocess + import sys + args = [sys.executable, str(Path(dataset_dir) / "5_validation" / "verify_dataset.py")] + return subprocess.call(args + (["--full"] if full else [])) + diff --git a/moldynx/verify_template.py b/moldynx/verify_template.py new file mode 100644 index 0000000..a645a7f --- /dev/null +++ b/moldynx/verify_template.py @@ -0,0 +1,97 @@ +""" +verify_dataset.py -- shipped inside every MolDynX Tools dataset. + + python 5_validation/verify_dataset.py [--full] + +1. layout: the expected folders exist and are not empty; +2. references: every Markdown link and image resolves; +3. raw files: 0_raw_production/ matches MANIFEST.sha256.tsv (sizes; SHA-256 with --full), + or is documented as not included; +4. PBC proof: the recorded checks passed; +5. smoke test (if MDAnalysis is installed): Rg recomputed from 4_trajectory/ matches + 3_results/rog.csv. + +Exit code 0 only if every check passes. Standard library only (MDAnalysis optional). +""" + +import csv +import hashlib +import json +import re +import sys +from pathlib import Path + +ROOT = Path(__file__).resolve().parents[1] +FAIL = [] + + +def check(ok, msg): + print((" ok " if ok else " FAIL ") + msg) + if not ok: + FAIL.append(msg) + + +def main(): + full = "--full" in sys.argv + print(f"MolDynX Tools dataset verification: {ROOT}") + print("[1] layout") + for d in ("1_report", "3_results", "4_trajectory", "5_validation", "6_workflow"): + p = ROOT / d + check(p.is_dir() and any(p.iterdir()), f"{d}/ present and not empty") + print("[2] references") + n = 0 + for md in ROOT.rglob("*.md"): + for target in re.findall(r"!?\[[^\]]*\]\(([^)\s]+)\)", md.read_text(encoding="utf-8")): + if re.match(r"^(https?:|mailto:|#)", target): + continue + n += 1 + ok = (md.parent / target.split("#")[0]).exists() + if not ok: + check(False, f"{md.relative_to(ROOT)} -> {target}") + check(True, f"{n} internal references checked") + print("[3] raw files") + man = ROOT / "0_raw_production" / "MANIFEST.sha256.tsv" + stub = ROOT / "0_raw_production" / "RAW_DATA_NOT_INCLUDED.md" + if man.exists(): + rows = list(csv.DictReader(man.open(encoding="utf-8"), delimiter="\t")) + present = [r for r in rows if (ROOT / "0_raw_production" / r["file"]).exists()] + if not present and stub.exists(): + check(True, f"{len(rows)} raw files listed; not included in this copy (see {stub.name})") + for r in present: + p = ROOT / "0_raw_production" / r["file"] + ok = p.stat().st_size == int(r["bytes"]) + if ok and full and r["sha256"] != "not-computed": + h = hashlib.sha256() + with open(p, "rb") as fh: + while b := fh.read(1 << 24): + h.update(b) + ok = h.hexdigest() == r["sha256"] + check(ok, f"raw file {r['file']}") + print("[4] PBC proof") + s = ROOT / "5_validation" / "pbc_summary.json" + if s.exists(): + checks = json.loads(s.read_text(encoding="utf-8")).get("checks", {}) + check(all(checks.values()), f"recorded PBC checks: {checks}") + print("[5] smoke test") + try: + import MDAnalysis as mda + import numpy as np + pdb, xtc = ROOT / "4_trajectory" / "core.pdb", ROOT / "4_trajectory" / "core.xtc" + rog = ROOT / "3_results" / "rog.csv" + if pdb.exists() and xtc.exists() and rog.exists(): + u = mda.Universe(str(pdb), str(xtc)) + ag = u.select_atoms("protein") or u.atoms + ref = [float(line.split(",")[1]) for line in rog.read_text().splitlines()[1:]] + k = len(u.trajectory) // 2 + u.trajectory[k] + check(abs(ag.radius_of_gyration() / 10 - ref[k]) < 1e-3, "Rg reproduced at the middle frame") + else: + print(" skip (no trajectory or rog.csv)") + except ImportError: + print(" skip (MDAnalysis not installed)") + print("RESULT:", "all checks passed" if not FAIL else f"{len(FAIL)} check(s) FAILED") + return 0 if not FAIL else 1 + + +if __name__ == "__main__": + sys.exit(main()) diff --git a/tests/test_dataset.py b/tests/test_dataset.py new file mode 100644 index 0000000..2356837 --- /dev/null +++ b/tests/test_dataset.py @@ -0,0 +1,55 @@ +"""Dataset assembly, the shipped verifier, and packaging.""" + +from __future__ import annotations + +import json +import subprocess +import sys + +from moldynx.dataset import assemble, package + + +def _fake_run(root): + (root / "report").mkdir(parents=True) + (root / "figures").mkdir() + (root / "results").mkdir() + (root / "intake").mkdir() + (root / "figures" / "pbc_validation.png").write_bytes(b"\x89PNG\r\n\x1a\n") + (root / "report" / "PBC_VALIDATION.md").write_text( + "# PBC\n\n## Observations\n\n![evidence](../figures/pbc_validation.png)\n", encoding="utf-8") + (root / "report" / "report.md").write_text("# Report\n\n[PBC](PBC_VALIDATION.md)\n", + encoding="utf-8") + (root / "results" / "pbc_summary.json").write_text(json.dumps({"checks": {"a": True}})) + (root / "results" / "rmsd.csv").write_text("time_ns,rmsd_nm\n0,0.1\n") + (root / "intake" / "INTAKE_REPORT.md").write_text("# Intake\n") + (root / "data").mkdir() + (root / "data" / "core_meta.json").write_text(json.dumps({"chains": []})) + (root / "manifest.json").write_text(json.dumps({"moldynx_version": "0.3.0", "system": {}, + "analyses": {"pbc_validation": { + "status": "ok", + "summary": {"all_checks_pass": True}}}})) + return root + + +def test_assemble_verify_package(tmp_path): + run = _fake_run(tmp_path / "run") + ds = assemble(run, tmp_path / "ds", title="Demo dataset") + assert (ds / "1_report" / "PBC_VALIDATION.html").exists() + assert "](figures/pbc_validation.png)" in (ds / "1_report" / "PBC_VALIDATION.md").read_text() + assert (ds / "5_validation" / "pbc_summary.json").exists() + assert (ds / "0_raw_production" / "RAW_DATA_NOT_INCLUDED.md").exists() + assert "pbc_validation | ok | all proof checks pass" in (ds / "README.md").read_text(encoding="utf-8") + r = subprocess.run([sys.executable, str(ds / "5_validation" / "verify_dataset.py")], + capture_output=True, text=True) + assert r.returncode == 0, r.stdout + z = package(ds, tmp_path / "ds.zip") + assert z["crc_ok"] and z["verify_exit"] == 0 + + +def test_verifier_catches_a_broken_link(tmp_path): + run = _fake_run(tmp_path / "run") + ds = assemble(run, tmp_path / "ds") + (ds / "1_report" / "report.md").write_text("# R\n\n![x](figures/missing.png)\n") + r = subprocess.run([sys.executable, str(ds / "5_validation" / "verify_dataset.py")], + capture_output=True, text=True) + assert r.returncode == 1 and "missing.png" in r.stdout From ebe5487f69844d1eff24e0f2f9d4a7d6795dd7c4 Mon Sep 17 00:00:00 2001 From: Hossam Mahmoud Date: Sat, 26 Sep 2026 06:22:51 +0300 Subject: [PATCH 5/6] Docs for MolDynX Tools 0.3: README, what's new vs mdforge, quickstart, next steps - README rewritten for MolDynX Tools: run stages, input-file contract, annotations schema, gmx_MMPBSA credit and citation, verification. - docs/WHATS_NEW_0.3.md: every change compared with mdforge 0.2. - QUICKSTART: intake, binding energy, datasets; conda env renamed to moldynx. - Example config: pbc, annotations, binding_energy. - CHANGELOG and NEXT_STEPS updated. Co-Authored-By: Claude Opus 5.5 --- CHANGELOG.md | 35 ++++- README.md | 190 ++++++++++++++----------- docs/NEXT_STEPS.md | 117 +++++---------- docs/QUICKSTART.md | 32 ++++- docs/WHATS_NEW_0.3.md | 72 ++++++++++ environment.yml | 6 +- examples/alpha_zein_A8HNE1/config.yaml | 15 ++ 7 files changed, 297 insertions(+), 170 deletions(-) create mode 100644 docs/WHATS_NEW_0.3.md diff --git a/CHANGELOG.md b/CHANGELOG.md index e893288..15b5870 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -84,4 +84,37 @@ `moldynx.core.surface.shrake_rupley` computes frame by frame; `sasa` and `interface` buried area now stride by 10 frames by default (parameters `stride` / `bsa_stride`). -Remaining 0.3 work: `docs/NEXT_STEPS.md`. +### Complex analyses, annotations, stationarity + +- `contact_lifetime`, `water_bridges`, `porcupine` (ported from a validated legacy pipeline, + generalised; porcupine uses an SVD and its own aligned copy of the trajectory). +- `annotation`: biological numbering, domains transferred by global alignment (BLOSUM62, + gap −10/−0.5) with uncertain edges flagged, motif location with partial matches, docking-site + involvement — only from the `annotations:` block of the configuration. +- `analysis_window`: Chodera equilibration detection, drift and half-vs-half tests; the + averaging window is a user decision (`binding_energy.primary_window`) when nothing is + stationary. +- `moldynx.statistics`: statistical inefficiency, corrected SEM, equilibration detection, drift. + +### MM-GBSA / MM-PBSA with gmx_MMPBSA + +- **Replaced** the placeholder `mmpbsa` (fixed groups `r 1-100`/`r 101-9999`, an Amber + force-field line in a CHARMM/GROMACS workflow, the full solvated system, never run; the + contact-count proxy is dropped — interface occupancy covers it) with a complete workflow built + on [gmx_MMPBSA](https://github.com/Valdes-Tresanca-MS/gmx_MMPBSA): protein-only system, + groups from chain identity, derived ionic strength, explicit decomposition residues, probe gate, + GB + PB, analysis with autocorrelation-corrected errors, drift, GB vs PB, hotspots, closure and + an entropy validity gate. New `moldynx binding-energy [--execute]`. Reproduces a manual + analysis of real outputs exactly. + +### Documents and datasets + +- `PBC_VALIDATION`, `EQUILIBRATION`, `BINDING_ENERGY` documents written from the results files; + Markdown rendered to self-contained HTML (the old renderer only deleted `**`); report wording + no longer claims stability from an RMSD plateau. +- `moldynx dataset` / `verify` / `package`: numbered dataset, raw-file manifest, + `verify_dataset.py` shipped inside, zip extracted and re-verified. +- New dependency: `markdown`; optional `biopython` (`[annotation]`). + +See `docs/WHATS_NEW_0.3.md` for the comparison with mdforge 0.2 and `docs/NEXT_STEPS.md` for +what is left. diff --git a/README.md b/README.md index abfd629..e3c932e 100644 --- a/README.md +++ b/README.md @@ -1,104 +1,130 @@ -# MolDynX Tools — reusable, reproducible, audited GROMACS MD analysis +# MolDynX Tools — audited, reproducible GROMACS MD analysis [![CI](https://github.com/SamDozer/molecular-dynamics-forge/actions/workflows/ci.yml/badge.svg)](https://github.com/SamDozer/molecular-dynamics-forge/actions) [![License: MIT](https://img.shields.io/badge/License-MIT-yellow.svg)](LICENSE) [![Python](https://img.shields.io/badge/python-3.10%2B-blue.svg)](pyproject.toml) [![DOI](https://zenodo.org/badge/DOI/10.5281/zenodo.21265946.svg)](https://doi.org/10.5281/zenodo.21265946) -**MolDynX Tools** (`moldynx`; formerly *mdforge*) turns a GROMACS simulation directory into a complete, publication-quality, -fully reproducible analysis — with minimal input. Point it at a folder; it discovers -the files, **detects the system** (protein / ligand / DNA / RNA / membrane / ions / -multi-chain / …), **auto-selects the right analyses**, runs them with streaming-friendly -performance, and produces figures, tables, a provenance manifest, and a report. +**MolDynX Tools** (`moldynx`; formerly *mdforge*) turns a GROMACS simulation folder into a +complete, **audited** analysis: it finds the right files by evidence, proves that trajectory +processing changed nothing it should not, audits how the system was prepared and equilibrated, +analyses the structure, dynamics and interfaces, computes MM-GBSA/MM-PBSA binding energies with +gmx_MMPBSA, and writes documents and a self-verifying dataset in which every number traces back +to a file. -> This repository began as a project-specific pipeline for an α-zein simulation -> (now preserved in [`legacy/`](legacy) and [`examples/`](examples)) and was -> refactored into this general toolkit. +> What is new compared with mdforge 0.2: [docs/WHATS_NEW_0.3.md](docs/WHATS_NEW_0.3.md). --- -## Highlights - -- **Zero-config detection** — recursively finds `*.xtc/.trr/.tpr/.gro/.edr/.ndx/.top`, - validates them (clear messages for anything missing), and classifies the system. -- **Automatic module selection** — each analysis declares the system types and files - it supports; the pipeline runs exactly what applies (`--plan` shows *why*). -- **Extensible via plugins** — drop a `BaseAnalysis` subclass into - `moldynx/analysis/plugins/` (or a `--plugin-dir`) and it is auto-discovered. -- **Config-driven** — describe a whole run in `config.yaml` and re-run with one command - (ideal for HPC/batch). -- **Reproducible by construction** — every run writes `manifest.json/yaml` with library - versions, git commit, input-file fingerprints, seeds, parameters and runtimes. -- **Publication-quality output** — 300-dpi PNG **and** vector PDF, consistent - Nature-like style, plus Markdown + self-contained HTML (+ optional PDF) reports. -- **Scales** — streams frame-by-frame and caches a solute-only trajectory, so large - (100 GB+) explicit-solvent runs stay tractable. - -## Supported systems - -protein-only · protein–protein · protein–peptide · protein–ligand · protein–DNA · -protein–RNA · protein–membrane · multi-chain · ions · cofactors · mixed biomolecular -systems (detected automatically; override with `--system-type`). - -## Install +## Quick start ```bash git clone https://github.com/SamDozer/molecular-dynamics-forge cd molecular-dynamics-forge -python -m pip install -e ".[all]" # or ".[dev]" for tests +python -m pip install -e ".[all]" # or ".[dev]" for tests + +moldynx intake --input /path/to/sim_dir # which files form the run, what is missing, why +moldynx analyze --input /path/to/sim_dir # everything applicable -> ./moldynx_results/ +moldynx binding-energy --input /path/to/sim_dir [--execute] # MM-GBSA/PBSA (gmx_MMPBSA) +moldynx dataset --run moldynx_results/ --out my_dataset +moldynx package my_dataset # zip, extract, verify ``` -## Usage +MolDynX never writes into the simulation folder. + +## What a run does + +| Stage | What happens | Output | +|---|---|---| +| **Intake** | files classified by simulation stage (setup / minimisation / NVT / NPT / production); the production run chosen by evidence (atom counts, time span, log completion) — never by size or name; earlier analysis folders, backups and crash dumps ignored; extensions, clock offsets and inputs referenced by job scripts but absent are reported | `intake/INTAKE_REPORT.md`, capability matrix | +| **PBC: diagnose → treat → prove** | each raw frame is measured (split molecules, continuity, PBC-aware partner distance), treated (`--pbc auto\|none\|whole\|nojump`), then checked against the raw frame | `PBC_VALIDATION.md`, `pbc_summary.json` | +| **Preparation audit** | minimisation outcome (and where the largest force sits), NVT/NPT statistics and residual drift, position restraints, protonation states, chain of custody, timeline | `EQUILIBRATION.md` | +| **Analyses** | RMSD, RMSF, Rg, SASA, PCA, DCCM, clustering, secondary structure, H-bonds, salt bridges, RIN, ProLIF, … plus for complexes: interface suite, contact lifetimes, water bridges, porcupine | figures, CSVs | +| **Annotations** | display names, biological numbering, domains transferred by alignment, motifs, docking site — only from your config | `annotation.json`, `region_summary.csv` | +| **Stationarity** | Chodera equilibration detection and drift tests; no silent choice of averaging window | `window_selection.json` | +| **Binding energy** | gmx_MMPBSA package (protein-only system, derived ionic strength, explicit decomposition residues, probe gate) → analysis with autocorrelation-corrected errors | `BINDING_ENERGY.md` | +| **Documents & dataset** | report + companion documents (Markdown and self-contained HTML), numbered dataset with `verify_dataset.py` | `report/`, dataset folder, zip | + +## The input-file contract + +| If this is present… | …this becomes possible | Without it | +|---|---|---| +| production run input (`.tpr`) + trajectory with equal atom counts | every analysis | **stop** | +| production `.edr` | energy analysis | skipped, with reason | +| production `.log` | Methods, sessions, extension audit, completion proof | completion reported as not proven | +| `topol.top` + `toppar/` | MM-GBSA/PBSA | skipped; package still prepared | +| minimisation / NVT / NPT `.log` + `.edr` | preparation audit | listed as missing | +| stage `.tpr` files | restraint audit | "not recoverable" | +| stage `.mdp` files | MDP-only settings (`gen-vel`, `define`) | "not recorded" | +| job scripts | explicit chain of custody | inferred, labelled inferred | + +## Configuration + +Everything can be driven from one YAML file (`moldynx analyze --config config.yaml`); see +[`examples/alpha_zein_A8HNE1/config.yaml`](examples/alpha_zein_A8HNE1/config.yaml). Facts that +cannot be read from files are supplied — never guessed: + +```yaml +annotations: + chains: + - {segid: seg_0_PROA, display: "Partner A", role: ligand} + - {segid: seg_1_PROB, display: "Partner B", role: receptor, numbering_offset: 187} + domains: + seg_1_PROB: {reference_name: "UniProt P11021", reference_sequence: "MKLS...", + regions: {NBD: [26, 405]}} + motifs: + - {name: "Motif 1", sequence: "CSQAPIASLLPPYLSPAVSSVC", chain: seg_0_PROA} + docking_site: {seg_1_PROB: [405, 434, 435]} + unresolved_metadata: {force_field_variant: null, salt_concentration_M: null} +binding_energy: + primary_window: [80, 100] # ns; omit and MolDynX shows the full run + final 20 % +``` + +## Binding energy (gmx_MMPBSA) + +MM-GBSA/MM-PBSA calculations are performed with +[**gmx_MMPBSA**](https://github.com/Valdes-Tresanca-MS/gmx_MMPBSA) (Valdés-Tresanca et al., +*J. Chem. Theory Comput.* 2021, 17, 6281) on top of AmberTools MMPBSA.py (Miller et al., 2012); +please cite both. MolDynX prepares the calculation (protein-only system; ionic strength derived +from the ions in the run input; decomposition of every residue that contacted the partner in any +frame), gates it (a short probe run must show zero bonded Δ terms), runs it on Linux/WSL, and +analyses it (windows, drift, GB vs PB, hotspots, decomposition closure, entropy validity). +Results are end-point estimates, not experimental affinities. + +## Reproducibility and verification + +Every run writes `manifest.json/yaml` (versions, git commit, parameters, input fingerprints and +the evidence for every chosen file). Every dataset ships `5_validation/verify_dataset.py` +(layout, links, raw-file manifest, PBC proof, Rg smoke test): ```bash -moldynx detect --input /path/to/sim_dir # what's in my system? -moldynx analyze --input /path/to/sim_dir --plan # what would run, and why? -moldynx analyze --input /path/to/sim_dir -o results # run everything applicable -moldynx analyze --config examples/alpha_zein_A8HNE1/config.yaml # reproducible -moldynx list-analyses # registered analyses (incl. plugins) +moldynx verify my_dataset [--full] ``` -See [`docs/QUICKSTART.md`](docs/QUICKSTART.md) for all options and the plugin template. - ## Architecture ``` moldynx/ - core/ system.py (detection) · base.py (BaseAnalysis) · registry.py (+plugins) - context.py · config.py (YAML+CLI) · provenance.py · pipeline.py - io/ discovery.py · validation.py - statistics/ descriptive · timeseries · correlation · bootstrap - plotting/ style · figures (PNG+PDF, 300 dpi) - analysis/ rmsd · rmsf · rog · sasa · … + plugins/ (auto-discovered) - report/ generator (Markdown + HTML + PDF) - cli/ main (analyze/detect/list) · interactive -tests/ · examples/ · docs/ · legacy/ · Dockerfile · .github/workflows/ + core/ system (detection, chain records) · pbc · surface · annotations · context · pipeline + registry (+plugins) · config (YAML+CLI) · provenance + io/ discovery (stage classification, evidence) · validation (capability matrix) + gromacs (log/mdp/xtc/tpr readers, gmx runner with WSL fallback) · jobscripts · intake + analysis/ 32 analyses (incl. the bundled plugin) + plugins/ (auto-discovered) + binding/ gmx_MMPBSA prepare + analyse + report/ generator + documents (Markdown → self-contained HTML) + dataset.py · verify_template.py +mdforge/ deprecated import shim (removed in 0.4) ``` -Every analysis subclasses `BaseAnalysis`, declaring `required_files`, -`supported_systems`, `outputs` and `default_params`; subclasses auto-register, so -**detection → selection → run → report** is entirely data-driven. - -## Analyses (24 built-in) - -| System scope | Analyses | -|---|---| -| **Any system** | RMSD, radius of gyration, SASA, COM, H-bonds, energies (`.edr`), ProLIF, MM/PBSA workflow, statistics | -| **Protein** | RMSF, structural descriptors (Dmax/κ²/volume), native contacts (Q), contact map, DSSP secondary structure, RIN, PCA + free-energy landscape, DCCM, clustering, salt bridges, convergence (RMSIP/block-avg), end-to-end *(plugin)* | -| **Complex** (protein–protein/–nucleic) | interface (BSA, contacts, iRMSD) | -| **Protein–ligand** | ligand RMSD, ligand contacts, binding pocket | -| **Protein–DNA/RNA** | protein–nucleic contacts, nucleic RMSD | +Every analysis subclasses `BaseAnalysis` and declares `required_files`, `supported_systems`, +`outputs` and `default_params`; `moldynx analyze --plan` shows what runs and why. -Each is a drop-in `BaseAnalysis`; the pipeline runs only those applicable to the -detected system (`moldynx list-analyses` shows all; `--plan` shows what runs and why). -The complex/ligand/nucleic modules are implemented and gate correctly but await -validation on a matching test trajectory. +## Tested -## Reproducibility - -```bash -cat results/manifest.json # versions, git commit, seeds, params, input hashes, runtimes -``` +66 tests run in CI on real, trimmed data from two 100 ns protein–protein simulations (GROMACS +logs, an energy file, gmx_MMPBSA outputs, raw-folder listings) — see +[`tests/fixtures/README.md`](tests/fixtures/README.md). Intake, PBC, preparation audit and +MM-GBSA/PBSA analysis reproduce an independent manual analysis of those datasets exactly. ## Container @@ -109,23 +135,15 @@ docker run --rm -v /data/sim:/sim moldynx analyze --input /sim --output /sim/res ## Roadmap -See **[ROADMAP.md](ROADMAP.md)** for the plan — the flagship being a -**comparison mode** (`moldynx compare`) that overlays control vs. protein–ligand / -protein–protein systems on shared axes (ΔRMSF maps, common-subspace PCA, ensemble -similarity), plus parallel execution, a functional API, membrane and multi-engine -support — drawing design influence from -[MDAnalysis](https://github.com/MDAnalysis/mdanalysis) and -[mdtraj](https://github.com/mdtraj/mdtraj). +See [ROADMAP.md](ROADMAP.md) — next: comparison mode (`moldynx compare`), parallel execution, +a process-pool SASA backend, membrane and multi-engine support. ## Citation -If you use MolDynX Tools, please cite it (concept DOI — always resolves to the latest version): - -> Mahmoud, H. *MolDynX Tools: a reusable, reproducible analysis framework for GROMACS -> molecular dynamics simulations.* Zenodo. https://doi.org/10.5281/zenodo.21265946 +> Mahmoud, H. *MolDynX Tools: a reusable, reproducible analysis framework for GROMACS molecular +> dynamics simulations.* Zenodo. https://doi.org/10.5281/zenodo.21265946 -A machine-readable [`CITATION.cff`](CITATION.cff) is included (GitHub shows a -"Cite this repository" button). +A machine-readable [`CITATION.cff`](CITATION.cff) is included. ## License diff --git a/docs/NEXT_STEPS.md b/docs/NEXT_STEPS.md index daa9efa..809c78e 100644 --- a/docs/NEXT_STEPS.md +++ b/docs/NEXT_STEPS.md @@ -1,79 +1,38 @@ -# MolDynX Tools 0.3 — remaining work - -Work paused after Phase 5 (part 1) on branch `feature/moldynx-dataset-pipeline`. The goal of -0.3: running MolDynX on any GROMACS/CHARMM-GUI simulation folder always yields the same audited, -reproducible dataset. The quality baseline is two hand-made analyses of protein–protein -complexes (α-zein A8HNE1–ZmBiP2, α-zein Q946V6–ZmBiP2). They are **test cases only**: no code -may name a system, chain, residue range or folder. The full specification is in the handoff -package `MolDynX_handoff_2026-09-23/PROMPT.md` (kept outside the repo); the essentials are below. - -## Done - -| Phase | Result | -|---|---| -| 1 Rename | mdforge → MolDynX Tools (`moldynx`), `mdforge` shim for one version | -| 2 Intake | evidence-based discovery, capability matrix, `moldynx intake`, INTAKE_REPORT | -| 3 PBC | one-pass diagnose → treat → prove (`core/pbc.py`), chain identity, cache key, `pbc_validation` | -| 4 Equilibration | `equilibration` analysis: stage parameters, minimisation outcome + Fmax location, crash dumps, energy statistics, position restraints (`gmx dump`), protonation, chain of custody, timeline, figure | -| 5 (part) | `interface` rewritten (full core lists, residues ever within 6 Å, buried area, iRMSD, trends); `core/surface.py` works around an mdtraj multi-frame SASA bug | - -Phases 2–4 reproduce the manual audit of both reference datasets exactly (see CHANGELOG). - -## Remaining - -### Phase 5 — finish complex analyses -- Run `interface` end to end on a real 100 ns complex; compare with the legacy values - (A8HNE1: 27.9 nm² buried, mean minimum distance 0.266 nm, 34 persistent contacts). -- Port from the legacy pipeline (`alpha-zein-md-analysis/scripts/`), generalised, as - `BaseAnalysis` classes: `contact_lifetime`, `dynamic_network`, `water_bridges`, `porcupine`, - eigenvalue spectrum into `pca`. Do **not** port `report_advanced.py` or legacy `mmgbsa.py`. -- SASA speed: single-frame calls are correct but serial (~4 s/frame for 13 k atoms); - `sasa` and `interface` stride by 10. Parallelise with a process pool, then restore stride 1. - -### Phase 6 — annotation and analysis window -- `annotations:` in the YAML config (template written by `moldynx intake --detect`): - chain display names and roles, biological numbering offsets, homology reference, motifs file, - docking-site residues, unresolved metadata (force-field variant, salt, protonation rationale — - never guessed). -- Domains by homology only (Needleman–Wunsch, BLOSUM62, gap −10/−0.5; CATH edges of a reference - structure through SIFTS; map author numbering to UniProt by alignment; DSSP cross-check; flag - uncertain edges). Motifs only from a user-supplied file; report partial matches honestly. -- Window selection: Chodera t₀ per binding-relevant observable, slope + Welch half-vs-half; - if nothing is stationary, report time-resolved values and ask for the primary window. - -### Phase 7 — MM-GBSA / MM-PBSA (rewrite `analysis/mmpbsa.py`) -- prepare → probe → run → analyse. Protein-only `complex.tpr` (convert-tpr on the Protein group), - `complex.top` with protein molecule types only, `complex.ndx` from chain identity - (0 = receptor, 1 = ligand); **no Amber `forcefields=` line**. -- Probe gate: 11 frames GB; abort unless every bonded Δ term is exactly 0. -- GB: igb 5, PBRadii 3, idecomp 2, IE + C2 entropy (reported only if σ < 3.6 / 6.0 kcal/mol). - PB: PBRadii 7, inp 1, fillratio 4. Same frames (1 ns spacing). Ionic strength **derived** from - ion counts (`SystemInfo.ion_counts`) and box volume, or labelled assumed. -- Decomposition: explicit `print_res` = `interface` residues within 6 Å ever - (`residues_within_6A_ever`), closure check Σ residues vs ΔG_bind (Q946V6 legacy gap +3.73 - kcal/mol at 80–100 ns should close). -- MPI ranks ≤ physical cores; PB ≈ 3.3 GB per rank. Analysis: autocorrelation-corrected SEM, - windows, drift, GB vs PB, hotspots, region sums. Test data: `tests/fixtures/gmx_mmpbsa/`. - -### Phase 8 — documents -- Data-driven README, AUDIT_PROVENANCE, PBC_VALIDATION, EQUILIBRATION, BINDING_ENERGY, - analysis report; Observations / Interpretation / Limitations / Conclusions; no hard-coded - system names or conclusions; never call a system "stable" from an RMSD plateau (current - `report/generator.py` `_interpretation` does — fix); render Markdown properly to - self-contained HTML; validator (pair `**` across the whole document). - -### Phase 9 — dataset, verify, package -- `moldynx dataset` (numbered layout 0_raw_production … 7_binding_energy with SHA-256 manifests), - `moldynx verify` (layout, references, manifests, Rg smoke test), `moldynx package` - (analysis-only zip when raw data > ~1 GB; stub injected into the zip only; extract + verify). - -### Phase 10 — wrap-up -- README rewrite (input-file contract, annotations schema, new commands), QUICKSTART, - CITATION/.zenodo titles (done), ROADMAP ticks, fix README claim about a `legacy/` folder. -- The user renames the GitHub repository and creates the release (Zenodo DOI); not automated. - -## Acceptance values (manual audit, for the desktop data) - -A8HNE1 — MM-GBSA 0–100 ns −55.2 ± 4.0, 80–100 −69.8 ± 3.7; MM-PBSA −57.2 ± 7.4 / −76.6 ± 4.6 -kcal/mol; r(GB,PB) 0.86. Q946V6 — MM-GBSA −33.4 ± 2.1 / −37.8 ± 2.4; MM-PBSA −40.7 ± 2.5 / -−47.1 ± 3.3; r 0.73; entropy σ 86.7 → invalid. +# MolDynX Tools — next steps after 0.3 + +All ten phases of the 0.3 plan are implemented (see [WHATS_NEW_0.3.md](WHATS_NEW_0.3.md) and +the CHANGELOG). The reference datasets used as the quality bar (two 100 ns protein–protein +complexes) are test cases only; nothing in the code names a system. + +## Left for the maintainer (cannot be automated) + +1. **Rename the GitHub repository** (Settings → General → Repository name), e.g. to + `moldynx-tools`; GitHub redirects the old URL. Then replace `molecular-dynamics-forge` with + the new name in `README.md` (badges, clone commands), `docs/QUICKSTART.md`, `CITATION.cff` + and `pyproject.toml` (7 URLs). +2. **Release 0.3.0** on GitHub when satisfied; Zenodo mints a new version DOI under the existing + concept DOI 10.5281/zenodo.21265946. +3. Decide whether the real test fixtures (trimmed logs and gmx_MMPBSA outputs from the reference + simulations) should stay public or be replaced by synthetic files. + +## Not yet validated end to end on real data + +- `interface`, `contact_lifetime`, `water_bridges`, `porcupine`, `annotation` and the full + `analyze` → `dataset` → `package` chain were tested on real geometry and synthetic + trajectories, not yet on a full 100 ns run. Reference values from the manual analysis of one + complex: 27.9 nm² buried area, 0.266 nm mean minimum distance, 34 persistent contacts. +- `moldynx binding-energy --execute` was prepared and gated on real data (complex.tpr, topology, + index, derived ionic strength 0.158 M); a full GB + PB run (≈ 2.5 h) through the new script is + still to be done. The explicit decomposition residue set should close the +3.7 kcal/mol gap + seen with `print_res = "within 6"` in one reference system. + +## Engineering backlog + +- SASA is computed frame by frame (mdtraj 1.11.1 multi-frame bug) and runs on one core; + `sasa` and the buried area stride by 10. A process-pool backend would allow stride 1. +- Comparison mode (`moldynx compare`, ROADMAP v0.3) should open with a table of preparation + differences (equilibration temperature, barostat, minimisation outcome, protonation) — the + equilibration audit already produces every field. +- Domain boundaries from CATH/SIFTS for a reference structure are supplied by the user today; + fetching and caching them would remove a manual step. +- Remove the `mdforge` compatibility shim in 0.4. diff --git a/docs/QUICKSTART.md b/docs/QUICKSTART.md index 36b4a21..28f0d96 100644 --- a/docs/QUICKSTART.md +++ b/docs/QUICKSTART.md @@ -9,10 +9,20 @@ python -m pip install -e ".[all]" # core + energy + fingerprints + pdf + ui ``` or with conda: ```bash -conda env create -f environment.yml && conda activate zein-md +conda env create -f environment.yml && conda activate moldynx python -m pip install -e . ``` +## 0. Check the inputs + +```bash +moldynx intake --input /path/to/simulation_dir [--detect] +``` +Writes `INTAKE_REPORT.md`: which files form the production run and the evidence for it, each +simulation stage, what every missing file disables, run extensions, temperature changes between +stages and inputs your job scripts reference but that are absent. Nothing is written into the +simulation folder. + ## 1. Detect your system ```bash @@ -59,6 +69,26 @@ Every run writes a `manifest.json` capturing library versions, git commit, input fingerprints, seeds, parameters and per-analysis runtimes — enough for another researcher to reproduce the analysis exactly. +## 5. Binding energy (MM-GBSA / MM-PBSA with gmx_MMPBSA) + +```bash +moldynx binding-energy --input /path/to/simulation_dir # prepare the package +moldynx binding-energy --input /path/to/simulation_dir --execute # run it (Linux/WSL, hours) +``` +Needs [gmx_MMPBSA](https://github.com/Valdes-Tresanca-MS/gmx_MMPBSA) in a conda environment +(default name `gmxMMPBSA`, override with `MOLDYNX_MMPBSA_ENV`), e.g. +`conda create -n gmxMMPBSA -c conda-forge --override-channels python=3.12 "ambertools>=24.8,<27" +"mpi4py>=4.0.1,<5" "numpy<2"` then `pip install gmx_MMPBSA`. When the outputs exist, `analyze` +reads them and writes `BINDING_ENERGY.md`. + +## 6. Share a dataset + +```bash +moldynx dataset --run moldynx_results/ --out my_dataset [--include-raw] +moldynx package my_dataset # zip -> extract -> run verify_dataset.py +moldynx verify my_dataset --full # re-check at any time +``` + ## 5. Add your own analysis (plugin) Drop a file into `moldynx/analysis/plugins/` (or any `--plugin-dir`): diff --git a/docs/WHATS_NEW_0.3.md b/docs/WHATS_NEW_0.3.md new file mode 100644 index 0000000..9e4f5da --- /dev/null +++ b/docs/WHATS_NEW_0.3.md @@ -0,0 +1,72 @@ +# MolDynX Tools 0.3 — what changed compared with mdforge 0.2 + +mdforge 0.2 was a good general framework: plugin registry, automatic system detection and +analysis selection, a run manifest, figures and a report. MolDynX Tools 0.3 keeps all of that +and adds what was missing to *trust* the results: evidence for every input choice, proof for +every trajectory transformation, an audit of how the system was prepared, binding energies +that are actually computed, and documents in which every number traces back to a file. + +The quality bar was two hand-made, fully audited analyses of 100 ns protein–protein complexes. +Where MolDynX now automates a step of those analyses, the tests check it reproduces them. + +## 1. Choosing the input files + +| mdforge 0.2 | MolDynX Tools 0.3 | +|---|---| +| Picked the **largest** file per role. On a real CHARMM-GUI folder this took an *equilibration* run input as the topology (it was larger than the production one) and a minimisation crash dump (`stepNc.pdb`) as the structure. | Files are classified by simulation stage. The production trajectory is the longest complete one whose atom count equals the run input's, read from file headers. Ties are reported as ambiguities, never resolved silently. | +| Scanned every subfolder, including earlier analysis outputs. | Ignores earlier analysis folders, GROMACS backups, caches and crash dumps, and says so. | +| Presence checks only. | Capability matrix: each missing file and what it disables. | +| — | `moldynx intake` / `INTAKE_REPORT.md`: extended runs (declared vs. executed steps), temperature changes between stages, the cluster-clock offset, inputs that job scripts reference but that are absent. | +| — | Read-only readers for mdrun logs, `.mdp`, XTC frame headers (no cache files written into the data folder) and `.tpr` headers; a GROMACS runner that falls back to WSL on Windows. | + +## 2. Periodic boundaries + +| mdforge 0.2 | MolDynX Tools 0.3 | +|---|---| +| `unwrap` inside `try/except: pass`: a topology without bonds silently gave a trajectory that was never made whole. | One pass per frame: **diagnose** the raw coordinates, **treat** them (`none`/`whole`/`nojump`, first-frame clustering for complexes), **prove** the result (only whole-box translations, bonds short, partner distances preserved). Failures are recorded, never swallowed. | +| Cached solute trajectory reused whenever the files existed. | Cache keyed on input fingerprints, PBC mode, frame slice and selection. | +| Chain identity lost: PDB truncates CHARMM-GUI segment IDs (`seg_0_PROA`, `seg_1_PROB` → `seg_`), so `interface` could not find two partners. | Chain identity persisted as atom-index ranges; every complex analysis uses it. | + +## 3. Preparation and equilibration (new) + +Minimisation outcome (states plainly when the force tolerance was *not* reached) and the residue +carrying the largest residual force; NVT/NPT statistics, settling times and residual drift; +position restraints read from each run input (and none in production); protonation states; +crash-dump accounting; chain of custody from job scripts; stage timeline — as +`EQUILIBRATION.md`. On one reference system it located the largest residual force on a +protonated aspartate, which the manual audit had not connected. + +## 4. Complex analyses + +| mdforge 0.2 | MolDynX Tools 0.3 | +|---|---| +| `interface` "implemented but not validated"; atom-pair contacts including hydrogens. | Rewritten from a validated suite: residue–residue heavy-atom contacts, full (never truncated) core interface lists, persistent contacts, residues within 6 Å at any time, buried area, interface Cα RMSD, trends. | +| — | `contact_lifetime`, `water_bridges`, `porcupine` (SVD, no 3N × 3N matrix). | +| — | `annotation`: biological numbering, domains transferred by sequence alignment (uncertain edges flagged), motif location with partial matches, docking-site involvement — only from user input. | +| — | `analysis_window`: Chodera equilibration detection and drift tests; the averaging window is a user decision when nothing is stationary. | +| SASA via one multi-frame `mdtraj.shrake_rupley` call. | **Bug found:** mdtraj 1.11.1 returns wrong values for some frames of a multi-frame call (±0.5 nm² on a real trajectory; negative buried areas in a rigid-body test). MolDynX computes frame by frame. | + +## 5. Binding energy + +| mdforge 0.2 | MolDynX Tools 0.3 | +|---|---| +| Wrote a placeholder script: groups `r 1-100` / `r 101-9999`, an Amber `forcefields = "oldff/leaprc.ff99SB, leaprc.gaff"` line in a CHARMM/GROMACS workflow, the full solvated system, never ran. A contact-count "proxy" was shown under an MM/PBSA label. | Full [gmx_MMPBSA](https://github.com/Valdes-Tresanca-MS/gmx_MMPBSA) workflow: protein-only system cut from the run input, groups from chain identity, **ionic strength derived** from the ions in the system, **explicit decomposition residues** (every residue that contacted the partner in any frame — `within 6` misses late binders), a **probe gate** (bonded Δ terms must be zero), GB and PB runs, and an analysis with autocorrelation-corrected errors, drift, GB vs PB, hotspots in biological numbering, decomposition closure and an entropy validity gate. Reproduces a manual analysis of real outputs exactly. gmx_MMPBSA and MMPBSA.py are credited and cited in every output. | + +## 6. Documents and datasets + +| mdforge 0.2 | MolDynX Tools 0.3 | +|---|---| +| One generic report; the HTML renderer deleted `**` and ignored inline code, links and subscripts. | Report plus `PBC_VALIDATION`, `EQUILIBRATION` and `BINDING_ENERGY` documents, written from the results files, with Observations / Interpretation / Limitations; proper Markdown → self-contained HTML. | +| Interpretation said "**Stability.** Backbone RMSD reached a plateau". | Wording rules: no stability claim from an RMSD plateau; binding energies are end-point estimates, not affinities; entropy only when valid; unknown metadata stated as not recorded. | +| — | `moldynx dataset` / `verify` / `package`: numbered dataset, raw-file manifest, `verify_dataset.py` shipped inside, zip extracted and re-verified. | + +## 7. Engineering + +- Renamed to **MolDynX Tools** (`moldynx`); `import mdforge`, the `mdforge` command and the + `mdforge.plugins` entry-point group keep working until 0.4 with a deprecation warning. +- Autocorrelation statistics (`statistical_inefficiency`, corrected SEM, equilibration + detection, drift) in `moldynx.statistics`. +- Tests: from 7 unit tests to a suite running on real, trimmed data from two production + simulations (logs, energy file, gmx_MMPBSA outputs, raw-folder listings). +- Default output `./moldynx_results/`; MolDynX never writes into the simulation + folder. diff --git a/environment.yml b/environment.yml index 6dc7e71..728544c 100644 --- a/environment.yml +++ b/environment.yml @@ -1,10 +1,10 @@ -# Conda environment for the alpha-zein (A8HNE1) MD analysis pipeline. +# Conda environment for MolDynX Tools. # conda env create -f environment.yml -# conda activate zein-md +# conda activate moldynx # # conda-forge provides pre-built wheels for rdkit / mdtraj / MDAnalysis, which is # the most reliable way to install the interaction-fingerprint stack. -name: zein-md +name: moldynx channels: - conda-forge dependencies: diff --git a/examples/alpha_zein_A8HNE1/config.yaml b/examples/alpha_zein_A8HNE1/config.yaml index c2bd1d0..e84e7fd 100644 --- a/examples/alpha_zein_A8HNE1/config.yaml +++ b/examples/alpha_zein_A8HNE1/config.yaml @@ -35,3 +35,18 @@ params: ref_frame: 0 rog: window: 20 + +# Periodic-boundary treatment: auto (nojump for complexes) | none | whole | nojump +pbc: auto + +# What only you can supply -- never guessed (all optional) +annotations: + chains: + - {segid: seg_0_PROA, display: "alpha-zein A8HNE1", role: ligand} + - {segid: seg_1_PROB, display: "ZmBiP2", role: receptor} # numbering_offset defaults to 1-based + motifs: [] # e.g. {name: "Motif 1", sequence: "...", chain: seg_0_PROA} + docking_site: {} # e.g. {seg_1_PROB: [405, 434, 435]} (biological numbering) + unresolved_metadata: {force_field_variant: null, salt_concentration_M: null} + +binding_energy: + primary_window: [80, 100] # ns; omit to show the full run and its final 20 % From b105eb739431f7d20043b74370b19c60453e77c0 Mon Sep 17 00:00:00 2001 From: Hossam Mahmoud Date: Sat, 26 Sep 2026 06:23:37 +0300 Subject: [PATCH 6/6] Point every link to the renamed repository SamDozer/MolDynX-Tools Co-Authored-By: Claude Opus 5.5 --- CITATION.cff | 4 ++-- README.md | 6 +++--- docs/NEXT_STEPS.md | 6 ++---- docs/QUICKSTART.md | 4 ++-- pyproject.toml | 4 ++-- 5 files changed, 11 insertions(+), 13 deletions(-) diff --git a/CITATION.cff b/CITATION.cff index 1df228b..94a0efb 100644 --- a/CITATION.cff +++ b/CITATION.cff @@ -11,8 +11,8 @@ authors: - family-names: Mahmoud given-names: Hossam orcid: "https://orcid.org/0009-0004-7804-4174" -repository-code: "https://github.com/SamDozer/molecular-dynamics-forge" -url: "https://github.com/SamDozer/molecular-dynamics-forge" +repository-code: "https://github.com/SamDozer/MolDynX-Tools" +url: "https://github.com/SamDozer/MolDynX-Tools" license: MIT version: 0.2.1 date-released: "2026-07-08" diff --git a/README.md b/README.md index e3c932e..de07a6c 100644 --- a/README.md +++ b/README.md @@ -1,6 +1,6 @@ # MolDynX Tools — audited, reproducible GROMACS MD analysis -[![CI](https://github.com/SamDozer/molecular-dynamics-forge/actions/workflows/ci.yml/badge.svg)](https://github.com/SamDozer/molecular-dynamics-forge/actions) +[![CI](https://github.com/SamDozer/MolDynX-Tools/actions/workflows/ci.yml/badge.svg)](https://github.com/SamDozer/MolDynX-Tools/actions) [![License: MIT](https://img.shields.io/badge/License-MIT-yellow.svg)](LICENSE) [![Python](https://img.shields.io/badge/python-3.10%2B-blue.svg)](pyproject.toml) [![DOI](https://zenodo.org/badge/DOI/10.5281/zenodo.21265946.svg)](https://doi.org/10.5281/zenodo.21265946) @@ -19,8 +19,8 @@ to a file. ## Quick start ```bash -git clone https://github.com/SamDozer/molecular-dynamics-forge -cd molecular-dynamics-forge +git clone https://github.com/SamDozer/MolDynX-Tools +cd MolDynX-Tools python -m pip install -e ".[all]" # or ".[dev]" for tests moldynx intake --input /path/to/sim_dir # which files form the run, what is missing, why diff --git a/docs/NEXT_STEPS.md b/docs/NEXT_STEPS.md index 809c78e..98d44cd 100644 --- a/docs/NEXT_STEPS.md +++ b/docs/NEXT_STEPS.md @@ -6,10 +6,8 @@ complexes) are test cases only; nothing in the code names a system. ## Left for the maintainer (cannot be automated) -1. **Rename the GitHub repository** (Settings → General → Repository name), e.g. to - `moldynx-tools`; GitHub redirects the old URL. Then replace `molecular-dynamics-forge` with - the new name in `README.md` (badges, clone commands), `docs/QUICKSTART.md`, `CITATION.cff` - and `pyproject.toml` (7 URLs). +1. ~~Rename the GitHub repository~~ — done: `SamDozer/MolDynX-Tools` (the old URL redirects); + every link in the repository points to the new name. 2. **Release 0.3.0** on GitHub when satisfied; Zenodo mints a new version DOI under the existing concept DOI 10.5281/zenodo.21265946. 3. Decide whether the real test fixtures (trimmed logs and gmx_MMPBSA outputs from the reference diff --git a/docs/QUICKSTART.md b/docs/QUICKSTART.md index 28f0d96..1facb0b 100644 --- a/docs/QUICKSTART.md +++ b/docs/QUICKSTART.md @@ -3,8 +3,8 @@ ## Install ```bash -git clone https://github.com/SamDozer/molecular-dynamics-forge -cd molecular-dynamics-forge +git clone https://github.com/SamDozer/MolDynX-Tools +cd MolDynX-Tools python -m pip install -e ".[all]" # core + energy + fingerprints + pdf + ui ``` or with conda: diff --git a/pyproject.toml b/pyproject.toml index 24c9e2f..6ea207c 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -43,8 +43,8 @@ all = ["panedr>=0.8", "biopython>=1.80", "prolif>=2.0", "rdkit>=2023.3", "weasyp "rich>=13.0"] [project.urls] -Homepage = "https://github.com/SamDozer/molecular-dynamics-forge" -Repository = "https://github.com/SamDozer/molecular-dynamics-forge" +Homepage = "https://github.com/SamDozer/MolDynX-Tools" +Repository = "https://github.com/SamDozer/MolDynX-Tools" [project.scripts] moldynx = "moldynx.cli.main:main"