Contact Map Explorer includes some tricks to study contact maps in protein dynamics, based on tools in MDTraj. This notebook shows examples and serves as documentation.
As an example, we'll use part of a trajectory of the KRas protein bound to GTP, which was provided by Sander Roet. KRas is a protein that plays a role in many cancers. For simplicity, the waters were removed from the trajectory (although ions are still included). To run this notebook, download the example files from https://figshare.com/s/453b1b215cf2f9270769 (total download size about 1.2 MB). Download all files, and extract in the same directory that you started Jupyter from (so that you have a directory called 5550217 in your current working directory).
In [1]:
%matplotlib inline
import matplotlib.pyplot as plt
import mdtraj as md
traj = md.load("5550217/kras.xtc", top="5550217/kras.pdb")
topology = traj.topology
In [2]:
from contact_map import ContactMap, ContactFrequency, ContactDifference
In [3]:
%%time
frame_contacts = ContactMap(traj[0])
The built-in plotter requires one of the matplotlib color maps the set the color scheme. If you select a divergent color map (useful if you want to look at contact differences), then you should give the parameters vmin=-1, and vmax=1. Otherwise, you should use vmin=0 and vmax=1.
In [4]:
%%time
(fig, ax) = frame_contacts.residue_contacts.plot(cmap='seismic', vmin=-1, vmax=1)
plt.xlabel("Residue")
plt.ylabel("Residue")
The plotting function return the matplotlib Figure and Axes objects, which allow you to make more manipulations to them later. I'll show an example of that in the "Changing the defaults" section.
We can also plot the atom-atom contacts, although it takes a little time. The built-in plotting function is best if there are not many contacts (if the matrix is sparse). If there are lots of contacts, sometimes other approaches can plot more quickly. See an example in the "Changing the defaults" section.
In [5]:
%%time
frame_contacts.atom_contacts.plot(cmap='seismic', vmin=-1, vmax=1);
Out[5]:
You'll notice that you don't see many points here. As the warning states: That is because the points are smaller than a single pixel at this resolution. To fix that, output or save this figure as a vector image (like PDF), increase the figure's size (by using the figsize parameter) or dpi (by using the dpi parameter), or both.
In [6]:
# If you want Jupyter to output pdf, instead of png:
# from IPython.display import set_matplotlib_formats
# set_matplotlib_formats('pdf')
# frame_contacts.atom_contacts.plot(cmap='seismic', vmin=-1, vmax=1);
# If you want to save as a vector image:
# frame_contacts.atom_contacts.plot(cmap='seismic', vmin=-1, vmax=1);
# plt.savefig('atom_contacts.pdf')
# If you want to increase the figure size:
# frame_contacts.atom_contacts.plot(cmap='seismic', vmin=-1, vmax=1, figsize=(38, 38));
# If you want to increase the dpi (and make the image square):
# frame_contacts.atom_contacts.plot(cmap='seismic', vmin=-1, vmax=1, figsize=(6, 6), dpi=454);
In [7]:
%%time
trajectory_contacts = ContactFrequency(traj)
In [8]:
# if you want to save this for later analysis
trajectory_contacts.save_to_file("traj_contacts.p")
# then load with ContactFrequency.from_file("traj_contacts.p")
In [9]:
%%time
fig, ax = trajectory_contacts.residue_contacts.plot(cmap='seismic', vmin=-1, vmax=1);
In [10]:
%%time
diff = trajectory_contacts - frame_contacts
A contact that appears in trajectory, but not in the frame, will be at +1 and will be shown in red below. A contact that appears in the frame, but not the trajectory, will be at -1 and will be shown in blue below. The values are the difference in the frequencies (of course, for a single frame, the frequency is always 0 or 1).
In [11]:
%%time
diff.residue_contacts.plot(cmap='seismic', vmin=-1, vmax=1);
Out[11]:
You could have created the same object with:
diff = ContactDifference(trajectory_contact, frame_contacts)
but the simple notation using - is much more straightforward. However, note that ContactDifference makes a difference between the frequencies in the two objects, not the absolute count. Otherwise the trajectory would swamp the single frame, and there would be no blue in that picture!
First we look at the contacts that are much more important in the trajectory than the frame. Then we look at the contacts that are more important in the frame than the trajectory.
The .most_common() method gives a list of the contact pairs and the frequency, sorted by frequency. See also collections.Counter.most_common() in the standard Python collections module.
Here we do this with the ContactDifference we created, although it works the same for ContactFrequency and ContactMap (with the single-frame contact map, the ordering is a bit nonsensical, since every entry is either 0 or 1).
In [12]:
%%time
# residue contact more important in trajectory than in frame (near +1)
diff.residue_contacts.most_common()[:10]
Out[12]:
In [13]:
# residue contact more important in frame than in trajectory (near -1)
list(reversed(diff.residue_contacts.most_common()))[:10]
# alternate: diff.residue_contacts.most_common()[:-10:-1] # (thanks Sander!)
Out[13]:
First let's select a few residues from the topology. Note that GTP has residue ID 201 in the PDB sequence, even though it is only residue 166 (counting from 0) in the topology. This is because some of the protein was removed, and therefore the PDB is missing those residues. The topology only counts the residues that are actually present.
In [14]:
val81 = topology.residue(80)
asn116 = topology.residue(115)
gtp201 = topology.residue(166)
print(val81, asn116, gtp201)
We extended the standard .most_common() to take an optional argument. When the argument is given, it will filter the output to only include the ones where that argument is part of the contact. For example, the following gives the residues most commonly in contact with GTP.
In [15]:
for contact in trajectory_contacts.residue_contacts.most_common(gtp201):
if contact[1] > 0.1:
print(contact)
We can also find all the atoms, for all residue contacts, that are in contact with a given residue, and return that sorted by frequency.
In [16]:
diff.most_common_atoms_for_residue(gtp201)[:15]
Out[16]:
Finally, we can look at which atoms are most commonly in contact within a given residue contact pair.
In [17]:
trajectory_contacts.most_common_atoms_for_contact([val81, asn116])
Out[17]:
This sections covers several options that you can modify to make the contact maps faster, and to focus on what you're most interested in.
The first three options change which atoms are included as possible contacts. We call these query and haystack, and although they are conceptually equivalent, the algorithm is designed such that the query should have fewer atoms than the haystack.
Both of these options take a list of atom index numbers. These are most easily created using MDTraj's atom selection language.
In [18]:
# the default selection is
default_selection = topology.select("not water and symbol != 'H'")
print(len(default_selection))
MDTraj allows queries based on different numbering systems: resid and resSeq. The resid is the internally-used residue number, and starts from 0. On the other hand, resSeq is the residue number given in the PDB, which usually starts from 1 (and is the number we usually refer to in literature). More details can be found in the documentation on MDTraj's atom selection language.
In [19]:
switch1 = topology.select("resSeq 32 to 38 and symbol != 'H'")
switch2 = topology.select("resSeq 59 to 67 and symbol != 'H'")
gtp = topology.select("resname GTP and symbol != 'H'")
mg = topology.select("element Mg")
cations = topology.select("resname NA or resname MG")
sodium = topology.select("resname NA")
In [20]:
%%time
sw1_contacts = ContactFrequency(trajectory=traj, query=switch1)
In [21]:
sw1_contacts.residue_contacts.plot(cmap='seismic', vmin=-1, vmax=1)
Out[21]:
In [22]:
%%time
cations_switch1 = ContactFrequency(trajectory=traj, query=cations, haystack=switch1)
In [23]:
fig, ax = cations_switch1.residue_contacts.plot(cmap='seismic', vmin=-1, vmax=1)
Now we'll plot again, but we'll change the x and y axes so that we only see switch 1 (the haystack) along x and cations (the query) along y:
In [24]:
fig, ax = cations_switch1.residue_contacts.plot(cmap='seismic', vmin=-1, vmax=1)
ax.set_xlim(*cations_switch1.haystack_residue_range)
ax.set_ylim(*cations_switch1.query_residue_range);
Here, of course, the boxes are much larger, and are long rectangles instead of squares. The box represents the residue number (in the resid numbering system) that is to its left and under it. So the most significant contacts here are between residue 36 and the ion listed as residue 167. Let's see just how frequently that contact is made:
In [25]:
print(cations_switch1.residue_contacts.counter[frozenset([36, 167])])
So about half the time. Now, which residue/ion are these? Remember, these indices start at 0, even though the tradition in science (and the PDB) is to count from 1. Furthermore, the PDB residue numbers for the ions skip the section of the protein that has been removed. But we can easily obtain the relevant residues:
In [26]:
print(traj.topology.residue(36))
print(traj.topology.residue(167))
So this is a contact between the Glu37 and the magnesium ion (which is listed as residue 202 in the PDB).
By default, we ignore atoms from 2 residues on either side of the given residue (and in the same chain). This is easily changed. However, even when you say to ignore no neighbors, you still ignore the residue's interactions with itself.
Note: for non-protein contacts, the chain is often poorly defined. In this example, the GTP and the Mg are listed sequentially in residue order, and therefore they are considered "neighbors" and their contacts are ignored.
In [27]:
%%time
ignore_none = ContactFrequency(trajectory=traj, n_neighbors_ignored=0)
In [28]:
ignore_none.residue_contacts.plot(cmap='seismic', vmin=-1, vmax=1);
In [29]:
%%time
large_cutoff = ContactFrequency(trajectory=traj, cutoff=1.5)
The cost of the built-in plot function also depends strongly on the number of contacts that are made. It is designed to work well for sparse matrices; as the matrix gets less sparse, other approaches may be better. Here's an example:
In [30]:
%%time
large_cutoff.residue_contacts.plot(cmap='seismic', vmin=-1, vmax=1);
Out[30]:
In [31]:
%%time
import matplotlib
cmap = matplotlib.pyplot.get_cmap('seismic')
norm = matplotlib.colors.Normalize(vmin=-1, vmax=1)
plot = matplotlib.pyplot.pcolor(large_cutoff.residue_contacts.df, cmap='seismic', vmin=-1, vmax=1)
plot.cmap.set_under(cmap(norm(0)));
In this case, using the pandas.DataFrame representation (obtained using .df) is faster. On the other hand, try using this approach on the atom-atom picture at the top! That will take a while.
You'll notice that these may not be pixel-perfect copies. This is because the number of pixels doesn't evenly divide into the number of residues. You can improve this by increasing the resolution (dpi in matplotlib) or the figure size. However, in both versions you can see the overall structure quite clearly. In addition, the color bar is only shown in the built-in version.