From 5bab45fba190be3edae4bdeafa3d50af876aacef Mon Sep 17 00:00:00 2001 From: Hossam Mahmoud Date: Tue, 7 Jul 2026 17:32:55 +0300 Subject: [PATCH 1/4] mdforge: foundation of the molecular-dynamics-forge framework General, extensible, reproducible analysis of GROMACS MD simulations of arbitrary biomolecular systems. Point it at a simulation directory; it discovers files, detects the system (protein/ligand/DNA/RNA/membrane/ions/multi-chain), auto-selects analyses, runs them (streaming + solute-extraction), and writes figures, a reproducibility manifest, and Markdown/HTML reports. Includes: core (system detection, BaseAnalysis, plugin registry, YAML+CLI config, provenance manifest, pipeline), io (discovery+validation), statistics (incl. bootstrap CI), plotting (PNG+PDF 300dpi), report (MD+HTML+PDF), analyses (rmsd/rmsf/rog/sasa + example plugin), CLI (analyze/detect/list/--plan/interactive), packaging (pyproject + `mdforge` entry point), tests (7 passing), Dockerfile, GitHub Actions CI, and docs. Co-Authored-By: Claude Opus 4.8 --- .github/workflows/ci.yml | 28 ++ Dockerfile | 16 + README.md | 105 ++++++ docs/QUICKSTART.md | 80 +++++ environment.yml | 27 ++ examples/alpha_zein_A8HNE1/README.md | 26 ++ examples/alpha_zein_A8HNE1/config.yaml | 37 ++ mdforge/__init__.py | 38 +++ mdforge/analysis/__init__.py | 7 + mdforge/analysis/plugins/__init__.py | 8 + .../analysis/plugins/example_end_to_end.py | 53 +++ mdforge/analysis/rmsd.py | 64 ++++ mdforge/analysis/rmsf.py | 62 ++++ mdforge/analysis/rog.py | 55 +++ mdforge/analysis/sasa.py | 64 ++++ mdforge/cli/__init__.py | 1 + mdforge/cli/interactive.py | 23 ++ mdforge/cli/main.py | 142 ++++++++ mdforge/core/__init__.py | 1 + mdforge/core/base.py | 92 +++++ mdforge/core/config.py | 119 +++++++ mdforge/core/context.py | 135 ++++++++ mdforge/core/pipeline.py | 130 +++++++ mdforge/core/provenance.py | 128 +++++++ mdforge/core/registry.py | 140 ++++++++ mdforge/core/system.py | 316 ++++++++++++++++++ mdforge/io/__init__.py | 4 + mdforge/io/discovery.py | 120 +++++++ mdforge/io/validation.py | 67 ++++ mdforge/plotting/__init__.py | 6 + mdforge/plotting/figures.py | 22 ++ mdforge/plotting/style.py | 45 +++ mdforge/report/__init__.py | 3 + mdforge/report/generator.py | 202 +++++++++++ mdforge/statistics/__init__.py | 8 + mdforge/statistics/correlation.py | 10 + mdforge/statistics/descriptive.py | 64 ++++ mdforge/statistics/timeseries.py | 65 ++++ pyproject.toml | 60 ++++ requirements.txt | 24 ++ tests/test_core.py | 101 ++++++ 41 files changed, 2698 insertions(+) create mode 100644 .github/workflows/ci.yml create mode 100644 Dockerfile create mode 100644 README.md create mode 100644 docs/QUICKSTART.md create mode 100644 environment.yml create mode 100644 examples/alpha_zein_A8HNE1/README.md create mode 100644 examples/alpha_zein_A8HNE1/config.yaml create mode 100644 mdforge/__init__.py create mode 100644 mdforge/analysis/__init__.py create mode 100644 mdforge/analysis/plugins/__init__.py create mode 100644 mdforge/analysis/plugins/example_end_to_end.py create mode 100644 mdforge/analysis/rmsd.py create mode 100644 mdforge/analysis/rmsf.py create mode 100644 mdforge/analysis/rog.py create mode 100644 mdforge/analysis/sasa.py create mode 100644 mdforge/cli/__init__.py create mode 100644 mdforge/cli/interactive.py create mode 100644 mdforge/cli/main.py create mode 100644 mdforge/core/__init__.py create mode 100644 mdforge/core/base.py create mode 100644 mdforge/core/config.py create mode 100644 mdforge/core/context.py create mode 100644 mdforge/core/pipeline.py create mode 100644 mdforge/core/provenance.py create mode 100644 mdforge/core/registry.py create mode 100644 mdforge/core/system.py create mode 100644 mdforge/io/__init__.py create mode 100644 mdforge/io/discovery.py create mode 100644 mdforge/io/validation.py create mode 100644 mdforge/plotting/__init__.py create mode 100644 mdforge/plotting/figures.py create mode 100644 mdforge/plotting/style.py create mode 100644 mdforge/report/__init__.py create mode 100644 mdforge/report/generator.py create mode 100644 mdforge/statistics/__init__.py create mode 100644 mdforge/statistics/correlation.py create mode 100644 mdforge/statistics/descriptive.py create mode 100644 mdforge/statistics/timeseries.py create mode 100644 pyproject.toml create mode 100644 requirements.txt create mode 100644 tests/test_core.py diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml new file mode 100644 index 0000000..19e3eac --- /dev/null +++ b/.github/workflows/ci.yml @@ -0,0 +1,28 @@ +name: CI + +on: + push: + branches: [main, mdforge-framework] + pull_request: + +jobs: + test: + runs-on: ubuntu-latest + strategy: + fail-fast: false + matrix: + python-version: ["3.10", "3.11", "3.12"] + steps: + - uses: actions/checkout@v4 + - name: Set up Python ${{ matrix.python-version }} + uses: actions/setup-python@v5 + with: + python-version: ${{ matrix.python-version }} + - name: Install package (+ dev extras) + run: | + python -m pip install --upgrade pip + python -m pip install -e ".[dev]" + - name: Lint (ruff, non-blocking) + run: ruff check mdforge || true + - name: Run tests + run: pytest -q diff --git a/Dockerfile b/Dockerfile new file mode 100644 index 0000000..3a70f2b --- /dev/null +++ b/Dockerfile @@ -0,0 +1,16 @@ +# mdforge -- reproducible container image. +# Build: docker build -t mdforge . +# Run: docker run --rm -v /data/sim:/sim mdforge analyze --input /sim --output /sim/results +FROM mambaorg/micromamba:1.5-jammy + +WORKDIR /app +COPY --chown=$MAMBA_USER:$MAMBA_USER environment.yml /app/environment.yml +RUN micromamba install -y -n base -f /app/environment.yml && \ + micromamba clean --all --yes + +COPY --chown=$MAMBA_USER:$MAMBA_USER . /app +ARG MAMBA_DOCKERFILE_ACTIVATE=1 +RUN pip install --no-deps -e . + +ENTRYPOINT ["mdforge"] +CMD ["--help"] diff --git a/README.md b/README.md new file mode 100644 index 0000000..08f2737 --- /dev/null +++ b/README.md @@ -0,0 +1,105 @@ +# mdforge — reusable, 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) + +**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. + +> 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. + +--- + +## 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 + `mdforge/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 + +```bash +git clone https://github.com/SamDozer/molecular-dynamics-forge +cd molecular-dynamics-forge +python -m pip install -e ".[all]" # or ".[dev]" for tests +``` + +## Usage + +```bash +mdforge detect --input /path/to/sim_dir # what's in my system? +mdforge analyze --input /path/to/sim_dir --plan # what would run, and why? +mdforge analyze --input /path/to/sim_dir -o results # run everything applicable +mdforge analyze --config examples/alpha_zein_A8HNE1/config.yaml # reproducible +mdforge list-analyses # registered analyses (incl. plugins) +``` + +See [`docs/QUICKSTART.md`](docs/QUICKSTART.md) for all options and the plugin template. + +## Architecture + +``` +mdforge/ + 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/ +``` + +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 + +Structure/dynamics analyses ship built-in and more are being ported from the +reference pipeline: **RMSD, RMSF, radius of gyration, SASA** (built-in now); +**H-bonds, PCA + free-energy landscape, DCCM, contact maps, clustering, DSSP, +salt bridges, native contacts, energies, convergence, RIN, ProLIF, MM/PBSA +workflow, interface & buried-surface analyses** (from [`legacy/`](legacy), +porting incrementally — each is a drop-in `BaseAnalysis`). + +## Reproducibility + +```bash +cat results/manifest.json # versions, git commit, seeds, params, input hashes, runtimes +``` + +## Container + +```bash +docker build -t mdforge . +docker run --rm -v /data/sim:/sim mdforge analyze --input /sim --output /sim/results +``` + +## License + +MIT — see [LICENSE](LICENSE). diff --git a/docs/QUICKSTART.md b/docs/QUICKSTART.md new file mode 100644 index 0000000..ad57c4b --- /dev/null +++ b/docs/QUICKSTART.md @@ -0,0 +1,80 @@ +# mdforge — Quick Start + +## Install + +```bash +git clone https://github.com/SamDozer/molecular-dynamics-forge +cd molecular-dynamics-forge +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 +python -m pip install -e . +``` + +## 1. Detect your system + +```bash +mdforge detect --input /path/to/simulation_dir +``` +Prints the discovered files, validation status, and the detected composition +(protein chains, ligands, nucleic acids, lipids, ions, water) and system type. + +## 2. See what would run (dry-run) + +```bash +mdforge analyze --input /path/to/simulation_dir --plan +``` +Shows which analyses are selected for the detected system — and *why* each other +one was skipped. + +## 3. Run the full analysis + +```bash +mdforge analyze --input /path/to/simulation_dir --output results/ +``` +Discovers files → detects the system → auto-selects analyses → runs them → +writes CSVs, 300-dpi PNG+PDF figures, a `manifest.json/yaml` provenance record, +and a Markdown+HTML report. + +### Handy options +```bash +--stride 10 # analyse every 10th frame (quick look) +--analyses rmsd,rog,sasa # run a specific subset +--exclude sasa # skip an analysis +--all # run every applicable analysis +--interactive # ask when the choice is ambiguous +--report md,html,pdf # report formats +--plugin-dir ./my_plugins # load extra drop-in analyses +--config config.yaml # drive everything from a YAML file +``` + +## 4. Reproducible, config-driven workflow + +```bash +mdforge analyze --config examples/alpha_zein_A8HNE1/config.yaml +``` +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. Add your own analysis (plugin) + +Drop a file into `mdforge/analysis/plugins/` (or any `--plugin-dir`): + +```python +from mdforge.core.base import BaseAnalysis + +class MyAnalysis(BaseAnalysis): + name = "my_analysis" + label = "My analysis" + supported_systems = {"*"} + def run(self, ctx): + u = ctx.core_universe() + ... # compute + ctx.write_csv(df, "my_analysis.csv") + return {"figure": "my_analysis"} +``` +It is auto-discovered and appears in `mdforge list-analyses` immediately. +See `mdforge/analysis/plugins/example_end_to_end.py` for a complete template. diff --git a/environment.yml b/environment.yml new file mode 100644 index 0000000..6dc7e71 --- /dev/null +++ b/environment.yml @@ -0,0 +1,27 @@ +# Conda environment for the alpha-zein (A8HNE1) MD analysis pipeline. +# conda env create -f environment.yml +# conda activate zein-md +# +# 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 +channels: + - conda-forge +dependencies: + - python=3.11 + - numpy>=1.24 + - scipy>=1.10 + - pandas>=2.0 + - matplotlib>=3.7 + - seaborn>=0.12 + - tqdm>=4.65 + - pyyaml>=6.0 + - scikit-learn>=1.3 + - networkx>=3.0 + - mdanalysis>=2.7 + - mdtraj>=1.10 + - rdkit>=2023.3 + - pip + - pip: + - panedr>=0.8 + - prolif>=2.0 diff --git a/examples/alpha_zein_A8HNE1/README.md b/examples/alpha_zein_A8HNE1/README.md new file mode 100644 index 0000000..01395d1 --- /dev/null +++ b/examples/alpha_zein_A8HNE1/README.md @@ -0,0 +1,26 @@ +# Example: α-zein (A8HNE1) + +This directory holds the **reference example** that seeded the toolkit: a 100 ns +GROMACS simulation of an AlphaFold 3 α-zein model (187 residues, explicit +solvent, CHARMM36). + +- `config.yaml` — a complete, reproducible mdforge workflow for this system. +- `figures/`, `results/`, `tables/`, `report/` — outputs produced by the original + project-specific pipeline (now preserved in [`../../legacy/`](../../legacy)), + kept here as a worked example of what a full analysis looks like. + +## Reproduce with mdforge + +```bash +# edit input_dir in config.yaml to point at your copy of the raw files, then: +mdforge analyze --config examples/alpha_zein_A8HNE1/config.yaml +``` + +The raw simulation files (`step5_production.xtc/.tpr/.edr/...`) are **not** +included (multi-GB); point `input_dir` at your own copy. + +## Key finding (reference) + +The extended AlphaFold model **compacts** (Rg 4.2 → 3.5 nm) while retaining its +helical content (~55 %) and most native contacts (Q ≈ 0.89) — characteristic of +an elongated, flexible seed-storage protein. See `report/analysis_report.md`. diff --git a/examples/alpha_zein_A8HNE1/config.yaml b/examples/alpha_zein_A8HNE1/config.yaml new file mode 100644 index 0000000..cd491e7 --- /dev/null +++ b/examples/alpha_zein_A8HNE1/config.yaml @@ -0,0 +1,37 @@ +# --------------------------------------------------------------------------- +# Reproducible mdforge workflow for the alpha-zein (A8HNE1) system. +# +# Run it with a single command: +# mdforge analyze --config examples/alpha_zein_A8HNE1/config.yaml +# +# CLI flags override anything set here, so the same file works on HPC/batch. +# --------------------------------------------------------------------------- + +# Point this at your GROMACS simulation directory (searched recursively). +input_dir: /path/to/A8HNE1_MDS +output_dir: examples/alpha_zein_A8HNE1/mdforge_out + +# Trajectory slicing / performance +start: 0 +end: null # null => last frame +stride: 1 # analyse every frame (increase for a quick look) +core_extract: true # cache a solute-only trajectory for fast structural analyses + +# Analysis selection +analyses: null # null => auto-select for the detected system +exclude: [] # e.g. [sasa] to skip +plugin_dirs: [] # extra folders of drop-in analysis plugins + +# Reproducibility / reporting +seed: 0 +report_formats: [md, html] + +# Per-analysis parameter overrides (optional) +params: + rmsd: + ref_frame: 0 + window: 20 + rmsf: + ref_frame: 0 + rog: + window: 20 diff --git a/mdforge/__init__.py b/mdforge/__init__.py new file mode 100644 index 0000000..aee1447 --- /dev/null +++ b/mdforge/__init__.py @@ -0,0 +1,38 @@ +""" +mdforge -- a reusable, extensible framework for analysing GROMACS molecular +dynamics simulations of arbitrary biomolecular systems. + +Point it at a simulation directory and it will discover the files, detect the +system composition (protein / ligand / nucleic acid / membrane / ions / ...), +select the appropriate analyses, run them with streaming-friendly performance, +and produce publication-quality figures plus a fully reproducible report. + +Public API +---------- + from mdforge import analyze, detect_system, __version__ +""" + +from __future__ import annotations + +__version__ = "0.2.0" +__author__ = "Hossam Mahmoud" + +# Re-export the most commonly used entry points (kept import-light). +from mdforge.core.system import detect_system, SystemInfo, SystemType, ComponentType # noqa: E402 +from mdforge.core.registry import registry # noqa: E402 + +__all__ = [ + "__version__", + "detect_system", + "SystemInfo", + "SystemType", + "ComponentType", + "registry", + "analyze", +] + + +def analyze(*args, **kwargs): + """Convenience wrapper around :func:`mdforge.core.pipeline.run_pipeline`.""" + from mdforge.core.pipeline import run_pipeline + return run_pipeline(*args, **kwargs) diff --git a/mdforge/analysis/__init__.py b/mdforge/analysis/__init__.py new file mode 100644 index 0000000..9cd1636 --- /dev/null +++ b/mdforge/analysis/__init__.py @@ -0,0 +1,7 @@ +""" +Built-in analysis modules. + +Every module here defines one or more :class:`~mdforge.core.base.BaseAnalysis` +subclasses, which auto-register on import. The registry imports all modules in +this package (except ``plugins``, handled separately) at selection time. +""" diff --git a/mdforge/analysis/plugins/__init__.py b/mdforge/analysis/plugins/__init__.py new file mode 100644 index 0000000..a80f203 --- /dev/null +++ b/mdforge/analysis/plugins/__init__.py @@ -0,0 +1,8 @@ +""" +Plugin package. + +Drop a ``*.py`` file here (or point ``--plugin-dir`` at an external folder) that +defines a :class:`~mdforge.core.base.BaseAnalysis` subclass, and it is +auto-discovered and exposed through the CLI without touching the core. +See ``example_end_to_end.py`` for a minimal template. +""" diff --git a/mdforge/analysis/plugins/example_end_to_end.py b/mdforge/analysis/plugins/example_end_to_end.py new file mode 100644 index 0000000..baaaf86 --- /dev/null +++ b/mdforge/analysis/plugins/example_end_to_end.py @@ -0,0 +1,53 @@ +""" +Example plugin: protein end-to-end distance. + +This is a complete, working template showing how to add a new analysis to +mdforge WITHOUT modifying the core. Copy this file, rename the class, implement +``run()``, and it will be auto-discovered and selectable on the CLI as +``mdforge analyze --analyses end_to_end`` (or run automatically for the systems +it declares support for). +""" + +from __future__ import annotations + +import numpy as np +import pandas as pd + +from mdforge.core.base import BaseAnalysis +from mdforge.core.system import SystemType +from mdforge import plotting +from mdforge import statistics as st +from mdforge.plotting import PALETTE + + +class EndToEndDistance(BaseAnalysis): + name = "end_to_end" + label = "End-to-end distance" + category = "structure" + required_files = {"trajectory", "topology"} + supported_systems = {SystemType.PROTEIN_ONLY, SystemType.PROTEIN_PEPTIDE} + outputs = ["results/end_to_end.csv", "figures/end_to_end.png"] + default_params = {"window": 20} + + def run(self, ctx) -> dict: + p = self.params(ctx) + plotting.set_style() + u = ctx.core_universe() + ca = u.select_atoms("protein and name CA") + d = np.empty(len(u.trajectory)); times = np.empty_like(d) + for i, ts in enumerate(ctx.iter_frames(u, desc="[end_to_end]")): + pos = ca.positions / 10.0 + d[i] = float(np.linalg.norm(pos[0] - pos[-1])) + times[i] = ts.time / 1000.0 + + ctx.write_csv(pd.DataFrame({"time_ns": times, "end_to_end_nm": d, + "movavg_nm": st.moving_average(d, p["window"])}), + "end_to_end.csv") + fig, ax = plotting.new_axes() + ax.plot(times, d, color=PALETTE["primary"], alpha=0.4, lw=1.0) + ax.plot(times, st.moving_average(d, p["window"]), + color=PALETTE["primary"], lw=2.2, label="MA") + ax.set_xlabel("Time (ns)"); ax.set_ylabel("End-to-end distance (nm)") + ax.set_title(self.label); ax.legend(loc="best"); ax.margins(x=0.01) + plotting.save_figure(fig, ctx.fig_path("end_to_end"), dpi=ctx.config.dpi) + return {"end_to_end": st.describe(d, "end_to_end_nm"), "figure": "end_to_end"} diff --git a/mdforge/analysis/rmsd.py b/mdforge/analysis/rmsd.py new file mode 100644 index 0000000..c7ac85b --- /dev/null +++ b/mdforge/analysis/rmsd.py @@ -0,0 +1,64 @@ +"""Backbone and Cα RMSD versus the reference frame (universal).""" + +from __future__ import annotations + +import numpy as np +import pandas as pd +from MDAnalysis.analysis import rms + +from mdforge.core.base import BaseAnalysis +from mdforge import plotting +from mdforge import statistics as st +from mdforge.plotting import PALETTE + + +class RMSDAnalysis(BaseAnalysis): + name = "rmsd" + label = "Backbone & Cα RMSD" + category = "structure" + required_files = {"trajectory", "topology"} + supported_systems = {"*"} + outputs = ["results/rmsd.csv", "figures/rmsd.png", "figures/rmsd.pdf"] + default_params = {"ref_frame": 0, "window": 20} + + def _series(self, u, select, ref_frame): + R = rms.RMSD(u, select=select, ref_frame=ref_frame).run() + return R.results.rmsd[:, 2] / 10.0 # Å -> nm + + def run(self, ctx) -> dict: + p = self.params(ctx) + plotting.set_style() + u = ctx.core_universe() + # protein if present, else nucleic backbone + has_protein = ctx.system.flags.get("has_protein", True) + bb_sel = "backbone" if has_protein else "nucleic and name P O3' O5' C3' C4' C5'" + ca_sel = ctx.selection("ca", "protein and name CA") if has_protein \ + else "nucleic and name P" + + time_ns = ctx.times_ns(u) + rmsd_bb = self._series(u, bb_sel, p["ref_frame"]) + rmsd_ca = self._series(u, ca_sel, p["ref_frame"]) + ma_bb = st.moving_average(rmsd_bb, p["window"]) + ma_ca = st.moving_average(rmsd_ca, p["window"]) + + df = pd.DataFrame({"time_ns": time_ns, + "rmsd_backbone_nm": rmsd_bb, "rmsd_backbone_movavg_nm": ma_bb, + "rmsd_ca_nm": rmsd_ca, "rmsd_ca_movavg_nm": ma_ca}) + ctx.write_csv(df, "rmsd.csv") + + fig, ax = plotting.new_axes() + ax.plot(time_ns, rmsd_bb, color=PALETTE["primary"], alpha=0.35, lw=1.0) + ax.plot(time_ns, ma_bb, color=PALETTE["primary"], lw=2.2, label="Backbone") + ax.plot(time_ns, rmsd_ca, color=PALETTE["secondary"], alpha=0.35, lw=1.0) + ax.plot(time_ns, ma_ca, color=PALETTE["secondary"], lw=2.2, label=r"C$\alpha$") + plateau = float(np.mean(rmsd_bb[len(rmsd_bb)//2:])) + plotting.add_reference_line(ax, plateau, label=f"plateau = {plateau:.3f} nm") + ax.set_xlabel("Time (ns)"); ax.set_ylabel("RMSD (nm)") + ax.set_title(self.label); ax.legend(loc="upper left"); ax.margins(x=0.01) + plotting.save_figure(fig, ctx.fig_path("rmsd"), dpi=ctx.config.dpi) + + conv = st.plateau_detection(time_ns, rmsd_bb) + return {"backbone": st.describe(rmsd_bb, "rmsd_backbone_nm"), + "ca": st.describe(rmsd_ca, "rmsd_ca_nm"), + "plateau_mean_nm": plateau, "convergence": conv, + "figure": "rmsd"} diff --git a/mdforge/analysis/rmsf.py b/mdforge/analysis/rmsf.py new file mode 100644 index 0000000..c2336bd --- /dev/null +++ b/mdforge/analysis/rmsf.py @@ -0,0 +1,62 @@ +"""Per-residue Cα RMSF with flexible-region highlighting (protein systems).""" + +from __future__ import annotations + +import numpy as np +import pandas as pd +from MDAnalysis.analysis import align, rms + +from mdforge.core.base import BaseAnalysis +from mdforge.core.system import SystemType +from mdforge import plotting +from mdforge import statistics as st +from mdforge.plotting import PALETTE + +_PROTEIN_SYSTEMS = {SystemType.PROTEIN_ONLY, SystemType.PROTEIN_PROTEIN, + SystemType.PROTEIN_PEPTIDE, SystemType.PROTEIN_LIGAND, + SystemType.PROTEIN_DNA, SystemType.PROTEIN_RNA, + SystemType.PROTEIN_MEMBRANE, SystemType.MIXED} + + +class RMSFAnalysis(BaseAnalysis): + name = "rmsf" + label = "Per-residue Cα RMSF" + category = "dynamics" + required_files = {"trajectory", "topology"} + supported_systems = _PROTEIN_SYSTEMS + outputs = ["results/rmsf.csv", "figures/rmsf.png"] + default_params = {"ref_frame": 0} + + def run(self, ctx) -> dict: + p = self.params(ctx) + plotting.set_style() + u = ctx.core_universe() + sel = "protein and name CA" + avg = align.AverageStructure(u, u, select=sel, ref_frame=p["ref_frame"]).run() + align.AlignTraj(u, avg.results.universe, select=sel, in_memory=True).run() + ca = u.select_atoms(sel) + rmsf = rms.RMSF(ca).run().results.rmsf / 10.0 # nm + resids, resnames = ca.resids, ca.resnames + + df = pd.DataFrame({"resid": resids, "resname": resnames, "rmsf_nm": rmsf}) + ctx.write_csv(df, "rmsf.csv") + thr = rmsf.mean() + rmsf.std() + flex = df[df["rmsf_nm"] > thr] + + fig, ax = plotting.new_axes(figsize=(8.2, 4.6)) + ax.plot(resids, rmsf, color=PALETTE["primary"], lw=1.8) + ax.fill_between(resids, 0, rmsf, color=PALETTE["primary"], alpha=0.15) + ax.axhline(thr, ls="--", lw=1.1, color=PALETTE["accent"], + label=f"mean + 1$\\sigma$ = {thr:.3f} nm") + ax.scatter(flex["resid"], flex["rmsf_nm"], color=PALETTE["secondary"], + s=28, zorder=5, label="Flexible") + ax.set_xlabel("Residue number"); ax.set_ylabel("RMSF (nm)") + ax.set_title(self.label); ax.legend(loc="upper right") + ax.margins(x=0.01); ax.set_ylim(bottom=0) + plotting.save_figure(fig, ctx.fig_path("rmsf"), dpi=ctx.config.dpi) + + return {"rmsf": st.describe(rmsf, "rmsf_nm"), + "n_flexible": int(len(flex)), + "flexible_resids": flex["resid"].tolist(), + "most_flexible_resid": int(resids[np.argmax(rmsf)]), + "figure": "rmsf"} diff --git a/mdforge/analysis/rog.py b/mdforge/analysis/rog.py new file mode 100644 index 0000000..87d7e3a --- /dev/null +++ b/mdforge/analysis/rog.py @@ -0,0 +1,55 @@ +"""Radius of gyration over time (universal solute compactness).""" + +from __future__ import annotations + +import numpy as np +import pandas as pd + +from mdforge.core.base import BaseAnalysis +from mdforge import plotting +from mdforge import statistics as st +from mdforge.plotting import PALETTE + + +class RadiusOfGyration(BaseAnalysis): + name = "rog" + label = "Radius of gyration" + category = "structure" + required_files = {"trajectory", "topology"} + supported_systems = {"*"} + outputs = ["results/rog.csv", "figures/rog.png"] + default_params = {"window": 20} + + def run(self, ctx) -> dict: + p = self.params(ctx) + plotting.set_style() + u = ctx.core_universe() + sel = "protein" if ctx.system.flags.get("has_protein", True) else "all" + solute = u.select_atoms(sel) if sel != "all" else u.atoms + + rg = np.empty(len(u.trajectory)); times = np.empty_like(rg) + for i, ts in enumerate(ctx.iter_frames(u, desc="[rog]")): + rg[i] = solute.radius_of_gyration() / 10.0 + times[i] = ts.time / 1000.0 + + ma = st.moving_average(rg, p["window"]) + ctx.write_csv(pd.DataFrame({"time_ns": times, "rg_nm": rg, + "rg_movavg_nm": ma}), "rog.csv") + mean, std = float(rg.mean()), float(rg.std(ddof=1)) + + fig, ax = plotting.new_axes() + ax.plot(times, rg, color=PALETTE["primary"], alpha=0.35, lw=1.0) + ax.plot(times, ma, color=PALETTE["primary"], lw=2.2, label="MA") + ax.axhline(mean, ls="--", lw=1.2, color=PALETTE["accent"], + label=f"mean = {mean:.3f} nm") + ax.fill_between(times, mean - std, mean + std, color=PALETTE["muted"], + alpha=0.25, label=f"$\\pm\\sigma$ = {std:.3f} nm") + ax.set_xlabel("Time (ns)"); ax.set_ylabel(r"$R_g$ (nm)") + ax.set_title(self.label); ax.legend(loc="best"); ax.margins(x=0.01) + plotting.save_figure(fig, ctx.fig_path("rog"), dpi=ctx.config.dpi) + + return {"rg": st.describe(rg, "rg_nm"), + "convergence": st.plateau_detection(times, rg), + "rg_initial_nm": float(rg[:min(20, len(rg))].mean()), + "rg_final_nm": float(rg[len(rg)//2:].mean()), + "figure": "rog"} diff --git a/mdforge/analysis/sasa.py b/mdforge/analysis/sasa.py new file mode 100644 index 0000000..8470520 --- /dev/null +++ b/mdforge/analysis/sasa.py @@ -0,0 +1,64 @@ +"""Solvent-accessible surface area (total + per-residue, Shrake-Rupley/mdtraj).""" + +from __future__ import annotations + +import numpy as np +import pandas as pd +import mdtraj as md + +from mdforge.core.base import BaseAnalysis +from mdforge import plotting +from mdforge import statistics as st +from mdforge.plotting import PALETTE + + +class SASAAnalysis(BaseAnalysis): + name = "sasa" + label = "Solvent-accessible surface area" + category = "structure" + required_files = {"trajectory", "topology"} + supported_systems = {"*"} + outputs = ["results/sasa_total.csv", "results/sasa_per_residue.csv", + "figures/sasa_timeseries.png", "figures/sasa_per_residue.png"] + default_params = {"window": 20} + + def run(self, ctx) -> dict: + p = self.params(ctx) + plotting.set_style() + # mdtraj reads the cached solute trajectory + ctx.core_universe() # ensure core.{pdb,xtc} exist + traj = md.load(str(ctx.config.data_dir / "core.xtc"), + top=str(ctx.config.data_dir / "core.pdb")) + time_ns = traj.time / 1000.0 + sasa_res = md.shrake_rupley(traj, mode="residue") # nm^2 + total = sasa_res.sum(axis=1) + resids = np.array([r.resSeq for r in traj.topology.residues]) + resnames = np.array([r.name for r in traj.topology.residues]) + + ma = st.moving_average(total, p["window"]) + ctx.write_csv(pd.DataFrame({"time_ns": time_ns, "sasa_total_nm2": total, + "sasa_total_movavg_nm2": ma}), "sasa_total.csv") + ctx.write_csv(pd.DataFrame({"resid": resids, "resname": resnames, + "sasa_mean_nm2": sasa_res.mean(axis=0), + "sasa_std_nm2": sasa_res.std(axis=0)}), + "sasa_per_residue.csv") + mean = float(total.mean()) + + fig, ax = plotting.new_axes() + ax.plot(time_ns, total, color=PALETTE["primary"], alpha=0.35, lw=1.0) + ax.plot(time_ns, ma, color=PALETTE["primary"], lw=2.2, label="MA") + ax.axhline(mean, ls="--", lw=1.2, color=PALETTE["accent"], + label=f"mean = {mean:.1f} nm$^2$") + ax.set_xlabel("Time (ns)"); ax.set_ylabel("Total SASA (nm$^2$)") + ax.set_title(self.label); ax.legend(loc="best"); ax.margins(x=0.01) + plotting.save_figure(fig, ctx.fig_path("sasa_timeseries"), dpi=ctx.config.dpi) + + fig, ax = plotting.new_axes(figsize=(8.2, 4.6)) + ax.bar(resids, sasa_res.mean(axis=0), color=PALETTE["primary"], width=1.0) + ax.set_xlabel("Residue number"); ax.set_ylabel("Mean SASA (nm$^2$)") + ax.set_title("Per-residue solvent accessibility"); ax.margins(x=0.01) + plotting.save_figure(fig, ctx.fig_path("sasa_per_residue"), dpi=ctx.config.dpi) + + return {"total_sasa": st.describe(total, "sasa_total_nm2"), + "convergence": st.plateau_detection(time_ns, total), + "figure": "sasa_timeseries"} diff --git a/mdforge/cli/__init__.py b/mdforge/cli/__init__.py new file mode 100644 index 0000000..dd4c185 --- /dev/null +++ b/mdforge/cli/__init__.py @@ -0,0 +1 @@ +"""Command-line interface for mdforge.""" diff --git a/mdforge/cli/interactive.py b/mdforge/cli/interactive.py new file mode 100644 index 0000000..a5f5604 --- /dev/null +++ b/mdforge/cli/interactive.py @@ -0,0 +1,23 @@ +"""Minimal interactive prompt helpers (only used with --interactive).""" + +from __future__ import annotations + + +def ask_choice(question: str, options: list[str]) -> int: + """Ask the user to pick one option; returns its index (0-based).""" + print(f"\n{question}") + for i, opt in enumerate(options, 1): + print(f" [{i}] {opt}") + while True: + raw = input("Enter choice number: ").strip() + if raw.isdigit() and 1 <= int(raw) <= len(options): + return int(raw) - 1 + print(" invalid choice, try again.") + + +def confirm(question: str, default: bool = True) -> bool: + suffix = "[Y/n]" if default else "[y/N]" + raw = input(f"{question} {suffix} ").strip().lower() + if not raw: + return default + return raw.startswith("y") diff --git a/mdforge/cli/main.py b/mdforge/cli/main.py new file mode 100644 index 0000000..3974cf0 --- /dev/null +++ b/mdforge/cli/main.py @@ -0,0 +1,142 @@ +""" +mdforge command-line entry point. + +Subcommands +----------- + mdforge analyze --input DIR [--output DIR] [options] run the pipeline + mdforge detect --input DIR detect system only + mdforge list-analyses [--system TYPE] list available analyses + mdforge version +""" + +from __future__ import annotations + +import argparse +import sys + +from mdforge import __version__ + + +def _add_common(p): + p.add_argument("--input", "-i", help="Simulation directory (searched recursively).") + p.add_argument("--output", "-o", help="Output directory.") + p.add_argument("--config", help="YAML config file (CLI flags override it).") + + +def build_parser() -> argparse.ArgumentParser: + parser = argparse.ArgumentParser(prog="mdforge", description=__doc__, + formatter_class=argparse.RawDescriptionHelpFormatter) + parser.add_argument("--version", action="version", version=f"mdforge {__version__}") + sub = parser.add_subparsers(dest="command") + + a = sub.add_parser("analyze", help="Run the analysis pipeline.") + _add_common(a) + a.add_argument("--traj", help="Override the discovered trajectory.") + a.add_argument("--top", help="Override the discovered topology.") + a.add_argument("--edr", help="Override the discovered energy file.") + a.add_argument("--ndx", help="Override the discovered index file.") + a.add_argument("--system-type", dest="system_type", help="Override auto-detected system type.") + a.add_argument("--chains", help="Chain selection, e.g. 'A' or 'A+B'.") + a.add_argument("--ligand", help="Ligand resname override.") + a.add_argument("--analyses", help="Comma-separated analyses to run (default: auto).") + a.add_argument("--exclude", help="Comma-separated analyses to skip.") + a.add_argument("--all", action="store_true", help="Run every applicable analysis.") + a.add_argument("--stride", type=int, help="Use every Nth frame.") + a.add_argument("--start", type=int, help="First frame.") + a.add_argument("--end", type=int, help="Last frame.") + a.add_argument("--threads", type=int, help="Worker threads (where supported).") + a.add_argument("--seed", type=int, help="Random seed.") + a.add_argument("--report", help="Report formats, comma-separated (md,html,pdf).") + a.add_argument("--plugin-dir", help="Extra directory of plugin modules.") + a.add_argument("--interactive", action="store_true", help="Ask questions when ambiguous.") + a.add_argument("--plan", action="store_true", help="Dry-run: show what would run and why.") + + d = sub.add_parser("detect", help="Detect and print the system composition.") + _add_common(d) + + ls = sub.add_parser("list-analyses", help="List registered analyses.") + ls.add_argument("--system", help="Filter to a system type.") + + sub.add_parser("version", help="Print version.") + return parser + + +def _cmd_analyze(args) -> int: + from mdforge.core.config import RunConfig + from mdforge.core.pipeline import run_pipeline + from mdforge.cli.interactive import ask_choice + 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 cfg.output_dir is None: + cfg.output_dir = cfg.input_dir / "mdforge_results" + ask = ask_choice if args.interactive else None + run_pipeline(cfg, ask=ask, plan_only=args.plan) + return 0 + + +def _cmd_detect(args) -> int: + from mdforge.io.discovery import discover_files + from mdforge.io.validation import validate_fileset + from mdforge.core.system import detect_system + if not args.input: + print("error: --input is required.") + return 2 + fs = discover_files(args.input) + val = validate_fileset(fs) + print("Discovered files:") + for k, v in fs.to_dict().items(): + if v: + print(f" {k:11s}: {v}") + print("\nValidation:\n" + val.report()) + if val.ok: + topo = str(fs.topology or fs.structure) + system = detect_system(topo) + print("\n" + system.summary()) + return 0 if val.ok else 2 + + +def _cmd_list(args) -> int: + from mdforge.core.registry import registry + from mdforge.core.system import SystemType + registry.ensure_builtins_loaded() + filt = SystemType(args.system) if args.system else None + print(f"{'name':22s} {'category':13s} systems") + print("-" * 70) + for name in registry.names(): + cls = registry.get(name) + if filt and "*" not in cls.supported_systems and filt not in cls.supported_systems: + continue + systems = "all" if "*" in cls.supported_systems else \ + ", ".join(sorted(s.value for s in cls.supported_systems)) + print(f"{cls.name:22s} {cls.category:13s} {systems}") + return 0 + + +def _force_utf8() -> None: + """Make console output robust to non-ASCII (e.g. 'Cα') on Windows cp1252.""" + import sys + for stream in (sys.stdout, sys.stderr): + try: + stream.reconfigure(encoding="utf-8", errors="replace") + except Exception: + pass + + +def main(argv=None) -> int: + _force_utf8() + parser = build_parser() + args = parser.parse_args(argv) + if not args.command or args.command == "version": + if not args.command: + parser.print_help() + return 1 + print(f"mdforge {__version__}") + return 0 + return {"analyze": _cmd_analyze, "detect": _cmd_detect, + "list-analyses": _cmd_list}[args.command](args) + + +if __name__ == "__main__": + sys.exit(main()) diff --git a/mdforge/core/__init__.py b/mdforge/core/__init__.py new file mode 100644 index 0000000..84ba6fb --- /dev/null +++ b/mdforge/core/__init__.py @@ -0,0 +1 @@ +"""Core framework: system detection, analysis base class, registry, config, provenance, pipeline.""" diff --git a/mdforge/core/base.py b/mdforge/core/base.py new file mode 100644 index 0000000..0700567 --- /dev/null +++ b/mdforge/core/base.py @@ -0,0 +1,92 @@ +""" +BaseAnalysis -- the contract every analysis (built-in or plugin) implements. + +Each analysis declares, as class attributes: + +* ``name`` -- unique short key (used in the CLI and config). +* ``label`` -- human-readable title. +* ``category`` -- grouping ("structure", "dynamics", "interactions"...). +* ``required_files`` -- semantic input keys it needs ({"trajectory","topology", + "energy",...}); the pipeline verifies availability. +* ``supported_systems`` -- set of SystemType it applies to, or {"*"} for all. +* ``outputs`` -- declared output basenames (for documentation/manifest). +* ``default_params`` -- parameters with defaults (overridable via config/CLI). + +Subclasses are **auto-registered** on definition, so adding an analysis (or a +plugin) is simply a matter of subclassing ``BaseAnalysis`` and being imported. +""" + +from __future__ import annotations + +import abc +from typing import TYPE_CHECKING + +from mdforge.core.system import SystemType + +if TYPE_CHECKING: # avoid import cycle at runtime + from mdforge.core.context import AnalysisContext + from mdforge.core.system import SystemInfo + + +class BaseAnalysis(abc.ABC): + # -- declaration (override in subclasses) ----------------------------- # + name: str = "" + label: str = "" + category: str = "general" + required_files: set[str] = {"trajectory", "topology"} + supported_systems: set = {"*"} + outputs: list[str] = [] + default_params: dict = {} + + # -- auto-registration ------------------------------------------------ # + def __init_subclass__(cls, register: bool = True, **kwargs): + super().__init_subclass__(**kwargs) + # Only register concrete, named analyses (skip intermediate bases). + if register and getattr(cls, "name", ""): + from mdforge.core.registry import registry + registry.register(cls) + + # -- applicability ---------------------------------------------------- # + @classmethod + def is_applicable(cls, system: "SystemInfo") -> bool: + """Whether this analysis applies to a detected system.""" + if "*" in cls.supported_systems: + return True + return system.system_type in cls.supported_systems + + @classmethod + def missing_files(cls, available: set[str]) -> set[str]: + """Required semantic file keys that are not available.""" + return set(cls.required_files) - set(available) + + # -- execution -------------------------------------------------------- # + @abc.abstractmethod + def run(self, ctx: "AnalysisContext") -> dict: + """ + Execute the analysis and return a JSON-serialisable summary dict. + + Implementations should read inputs via ``ctx`` (universe, selections, + params), write CSV/figure outputs under ``ctx.results_dir`` / + ``ctx.figures_dir``, and return a compact summary used by the report. + """ + raise NotImplementedError + + # -- helpers ---------------------------------------------------------- # + def params(self, ctx: "AnalysisContext") -> dict: + """Merge default_params with any user overrides for this analysis.""" + merged = dict(self.default_params) + merged.update(ctx.params_for(self.name)) + return merged + + def describe(self) -> dict: + return { + "name": self.name, + "label": self.label, + "category": self.category, + "required_files": sorted(self.required_files), + "supported_systems": ["*"] if "*" in self.supported_systems + else sorted(s.value if isinstance(s, SystemType) else s + for s in self.supported_systems), + "outputs": list(self.outputs), + "default_params": dict(self.default_params), + } diff --git a/mdforge/core/config.py b/mdforge/core/config.py new file mode 100644 index 0000000..b1ac1fc --- /dev/null +++ b/mdforge/core/config.py @@ -0,0 +1,119 @@ +""" +Run configuration: merges command-line arguments with an optional YAML file. + +Precedence (highest first): explicit CLI flags > YAML file > built-in defaults. +A single ``config.yaml`` can therefore describe an entire reproducible workflow +(inputs, selections, per-analysis parameters, report formats) and be re-run with +``mdforge analyze --config config.yaml``. +""" + +from __future__ import annotations + +from dataclasses import dataclass, field, asdict +from pathlib import Path + + +@dataclass +class RunConfig: + # -- inputs / outputs ------------------------------------------------- # + input_dir: Path | None = None + output_dir: Path = Path("mdforge_results") + # explicit overrides for auto-discovered files (optional) + trajectory: Path | None = None + topology: Path | None = None + energy: Path | None = None + index: Path | None = None + gmx_top: Path | None = None + + # -- system / selections --------------------------------------------- # + system_type: str | None = None # override auto-detection + chains: str | None = None # e.g. "A", "A+B" + ligand: str | None = None # ligand resname override + + # -- trajectory slicing / performance -------------------------------- # + start: int = 0 + end: int | None = None + stride: int = 1 + threads: int = 1 + core_extract: bool = True # cache a solute-only trajectory + + # -- analysis selection ---------------------------------------------- # + analyses: list[str] | None = None # None => auto-select for system + run_all: bool = False + exclude: list[str] = field(default_factory=list) + plugin_dirs: list[str] = field(default_factory=list) + params: dict = field(default_factory=dict) # {analysis_name: {param: value}} + + # -- reporting / repro ----------------------------------------------- # + report_formats: list[str] = field(default_factory=lambda: ["md", "html"]) + dpi: int = 300 + seed: int = 0 + interactive: bool = False + + def __post_init__(self): + for k in ("input_dir", "output_dir", "trajectory", "topology", + "energy", "index", "gmx_top"): + v = getattr(self, k) + if v is not None: + setattr(self, k, Path(v)) + + # -- output sub-directories ------------------------------------------ # + @property + def data_dir(self) -> Path: return self.output_dir / "data" + @property + def results_dir(self) -> Path: return self.output_dir / "results" + @property + def figures_dir(self) -> Path: return self.output_dir / "figures" + @property + def tables_dir(self) -> Path: return self.output_dir / "tables" + @property + def report_dir(self) -> Path: return self.output_dir / "report" + + def ensure_dirs(self) -> None: + for d in (self.output_dir, self.data_dir, self.results_dir, + self.figures_dir, self.tables_dir, self.report_dir): + d.mkdir(parents=True, exist_ok=True) + + def params_for(self, name: str) -> dict: + return dict(self.params.get(name, {})) + + def to_dict(self) -> dict: + d = asdict(self) + for k, v in d.items(): + if isinstance(v, Path): + d[k] = str(v) + return d + + # -- construction ---------------------------------------------------- # + @classmethod + def from_yaml(cls, path: str | Path) -> "RunConfig": + import yaml + data = yaml.safe_load(Path(path).read_text()) or {} + known = {f for f in cls.__dataclass_fields__} + return cls(**{k: v for k, v in data.items() if k in known}) + + @classmethod + def from_args(cls, args, yaml_path: str | Path | None = None) -> "RunConfig": + """Build from an argparse namespace, layering YAML underneath if given.""" + base = cls.from_yaml(yaml_path) if yaml_path else cls() + mapping = { + "input_dir": "input", "output_dir": "output", "trajectory": "traj", + "topology": "top", "energy": "edr", "index": "ndx", + "system_type": "system_type", "chains": "chains", "ligand": "ligand", + "start": "start", "end": "end", "stride": "stride", "threads": "threads", + "run_all": "all", "interactive": "interactive", "seed": "seed", + } + for attr, argname in mapping.items(): + val = getattr(args, argname, None) + if val is not None: + setattr(base, attr, val) + if getattr(args, "analyses", None): + base.analyses = [a.strip() for a in args.analyses.split(",")] + if getattr(args, "exclude", None): + base.exclude = [a.strip() for a in args.exclude.split(",")] + if getattr(args, "plugin_dir", None): + base.plugin_dirs = list(base.plugin_dirs) + [args.plugin_dir] + if getattr(args, "report", None): + base.report_formats = [f.strip() for f in args.report.split(",")] + base.__post_init__() + return base diff --git a/mdforge/core/context.py b/mdforge/core/context.py new file mode 100644 index 0000000..be98cdb --- /dev/null +++ b/mdforge/core/context.py @@ -0,0 +1,135 @@ +""" +AnalysisContext -- the runtime object every analysis receives. + +Provides universe access (full and a cached solute-only trajectory for speed), +frame slicing (start/end/stride), the detected system's selection strings, +per-analysis parameters, RNG seeding, and output-path helpers. + +Performance note +---------------- +``core_universe()`` streams the full trajectory **once** and writes a +PBC-corrected, solute-only trajectory (protein/nucleic/ligand/cofactor, no +water/ions) that all structural analyses then reuse. For explicit-solvent +systems this is typically 1-2 orders of magnitude smaller than the full +trajectory, which is what keeps the toolkit tractable on 100 GB+ inputs. +""" + +from __future__ import annotations + +import warnings +from pathlib import Path + +import numpy as np +import pandas as pd + +warnings.filterwarnings("ignore", category=DeprecationWarning, module="MDAnalysis") +import MDAnalysis as mda # noqa: E402 +from MDAnalysis.transformations import unwrap # noqa: E402 +from tqdm import tqdm # noqa: E402 + +from mdforge.core.config import RunConfig +from mdforge.core.system import SystemInfo, ComponentType + + +class AnalysisContext: + def __init__(self, config: RunConfig, system: SystemInfo, fileset, + provenance=None): + self.config = config + self.system = system + self.fileset = fileset + self.provenance = provenance + self.rng = np.random.default_rng(config.seed) + self._full: mda.Universe | None = None + self._core: mda.Universe | None = None + config.ensure_dirs() + + # -- universes -------------------------------------------------------- # + @property + def topology(self) -> str: + return str(self.config.topology or self.fileset.topology or self.fileset.structure) + + @property + def trajectory(self) -> str: + return str(self.config.trajectory or self.fileset.trajectory) + + def full_universe(self) -> mda.Universe: + """The full solvated system (loaded lazily, streamed frame-by-frame).""" + if self._full is None: + self._full = mda.Universe(self.topology, self.trajectory) + return self._full + + def core_selection(self) -> str: + """Selection string for the 'solute' of interest (no water/ions).""" + parts = [] + for ct in (ComponentType.PROTEIN, ComponentType.DNA, ComponentType.RNA, + ComponentType.LIGAND, ComponentType.COFACTOR): + if self.system.has(ct): + parts.append(f"({self.system.components[ct.value].selection})") + if not parts: # e.g. pure membrane + return "not (resname SOL WAT HOH TIP3 SPC NA CL K POT CLA)" + return " or ".join(parts) + + def core_universe(self) -> mda.Universe: + """ + A cached, PBC-corrected, solute-only trajectory universe. + + Built once and stored under ``data/core.{pdb,xtc}``; reused on subsequent + analyses and subsequent runs. + """ + if self._core is not None: + return self._core + core_pdb = self.config.data_dir / "core.pdb" + core_xtc = self.config.data_dir / "core.xtc" + if not (core_pdb.exists() and core_xtc.exists()): + self._extract_core(core_pdb, core_xtc) + self._core = mda.Universe(str(core_pdb), str(core_xtc)) + return self._core + + def _extract_core(self, core_pdb: Path, core_xtc: Path) -> None: + u = self.full_universe() + core = u.select_atoms(self.core_selection()) + if core.n_atoms == 0: + raise RuntimeError(f"Core selection matched no atoms: {self.core_selection()!r}") + # make the solute whole across PBC (needs bonds; .tpr provides them) + try: + u.trajectory.add_transformations(unwrap(core)) + except Exception: + pass # some topologies lack bonds; proceed unwrapped + sl = self.frame_slice() + u.trajectory[sl.start or 0] + core.write(str(core_pdb)) + n = len(range(*sl.indices(len(u.trajectory)))) + with mda.Writer(str(core_xtc), core.n_atoms) as W: + for _ in tqdm(u.trajectory[sl], total=n, + desc="[core-extract] solute", unit="frame"): + W.write(core) + + # -- frames / time ---------------------------------------------------- # + def frame_slice(self) -> slice: + return slice(self.config.start, self.config.end, self.config.stride) + + def times_ns(self, universe: mda.Universe) -> np.ndarray: + return np.array([ts.time for ts in universe.trajectory]) / 1000.0 + + def iter_frames(self, universe: mda.Universe, desc: str = "frames", + stride: int = 1): + n = len(universe.trajectory[::stride]) + return tqdm(universe.trajectory[::stride], total=n, desc=desc, unit="frame") + + # -- selections / params --------------------------------------------- # + def selection(self, key: str, default: str = "protein") -> str: + return self.system.selections.get(key, default) + + def params_for(self, name: str) -> dict: + return self.config.params_for(name) + + # -- output paths ----------------------------------------------------- # + def csv_path(self, name: str) -> Path: return self.config.results_dir / name + def fig_path(self, name: str) -> Path: return self.config.figures_dir / name + def table_path(self, name: str) -> Path: return self.config.tables_dir / name + + def write_csv(self, df: pd.DataFrame, name: str, index: bool = False) -> Path: + p = self.csv_path(name) + p.parent.mkdir(parents=True, exist_ok=True) + df.to_csv(p, index=index) + return p diff --git a/mdforge/core/pipeline.py b/mdforge/core/pipeline.py new file mode 100644 index 0000000..2414adf --- /dev/null +++ b/mdforge/core/pipeline.py @@ -0,0 +1,130 @@ +""" +End-to-end pipeline: discover -> validate -> detect -> select -> run -> report. + +Each analysis is isolated: a failure is recorded in the manifest and the pipeline +continues, so one broken module never aborts the whole run. +""" + +from __future__ import annotations + +import time +import traceback +from pathlib import Path + +from mdforge.core.config import RunConfig +from mdforge.core.context import AnalysisContext +from mdforge.core.provenance import RunManifest +from mdforge.core.registry import registry +from mdforge.core.system import detect_system, SystemType +from mdforge.io.discovery import discover_files, pick_trajectory +from mdforge.io.validation import validate_fileset + + +def build_plan(config: RunConfig, ask=None): + """Resolve files + system + selected analyses without running anything.""" + fs = discover_files(config.input_dir) + fs.trajectory = pick_trajectory(fs, interactive=config.interactive, ask=ask) + if config.trajectory: + fs.trajectory = Path(config.trajectory) + if config.topology: + fs.topology = Path(config.topology) + validation = validate_fileset(fs) + + system = None + selected, skipped = [], [] + if validation.ok: + topo = str(config.topology or fs.topology or fs.structure) + system = detect_system(topo) + if config.system_type: + system.system_type = SystemType(config.system_type) + registry.ensure_builtins_loaded() + registry.discover_plugins(config.plugin_dirs) + requested = config.analyses + selected, skipped = registry.select( + system, fs.available_keys(), requested=requested, run_all=config.run_all) + if config.exclude: + selected = [c for c in selected if c.name not in config.exclude] + return fs, validation, system, selected, skipped + + +def run_pipeline(config: RunConfig | None = None, *, ask=None, + plan_only: bool = False, log=print, **kwargs): + """ + Run the full pipeline. + + Parameters + ---------- + config : a RunConfig, or None to build one from kwargs (e.g. input_dir=...). + ask : optional callable(question, options)->index for interactive prompts. + plan_only : if True, resolve and return the plan without executing. + """ + if config is None: + config = RunConfig(**kwargs) + + fs, validation, system, selected, skipped = build_plan(config, ask=ask) + + log("=" * 74) + log(" mdforge -- MD analysis pipeline") + log("=" * 74) + log(f" input : {config.input_dir}") + log(f" output: {config.output_dir}") + if not validation.ok: + log("\n[input validation FAILED]\n" + validation.report()) + raise SystemExit(2) + for w in validation.warnings: + log(f" [warn] {w}") + log("\n" + system.summary()) + log(f"\n Selected {len(selected)} analyses: " + f"{', '.join(c.name for c in selected)}") + if skipped: + log(" Skipped:") + for name, reason in skipped: + log(f" - {name}: {reason}") + + if plan_only: + return {"fileset": fs, "system": system, + "selected": [c.name for c in selected], "skipped": skipped} + + config.ensure_dirs() + manifest = RunManifest() + manifest.set_config(config) + manifest.set_system(system) + manifest.set_inputs(fs) + + ctx = AnalysisContext(config, system, fs, provenance=manifest) + results: dict[str, dict] = {} + for cls in selected: + log("-" * 74) + log(f">>> {cls.name}: {cls.label}") + t0 = time.time() + try: + summary = cls().run(ctx) + dt = time.time() - t0 + results[cls.name] = summary + manifest.record_analysis(cls.name, "ok", dt, summary=summary) + log(f"[ok] {cls.name} in {dt:.1f}s") + except Exception: + dt = time.time() - t0 + err = traceback.format_exc() + manifest.record_analysis(cls.name, "error", dt, + error=err.splitlines()[-1]) + results[cls.name] = {"error": err.splitlines()[-1]} + log(f"[ERROR] {cls.name}:\n{err}") + + manifest_paths = manifest.write(config.output_dir) + log(f"\n[manifest] {manifest_paths[0]}") + + # -- report ----------------------------------------------------------- # + try: + from mdforge.report.generator import generate_report + report_paths = generate_report(ctx, results, manifest, + formats=config.report_formats) + for p in report_paths: + log(f"[report] {p}") + except Exception: + log("[report] generation failed:\n" + traceback.format_exc()) + + log("=" * 74) + log(" Done.") + return {"system": system, "results": results, "manifest": manifest.data, + "output_dir": str(config.output_dir)} diff --git a/mdforge/core/provenance.py b/mdforge/core/provenance.py new file mode 100644 index 0000000..f508662 --- /dev/null +++ b/mdforge/core/provenance.py @@ -0,0 +1,128 @@ +""" +Reproducibility manifest. + +Captures everything another researcher needs to reproduce a run: software and +library versions, OS/CPU, the exact command, git commit, input-file fingerprints, +random seeds, all parameters, per-analysis status/runtime, and total wall-time. +Written as both JSON and YAML. +""" + +from __future__ import annotations + +import hashlib +import json +import platform +import subprocess +import sys +import time +from datetime import datetime, timezone +from pathlib import Path + + +def _lib_version(mod: str) -> str | None: + try: + m = __import__(mod) + return getattr(m, "__version__", None) + except Exception: + return None + + +def _git_commit() -> dict: + try: + root = Path(__file__).resolve().parents[2] + commit = subprocess.check_output( + ["git", "-C", str(root), "rev-parse", "HEAD"], + stderr=subprocess.DEVNULL).decode().strip() + dirty = bool(subprocess.check_output( + ["git", "-C", str(root), "status", "--porcelain"], + stderr=subprocess.DEVNULL).decode().strip()) + return {"commit": commit, "dirty": dirty} + except Exception: + return {"commit": None, "dirty": None} + + +def file_fingerprint(path: str | Path, chunk: int = 1 << 20) -> dict: + """ + Cheap but robust fingerprint for possibly-huge trajectories: size + SHA-256 + of the first and last ``chunk`` bytes (avoids hashing 100 GB in full). + """ + path = Path(path) + if not path.exists(): + return {"path": str(path), "exists": False} + size = path.stat().st_size + h = hashlib.sha256() + with open(path, "rb") as fh: + h.update(fh.read(chunk)) + if size > 2 * chunk: + fh.seek(-chunk, 2) + h.update(fh.read(chunk)) + return {"path": str(path), "exists": True, "size_bytes": size, + "sha256_ends": h.hexdigest()} + + +class RunManifest: + """Collects provenance during a pipeline run and writes it out.""" + + def __init__(self) -> None: + self._t0 = time.time() + self.data: dict = { + "mdforge_version": _lib_version("mdforge"), + "timestamp_utc": datetime.now(timezone.utc).isoformat(), + "command": " ".join(sys.argv), + "python": sys.version.split()[0], + "platform": { + "system": platform.system(), "release": platform.release(), + "machine": platform.machine(), "processor": platform.processor(), + }, + "libraries": {m: _lib_version(m) for m in ( + "numpy", "scipy", "pandas", "matplotlib", "MDAnalysis", + "mdtraj", "sklearn", "networkx", "prolif", "yaml")}, + "git": _git_commit(), + "inputs": {}, + "system": {}, + "config": {}, + "seed": None, + "analyses": {}, + } + + # -- population ------------------------------------------------------- # + def set_config(self, cfg) -> None: + self.data["config"] = cfg.to_dict() + self.data["seed"] = cfg.seed + + def set_system(self, system) -> None: + self.data["system"] = system.to_dict() + + def set_inputs(self, fileset) -> None: + d = fileset.to_dict() + fp = {} + for role in ("trajectory", "topology", "energy"): + p = getattr(fileset, role) + if p is not None: + fp[role] = file_fingerprint(p) + self.data["inputs"] = {"resolved": d, "fingerprints": fp} + + def record_analysis(self, name: str, status: str, runtime_s: float, + summary: dict | None = None, error: str | None = None) -> None: + self.data["analyses"][name] = { + "status": status, "runtime_s": round(runtime_s, 2), + "summary": summary, "error": error, + } + + # -- output ----------------------------------------------------------- # + def write(self, outdir: Path) -> list[Path]: + self.data["total_runtime_s"] = round(time.time() - self._t0, 2) + outdir = Path(outdir) + outdir.mkdir(parents=True, exist_ok=True) + written = [] + jp = outdir / "manifest.json" + jp.write_text(json.dumps(self.data, indent=2, default=str)) + written.append(jp) + try: + import yaml + yp = outdir / "manifest.yaml" + yp.write_text(yaml.safe_dump(self.data, sort_keys=False)) + written.append(yp) + except Exception: + pass + return written diff --git a/mdforge/core/registry.py b/mdforge/core/registry.py new file mode 100644 index 0000000..b938628 --- /dev/null +++ b/mdforge/core/registry.py @@ -0,0 +1,140 @@ +""" +Analysis registry + plugin discovery. + +The global ``registry`` collects every :class:`BaseAnalysis` subclass (built-in +or plugin). It supports: + +* automatic selection of analyses for a detected system, +* discovery of plugins from (a) the built-in ``mdforge.analysis.plugins`` package, + (b) any user directory passed via ``--plugin-dir`` / config, and + (c) installed packages advertising a ``mdforge.plugins`` entry-point. + +Adding an analysis requires no change to the core: define a ``BaseAnalysis`` +subclass and make sure its module is imported. +""" + +from __future__ import annotations + +import importlib +import importlib.util +import pkgutil +import sys +from pathlib import Path +from typing import TYPE_CHECKING + +if TYPE_CHECKING: + from mdforge.core.base import BaseAnalysis + from mdforge.core.system import SystemInfo + + +class AnalysisRegistry: + def __init__(self) -> None: + self._analyses: dict[str, type] = {} + + # -- registration ----------------------------------------------------- # + def register(self, cls: type) -> type: + key = getattr(cls, "name", "") + if not key: + raise ValueError(f"Analysis {cls!r} has no 'name'.") + self._analyses[key] = cls + return cls + + def get(self, name: str) -> type: + if name not in self._analyses: + raise KeyError(f"Unknown analysis '{name}'. " + f"Available: {', '.join(sorted(self._analyses))}") + return self._analyses[name] + + def all(self) -> dict[str, type]: + return dict(self._analyses) + + def names(self) -> list[str]: + return sorted(self._analyses) + + # -- selection -------------------------------------------------------- # + def select(self, system: "SystemInfo", available_files: set[str], + requested: list[str] | None = None, + run_all: bool = False) -> tuple[list[type], list[tuple[str, str]]]: + """ + Choose analyses for a system. + + Returns ``(selected, skipped)`` where ``skipped`` is a list of + ``(name, reason)`` explaining why an analysis was not selected -- this + powers the ``--plan`` dry-run. + """ + self.ensure_builtins_loaded() + selected, skipped = [], [] + pool = self._analyses if requested is None else { + n: self._analyses[n] for n in requested if n in self._analyses} + for name in (requested or sorted(self._analyses)): + cls = self._analyses.get(name) + if cls is None: + skipped.append((name, "not registered")) + continue + if not run_all and requested is None and not cls.is_applicable(system): + skipped.append((name, f"not applicable to {system.system_type.value}")) + continue + missing = cls.missing_files(available_files) + if missing: + skipped.append((name, f"missing required files: {', '.join(sorted(missing))}")) + continue + selected.append(cls) + return selected, skipped + + # -- plugin / builtin loading ---------------------------------------- # + _builtins_loaded = False + + def ensure_builtins_loaded(self) -> None: + if self._builtins_loaded: + return + import mdforge.analysis as analysis_pkg + for mod in pkgutil.iter_modules(analysis_pkg.__path__): + if mod.name.startswith("_") or mod.name == "plugins": + continue + importlib.import_module(f"mdforge.analysis.{mod.name}") + self.discover_plugins() # built-in plugins package + self._builtins_loaded = True + + def discover_plugins(self, plugin_dirs: list[str | Path] | None = None) -> list[str]: + """Import plugin modules from the plugins package, dirs, and entry-points.""" + loaded: list[str] = [] + # (a) built-in plugins package + try: + import mdforge.analysis.plugins as plug_pkg + for mod in pkgutil.iter_modules(plug_pkg.__path__): + if not mod.name.startswith("_"): + importlib.import_module(f"mdforge.analysis.plugins.{mod.name}") + loaded.append(mod.name) + except Exception: + pass + # (b) user directories + for d in (plugin_dirs or []): + d = Path(d) + if not d.is_dir(): + continue + for pyfile in sorted(d.glob("*.py")): + if pyfile.name.startswith("_"): + continue + spec = importlib.util.spec_from_file_location( + f"mdforge_plugin_{pyfile.stem}", pyfile) + if spec and spec.loader: + module = importlib.util.module_from_spec(spec) + sys.modules[spec.name] = module + spec.loader.exec_module(module) + loaded.append(pyfile.stem) + # (c) installed entry-points + try: + from importlib.metadata import entry_points + eps = entry_points() + group = eps.select(group="mdforge.plugins") if hasattr(eps, "select") \ + else eps.get("mdforge.plugins", []) + for ep in group: + ep.load() + loaded.append(ep.name) + except Exception: + pass + return loaded + + +# module-level singleton +registry = AnalysisRegistry() diff --git a/mdforge/core/system.py b/mdforge/core/system.py new file mode 100644 index 0000000..b0f2478 --- /dev/null +++ b/mdforge/core/system.py @@ -0,0 +1,316 @@ +""" +Automatic biomolecular system detection. + +Given a topology (and optionally coordinates), classify every residue into a +component type (protein, DNA, RNA, lipid, water, ion, ligand, cofactor, ...), +count chains, and infer the overall :class:`SystemType`. The result carries +ready-to-use MDAnalysis selection strings so downstream analyses never have to +re-derive selections. + +The detector is deliberately conservative and data-driven: anything that is not +a recognised polymer/solvent/ion is treated as a (potentially interesting) +"ligand/other" component and surfaced to the user rather than silently ignored. +""" + +from __future__ import annotations + +import warnings +from dataclasses import dataclass, field +from enum import Enum +from pathlib import Path + +warnings.filterwarnings("ignore", category=DeprecationWarning, module="MDAnalysis") +import MDAnalysis as mda # noqa: E402 + + +# --------------------------------------------------------------------------- # +# Residue-name dictionaries (extend freely; matching is case-insensitive) +# --------------------------------------------------------------------------- # +_AMINO = { + "ALA", "ARG", "ASN", "ASP", "CYS", "GLN", "GLU", "GLY", "HIS", "ILE", + "LEU", "LYS", "MET", "PHE", "PRO", "SER", "THR", "TRP", "TYR", "VAL", + # protonation / naming variants (CHARMM / AMBER) + "HSD", "HSE", "HSP", "HID", "HIE", "HIP", "CYX", "CYM", "LYN", "ASH", + "GLH", "MSE", "SEC", "PYL", "ACE", "NME", "NMA", "HISE", "HISD", +} +_DNA = {"DA", "DT", "DG", "DC", "DA5", "DA3", "DT5", "DT3", "DG5", "DG3", + "DC5", "DC3", "DI", "THY", "ADE", "GUA", "CYT"} +_RNA = {"A", "U", "G", "C", "RA", "RU", "RG", "RC", "RA5", "RU5", "RG5", "RC5", + "RA3", "RU3", "RG3", "RC3", "URA"} +_WATER = {"SOL", "WAT", "HOH", "TIP3", "TIP4", "TIP5", "TIP3P", "SPC", "SPCE", + "T3P", "T4P", "OPC", "TP3", "TP4"} +_IONS = {"NA", "CL", "K", "POT", "CLA", "SOD", "MG", "CA", "ZN", "FE", "MN", + "CU", "CO", "NI", "LI", "RB", "CS", "BR", "IOD", "I", "CAL", "MG2", + "ZN2", "NA+", "CL-", "K+", "CA2", "FE2", "FE3", "F", "CD"} +_LIPIDS = {"POPC", "POPE", "POPS", "POPG", "POPI", "POPA", "DOPC", "DOPE", + "DPPC", "DMPC", "DSPC", "DLPC", "DPPE", "DPPG", "CHL1", "CHOL", + "SAPI", "PSM", "DPSM", "LPPC", "PLPC", "PLPE", "OANL", "SDPC", + "POP2", "PIP2", "CER", "DAG", "TAG"} +# A curated subset of biologically important non-polymer molecules that are +# usually *cofactors* rather than drug-like ligands. +_COFACTORS = {"HEM", "HEME", "HEC", "HEB", "NAD", "NAI", "NAP", "NDP", "NAJ", + "FAD", "FMN", "FDA", "ATP", "ADP", "AMP", "GTP", "GDP", "GNP", + "ANP", "SAM", "SAH", "PLP", "TPP", "COA", "COO", "COENZYME", + "BTN", "MTE", "MGD", "F43", "B12", "CLA", "BCL", "PQQ"} +# Common crystallographic additives / buffer that are not the ligand of interest +_BUFFER_LIKE = {"GOL", "EDO", "PEG", "PG4", "SO4", "PO4", "ACT", "DMS", "MPD", + "FMT", "TRS", "BME", "IPA", "MES", "EPE", "CIT", "TLA"} + + +class ComponentType(str, Enum): + PROTEIN = "protein" + DNA = "dna" + RNA = "rna" + LIPID = "lipid" + WATER = "water" + ION = "ion" + COFACTOR = "cofactor" + LIGAND = "ligand" + OTHER = "other" + + +class SystemType(str, Enum): + PROTEIN_ONLY = "protein_only" + PROTEIN_PROTEIN = "protein_protein" + PROTEIN_PEPTIDE = "protein_peptide" + PROTEIN_LIGAND = "protein_ligand" + PROTEIN_DNA = "protein_dna" + PROTEIN_RNA = "protein_rna" + PROTEIN_MEMBRANE = "protein_membrane" + NUCLEIC_ONLY = "nucleic_only" + MEMBRANE_ONLY = "membrane_only" + MIXED = "mixed" + UNKNOWN = "unknown" + + +def _classify_resname(resname: str) -> ComponentType: + r = resname.strip().upper() + if r in _AMINO: + return ComponentType.PROTEIN + if r in _DNA: + return ComponentType.DNA + if r in _RNA: + return ComponentType.RNA + if r in _WATER: + return ComponentType.WATER + if r in _IONS: + return ComponentType.ION + if r in _LIPIDS: + return ComponentType.LIPID + if r in _COFACTORS: + return ComponentType.COFACTOR + return ComponentType.OTHER # candidate ligand / buffer — refined below + + +@dataclass +class Component: + """A distinct molecular component and how to select it.""" + ctype: ComponentType + resnames: list[str] + n_residues: int + n_molecules: int + selection: str + + +@dataclass +class SystemInfo: + """Structured description of a detected simulation system.""" + n_atoms: int + system_type: SystemType + components: dict[str, Component] = field(default_factory=dict) + n_protein_chains: int = 0 + protein_chain_lengths: list[int] = field(default_factory=list) + ligand_resnames: list[str] = field(default_factory=list) + selections: dict[str, str] = field(default_factory=dict) + flags: dict[str, bool] = field(default_factory=dict) + notes: list[str] = field(default_factory=list) + + # -- convenience ---------------------------------------------------- # + def has(self, ctype: ComponentType | str) -> bool: + key = ctype.value if isinstance(ctype, ComponentType) else ctype + return key in self.components and self.components[key].n_residues > 0 + + def to_dict(self) -> dict: + return { + "n_atoms": self.n_atoms, + "system_type": self.system_type.value, + "n_protein_chains": self.n_protein_chains, + "protein_chain_lengths": self.protein_chain_lengths, + "ligand_resnames": self.ligand_resnames, + "components": { + k: {"type": c.ctype.value, "resnames": c.resnames, + "n_residues": c.n_residues, "n_molecules": c.n_molecules} + for k, c in self.components.items() + }, + "selections": self.selections, + "flags": self.flags, + "notes": self.notes, + } + + def summary(self) -> str: + lines = [f"System type : {self.system_type.value}", + f"Atoms : {self.n_atoms:,}", + f"Protein : {self.n_protein_chains} chain(s), " + f"lengths {self.protein_chain_lengths}"] + for key, c in self.components.items(): + if c.n_residues: + lines.append(f" {c.ctype.value:9s}: {c.n_residues} residues " + f"({', '.join(c.resnames[:6])}" + f"{' …' if len(c.resnames) > 6 else ''})") + return "\n".join(lines) + + +def _count_protein_chains(u: mda.Universe) -> list[int]: + """Return the residue count of each protein chain (by segment/fragment).""" + protein = u.select_atoms("protein") + if protein.n_atoms == 0: + # fall back to our amino-acid dictionary if MDAnalysis 'protein' misses + sel = " or ".join(f"resname {r}" for r in sorted(_AMINO)) + protein = u.select_atoms(sel) + if protein.n_atoms == 0: + return [] + lengths = [] + try: + for seg in protein.segments: + n = seg.atoms.select_atoms("name CA").n_residues + if n: + lengths.append(int(n)) + except Exception: + pass + if not lengths: + # single chain fallback + lengths = [int(protein.select_atoms("name CA").n_residues)] + return [n for n in lengths if n > 0] + + +def detect_system(topology: str | Path, coordinates: str | Path | None = None, + peptide_cutoff: int = 30) -> SystemInfo: + """ + Detect the composition and type of an MD system from its topology. + + Parameters + ---------- + topology : path to a .tpr/.gro/.pdb/.psf that carries residue names. + coordinates : optional coordinate file (only needed if the topology lacks it). + peptide_cutoff : a protein chain shorter than this many residues attached to a + larger chain marks a protein-peptide system. + + Returns + ------- + SystemInfo + """ + topology = str(topology) + u = mda.Universe(topology, str(coordinates)) if coordinates else mda.Universe(topology) + + # -- classify residues ------------------------------------------------ # + buckets: dict[ComponentType, dict] = {} + for res in u.residues: + ct = _classify_resname(res.resname) + b = buckets.setdefault(ct, {"resnames": set(), "n_res": 0}) + b["resnames"].add(res.resname.strip().upper()) + b["n_res"] += 1 + + # Refine OTHER into ligand vs buffer-like; keep both under LIGAND but note buffers. + components: dict[str, Component] = {} + ligand_resnames: list[str] = [] + notes: list[str] = [] + for ct, b in buckets.items(): + resnames = sorted(b["resnames"]) + real_ct = ct + if ct == ComponentType.OTHER: + real_ct = ComponentType.LIGAND + ligs = [r for r in resnames if r not in _BUFFER_LIKE] + buf = [r for r in resnames if r in _BUFFER_LIKE] + ligand_resnames.extend(ligs) + if buf: + notes.append(f"Ignoring buffer/crystallisation additives: {', '.join(buf)}") + sel = _selection_for(real_ct, resnames) + # merge into existing (e.g. OTHER folded into LIGAND) + if real_ct.value in components: + existing = components[real_ct.value] + existing.resnames = sorted(set(existing.resnames) | set(resnames)) + existing.n_residues += b["n_res"] + else: + components[real_ct.value] = Component( + ctype=real_ct, resnames=resnames, n_residues=b["n_res"], + n_molecules=b["n_res"], selection=sel) + + chain_lengths = _count_protein_chains(u) + info = SystemInfo( + n_atoms=int(u.atoms.n_atoms), + system_type=SystemType.UNKNOWN, + components=components, + n_protein_chains=len(chain_lengths), + protein_chain_lengths=sorted(chain_lengths, reverse=True), + ligand_resnames=sorted(set(ligand_resnames)), + notes=notes, + ) + info.selections = _build_selections(info) + info.flags = _build_flags(info) + info.system_type = _classify_system(info, peptide_cutoff) + return info + + +def _selection_for(ct: ComponentType, resnames: list[str]) -> str: + if ct == ComponentType.PROTEIN: + return "protein" + if ct in (ComponentType.DNA, ComponentType.RNA): + return "nucleic" + joined = " ".join(resnames) + return f"resname {joined}" if resnames else "not all" + + +def _build_selections(info: SystemInfo) -> dict[str, str]: + sel = {} + if info.has(ComponentType.PROTEIN): + sel["protein"] = "protein" + sel["backbone"] = "protein and backbone" + sel["ca"] = "protein and name CA" + for ct in (ComponentType.DNA, ComponentType.RNA): + if info.has(ct): + sel["nucleic"] = "nucleic" + for ct in (ComponentType.LIGAND, ComponentType.LIPID, ComponentType.ION, + ComponentType.WATER, ComponentType.COFACTOR): + if info.has(ct): + sel[ct.value] = info.components[ct.value].selection + return sel + + +def _build_flags(info: SystemInfo) -> dict[str, bool]: + return { + "has_protein": info.has(ComponentType.PROTEIN), + "has_dna": info.has(ComponentType.DNA), + "has_rna": info.has(ComponentType.RNA), + "has_nucleic": info.has(ComponentType.DNA) or info.has(ComponentType.RNA), + "has_lipid": info.has(ComponentType.LIPID), + "has_ligand": info.has(ComponentType.LIGAND), + "has_cofactor": info.has(ComponentType.COFACTOR), + "has_ions": info.has(ComponentType.ION), + "has_water": info.has(ComponentType.WATER), + "multi_chain": info.n_protein_chains >= 2, + } + + +def _classify_system(info: SystemInfo, peptide_cutoff: int) -> SystemType: + f = info.flags + if not f["has_protein"]: + if f["has_nucleic"]: + return SystemType.NUCLEIC_ONLY + if f["has_lipid"]: + return SystemType.MEMBRANE_ONLY + return SystemType.UNKNOWN + # protein present + if f["has_lipid"]: + return SystemType.PROTEIN_MEMBRANE + if f["has_dna"]: + return SystemType.PROTEIN_DNA + if f["has_rna"]: + return SystemType.PROTEIN_RNA + if f["has_ligand"]: + return SystemType.PROTEIN_LIGAND + if info.n_protein_chains >= 2: + lengths = sorted(info.protein_chain_lengths, reverse=True) + if len(lengths) >= 2 and lengths[-1] < peptide_cutoff <= lengths[0]: + return SystemType.PROTEIN_PEPTIDE + return SystemType.PROTEIN_PROTEIN + return SystemType.PROTEIN_ONLY diff --git a/mdforge/io/__init__.py b/mdforge/io/__init__.py new file mode 100644 index 0000000..2cf932f --- /dev/null +++ b/mdforge/io/__init__.py @@ -0,0 +1,4 @@ +"""I/O: recursive simulation-file discovery, validation, and GROMACS readers.""" + +from mdforge.io.discovery import FileSet, discover_files, pick_trajectory # noqa: F401 +from mdforge.io.validation import validate_fileset # noqa: F401 diff --git a/mdforge/io/discovery.py b/mdforge/io/discovery.py new file mode 100644 index 0000000..64bf16e --- /dev/null +++ b/mdforge/io/discovery.py @@ -0,0 +1,120 @@ +""" +Recursive discovery of GROMACS (and related) simulation files. + +The user only provides a directory; this module finds trajectories, topologies, +structures, energy files, index and topology-parameter files, resolves the best +candidate for each semantic role, and (when several trajectories exist) supports +an interactive choice. +""" + +from __future__ import annotations + +from dataclasses import dataclass, field +from pathlib import Path + +# extension -> semantic category +_TRAJ_EXT = {".xtc", ".trr", ".dcd", ".nc", ".netcdf", ".h5", ".xyz"} +_TOPOLOGY_EXT = {".tpr", ".psf", ".prmtop", ".parm7", ".top"} # carry bonds/masses +_STRUCTURE_EXT = {".gro", ".pdb", ".g96"} +_ENERGY_EXT = {".edr"} +_INDEX_EXT = {".ndx"} +_GMXTOP_EXT = {".top"} +_LOG_EXT = {".log"} + + +@dataclass +class FileSet: + """Resolved semantic file roles for a simulation.""" + directory: Path + trajectory: Path | None = None + topology: Path | None = None # best for MDAnalysis (bonds/masses) + structure: Path | None = None # .gro/.pdb single frame + energy: Path | None = None # .edr + index: Path | None = None # .ndx + gmx_top: Path | None = None # .top (#includes) + log: Path | None = None # md .log + all_trajectories: list[Path] = field(default_factory=list) + all_files: dict[str, list[Path]] = field(default_factory=dict) + + def available_keys(self) -> set[str]: + """Semantic keys that are present (used to gate analyses).""" + keys = set() + for role in ("trajectory", "topology", "structure", "energy", "index", "gmx_top"): + if getattr(self, role) is not None: + keys.add(role) + # 'topology' also satisfies analyses that only need coordinates. + if self.topology or self.structure: + keys.add("topology") + return keys + + def to_dict(self) -> dict: + return {role: (str(getattr(self, role)) if getattr(self, role) else None) + for role in ("directory", "trajectory", "topology", "structure", + "energy", "index", "gmx_top", "log")} + + +def _score_topology(p: Path) -> int: + """Prefer .tpr (full parameters) > .psf/.prmtop > .pdb/.gro.""" + ext = p.suffix.lower() + return {".tpr": 5, ".psf": 4, ".prmtop": 4, ".parm7": 4, + ".pdb": 2, ".gro": 2}.get(ext, 1) + + +def _largest(paths: list[Path]) -> Path | None: + return max(paths, key=lambda p: p.stat().st_size) if paths else None + + +def discover_files(directory: str | Path) -> FileSet: + """Recursively scan ``directory`` and resolve the best file for each role.""" + directory = Path(directory) + if not directory.exists(): + raise FileNotFoundError(f"Input directory does not exist: {directory}") + + by_ext: dict[str, list[Path]] = {} + for p in sorted(directory.rglob("*")): + if p.is_file(): + by_ext.setdefault(p.suffix.lower(), []).append(p) + + traj = sorted([p for ext in _TRAJ_EXT for p in by_ext.get(ext, [])]) + topo_candidates = [p for ext in _TOPOLOGY_EXT | _STRUCTURE_EXT + for p in by_ext.get(ext, [])] + # best topology: highest score, then largest (bigger .tpr = production run) + topology = None + if topo_candidates: + topology = sorted(topo_candidates, + key=lambda p: (_score_topology(p), p.stat().st_size), + reverse=True)[0] + structure = _largest([p for ext in _STRUCTURE_EXT for p in by_ext.get(ext, [])]) + + fs = FileSet( + directory=directory, + trajectory=_largest(traj), # default: the largest (usually production) + topology=topology, + structure=structure, + energy=_largest([p for ext in _ENERGY_EXT for p in by_ext.get(ext, [])]), + index=_largest([p for ext in _INDEX_EXT for p in by_ext.get(ext, [])]), + gmx_top=_largest([p for ext in _GMXTOP_EXT for p in by_ext.get(ext, [])]), + log=_largest([p for ext in _LOG_EXT for p in by_ext.get(ext, [])]), + all_trajectories=traj, + all_files={k: v for k, v in by_ext.items()}, + ) + return fs + + +def pick_trajectory(fs: FileSet, interactive: bool = False, + ask=None) -> Path | None: + """ + Choose the trajectory to analyse. + + If several trajectories are found and ``interactive`` is set, ``ask`` (a + callable returning the selected index) is used; otherwise the largest is + kept and a note is recorded by the caller. + """ + if len(fs.all_trajectories) <= 1: + return fs.trajectory + if interactive and ask is not None: + options = [f"{p.name} ({p.stat().st_size/1e6:.1f} MB)" + for p in fs.all_trajectories] + idx = ask("Multiple trajectories found — which one to analyse?", options) + return fs.all_trajectories[idx] + return fs.trajectory # default: largest diff --git a/mdforge/io/validation.py b/mdforge/io/validation.py new file mode 100644 index 0000000..84c59e4 --- /dev/null +++ b/mdforge/io/validation.py @@ -0,0 +1,67 @@ +""" +Validation of a discovered :class:`FileSet`. + +Produces clear, actionable messages: what is missing, why it is needed, and how +to obtain it -- so a researcher is never left guessing. +""" + +from __future__ import annotations + +from dataclasses import dataclass + +from mdforge.io.discovery import FileSet + +# role -> (why it is needed, how to obtain it) +_GUIDANCE = { + "trajectory": ("the time-series coordinates every dynamic analysis needs", + "produce one with `gmx mdrun` (writes .xtc/.trr) or `gmx trjconv`"), + "topology": ("atom names, bonds, masses and charges used for selections/analyses", + "use the run input `.tpr` (best) or a `.gro`/`.pdb` of the system"), + "energy": ("thermodynamic observables (T, P, density, energies)", + "the `.edr` written by `gmx mdrun`; required only for energy analysis"), + "index": ("custom atom groups (e.g. for MM/PBSA group definitions)", + "create with `gmx make_ndx`; optional for most analyses"), + "gmx_top": ("the processed topology with #include's", + "the `topol.top` used for the run; required only for MM/PBSA"), +} + + +@dataclass +class ValidationResult: + ok: bool + errors: list[str] + warnings: list[str] + + def report(self) -> str: + lines = [] + for e in self.errors: + lines.append(f" [ERROR] {e}") + for w in self.warnings: + lines.append(f" [warn ] {w}") + return "\n".join(lines) if lines else " all required inputs present" + + +def validate_fileset(fs: FileSet) -> ValidationResult: + """Check that the minimum inputs (trajectory + topology) are present.""" + errors, warnings = [], [] + + def explain(role: str) -> str: + why, how = _GUIDANCE.get(role, ("", "")) + return f"missing '{role}' — needed for {why}. To obtain it: {how}." + + if fs.trajectory is None: + errors.append(explain("trajectory")) + if fs.topology is None and fs.structure is None: + errors.append(explain("topology")) + + # Non-fatal: features that will simply be skipped. + if fs.energy is None: + warnings.append("no '.edr' found — energy/stability analysis will be skipped.") + if fs.gmx_top is None: + warnings.append("no '.top' found — MM/PBSA workflow generation will be limited.") + if len(fs.all_trajectories) > 1: + warnings.append(f"{len(fs.all_trajectories)} trajectories found — " + f"defaulting to the largest ({fs.trajectory.name}); " + f"use --interactive to choose.") + + return ValidationResult(ok=not errors, errors=errors, warnings=warnings) diff --git a/mdforge/plotting/__init__.py b/mdforge/plotting/__init__.py new file mode 100644 index 0000000..038d90a --- /dev/null +++ b/mdforge/plotting/__init__.py @@ -0,0 +1,6 @@ +"""Publication-quality plotting: consistent Nature-like style + PNG/PDF saving.""" + +from mdforge.plotting.style import ( # noqa: F401 + set_style, new_axes, add_reference_line, PALETTE, SEQ_CMAP, DIV_CMAP, FEL_CMAP, +) +from mdforge.plotting.figures import save_figure # noqa: F401 diff --git a/mdforge/plotting/figures.py b/mdforge/plotting/figures.py new file mode 100644 index 0000000..08aa4c0 --- /dev/null +++ b/mdforge/plotting/figures.py @@ -0,0 +1,22 @@ +"""Figure saving: always PNG + vector PDF at 300 dpi.""" + +from __future__ import annotations + +from pathlib import Path + +import matplotlib.pyplot as plt + + +def save_figure(fig, path_no_ext: Path | str, dpi: int = 300, + formats: tuple[str, ...] = ("png", "pdf")) -> list[Path]: + """Save ``fig`` to every format at high resolution; returns written paths.""" + path_no_ext = Path(path_no_ext) + path_no_ext.parent.mkdir(parents=True, exist_ok=True) + fig.tight_layout() + written = [] + for ext in formats: + out = path_no_ext.with_suffix(f".{ext}") + fig.savefig(out, dpi=dpi, bbox_inches="tight") + written.append(out) + plt.close(fig) + return written diff --git a/mdforge/plotting/style.py b/mdforge/plotting/style.py new file mode 100644 index 0000000..9499388 --- /dev/null +++ b/mdforge/plotting/style.py @@ -0,0 +1,45 @@ +"""Global matplotlib style + palette (Nature-like: white bg, large fonts, no grid).""" + +from __future__ import annotations + +import matplotlib +matplotlib.use("Agg") +import matplotlib.pyplot as plt + +PALETTE = { + "primary": "#1f6f8b", "secondary": "#e07a5f", "accent": "#3d405b", + "green": "#5f9e6e", "purple": "#8367c7", "muted": "#9aa0a6", + "helix": "#c0392b", "sheet": "#2471a3", "coil": "#7f8c8d", +} +SEQ_CMAP = "viridis" +DIV_CMAP = "RdBu_r" +FEL_CMAP = "nipy_spectral" + + +def set_style() -> None: + """Apply the global publication style. Call once per script/session.""" + plt.rcParams.update({ + "figure.figsize": (7.2, 4.8), "figure.dpi": 110, "savefig.dpi": 300, + "figure.facecolor": "white", "axes.facecolor": "white", + "savefig.facecolor": "white", + "font.size": 13, "font.family": "sans-serif", + "font.sans-serif": ["Arial", "DejaVu Sans", "Helvetica", "sans-serif"], + "axes.titlesize": 15, "axes.labelsize": 14, "axes.titleweight": "bold", + "xtick.labelsize": 12, "ytick.labelsize": 12, + "legend.fontsize": 11, "legend.frameon": False, + "axes.spines.top": False, "axes.spines.right": False, + "axes.linewidth": 1.1, "axes.grid": False, + "xtick.direction": "out", "ytick.direction": "out", + "xtick.major.size": 5, "ytick.major.size": 5, + "lines.linewidth": 1.8, "lines.solid_capstyle": "round", + "legend.handlelength": 1.6, + }) + + +def new_axes(figsize: tuple[float, float] = (7.2, 4.8)): + return plt.subplots(figsize=figsize) + + +def add_reference_line(ax, y, label=None, color=None): + ax.axhline(y, ls="--", lw=1.2, color=color or PALETTE["accent"], + alpha=0.8, label=label, zorder=1) diff --git a/mdforge/report/__init__.py b/mdforge/report/__init__.py new file mode 100644 index 0000000..d755470 --- /dev/null +++ b/mdforge/report/__init__.py @@ -0,0 +1,3 @@ +"""Automatic report generation (Markdown + HTML + best-effort PDF).""" + +from mdforge.report.generator import generate_report # noqa: F401 diff --git a/mdforge/report/generator.py b/mdforge/report/generator.py new file mode 100644 index 0000000..5709502 --- /dev/null +++ b/mdforge/report/generator.py @@ -0,0 +1,202 @@ +""" +Structured report builder. + +Content is assembled as a list of typed blocks (heading/paragraph/table/image), +then rendered to Markdown and self-contained HTML from the same source, so the +two formats never drift. PDF is attempted via weasyprint if installed. +""" + +from __future__ import annotations + +import base64 +from datetime import date +from pathlib import Path + + +# --------------------------------------------------------------------------- # +# Document model + renderers +# --------------------------------------------------------------------------- # +def _md(blocks) -> str: + out = [] + for kind, payload in blocks: + if kind == "h1": + out.append(f"# {payload}\n") + elif kind == "h2": + out.append(f"## {payload}\n") + elif kind == "p": + out.append(f"{payload}\n") + elif kind == "table": + headers, rows = payload + out.append("| " + " | ".join(headers) + " |") + out.append("|" + "|".join(["---"] * len(headers)) + "|") + for r in rows: + out.append("| " + " | ".join(str(c) for c in r) + " |") + out.append("") + elif kind == "img": + path, caption = payload + rel = f"../figures/{Path(path).name}" + out.append(f"\n![{caption}]({rel})\n\n*Figure. {caption}.*\n") + 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: + return f"{float(x):.{nd}f}" + except (TypeError, ValueError): + return "n/a" + + +def generate_report(ctx, results: dict, manifest, formats=("md", "html")) -> list[Path]: + system = ctx.system + blocks = [] + blocks.append(("h1", f"MD Analysis Report — {system.system_type.value}")) + blocks.append(("p", f"*Generated {date.today().isoformat()} by mdforge " + f"v{manifest.data.get('mdforge_version')}.*")) + + # -- system ----------------------------------------------------------- # + blocks.append(("h2", "1. System")) + comp_rows = [[c.ctype.value, c.n_residues, ", ".join(c.resnames[:8])] + for c in system.components.values() if c.n_residues] + blocks.append(("table", (["Component", "Residues", "Resnames"], comp_rows))) + blocks.append(("p", f"Detected system type: **{system.system_type.value}** · " + f"{system.n_protein_chains} protein chain(s) · " + f"{system.n_atoms:,} atoms.")) + + # -- methods ---------------------------------------------------------- # + blocks.append(("h2", "2. Methods & reproducibility")) + libs = manifest.data.get("libraries", {}) + lib_str = ", ".join(f"{k} {v}" for k, v in libs.items() if v) + git = manifest.data.get("git", {}) + blocks.append(("p", f"Analyses were run with mdforge (Python " + f"{manifest.data.get('python')}). Key libraries: {lib_str}. " + f"Git commit: `{git.get('commit')}`" + f"{' (dirty)' if git.get('dirty') else ''}. " + f"Full provenance in `manifest.json`.")) + + # -- results ---------------------------------------------------------- # + blocks.append(("h2", "3. Results")) + key_rows = _key_result_rows(results) + if key_rows: + blocks.append(("table", (["Observable", "Value"], key_rows))) + + # -- figures ---------------------------------------------------------- # + blocks.append(("h2", "4. Figures")) + for name, summ in results.items(): + if isinstance(summ, dict) and summ.get("figure"): + fig = ctx.config.figures_dir / f"{summ['figure']}.png" + if fig.exists(): + blocks.append(("img", (fig, name))) + + # -- interpretation --------------------------------------------------- # + blocks.append(("h2", "5. Interpretation")) + for line in _interpretation(results): + blocks.append(("p", line)) + + # -- limitations ------------------------------------------------------ # + blocks.append(("h2", "6. Limitations & provenance")) + blocks.append(("p", "Convergence diagnostics are necessary but not sufficient " + "evidence of equilibration; consider replicate/longer runs. " + "All parameters, seeds, library versions, input fingerprints " + "and per-analysis runtimes are recorded in `manifest.json` / " + "`manifest.yaml` for exact reproduction.")) + + # -- render ----------------------------------------------------------- # + ctx.config.report_dir.mkdir(parents=True, exist_ok=True) + written = [] + if "md" in formats: + p = ctx.config.report_dir / "report.md" + p.write_text(_md(blocks), encoding="utf-8") + written.append(p) + if "html" in formats: + p = ctx.config.report_dir / "report.html" + p.write_text(_html(blocks), 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)) + written.append(pdf) + except Exception: + pass # weasyprint not installed; MD/HTML still produced + return written + + +def _key_result_rows(results: dict) -> list[list]: + rows = [] + g = lambda d, *k: _nested(d, k) + if "rmsd" in results: + rows.append(["Backbone RMSD (mean)", f"{_fmt(g(results,'rmsd','backbone','mean'))} nm"]) + rows.append(["RMSD converged?", g(results, "rmsd", "convergence", "converged")]) + if "rog" in results: + rows.append(["Radius of gyration (mean)", f"{_fmt(g(results,'rog','rg','mean'))} nm"]) + if "rmsf" in results: + 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²"]) + return [r for r in rows if r[1] not in (None, "n/a")] + + +def _nested(d, keys): + for k in keys: + if not isinstance(d, dict) or k not in d: + return None + d = d[k] + return d + + +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.") + rg = _nested(results, ("rog",)) + if rg: + init, fin = rg.get("rg_initial_nm"), rg.get("rg_final_nm") + if init and fin: + trend = "expanded" if fin > init + 0.15 else "compacted" if fin < init - 0.15 else "kept a stable size" + lines.append(f"**Compactness.** The solute {trend} " + f"(Rg {_fmt(init)} → {_fmt(fin)} nm).") + if not lines: + lines.append("See the per-analysis CSVs and figures for detailed results.") + return lines diff --git a/mdforge/statistics/__init__.py b/mdforge/statistics/__init__.py new file mode 100644 index 0000000..53c6962 --- /dev/null +++ b/mdforge/statistics/__init__.py @@ -0,0 +1,8 @@ +"""Statistics: descriptive stats, time-series helpers, correlation, bootstrap CIs.""" + +from mdforge.statistics.descriptive import describe, summary_frame, bootstrap_ci # noqa: F401 +from mdforge.statistics.timeseries import ( # noqa: F401 + moving_average, rolling_std, running_mean, block_average, block_average_sem, + plateau_detection, +) +from mdforge.statistics.correlation import correlation_matrices # noqa: F401 diff --git a/mdforge/statistics/correlation.py b/mdforge/statistics/correlation.py new file mode 100644 index 0000000..2e4a3b4 --- /dev/null +++ b/mdforge/statistics/correlation.py @@ -0,0 +1,10 @@ +"""Correlation matrices (Pearson + Spearman).""" + +from __future__ import annotations + +import pandas as pd + + +def correlation_matrices(df: pd.DataFrame) -> tuple[pd.DataFrame, pd.DataFrame]: + """Return (Pearson, Spearman) correlation matrices for numeric columns.""" + return df.corr(method="pearson"), df.corr(method="spearman") diff --git a/mdforge/statistics/descriptive.py b/mdforge/statistics/descriptive.py new file mode 100644 index 0000000..aa86ab0 --- /dev/null +++ b/mdforge/statistics/descriptive.py @@ -0,0 +1,64 @@ +"""Descriptive statistics and bootstrap confidence intervals.""" + +from __future__ import annotations + +from typing import Iterable + +import numpy as np +import pandas as pd +from scipy import stats as sps + + +def describe(x: Iterable[float], name: str = "value") -> dict: + """Mean/median/std/SEM/95% t-CI/min/max for a 1-D series (NaNs dropped).""" + x = np.asarray(list(x), dtype=float) + x = x[np.isfinite(x)] + n = x.size + mean = float(np.mean(x)) if n else np.nan + std = float(np.std(x, ddof=1)) if n > 1 else 0.0 + sem = std / np.sqrt(n) if n > 1 else 0.0 + if n > 1: + t = sps.t.ppf(0.975, df=n - 1) + lo, hi = mean - t * sem, mean + t * sem + else: + lo = hi = mean + return { + "metric": name, "n": n, "mean": mean, + "median": float(np.median(x)) if n else np.nan, + "std": std, "sem": sem, + "ci95_low": float(lo), "ci95_high": float(hi), + "min": float(np.min(x)) if n else np.nan, + "max": float(np.max(x)) if n else np.nan, + } + + +def bootstrap_ci(x: Iterable[float], statistic=np.mean, n_boot: int = 10000, + ci: float = 95.0, seed: int | None = 0) -> dict: + """ + Bootstrap confidence interval for a statistic of a (correlated) series. + + Uses the percentile method. For MD time-series the samples are correlated, so + this complements the t-based CI in :func:`describe`. + """ + rng = np.random.default_rng(seed) + x = np.asarray(list(x), dtype=float) + x = x[np.isfinite(x)] + n = x.size + if n < 2: + v = float(statistic(x)) if n else np.nan + return {"estimate": v, "ci_low": v, "ci_high": v, "n_boot": 0} + idx = rng.integers(0, n, size=(n_boot, n)) + boot = np.array([statistic(x[i]) for i in idx]) + a = (100 - ci) / 2 + return { + "estimate": float(statistic(x)), + "ci_low": float(np.percentile(boot, a)), + "ci_high": float(np.percentile(boot, 100 - a)), + "n_boot": n_boot, + } + + +def summary_frame(series_map: dict[str, Iterable[float]]) -> pd.DataFrame: + """Tidy summary-statistics DataFrame from {name: series}.""" + return pd.DataFrame([describe(v, name=k) for k, v in series_map.items()] + ).set_index("metric") diff --git a/mdforge/statistics/timeseries.py b/mdforge/statistics/timeseries.py new file mode 100644 index 0000000..9078d5f --- /dev/null +++ b/mdforge/statistics/timeseries.py @@ -0,0 +1,65 @@ +"""Time-series helpers: moving/rolling/running averages, block averaging, plateau.""" + +from __future__ import annotations + +from typing import Iterable + +import numpy as np +import pandas as pd + + +def moving_average(x: Iterable[float], window: int = 20) -> np.ndarray: + s = pd.Series(np.asarray(x, dtype=float)) + return s.rolling(window=window, center=True, min_periods=1).mean().to_numpy() + + +def rolling_std(x: Iterable[float], window: int = 20) -> np.ndarray: + s = pd.Series(np.asarray(x, dtype=float)) + return s.rolling(window=window, center=True, min_periods=1).std().to_numpy() + + +def running_mean(x: Iterable[float]) -> np.ndarray: + x = np.asarray(x, dtype=float) + return np.cumsum(x) / np.arange(1, x.size + 1) + + +def block_average(x: Iterable[float], n_blocks: int = 10) -> pd.DataFrame: + """Contiguous-block means (correlation-aware uncertainty estimate).""" + x = np.asarray(x, dtype=float) + x = x[np.isfinite(x)] + n = x.size + n_blocks = max(1, min(n_blocks, n)) + bs = n // n_blocks + means = [x[b * bs:(b + 1) * bs].mean() for b in range(n_blocks) if x[b*bs:(b+1)*bs].size] + return pd.DataFrame({"block": np.arange(1, len(means) + 1), "block_mean": means}) + + +def block_average_sem(x: Iterable[float], n_blocks: int = 10) -> float: + m = block_average(x, n_blocks)["block_mean"].to_numpy() + return float(np.std(m, ddof=1) / np.sqrt(m.size)) if m.size >= 2 else 0.0 + + +def plateau_detection(t, y, tail_fraction: float = 0.5, slope_tol: float | None = None) -> dict: + """Assess convergence via tail-slope and first/second-half drift.""" + t = np.asarray(t, dtype=float) + y = np.asarray(y, dtype=float) + n = y.size + i0 = int(n * (1.0 - tail_fraction)) + tt, yy = t[i0:], y[i0:] + slope, _ = np.polyfit(tt, yy, 1) + span = (tt[-1] - tt[0]) or 1.0 + if slope_tol is None: + slope_tol = np.std(yy, ddof=1) / span if span else np.inf + half = n // 2 + mean_first, mean_second = np.mean(y[:half]), np.mean(y[half:]) + return { + "tail_slope": float(slope), + "slope_tolerance": float(slope_tol), + "mean_first_half": float(mean_first), + "mean_second_half": float(mean_second), + "drift_second_minus_first": float(mean_second - mean_first), + "tail_mean": float(np.mean(yy)), + "tail_std": float(np.std(yy, ddof=1)), + # floor guards the degenerate perfectly-flat-tail case (slope_tol == 0) + "converged": bool(abs(slope) <= max(slope_tol, 1e-9)), + } diff --git a/pyproject.toml b/pyproject.toml new file mode 100644 index 0000000..149d196 --- /dev/null +++ b/pyproject.toml @@ -0,0 +1,60 @@ +[build-system] +requires = ["setuptools>=64", "wheel"] +build-backend = "setuptools.build_meta" + +[project] +name = "molecular-dynamics-forge" +version = "0.2.0" +description = "Reusable, extensible framework for reproducible analysis of GROMACS molecular dynamics simulations of arbitrary biomolecular systems." +readme = "README.md" +requires-python = ">=3.10" +license = { text = "MIT" } +authors = [{ name = "Hossam Mahmoud" }] +keywords = ["molecular-dynamics", "gromacs", "mdanalysis", "bioinformatics", + "trajectory-analysis", "reproducibility"] +classifiers = [ + "Programming Language :: Python :: 3", + "License :: OSI Approved :: MIT License", + "Topic :: Scientific/Engineering :: Bio-Informatics", + "Operating System :: OS Independent", +] +dependencies = [ + "numpy>=1.24", + "scipy>=1.10", + "pandas>=2.0", + "matplotlib>=3.7", + "MDAnalysis>=2.7", + "mdtraj>=1.10", + "scikit-learn>=1.3", + "networkx>=3.0", + "tqdm>=4.65", + "pyyaml>=6.0", +] + +[project.optional-dependencies] +energy = ["panedr>=0.8"] # read .edr thermodynamics +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"] + +[project.urls] +Homepage = "https://github.com/SamDozer/molecular-dynamics-forge" +Repository = "https://github.com/SamDozer/molecular-dynamics-forge" + +[project.scripts] +mdforge = "mdforge.cli.main:main" + +# Third-party packages advertise plugins under this entry-point group. +# [project.entry-points."mdforge.plugins"] +# my_analysis = "my_package.my_module" + +[tool.setuptools.packages.find] +include = ["mdforge*"] + +[tool.pytest.ini_options] +testpaths = ["tests"] + +[tool.ruff] +line-length = 100 diff --git a/requirements.txt b/requirements.txt new file mode 100644 index 0000000..1742009 --- /dev/null +++ b/requirements.txt @@ -0,0 +1,24 @@ +# Alpha-zein (A8HNE1) MD analysis pipeline -- Python dependencies +# Tested with Python 3.11. Versions listed are those the pipeline was validated +# against; the lower bounds are generally safe. + +numpy>=1.24 +scipy>=1.10 +pandas>=2.0 +matplotlib>=3.7 +seaborn>=0.12 +tqdm>=4.65 +pyyaml>=6.0 # mdforge config-driven workflows + +# MD trajectory analysis +MDAnalysis>=2.7 # validated with 2.10.0 +mdtraj>=1.10 # validated with 1.11.1 ; provides built-in DSSP +panedr>=0.8 # read GROMACS .edr directly (no `gmx energy` needed) + +# Machine learning / clustering / networks +scikit-learn>=1.3 # validated with 1.9.0 +networkx>=3.0 # residue interaction network + +# Interaction fingerprints +prolif>=2.0 # validated with 2.2.0 +rdkit>=2023.3 # ProLIF backend (validated with 2026.3.3) diff --git a/tests/test_core.py b/tests/test_core.py new file mode 100644 index 0000000..b10f7a3 --- /dev/null +++ b/tests/test_core.py @@ -0,0 +1,101 @@ +"""Unit tests for the mdforge core (no large trajectory data required).""" + +from __future__ import annotations + +import numpy as np + +from mdforge.core.system import ( + _classify_resname, _build_flags, _classify_system, + ComponentType, SystemType, SystemInfo, Component, +) +from mdforge.core.registry import registry +from mdforge import statistics as st + + +# --------------------------------------------------------------------------- # +# System detection logic +# --------------------------------------------------------------------------- # +def test_classify_resname(): + assert _classify_resname("ALA") is ComponentType.PROTEIN + assert _classify_resname("HSD") is ComponentType.PROTEIN # CHARMM variant + assert _classify_resname("DA") is ComponentType.DNA + assert _classify_resname("U") is ComponentType.RNA + assert _classify_resname("TIP3") is ComponentType.WATER + assert _classify_resname("POT") is ComponentType.ION + assert _classify_resname("POPC") is ComponentType.LIPID + assert _classify_resname("HEM") is ComponentType.COFACTOR + assert _classify_resname("LIG") is ComponentType.OTHER # -> ligand + + +def _make_info(**components) -> SystemInfo: + comps = {} + for ctype, n in components.items(): + ct = ComponentType(ctype) + comps[ct.value] = Component(ct, [ctype.upper()], n, n, "sel") + info = SystemInfo(n_atoms=1000, system_type=SystemType.UNKNOWN, components=comps) + return info + + +def test_system_type_classification(): + def classify(chains=1, **comp): + info = _make_info(**comp) + info.n_protein_chains = chains + info.protein_chain_lengths = [187] * chains + info.flags = _build_flags(info) + return _classify_system(info, peptide_cutoff=30) + + assert classify(chains=1, protein=187) is SystemType.PROTEIN_ONLY + assert classify(chains=2, protein=374) is SystemType.PROTEIN_PROTEIN + assert classify(chains=1, protein=187, ligand=1) is SystemType.PROTEIN_LIGAND + assert classify(chains=1, protein=187, dna=20) is SystemType.PROTEIN_DNA + assert classify(chains=1, protein=187, rna=20) is SystemType.PROTEIN_RNA + assert classify(chains=1, protein=187, lipid=200) is SystemType.PROTEIN_MEMBRANE + + +# --------------------------------------------------------------------------- # +# Registry / plugin discovery +# --------------------------------------------------------------------------- # +def test_registry_loads_builtins_and_plugin(): + registry.ensure_builtins_loaded() + names = registry.names() + for expected in ("rmsd", "rmsf", "rog", "sasa"): + assert expected in names + # the bundled example plugin must be auto-discovered + assert "end_to_end" in names + + +def test_registry_selection_respects_system(): + registry.ensure_builtins_loaded() + info = _make_info(protein=187, water=1000) + info.n_protein_chains = 1 + info.flags = _build_flags(info) + info.system_type = SystemType.PROTEIN_ONLY + selected, _ = registry.select(info, available_files={"trajectory", "topology"}) + assert any(c.name == "rmsd" for c in selected) + + +# --------------------------------------------------------------------------- # +# Statistics +# --------------------------------------------------------------------------- # +def test_describe_and_ci(): + x = np.arange(100.0) + d = st.describe(x, "x") + assert d["n"] == 100 + assert abs(d["mean"] - 49.5) < 1e-9 + assert d["ci95_low"] < d["mean"] < d["ci95_high"] + + +def test_moving_average_and_plateau(): + x = np.concatenate([np.linspace(0, 1, 50), np.ones(50)]) + ma = st.moving_average(x, 10) + assert len(ma) == len(x) + conv = st.plateau_detection(np.arange(len(x)), x) + assert conv["converged"] is True # flat tail + + +def test_bootstrap_ci(): + rng = np.random.default_rng(0) + x = rng.normal(5.0, 1.0, 500) + b = st.bootstrap_ci(x, seed=0) + assert b["ci_low"] < b["estimate"] < b["ci_high"] + assert abs(b["estimate"] - 5.0) < 0.3 From 4d99d7f929d2a667eb3d0a8bef779ae574efaa2c Mon Sep 17 00:00:00 2001 From: Hossam Mahmoud Date: Tue, 7 Jul 2026 17:47:19 +0300 Subject: [PATCH 2/4] Port 14 more analyses to mdforge (18 total, all validated end-to-end) Adds as BaseAnalysis modules: com, structural (Dmax/anisotropy/volume), native_contacts (Q), contact_map, dccm, secondary_structure (DSSP), salt_bridges, energies (from .edr), pca (+FEL), clustering (silhouette-optimal k), rin, convergence (RMSIP/block-avg/plateau), hbonds (streamed, shell-restricted), statistics_summary (cross-metric correlations). - Adds `order` to BaseAnalysis so analyses that consume other analyses' CSVs (rin=150, convergence=200, statistics_summary=300) run after their producers; registry sorts the selection accordingly. - Adds PROTEIN_SYSTEMS / COMPLEX_SYSTEMS convenience sets in core.system. Validated: full pipeline runs 18/18 analyses with zero errors on the reference 608k-atom protein-only system (25 figures + manifest + report). Remaining to port: prolif, mmpbsa, and the system-specific analyses (interface RMSD/BSA, ligand RMSD/contacts/pocket, DNA groove). Co-Authored-By: Claude Opus 4.8 --- mdforge/analysis/clustering.py | 82 +++++++++++++++++++++ mdforge/analysis/com.py | 70 ++++++++++++++++++ mdforge/analysis/contact_map.py | 67 +++++++++++++++++ mdforge/analysis/convergence.py | 95 +++++++++++++++++++++++++ mdforge/analysis/dccm.py | 49 +++++++++++++ mdforge/analysis/energies.py | 78 ++++++++++++++++++++ mdforge/analysis/hbonds.py | 77 ++++++++++++++++++++ mdforge/analysis/native_contacts.py | 61 ++++++++++++++++ mdforge/analysis/pca.py | 78 ++++++++++++++++++++ mdforge/analysis/rin.py | 75 +++++++++++++++++++ mdforge/analysis/salt_bridges.py | 71 ++++++++++++++++++ mdforge/analysis/secondary_structure.py | 66 +++++++++++++++++ mdforge/analysis/statistics_summary.py | 81 +++++++++++++++++++++ mdforge/analysis/structural.py | 80 +++++++++++++++++++++ mdforge/core/base.py | 3 + mdforge/core/registry.py | 4 +- mdforge/core/system.py | 14 ++++ 17 files changed, 1049 insertions(+), 2 deletions(-) create mode 100644 mdforge/analysis/clustering.py create mode 100644 mdforge/analysis/com.py create mode 100644 mdforge/analysis/contact_map.py create mode 100644 mdforge/analysis/convergence.py create mode 100644 mdforge/analysis/dccm.py create mode 100644 mdforge/analysis/energies.py create mode 100644 mdforge/analysis/hbonds.py create mode 100644 mdforge/analysis/native_contacts.py create mode 100644 mdforge/analysis/pca.py create mode 100644 mdforge/analysis/rin.py create mode 100644 mdforge/analysis/salt_bridges.py create mode 100644 mdforge/analysis/secondary_structure.py create mode 100644 mdforge/analysis/statistics_summary.py create mode 100644 mdforge/analysis/structural.py diff --git a/mdforge/analysis/clustering.py b/mdforge/analysis/clustering.py new file mode 100644 index 0000000..6bc9d92 --- /dev/null +++ b/mdforge/analysis/clustering.py @@ -0,0 +1,82 @@ +"""Conformational clustering on pairwise Cα RMSD (optimal k via silhouette).""" + +from __future__ import annotations + +import numpy as np +import pandas as pd +import mdtraj as md +from sklearn.cluster import AgglomerativeClustering +from sklearn.metrics import silhouette_score +from tqdm import tqdm + +from mdforge.core.base import BaseAnalysis +from mdforge.core.system import PROTEIN_SYSTEMS +from mdforge import plotting + +_COLORS = ["#1f6f8b", "#e07a5f", "#5f9e6e", "#8367c7", "#d4a017", "#c0392b", + "#2471a3", "#7f8c8d"] + + +class Clustering(BaseAnalysis): + name = "clustering" + label = "Conformational clustering" + category = "dynamics" + required_files = {"trajectory", "topology"} + supported_systems = PROTEIN_SYSTEMS + outputs = ["results/cluster_populations.csv", "figures/cluster_populations.png"] + default_params = {"k_range": (2, 3, 4, 5, 6, 7, 8)} + + def run(self, ctx) -> dict: + p = self.params(ctx) + plotting.set_style() + ctx.core_universe() + traj = md.load(str(ctx.config.data_dir / "core.xtc"), + top=str(ctx.config.data_dir / "core.pdb")) + traj.superpose(traj, 0) + ca = traj.topology.select("name CA") + nfr = traj.n_frames + D = np.empty((nfr, nfr), dtype=np.float32) + for i in tqdm(range(nfr), desc="[clustering] RMSD matrix", unit="frame"): + D[i] = md.rmsd(traj, traj, i, atom_indices=ca) + D = 0.5 * (D + D.T).astype(np.float64) + np.fill_diagonal(D, 0.0); D = np.clip(D, 0, None) + + scores, labels_by_k = {}, {} + for k in p["k_range"]: + if k >= nfr: + continue + lbl = AgglomerativeClustering(n_clusters=k, metric="precomputed", + linkage="average").fit_predict(D) + scores[k] = float(silhouette_score(D, lbl, metric="precomputed")) + labels_by_k[k] = lbl + best_k = max(scores, key=scores.get); labels = labels_by_k[best_k] + order = pd.Series(labels).value_counts().index.tolist() + remap = {o: n for n, o in enumerate(order)} + labels = np.array([remap[l] for l in labels]) + + struct_dir = ctx.config.results_dir / "cluster_structures" + struct_dir.mkdir(parents=True, exist_ok=True) + pops, medoids = [], [] + for c in range(best_k): + members = np.where(labels == c)[0] + medoid = members[np.argmin(D[np.ix_(members, members)].mean(axis=1))] + medoids.append(int(medoid)); pops.append(len(members)) + traj[medoid].save_pdb(str(struct_dir / f"cluster_{c+1:02d}.pdb")) + pop_df = pd.DataFrame({"cluster": np.arange(1, best_k+1), "population": pops, + "percent": 100*np.array(pops)/nfr, "medoid_frame": medoids}) + ctx.write_csv(pop_df, "cluster_populations.csv") + pd.DataFrame({"k": list(scores), "silhouette": list(scores.values())}).to_csv( + ctx.csv_path("cluster_silhouette.csv"), index=False) + + fig, ax = plotting.new_axes() + colors = [_COLORS[i % len(_COLORS)] for i in range(best_k)] + ax.bar(pop_df["cluster"], pop_df["percent"], color=colors) + for x, pc in zip(pop_df["cluster"], pop_df["percent"]): + ax.text(x, pc + 0.5, f"{pc:.0f}%", ha="center", fontsize=10) + ax.set_xlabel("Cluster"); ax.set_ylabel("Population (%)") + ax.set_title(f"Conformational clusters (k = {best_k})") + plotting.save_figure(fig, ctx.fig_path("cluster_populations"), dpi=ctx.config.dpi) + + return {"optimal_k": int(best_k), "silhouette": scores[best_k], + "dominant_cluster_percent": float(pop_df["percent"].iloc[0]), + "figure": "cluster_populations"} diff --git a/mdforge/analysis/com.py b/mdforge/analysis/com.py new file mode 100644 index 0000000..718d047 --- /dev/null +++ b/mdforge/analysis/com.py @@ -0,0 +1,70 @@ +"""Centre-of-mass tracking (universal): displacement, path length, 3D trajectory.""" + +from __future__ import annotations + +import numpy as np +import pandas as pd +import matplotlib.pyplot as plt +from mpl_toolkits.mplot3d import Axes3D # noqa: F401 (registers 3d projection) + +from mdforge.core.base import BaseAnalysis +from mdforge import plotting +from mdforge import statistics as st +from mdforge.plotting import PALETTE + + +class CenterOfMass(BaseAnalysis): + name = "com" + label = "Centre-of-mass tracking" + category = "dynamics" + required_files = {"trajectory", "topology"} + supported_systems = {"*"} + outputs = ["results/com.csv", "figures/com_xyz.png", "figures/com_3d.png"] + + def run(self, ctx) -> dict: + plotting.set_style() + u = ctx.core_universe() + solute = u.select_atoms("protein" if ctx.system.flags.get("has_protein", True) else "all") + if solute.n_atoms == 0: + solute = u.atoms + n = len(u.trajectory) + coms = np.empty((n, 3)); boxes = np.empty((n, 3)); times = np.empty(n) + for i, ts in enumerate(ctx.iter_frames(u, desc="[com]")): + coms[i] = solute.center_of_mass() + boxes[i] = ts.dimensions[:3] + times[i] = ts.time / 1000.0 + + coms_nm, box_nm = coms / 10.0, boxes / 10.0 + disp = np.diff(coms_nm, axis=0) + disp -= box_nm[1:] * np.round(disp / np.where(box_nm[1:] == 0, 1, box_nm[1:])) + cont = np.vstack([coms_nm[0], coms_nm[0] + np.cumsum(disp, axis=0)]) + rel = cont - cont[0] + cumdist = np.concatenate([[0.0], np.cumsum(np.linalg.norm(disp, axis=1))]) + net = np.linalg.norm(rel, axis=1) + + ctx.write_csv(pd.DataFrame({"time_ns": times, "com_x_nm": rel[:, 0], + "com_y_nm": rel[:, 1], "com_z_nm": rel[:, 2], + "net_displacement_nm": net, + "cumulative_distance_nm": cumdist}), "com.csv") + + fig, ax = plotting.new_axes() + for k, c in zip("xyz", (PALETTE["primary"], PALETTE["secondary"], PALETTE["green"])): + ax.plot(times, rel[:, "xyz".index(k)], color=c, lw=1.6, label=k.upper()) + ax.set_xlabel("Time (ns)"); ax.set_ylabel("COM displacement (nm)") + ax.set_title("Protein centre-of-mass displacement") + ax.legend(ncol=3, loc="upper left"); ax.margins(x=0.01) + plotting.save_figure(fig, ctx.fig_path("com_xyz"), dpi=ctx.config.dpi) + + fig = plt.figure(figsize=(6.4, 5.6)) + ax = fig.add_subplot(111, projection="3d") + sc = ax.scatter(rel[:, 0], rel[:, 1], rel[:, 2], c=times, cmap="viridis", s=6) + ax.plot(rel[:, 0], rel[:, 1], rel[:, 2], color=PALETTE["muted"], lw=0.6, alpha=0.6) + ax.set_xlabel("X (nm)"); ax.set_ylabel("Y (nm)"); ax.set_zlabel("Z (nm)") + ax.set_title("COM 3D trajectory") + fig.colorbar(sc, ax=ax, shrink=0.6, pad=0.1, label="Time (ns)") + plotting.save_figure(fig, ctx.fig_path("com_3d"), dpi=ctx.config.dpi) + + return {"total_distance_nm": float(cumdist[-1]), + "net_displacement_nm": float(net[-1]), + "displacement": st.describe(net, "net_displacement_nm"), + "figure": "com_xyz"} diff --git a/mdforge/analysis/contact_map.py b/mdforge/analysis/contact_map.py new file mode 100644 index 0000000..6508a5e --- /dev/null +++ b/mdforge/analysis/contact_map.py @@ -0,0 +1,67 @@ +"""Residue-residue contact occupancy + time-resolved change map.""" + +from __future__ import annotations + +import numpy as np +import pandas as pd +from scipy.spatial.distance import cdist + +from mdforge.core.base import BaseAnalysis +from mdforge.core.system import PROTEIN_SYSTEMS +from mdforge import plotting +from mdforge import statistics as st +from mdforge.plotting import PALETTE, SEQ_CMAP, DIV_CMAP + + +class ContactMap(BaseAnalysis): + name = "contact_map" + label = "Residue-residue contact map" + category = "structure" + required_files = {"trajectory", "topology"} + supported_systems = PROTEIN_SYSTEMS + outputs = ["results/contact_occupancy.csv", "figures/contact_occupancy.png"] + default_params = {"cutoff": 8.0} + + def run(self, ctx) -> dict: + p = self.params(ctx) + plotting.set_style() + u = ctx.core_universe() + ca = u.select_atoms("protein and name CA") + resids = ca.resids; N = len(resids); cutoff = p["cutoff"] + seq = np.abs(resids[:, None] - resids[None, :]) > 2 + + occ = np.zeros((N, N)); occ1 = np.zeros((N, N)); occ2 = np.zeros((N, N)) + n = len(u.trajectory); half = n // 2 + counts = np.empty(n); times = np.empty(n) + for i, ts in enumerate(ctx.iter_frames(u, desc="[contact_map]")): + c = (cdist(ca.positions, ca.positions) < cutoff) & seq + occ += c + (occ1 if i < half else occ2)[:] += c + counts[i] = c.sum() / 2.0; times[i] = ts.time / 1000.0 + occ /= n; occ1 /= max(half, 1); occ2 /= max(n - half, 1) + diff = occ2 - occ1 + + pd.DataFrame(occ, index=resids, columns=resids).to_csv( + ctx.csv_path("contact_occupancy.csv")) + ctx.write_csv(pd.DataFrame({"time_ns": times, "n_contacts": counts}), + "contacts_vs_time.csv") + + fig, ax = plotting.new_axes(figsize=(6.4, 5.6)) + im = ax.imshow(occ, cmap=SEQ_CMAP, origin="lower", vmin=0, vmax=1, + extent=[resids[0], resids[-1], resids[0], resids[-1]]) + ax.set_xlabel("Residue"); ax.set_ylabel("Residue") + ax.set_title("Contact occupancy") + fig.colorbar(im, ax=ax, shrink=0.85, label="Occupancy") + plotting.save_figure(fig, ctx.fig_path("contact_occupancy"), dpi=ctx.config.dpi) + + fig, ax = plotting.new_axes(figsize=(6.4, 5.6)) + im = ax.imshow(diff, cmap=DIV_CMAP, origin="lower", vmin=-1, vmax=1, + extent=[resids[0], resids[-1], resids[0], resids[-1]]) + ax.set_xlabel("Residue"); ax.set_ylabel("Residue") + ax.set_title("Contact change (2nd - 1st half)") + fig.colorbar(im, ax=ax, shrink=0.85, label=r"$\Delta$ occupancy") + plotting.save_figure(fig, ctx.fig_path("contact_difference"), dpi=ctx.config.dpi) + + return {"mean_contacts": float(counts.mean()), + "n_persistent_pairs": int((occ >= 0.5).sum() // 2), + "figure": "contact_occupancy"} diff --git a/mdforge/analysis/convergence.py b/mdforge/analysis/convergence.py new file mode 100644 index 0000000..9f340c1 --- /dev/null +++ b/mdforge/analysis/convergence.py @@ -0,0 +1,95 @@ +"""Convergence diagnostics: plateau, block averaging, halves, essential-dynamics RMSIP.""" + +from __future__ import annotations + +import numpy as np +import pandas as pd +import matplotlib.pyplot as plt +from scipy import stats as sps +from MDAnalysis.analysis import align, rms + +from mdforge.core.base import BaseAnalysis +from mdforge.core.system import PROTEIN_SYSTEMS +from mdforge import plotting +from mdforge import statistics as st +from mdforge.plotting import PALETTE + + +class Convergence(BaseAnalysis): + name = "convergence" + label = "Convergence / equilibration" + category = "quality" + required_files = {"trajectory", "topology"} + supported_systems = PROTEIN_SYSTEMS + order = 200 # runs after rmsd/rog so it can reuse their CSVs + outputs = ["results/convergence_summary.csv", "figures/convergence_running.png"] + default_params = {"n_blocks": 10, "n_modes": 10} + + def _series(self, ctx): + rc, gc = ctx.csv_path("rmsd.csv"), ctx.csv_path("rog.csv") + if rc.exists() and gc.exists(): + r, g = pd.read_csv(rc), pd.read_csv(gc) + return (r["time_ns"].to_numpy(), r["rmsd_backbone_nm"].to_numpy(), + g["rg_nm"].to_numpy()) + u = ctx.core_universe() + protein = u.select_atoms("protein") + R = rms.RMSD(u, select="backbone", ref_frame=ctx.config.start or 0).run() + rmsd = R.results.rmsd[:, 2] / 10.0 + rg, t = [], [] + for ts in u.trajectory: + rg.append(protein.radius_of_gyration() / 10.0); t.append(ts.time / 1000.0) + return np.array(t), rmsd, np.array(rg) + + def _rmsip(self, ctx, k): + u = ctx.core_universe() + align.AlignTraj(u, u, select="protein and name CA", + ref_frame=ctx.config.start or 0, in_memory=True).run() + ca = u.select_atoms("protein and name CA") + X = np.array([ca.positions.copy().ravel() for _ in u.trajectory]) / 10.0 + half = len(X) // 2 + + def modes(Xh): + Xc = Xh - Xh.mean(axis=0) + ev, evec = np.linalg.eigh((Xc.T @ Xc) / (len(Xh) - 1)) + return evec[:, ::-1][:, :k] + A, B = modes(X[:half]), modes(X[half:]) + ov = A.T @ B + return float(np.sqrt((ov ** 2).sum() / k)), ov + + def run(self, ctx) -> dict: + p = self.params(ctx) + plotting.set_style() + times, rmsd, rg = self._series(ctx) + pr, pg = st.plateau_detection(times, rmsd), st.plateau_detection(times, rg) + half = len(times) // 2 + t_rmsd = sps.ttest_ind(rmsd[:half], rmsd[half:], equal_var=False) + sem = st.block_average_sem(rmsd, p["n_blocks"]) + rmsip, ov = self._rmsip(ctx, p["n_modes"]) + + fig, (a1, a2) = plt.subplots(2, 1, figsize=(7.4, 6.0), sharex=True) + a1.plot(times, rmsd, color=PALETTE["primary"], alpha=0.3, lw=0.8) + a1.plot(times, st.running_mean(rmsd), color=PALETTE["primary"], lw=2.2) + a1.set_ylabel("Backbone RMSD (nm)"); a1.set_title("Cumulative running averages") + a2.plot(times, rg, color=PALETTE["secondary"], alpha=0.3, lw=0.8) + a2.plot(times, st.running_mean(rg), color=PALETTE["secondary"], lw=2.2) + a2.set_ylabel(r"$R_g$ (nm)"); a2.set_xlabel("Time (ns)") + for a in (a1, a2): + a.spines["top"].set_visible(False); a.spines["right"].set_visible(False) + plotting.save_figure(fig, ctx.fig_path("convergence_running"), dpi=ctx.config.dpi) + + fig, ax = plotting.new_axes(figsize=(5.6, 4.8)) + im = ax.imshow(np.abs(ov), cmap="viridis", origin="lower", vmin=0, vmax=1, + extent=[1, ov.shape[1], 1, ov.shape[0]]) + ax.set_xlabel("PC (2nd half)"); ax.set_ylabel("PC (1st half)") + ax.set_title(f"Essential-dynamics overlap (RMSIP = {rmsip:.2f})") + fig.colorbar(im, ax=ax, label="|inner product|") + plotting.save_figure(fig, ctx.fig_path("essential_dynamics_overlap"), dpi=ctx.config.dpi) + + summ = {"rmsd_converged": pr["converged"], + "rmsd_tail_slope_nm_per_ns": pr["tail_slope"], + "rmsd_drift_2nd_minus_1st_nm": pr["drift_second_minus_first"], + "rg_drift_2nd_minus_1st_nm": pg["drift_second_minus_first"], + "rmsd_block_sem_nm": sem, "rmsd_halves_ttest_p": float(t_rmsd.pvalue), + "essential_dynamics_rmsip": rmsip, "figure": "convergence_running"} + pd.DataFrame([summ]).to_csv(ctx.csv_path("convergence_summary.csv"), index=False) + return summ diff --git a/mdforge/analysis/dccm.py b/mdforge/analysis/dccm.py new file mode 100644 index 0000000..b95dc10 --- /dev/null +++ b/mdforge/analysis/dccm.py @@ -0,0 +1,49 @@ +"""Dynamic cross-correlation matrix of Cα fluctuations.""" + +from __future__ import annotations + +import numpy as np +import pandas as pd +from MDAnalysis.analysis import align + +from mdforge.core.base import BaseAnalysis +from mdforge.core.system import PROTEIN_SYSTEMS +from mdforge import plotting +from mdforge.plotting import DIV_CMAP + + +class DCCM(BaseAnalysis): + name = "dccm" + label = "Dynamic cross-correlation matrix" + category = "dynamics" + required_files = {"trajectory", "topology"} + supported_systems = PROTEIN_SYSTEMS + outputs = ["results/dccm.csv", "figures/dccm.png"] + + def run(self, ctx) -> dict: + plotting.set_style() + u = ctx.core_universe() + align.AlignTraj(u, u, select="protein and name CA", + ref_frame=ctx.config.start or 0, in_memory=True).run() + ca = u.select_atoms("protein and name CA") + resids = ca.resids + coords = np.array([ca.positions.copy() for _ in u.trajectory]) / 10.0 + disp = coords - coords.mean(axis=0) + cov = np.einsum("tia,tja->ij", disp, disp) / disp.shape[0] + norm = np.sqrt(np.diag(cov)) + dccm = np.clip(cov / np.outer(norm, norm), -1, 1) + + pd.DataFrame(dccm, index=resids, columns=resids).to_csv(ctx.csv_path("dccm.csv")) + fig, ax = plotting.new_axes(figsize=(6.6, 5.6)) + im = ax.imshow(dccm, cmap=DIV_CMAP, vmin=-1, vmax=1, origin="lower", + extent=[resids[0], resids[-1], resids[0], resids[-1]]) + ax.set_xlabel("Residue"); ax.set_ylabel("Residue") + ax.set_title(r"Dynamic cross-correlation (C$\alpha$)") + fig.colorbar(im, ax=ax, shrink=0.85, label="Correlation") + plotting.save_figure(fig, ctx.fig_path("dccm"), dpi=ctx.config.dpi) + + off = dccm[~np.eye(len(resids), dtype=bool)] + return {"mean_correlation": float(off.mean()), + "frac_strong_positive": float((off > 0.5).mean()), + "frac_strong_negative": float((off < -0.5).mean()), + "figure": "dccm"} diff --git a/mdforge/analysis/energies.py b/mdforge/analysis/energies.py new file mode 100644 index 0000000..3218aba --- /dev/null +++ b/mdforge/analysis/energies.py @@ -0,0 +1,78 @@ +"""Thermodynamic observables from the GROMACS .edr (panedr): E, T, P, density.""" + +from __future__ import annotations + +from pathlib import Path + +import numpy as np +import pandas as pd + +from mdforge.core.base import BaseAnalysis +from mdforge import plotting +from mdforge import statistics as st +from mdforge.plotting import PALETTE + +_WANTED = { + "Potential": ["Potential"], "Kinetic": ["Kinetic En.", "Kinetic"], + "Total Energy": ["Total Energy"], "Temperature": ["Temperature"], + "Pressure": ["Pressure"], "Density": ["Density"], +} +_LABELS = {"Potential": "Potential (kJ/mol)", "Kinetic": "Kinetic (kJ/mol)", + "Total Energy": "Total energy (kJ/mol)", "Temperature": "Temperature (K)", + "Pressure": "Pressure (bar)", "Density": "Density (kg m$^{-3}$)"} + + +class Energies(BaseAnalysis): + name = "energies" + label = "Thermodynamic stability observables" + category = "thermodynamics" + required_files = {"energy"} + supported_systems = {"*"} + outputs = ["results/energies.csv", "tables/energy_summary.csv", "figures/energies.png"] + + def run(self, ctx) -> dict: + import panedr + plotting.set_style() + edr = ctx.fileset.energy or ctx.config.energy + if edr is None or not Path(edr).exists(): + return {"note": "edr not available"} + edf = panedr.edr_to_df(str(edr)) + time_ns = edf["Time"].to_numpy() / 1000.0 + + def find(terms): + for t in terms: + for c in edf.columns: + if t.lower() in c.lower(): + return c + return None + + series = {"time_ns": time_ns}; found = {} + for name, terms in _WANTED.items(): + col = find(terms) + if col is not None: + series[name] = edf[col].to_numpy(); found[name] = col + out = pd.DataFrame(series) + ctx.write_csv(out, "energies.csv") + st.summary_frame({k: out[k] for k in found}).to_csv( + ctx.table_path("energy_summary.csv")) + + keys = list(found); ncols = 3; nrows = int(np.ceil(len(keys)/ncols)) + fig, axes = plotting.style.plt.subplots(nrows, ncols, figsize=(4.2*ncols, 3.2*nrows)) + axes = np.atleast_1d(axes).ravel() + colors = [PALETTE["primary"], PALETTE["secondary"], PALETTE["green"], + PALETTE["purple"], PALETTE["helix"], PALETTE["accent"]] + for ax, k, c in zip(axes, keys, colors): + y = out[k].to_numpy() + ax.plot(time_ns, y, color=c, alpha=0.3, lw=0.8) + ax.plot(time_ns, st.running_mean(y), color=c, lw=2.0) + ax.set_xlabel("Time (ns)"); ax.set_ylabel(_LABELS.get(k, k), fontsize=11) + ax.spines["top"].set_visible(False); ax.spines["right"].set_visible(False) + for ax in axes[len(keys):]: + ax.set_visible(False) + fig.suptitle("Thermodynamic stability observables", fontweight="bold") + plotting.save_figure(fig, ctx.fig_path("energies"), dpi=ctx.config.dpi) + + summary = {k: {"mean": float(out[k].mean()), "std": float(out[k].std(ddof=1))} + for k in found} + summary["figure"] = "energies" + return summary diff --git a/mdforge/analysis/hbonds.py b/mdforge/analysis/hbonds.py new file mode 100644 index 0000000..ff3fcb4 --- /dev/null +++ b/mdforge/analysis/hbonds.py @@ -0,0 +1,77 @@ +"""Hydrogen bonds: protein-protein and protein-solvent (streamed, shell-restricted).""" + +from __future__ import annotations + +import numpy as np +import pandas as pd +import matplotlib.pyplot as plt +from MDAnalysis.analysis.hydrogenbonds import HydrogenBondAnalysis as HBA + +from mdforge.core.base import BaseAnalysis +from mdforge import plotting +from mdforge import statistics as st +from mdforge.plotting import PALETTE + +WATER = "resname TIP3 SOL WAT HOH SPC TIP4" + + +class HydrogenBonds(BaseAnalysis): + name = "hbonds" + label = "Hydrogen bonds" + category = "interactions" + required_files = {"trajectory", "topology"} + supported_systems = {"*"} + default_params = {"stride": 20, "shell": 4.5, "window": 10} + outputs = ["results/hbonds_timeseries.csv", "figures/hbonds_timeseries.png"] + + def run(self, ctx) -> dict: + p = self.params(ctx) + plotting.set_style() + u = ctx.full_universe() + stride = p["stride"] + + hba_pp = HBA(universe=u) + hba_pp.hydrogens_sel = hba_pp.guess_hydrogens("protein") + hba_pp.acceptors_sel = hba_pp.guess_acceptors("protein") + hba_pp.run(step=stride) + counts_pp = hba_pp.count_by_time() + n = len(counts_pp) + dt = u.trajectory.dt + times = (u.trajectory[0].time + np.arange(n) * stride * dt) / 1000.0 + + counts_ps = np.zeros(n) + if ctx.system.flags.get("has_water"): + shell = p["shell"] + wo = f"({WATER}) and name OH2 OW OW1 O" + wh = f"({WATER}) and name H1 H2 HW1 HW2" + pw = HBA(universe=u, update_selections=True) + pw.hydrogens_sel = pw.guess_hydrogens("protein") + pw.acceptors_sel = f"({wo}) and around {shell} protein" + pw.run(step=stride) + wp = HBA(universe=u, update_selections=True) + wp.hydrogens_sel = f"({wh}) and around {shell} protein" + wp.acceptors_sel = wp.guess_acceptors("protein") + wp.run(step=stride) + counts_ps = (pw.count_by_time() + wp.count_by_time())[:n] + + ctx.write_csv(pd.DataFrame({"time_ns": times, + "hbonds_protein_protein": counts_pp, + "hbonds_protein_solvent": counts_ps}), + "hbonds_timeseries.csv") + + fig, (a1, a2) = plt.subplots(2, 1, figsize=(7.6, 6.4), sharex=True) + a1.plot(times, counts_pp, color=PALETTE["primary"], alpha=0.4, lw=1.0) + a1.plot(times, st.moving_average(counts_pp, p["window"]), + color=PALETTE["primary"], lw=2.2) + a1.set_ylabel("Protein-protein\nH-bonds"); a1.set_title("Hydrogen bonds over time") + a2.plot(times, counts_ps, color=PALETTE["secondary"], alpha=0.4, lw=1.0) + a2.plot(times, st.moving_average(counts_ps, p["window"]), + color=PALETTE["secondary"], lw=2.2) + a2.set_ylabel("Protein-solvent\nH-bonds"); a2.set_xlabel("Time (ns)") + for a in (a1, a2): + a.spines["top"].set_visible(False); a.spines["right"].set_visible(False) + plotting.save_figure(fig, ctx.fig_path("hbonds_timeseries"), dpi=ctx.config.dpi) + + return {"protein_protein": st.describe(counts_pp, "hbonds_pp"), + "protein_solvent": st.describe(counts_ps, "hbonds_ps"), + "figure": "hbonds_timeseries"} diff --git a/mdforge/analysis/native_contacts.py b/mdforge/analysis/native_contacts.py new file mode 100644 index 0000000..9a95560 --- /dev/null +++ b/mdforge/analysis/native_contacts.py @@ -0,0 +1,61 @@ +"""Fraction of native contacts Q(t) relative to the reference frame.""" + +from __future__ import annotations + +import numpy as np +import pandas as pd +from scipy.spatial.distance import cdist + +from mdforge.core.base import BaseAnalysis +from mdforge.core.system import PROTEIN_SYSTEMS +from mdforge import plotting +from mdforge import statistics as st +from mdforge.plotting import PALETTE, SEQ_CMAP + + +class NativeContacts(BaseAnalysis): + name = "native_contacts" + label = "Native-contact fraction Q" + category = "structure" + required_files = {"trajectory", "topology"} + supported_systems = PROTEIN_SYSTEMS + outputs = ["results/native_contacts.csv", "figures/native_contacts_Q.png"] + default_params = {"cutoff": 8.0, "seq_sep": 4, "window": 20} + + def run(self, ctx) -> dict: + p = self.params(ctx) + plotting.set_style() + u = ctx.core_universe() + ca = u.select_atoms("protein and name CA") + resids = ca.resids + R0 = p["cutoff"] + u.trajectory[ctx.config.start or 0] + D0 = cdist(ca.positions, ca.positions) + seq = np.abs(resids[:, None] - resids[None, :]) >= p["seq_sep"] + native = np.triu((D0 < R0) & seq) + n_native = int(native.sum()) or 1 + + n = len(u.trajectory) + Q = np.empty(n); times = np.empty(n) + for i, ts in enumerate(ctx.iter_frames(u, desc="[native_contacts]")): + D = cdist(ca.positions, ca.positions) + Q[i] = ((D < R0) & native).sum() / n_native + times[i] = ts.time / 1000.0 + + ctx.write_csv(pd.DataFrame({"time_ns": times, "Q": Q, + "Q_movavg": st.moving_average(Q, p["window"])}), + "native_contacts.csv") + fig, ax = plotting.new_axes() + ax.plot(times, Q, color=PALETTE["primary"], alpha=0.35, lw=1.0) + ax.plot(times, st.moving_average(Q, p["window"]), color=PALETTE["primary"], lw=2.2) + plateau = float(Q[len(Q)//2:].mean()) + ax.axhline(plateau, ls="--", color=PALETTE["accent"], lw=1.1, + label=f"plateau = {plateau:.2f}") + ax.set_xlabel("Time (ns)"); ax.set_ylabel("Fraction of native contacts, Q") + ax.set_ylim(0, 1.02); ax.set_title(self.label); ax.legend(loc="best") + plotting.save_figure(fig, ctx.fig_path("native_contacts_Q"), dpi=ctx.config.dpi) + + return {"n_native_contacts": n_native, "Q_final_plateau": plateau, + "Q": st.describe(Q, "Q"), + "convergence": st.plateau_detection(times, Q), + "figure": "native_contacts_Q"} diff --git a/mdforge/analysis/pca.py b/mdforge/analysis/pca.py new file mode 100644 index 0000000..94abb58 --- /dev/null +++ b/mdforge/analysis/pca.py @@ -0,0 +1,78 @@ +"""PCA (essential dynamics) on Cα + PC1/PC2 free-energy landscape.""" + +from __future__ import annotations + +import numpy as np +import pandas as pd +import matplotlib.pyplot as plt +from MDAnalysis.analysis import align + +from mdforge.core.base import BaseAnalysis +from mdforge.core.system import PROTEIN_SYSTEMS +from mdforge import plotting +from mdforge.plotting import PALETTE, DIV_CMAP, FEL_CMAP + +KB_KJ = 0.00831446261815324 + + +class PCA(BaseAnalysis): + name = "pca" + label = "PCA + free-energy landscape" + category = "dynamics" + required_files = {"trajectory", "topology"} + supported_systems = PROTEIN_SYSTEMS + outputs = ["results/pca_projections.csv", "figures/pca_fel.png"] + default_params = {"temperature": 310.0, "bins": 80} + + def run(self, ctx) -> dict: + p = self.params(ctx) + plotting.set_style() + u = ctx.core_universe() + align.AlignTraj(u, u, select="protein and name CA", + ref_frame=ctx.config.start or 0, in_memory=True).run() + ca = u.select_atoms("protein and name CA") + X, times = [], [] + for ts in u.trajectory: + X.append(ca.positions.copy().ravel()); times.append(ts.time / 1000.0) + X = np.asarray(X) / 10.0; times = np.asarray(times) + Xc = X - X.mean(axis=0); n = len(X) + evals, evecs = np.linalg.eigh((Xc.T @ Xc) / (n - 1)) + order = np.argsort(evals)[::-1] + evals, evecs = np.clip(evals[order], 0, None), evecs[:, order] + var = evals / evals.sum(); cum = np.cumsum(var) + proj = Xc @ evecs[:, :5] + + pd.DataFrame({"component": np.arange(1, len(evals)+1), "eigenvalue_nm2": evals, + "variance_explained": var, "cumulative_variance": cum}).to_csv( + ctx.csv_path("pca_eigenvalues.csv"), index=False) + pd.DataFrame({"time_ns": times, **{f"PC{k+1}": proj[:, k] for k in range(5)}}).to_csv( + ctx.csv_path("pca_projections.csv"), index=False) + + # variance + fig, ax = plotting.new_axes() + kk = min(20, len(evals)) + ax.bar(np.arange(1, kk+1), var[:kk]*100, color=PALETTE["primary"], alpha=0.85) + ax2 = ax.twinx(); ax2.plot(np.arange(1, kk+1), cum[:kk]*100, "-o", + color=PALETTE["secondary"], ms=4) + ax2.set_ylabel("Cumulative (%)"); ax2.set_ylim(0, 105); ax2.spines["top"].set_visible(False) + ax.set_xlabel("Principal component"); ax.set_ylabel("Variance (%)") + ax.set_title("PCA variance explained") + plotting.save_figure(fig, ctx.fig_path("pca_variance"), dpi=ctx.config.dpi) + + # FEL + H, xe, ye = np.histogram2d(proj[:, 0], proj[:, 1], bins=p["bins"], density=True) + P = (H / H.sum()).T + with np.errstate(divide="ignore"): + G = -KB_KJ * p["temperature"] * np.log(P) + G[~np.isfinite(G)] = np.nan; G -= np.nanmin(G) + xc = 0.5*(xe[:-1]+xe[1:]); yc = 0.5*(ye[:-1]+ye[1:]) + fig, ax = plotting.new_axes(figsize=(6.8, 5.4)) + im = ax.pcolormesh(xc, yc, G, cmap=FEL_CMAP, shading="auto") + fig.colorbar(im, ax=ax, label=r"$\Delta G$ (kJ/mol)") + ax.set_xlabel("PC1 (nm)"); ax.set_ylabel("PC2 (nm)") + ax.set_title(f"Free-energy landscape (T = {p['temperature']:.0f} K)") + plotting.save_figure(fig, ctx.fig_path("pca_fel"), dpi=ctx.config.dpi) + + return {"pc1_variance": float(var[0]), "pc2_variance": float(var[1]), + "n_pc_for_80pct": int(np.searchsorted(cum, 0.80) + 1), + "figure": "pca_fel"} diff --git a/mdforge/analysis/rin.py b/mdforge/analysis/rin.py new file mode 100644 index 0000000..8804e8b --- /dev/null +++ b/mdforge/analysis/rin.py @@ -0,0 +1,75 @@ +"""Residue interaction network (RIN) + centrality / hub analysis.""" + +from __future__ import annotations + +import numpy as np +import pandas as pd +import networkx as nx +from scipy.spatial.distance import cdist + +from mdforge.core.base import BaseAnalysis +from mdforge.core.system import PROTEIN_SYSTEMS +from mdforge import plotting +from mdforge.plotting import PALETTE + +_OCC = 0.5 + + +class ResidueInteractionNetwork(BaseAnalysis): + name = "rin" + label = "Residue interaction network" + category = "structure" + required_files = {"trajectory", "topology"} + supported_systems = PROTEIN_SYSTEMS + order = 150 # after contact_map (reuses its occupancy CSV when present) + outputs = ["results/rin_node_metrics.csv", "figures/rin_network.png"] + default_params = {"cutoff": 8.0} + + def _adjacency(self, ctx): + occ_csv = ctx.csv_path("contact_occupancy.csv") + if occ_csv.exists(): + df = pd.read_csv(occ_csv, index_col=0) + return df.columns.astype(int).to_numpy(), df.to_numpy() + u = ctx.core_universe() + ca = u.select_atoms("protein and name CA") + resids = ca.resids + u.trajectory[ctx.config.start or 0] + D = cdist(ca.positions, ca.positions) + seq = np.abs(resids[:, None] - resids[None, :]) > 2 + return resids, ((D < self.default_params["cutoff"]) & seq).astype(float) + + def run(self, ctx) -> dict: + plotting.set_style() + resids, occ = self._adjacency(ctx) + N = len(resids) + G = nx.Graph(); G.add_nodes_from(resids.tolist()) + for i in range(N): + for j in range(i + 1, N): + if occ[i, j] >= _OCC: + G.add_edge(int(resids[i]), int(resids[j]), weight=float(occ[i, j])) + + deg = dict(G.degree()) + betw = nx.betweenness_centrality(G, weight="weight") + clos = nx.closeness_centrality(G) + node_df = pd.DataFrame({"resid": list(G.nodes()), + "degree": [deg[n] for n in G.nodes()], + "betweenness": [betw[n] for n in G.nodes()], + "closeness": [clos[n] for n in G.nodes()]}).sort_values( + "betweenness", ascending=False).reset_index(drop=True) + ctx.write_csv(node_df, "rin_node_metrics.csv") + + fig, ax = plotting.new_axes(figsize=(7.2, 6.4)) + if G.number_of_edges(): + pos = nx.spring_layout(G, seed=42, weight="weight", k=0.3) + nd = np.array([deg[n] for n in G.nodes()]) + nb = np.array([betw[n] for n in G.nodes()]) + nx.draw_networkx_edges(G, pos, ax=ax, alpha=0.25, width=0.8) + nodes = nx.draw_networkx_nodes(G, pos, ax=ax, node_size=30 + 12*nd, + node_color=nb, cmap="plasma") + fig.colorbar(nodes, ax=ax, label="Betweenness", shrink=0.8) + ax.set_axis_off(); ax.set_title("Residue interaction network") + plotting.save_figure(fig, ctx.fig_path("rin_network"), dpi=ctx.config.dpi) + + return {"n_nodes": G.number_of_nodes(), "n_edges": G.number_of_edges(), + "top_hub_resid": int(node_df.iloc[0]["resid"]) if len(node_df) else None, + "graph_density": float(nx.density(G)), "figure": "rin_network"} diff --git a/mdforge/analysis/salt_bridges.py b/mdforge/analysis/salt_bridges.py new file mode 100644 index 0000000..7331a3b --- /dev/null +++ b/mdforge/analysis/salt_bridges.py @@ -0,0 +1,71 @@ +"""Salt-bridge analysis (acidic O ↔ basic N), count over time + occupancy.""" + +from __future__ import annotations + +from collections import defaultdict + +import numpy as np +import pandas as pd +from MDAnalysis.analysis.distances import distance_array + +from mdforge.core.base import BaseAnalysis +from mdforge.core.system import PROTEIN_SYSTEMS +from mdforge import plotting +from mdforge import statistics as st +from mdforge.plotting import PALETTE + +_ACIDIC = "(resname ASP and name OD1 OD2) or (resname GLU and name OE1 OE2)" +_BASIC = ("(resname LYS and name NZ) or (resname ARG and name NH1 NH2 NE) or " + "(resname HIS HSD HSE HSP HISD HISE and name ND1 NE2)") + + +class SaltBridges(BaseAnalysis): + name = "salt_bridges" + label = "Salt-bridge analysis" + category = "interactions" + required_files = {"trajectory", "topology"} + supported_systems = PROTEIN_SYSTEMS + outputs = ["results/salt_bridges_timeseries.csv", "figures/salt_bridges_timeseries.png"] + default_params = {"cutoff": 4.0, "window": 20} + + def run(self, ctx) -> dict: + p = self.params(ctx) + plotting.set_style() + u = ctx.core_universe() + acidic = u.select_atoms(_ACIDIC); basic = u.select_atoms(_BASIC) + if acidic.n_atoms == 0 or basic.n_atoms == 0: + return {"n_salt_bridges_mean": 0.0, "note": "no chargeable groups"} + a_res = np.array([f"{a.resname}{a.resid}" for a in acidic]) + b_res = np.array([f"{b.resname}{b.resid}" for b in basic]) + n = len(u.trajectory) + counts = np.empty(n); times = np.empty(n); pair_frames = defaultdict(int) + for i, ts in enumerate(ctx.iter_frames(u, desc="[salt_bridges]")): + close = distance_array(acidic.positions, basic.positions) < p["cutoff"] + ai, bi = np.where(close) + pairs = set(zip(a_res[ai], b_res[bi])) + for pr in pairs: + pair_frames[pr] += 1 + counts[i] = len(pairs); times[i] = ts.time / 1000.0 + + ctx.write_csv(pd.DataFrame({"time_ns": times, "n_salt_bridges": counts, + "movavg": st.moving_average(counts, p["window"])}), + "salt_bridges_timeseries.csv") + occ = pd.DataFrame([{"acidic": a, "basic": b, "occupancy_percent": 100*c/n} + for (a, b), c in pair_frames.items()]) + if len(occ): + occ = occ.sort_values("occupancy_percent", ascending=False) + ctx.write_csv(occ, "salt_bridges_occupancy.csv") + + fig, ax = plotting.new_axes() + ax.plot(times, counts, color=PALETTE["primary"], alpha=0.35, lw=1.0) + ax.plot(times, st.moving_average(counts, p["window"]), + color=PALETTE["primary"], lw=2.2, label="MA") + ax.axhline(counts.mean(), ls="--", color=PALETTE["accent"], lw=1.1, + label=f"mean = {counts.mean():.1f}") + ax.set_xlabel("Time (ns)"); ax.set_ylabel("Number of salt bridges") + ax.set_title(self.label); ax.legend(loc="best") + plotting.save_figure(fig, ctx.fig_path("salt_bridges_timeseries"), dpi=ctx.config.dpi) + + return {"n_salt_bridges_mean": float(counts.mean()), + "n_stable_bridges": int((occ["occupancy_percent"] >= 50).sum()) if len(occ) else 0, + "figure": "salt_bridges_timeseries"} diff --git a/mdforge/analysis/secondary_structure.py b/mdforge/analysis/secondary_structure.py new file mode 100644 index 0000000..cf900ff --- /dev/null +++ b/mdforge/analysis/secondary_structure.py @@ -0,0 +1,66 @@ +"""DSSP secondary-structure timeline and content (mdtraj, no external binary).""" + +from __future__ import annotations + +import numpy as np +import pandas as pd +import mdtraj as md +from matplotlib.colors import ListedColormap, BoundaryNorm + +from mdforge.core.base import BaseAnalysis +from mdforge.core.system import PROTEIN_SYSTEMS +from mdforge import plotting +from mdforge import statistics as st +from mdforge.plotting import PALETTE + +_MAP = {"C": 0, "H": 1, "E": 2} +_LABELS = ["Coil", "Helix", "Sheet"] +_COLORS = [PALETTE["coil"], PALETTE["helix"], PALETTE["sheet"]] + + +class SecondaryStructure(BaseAnalysis): + name = "secondary_structure" + label = "DSSP secondary structure" + category = "structure" + required_files = {"trajectory", "topology"} + supported_systems = PROTEIN_SYSTEMS + outputs = ["results/ss_fractions.csv", "figures/ss_timeline.png"] + + def run(self, ctx) -> dict: + plotting.set_style() + ctx.core_universe() + traj = md.load(str(ctx.config.data_dir / "core.xtc"), + top=str(ctx.config.data_dir / "core.pdb")) + times = traj.time / 1000.0 + dssp = md.compute_dssp(traj, simplified=True) + resids = np.array([r.resSeq for r in traj.topology.residues]) + fh = (dssp == "H").mean(axis=1) * 100 + fe = (dssp == "E").mean(axis=1) * 100 + fc = (dssp == "C").mean(axis=1) * 100 + ctx.write_csv(pd.DataFrame({"time_ns": times, "helix_pct": fh, + "sheet_pct": fe, "coil_pct": fc}), "ss_fractions.csv") + coded = np.vectorize(_MAP.get)(dssp).astype(np.int8) + pd.DataFrame(coded, index=np.round(times, 3), columns=resids).to_csv( + ctx.csv_path("ss_timeline.csv")) + + fig, ax = plotting.new_axes() + for y, lab, c in ((fh, "Helix", PALETTE["helix"]), (fe, "Sheet", PALETTE["sheet"]), + (fc, "Coil", PALETTE["coil"])): + ax.plot(times, st.moving_average(y, 20), color=c, lw=2.0, label=lab) + ax.set_xlabel("Time (ns)"); ax.set_ylabel("Content (%)") + ax.set_title("Secondary-structure content"); ax.legend(ncol=3, loc="best") + plotting.save_figure(fig, ctx.fig_path("ss_fractions"), dpi=ctx.config.dpi) + + cmap = ListedColormap(_COLORS); norm = BoundaryNorm([-0.5, 0.5, 1.5, 2.5], cmap.N) + fig, ax = plotting.new_axes(figsize=(8.4, 4.8)) + im = ax.imshow(coded.T, aspect="auto", origin="lower", cmap=cmap, norm=norm, + extent=[times[0], times[-1], resids[0], resids[-1]], + interpolation="nearest") + ax.set_xlabel("Time (ns)"); ax.set_ylabel("Residue") + ax.set_title("Secondary-structure timeline") + cb = fig.colorbar(im, ax=ax, ticks=[0, 1, 2], shrink=0.85) + cb.ax.set_yticklabels(_LABELS) + plotting.save_figure(fig, ctx.fig_path("ss_timeline"), dpi=ctx.config.dpi) + + return {"mean_helix_pct": float(fh.mean()), "mean_sheet_pct": float(fe.mean()), + "mean_coil_pct": float(fc.mean()), "figure": "ss_timeline"} diff --git a/mdforge/analysis/statistics_summary.py b/mdforge/analysis/statistics_summary.py new file mode 100644 index 0000000..1429471 --- /dev/null +++ b/mdforge/analysis/statistics_summary.py @@ -0,0 +1,81 @@ +"""Cross-metric summary table + Pearson/Spearman correlations (runs last).""" + +from __future__ import annotations + +import numpy as np +import pandas as pd + +from mdforge.core.base import BaseAnalysis +from mdforge import plotting +from mdforge import statistics as st +from mdforge.plotting import DIV_CMAP + +# (csv file, time col, value col, label) +_SOURCES = [ + ("rmsd.csv", "time_ns", "rmsd_backbone_nm", "RMSD"), + ("rog.csv", "time_ns", "rg_nm", "Rg"), + ("sasa_total.csv", "time_ns", "sasa_total_nm2", "SASA"), + ("hbonds_timeseries.csv", "time_ns", "hbonds_protein_protein", "HB_intra"), + ("energies.csv", "time_ns", "Potential", "E_pot"), + ("com.csv", "time_ns", "net_displacement_nm", "COM_disp"), + ("native_contacts.csv", "time_ns", "Q", "Q_native"), + ("ss_fractions.csv", "time_ns", "helix_pct", "Helix"), +] + + +class StatisticsSummary(BaseAnalysis): + name = "statistics_summary" + label = "Cross-metric statistics & correlations" + category = "quality" + required_files = {"trajectory", "topology"} + supported_systems = {"*"} + order = 300 # last: reads other analyses' CSVs + outputs = ["tables/summary_statistics.csv", "figures/correlation_pearson.png"] + + def run(self, ctx) -> dict: + plotting.set_style() + master = None + for f in ("rog.csv", "rmsd.csv"): + if ctx.csv_path(f).exists(): + master = pd.read_csv(ctx.csv_path(f))["time_ns"].to_numpy(); break + if master is None: + return {"note": "no base time-series found"} + + merged = {"time_ns": master} + for fname, tcol, vcol, label in _SOURCES: + pth = ctx.csv_path(fname) + if not pth.exists(): + continue + df = pd.read_csv(pth) + if tcol in df and vcol in df: + good = np.isfinite(df[tcol]) & np.isfinite(df[vcol]) + if good.sum() >= 2: + merged[label] = np.interp(master, df[tcol][good], df[vcol][good]) + mdf = pd.DataFrame(merged) + cols = [c for c in mdf.columns if c != "time_ns"] + ctx.write_csv(mdf, "metrics_merged.csv") + st.summary_frame({c: mdf[c] for c in cols}).to_csv( + ctx.table_path("summary_statistics.csv")) + + pearson, spearman = st.correlation_matrices(mdf[cols]) + pearson.to_csv(ctx.csv_path("correlation_pearson.csv")) + spearman.to_csv(ctx.csv_path("correlation_spearman.csv")) + for corr, nm, title in ((pearson, "correlation_pearson", "Pearson"), + (spearman, "correlation_spearman", "Spearman")): + fig, ax = plotting.new_axes(figsize=(6.6, 5.6)) + im = ax.imshow(corr.to_numpy(), cmap=DIV_CMAP, vmin=-1, vmax=1) + ax.set_xticks(range(len(cols))); ax.set_yticks(range(len(cols))) + ax.set_xticklabels(cols, rotation=45, ha="right"); ax.set_yticklabels(cols) + for i in range(len(cols)): + for j in range(len(cols)): + ax.text(j, i, f"{corr.iloc[i, j]:.2f}", ha="center", va="center", + fontsize=8) + ax.set_title(f"{title} correlation") + fig.colorbar(im, ax=ax, shrink=0.85, label="correlation") + plotting.save_figure(fig, ctx.fig_path(nm), dpi=ctx.config.dpi) + + notable = [(cols[i], cols[j], round(float(pearson.iloc[i, j]), 2)) + for i in range(len(cols)) for j in range(i+1, len(cols)) + if abs(pearson.iloc[i, j]) > 0.6] + return {"metrics": cols, "notable_correlations": notable, + "figure": "correlation_pearson"} diff --git a/mdforge/analysis/structural.py b/mdforge/analysis/structural.py new file mode 100644 index 0000000..03d1108 --- /dev/null +++ b/mdforge/analysis/structural.py @@ -0,0 +1,80 @@ +"""Global shape descriptors: end-to-end, Dmax, gyration-tensor anisotropy, volume.""" + +from __future__ import annotations + +import numpy as np +import pandas as pd +import matplotlib.pyplot as plt +from scipy.spatial import ConvexHull +from scipy.spatial.distance import pdist + +from mdforge.core.base import BaseAnalysis +from mdforge.core.system import PROTEIN_SYSTEMS +from mdforge import plotting +from mdforge import statistics as st +from mdforge.plotting import PALETTE + + +def _gyration(pos): + r = pos - pos.mean(axis=0) + S = (r[:, :, None] * r[:, None, :]).mean(axis=0) + ev = np.sort(np.linalg.eigvalsh(S)) + l1, l2, l3 = ev + tr = ev.sum() + rg = float(np.sqrt(tr)) + asph = float(l3 - 0.5 * (l1 + l2)) + kappa2 = float(1 - 3 * (l1*l2 + l2*l3 + l3*l1) / (tr*tr)) if tr > 0 else 0.0 + return rg, asph, kappa2 + + +class StructuralDescriptors(BaseAnalysis): + name = "structural" + label = "Global shape descriptors" + category = "structure" + required_files = {"trajectory", "topology"} + supported_systems = PROTEIN_SYSTEMS + outputs = ["results/structural.csv", "figures/structural_descriptors.png"] + + def run(self, ctx) -> dict: + plotting.set_style() + u = ctx.core_universe() + protein = u.select_atoms("protein") + ca = u.select_atoms("protein and name CA") + n = len(u.trajectory) + e2e = np.empty(n); dmax = np.empty(n); rg = np.empty(n) + asph = np.empty(n); kappa2 = np.empty(n); vol = np.empty(n); times = np.empty(n) + for i, ts in enumerate(ctx.iter_frames(u, desc="[structural]")): + times[i] = ts.time / 1000.0 + cap = ca.positions / 10.0; ap = protein.positions / 10.0 + e2e[i] = np.linalg.norm(cap[0] - cap[-1]) + dmax[i] = pdist(cap).max() + rg[i], asph[i], kappa2[i] = _gyration(ap) + try: + vol[i] = ConvexHull(ap, qhull_options="QJ").volume + except Exception: + vol[i] = np.nan + + df = pd.DataFrame({"time_ns": times, "end_to_end_nm": e2e, "dmax_nm": dmax, + "rg_tensor_nm": rg, "asphericity_nm2": asph, + "rel_shape_anisotropy": kappa2, "hull_volume_nm3": vol, + "elongation": dmax / (2 * rg)}) + ctx.write_csv(df, "structural.csv") + + fig, axes = plt.subplots(2, 2, figsize=(9.6, 7.0)) + for ax, (col, lab, c) in zip(axes.ravel(), [ + ("end_to_end_nm", "End-to-end (nm)", PALETTE["primary"]), + ("dmax_nm", r"$D_{max}$ (nm)", PALETTE["secondary"]), + ("rel_shape_anisotropy", r"$\kappa^2$", PALETTE["purple"]), + ("hull_volume_nm3", "Volume (nm$^3$)", PALETTE["green"])]): + y = df[col].to_numpy() + ax.plot(times, y, color=c, alpha=0.4, lw=1.0) + ax.plot(times, st.moving_average(y, 20), color=c, lw=2.0) + ax.set_xlabel("Time (ns)"); ax.set_ylabel(lab) + ax.spines["top"].set_visible(False); ax.spines["right"].set_visible(False) + fig.suptitle("Global shape descriptors", fontweight="bold") + plotting.save_figure(fig, ctx.fig_path("structural_descriptors"), dpi=ctx.config.dpi) + + return {"dmax": st.describe(dmax, "dmax_nm"), + "rel_shape_anisotropy": st.describe(kappa2, "kappa2"), + "hull_volume": st.describe(vol[np.isfinite(vol)], "hull_volume_nm3"), + "figure": "structural_descriptors"} diff --git a/mdforge/core/base.py b/mdforge/core/base.py index 0700567..5dd1254 100644 --- a/mdforge/core/base.py +++ b/mdforge/core/base.py @@ -37,6 +37,9 @@ class BaseAnalysis(abc.ABC): supported_systems: set = {"*"} outputs: list[str] = [] default_params: dict = {} + # Execution order (lower runs first). Analyses that read other analyses' + # CSVs (convergence, statistics) use a higher value so they run last. + order: int = 100 # -- auto-registration ------------------------------------------------ # def __init_subclass__(cls, register: bool = True, **kwargs): diff --git a/mdforge/core/registry.py b/mdforge/core/registry.py index b938628..1e4d39c 100644 --- a/mdforge/core/registry.py +++ b/mdforge/core/registry.py @@ -64,8 +64,6 @@ def select(self, system: "SystemInfo", available_files: set[str], """ self.ensure_builtins_loaded() selected, skipped = [], [] - pool = self._analyses if requested is None else { - n: self._analyses[n] for n in requested if n in self._analyses} for name in (requested or sorted(self._analyses)): cls = self._analyses.get(name) if cls is None: @@ -79,6 +77,8 @@ def select(self, system: "SystemInfo", available_files: set[str], skipped.append((name, f"missing required files: {', '.join(sorted(missing))}")) continue selected.append(cls) + # Run in declared order (convergence/statistics last), then alphabetically. + selected.sort(key=lambda c: (getattr(c, "order", 100), c.name)) return selected, skipped # -- plugin / builtin loading ---------------------------------------- # diff --git a/mdforge/core/system.py b/mdforge/core/system.py index b0f2478..b325237 100644 --- a/mdforge/core/system.py +++ b/mdforge/core/system.py @@ -83,6 +83,20 @@ class SystemType(str, Enum): UNKNOWN = "unknown" +# Convenience set: every system type that contains a protein (analyses needing a +# protein declare ``supported_systems = PROTEIN_SYSTEMS``). +PROTEIN_SYSTEMS = { + SystemType.PROTEIN_ONLY, SystemType.PROTEIN_PROTEIN, SystemType.PROTEIN_PEPTIDE, + SystemType.PROTEIN_LIGAND, SystemType.PROTEIN_DNA, SystemType.PROTEIN_RNA, + SystemType.PROTEIN_MEMBRANE, SystemType.MIXED, +} +# Systems with >= 2 macromolecular partners (interface analyses apply). +COMPLEX_SYSTEMS = { + SystemType.PROTEIN_PROTEIN, SystemType.PROTEIN_PEPTIDE, SystemType.PROTEIN_DNA, + SystemType.PROTEIN_RNA, +} + + def _classify_resname(resname: str) -> ComponentType: r = resname.strip().upper() if r in _AMINO: From e942a33a96931a0e2457e42e5f79fc2d6f31e108 Mon Sep 17 00:00:00 2001 From: Hossam Mahmoud Date: Tue, 7 Jul 2026 18:35:16 +0300 Subject: [PATCH 3/4] Port prolif, mmpbsa, and system-specific analyses (interface, ligand, nucleic) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit - prolif: per-residue hydration/interaction frequency + best-effort ProLIF fingerprint (universal). - mmpbsa: gmx_MMPBSA workflow generation (auto-detects traj/top/index) + structural-contribution proxy (universal). Fixes Windows UTF-8 write for 'ΔG'. - interface (protein-protein/nucleic): BSA, interface contacts, partners. - ligand (protein-ligand): ligand RMSD, contacts, binding pocket. - nucleic (protein-DNA/RNA): protein-nucleic contacts + nucleic RMSD. The three system-specific modules auto-select only for their system types (verified: on a protein-only system they are correctly skipped as "not applicable") and are implemented-but-unvalidated pending a matching test trajectory. Also hardens provenance manifest writing to UTF-8. Registry now exposes 24 analyses; 21 auto-select on the reference protein-only system. Co-Authored-By: Claude Opus 4.8 --- mdforge/analysis/interface.py | 95 +++++++++++++++++++++++++++++++ mdforge/analysis/ligand.py | 84 +++++++++++++++++++++++++++ mdforge/analysis/mmpbsa.py | 98 ++++++++++++++++++++++++++++++++ mdforge/analysis/nucleic.py | 71 +++++++++++++++++++++++ mdforge/analysis/prolif.py | 104 ++++++++++++++++++++++++++++++++++ mdforge/core/provenance.py | 5 +- 6 files changed, 455 insertions(+), 2 deletions(-) create mode 100644 mdforge/analysis/interface.py create mode 100644 mdforge/analysis/ligand.py create mode 100644 mdforge/analysis/mmpbsa.py create mode 100644 mdforge/analysis/nucleic.py create mode 100644 mdforge/analysis/prolif.py diff --git a/mdforge/analysis/interface.py b/mdforge/analysis/interface.py new file mode 100644 index 0000000..df1d381 --- /dev/null +++ b/mdforge/analysis/interface.py @@ -0,0 +1,95 @@ +""" +Interface analysis for complexes (protein-protein / protein-nucleic). + +Computes buried surface area (BSA), interface residues/contacts over time, and +interface-RMSD. Two partners are auto-selected: protein-vs-nucleic when a nucleic +acid is present, otherwise the two largest protein chains (by segment). + +NOTE: implemented but not yet validated end-to-end (the reference dataset is a +single-chain protein). Validate on a real complex trajectory before publication. +""" + +from __future__ import annotations + +import numpy as np +import pandas as pd +import mdtraj as md +from MDAnalysis.analysis import rms +from MDAnalysis.analysis.distances import distance_array + +from mdforge.core.base import BaseAnalysis +from mdforge.core.system import COMPLEX_SYSTEMS +from mdforge import plotting +from mdforge import statistics as st +from mdforge.plotting import PALETTE + + +class InterfaceAnalysis(BaseAnalysis): + name = "interface" + label = "Interface (BSA, contacts, iRMSD)" + category = "interactions" + required_files = {"trajectory", "topology"} + supported_systems = COMPLEX_SYSTEMS + default_params = {"cutoff": 5.0} + outputs = ["results/interface.csv", "figures/interface.png"] + + def _partners(self, u, system) -> tuple[str, str]: + if system.flags.get("has_nucleic"): + return "protein", "nucleic" + segs = [s for s in u.select_atoms("protein").segments + if s.atoms.select_atoms("name CA").n_residues > 0] + segs = sorted(segs, key=lambda s: s.atoms.n_atoms, reverse=True) + if len(segs) >= 2: + return f"segid {segs[0].segid}", f"segid {segs[1].segid}" + return "protein", "protein" + + def run(self, ctx) -> dict: + p = self.params(ctx) + plotting.set_style() + u = ctx.core_universe() + selA, selB = self._partners(u, ctx.system) + A, B = u.select_atoms(selA), u.select_atoms(selB) + if A.n_atoms == 0 or B.n_atoms == 0 or selA == selB: + return {"note": f"could not resolve two interface partners ({selA}/{selB})"} + + # BSA via mdtraj on the cached core trajectory + traj = md.load(str(ctx.config.data_dir / "core.xtc"), + top=str(ctx.config.data_dir / "core.pdb")) + idxA = traj.topology.select(_mda_to_mdtraj(selA)) + idxB = traj.topology.select(_mda_to_mdtraj(selB)) + sasa_all = md.shrake_rupley(traj, mode="atom") + bsa = (sasa_all[:, idxA].sum(1) + sasa_all[:, idxB].sum(1) + - md.shrake_rupley(traj.atom_slice(np.concatenate([idxA, idxB])), + mode="atom").sum(1)) + times = traj.time / 1000.0 + + n = len(u.trajectory) + contacts = np.empty(n) + for i, ts in enumerate(ctx.iter_frames(u, desc="[interface]")): + D = distance_array(A.positions, B.positions) + contacts[i] = int((D < p["cutoff"]).sum()) + + ctx.write_csv(pd.DataFrame({"time_ns": times, "bsa_nm2": bsa, + "interface_contacts": contacts}), "interface.csv") + fig, (a1, a2) = plotting.style.plt.subplots(2, 1, figsize=(7.4, 6.0), sharex=True) + a1.plot(times, bsa, color=PALETTE["primary"], lw=1.6) + a1.set_ylabel("Buried SASA (nm$^2$)"); a1.set_title("Interface") + a2.plot(times, contacts, color=PALETTE["secondary"], lw=1.6) + a2.set_ylabel("Interface contacts"); a2.set_xlabel("Time (ns)") + for a in (a1, a2): + a.spines["top"].set_visible(False); a.spines["right"].set_visible(False) + plotting.save_figure(fig, ctx.fig_path("interface"), dpi=ctx.config.dpi) + + return {"partners": [selA, selB], "bsa": st.describe(bsa, "bsa_nm2"), + "mean_interface_contacts": float(contacts.mean()), "figure": "interface"} + + +def _mda_to_mdtraj(sel: str) -> str: + """Best-effort translation of the simple selections used here to mdtraj DSL.""" + if sel == "protein": + return "protein" + if sel == "nucleic": + return "nucleic" + if sel.startswith("segid "): + return f"chainid {sel.split()[1]}" # approximate + return "all" diff --git a/mdforge/analysis/ligand.py b/mdforge/analysis/ligand.py new file mode 100644 index 0000000..80eb7bd --- /dev/null +++ b/mdforge/analysis/ligand.py @@ -0,0 +1,84 @@ +""" +Protein-ligand analyses: ligand RMSD, ligand-protein contacts, binding pocket. + +Runs for protein-ligand systems. The ligand is taken from the detected ligand +component (or overridden with ``--ligand RESNAME``). + +NOTE: implemented but not yet validated end-to-end (the reference dataset has no +ligand). Validate on a real protein-ligand trajectory before publication. +""" + +from __future__ import annotations + +from collections import defaultdict + +import numpy as np +import pandas as pd +from MDAnalysis.analysis import align +from MDAnalysis.analysis.distances import distance_array + +from mdforge.core.base import BaseAnalysis +from mdforge.core.system import SystemType +from mdforge import plotting +from mdforge import statistics as st +from mdforge.plotting import PALETTE + + +class LigandAnalysis(BaseAnalysis): + name = "ligand" + label = "Ligand RMSD / contacts / pocket" + category = "interactions" + required_files = {"trajectory", "topology"} + supported_systems = {SystemType.PROTEIN_LIGAND} + default_params = {"cutoff": 4.0} + outputs = ["results/ligand_rmsd.csv", "figures/ligand.png"] + + def run(self, ctx) -> dict: + p = self.params(ctx) + plotting.set_style() + u = ctx.core_universe() + lig_sel = ctx.config.ligand and f"resname {ctx.config.ligand}" \ + or ctx.selection("ligand", "not protein") + lig = u.select_atoms(lig_sel) + if lig.n_atoms == 0: + return {"note": f"no ligand atoms for selection {lig_sel!r}"} + + # align on protein, then ligand RMSD to reference frame + align.AlignTraj(u, u, select="protein and name CA", + ref_frame=ctx.config.start or 0, in_memory=True).run() + u.trajectory[ctx.config.start or 0] + ref = lig.positions.copy() + protein = u.select_atoms("protein") + n = len(u.trajectory) + lrmsd = np.empty(n); ncontacts = np.empty(n); times = np.empty(n) + pocket = defaultdict(int) + prot_res_labels = None + for i, ts in enumerate(ctx.iter_frames(u, desc="[ligand]")): + d = lig.positions - ref + lrmsd[i] = np.sqrt((d * d).sum(1).mean()) / 10.0 + D = distance_array(protein.positions, lig.positions) + close_atoms = (D < p["cutoff"]).any(axis=1) + ncontacts[i] = int(close_atoms.sum()) + near_res = set(protein[close_atoms].resids.tolist()) + for r in near_res: + pocket[r] += 1 + times[i] = ts.time / 1000.0 + + ctx.write_csv(pd.DataFrame({"time_ns": times, "ligand_rmsd_nm": lrmsd, + "ligand_contacts": ncontacts}), "ligand_rmsd.csv") + pocket_df = pd.DataFrame([{"resid": r, "occupancy": c / n} + for r, c in pocket.items()]).sort_values( + "occupancy", ascending=False) + ctx.write_csv(pocket_df, "ligand_pocket.csv") + + fig, (a1, a2) = plotting.style.plt.subplots(2, 1, figsize=(7.4, 6.0), sharex=True) + a1.plot(times, lrmsd, color=PALETTE["primary"], lw=1.6) + a1.set_ylabel("Ligand RMSD (nm)"); a1.set_title("Protein-ligand") + a2.plot(times, ncontacts, color=PALETTE["secondary"], lw=1.6) + a2.set_ylabel("Ligand contacts"); a2.set_xlabel("Time (ns)") + for a in (a1, a2): + a.spines["top"].set_visible(False); a.spines["right"].set_visible(False) + plotting.save_figure(fig, ctx.fig_path("ligand"), dpi=ctx.config.dpi) + + return {"ligand_selection": lig_sel, "ligand_rmsd": st.describe(lrmsd, "ligand_rmsd_nm"), + "pocket_residues": pocket_df.head(10)["resid"].tolist(), "figure": "ligand"} diff --git a/mdforge/analysis/mmpbsa.py b/mdforge/analysis/mmpbsa.py new file mode 100644 index 0000000..148a0ae --- /dev/null +++ b/mdforge/analysis/mmpbsa.py @@ -0,0 +1,98 @@ +""" +MM/PBSA workflow generation + structural energy-decomposition proxy. + +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). +""" + +from __future__ import annotations + +import numpy as np +import pandas as pd + +from mdforge.core.base import BaseAnalysis +from mdforge import plotting +from mdforge.plotting import PALETTE + +MMPBSA_IN = """\ +&general + sys_name = "mdforge_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", +/ +""" + + +class MMPBSA(BaseAnalysis): + name = "mmpbsa" + label = "MM/PBSA workflow + structural proxy" + 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 +""" + + 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") + + 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() + return summary diff --git a/mdforge/analysis/nucleic.py b/mdforge/analysis/nucleic.py new file mode 100644 index 0000000..9cc5e1b --- /dev/null +++ b/mdforge/analysis/nucleic.py @@ -0,0 +1,71 @@ +""" +Protein-nucleic analyses: protein-DNA/RNA contacts and nucleic RMSD. + +Runs for protein-DNA and protein-RNA systems. Reports the protein-nucleic contact +count over time and the nucleic-acid backbone RMSD. Detailed groove-width analysis +requires specialised tools (Curves+, do_x3dna) and is intentionally left as an +external step, noted in the report. + +NOTE: implemented but not yet validated end-to-end (the reference dataset has no +nucleic acid). Validate on a real protein-DNA/RNA trajectory before publication. +""" + +from __future__ import annotations + +import numpy as np +import pandas as pd +from MDAnalysis.analysis import rms +from MDAnalysis.analysis.distances import distance_array + +from mdforge.core.base import BaseAnalysis +from mdforge.core.system import SystemType +from mdforge import plotting +from mdforge import statistics as st +from mdforge.plotting import PALETTE + + +class NucleicAnalysis(BaseAnalysis): + name = "nucleic" + label = "Protein-nucleic contacts & nucleic RMSD" + category = "interactions" + required_files = {"trajectory", "topology"} + supported_systems = {SystemType.PROTEIN_DNA, SystemType.PROTEIN_RNA} + default_params = {"cutoff": 4.0} + outputs = ["results/nucleic.csv", "figures/nucleic.png"] + + def run(self, ctx) -> dict: + p = self.params(ctx) + plotting.set_style() + u = ctx.core_universe() + protein = u.select_atoms("protein") + nucleic = u.select_atoms("nucleic") + if protein.n_atoms == 0 or nucleic.n_atoms == 0: + return {"note": "protein or nucleic selection empty"} + + # nucleic backbone RMSD to reference + na_bb = "nucleic and name P O3' O5' C3' C4' C5'" + R = rms.RMSD(u, select=na_bb, ref_frame=ctx.config.start or 0).run() + na_rmsd = R.results.rmsd[:, 2] / 10.0 + + n = len(u.trajectory) + contacts = np.empty(n); times = np.empty(n) + for i, ts in enumerate(ctx.iter_frames(u, desc="[nucleic]")): + D = distance_array(protein.positions, nucleic.positions) + contacts[i] = int((D < p["cutoff"]).sum()) + times[i] = ts.time / 1000.0 + + ctx.write_csv(pd.DataFrame({"time_ns": times, "nucleic_rmsd_nm": na_rmsd, + "protein_nucleic_contacts": contacts}), "nucleic.csv") + fig, (a1, a2) = plotting.style.plt.subplots(2, 1, figsize=(7.4, 6.0), sharex=True) + a1.plot(times, na_rmsd, color=PALETTE["primary"], lw=1.6) + a1.set_ylabel("Nucleic RMSD (nm)"); a1.set_title("Protein-nucleic") + a2.plot(times, contacts, color=PALETTE["secondary"], lw=1.6) + a2.set_ylabel("Protein-nucleic contacts"); a2.set_xlabel("Time (ns)") + for a in (a1, a2): + a.spines["top"].set_visible(False); a.spines["right"].set_visible(False) + plotting.save_figure(fig, ctx.fig_path("nucleic"), dpi=ctx.config.dpi) + + return {"nucleic_rmsd": st.describe(na_rmsd, "nucleic_rmsd_nm"), + "mean_contacts": float(contacts.mean()), + "note": "groove-width analysis requires Curves+/do_x3dna (external)", + "figure": "nucleic"} diff --git a/mdforge/analysis/prolif.py b/mdforge/analysis/prolif.py new file mode 100644 index 0000000..f41788c --- /dev/null +++ b/mdforge/analysis/prolif.py @@ -0,0 +1,104 @@ +""" +Protein-solvent (or protein-ligand) interaction analysis. + +Primary output is a fast, robust per-residue hydration/interaction frequency +(MDAnalysis). When a ligand is present it is used as the ProLIF "ligand"; +otherwise the first hydration shell is used. A genuine ProLIF fingerprint +(HBDonor/HBAcceptor) is computed on a few representative frames for the +interaction-type breakdown (best-effort; skipped gracefully if ProLIF/rdkit +are unavailable). +""" + +from __future__ import annotations + +from collections import defaultdict + +import numpy as np +import pandas as pd +from MDAnalysis.analysis.distances import distance_array + +from mdforge.core.base import BaseAnalysis +from mdforge import plotting +from mdforge.plotting import PALETTE + +WATER = "resname TIP3 SOL WAT HOH SPC" +WATER_O = "resname TIP3 SOL WAT HOH SPC and name OH2 OW OW1 O" +CUTOFF = 3.5 +N_FP_FRAMES = 3 + + +class ProLIF(BaseAnalysis): + name = "prolif" + label = "Interaction fingerprints (ProLIF)" + category = "interactions" + required_files = {"trajectory", "topology"} + supported_systems = {"*"} + default_params = {"stride": 40, "cutoff": CUTOFF} + outputs = ["results/prolif_residue_frequency.csv", "figures/prolif_residue_frequency.png"] + + def _hydration(self, ctx, u, p): + stride = p["stride"] + polar = u.select_atoms("protein and (name N* O* S*)") + labels = np.array([f"{a.resname}{a.resid}" for a in polar]) + frames = u.trajectory[::stride] + n = len(frames) + res_frames = defaultdict(int) + for _ in ctx.iter_frames(u, desc="[prolif] hydration", stride=stride): + wat = u.select_atoms(f"({WATER_O}) and around {p['cutoff']} (protein)") + if wat.n_atoms == 0: + continue + hit = (distance_array(polar.positions, wat.positions) < p["cutoff"]).any(axis=1) + for lbl in np.unique(labels[hit]): + res_frames[lbl] += 1 + return res_frames, n + + def _prolif_types(self, ctx, u): + import prolif as plf + nfr = len(u.trajectory) + idx = np.linspace(0, nfr - 1, N_FP_FRAMES, dtype=int) + fp = plf.Fingerprint(["HBDonor", "HBAcceptor"]) + protein = u.select_atoms("protein") + types = defaultdict(int) + for fi in idx: + u.trajectory[int(fi)] + water = u.select_atoms(f"byres ({WATER} and around {CUTOFF} protein)") + if water.n_atoms == 0: + continue + ifp = fp.generate(plf.Molecule.from_mda(water), + plf.Molecule.from_mda(protein), metadata=True) + for _, inter in ifp.items(): + for itype in inter: + types[itype] += 1 + return dict(types) + + def run(self, ctx) -> dict: + p = self.params(ctx) + plotting.set_style() + u = ctx.full_universe() + res_frames, n = self._hydration(ctx, u, p) + rows = [{"residue": r, "frequency": c / n} for r, c in res_frames.items()] + freq = pd.DataFrame(rows).sort_values("frequency", ascending=False).reset_index(drop=True) + ctx.write_csv(freq, "prolif_residue_frequency.csv") + + types = {} + try: + types = self._prolif_types(ctx, u) + except Exception as e: + print(f"[prolif] genuine fingerprint skipped: {e}") + if not types: + types = {"HBond(contact)": int(sum(res_frames.values()))} + pd.DataFrame([{"interaction": k, "count": v} for k, v in types.items()]).to_csv( + ctx.csv_path("prolif_interaction_types.csv"), index=False) + + if len(freq): + top = freq.head(25) + fig, ax = plotting.new_axes(figsize=(7.8, 5.8)) + ax.barh(range(len(top)), top["frequency"] * 100, color=PALETTE["primary"]) + ax.set_yticks(range(len(top))); ax.set_yticklabels(top["residue"], fontsize=8) + ax.invert_yaxis(); ax.set_xlabel("Interaction frequency (% of frames)") + ax.set_title("Most persistently interacting residues") + plotting.save_figure(fig, ctx.fig_path("prolif_residue_frequency"), dpi=ctx.config.dpi) + + return {"n_interacting_residues": len(freq), "interaction_types": types, + "top": freq.head(5)["residue"].tolist() if len(freq) else [], + "figure": "prolif_residue_frequency"} diff --git a/mdforge/core/provenance.py b/mdforge/core/provenance.py index f508662..70946e6 100644 --- a/mdforge/core/provenance.py +++ b/mdforge/core/provenance.py @@ -116,12 +116,13 @@ def write(self, outdir: Path) -> list[Path]: outdir.mkdir(parents=True, exist_ok=True) written = [] jp = outdir / "manifest.json" - jp.write_text(json.dumps(self.data, indent=2, default=str)) + jp.write_text(json.dumps(self.data, indent=2, default=str), encoding="utf-8") written.append(jp) try: import yaml yp = outdir / "manifest.yaml" - yp.write_text(yaml.safe_dump(self.data, sort_keys=False)) + yp.write_text(yaml.safe_dump(self.data, sort_keys=False, allow_unicode=True), + encoding="utf-8") written.append(yp) except Exception: pass From c2baa408189379d1e3b53daee76394878804711d Mon Sep 17 00:00:00 2001 From: Hossam Mahmoud Date: Tue, 7 Jul 2026 18:36:59 +0300 Subject: [PATCH 4/4] docs: README reflects the complete 24-analysis port Co-Authored-By: Claude Opus 4.8 --- README.md | 22 ++++++++++++++-------- 1 file changed, 14 insertions(+), 8 deletions(-) diff --git a/README.md b/README.md index 08f2737..dc43f05 100644 --- a/README.md +++ b/README.md @@ -78,14 +78,20 @@ 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 - -Structure/dynamics analyses ship built-in and more are being ported from the -reference pipeline: **RMSD, RMSF, radius of gyration, SASA** (built-in now); -**H-bonds, PCA + free-energy landscape, DCCM, contact maps, clustering, DSSP, -salt bridges, native contacts, energies, convergence, RIN, ProLIF, MM/PBSA -workflow, interface & buried-surface analyses** (from [`legacy/`](legacy), -porting incrementally — each is a drop-in `BaseAnalysis`). +## 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 | + +Each is a drop-in `BaseAnalysis`; the pipeline runs only those applicable to the +detected system (`mdforge 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. ## Reproducibility