%load_ext autoreload
%autoreload 2

import molsysmt as msm
from molsysmt import pyunitwizard as puw
import numpy as np
import matplotlib.pyplot as plt

Get contacts

A contact map is a logic matrix where the element (i-th,j-th) is True if the distance between i-th and j-th is lower, or equal, than a certain threshold. The contact map is a common tool, simple but effective, used to represent structural motifs from the protein-protein interaface of the protein folding, for example. MolSysMT includes a method, based on molsysmt.distance(), to obtain contact maps: molsysmt.contact_map(). As such, molsysmt.contact_map inherits many input arguments from molsysmt.distance(). Lets see a couple of examples in this section.

Lets get the contact map obtained with the threshold 1.2 nm and the CA atoms from the dimer found in the system with pdb id 1TCD:

molecular_system = msm.convert('1TCD', 'molsysmt.MolSys')

CA_atoms = msm.select(molecular_system, selection='atom_name=="CA"')

contact_map = msm.structure.get_contacts(molecular_system, selection=CA_atoms, threshold='1.2 nm')

contact_map

array([[[ True,  True,  True, ..., False, False, False],
        [ True,  True,  True, ..., False, False, False],
        [ True,  True,  True, ..., False, False, False],
        ...,
        [False, False, False, ...,  True,  True,  True],
        [False, False, False, ...,  True,  True,  True],
        [False, False, False, ...,  True,  True,  True]]])

Our molecular system has a single frame. That’s why the shape of the output is the following:

contact_map.shape

(1, 497, 497)

Where 497 is the total number of CA atoms in our system:

print(msm.get(molecular_system, element='atom', selection="atom_name=='CA'", n_atoms=True))

497

Usually the contact map is represented as a plot. The structural motifs can be then easily identify:

plt.imshow(contact_map[0], cmap='Greys', origin='lower')
plt.show()

The contact map can also be computed among elements from different sets. Lets get the contact map with the same threshold but this time between the CA atoms from both chains:

CA_atoms_chain_0 = msm.get(molecular_system, element='atom', selection="atom_name=='CA' and chain_index==0", atom_index=True)
CA_atoms_chain_1 = msm.get(molecular_system, element='atom', selection="atom_name=='CA' and chain_index==1", atom_index=True)

contact_map = msm.structure.get_contacts(molecular_system, selection=CA_atoms_chain_0, selection_2=CA_atoms_chain_1,
                              threshold='1.2 nm')

plt.imshow(contact_map[0], cmap='Greys', origin='lower')
plt.show()

The molecular system is homodimeric, thereby the interface is symmetric as it can be seen in this last plot.

Finnally, molsysmt.contact_map() can also work with atoms groups. Instead of getting the contact map between CA atoms, lets do it now with the geometric centers of the residues:

atoms_in_residues_chain_0 = msm.get(molecular_system, element='group',
                                    selection="molecule_type=='protein' and chain_index==0", atom_index=True)
atoms_in_residues_chain_1 = msm.get(molecular_system, element='group',
                                    selection="molecule_type=='protein' and chain_index==1", atom_index=True)

contact_map = msm.structure.get_contacts(molecular_system,
                              selection=atoms_in_residues_chain_0, center_of_atoms=True,
                              selection_2=atoms_in_residues_chain_1, center_of_atoms_2=True,
                              threshold='1.2 nm')

contact_map.shape

(1, 248, 249)

plt.imshow(contact_map[0], cmap='Greys', origin='lower')
plt.show()

molsys = msm.convert(msm.systems['chicken villin HP35']['chicken_villin_HP35.h5msm'])

contacts = msm.structure.get_contacts(molsys, selection='atom_name=="CA"', threshold='5 angstroms',
                                      output_type='pairs', output_indices='atom')

contacts

[[[8, 25],
  [25, 44],
  [44, 55],
  [55, 67],
  [55, 94],
  [67, 82],
  [82, 94],
  [82, 136],
  [94, 114],
  [114, 136],
  [114, 182],
  [136, 146],
  [146, 162],
  [162, 182],
  [182, 189],
  [189, 206],
  [206, 220],
  [220, 244],
  [244, 255],
  [255, 265],
  [265, 285],
  [285, 295],
  [295, 309],
  [309, 327],
  [327, 342],
  [342, 361],
  [342, 407],
  [361, 385],
  [385, 407],
  [407, 424],
  [424, 441],
  [424, 474],
  [441, 455],
  [455, 474],
  [455, 518],
  [474, 496],
  [496, 518],
  [518, 533],
  [533, 555],
  [555, 562],
  [555, 581],
  [562, 581]]]