2. How-to

This notebook serves as a practical guide to common questions users might have.

[1]:
import MDAnalysis as mda
import prolif as plf
import numpy as np
[2]:
u = mda.Universe(plf.datafiles.TOP, plf.datafiles.TRAJ)
lig = u.select_atoms("resname LIG")
prot = u.select_atoms("protein")
lmol = plf.Molecule.from_mda(lig)
pmol = plf.Molecule.from_mda(prot)

2.1. 1. Changing the parameters for an interaction

You can list all the available interactions as follow:

[3]:
plf.Fingerprint.list_available(show_hidden=True)
[3]:
['Interaction',
 '_Distance',
 'Hydrophobic',
 '_BaseHBond',
 'HBDonor',
 'HBAcceptor',
 '_BaseXBond',
 'XBAcceptor',
 'XBDonor',
 '_BaseIonic',
 'Cationic',
 'Anionic',
 '_BaseCationPi',
 'PiCation',
 'CationPi',
 'PiStacking',
 'FaceToFace',
 'EdgeToFace',
 '_BaseMetallic',
 'MetalDonor',
 'MetalAcceptor']

In this example, we’ll redefine the hydrophobic interaction with a shorter distance.

You have the choice between overwriting the original hydrophobic interaction with the new one, or giving it an original name.

Let’s start with a test case: with the default parameters, Y109 is interacting with our ligand.

[4]:
fp = plf.Fingerprint()
fp.hydrophobic(lmol, pmol["TYR109.A"])
[4]:
True

2.1.1. 1.a Overwriting the original interaction

You have to define a class that inherits one of the classes listed in the prolif.interactions module.

[5]:
class Hydrophobic(plf.interactions.Hydrophobic):
    def __init__(self):
        super().__init__(distance=4.0)
