Extract autodocked protein-ligand connections programatically

Question

I have 2 crystal structures, one for my protein and one for my ligand (I have several protein-ligand pairs). I am using AutoDock Vina to simulate docking, which returns another file of the ligand with updated 3D coordinates. Everything is in pdbqt format, but I can use OpenBabel to convert them if needed.

I am looking for a programmatic way (mainly python, but other languages work as well) to load both the protein and the docked ligand, then extract only the acting atoms which connect the 2 structures together (and the bonds between them).

If a direct answer is too complex, I would be happy if someone could suggest a tool/framework I could use to manipulate this sort of data easily so I can implement a solution on my own.

Answer 1

With pymol as a python module (conda install -c schrodinger pymol, not the GUI one), it is very easy.

import pymol2
with pymol2.PyMOL() as pymol:
     pymol.cmd.load('protein.pdbqt', 'protein') # the second is the name within pymol instance.
     pymol.cmd.load('ligand.sdf', 'ligand')
     pymol.cmd.save('neighbours.pdb','ligand expand 3') # get atoms 3 Å around

Done! However, if you are doing stuff in python in the first place you might want the last line as a string: pymol.cmd.get_pdbstr('ligand expand 3'). Changing the selection 'ligand expand 3' to 'byres ligand expand 3' will save the full residues.

Alternatively you can get the atom coordinates themselves etc. with it. To see what attributes Atom has you have you can do this:

for atom in pymol.cmd.get_model('ligand expand 3').atom:
    print(dir(atom))
    break

Which will give a large amount of attributes, such as atom.coord, which is a 3-element list.