Merge pull request #1 from pegerto/basic_freesasa

Calcuate SASA using free sasa

devtools/conda-envs/test_env.yaml

@@ -7,8 +7,8 @@ dependencies:
   - python
   - pip
 
-  # MDAnalysis
   - MDAnalysis
+  - freesasa
 
   # Testing
   - pytest

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
@@ -17,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
     ----------
@@ -58,58 +60,31 @@ 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)
 
     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()

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

mdakit_sasa/tests/test_mdakit_sasa.py

@@ -2,18 +2,11 @@
 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
 
 
 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

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

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",