Back to index

python-biopython  1.60
ResidueDepth.py
Go to the documentation of this file.
00001 # Copyright (C) 2002, Thomas Hamelryck (thamelry@binf.ku.dk)
00002 # This code is part of the Biopython distribution and governed by its
00003 # license.  Please see the LICENSE file that should have been included
00004 # as part of this package.
00005 
00006 """Calculation of residue depth using command line tool MSMS.
00007 
00008 This module uses Michel Sanner's MSMS program for the surface calculation
00009 (specifically commands msms and pdb_to_xyzr). See:
00010 http://mgltools.scripps.edu/packages/MSMS
00011 
00012 Residue depth is the average distance of the atoms of a residue from 
00013 the solvent accessible surface.
00014 
00015 Residue Depth:
00016 
00017     >>> rd = ResidueDepth(model, pdb_file)
00018     >>> print rd[(chain_id, res_id)]
00019 
00020 Direct MSMS interface:
00021 
00022     Typical use:
00023 
00024         >>> surface = get_surface("1FAT.pdb")
00025 
00026     Surface is a Numeric array with all the surface 
00027     vertices.  
00028 
00029     Distance to surface:
00030 
00031         >>> dist = min_dist(coord, surface)
00032 
00033     where coord is the coord of an atom within the volume
00034     bound by the surface (ie. atom depth).
00035 
00036     To calculate the residue depth (average atom depth
00037     of the atoms in a residue):
00038 
00039         >>> rd = residue_depth(residue, surface)
00040 """
00041 
00042 import os
00043 import tempfile
00044 
00045 import numpy
00046 
00047 from Bio.PDB import Selection
00048 from Bio.PDB.AbstractPropertyMap import AbstractPropertyMap
00049 from Bio.PDB.Polypeptide import is_aa
00050 
00051 
00052 def _read_vertex_array(filename):
00053     """
00054     Read the vertex list into a Numeric array.
00055     """
00056     fp=open(filename, "r")
00057     vertex_list=[]
00058     for l in fp.readlines():
00059         sl=l.split()
00060         if not len(sl)==9:
00061             # skip header
00062             continue
00063         vl=map(float, sl[0:3])
00064         vertex_list.append(vl)
00065     fp.close()
00066     return numpy.array(vertex_list)
00067 
00068 def get_surface(pdb_file, PDB_TO_XYZR="pdb_to_xyzr", MSMS="msms"):
00069     """
00070     Return a Numeric array that represents 
00071     the vertex list of the molecular surface.
00072 
00073     PDB_TO_XYZR --- pdb_to_xyzr executable (arg. to os.system)
00074     MSMS --- msms executable (arg. to os.system)
00075     """
00076     # extract xyz and set radii
00077     xyz_tmp=tempfile.mktemp()
00078     PDB_TO_XYZR=PDB_TO_XYZR+" %s > %s"
00079     make_xyz=PDB_TO_XYZR % (pdb_file, xyz_tmp)
00080     os.system(make_xyz)
00081     assert os.path.isfile(xyz_tmp), \
00082         "Failed to generate XYZR file using command:\n%s" % make_xyz
00083     # make surface
00084     surface_tmp=tempfile.mktemp()
00085     MSMS=MSMS+" -probe_radius 1.5 -if %s -of %s > "+tempfile.mktemp()
00086     make_surface=MSMS % (xyz_tmp, surface_tmp)
00087     os.system(make_surface)
00088     surface_file=surface_tmp+".vert"
00089     assert os.path.isfile(surface_file), \
00090         "Failed to generate surface file using command:\n%s" % make_surface
00091     # read surface vertices from vertex file
00092     surface=_read_vertex_array(surface_file)
00093     # clean up tmp files
00094     # ...this is dangerous
00095     #os.system("rm "+xyz_tmp)
00096     #os.system("rm "+surface_tmp+".vert")
00097     #os.system("rm "+surface_tmp+".face")
00098     return surface
00099 
00100 
00101 def min_dist(coord, surface):
00102     """
00103     Return minimum distance between coord
00104     and surface.
00105     """
00106     d=surface-coord
00107     d2=numpy.sum(d*d, 1)
00108     return numpy.sqrt(min(d2))
00109 
00110 def residue_depth(residue, surface):
00111     """
00112     Return average distance to surface for all
00113     atoms in a residue, ie. the residue depth.
00114     """
00115     atom_list=residue.get_unpacked_list()
00116     length=len(atom_list)
00117     d=0
00118     for atom in atom_list:
00119         coord=atom.get_coord()
00120         d=d+min_dist(coord, surface)
00121     return d/length
00122 
00123 def ca_depth(residue, surface):
00124     if not residue.has_id("CA"):
00125         return None
00126     ca=residue["CA"]
00127     coord=ca.get_coord()
00128     return min_dist(coord, surface)
00129 
00130 class ResidueDepth(AbstractPropertyMap):
00131     """
00132     Calculate residue and CA depth for all residues.
00133     """
00134     def __init__(self, model, pdb_file):
00135         depth_dict={}
00136         depth_list=[]
00137         depth_keys=[]
00138         # get_residue
00139         residue_list=Selection.unfold_entities(model, 'R')
00140         # make surface from PDB file
00141         surface=get_surface(pdb_file)
00142         # calculate rdepth for each residue
00143         for residue in residue_list:
00144             if not is_aa(residue):
00145                 continue
00146             rd=residue_depth(residue, surface)
00147             ca_rd=ca_depth(residue, surface)
00148             # Get the key
00149             res_id=residue.get_id()
00150             chain_id=residue.get_parent().get_id()
00151             depth_dict[(chain_id, res_id)]=(rd, ca_rd)
00152             depth_list.append((residue, (rd, ca_rd)))
00153             depth_keys.append((chain_id, res_id))
00154             # Update xtra information
00155             residue.xtra['EXP_RD']=rd
00156             residue.xtra['EXP_RD_CA']=ca_rd
00157         AbstractPropertyMap.__init__(self, depth_dict, depth_keys, depth_list)
00158 
00159 
00160 if __name__=="__main__":
00161 
00162     import sys
00163     from Bio.PDB import PDBParser
00164 
00165     p=PDBParser()
00166     s=p.get_structure("X", sys.argv[1])
00167     model=s[0]
00168 
00169     rd=ResidueDepth(model, sys.argv[1])
00170 
00171 
00172     for item in rd:
00173         print item
00174