schrodinger / pymol-open-source

Open-source foundation of the user-sponsored PyMOL molecular visualization system.
https://pymol.org/
Other
1.15k stars 275 forks source link

get_area result does not correspond with SASA calculated using other libraries #302

Closed Glucksistemi closed 1 year ago

Glucksistemi commented 1 year ago

There is a script to test the situation described in header witch uses Pymol, Biopython and Biotite to get the same information:

from pymol import cmd

from Bio.PDB import PDBParser, SASA

from biotite import structure
from biotite.structure.io import load_structure

PROBE_FILE = '/home/gluck/PTMv2/aligned/wt/2UZB_C-2WIP_A.pdb' # replace with your local PDB file here

# Using PyMol
cmd.load(PROBE_FILE)
pm_sasa = cmd.get_area()
print(f"Pymol SASA: {pm_sasa}")
# Using Biopython
parser = PDBParser(QUIET=True)
struct = parser.get_structure('test', PROBE_FILE)
sasa_conf = SASA.ShrakeRupley()
sasa_conf.compute(struct, level='C')
bp_sasa = struct[0].get_chains().__next__().sasa
print(f"Biopython SASA: {bp_sasa}")
# Using biotite
bt_atoms = load_structure(PROBE_FILE)
bt_sasa = structure.sasa(bt_atoms).sum()
print(f"Biotite SASA: {bt_sasa}")

The output is:

Pymol SASA: 30262.603515625
Biopython SASA: 14397.994692988526
Biotite SASA: 14210.21484375

I tried it with several files, all containing a single chain without ligands and solvent (previously cleaned using remove not polymer), and always the result given py PyMol is greater than with other two ways to calculate. Commonly it is about two times higher, but I found no real dependency, for some structures difference is lower than two times.

If there were smaller difference I would belive this is a question of accuracy or default settings, like in case of biopython vs biotite, but this difference I got has no excuse about accuracy.

Or I assume I can be wrong and result of get_area is not exaxtly SASA (solvent-accessible surface area), so please fix me and recommend a way to get exactly SASA value for a whole structure/selection using PyMol

speleo3 commented 1 year ago

Short answer: You need to set dot_solvent.

cmd.set("dot_solvent")
pm_sasa = cmd.get_area()

From the help get_area output:

Get the surface area of an selection. Depends on the "dot_solvent" setting. With "dot_solvent=off" (default) it calculates the solvent excluded surface area, else the surface accessible surface.

JarrettSJohnson commented 1 year ago

Closing since OP reacted positively to the last comment--assuming this fixes your issue.

Thanks, Thomas.