From 8044ecd00da4bcdd6e1a1d382c4886b1c8980c18 Mon Sep 17 00:00:00 2001 From: Derrick K Date: Wed, 11 Feb 2026 13:22:40 -0800 Subject: [PATCH 1/8] Fix dihedral residue option aliases and add coverage - Accept residue/residue_selection aliases for phi/psi/omega and combined dihedrals without warnings - Preserve residue filtering in combined dihedrals and align warning handling with other analyses - Add tests for residue alias parsing and filtered outputs across dihedral analyses --- src/fastmdanalysis/analysis/dihedrals.py | 28 +++++++++++++----- tests/test_dihedrals_residue_options.py | 36 ++++++++++++++++++++++++ 2 files changed, 57 insertions(+), 7 deletions(-) create mode 100644 tests/test_dihedrals_residue_options.py diff --git a/src/fastmdanalysis/analysis/dihedrals.py b/src/fastmdanalysis/analysis/dihedrals.py index b4b9020..a0cb5f7 100644 --- a/src/fastmdanalysis/analysis/dihedrals.py +++ b/src/fastmdanalysis/analysis/dihedrals.py @@ -36,7 +36,8 @@ class PhiAnalysis(BaseAnalysis): """ _ALIASES = { - "residues": "residue_selection", + "residue": "residues", + "residue_selection": "residues", } def __init__( @@ -241,7 +242,8 @@ class PsiAnalysis(BaseAnalysis): # Similar structure to PhiAnalysis, but for psi angles _ALIASES = { - "residues": "residue_selection", + "residue": "residues", + "residue_selection": "residues", } def __init__( @@ -253,13 +255,15 @@ def __init__( **kwargs ): logger.info("Initializing Psi analysis") + warn_unknown = kwargs.pop("_warn_unknown", False) + analysis_opts = {"residues": residues, "units": units, "strict": strict} analysis_opts.update(kwargs) forwarder = OptionsForwarder(aliases=self._ALIASES, strict=strict) resolved = forwarder.apply_aliases(analysis_opts) resolved = forwarder.filter_known( - resolved, {"residues", "units", "strict", "output"}, context="psi" + resolved, {"residues", "units", "strict", "output"}, context="psi", warn=warn_unknown ) residues = resolved.get("residues", None) @@ -366,7 +370,8 @@ class OmegaAnalysis(BaseAnalysis): """ _ALIASES = { - "residues": "residue_selection", + "residue": "residues", + "residue_selection": "residues", } def __init__( @@ -378,13 +383,15 @@ def __init__( **kwargs ): logger.info("Initializing Omega analysis") + warn_unknown = kwargs.pop("_warn_unknown", False) + analysis_opts = {"residues": residues, "units": units, "strict": strict} analysis_opts.update(kwargs) forwarder = OptionsForwarder(aliases=self._ALIASES, strict=strict) resolved = forwarder.apply_aliases(analysis_opts) resolved = forwarder.filter_known( - resolved, {"residues", "units", "strict", "output"}, context="omega" + resolved, {"residues", "units", "strict", "output"}, context="omega", warn=warn_unknown ) residues = resolved.get("residues", None) @@ -489,6 +496,11 @@ class DihedralsAnalysis(BaseAnalysis): Combined dihedral analysis for phi, psi, omega with Ramachandran plotting. """ + _ALIASES = { + "residue": "residues", + "residue_selection": "residues", + } + def __init__( self, trajectory: md.Trajectory, @@ -499,13 +511,15 @@ def __init__( **kwargs ): logger.info("Initializing Dihedrals analysis") + warn_unknown = kwargs.pop("_warn_unknown", False) + analysis_opts = {"types": types, "residues": residues, "units": units, "strict": strict} analysis_opts.update(kwargs) - forwarder = OptionsForwarder(strict=strict) + forwarder = OptionsForwarder(aliases=self._ALIASES, strict=strict) resolved = forwarder.apply_aliases(analysis_opts) resolved = forwarder.filter_known( - resolved, {"types", "residues", "units", "strict", "output"}, context="dihedrals" + resolved, {"types", "residues", "units", "strict", "output"}, context="dihedrals", warn=warn_unknown ) types = resolved.get("types", ["phi", "psi", "omega"]) diff --git a/tests/test_dihedrals_residue_options.py b/tests/test_dihedrals_residue_options.py new file mode 100644 index 0000000..36728cf --- /dev/null +++ b/tests/test_dihedrals_residue_options.py @@ -0,0 +1,36 @@ +import warnings +import numpy as np + + +def _unknown_warnings(messages): + return [m for m in messages if "Unknown options" in m or "Unsupported options" in m] + + +def test_phi_residue_alias_no_unknown_warning(fastmda): + with warnings.catch_warnings(record=True) as w: + warnings.simplefilter("always") + analysis = fastmda.phi(residue=0) + + messages = [str(wi.message) for wi in w] + assert not _unknown_warnings(messages) + + assert "phi_avg_filtered" in analysis.results + filtered = analysis.results["phi_avg_filtered"] + assert filtered.shape[0] == 1 + np.testing.assert_allclose(filtered[0], analysis.data[0], atol=1e-6) + + +def test_dihedrals_residue_selection_alias_propagates(fastmda): + with warnings.catch_warnings(record=True) as w: + warnings.simplefilter("always") + analysis = fastmda.dihedrals(residue_selection=[0, 1]) + + messages = [str(wi.message) for wi in w] + assert not _unknown_warnings(messages) + + for key in ("phi_avg_filtered", "psi_avg_filtered", "omega_avg_filtered"): + assert key in analysis.results + filtered = analysis.results[key] + assert filtered.shape[0] == 2 + base_key = key.replace("_filtered", "") + np.testing.assert_allclose(filtered[:, 0], analysis.results[base_key][:2, 0], atol=1e-6) From dadb281d3703a49de7628837b5e6a2a1cf78d36d Mon Sep 17 00:00:00 2001 From: Derrick K Date: Wed, 11 Feb 2026 17:58:11 -0800 Subject: [PATCH 2/8] Fix dihedrals --residues filtering and add tests --- src/fastmdanalysis/analysis/dihedrals.py | 98 ++++++++++++++++-------- tests/test_dihedrals_residue_options.py | 14 ++++ 2 files changed, 82 insertions(+), 30 deletions(-) diff --git a/src/fastmdanalysis/analysis/dihedrals.py b/src/fastmdanalysis/analysis/dihedrals.py index a0cb5f7..eca2762 100644 --- a/src/fastmdanalysis/analysis/dihedrals.py +++ b/src/fastmdanalysis/analysis/dihedrals.py @@ -94,6 +94,9 @@ def __init__( self.units: str = units self.strict = strict + # Residue indices corresponding to the rows in self.data + self.residue_indices: Optional[np.ndarray] = None + # Populated during run() self.data: Optional[np.ndarray] = None self.results: Dict[str, np.ndarray] = {} @@ -115,6 +118,15 @@ def run(self) -> Dict[str, np.ndarray]: if angles.size == 0: raise AnalysisError("No phi angles found in trajectory (no protein?)") + # Restrict computation to selected residues (0-based dihedral/residue index) + full_res_indices = np.arange(angles.shape[1]) + if self.residues is not None: + res_list = [self.residues] if isinstance(self.residues, int) else list(self.residues) + angles = angles[:, res_list] + self.residue_indices = np.asarray(res_list, dtype=int) + else: + self.residue_indices = full_res_indices + # Circular mean per residue n_residues = angles.shape[1] avg_angles = np.zeros(n_residues) @@ -132,16 +144,12 @@ def run(self) -> Dict[str, np.ndarray]: avg_angles = np.degrees(avg_angles) self.data = avg_angles.reshape(-1, 1) - self.results = {"phi_avg": self.data} - - # Filter by residues if specified - if self.residues is not None: - if isinstance(self.residues, int): - res_list = [self.residues] - else: - res_list = list(self.residues) - filtered_data = self.data[res_list] - self.results["phi_avg_filtered"] = filtered_data.reshape(-1, 1) + # If residues were provided, self.data is already filtered. + self.results = { + "phi_avg": self.data, + "phi_avg_filtered": self.data, + "phi_residues": self.residue_indices, + } # Save data self._save_data(self.data, "phi_avg", header=f"phi_avg_{self.units}") @@ -197,8 +205,12 @@ def plot( raise AnalysisError("No phi data available to plot.") y = np.asarray(data, dtype=float).flatten() - n = len(y) - x = np.arange(n) + + # X-axis should reflect residue indices of the computed data (not 0..N-1) + if self.residue_indices is not None and len(self.residue_indices) == len(y): + x = self.residue_indices.astype(int) + else: + x = np.arange(len(y)) # Filter residues if residues is not None: @@ -278,6 +290,8 @@ def __init__( self.data = None self.results = {} + self.residue_indices: Optional[np.ndarray] = None + def run(self) -> Dict[str, np.ndarray]: logger.info("Starting Psi analysis") try: @@ -285,6 +299,14 @@ def run(self) -> Dict[str, np.ndarray]: if angles.size == 0: raise AnalysisError("No psi angles found in trajectory") + full_res_indices = np.arange(angles.shape[1]) + if self.residues is not None: + res_list = [self.residues] if isinstance(self.residues, int) else list(self.residues) + angles = angles[:, res_list] + self.residue_indices = np.asarray(res_list, dtype=int) + else: + self.residue_indices = full_res_indices + n_residues = angles.shape[1] avg_angles = np.zeros(n_residues) for i in range(n_residues): @@ -299,12 +321,11 @@ def run(self) -> Dict[str, np.ndarray]: avg_angles = np.degrees(avg_angles) self.data = avg_angles.reshape(-1, 1) - self.results = {"psi_avg": self.data} - - if self.residues is not None: - res_list = [self.residues] if isinstance(self.residues, int) else list(self.residues) - filtered_data = self.data[res_list] - self.results["psi_avg_filtered"] = filtered_data.reshape(-1, 1) + self.results = { + "psi_avg": self.data, + "psi_avg_filtered": self.data, + "psi_residues": self.residue_indices, + } self._save_data(self.data, "psi_avg", header=f"psi_avg_{self.units}") plot_path = self.plot() @@ -327,8 +348,10 @@ def plot(self, **kwargs) -> str: raise AnalysisError("No psi data available to plot.") y = np.asarray(kwargs["data"], dtype=float).flatten() - n = len(y) - x = np.arange(n) + if self.residue_indices is not None and len(self.residue_indices) == len(y): + x = self.residue_indices.astype(int) + else: + x = np.arange(len(y)) # Filter residues residues = kwargs.get("residues") @@ -406,6 +429,8 @@ def __init__( self.data = None self.results = {} + self.residue_indices: Optional[np.ndarray] = None + def run(self) -> Dict[str, np.ndarray]: logger.info("Starting Omega analysis") try: @@ -413,6 +438,14 @@ def run(self) -> Dict[str, np.ndarray]: if angles.size == 0: raise AnalysisError("No omega angles found in trajectory") + full_res_indices = np.arange(angles.shape[1]) + if self.residues is not None: + res_list = [self.residues] if isinstance(self.residues, int) else list(self.residues) + angles = angles[:, res_list] + self.residue_indices = np.asarray(res_list, dtype=int) + else: + self.residue_indices = full_res_indices + n_residues = angles.shape[1] avg_angles = np.zeros(n_residues) for i in range(n_residues): @@ -427,12 +460,11 @@ def run(self) -> Dict[str, np.ndarray]: avg_angles = np.degrees(avg_angles) self.data = avg_angles.reshape(-1, 1) - self.results = {"omega_avg": self.data} - - if self.residues is not None: - res_list = [self.residues] if isinstance(self.residues, int) else list(self.residues) - filtered_data = self.data[res_list] - self.results["omega_avg_filtered"] = filtered_data.reshape(-1, 1) + self.results = { + "omega_avg": self.data, + "omega_avg_filtered": self.data, + "omega_residues": self.residue_indices, + } self._save_data(self.data, "omega_avg", header=f"omega_avg_{self.units}") plot_path = self.plot() @@ -454,8 +486,10 @@ def plot(self, **kwargs) -> str: raise AnalysisError("No omega data available to plot.") y = np.asarray(kwargs["data"], dtype=float).flatten() - n = len(y) - x = np.arange(n) + if self.residue_indices is not None and len(self.residue_indices) == len(y): + x = self.residue_indices.astype(int) + else: + x = np.arange(len(y)) # Filter residues residues = kwargs.get("residues") @@ -602,8 +636,12 @@ def plot_ramachandran( x = phi_data y = psi_data - n = len(x) - res_indices = np.arange(n) + # Use residue indices when the per-angle analyses were residue-filtered + res_indices = self.results.get("phi_residues") + if res_indices is None: + res_indices = np.arange(len(x)) + else: + res_indices = np.asarray(res_indices, dtype=int) # Filter residues if residues is not None: diff --git a/tests/test_dihedrals_residue_options.py b/tests/test_dihedrals_residue_options.py index 36728cf..41a6e7e 100644 --- a/tests/test_dihedrals_residue_options.py +++ b/tests/test_dihedrals_residue_options.py @@ -14,6 +14,11 @@ def test_phi_residue_alias_no_unknown_warning(fastmda): messages = [str(wi.message) for wi in w] assert not _unknown_warnings(messages) + assert analysis.data.shape[0] == 1 + + assert "phi_residues" in analysis.results + assert list(np.asarray(analysis.results["phi_residues"]).astype(int)) == [0] + assert "phi_avg_filtered" in analysis.results filtered = analysis.results["phi_avg_filtered"] assert filtered.shape[0] == 1 @@ -28,6 +33,15 @@ def test_dihedrals_residue_selection_alias_propagates(fastmda): messages = [str(wi.message) for wi in w] assert not _unknown_warnings(messages) + # Ensure the combined analysis truly computed only 2 residues worth of data + assert analysis.results["phi_avg"].shape[0] == 2 + assert analysis.results["psi_avg"].shape[0] == 2 + assert analysis.results["omega_avg"].shape[0] == 2 + + assert list(np.asarray(analysis.results["phi_residues"]).astype(int)) == [0, 1] + assert list(np.asarray(analysis.results["psi_residues"]).astype(int)) == [0, 1] + assert list(np.asarray(analysis.results["omega_residues"]).astype(int)) == [0, 1] + for key in ("phi_avg_filtered", "psi_avg_filtered", "omega_avg_filtered"): assert key in analysis.results filtered = analysis.results[key] From 8df9395f42b12dc3652fb6b2fc6456110ebc407f Mon Sep 17 00:00:00 2001 From: Derrick K Date: Wed, 11 Feb 2026 18:53:13 -0800 Subject: [PATCH 3/8] Update dihedrals.py --- src/fastmdanalysis/analysis/dihedrals.py | 148 +++++++++++++++++++++-- 1 file changed, 135 insertions(+), 13 deletions(-) diff --git a/src/fastmdanalysis/analysis/dihedrals.py b/src/fastmdanalysis/analysis/dihedrals.py index eca2762..e533f56 100644 --- a/src/fastmdanalysis/analysis/dihedrals.py +++ b/src/fastmdanalysis/analysis/dihedrals.py @@ -13,11 +13,12 @@ from __future__ import annotations from typing import Dict, Optional, Sequence, Union, List +from pathlib import Path import logging import numpy as np import mdtraj as md -from scipy.stats import circmean +from scipy.stats import circmean, circstd import matplotlib matplotlib.use("Agg") @@ -127,21 +128,25 @@ def run(self) -> Dict[str, np.ndarray]: else: self.residue_indices = full_res_indices - # Circular mean per residue + # Circular mean and std per residue n_residues = angles.shape[1] avg_angles = np.zeros(n_residues) + std_angles = np.zeros(n_residues) for i in range(n_residues): res_angles = angles[:, i] # Remove NaN values (missing angles) valid_angles = res_angles[~np.isnan(res_angles)] if len(valid_angles) > 0: avg_angles[i] = circmean(valid_angles, high=np.pi, low=-np.pi) + std_angles[i] = circstd(valid_angles, high=np.pi, low=-np.pi) else: avg_angles[i] = np.nan + std_angles[i] = np.nan # Convert units if self.units == "degrees": avg_angles = np.degrees(avg_angles) + std_angles = np.degrees(std_angles) self.data = avg_angles.reshape(-1, 1) # If residues were provided, self.data is already filtered. @@ -149,6 +154,7 @@ def run(self) -> Dict[str, np.ndarray]: "phi_avg": self.data, "phi_avg_filtered": self.data, "phi_residues": self.residue_indices, + "phi_std": std_angles.reshape(-1, 1), } # Save data @@ -172,6 +178,7 @@ def plot( self, data: Optional[Union[Sequence[float], np.ndarray]] = None, *, + std_data: Optional[Union[Sequence[float], np.ndarray]] = None, residues: Optional[Union[int, Sequence[int]]] = None, highlight_residues: Optional[Union[int, Sequence[int]]] = None, figsize=(12, 6), @@ -201,10 +208,15 @@ def plot( """ if data is None: data = self.data + if std_data is None: + std_data = self.results.get("phi_std") if data is None: raise AnalysisError("No phi data available to plot.") y = np.asarray(data, dtype=float).flatten() + yerr = None + if std_data is not None: + yerr = np.asarray(std_data, dtype=float).flatten() # X-axis should reflect residue indices of the computed data (not 0..N-1) if self.residue_indices is not None and len(self.residue_indices) == len(y): @@ -222,7 +234,16 @@ def plot( # Plot fig, ax = plt.subplots(figsize=figsize) - ax.bar(x, y, color=color, **kwargs) + ax.errorbar( + x, + y, + yerr=yerr, + fmt="o", + color=color, + ecolor=kwargs.pop("ecolor", "gray"), + capsize=kwargs.pop("capsize", 3), + **kwargs, + ) ax.set_title(title) ax.set_xlabel(xlabel) @@ -237,7 +258,7 @@ def plot( for res in highlight_residues: if res in x: idx = np.where(x == res)[0][0] - ax.bar([x[idx]], [y[idx]], color="red", alpha=0.7) + ax.scatter([x[idx]], [y[idx]], color="red", alpha=0.7, zorder=3) apply_slide_style(ax, x_values=x, y_values=y) fig.tight_layout() @@ -309,22 +330,27 @@ def run(self) -> Dict[str, np.ndarray]: n_residues = angles.shape[1] avg_angles = np.zeros(n_residues) + std_angles = np.zeros(n_residues) for i in range(n_residues): res_angles = angles[:, i] valid_angles = res_angles[~np.isnan(res_angles)] if len(valid_angles) > 0: avg_angles[i] = circmean(valid_angles, high=np.pi, low=-np.pi) + std_angles[i] = circstd(valid_angles, high=np.pi, low=-np.pi) else: avg_angles[i] = np.nan + std_angles[i] = np.nan if self.units == "degrees": avg_angles = np.degrees(avg_angles) + std_angles = np.degrees(std_angles) self.data = avg_angles.reshape(-1, 1) self.results = { "psi_avg": self.data, "psi_avg_filtered": self.data, "psi_residues": self.residue_indices, + "psi_std": std_angles.reshape(-1, 1), } self._save_data(self.data, "psi_avg", header=f"psi_avg_{self.units}") @@ -344,10 +370,15 @@ def plot(self, **kwargs) -> str: if "data" not in kwargs: kwargs["data"] = self.data + if "std_data" not in kwargs: + kwargs["std_data"] = self.results.get("psi_std") if kwargs["data"] is None: raise AnalysisError("No psi data available to plot.") y = np.asarray(kwargs["data"], dtype=float).flatten() + yerr = None + if kwargs.get("std_data") is not None: + yerr = np.asarray(kwargs["std_data"], dtype=float).flatten() if self.residue_indices is not None and len(self.residue_indices) == len(y): x = self.residue_indices.astype(int) else: @@ -363,7 +394,33 @@ def plot(self, **kwargs) -> str: y = y[mask] fig, ax = plt.subplots(figsize=kwargs.get("figsize", (12, 6))) - ax.bar(x, y, color=kwargs.get("color", "green"), **{k: v for k, v in kwargs.items() if k not in ["data", "residues", "figsize", "title", "xlabel", "ylabel", "color", "highlight_residues"]}) + ax.errorbar( + x, + y, + yerr=yerr, + fmt="o", + color=kwargs.get("color", "green"), + ecolor=kwargs.get("ecolor", "gray"), + capsize=kwargs.get("capsize", 3), + **{ + k: v + for k, v in kwargs.items() + if k + not in [ + "data", + "std_data", + "residues", + "figsize", + "title", + "xlabel", + "ylabel", + "color", + "highlight_residues", + "ecolor", + "capsize", + ] + }, + ) ax.set_title(kwargs.get("title", "Average Psi Angles per Residue")) ax.set_xlabel(kwargs.get("xlabel", "Residue Index")) @@ -377,7 +434,7 @@ def plot(self, **kwargs) -> str: for res in highlight_residues: if res in x: idx = np.where(x == res)[0][0] - ax.bar([x[idx]], [y[idx]], color="red", alpha=0.7) + ax.scatter([x[idx]], [y[idx]], color="red", alpha=0.7, zorder=3) apply_slide_style(ax, x_values=x, y_values=y) fig.tight_layout() @@ -448,22 +505,27 @@ def run(self) -> Dict[str, np.ndarray]: n_residues = angles.shape[1] avg_angles = np.zeros(n_residues) + std_angles = np.zeros(n_residues) for i in range(n_residues): res_angles = angles[:, i] valid_angles = res_angles[~np.isnan(res_angles)] if len(valid_angles) > 0: avg_angles[i] = circmean(valid_angles, high=np.pi, low=-np.pi) + std_angles[i] = circstd(valid_angles, high=np.pi, low=-np.pi) else: avg_angles[i] = np.nan + std_angles[i] = np.nan if self.units == "degrees": avg_angles = np.degrees(avg_angles) + std_angles = np.degrees(std_angles) self.data = avg_angles.reshape(-1, 1) self.results = { "omega_avg": self.data, "omega_avg_filtered": self.data, "omega_residues": self.residue_indices, + "omega_std": std_angles.reshape(-1, 1), } self._save_data(self.data, "omega_avg", header=f"omega_avg_{self.units}") @@ -482,10 +544,15 @@ def plot(self, **kwargs) -> str: if "data" not in kwargs: kwargs["data"] = self.data + if "std_data" not in kwargs: + kwargs["std_data"] = self.results.get("omega_std") if kwargs["data"] is None: raise AnalysisError("No omega data available to plot.") y = np.asarray(kwargs["data"], dtype=float).flatten() + yerr = None + if kwargs.get("std_data") is not None: + yerr = np.asarray(kwargs["std_data"], dtype=float).flatten() if self.residue_indices is not None and len(self.residue_indices) == len(y): x = self.residue_indices.astype(int) else: @@ -501,7 +568,33 @@ def plot(self, **kwargs) -> str: y = y[mask] fig, ax = plt.subplots(figsize=kwargs.get("figsize", (12, 6))) - ax.bar(x, y, color=kwargs.get("color", "orange"), **{k: v for k, v in kwargs.items() if k not in ["data", "residues", "figsize", "title", "xlabel", "ylabel", "color", "highlight_residues"]}) + ax.errorbar( + x, + y, + yerr=yerr, + fmt="o", + color=kwargs.get("color", "orange"), + ecolor=kwargs.get("ecolor", "gray"), + capsize=kwargs.get("capsize", 3), + **{ + k: v + for k, v in kwargs.items() + if k + not in [ + "data", + "std_data", + "residues", + "figsize", + "title", + "xlabel", + "ylabel", + "color", + "highlight_residues", + "ecolor", + "capsize", + ] + }, + ) ax.set_title(kwargs.get("title", "Average Omega Angles per Residue")) ax.set_xlabel(kwargs.get("xlabel", "Residue Index")) @@ -515,7 +608,7 @@ def plot(self, **kwargs) -> str: for res in highlight_residues: if res in x: idx = np.where(x == res)[0][0] - ax.bar([x[idx]], [y[idx]], color="red", alpha=0.7) + ax.scatter([x[idx]], [y[idx]], color="red", alpha=0.7, zorder=3) apply_slide_style(ax, x_values=x, y_values=y) fig.tight_layout() @@ -618,6 +711,8 @@ def plot_ramachandran( self, phi_data: Optional[np.ndarray] = None, psi_data: Optional[np.ndarray] = None, + phi_std: Optional[np.ndarray] = None, + psi_std: Optional[np.ndarray] = None, residues: Optional[Union[int, Sequence[int]]] = None, figsize=(8, 8), title: str = "Ramachandran Plot (Average Angles)", @@ -630,12 +725,17 @@ def plot_ramachandran( phi_data = self.results["phi_avg"].flatten() if psi_data is None and "psi_avg" in self.results: psi_data = self.results["psi_avg"].flatten() + if phi_std is None and "phi_std" in self.results: + phi_std = self.results["phi_std"].flatten() + if psi_std is None and "psi_std" in self.results: + psi_std = self.results["psi_std"].flatten() if phi_data is None or psi_data is None: raise AnalysisError("Phi and psi data required for Ramachandran plot") x = phi_data y = psi_data + # Use residue indices when the per-angle analyses were residue-filtered res_indices = self.results.get("phi_residues") if res_indices is None: @@ -643,7 +743,8 @@ def plot_ramachandran( else: res_indices = np.asarray(res_indices, dtype=int) - # Filter residues + if residues is None: + residues = self.residues if residues is not None: if isinstance(residues, int): residues = [residues] @@ -651,17 +752,38 @@ def plot_ramachandran( x = x[mask] y = y[mask] res_indices = res_indices[mask] + if phi_std is not None: + phi_std = phi_std[mask] + if psi_std is not None: + psi_std = psi_std[mask] fig, ax = plt.subplots(figsize=figsize) - scatter = ax.scatter(x, y, c=res_indices, cmap="viridis", **kwargs) + cmap = plt.get_cmap("viridis") + norm = plt.Normalize(vmin=res_indices.min(), vmax=res_indices.max()) if len(res_indices) else None + + for i, res in enumerate(res_indices): + color = cmap(norm(res)) if norm is not None else "blue" + ax.errorbar( + x[i], + y[i], + xerr=None if phi_std is None else phi_std[i], + yerr=None if psi_std is None else psi_std[i], + fmt="o", + color=color, + ecolor=color, + capsize=3, + **{k: v for k, v in kwargs.items() if k not in ["capsize"]}, + ) + ax.set_title(title) ax.set_xlabel(f"Phi ({self.units})") ax.set_ylabel(f"Psi ({self.units})") ax.grid(True, alpha=0.3) - # Add colorbar - cbar = plt.colorbar(scatter, ax=ax) - cbar.set_label("Residue Index") + if len(res_indices): + mappable = matplotlib.cm.ScalarMappable(norm=norm, cmap=cmap) + cbar = plt.colorbar(mappable, ax=ax) + cbar.set_label("Residue Index") fig.tight_layout() outpath = self._save_plot(fig, "ramachandran") From f0096879343321c2af92501990b7f0b705eec028 Mon Sep 17 00:00:00 2001 From: Derrick K Date: Wed, 11 Feb 2026 19:12:37 -0800 Subject: [PATCH 4/8] Update dihedrals.py --- src/fastmdanalysis/analysis/dihedrals.py | 66 +++++++++++++++--------- 1 file changed, 43 insertions(+), 23 deletions(-) diff --git a/src/fastmdanalysis/analysis/dihedrals.py b/src/fastmdanalysis/analysis/dihedrals.py index d15b7bf..fa40b02 100644 --- a/src/fastmdanalysis/analysis/dihedrals.py +++ b/src/fastmdanalysis/analysis/dihedrals.py @@ -128,11 +128,7 @@ def run(self) -> Dict[str, np.ndarray]: else: self.residue_indices = full_res_indices -<<<<<<< HEAD # Circular mean and std per residue -======= - # Circular mean per residue ->>>>>>> upstream/main n_residues = angles.shape[1] avg_angles = np.zeros(n_residues) std_angles = np.zeros(n_residues) @@ -158,10 +154,7 @@ def run(self) -> Dict[str, np.ndarray]: "phi_avg": self.data, "phi_avg_filtered": self.data, "phi_residues": self.residue_indices, -<<<<<<< HEAD "phi_std": std_angles.reshape(-1, 1), -======= ->>>>>>> upstream/main } # Save data @@ -221,12 +214,9 @@ def plot( raise AnalysisError("No phi data available to plot.") y = np.asarray(data, dtype=float).flatten() -<<<<<<< HEAD yerr = None if std_data is not None: yerr = np.asarray(std_data, dtype=float).flatten() -======= ->>>>>>> upstream/main # X-axis should reflect residue indices of the computed data (not 0..N-1) if self.residue_indices is not None and len(self.residue_indices) == len(y): @@ -360,10 +350,7 @@ def run(self) -> Dict[str, np.ndarray]: "psi_avg": self.data, "psi_avg_filtered": self.data, "psi_residues": self.residue_indices, -<<<<<<< HEAD "psi_std": std_angles.reshape(-1, 1), -======= ->>>>>>> upstream/main } self._save_data(self.data, "psi_avg", header=f"psi_avg_{self.units}") @@ -389,12 +376,9 @@ def plot(self, **kwargs) -> str: raise AnalysisError("No psi data available to plot.") y = np.asarray(kwargs["data"], dtype=float).flatten() -<<<<<<< HEAD yerr = None if kwargs.get("std_data") is not None: yerr = np.asarray(kwargs["std_data"], dtype=float).flatten() -======= ->>>>>>> upstream/main if self.residue_indices is not None and len(self.residue_indices) == len(y): x = self.residue_indices.astype(int) else: @@ -541,10 +525,7 @@ def run(self) -> Dict[str, np.ndarray]: "omega_avg": self.data, "omega_avg_filtered": self.data, "omega_residues": self.residue_indices, -<<<<<<< HEAD "omega_std": std_angles.reshape(-1, 1), -======= ->>>>>>> upstream/main } self._save_data(self.data, "omega_avg", header=f"omega_avg_{self.units}") @@ -569,12 +550,9 @@ def plot(self, **kwargs) -> str: raise AnalysisError("No omega data available to plot.") y = np.asarray(kwargs["data"], dtype=float).flatten() -<<<<<<< HEAD yerr = None if kwargs.get("std_data") is not None: yerr = np.asarray(kwargs["std_data"], dtype=float).flatten() -======= ->>>>>>> upstream/main if self.residue_indices is not None and len(self.residue_indices) == len(y): x = self.residue_indices.astype(int) else: @@ -767,10 +745,12 @@ def plot_ramachandran( if residues is None: residues = self.residues + res_list = None if residues is not None: if isinstance(residues, int): residues = [residues] - mask = np.isin(res_indices, residues) + res_list = list(residues) + mask = np.isin(res_indices, res_list) x = x[mask] y = y[mask] res_indices = res_indices[mask] @@ -811,6 +791,46 @@ def plot_ramachandran( outpath = self._save_plot(fig, "ramachandran") plt.close(fig) + # If residues are selected, also generate per-residue frame-level plots + if res_list: + try: + _, phi_angles = md.compute_phi(self.traj) + _, psi_angles = md.compute_psi(self.traj) + if phi_angles.size and psi_angles.size: + phi_angles = phi_angles[:, res_list] + psi_angles = psi_angles[:, res_list] + + if self.units == "degrees": + phi_angles = np.degrees(phi_angles) + psi_angles = np.degrees(psi_angles) + + per_residue: Dict[int, Path] = {} + for idx, res in enumerate(res_list): + fig_res, ax_res = plt.subplots(figsize=figsize) + ax_res.scatter( + phi_angles[:, idx], + psi_angles[:, idx], + color=kwargs.get("color", "blue"), + alpha=kwargs.get("alpha", 0.7), + ) + ax_res.set_title(f"Ramachandran Plot — Residue {res}") + ax_res.set_xlabel(f"Phi ({self.units})") + ax_res.set_ylabel(f"Psi ({self.units})") + ax_res.grid(True, alpha=0.3) + fig_res.tight_layout() + + per_path = self._save_plot( + fig_res, + "ramachandran", + filename=f"ramachandran_res{res}", + ) + plt.close(fig_res) + per_residue[res] = per_path + + self.results["ramachandran_per_residue"] = per_residue + except Exception as exc: + logger.warning("Per-residue Ramachandran plots failed: %s", exc) + return outpath def plot(self, **kwargs) -> str: From 0570775e5e1576285f447ca7835d4e3c7199fa0b Mon Sep 17 00:00:00 2001 From: Derrick K Date: Wed, 11 Feb 2026 19:50:01 -0800 Subject: [PATCH 5/8] Update dihedrals.py --- src/fastmdanalysis/analysis/dihedrals.py | 48 ++++++++++++++++++++++++ 1 file changed, 48 insertions(+) diff --git a/src/fastmdanalysis/analysis/dihedrals.py b/src/fastmdanalysis/analysis/dihedrals.py index fa40b02..5e33386 100644 --- a/src/fastmdanalysis/analysis/dihedrals.py +++ b/src/fastmdanalysis/analysis/dihedrals.py @@ -231,6 +231,15 @@ def plot( mask = np.isin(x, residues) x = x[mask] y = y[mask] + if yerr is not None: + yerr = yerr[mask] + + plot_matrix = np.column_stack([x, y]) + header = f"residue_index phi_mean_{self.units}" + if yerr is not None: + plot_matrix = np.column_stack([plot_matrix, yerr]) + header = f"{header} phi_std_{self.units}" + self._save_data(plot_matrix, "phi_avg_plot", header=header) # Plot fig, ax = plt.subplots(figsize=figsize) @@ -392,6 +401,15 @@ def plot(self, **kwargs) -> str: mask = np.isin(x, residues) x = x[mask] y = y[mask] + if yerr is not None: + yerr = yerr[mask] + + plot_matrix = np.column_stack([x, y]) + header = f"residue_index psi_mean_{self.units}" + if yerr is not None: + plot_matrix = np.column_stack([plot_matrix, yerr]) + header = f"{header} psi_std_{self.units}" + self._save_data(plot_matrix, "psi_avg_plot", header=header) fig, ax = plt.subplots(figsize=kwargs.get("figsize", (12, 6))) ax.errorbar( @@ -566,6 +584,15 @@ def plot(self, **kwargs) -> str: mask = np.isin(x, residues) x = x[mask] y = y[mask] + if yerr is not None: + yerr = yerr[mask] + + plot_matrix = np.column_stack([x, y]) + header = f"residue_index omega_mean_{self.units}" + if yerr is not None: + plot_matrix = np.column_stack([plot_matrix, yerr]) + header = f"{header} omega_std_{self.units}" + self._save_data(plot_matrix, "omega_avg_plot", header=header) fig, ax = plt.subplots(figsize=kwargs.get("figsize", (12, 6))) ax.errorbar( @@ -759,6 +786,13 @@ def plot_ramachandran( if psi_std is not None: psi_std = psi_std[mask] + avg_matrix = np.column_stack([res_indices, x, y]) + header = f"residue_index phi_mean_{self.units} psi_mean_{self.units}" + if phi_std is not None and psi_std is not None: + avg_matrix = np.column_stack([avg_matrix, phi_std, psi_std]) + header = f"{header} phi_std_{self.units} psi_std_{self.units}" + self._save_data(avg_matrix, "ramachandran_avg", header=header) + fig, ax = plt.subplots(figsize=figsize) cmap = plt.get_cmap("viridis") norm = plt.Normalize(vmin=res_indices.min(), vmax=res_indices.max()) if len(res_indices) else None @@ -782,6 +816,10 @@ def plot_ramachandran( ax.set_ylabel(f"Psi ({self.units})") ax.grid(True, alpha=0.3) + limit = 180.0 if self.units == "degrees" else np.pi + ax.set_xlim(-limit, limit) + ax.set_ylim(-limit, limit) + if len(res_indices): mappable = matplotlib.cm.ScalarMappable(norm=norm, cmap=cmap) cbar = plt.colorbar(mappable, ax=ax) @@ -817,8 +855,18 @@ def plot_ramachandran( ax_res.set_xlabel(f"Phi ({self.units})") ax_res.set_ylabel(f"Psi ({self.units})") ax_res.grid(True, alpha=0.3) + ax_res.set_xlim(-limit, limit) + ax_res.set_ylim(-limit, limit) fig_res.tight_layout() + frame_matrix = np.column_stack([phi_angles[:, idx], psi_angles[:, idx]]) + frame_header = f"phi_{self.units} psi_{self.units}" + self._save_data( + frame_matrix, + f"ramachandran_res{res}", + header=frame_header, + ) + per_path = self._save_plot( fig_res, "ramachandran", From 55df74ddbb34feae16a892dff186c4ecc3127c49 Mon Sep 17 00:00:00 2001 From: Derrick K Date: Wed, 11 Feb 2026 13:22:40 -0800 Subject: [PATCH 6/8] Fix dihedral residue option aliases and add coverage - Accept residue/residue_selection aliases for phi/psi/omega and combined dihedrals without warnings - Preserve residue filtering in combined dihedrals and align warning handling with other analyses - Add tests for residue alias parsing and filtered outputs across dihedral analyses --- tests/test_dihedrals_residue_options.py | 4 ---- 1 file changed, 4 deletions(-) diff --git a/tests/test_dihedrals_residue_options.py b/tests/test_dihedrals_residue_options.py index 41a6e7e..dad6ceb 100644 --- a/tests/test_dihedrals_residue_options.py +++ b/tests/test_dihedrals_residue_options.py @@ -13,12 +13,10 @@ def test_phi_residue_alias_no_unknown_warning(fastmda): messages = [str(wi.message) for wi in w] assert not _unknown_warnings(messages) - assert analysis.data.shape[0] == 1 assert "phi_residues" in analysis.results assert list(np.asarray(analysis.results["phi_residues"]).astype(int)) == [0] - assert "phi_avg_filtered" in analysis.results filtered = analysis.results["phi_avg_filtered"] assert filtered.shape[0] == 1 @@ -32,7 +30,6 @@ def test_dihedrals_residue_selection_alias_propagates(fastmda): messages = [str(wi.message) for wi in w] assert not _unknown_warnings(messages) - # Ensure the combined analysis truly computed only 2 residues worth of data assert analysis.results["phi_avg"].shape[0] == 2 assert analysis.results["psi_avg"].shape[0] == 2 @@ -41,7 +38,6 @@ def test_dihedrals_residue_selection_alias_propagates(fastmda): assert list(np.asarray(analysis.results["phi_residues"]).astype(int)) == [0, 1] assert list(np.asarray(analysis.results["psi_residues"]).astype(int)) == [0, 1] assert list(np.asarray(analysis.results["omega_residues"]).astype(int)) == [0, 1] - for key in ("phi_avg_filtered", "psi_avg_filtered", "omega_avg_filtered"): assert key in analysis.results filtered = analysis.results[key] From b5534a5670a2ca15d26226b27d364a2600e405e3 Mon Sep 17 00:00:00 2001 From: Derrick K Date: Wed, 11 Feb 2026 17:58:11 -0800 Subject: [PATCH 7/8] Fix dihedrals --residues filtering and add tests --- tests/test_dihedrals_residue_options.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/tests/test_dihedrals_residue_options.py b/tests/test_dihedrals_residue_options.py index dad6ceb..b28fea5 100644 --- a/tests/test_dihedrals_residue_options.py +++ b/tests/test_dihedrals_residue_options.py @@ -14,7 +14,6 @@ def test_phi_residue_alias_no_unknown_warning(fastmda): messages = [str(wi.message) for wi in w] assert not _unknown_warnings(messages) assert analysis.data.shape[0] == 1 - assert "phi_residues" in analysis.results assert list(np.asarray(analysis.results["phi_residues"]).astype(int)) == [0] assert "phi_avg_filtered" in analysis.results @@ -34,7 +33,6 @@ def test_dihedrals_residue_selection_alias_propagates(fastmda): assert analysis.results["phi_avg"].shape[0] == 2 assert analysis.results["psi_avg"].shape[0] == 2 assert analysis.results["omega_avg"].shape[0] == 2 - assert list(np.asarray(analysis.results["phi_residues"]).astype(int)) == [0, 1] assert list(np.asarray(analysis.results["psi_residues"]).astype(int)) == [0, 1] assert list(np.asarray(analysis.results["omega_residues"]).astype(int)) == [0, 1] From 595458ede3b4c44f61cf6a02b0228fd54cadd565 Mon Sep 17 00:00:00 2001 From: Derrick K Date: Wed, 11 Feb 2026 19:50:01 -0800 Subject: [PATCH 8/8] Update dihedrals.py --- src/fastmdanalysis/analysis/dihedrals.py | 48 ++++++++++++++++++++++++ 1 file changed, 48 insertions(+) diff --git a/src/fastmdanalysis/analysis/dihedrals.py b/src/fastmdanalysis/analysis/dihedrals.py index fa40b02..5e33386 100644 --- a/src/fastmdanalysis/analysis/dihedrals.py +++ b/src/fastmdanalysis/analysis/dihedrals.py @@ -231,6 +231,15 @@ def plot( mask = np.isin(x, residues) x = x[mask] y = y[mask] + if yerr is not None: + yerr = yerr[mask] + + plot_matrix = np.column_stack([x, y]) + header = f"residue_index phi_mean_{self.units}" + if yerr is not None: + plot_matrix = np.column_stack([plot_matrix, yerr]) + header = f"{header} phi_std_{self.units}" + self._save_data(plot_matrix, "phi_avg_plot", header=header) # Plot fig, ax = plt.subplots(figsize=figsize) @@ -392,6 +401,15 @@ def plot(self, **kwargs) -> str: mask = np.isin(x, residues) x = x[mask] y = y[mask] + if yerr is not None: + yerr = yerr[mask] + + plot_matrix = np.column_stack([x, y]) + header = f"residue_index psi_mean_{self.units}" + if yerr is not None: + plot_matrix = np.column_stack([plot_matrix, yerr]) + header = f"{header} psi_std_{self.units}" + self._save_data(plot_matrix, "psi_avg_plot", header=header) fig, ax = plt.subplots(figsize=kwargs.get("figsize", (12, 6))) ax.errorbar( @@ -566,6 +584,15 @@ def plot(self, **kwargs) -> str: mask = np.isin(x, residues) x = x[mask] y = y[mask] + if yerr is not None: + yerr = yerr[mask] + + plot_matrix = np.column_stack([x, y]) + header = f"residue_index omega_mean_{self.units}" + if yerr is not None: + plot_matrix = np.column_stack([plot_matrix, yerr]) + header = f"{header} omega_std_{self.units}" + self._save_data(plot_matrix, "omega_avg_plot", header=header) fig, ax = plt.subplots(figsize=kwargs.get("figsize", (12, 6))) ax.errorbar( @@ -759,6 +786,13 @@ def plot_ramachandran( if psi_std is not None: psi_std = psi_std[mask] + avg_matrix = np.column_stack([res_indices, x, y]) + header = f"residue_index phi_mean_{self.units} psi_mean_{self.units}" + if phi_std is not None and psi_std is not None: + avg_matrix = np.column_stack([avg_matrix, phi_std, psi_std]) + header = f"{header} phi_std_{self.units} psi_std_{self.units}" + self._save_data(avg_matrix, "ramachandran_avg", header=header) + fig, ax = plt.subplots(figsize=figsize) cmap = plt.get_cmap("viridis") norm = plt.Normalize(vmin=res_indices.min(), vmax=res_indices.max()) if len(res_indices) else None @@ -782,6 +816,10 @@ def plot_ramachandran( ax.set_ylabel(f"Psi ({self.units})") ax.grid(True, alpha=0.3) + limit = 180.0 if self.units == "degrees" else np.pi + ax.set_xlim(-limit, limit) + ax.set_ylim(-limit, limit) + if len(res_indices): mappable = matplotlib.cm.ScalarMappable(norm=norm, cmap=cmap) cbar = plt.colorbar(mappable, ax=ax) @@ -817,8 +855,18 @@ def plot_ramachandran( ax_res.set_xlabel(f"Phi ({self.units})") ax_res.set_ylabel(f"Psi ({self.units})") ax_res.grid(True, alpha=0.3) + ax_res.set_xlim(-limit, limit) + ax_res.set_ylim(-limit, limit) fig_res.tight_layout() + frame_matrix = np.column_stack([phi_angles[:, idx], psi_angles[:, idx]]) + frame_header = f"phi_{self.units} psi_{self.units}" + self._save_data( + frame_matrix, + f"ramachandran_res{res}", + header=frame_header, + ) + per_path = self._save_plot( fig_res, "ramachandran",