From f0722d32ac7331c581487d7bd48a03bd23425e1a Mon Sep 17 00:00:00 2001 From: pegerto-fernandez-ck Date: Wed, 4 Oct 2023 23:08:51 +0100 Subject: [PATCH 1/3] Calcuate SASA using free sasa --- mdakit_sasa/analysis/sasaanalysis.py | 44 ++++++------------- .../tests/analysis/test_sasaanalysis.py | 35 +++++---------- mdakit_sasa/tests/test_mdakit_sasa.py | 9 +--- mdakit_sasa/tests/utils.py | 11 +---- pyproject.toml | 1 + 5 files changed, 26 insertions(+), 74 deletions(-) diff --git a/mdakit_sasa/analysis/sasaanalysis.py b/mdakit_sasa/analysis/sasaanalysis.py index 40cf789..33854be 100644 --- a/mdakit_sasa/analysis/sasaanalysis.py +++ b/mdakit_sasa/analysis/sasaanalysis.py @@ -9,6 +9,7 @@ from MDAnalysis.analysis.base import AnalysisBase import numpy as np +import freesasa if TYPE_CHECKING: from MDAnalysis.core.universe import Universe, AtomGroup @@ -74,42 +75,23 @@ def __init__( def _prepare(self): """Set things up before the analysis loop begins""" - # This is an optional method that runs before - # _single_frame loops over the trajectory. - # It is useful for setting up results arrays - # For example, below we create an array to store - # the number of atoms with negative coordinates - # in each frame. - self.results.is_negative = np.zeros( + self.results.total_area = np.zeros( (self.n_frames, self.atomgroup.n_atoms), dtype=bool, ) def _single_frame(self): """Calculate data from a single frame of trajectory""" - # This runs once for each frame of the trajectory - # It can contain the main analysis method, or just collect data - # so that analysis can be done over the aggregate data - # in _conclude. - - # The trajectory positions update automatically - negative = self.atomgroup.positions < 0 - # You can access the frame number using self._frame_index - self.results.is_negative[self._frame_index] = negative.any(axis=1) + + structure = freesasa.Structure() + # FreeSasa structure accepts PDBS if not available requires to reconstruct the structure using `adAtom` + for a in self.atomgroup: + x,y,z = a.position + structure.addAtom(a.name, a.resname, a.resnum.item(), "", x, y, z) + + result = freesasa.calc(structure) + + self.results.total_area[self._frame_index] = result.totalArea() def _conclude(self): - """Calculate the final results of the analysis""" - # This is an optional method that runs after - # _single_frame loops over the trajectory. - # It is useful for calculating the final results - # of the analysis. - # For example, below we determine the - # which atoms always have negative coordinates. - self.results.always_negative = self.results.is_negative.all(axis=0) - always_negative_atoms = self.atomgroup[self.results.always_negative] - self.results.always_negative_atoms = always_negative_atoms - self.results.always_negative_atom_names = always_negative_atoms.names - - # results don't have to be arrays -- they can be any value, e.g. floats - self.results.n_negative_atoms = self.results.is_negative.sum(axis=1) - self.results.mean_negative_atoms = self.results.n_negative_atoms.mean() + self.results.mean_total_area= self.results.total_area.mean() diff --git a/mdakit_sasa/tests/analysis/test_sasaanalysis.py b/mdakit_sasa/tests/analysis/test_sasaanalysis.py index 07bb366..daad37a 100644 --- a/mdakit_sasa/tests/analysis/test_sasaanalysis.py +++ b/mdakit_sasa/tests/analysis/test_sasaanalysis.py @@ -3,7 +3,7 @@ from mdakit_sasa.analysis.sasaanalysis import SASAAnalysis from mdakit_sasa.tests.utils import make_Universe - +from MDAnalysis.core.topologyattrs import Atomnames, Resnames, Resids, Resnums class TestSASAAnalysis: @@ -12,12 +12,16 @@ class TestSASAAnalysis: @pytest.fixture def universe(self): u = make_Universe( - extras=("names", "resnames",), n_frames=3, ) - # create toy data to test assumptions + for ts in u.trajectory: ts.positions[:ts.frame] *= -1 + + u.add_TopologyAttr(Atomnames(["H"] * len(u.atoms))) + u.add_TopologyAttr(Resnames(["GLY"] * len(u.residues))) + u.add_TopologyAttr(Resids(list(range(0, len(u.residues))))) + u.add_TopologyAttr(Resnums(list(range(0, len(u.residues))))) return u @pytest.fixture @@ -38,26 +42,7 @@ def test_atom_selection(self, universe, select, n_atoms): universe, select=select) assert analysis.atomgroup.n_atoms == n_atoms - @pytest.mark.parametrize( - "stop, expected_mean", - [ - (1, 0), - (2, 0.5), - (3, 1) - ] - ) - def test_mean_negative_atoms(self, analysis, stop, expected_mean): - # assert we haven't run yet and the result doesn't exist yet + def test_residue_sasa_calculation(self, analysis): assert "mean_negative_atoms" not in analysis.results - analysis.run(stop=stop) - assert analysis.n_frames == stop - - # when comparing floating point values, it's best to use assert_allclose - # to allow for floating point precision differences - assert_allclose( - analysis.results.mean_negative_atoms, # computed data - expected_mean, # reference / desired data - rtol=1e-07, # relative tolerance - atol=0, # absolute tolerance - err_msg="mean_negative_atoms is not correct", - ) + analysis.run(stop=3) + assert analysis.n_frames == 3 \ No newline at end of file diff --git a/mdakit_sasa/tests/test_mdakit_sasa.py b/mdakit_sasa/tests/test_mdakit_sasa.py index 1485240..6a191f7 100644 --- a/mdakit_sasa/tests/test_mdakit_sasa.py +++ b/mdakit_sasa/tests/test_mdakit_sasa.py @@ -2,7 +2,6 @@ Unit and regression test for the mdakit_sasa package. """ -# Import package, test suite, and other packages as needed import mdakit_sasa import pytest import sys @@ -10,10 +9,4 @@ def test_mdakit_sasa_imported(): """Sample test, will always pass so long as import statement worked""" - assert "mdakit_sasa" in sys.modules - - -def test_mdanalysis_logo_length(mdanalysis_logo_text): - """Example test using a fixture defined in conftest.py""" - logo_lines = mdanalysis_logo_text.split("\n") - assert len(logo_lines) == 46, "Logo file does not have 46 lines!" + assert "mdakit_sasa" in sys.modules \ No newline at end of file diff --git a/mdakit_sasa/tests/utils.py b/mdakit_sasa/tests/utils.py index 63a2e16..f5eeb31 100644 --- a/mdakit_sasa/tests/utils.py +++ b/mdakit_sasa/tests/utils.py @@ -6,7 +6,6 @@ def make_Universe( - extras: Tuple[str] = tuple(), size: Tuple[int, int, int] = (125, 25, 5), n_frames: int = 0, velocities: bool = False, @@ -23,10 +22,6 @@ def make_Universe( Parameters ---------- - extras : tuple of strings, optional - extra attributes to add to Universe: - u = make_Universe(('masses', 'charges')) - Creates a lightweight Universe with only masses and charges. size : tuple of int, optional number of elements of the Universe (n_atoms, n_residues, n_segments) n_frames : int @@ -58,11 +53,7 @@ def make_Universe( velocities=velocities, forces=forces, ) - if extras is None: - extras = [] - for ex in extras: - u.add_TopologyAttr(ex) - + if trajectory: pos = np.arange(3 * n_atoms * n_frames).reshape(n_frames, n_atoms, 3) vel = pos + 100 if velocities else None diff --git a/pyproject.toml b/pyproject.toml index f220339..11b9179 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -19,6 +19,7 @@ readme = "README.md" requires-python = ">=3.9" dependencies = [ "MDAnalysis>=2.0.0", + "freesasa>= 2.2.0" ] keywords = [ "molecular simulations", From bdeb20a9b7d6187fb293088244727db771ab3c07 Mon Sep 17 00:00:00 2001 From: pegerto-fernandez-ck Date: Wed, 4 Oct 2023 23:16:53 +0100 Subject: [PATCH 2/3] wip --- mdakit_sasa/analysis/sasaanalysis.py | 11 ++--------- 1 file changed, 2 insertions(+), 9 deletions(-) diff --git a/mdakit_sasa/analysis/sasaanalysis.py b/mdakit_sasa/analysis/sasaanalysis.py index 33854be..7278693 100644 --- a/mdakit_sasa/analysis/sasaanalysis.py +++ b/mdakit_sasa/analysis/sasaanalysis.py @@ -18,7 +18,8 @@ class SASAAnalysis(AnalysisBase): """SASAAnalysis class. - This class is used to perform analysis on a trajectory. + This class is used to compute the solvant accessible area of a trajectory. +. Parameters ---------- @@ -59,17 +60,9 @@ def __init__( self, universe_or_atomgroup: Union["Universe", "AtomGroup"], select: str = "all", - # TODO: add your own parameters here **kwargs ): - # the below line must be kept to initialize the AnalysisBase class! super().__init__(universe_or_atomgroup.trajectory, **kwargs) - # after this you will be able to access `self.results` - # `self.results` is a dictionary-like object - # that can should used to store and retrieve results - # See more at the MDAnalysis documentation: - # https://docs.mdanalysis.org/stable/documentation_pages/analysis/base.html?highlight=results#MDAnalysis.analysis.base.Results - self.universe = universe_or_atomgroup.universe self.atomgroup = universe_or_atomgroup.select_atoms(select) From b599976002fcf87a282689341a9e41b3b49fe535 Mon Sep 17 00:00:00 2001 From: pegerto-fernandez-ck Date: Wed, 4 Oct 2023 23:27:47 +0100 Subject: [PATCH 3/3] wip --- devtools/conda-envs/test_env.yaml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/devtools/conda-envs/test_env.yaml b/devtools/conda-envs/test_env.yaml index bea319c..52fa0db 100644 --- a/devtools/conda-envs/test_env.yaml +++ b/devtools/conda-envs/test_env.yaml @@ -7,8 +7,8 @@ dependencies: - python - pip - # MDAnalysis - MDAnalysis + - freesasa # Testing - pytest