text: user5118: 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.
  
  user6322: 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.
  
original_text: 5118: <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>
   6322: <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>
  