-
Notifications
You must be signed in to change notification settings - Fork 10
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
Volume vs. Time plot for a particular pocket #16
Comments
Hi @mukherjeesutanu,
To plot a volume/area of a particular cryptic pocket, you must perform cavity detection to ensure you are only getting the cavity you want in each frame. For this, you could use a space segmentation approach (box adjustment or ligand adjustment). So you could use the maximum volume or the sum of the volumes, depending on your problem. # Get all frames of the molecular dynamics simulation
frames = [f for f in sorted(os.listdir('./frames')) if f.endswith('.pdb')]
# Calculating volume and area
volarea = numpy.zeros((2, len(frames)))
# Check if cavities file already exists
if os.path.exists('cavities.pdb'):
os.remove('cavities.pdb')
# Define a common custom box using parKVFinder's PyMOL plugin
box = {
'box':{
'p1': [ 2.03, -3.97, 7.72,],
'p2': [ 18.93, -3.97, 7.72,],
'p3': [ 2.03, 16.23, 7.72,],
'p4': [ 2.03, -3.97, 33.72,],
}
}
# Write common custom box to file
with open('box.toml', 'w') as f:
toml.dump(box, f)
# Cavity detection with box adjustment
for frame in frames:
# Get model number
num = (int(frame.replace('.pdb', '')))
# Load atomic data
atomic = pyKVFinder.read_pdb(os.path.join('./frames', frame))
# Get vertices from file
vertices, atomic = pyKVFinder.get_vertices_from_file('./box.toml', atomic, probe_out=12.0)
# Detect biomolecular cavities
ncav, cavities = pyKVFinder.detect(atomic, vertices, probe_out=12.0, box_adjustment=True, volume_cutoff=50.0)
# Spatial characterization
surface, volume, area = pyKVFinder.spatial(cavities)
# Accumulate volume and area
volarea[0, num-1] = sum(volume.values())
volarea[1, num-1] = sum(area.values())
# Export cavity over time
pyKVFinder.export('./cavities.pdb', cavities, surface, vertices, append=True, model=num)
Unfortunately, it is still not possible to see the pocket volume change with nglview. Since cavities have a different number of points to represent their volume, nglview does not handle it appropriately. However, you can easily visualize them in pymol. Here is an attachment with a Jupyter Notebook with a complete example of cavity detection using box adjustment: Besides the above, we provide an example of pyKVFinder in molecular dynamics simulation at https://github.com/LBC-LNBio/pyKVFinder/blob/master/examples/md-analysis/md-analysis.ipynb and other examples at https://github.com/LBC-LNBio/pyKVFinder/blob/master/examples. Please check them too. Hope it helps you! |
Thanks for answering. I have another question related to this: if a cryptic pocket is not present in the reference.pdb structure, but it opens during the course of the simulation; how to track these kind of transient pockets using pyKVFinder? Thanks and Regards, |
In this scenario, you can perform two analyses:
First, you could detect cavities over time and save the occurrence (number of times a voxel - grid point - appears in a cavity). With it, you can export a PDB file with all cavity points over time and the occurrence of each point as a B-factor. # Create an empty array
occurrence = None
for frame in frames:
# Load atomic data
atomic = pyKVFinder.read_pdb(os.path.join('./frames', frame))
# Get vertices from file
vertices, atomic = pyKVFinder.get_vertices_from_file('./box.toml', atomic, probe_out=12.0)
# Detect biomolecular cavities
ncav, cavities = pyKVFinder.detect(atomic, vertices, probe_out=12.0, volume_cutoff=100.0, box_adjustment=True)
if occurrence is None:
occurrence = (cavities > 1).astype(int)
else:
occurrence += (cavities > 1).astype(int)
# Get cavity points
noise_cutoff = 2 # minimum number of frames that a cavity point must appear
cavities = 2 * ((occurrence >= noise_cutoff).astype('int32'))
# Export cavities with the percentage of occurrence in the B-factor column
pyKVFinder.export('./occurrence.pdb', cavities, None, vertices, B=occurrence / len(frames)) This analysis would help identify transient pockets across the biomolecular surface. Examples are available at:
After identifying the transient pockets, you can perform cavity detection around regions of interest and save the number of cavities in each frame. # Create an empty array
ncavities = []
for pdb in pdbs:
# Load atomic data
atomic = pyKVFinder.read_pdb(os.path.join('./data', pdb))
# Get vertices from file
vertices, atomic = pyKVFinder.get_vertices_from_file('./box.toml', atomic, probe_out=12.0)
# Detect biomolecular cavities
ncav, cavities = pyKVFinder.detect(atomic, vertices, probe_out=12.0, volume_cutoff=100.0, box_adjustment=True)
# Append
ncavities.append(ncav) Then, you could plot the number of cavities over time. print(ncavities)
# [1, 1, 1, 1, 2, 2, 2, 1, 1, 2]
ax = pandas.DataFrame(ncavities).plot.line()
# Customize axis
ax.set_ylim(0, 3)
ax.set_xlim(0, 9)
ax.set_ylabel('Number of cavities')
ax.set_xlabel('Time (ns)')
ax.grid(True)
# Save to file
fig = ax.get_figure()
fig.savefig('./ncavities.png', dpi=300) |
@mukherjeesutanu, I am closing this issue for now. |
How to plot the volume/area of a particular cryptic pocket(s) with respect to time? Can we see the changing pocket volume interactively on Jupyter notebook (using nglview)? Thanks.
The text was updated successfully, but these errors were encountered: