/usr/share/pyshared/Bio/PDB/ResidueDepth.py is in python-biopython 1.58-1.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 | # Copyright (C) 2002, Thomas Hamelryck (thamelry@binf.ku.dk)
# This code is part of the Biopython distribution and governed by its
# license. Please see the LICENSE file that should have been included
# as part of this package.
"""Calculation of residue depth using command line tool MSMS.
This module uses Michel Sanner's MSMS program for the surface calculation
(specifically commands msms and pdb_to_xyzr). See:
http://mgltools.scripps.edu/packages/MSMS
Residue depth is the average distance of the atoms of a residue from
the solvent accessible surface.
Residue Depth:
>>> rd = ResidueDepth(model, pdb_file)
>>> print rd[(chain_id, res_id)]
Direct MSMS interface:
Typical use:
>>> surface = get_surface("1FAT.pdb")
Surface is a Numeric array with all the surface
vertices.
Distance to surface:
>>> dist = min_dist(coord, surface)
where coord is the coord of an atom within the volume
bound by the surface (ie. atom depth).
To calculate the residue depth (average atom depth
of the atoms in a residue):
>>> rd = residue_depth(residue, surface)
"""
import os
import tempfile
import numpy
from Bio.PDB import Selection
from Bio.PDB.AbstractPropertyMap import AbstractPropertyMap
from Bio.PDB.Polypeptide import is_aa
def _read_vertex_array(filename):
"""
Read the vertex list into a Numeric array.
"""
fp=open(filename, "r")
vertex_list=[]
for l in fp.readlines():
sl=l.split()
if not len(sl)==9:
# skip header
continue
vl=map(float, sl[0:3])
vertex_list.append(vl)
fp.close()
return numpy.array(vertex_list)
def get_surface(pdb_file, PDB_TO_XYZR="pdb_to_xyzr", MSMS="msms"):
"""
Return a Numeric array that represents
the vertex list of the molecular surface.
PDB_TO_XYZR --- pdb_to_xyzr executable (arg. to os.system)
MSMS --- msms executable (arg. to os.system)
"""
# extract xyz and set radii
xyz_tmp=tempfile.mktemp()
PDB_TO_XYZR=PDB_TO_XYZR+" %s > %s"
make_xyz=PDB_TO_XYZR % (pdb_file, xyz_tmp)
os.system(make_xyz)
assert os.path.isfile(xyz_tmp), \
"Failed to generate XYZR file using command:\n%s" % make_xyz
# make surface
surface_tmp=tempfile.mktemp()
MSMS=MSMS+" -probe_radius 1.5 -if %s -of %s > "+tempfile.mktemp()
make_surface=MSMS % (xyz_tmp, surface_tmp)
os.system(make_surface)
surface_file=surface_tmp+".vert"
assert os.path.isfile(surface_file), \
"Failed to generate surface file using command:\n%s" % make_surface
# read surface vertices from vertex file
surface=_read_vertex_array(surface_file)
# clean up tmp files
# ...this is dangerous
#os.system("rm "+xyz_tmp)
#os.system("rm "+surface_tmp+".vert")
#os.system("rm "+surface_tmp+".face")
return surface
def min_dist(coord, surface):
"""
Return minimum distance between coord
and surface.
"""
d=surface-coord
d2=numpy.sum(d*d, 1)
return numpy.sqrt(min(d2))
def residue_depth(residue, surface):
"""
Return average distance to surface for all
atoms in a residue, ie. the residue depth.
"""
atom_list=residue.get_unpacked_list()
length=len(atom_list)
d=0
for atom in atom_list:
coord=atom.get_coord()
d=d+min_dist(coord, surface)
return d/length
def ca_depth(residue, surface):
if not residue.has_id("CA"):
return None
ca=residue["CA"]
coord=ca.get_coord()
return min_dist(coord, surface)
class ResidueDepth(AbstractPropertyMap):
"""
Calculate residue and CA depth for all residues.
"""
def __init__(self, model, pdb_file):
depth_dict={}
depth_list=[]
depth_keys=[]
# get_residue
residue_list=Selection.unfold_entities(model, 'R')
# make surface from PDB file
surface=get_surface(pdb_file)
# calculate rdepth for each residue
for residue in residue_list:
if not is_aa(residue):
continue
rd=residue_depth(residue, surface)
ca_rd=ca_depth(residue, surface)
# Get the key
res_id=residue.get_id()
chain_id=residue.get_parent().get_id()
depth_dict[(chain_id, res_id)]=(rd, ca_rd)
depth_list.append((residue, (rd, ca_rd)))
depth_keys.append((chain_id, res_id))
# Update xtra information
residue.xtra['EXP_RD']=rd
residue.xtra['EXP_RD_CA']=ca_rd
AbstractPropertyMap.__init__(self, depth_dict, depth_keys, depth_list)
if __name__=="__main__":
import sys
from Bio.PDB import PDBParser
p=PDBParser()
s=p.get_structure("X", sys.argv[1])
model=s[0]
rd=ResidueDepth(model, sys.argv[1])
for item in rd:
print item
|