4. Visualisation

This notebook showcases different ways of visualizing lig-prot and prot-prot interactions, either with atomistic details or simply at the residue level.

This is a work in progress…

import MDAnalysis as mda
import prolif as plf
# load topology
u = mda.Universe(plf.datafiles.TOP)
lig = u.select_atoms("resname LIG")
prot = u.select_atoms("protein")

# create RDKit-like molecules for visualisation
lmol = plf.Molecule.from_mda(lig)
pmol = plf.Molecule.from_mda(prot)

# get lig-prot interactions with atom info
fp = plf.Fingerprint(["HBDonor", "HBAcceptor", "Cationic", "PiStacking"])
fp.run(u.trajectory, lig, prot, return_atoms=True)
df = fp.to_dataframe()
df.T

		Frame	0
ligand	protein	interaction
LIG1.G	ASP129.A	HBDonor	[52, 10]
		Cationic	[13, 10]
	VAL201.A	HBAcceptor	[28, 6]
	SER212.A	HBDonor	[44, 2]
	PHE331.B	PiStacking	[3, 9]

4.1. py3Dmol (3Dmol.js)

With py3Dmol we can easily display the interactions.

For interactions involving a ring (pi-cation, pi-stacking…etc.) ProLIF returns the index of one of the ring atoms, but for visualisation having the centroid of the ring looks nicer. We’ll start by writing a function to find the centroid, given the index of one of the ring atoms.

from rdkit import Chem
from rdkit import Geometry

def get_ring_centroid(mol, index):
    # find ring using the atom index
    Chem.SanitizeMol(mol, Chem.SanitizeFlags.SANITIZE_SETAROMATICITY)
    ri = mol.GetRingInfo()
    for r in ri.AtomRings():
        if index in r:
            break
    else:
        raise ValueError("No ring containing this atom index was found in the given molecule")
    # get centroid
    coords = mol.xyz[list(r)]
    ctd = plf.utils.get_centroid(coords)
    return Geometry.Point3D(*ctd)

Finally, the actual visualisation code. The API of py3Dmol is exactly the same as the GLViewer class of 3Dmol.js, for which the documentation can be found here.

import py3Dmol

colors = {
    "HBAcceptor": "blue",
    "HBDonor": "red",
    "Cationic": "green",
    "PiStacking": "purple",
}

hover_func = """
function(atom,viewer) {
    if(!atom.label) {
        atom.label = viewer.addLabel(atom.interaction,
            {position: atom, backgroundColor: 'mintcream', fontColor:'black'});
    }
}"""
unhover_func = """
function(atom,viewer) {
    if(atom.label) {
        viewer.removeLabel(atom.label);
        delete atom.label;
    }
}"""

v = py3Dmol.view(650, 600)
v.removeAllModels()

models = {}
mid = -1
for i, row in df.T.iterrows():
    lresid, presid, interaction = i
    lindex, pindex = row[0]
    lres = lmol[lresid]
    pres = pmol[presid]
    # set model ids for reusing later
    if lresid not in models.keys():
        mid += 1
        v.addModel(Chem.MolToMolBlock(lres), "sdf")
        model = v.getModel()
        model.setStyle({}, {"stick": {"colorscheme": "cyanCarbon"}})
        models[lresid] = mid
    if presid not in models.keys():
        mid += 1
        v.addModel(Chem.MolToMolBlock(pres), "sdf")
        model = v.getModel()
        model.setStyle({}, {"stick": {}})
        models[presid] = mid
    # get coordinates for both points of the interaction
    if interaction in ["PiStacking", "EdgeToFace", "FaceToFace", "PiCation"]:
        p1 = get_ring_centroid(lres, lindex)
    else:
        p1 = lres.GetConformer().GetAtomPosition(lindex)
    if interaction in ["PiStacking", "EdgeToFace", "FaceToFace", "CationPi"]:
        p2 = get_ring_centroid(pres, pindex)
    else:
        p2 = pres.GetConformer().GetAtomPosition(pindex)
    # add interaction line
    v.addCylinder({"start": dict(x=p1.x, y=p1.y, z=p1.z),
                   "end":   dict(x=p2.x, y=p2.y, z=p2.z),
                   "color": colors[interaction],
                   "radius": .15,
                   "dashed": True,
                   "fromCap": 1,
                   "toCap": 1,
                  })
    # add label when hovering the middle of the dashed line by adding a dummy atom
    c = Geometry.Point3D(*plf.utils.get_centroid([p1, p2]))
    mid = models[lresid]
    model = v.getModel(mid)
    model.addAtoms([{"elem": 'Z',
                     "x": c.x, "y": c.y, "z": c.z,
                     "interaction": interaction}])
    model.setStyle({"interaction": interaction}, {"clicksphere": {"radius": .5}})
    model.setHoverable(
        {"interaction": interaction}, True,
        hover_func, unhover_func)

v.zoomTo()
v.show()

You appear to be running in JupyterLab (or JavaScript failed to load for some other reason). You need to install the 3dmol extension:
jupyter labextension install jupyterlab_3dmol

4.2. RDKit

This script is given as a starting point for visualising the interactions with RDKit. It’s definitely not ideal for 3D visualisation as the drawing ends up being crowded and quite unreadable very quickly. I’m sure there’s a better way to do this but it will do as a code snippet.

from rdkit.Chem import Draw
from IPython import display

colors = {
    "HBAcceptor": (0, 1, 1),
    "HBDonor": (1, .7,.3),
    "Cationic": (0, 1, 0),
    "PiStacking": (.5, 0, .5),
}

d = Draw.MolDraw2DSVG(650, 600)
opts = Draw.MolDrawOptions()
opts.fixedBondLength = 25
d.SetDrawOptions(opts)

displayed = []
for i, row in df.T.iterrows():
    lresid, presid, interaction = i
    lindex, pindex = row[0]
    lres = lmol[lresid]
    pres = pmol[presid]
    if lresid not in displayed:
        displayed.append(lresid)
        d.DrawMolecule(lres)
    if presid not in displayed:
        displayed.append(presid)
        d.DrawMolecule(pres)
    if interaction in ["PiStacking", "EdgeToFace", "FaceToFace", "PiCation"]:
        p1 = get_ring_centroid(lres, lindex)
    else:
        p1 = lres.GetConformer().GetAtomPosition(lindex)
    if interaction in ["PiStacking", "EdgeToFace", "FaceToFace", "CationPi"]:
        p2 = get_ring_centroid(pres, pindex)
    else:
        p2 = pres.GetConformer().GetAtomPosition(pindex)
    p1 = Geometry.Point2D(p1)
    p2 = Geometry.Point2D(p2)
    d.DrawWavyLine(p1, p2, colors[interaction], colors[interaction])

# display
d.FinishDrawing()
display.SVG(d.GetDrawingText())