Interaction coverage of Protein-Ligand system

Recently, I read blog from iwantobipen and found a python package for drug design, and it is called Open Drug Discovery Toolkit (oddt). Based on RDkit and OpenBable, it has two amazing function:

• ligand-protein interaction finger print
• electronic shape similarity calculation

Here I want to show you about the function of interaction calculation. All interaction, even including halogen bond, pi-pi stacting and cation-pi intercation.

Preparation of protein and ligand

First of all, I have a 200 ns trjectory of a protein-ligand system in gromacs format. I dumped one frame as pdb format in every 20 ns, so that I have 11 frames including initial state.

Each pdb structure was processed as follow: 1. seperate protein and ligand into two file. 2. convert ligand into mol2 format(oddt can’t recognize ligand in pdb format…).

You can prepare any pdb file from RCSB PDB website. BTW, oddt does have some difficluty to deal any original pdb file from web… I would advise you should prepare protein by UCSF Chimera or Maestro after download structures. ;P

Import ODDT and set environment

from oddt import interactions
from oddt import toolkit

Then I packed some function into one function and a formatted output function.

def GetInteraction(pdb, lig):
    '''Only Accepte PDB format for protein and MOL2 for Ligand'''
    pdb = next(toolkit.readfile('pdb', pdb))
    lig = next(toolkit.readfile('mol2', lig))
    pdb.protein = True
    interac = []
    #Hydrogen Bond:
    hydrogen_bond = interactions.hbonds(pdb, lig, cutoff=4)
    interac.append(set(zip(hydrogen_bond[0]['resid'], hydrogen_bond[0]['resname'])))
    #Halogen Bond:
    halogen_bond = interactions.halogenbonds(pdb, lig)
    interac.append(set(zip(halogen_bond[0]['resid'], halogen_bond[0]['resname'])))
    #Pi-Pi Bond:
    pi_stack = interactions.pi_stacking(pdb, lig, cutoff=6)
    interac.append(set(zip(pi_stack[0]['resid'], pi_stack[0]['resname'])))
    #Hyrdophobic Interaion:
    hydrophobic = interactions.hydrophobic_contacts(pdb, lig)
    interac.append(set(zip(hydrophobic[0]['resid'], hydrophobic[0]['resname'])))
    
    return interac

def PrintInteration(frame, inp):
    t = frame * 20
    print('[+][+]Interactions for %d ns:' % t)
    print('      Hydrogen Bond:')
    print('      ', inp[0])
    print('      Halogen Bond:')
    print('      ', inp[1])
    print('      Pi-Pi Stacking:')
    print('      ', inp[2])
    print('      Hrdrophobic Interaction:')
    print('      ', inp[3])
    print('[+][+]......\n')

Calculate interaction between protein and ligand

Honestly, I use a GPCR sturcture here as an example, but I can’t tell you the ubpublished ligand. The GPCR family is one of most important drug targets at current stage. Another important target is kinease~

%%time

for i in range(0, 11):
    j = 1 + i
    print('[ ][ ]Processing NO.%d Frame' % j)
    pdb = "/home/dingkang/Documents/Project/SANT-1/SANT-MD/analysis/clean/0.%d.pdb" %j
    lig = "/home/dingkang/Documents/Project/SANT-1/SANT-MD/analysis/clean/0.%d.mol2" %j
    try:
        out = GetInteraction(pdb, lig)
        PrintInteration(i, out)
    except:
        print('[ ][ ]WARNNING: Process NO.%d Frame Failed!\n[+][+]...\n' % j)

OUTPUT:

[ ][ ]Processing NO.1 Frame
[+][+]Interactions for 0 ns:
      Hydrogen Bond:
       {(281, 'TRP'), (470, 'HIS')}
      Halogen Bond:
       {(394, 'TYR')}
      Pi-Pi Stacking:
       {(281, 'TRP'), (470, 'HIS'), (391, 'PHE')}
      Hrdrophobic Interaction:
       {(281, 'TRP'), (325, 'LEU'), (274, 'PHE'), (522, 'LEU'), (329, 'VAL'), (391, 'PHE'), (466, 'THR'), (521, 'ASN')}
[+][+]......

[ ][ ]Processing NO.2 Frame
[+][+]Interactions for 20 ns:
      Hydrogen Bond:
       {(281, 'TRP'), (470, 'HIS')}
      Halogen Bond:
       {(387, 'SER')}
      Pi-Pi Stacking:
       {(281, 'TRP'), (470, 'HIS'), (391, 'PHE')}
      Hrdrophobic Interaction:
       {(281, 'TRP'), (325, 'LEU'), (274, 'PHE'), (522, 'LEU'), (329, 'VAL'), (391, 'PHE')}
[+][+]......

[ ][ ]Processing NO.3 Frame
[+][+]Interactions for 40 ns:
      Hydrogen Bond:
       {(281, 'TRP'), (470, 'HIS')}
      Halogen Bond:
       set()
      Pi-Pi Stacking:
       {(281, 'TRP'), (470, 'HIS'), (391, 'PHE')}
      Hrdrophobic Interaction:
       {(463, 'VAL'), (281, 'TRP'), (325, 'LEU'), (274, 'PHE'), (522, 'LEU'), (329, 'VAL'), (391, 'PHE'), (466, 'THR')}
