Binding Pocket Analysis#
ProteoTwirl can be used not only to study the exterior of a protein as a whole, but also to analyze internal regions of a protein, such as binding pockets.
In this tutorial, we will analyze the hydrophobicity of the binding-pocket surface in human carbonic anhydrase II (HCA II). The HCA II binding pocket contains notable hydrophobic and hydrophilic patches of residues that are important for the function of the protein (see this article for more information).
We will use ProteoTwirl to analyze the binding-pocket surface and visualize
the hydrophobicity of the surface residues. For this tutorial, we will use
the 6GDC structure, which can be downloaded from here.
Loading and Centering the Protein#
We need to load the protein structure into ProteoTwirl; however, we want to center the protein structure on the binding pocket. To do this, we first need to find the center of the binding pocket. Since a ligand is present in the structure, we can use its location as the binding-pocket center.
PyMOL was used here to visually identify the ligand and then read the residue name and residue number. We need the residue number because there are two ligands present in the structure, but only one of them is in the binding pocket.
import MDAnalysis as mda
import numpy as np
# Load in the structure in MDAnalysis
universe = mda.Universe("6GDC.pdb")
# Select the ligand
ligand = universe.select_atoms("resname FB2 and resid 305")
# Get the position of ligand atoms
ligand_pos = ligand.positions
# Calculate the average position (center of geometry (COG))
ligand_cog = np.mean(ligand_pos, axis = 0)
print(ligand_cog)
>>> [-26.5358 2.6104002 15.9722 ]
Then we can load the protein centered on this position by disabling
center_selection and setting center = ligand_cog.
import ProteoTwirl.ProteinRotation as pt
protein_structure = pt.ProteinRotation(
protein = "6GDC.pdb",
selection = "protein and not element H",
base_rotation = None,
center_selection = False,
center = ligand_cog,
)
protein_structure.show()
Creating a Binding-Pocket Selection#
To select the surface of the binding pocket, we will combine three different selectors: the clipping-plane, field-of-view, and exposed-atom selectors.
The first one is the clipping-plane selector. We use this to “cut” the protein in half, so we do not see all the atoms in front of the binding pocket, blocking our view of it. The clipping plane will select anything behind the center (the binding pocket, because we have centered the protein on it) and looks like this.
from ProteoTwirl.Selectors.Selector import ClippingPlaneSelector
protein_structure.iter(
# Note we are using relatively few angles for the visualization here
rotation_res = 64,
selectors = [ClippingPlaneSelector((None,0))]
).show_3d(background = "#14181e")
Next, we will apply a field-of-view selector. This limits the atoms we are looking at to those within a small radius from the viewing direction, essentially selecting atoms within a cylinder. This is needed so that we do not select atoms outside the binding pocket.
from ProteoTwirl.Selectors.Selector import ClippingPlaneSelector, FieldOfViewSelector
protein_structure.iter(
# Note we are using relatively few angles for the visualization here
rotation_res = 64,
selectors = [
ClippingPlaneSelector((None,0)),
FieldOfViewSelector(7),
]
).show_3d(background = "#14181e")
Finally, we will use an exposed-atom selector to select only the atoms that are present on the surface of the binding pocket and exposed for a ligand to interact with. Combining all three selectors gives:
from ProteoTwirl.Selectors.Selector import ExposedAtomSelector, ClippingPlaneSelector, FieldOfViewSelector
protein_structure.iter(
# Note we are using relatively few angles for the visualization here
rotation_res = 64,
selectors = [
ClippingPlaneSelector((None,0)),
FieldOfViewSelector(7),
ExposedAtomSelector()
]
).show_3d(background = "#14181e")
Calculating Hydrophobicity#
Using the selector setup described above, we can now calculate the hydropathy index of the selected residues and calculate the mean to identify hydrophobic patches.
protein_structure.add_property(
property = "count.Hydropathy",
rotation_res = 256,
agg_method = "mean",
use_default = False,
rotation_range = None,
selectors = [
ClippingPlaneSelector((None,0)),
FieldOfViewSelector(7),
ExposedAtomSelector()
]
)
Visualizing Results#
We can visualize the results using plot_carto_contour to show which angles
select the most hydrophobic residues.
from ProteoTwirl.Visualization.CatographyPlot import plot_carto_contour
import cartopy.crs as ccrs
ax = plot_carto_contour(
angles = protein_structure._angle_spaces[0],
values = protein_structure.mean_Hydropathy,
projection = ccrs.Robinson(),
cmap = "BuPu",
figsize = (13,7),
)
We can alternatively make an interactive visualization where, by setting the protein
color to "Hydropathy Index", we can also color the residues according to their
hydropathy.
protein_structure.plot_property_3d(
property = "mean_Hydropathy",
protein_color = "Hydropathy Index",
cmap = "plasma",
background = "#14181e",
)
A more detailed and complex visualization of this can be found as an example in the gallery: gallery6gdc.
Confirming Results#
Does the hydrophobic patch we found agree with the residues that were identified in the article? To investigate this, we can check which residues are visible from each angle. We can then count how many of those residues match the ones noted in the article and check whether the angles with the most residue overlap are also the ones with the highest hydrophobicity.
First, we make a RotationCounter with the same settings as used for
calculating the hydropathy index. Then we can loop over each angle, which
provides positions, an atom selection, and atom indices for that angle, and
just keep the atom selection.
prot_it = protein_structure.iter(
rotation_res = 256,
selectors = [
ClippingPlaneSelector((None,0)),
FieldOfViewSelector(7),
ExposedAtomSelector()
]
)
md_atom_selection = []
for it_positions, it_atom, it_indices in prot_it:
md_atom_selection.append(it_atom)
After this, we will check how many of the residues in each angle selection match the residues noted in the article.
# Hydrophobic-patch residues from the article:
# Val121, Val143, Leu198, Val207, and Trp209
hydrophobic_patch = [121, 143, 198, 207, 209]
# Get residue indices for each angle
res_index = [mi.resids for mi in md_atom_selection]
# Get the residues in common between angle selection and article residues
overlap = [set(ri).intersection(hydrophobic_patch) for ri in res_index]
# Finally get the size of this overlap
overlap_count = [float(len(ov)) for ov in overlap]
We can then visualize the similarity between the overlap with the article and the calculated hydropathy. The 2D plot below shows good agreement between the hydropathy signal identified in the binding pocket and the residues identified in the article.
from ProteoTwirl.Visualization.CatographyPlot import plot_carto_contour
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
fig, axes = plt.subplots(
2, 1,
figsize = (5, 6),
subplot_kw = {"projection": ccrs.Robinson()}
)
ax1 = plot_carto_contour(
angles = protein_structure._angle_spaces[0],
values = protein_structure.mean_Hydropathy,
projection = ccrs.Robinson(),
cmap = "BuPu",
ax = axes[0]
)
ax1.set_title("Mean Hydropathy", color="#5555dd")
ax2 = plot_carto_contour(
angles = protein_structure._angle_spaces[0],
values = overlap_count,
projection = ccrs.Robinson(),
cmap = "BuPu",
ax = axes[1]
)
ax2.set_title("Number of Overlapping Residues", color="#5555dd")
