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

<pre><code>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
</code></pre>

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

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

<pre><code>for atom in pymol.cmd.get_model('ligand expand 3').atom:
    print(dir(atom))
    break
</code></pre>

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

CC BY-SA 4.0
2020-04-09T07:37:42.307Z
2020-04-11T16:27:21.323Z
2020-04-11T16:27:21.323Z
Matteo Ferla
<p>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.</p>

<p>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).</p>

<p>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.</p>

CC BY-SA 4.0
2020-04-07T15:24:31.677Z
See a related post at Biostars for examples of how to do this with BioStructures.jl and LigPlot+: https://www.biostars.org/p/425276
jgreener
CC BY-SA 4.0
2020-04-09T02:11:50.773Z
Thanks. That helps.
Zac R.
<python><protein-structure><docking>
Extract autodocked protein-ligand connections programatically
# Extract autodocked protein-ligand connections programatically

(Asked by: username_0 on 2020-04-07)

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.

**Tags:** python, protein-structure, docking
> **Comments:**
> **username_1**: See a related post at Biostars for examples of how to do this with BioStructures.jl and LigPlot+: https://www.biostars.org/p/425276  
> **username_0**: Thanks. That helps.  

---

## 1 answers

### Answer by username_2 on 2020-04-09. Score: 4 | Accepted answer

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.


---