/home/docs/checkouts/readthedocs.org/user_builds/prolif/conda/v0.3.0/lib/python3.9/site-packages/prolif/interactions.py:55: UserWarning: The 'Hydrophobic' interaction has been superseded by a new class with id 0x562970ee4770
  warnings.warn(f"The {name!r} interaction has been superseded by a "
[6]:
fp = plf.Fingerprint()
fp.hydrophobic(lmol, pmol["TYR109.A"])
[6]:
False

The interaction is not detected anymore. You can reset to the default interaction like so:

[7]:
class Hydrophobic(plf.interactions.Hydrophobic):
    pass

fp = plf.Fingerprint()
fp.hydrophobic(lmol, pmol["TYR109.A"])
/home/docs/checkouts/readthedocs.org/user_builds/prolif/conda/v0.3.0/lib/python3.9/site-packages/prolif/interactions.py:55: UserWarning: The 'Hydrophobic' interaction has been superseded by a new class with id 0x5629761db2f0
  warnings.warn(f"The {name!r} interaction has been superseded by a "
[7]:
True

2.1.2. 1.b Reparameterizing an interaction with another name

The steps are identical to above, just give the class a different name:

[8]:
class CustomHydrophobic(plf.interactions.Hydrophobic):
    def __init__(self):
        super().__init__(distance=4.0)

fp = plf.Fingerprint()
fp.hydrophobic(lmol, pmol["TYR109.A"])
[8]:
True
[9]:
fp.customhydrophobic(lmol, pmol["TYR109.A"])
[9]:
False
[10]:
fp = plf.Fingerprint(["Hydrophobic", "CustomHydrophobic"])
fp.bitvector(lmol, pmol["TYR109.A"])
[10]:
array([ True, False])

2.2. 2. Writing your own interaction

Before you dive into this section, make sure that there isn’t already an interaction that could just be reparameterized to do what you want!

For this, the best is to check the section of the documentation corresponding to the prolif.interactions module. There are some generic interactions, like the _Distance class, if you just need to define two chemical moieties within a certain distance. Both the Hydrophobic, Ionic, and Metallic interactions inherit from this class!

With that being said, there are a few rules that you must respect when writing your own interaction:

  • Inherit the ProLIF Interaction class

This class is located in prolif.interactions.Interaction. If for any reason you must inherit from another class, you can also define the prolif.interactions._InteractionMeta as a metaclass.

  • Naming convention

Your class name must not start with _ or be named Interaction. For non-symmetrical interactions, like hydrogen bonds or salt-bridges, the convention used here is to named the class after the function of the ligand. For example, the class HBDonor detects if a ligand acts as a hydrogen bond donor, and the class Cationic detects if a ligand acts as a cation.

  • Define a ``detect`` method

This method takes exactly two positional arguments (and as many named arguments as you need): a ligand Residue or Molecule and a protein Residue or Molecule (in this order).

  • Return value(s) for the ``detect`` method

There are two possibilities here, depending on whether or not you want to access the indices of atoms responsible for the interaction. If you don’t need this information, just return True if the interaction is detected, False otherwise. If you need to access atomic indices, you must return the following items in this order:

  • True or False for the detection of the interaction

  • The index of the ligand atom, or None if not detected

  • The index of the protein atom, or None if not detected

[11]:
from scipy.spatial import distance_matrix

# without atom indices
class CloseContact(plf.interactions.Interaction):
    def detect(self, res1, res2, threshold=2.0):
        dist_matrix = distance_matrix(res1.xyz, res2.xyz)
        if (dist_matrix <= threshold).any():
            return True
        return False

fp = plf.Fingerprint()
fp.closecontact(lmol, pmol["ASP129.A"])
[11]:
True
[12]:
# with atom indices
class CloseContact(plf.interactions.Interaction):
    def detect(self, res1, res2, threshold=2.0):
        dist_matrix = distance_matrix(res1.xyz, res2.xyz)
        mask = dist_matrix <= threshold
        if mask.any():
            res1_i, res2_i = np.where(mask)
            # return the first solution
            return True, res1_i[0], res2_i[0]
        return False, None, None

fp = plf.Fingerprint()
fp.closecontact(lmol, pmol["ASP129.A"])
/home/docs/checkouts/readthedocs.org/user_builds/prolif/conda/v0.3.0/lib/python3.9/site-packages/prolif/interactions.py:55: UserWarning: The 'CloseContact' interaction has been superseded by a new class with id 0x5629741118b0
  warnings.warn(f"The {name!r} interaction has been superseded by a "
[12]:
True

By default, the fingerprint will modify all interaction classes to only return the boolean value. To get the atom indices when using your custom class, you must choose one of the following options:

  • Use the __wrapped__ argument when calling the class as a fingerprint method:

[13]:
fp.closecontact.__wrapped__(lmol, pmol["ASP129.A"])
[13]:
(True, 52, 10)
  • Use the bitvector_atoms method instead of bitvector:

[14]:
fp = plf.Fingerprint(["CloseContact"])
bv, indices = fp.bitvector_atoms(lmol, pmol["ASP129.A"])
bv, indices
[14]:
(array([ True]), [[52, 10]])
  • Use the return_atoms=True argument when calling the run method:

[15]:
fp.run(u.trajectory[0:1], lig, prot, return_atoms=True)
fp.ifp

[15]:
[{'Frame': 0,
  (LIG1.G, MET337.B): [[None, None]],
  (LIG1.G, GLN41.A): [[None, None]],
  (LIG1.G, SER212.A): [[None, None]],
  (LIG1.G, SER106.A): [[None, None]],
  (LIG1.G, ALA349.B): [[None, None]],
  (LIG1.G, PHE331.B): [[None, None]],
  (LIG1.G, VAL210.A): [[None, None]],
  (LIG1.G, THR203.A): [[None, None]],
  (LIG1.G, TYR109.A): [[None, None]],
  (LIG1.G, ASP352.B): [[None, None]],
  (LIG1.G, VAL201.A): [[28, 6]],
  (LIG1.G, TYR208.A): [[None, None]],
  (LIG1.G, TRP356.B): [[None, None]],
  (LIG1.G, VAL102.A): [[None, None]],
  (LIG1.G, ILE333.B): [[None, None]],
  (LIG1.G, ILE137.A): [[None, None]],
  (LIG1.G, TRP125.A): [[None, None]],
  (LIG1.G, SER334.B): [[None, None]],
  (LIG1.G, VAL214.A): [[None, None]],
  (LIG1.G, THR110.A): [[None, None]],
  (LIG1.G, PHE330.B): [[None, None]],
  (LIG1.G, TRP327.B): [[None, None]],
  (LIG1.G, THR209.A): [[None, None]],
  (LIG1.G, TYR38.A): [[None, None]],
  (LIG1.G, THR134.A): [[None, None]],
  (LIG1.G, ASP129.A): [[52, 10]],
  (LIG1.G, PHE353.B): [[None, None]],
  (LIG1.G, THR213.A): [[None, None]],
  (LIG1.G, GLY215.A): [[None, None]],
  (LIG1.G, ILE130.A): [[None, None]],
  (LIG1.G, ILE350.B): [[None, None]],
  (LIG1.G, PHE354.B): [[None, None]],
  (LIG1.G, ALA216.A): [[None, None]],
  (LIG1.G, TYR359.B): [[None, None]],
  (LIG1.G, TYR40.A): [[None, None]],
  (LIG1.G, LEU126.A): [[None, None]],
  (LIG1.G, PRO338.B): [[None, None]],
  (LIG1.G, VAL200.A): [[None, None]],
  (LIG1.G, CYS199.A): [[None, None]],
  (LIG1.G, GLY358.B): [[None, None]],
  (LIG1.G, TRP115.A): [[None, None]],
  (LIG1.G, LEU335.B): [[None, None]],
  (LIG1.G, GLU198.A): [[None, None]],
  (LIG1.G, LEU348.B): [[None, None]],
  (LIG1.G, CYS133.A): [[None, None]],
  (LIG1.G, PHE217.A): [[None, None]],
  (LIG1.G, THR131.A): [[None, None]],
  (LIG1.G, PHE351.B): [[None, None]],
  (LIG1.G, TYR211.A): [[None, None]],
  (LIG1.G, ILE180.A): [[None, None]],
  (LIG1.G, ASN202.A): [[None, None]],
  (LIG1.G, THR355.B): [[None, None]]}]
  • Directly use your class:

[16]:
cc = CloseContact()
cc.detect(lmol, pmol["ASP129.A"])
[16]:
(True, 52, 10)

2.3. 3. Ignoring the protein backbone when computing interactions

[17]:
# not implemented yet