Today’s blog post is about using PLIP to extract information about interactions between a protein and ligand in a bound complex, using data from PDBbind. The blog post will cover how to combine the protein pdb file and the ligand mol2 file into a pdb file, and how to use PLIP in a high-throughput manner with python.

In order for PLIP to consider the ligand as one molecule interacting with the protein, we need to modify the mol2 file of the ligand. The 8th column of the atom portion of a mol2 file (the portion starts with @<TRIPOS>ATOM) includes the ID of the ligand that the atom belongs to. Most often all the atoms have the same ligand ID, but for peptides for instance, the atoms have the ID of the residue they’re part of. The following code snippet will make the required changes:

ligand_file = 'data/5oxm/5oxm_ligand.mol2'

with open(ligand_file, 'r') as f:
    ligand_lines = f.readlines()

mod = False
for i in range(len(ligand_lines)):
    line = ligand_lines[i]
    if line == '@<TRIPOS>BOND\n':
        mod = False
        
    if mod:
        ligand_lines[i] = line[:59] + 'ISK     ' + line[67:]
        
    if line == '@<TRIPOS>ATOM\n':
        mod = True

with open('data/5oxm/5oxm_ligand_mod.mol2', 'w') as g:
    for j in ligand_lines:
        g.write(j)

Next, we have to combine the modified mol2 file and the pdb file of the protein into a pdb file of the complex. We do this using the pymol package in python:

from pymol import cmd

ligand_file = 'data/5oxm/5oxm_ligand_mod.mol2'
protein_file = 'data/5oxm/5oxm_protein.pdb'
complex_file = 'data/5oxm/5oxm_complex.pdb'

cmd.reinitialize()
cmd.load(protein_file, "prot")
cmd.load(ligand_file, "ligand")
cmd.create("complex", "ligand, prot")
cmd.save(complex_file, "complex")

Finally, we can use PLIP in python to get the interactions. Notice that we extract the output from PLIP that has ISK as longname, that’s the renamed ligand ID in the mol2 files.

from plip.structure.preparation import PDBComplex

mol = PDBComplex()
mol.load_pdb(complex_file)
mol.analyze()

longnames = [x.longname for x in mol.ligands]
bsids = [":".join([x.hetid, x.chain, str(x.position)]) for x in mol.ligands]

indices = [j for j,x in enumerate(longnames) if x == 'ISK']

for idx in indices:
    bsid = bsids[idx]
    interactions = mol.interaction_sets[bsid]
    print(interactions.all_itypes)