Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Solvant Accesible surface area #4025

Closed
wants to merge 13 commits into from
Closed

Conversation

pegerto
Copy link

@pegerto pegerto commented Feb 13, 2023

Fixes #2439

Changes made in this Pull Request:

  • Add an implementation of Shrake-Rupley algorithm to calculate SASA per atom.

PR Checklist

  • Tests?
  • Docs?
  • CHANGELOG updated?
  • Issue raised/referenced?

Copy link

@github-actions github-actions bot left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Hello there first time contributor! Welcome to the MDAnalysis community! We ask that all contributors abide by our Code of Conduct and that first time contributors introduce themselves on the developer mailing list so we can get to know you. You can learn more about participating here. Please also add yourself to package/AUTHORS as part of this PR.

@pegerto pegerto changed the title basic sasa implementation Solvant Accesible surface area Feb 13, 2023
@codecov
Copy link

codecov bot commented Feb 13, 2023

Codecov Report

Base: 93.53% // Head: 93.54% // Increases project coverage by +0.01% 🎉

Coverage data is based on head (b9bef46) compared to base (d27a32a).
Patch coverage: 100.00% of modified lines in pull request are covered.

Additional details and impacted files
@@             Coverage Diff             @@
##           develop    #4025      +/-   ##
===========================================
+ Coverage    93.53%   93.54%   +0.01%     
===========================================
  Files          191      192       +1     
  Lines        25062    25108      +46     
  Branches      3547     3550       +3     
===========================================
+ Hits         23442    23488      +46     
  Misses        1099     1099              
  Partials       521      521              
Impacted Files Coverage Δ
package/MDAnalysis/analysis/accesibility.py 100.00% <100.00%> (ø)

Help us with your feedback. Take ten seconds to tell us how you rate us. Have a feature suggestion? Share it here.

☔ View full report at Codecov.
📢 Do you have feedback about the report comment? Let us know in this issue.

@orbeckst
Copy link
Member

Hello @pegerto , thank you for your contribution.

Because this project is listed as a potential GSOC project I have to ask if you intend to apply for Google Summer of Code with MDAnalysis (if MDAnalysis is selected to take part this year).

@pegerto
Copy link
Author

pegerto commented Feb 13, 2023

Google Summer of Code

Hi @orbeckst I saw the documentation of this request listed under GSOC 2021, are you planning to present this project for 2023 as organisation. It was not my intention to do so rather contribute to this analysis module in parallel to my research.

@orbeckst
Copy link
Member

Ok, good, thank you for letting us know!

@pegerto pegerto force-pushed the sasa branch 2 times, most recently from 001ae79 to 58d9cfc Compare February 18, 2023 11:35
@pegerto pegerto marked this pull request as ready for review February 18, 2023 19:55
Copy link
Member

