CaviTracer functionality - calculating surface cavities

[karolamik13]

I added several functionalities and improvements to CaviTracer.
Now, we are able to calculate surface cavities using calcSurfaceCavities() and display them using showSurfaceCavities().
I added improvements to calcChannels() and getVmdModel() to be able to compute and display surface cavities using the existing code.

To test the code:
from prody import *

p = parsePDB('1tqn')
protein = p.select('protein')

cavities, surface = calcSurfaceCavities(protein, output_path='test_surf_cav.pqr')
cavities2, surface2 = calcSurfaceCavities(protein, min_tetrahedra=15, output_path='test_surf_cav_min_tetra_5.pqr')
cavities3, surface3 = calcSurfaceCavities(protein, min_volume=200, max_volume=1500, output_path='test_surf_cav_MinMax_volume.pq

vmd_path = '/usr/local/bin/vmd'
model = getVmdModel(vmd_path, protein)

:::::
showSurfaceCavities(surface, model=model, show_surface=True)
showSurfaceCavities(surface2, model=model, show_surface=True)
showSurfaceCavities(surface3, model=model, show_surface=True)

showSurfaceCavities(surface3, model=model, show_surface=True, cavity_atoms='test_surf_cav_MinMax_volume.pqr')

The best visualization (visualization of the protein and cavities generated by VMD and displayed in ProDy):
model_protein = getVmdModel(vmd_path, protein)

cav_model = getVmdModel(
vmd_path,
parsePQR('test_surf_cav_MinMax_volume.pqr'),
representation='QuickSurf')

showSurfaceCavities(
surface3,
model=model_protein,
cavity_atoms=cav_model,
show_surface=True)

several other examples:
showSurfaceCavities(surface, model=model, show_surface=True)
showSurfaceCavities(surface, model=model, show_surface=True, mode='tetra')
showSurfaceCavities(surface, cavities=cavities, model=model, show_surface=True, mode='smooth', alpha=50.0, smoothing=1)
showSurfaceCavities(surface, model=model, mode='tetra', smoothing=0)

[jamesmkrieger]

I had a first try on my MacBook Pro laptop from the Crick that has macOS Tahoe 26.5.1 in my prody-github environment with python 3.10.19

Running the example code above gives me ImportError: To run getVmdModel, please install open3d but it doesn't tell me where from or what version.

I tried installing the default one from PyPI with pip, which gave me open3d-0.19.0, and it came with lots of other packages and I then get another error:

@> PDB file is found in working directory (1tqn.pdb.gz).
@> 3999 atoms and 1 coordinate set(s) were parsed in 0.01s.
@> Returning surface cavities
@> Saving surface cavities to test_surf_cav.pqr.
@> Returning surface cavities
@> Saving surface cavities to test_surf_cav_min_tetra_5.pqr.
@> Returning surface cavities
@> Saving surface cavities to test_surf_cav_MinMax_volume.pqr.
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In[1], line 11
      8 cavities3, surface3 = calcSurfaceCavities(protein, min_volume=200, max_volume=1500, output_path='test_surf_cav_MinMax_volume.pqr')
     10 vmd_path = "/Volumes/Macintosh\ HD.,juil.:\[\'/Applications/VMD\ 1.9.4a57-arm64-Rev12.app/Contents/vmd/vmd_MACOSXARM64"
---> 11 model = getVmdModel(vmd_path, protein)
     13 # :::::
     14 showSurfaceCavities(surface, model=model, show_surface=True)

File ~/scipion3/software/em/prody-github/ProDy/prody/proteins/channels.py:183, in getVmdModel(vmd_path, atoms, representation)
    180     os.unlink(temp_pdb_path)
    182 if not os.path.exists(output_path) or os.stat(output_path).st_size == 0:
--> 183     raise ValueError("STL file was not created or is empty.")
    185 stl_mesh = o3d.io.read_triangle_mesh(str(output_path))
    187 if stl_mesh.is_empty():

ValueError: STL file was not created or is empty.

[jamesmkrieger]

I had a first try on my MacBook Pro laptop from the Crick that has macOS Tahoe 26.5.1 in my prody-github environment with python 3.10.19

Running the example code above gives me ImportError: To run getVmdModel, please install open3d but it doesn't tell me where from or what version.

I tried installing the default one from PyPI with pip, which gave me open3d-0.19.0, and it came with lots of other packages and I then get another error:

@> PDB file is found in working directory (1tqn.pdb.gz).
@> 3999 atoms and 1 coordinate set(s) were parsed in 0.01s.
@> Returning surface cavities
@> Saving surface cavities to test_surf_cav.pqr.
@> Returning surface cavities
@> Saving surface cavities to test_surf_cav_min_tetra_5.pqr.
@> Returning surface cavities
@> Saving surface cavities to test_surf_cav_MinMax_volume.pqr.
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In[1], line 11
      8 cavities3, surface3 = calcSurfaceCavities(protein, min_volume=200, max_volume=1500, output_path='test_surf_cav_MinMax_volume.pqr')
     10 vmd_path = "/Volumes/Macintosh\ HD.,juil.:\[\'/Applications/VMD\ 1.9.4a57-arm64-Rev12.app/Contents/vmd/vmd_MACOSXARM64"
---> 11 model = getVmdModel(vmd_path, protein)
     13 # :::::
     14 showSurfaceCavities(surface, model=model, show_surface=True)

File ~/scipion3/software/em/prody-github/ProDy/prody/proteins/channels.py:183, in getVmdModel(vmd_path, atoms, representation)
    180     os.unlink(temp_pdb_path)
    182 if not os.path.exists(output_path) or os.stat(output_path).st_size == 0:
--> 183     raise ValueError("STL file was not created or is empty.")
    185 stl_mesh = o3d.io.read_triangle_mesh(str(output_path))
    187 if stl_mesh.is_empty():

ValueError: STL file was not created or is empty.

This error doesn't come up on a linux machine with vmd loaded from modules

[jamesmkrieger]

It all looks reasonable to me. I have few changes but they are optional

[karolamik13]

I forgot to add calcSurfaceCavitiesMultipleFrames(). It has now been added to the ProDy code.

Now, cavities can be analyzed using a multi-model PDB or a trajectory via one function.

Multi-model PDB:
from prody import *
p = parsePDB('multi.pdb').select('protein')
cavities, surface = calcSurfaceCavitiesMultipleFrames(p, output_path='_mm_test_CAV_1')

Trajectory:
from prody import *
PDBfile = 'protein.pdb'
DCDfile = 'traj.dcd'
atoms = parsePDB(PDBfile)
dcd = Trajectory(DCDfile)
dcd.link(atoms)
dcd.setCoords(atoms)
calcSurfaceCavitiesMultipleFrames(atoms, dcd, output_path='_m_test_CAV_1')

[karolamik13]

I added two functions similarly to the channels to provide information about cavities and residues.

To test it (and a few other functions before) execute the following commands:
from prody import *
p = parsePDB('1tqn')
protein = p.select('protein')
cavities, surface = calcSurfaceCavities(protein, output_path='test_surf_cav.pqr')
cavities3, surface3 = calcSurfaceCavities(protein, min_volume=500, max_volume=1000, output_path='test_surf_cav_MinMax_volume.pqr')

vmd_path = '/usr/local/bin/vmd'
model = getVmdModel(vmd_path, protein)

showSurfaceCavities(surface, model=model, show_surface=True)
showSurfaceCavities(surface3, model=model, show_surface=True)

model_protein = getVmdModel(vmd_path, protein)
cav_model = getVmdModel(
vmd_path,
parsePQR('test_surf_cav_MinMax_volume.pqr'),
representation='QuickSurf'
)

showSurfaceCavities(
surface3,
model=model_protein,
cavity_atoms=cav_model,
show_surface=True
)

results1 = getSurfaceCavityParameters(cavities3)
results2 = getSurfaceCavityResidueNames(protein, cavities3, surface3, residues_file_name='results')

You will obtain:

...
@> Model created successfully.
@> Cavity ID: Volume [ų] Depth [Å] Tetrahedra count
@> cavity 0: 928.09 8 97
@> cavity 1: 771.55 5 77

In [2]: results1
Out[2]: ([928.0878561950001, 771.5530158881667], [8, 5], [97, 77])

In [3]: results2
Out[3]:
['cavity0: LYS173, ASP174, GLY177, ALA178, SER195, LEU196, PRO199, LYS208, LYS209, LEU211, ARG212, PHE213, PHE219, PHE220, CYS239, VAL240, PHE241, PRO242, PHE304, TYR307, GLU308, SER311, LEU482',
'cavity1: LYS55, MET59, MET62, GLU63, TYR319, GLU320, THR323, HIS324, PHE367, TYR399, ARG403, GLU412, LYS413, PHE414, ILE473, PRO474, LEU475, LEU477, SER478, LEU479']

Provided residues are displayed in VMD (attached JPG file).

[karolamik13]

I am not sure what I can do with that. Maybe the version of something was different.
I am using Python 3.10.19 and open3d-0.19.0 as well, but on Linux (Ubuntu).

I had a first try on my MacBook Pro laptop from the Crick that has macOS Tahoe 26.5.1 in my prody-github environment with python 3.10.19

Running the example code above gives me ImportError: To run getVmdModel, please install open3d but it doesn't tell me where from or what version.

I tried installing the default one from PyPI with pip, which gave me open3d-0.19.0, and it came with lots of other packages and I then get another error:

@> PDB file is found in working directory (1tqn.pdb.gz).
@> 3999 atoms and 1 coordinate set(s) were parsed in 0.01s.
@> Returning surface cavities
@> Saving surface cavities to test_surf_cav.pqr.
@> Returning surface cavities
@> Saving surface cavities to test_surf_cav_min_tetra_5.pqr.
@> Returning surface cavities
@> Saving surface cavities to test_surf_cav_MinMax_volume.pqr.
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In[1], line 11
      8 cavities3, surface3 = calcSurfaceCavities(protein, min_volume=200, max_volume=1500, output_path='test_surf_cav_MinMax_volume.pqr')
     10 vmd_path = "/Volumes/Macintosh\ HD.,juil.:\[\'/Applications/VMD\ 1.9.4a57-arm64-Rev12.app/Contents/vmd/vmd_MACOSXARM64"
---> 11 model = getVmdModel(vmd_path, protein)
     13 # :::::
     14 showSurfaceCavities(surface, model=model, show_surface=True)

File ~/scipion3/software/em/prody-github/ProDy/prody/proteins/channels.py:183, in getVmdModel(vmd_path, atoms, representation)
    180     os.unlink(temp_pdb_path)
    182 if not os.path.exists(output_path) or os.stat(output_path).st_size == 0:
--> 183     raise ValueError("STL file was not created or is empty.")
    185 stl_mesh = o3d.io.read_triangle_mesh(str(output_path))
    187 if stl_mesh.is_empty():

ValueError: STL file was not created or is empty.

[karolamik13]

Tomorrow I will write a tutorial for the ProDy website to show the functionalities.

[karolamik13]

In the previous version of start_point for calcChannels(), I forgot to add an option to provide a ProDy selection instead of a 3D coordinate point. Now it is added.

from prody import *
p = parsePDB('1tqn')
protein = p.select('protein')

-> Start point via coordinates
calcChannels(protein, output_path='test7x.pdb', start_point=[-22.312, -20.065, -11.144])

-> Start point via protein selection (if more than one atom, a center of mass is computed):
start_sel = protein.select('resid 212 309 483')
calcChannels(protein, output_path='test7xx.pdb', start_point=start_sel)

I am attaching a comparison of those predictions (in orange: resid 212 309 483 are displayed).