[BioPython] protein-ligand interactions

[Peter Cock]

On Fri, Mar 20, 2009 at 12:18 PM, mitlox <mitlox at op.pl> wrote:
> Thank you very much for your code, it works and the output is exactly for
> what I was looking for.
>
> I try to get a structureCA object to write out the results in a PDB file
> (outCA.pdb) like this:
> ATOM   5275  CA  ILE A 881      17.242  57.141  22.062  1.00 38.49
> C
> ATOM   5283  CA  VAL A 882      16.292  57.880  25.678  1.00 38.90
> C ....
>
> Unfortunately I get this error with ... Here is the code:
> ...
> structureCA = []
> ...
> io=Bio.PDB.PDBIO()
> io.set_structure(structureCA)
> io.save('outCA.pdb')

Your structureCA object is just a python list, containing Residue objects.
Instead you need to create a new object with the partial chain - which
can be done by creating structure, model and chain objects manually.

However, I suggest you re-read pages 5 and 6 of the Bio.PDB
documentation for the recommend approach:
http://biopython.org/DIST/docs/cookbook/biopdb_faq.pdf
In your case, you'll want to write your own selection class using the
residue distance to the ligand.  I recognise this might seem rather
complicated for a python novice as you have to create your own
class - so here is my solution:

import Bio.PDB
import numpy

pdb_code = "1E8W"
pdb_filename = "1E8W.pdb"

structure = Bio.PDB.PDBParser().get_structure(pdb_code, pdb_filename)
model = structure[0]
chainA = model["A"]

def residue_dist_to_ligand(protein_residue, ligand_residue) :
    """Returns distance from the protein C-alpha to the closest ligand atom."""
    distances = []
    for atom in ligand_residue :
        diff_vector  = protein_residue["CA"].coord - atom.coord
        distances.append(numpy.sqrt(numpy.sum(diff_vector * diff_vector)))
    return min(distances)

class NearLigandSelect(Bio.PDB.Select):
    def __init__(self, distance_threshold, ligand_residue) :
        self.threshold = distance_threshold
        self.ligand_res = ligand_residue

    def accept_residue(self, residue):
        if residue == self.ligand_res :
            return True #change this to False if you don't want the ligand
        else :
            dist = residue_dist_to_ligand(residue, self.ligand_res)
            return dist < self.threshold

io=Bio.PDB.PDBIO()
io.set_structure(structure)
#From looking at the PDB file, ligand is last residue in chain A
ligand_res = chainA.child_list[-1]
#Going to use a distance theshold of 4A
io.save("near_ligand.pdb", NearLigandSelect(4, ligand_res))
print "Done"

Peter