@RMeli RMeli left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Hi @pegerto , a few initial comments from skimming through (I haven't looked at the core of the code yet).

I also think it might be better to call the module sasa.py or something similar, accessibility.py is too general IMO (but others might disagree).

Comment on lines +56 to +59
def atomic_radi(e):
"""van der Waals radii"""
DEFAULT_ATOMIC_RADI = 2.0
return ATOMIC_RADI[e] if e in ATOMIC_RADI else DEFAULT_ATOMIC_RADI
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This looks superfluous. You can use the dictionary's get method with a default argument, or maybe even a defaultdict.

Comment on lines +31 to +53
ATOMIC_RADI = {
"H": 1.200,
"HE": 1.400,
"C": 1.700,
"N": 1.550,
"O": 1.520,
"F": 1.470,
"NA": 2.270,
"MG": 1.730,
"P": 1.800,
"S": 1.800,
"CL": 1.750,
"K": 2.750,
"CA": 2.310,
"NI": 1.630,
"CU": 1.400,
"ZN": 1.390,
"SE": 1.900,
"BR": 1.850,
"CD": 1.580,
"I": 1.980,
"HG": 1.550,
}
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It would be good to add a comment on where you got these numbers from. Even better, you can use the ones already available in MDAnalysis:

TABLE_VDWRADII = r"""
# Van der Waals radii taken from
# [1] Bondi, A. (1964). "Van der Waals Volumes and Radii".
# J. Phys. Chem. 68 (3): 441-451. doi:10.1021/j100785a001.
# [2] Rowland and Taylor (1996). "Intermolecular Nonbonded Contact Distances in Organic Crystal Structures:
# Comparison with Distances Expected from van der Waals Radii".
# J. Phys. Chem., 1996, 100 (18), 7384.7391. doi:10.1021/jp953141+.
# [3] Mantina, et al. (2009). "Consistent van der Waals Radii for the Whole Main Group".
# J. Phys. Chem. A, 2009, 113 (19), 5806-5812. doi:10.1021/jp8111556.

n_points: int (optional)
Number of points used for the estimation of the atom surface area
probe_radius: double ( optional )
Radius of the probe atom, by default watter radious
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
Radius of the probe atom, by default watter radious
Radius of the probe atom, by default water radious

Atom group used for the calculation.
n_points: int (optional)
Number of points used for the estimation of the atom surface area
probe_radius: double ( optional )
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
probe_radius: double ( optional )
probe_radius: double (optional)

DEFAULT_ATOMIC_RADI = 2.0
return ATOMIC_RADI[e] if e in ATOMIC_RADI else DEFAULT_ATOMIC_RADI


Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Can you please add the DOI(s) of the relevant publication(s)/source(s) on which this analysis class is based? People need such references to double check the method and see where the default values come from.

You can find an example of how to do that here:

due.cite(Doi("10.21105/joss.00877"),
description="Mean Squared Displacements with tidynamics",
path="MDAnalysis.analysis.msd",
cite_module=True)
due.cite(Doi("10.1051/sfn/201112010"),
description="FCA fast correlation algorithm",
path="MDAnalysis.analysis.msd",
cite_module=True)

from MDAnalysis.analysis.accesibility import SASA


class TestSASA(object):
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Are there simple cases for which an analytical result is available?

self._atom_neighbours = KDTree(ag.positions, 10)
self._n_points = n_points
self._sphere = self._compute_sphere()
self._twice_max_radi = (max(list(ATOMIC_RADI.values())) + probe_radius) * 2
Copy link
Member

@RMeli RMeli Feb 18, 2023

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is there a better name for this variable, since it is not twice the maximum radius?


for atom in self._ag.atoms:
r_w = atomic_radi(atom.type)
raddi = r_w + self._probe_radius
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
raddi = r_w + self._probe_radius
radii = r_w + self._probe_radius

But maybe radius (singular) is more pertinent here? See also the other occurrences.

Copy link
Member

@IAlibay IAlibay left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks for this @pegerto

I however have two concerns:

  1. The _compute and _compute_sphere methods are very reminiscent of the way classes are defined in Biopython. I haven't had a chance to do a direct comparison of the class content here, so please forgive me rudely asking - how much, if any, of this code comes from Biopython?
  2. The AnalysisBase class methods aren't being used properly, I would suggest having a look at the docs for how to create such a child class: https://userguide.mdanalysis.org/stable/examples/analysis/custom_trajectory_analysis.html

return ATOMIC_RADI[e] if e in ATOMIC_RADI else DEFAULT_ATOMIC_RADI


class SASA(AnalysisBase):
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'm not seeing anything being defined for _prepare, _single_frame, or _conclude, and so it looks like you're not using the AnalysisBase class properly to allow it to iterate over a trajectory.

@pegerto
Copy link
Author

pegerto commented Feb 20, 2023

  1. The _compute and _compute_sphere methods are very reminiscent of the way classes are defined in Biopython. I haven't had a chance to do a direct comparison of the class content here, so please forgive me rudely asking - how much, if any, of this code comes from Biopython?

Hey @IAlibay not problem at all I definitely looked at the biopython implementation, the idea of precompute the sphere points and then stretch it and translate centroid of the atoms comes from this implementations.

Thanks for the feedback

@orbeckst
Copy link
Member

Do you have an sense of relative performance compared to eg freesasa?

orbeckst added a commit to MDAnalysis/MDAnalysis.github.io that referenced this pull request Feb 22, 2023
We have an active PR open for SASA. MDAnalysis/mdanalysis#4025
@orbeckst
Copy link
Member

@pegerto if your code is based on biopython's Bio.PDB.SASA, wouldn't it be easier to just import the code and work on moving data between MDA and Biopython? How much of a cost does one incur doing so?

orbeckst added a commit to MDAnalysis/MDAnalysis.github.io that referenced this pull request Feb 22, 2023
* GSoC 2023 blog post draft

* Apply suggestions from code review

* add gsoc list to config

* remove SASA project

We have an active PR open for SASA. MDAnalysis/mdanalysis#4025

* embed MDAnalysis 2021 trailer

* post updates

  - fix project numbering (now same as in https://github.com/MDAnalysis/mdanalysis/wiki/GSoC-2023-Project-Ideas)
  - fix broken link constructs to GSoC list
  - add final remarks subheading for better showing difference between requirements and comments
  - add Discord as contact to last paragraph
  - close with a friendly sign off
  - use GitHub and Discord team handles in signature

* fixed Discord gsoc-mentor handle

* removed unavailable mentors

---------

Co-authored-by: Rocco Meli <[email protected]>
Co-authored-by: Lily Wang <[email protected]>
@IAlibay
Copy link
Member

IAlibay commented Jul 11, 2023

@pegerto is this still something you are interested in contributing?

@IAlibay
Copy link
Member

IAlibay commented Jul 31, 2023

I'm going to close this PR as stale, please do go ahead and re-open if you still wish to contribute this @pegerto, it is still a very much needed feature.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
None yet
Development

Successfully merging this pull request may close these issues.

Molecular volume and (solvent-accessible) surface area
5 participants