Skip to content
This repository has been archived by the owner on Sep 11, 2023. It is now read-only.

how to import solvent accessible surface area as a feature in pyemma? #1586

Closed
satyajitkhatua09 opened this issue Jan 1, 2023 · 5 comments
Closed

Comments

@satyajitkhatua09
Copy link

Hi,

For one ligand unbinding, I want to use SASA as one of my CVs. Is there a way to import this feature in PyEMMA? If yes, can you tell me how to do that? This may be a simple question, but I am unable to include this feature. Please help.

WIth regards,
Satyajit Khatua

@thempel
Copy link
Member

thempel commented Jan 2, 2023

Hi Satyajit Khatua,
you can use mdtraj to compute the SASA using their implementation of the Shrake-Rupley algorithm. It's documented here. This function computes the SASA from an mdtraj-Trajectoryobject, which is exactly what you need in order to implement your own custom feature. There's an example on how to do that here #1566 (comment), can you use that as a template?

@satyajitkhatua09
Copy link
Author

satyajitkhatua09 commented Jan 2, 2023 via email

@satyajitkhatua09
Copy link
Author

Hi,

I have tried to calculate the SASA feature in these two different instances. The problems are attached herewith. I actually intended to calculate the SASA value for the ligand molecule, not all of the residues. That is where I am getting the errors. Can you help me with this? I have tried importing different options inside mdtraj.shrake_rupley() function. It seems nothing is working for me.

With regards,
Satyajit Khatua

first_instance
second_instance

@thempel
Copy link
Member

thempel commented Jan 3, 2023

The function call shrake_rupley(traj, mode='residue') returns an array of shape (n_frames, n_residues) that contains the SASA per residue for each time frame. You can select the SASA of your ligand if you know which residue-index your ligand has (make sure to use 0-based residue indexing as this is a numpy array). The code inside the function could look something like this: (CAUTION: I didn't check if it works! Please let me know if it does.)

def custom_function(traj, ligand_resid=FIXME):
    # compute ligand sasa and select ligand
    # returns array with length n_frames
    ligand_sasa = mdtraj.shrake_rupley(traj, mode='residue')[:, ligand_resid]
    
    # expand dimension to (n_frames, 1) as required by pyemma and return
    return ligand_sasa[:, None]

@satyajitkhatua09
Copy link
Author

satyajitkhatua09 commented Jan 7, 2023 via email

@thempel thempel closed this as completed Jan 9, 2023
Sign up for free to subscribe to this conversation on GitHub. Already have an account? Sign in.
Labels
None yet
Projects
None yet
Development

No branches or pull requests

2 participants