[+][+]......

[ ][ ]Processing NO.4 Frame
[+][+]Interactions for 60 ns:
      Hydrogen Bond:
       {(281, 'TRP'), (470, 'HIS')}
      Halogen Bond:
       {(394, 'TYR')}
      Pi-Pi Stacking:
       {(281, 'TRP'), (470, 'HIS'), (391, 'PHE')}
      Hrdrophobic Interaction:
       {(408, 'ILE'), (281, 'TRP'), (325, 'LEU'), (274, 'PHE'), (329, 'VAL'), (391, 'PHE')}
[+][+]......

[ ][ ]Processing NO.5 Frame
[+][+]Interactions for 80 ns:
      Hydrogen Bond:
       {(470, 'HIS')}
      Halogen Bond:
       {(518, 'GLU')}
      Pi-Pi Stacking:
       {(470, 'HIS'), (391, 'PHE')}
      Hrdrophobic Interaction:
       {(528, 'THR'), (325, 'LEU'), (274, 'PHE'), (522, 'LEU'), (329, 'VAL'), (391, 'PHE'), (521, 'ASN')}
[+][+]......

[ ][ ]Processing NO.6 Frame
[+][+]Interactions for 100 ns:
      Hydrogen Bond:
       {(470, 'HIS')}
      Halogen Bond:
       {(394, 'TYR')}
      Pi-Pi Stacking:
       {(281, 'TRP'), (470, 'HIS'), (391, 'PHE')}
      Hrdrophobic Interaction:
       {(463, 'VAL'), (281, 'TRP'), (329, 'VAL'), (466, 'THR'), (391, 'PHE')}
[+][+]......

[ ][ ]Processing NO.7 Frame
[+][+]Interactions for 120 ns:
      Hydrogen Bond:
       {(281, 'TRP'), (470, 'HIS')}
      Halogen Bond:
       {(518, 'GLU'), (394, 'TYR')}
      Pi-Pi Stacking:
       {(281, 'TRP'), (470, 'HIS'), (391, 'PHE')}
      Hrdrophobic Interaction:
       {(281, 'TRP'), (325, 'LEU'), (274, 'PHE'), (522, 'LEU'), (329, 'VAL'), (391, 'PHE'), (466, 'THR')}
[+][+]......

[ ][ ]Processing NO.8 Frame
[+][+]Interactions for 140 ns:
      Hydrogen Bond:
       {(281, 'TRP'), (470, 'HIS')}
      Halogen Bond:
       set()
      Pi-Pi Stacking:
       {(281, 'TRP'), (470, 'HIS'), (391, 'PHE')}
      Hrdrophobic Interaction:
       {(528, 'THR'), (463, 'VAL'), (325, 'LEU'), (274, 'PHE'), (522, 'LEU'), (329, 'VAL'), (391, 'PHE'), (466, 'THR')}
[+][+]......

[ ][ ]Processing NO.9 Frame
[+][+]Interactions for 160 ns:
      Hydrogen Bond:
       {(281, 'TRP'), (470, 'HIS')}
      Halogen Bond:
       set()
      Pi-Pi Stacking:
       {(281, 'TRP'), (470, 'HIS'), (274, 'PHE'), (391, 'PHE')}
      Hrdrophobic Interaction:
       {(325, 'LEU'), (463, 'VAL'), (329, 'VAL'), (274, 'PHE'), (391, 'PHE')}
[+][+]......

[ ][ ]Processing NO.10 Frame
[+][+]Interactions for 180 ns:
      Hydrogen Bond:
       {(281, 'TRP'), (470, 'HIS')}
      Halogen Bond:
       {(518, 'GLU')}
      Pi-Pi Stacking:
       {(281, 'TRP'), (470, 'HIS'), (391, 'PHE')}
      Hrdrophobic Interaction:
       {(281, 'TRP'), (325, 'LEU'), (274, 'PHE'), (329, 'VAL'), (391, 'PHE'), (466, 'THR')}
[+][+]......

[ ][ ]Processing NO.11 Frame
[+][+]Interactions for 200 ns:
      Hydrogen Bond:
       {(281, 'TRP'), (470, 'HIS')}
      Halogen Bond:
       {(518, 'GLU')}
      Pi-Pi Stacking:
       {(281, 'TRP'), (470, 'HIS'), (391, 'PHE')}
      Hrdrophobic Interaction:
       {(466, 'THR'), (281, 'TRP'), (329, 'VAL'), (274, 'PHE'), (391, 'PHE')}
[+][+]......

CPU times: user 1min 53s, sys: 36 ms, total: 1min 53s
Wall time: 1min 53s

That’s all.Then I collected these imfomation and recorded into excel to see the interation converage during this 200 ns trajectory. Next time I hope I can driectly generate an illustrater to show the interation converage.