IMP  2.0.0
The Integrative Modeling Platform
emagefit_dock.py
1 #!/usr/bin/env python
2 
3 import IMP
4 import IMP.algebra as alg
5 import IMP.atom as atom
6 import IMP.em2d
7 import IMP.em2d.imp_general.representation as representation
8 import IMP.em2d.imp_general.io as io
9 import IMP.em2d.csv_related as csv_related
10 import csv
11 import sys
12 import logging
13 import subprocess
14 import errno
15 log = logging.getLogger("emagefit_dock")
16 
17 
18 """
19  This script can be used for:
20  1) docking:
21  $IMP python docking_related.py --dock --hex hex_transforms.txt
22  --int internal_transforms.txt receptor.pdb ligand.pdb
23  2) Writing docking solutions
24  $IMP python docking_related.py --w 10 --int hex_transforms.txt
25  receptor.pdb ligand.pdb
26 
27  NOTE: The script requires the HEX docking program to perform the docking.
28 """
29 
30 
31 
32 def check_for_hexdock():
33  log.debug("Checking for HEXDOCK")
34  try:
35  subprocess.Popen(["hex"],stdin=subprocess.PIPE ,
36  stdout=subprocess.PIPE, stderr=subprocess.PIPE)
37  except OSError, e:
38  if e.errno == errno.ENOENT:
39  print "Docking requires the program HEXDOCK."
40  print "Make sure that the command hex is available in your path."
41  sys.exit()
42  else:
43  raise e
44 
45 
46 class HexDocking(object):
47 
48  def __init__(self, ):
49  ##### commands for HEX (the %s are filled when calling self.dock)
50  self.text_base = \
51  """
52  open_receptor %s
53  open_ligand %s
54  docking_correlation 1
55  docking_fft_device 2
56  docking_fft_type 5
57  max_docking_solutions 100000
58  receptor_range_angle 90
59  ligand_range_angle 90
60  docking_receptor_stepsize 10.0
61  docking_ligand_stepsize 10.0
62  docking_alpha_stepsize 6.0
63  docking_r12_range 60
64  docking_r12_step 2.0
65  docking_main_scan 20
66  activate_docking
67  save_transform %s
68  save_ligand %s
69  """
70 
71  def dock(self, fn_receptor, fn_ligand, fn_transforms,
72  fn_docked="docked.pdb",
73  write_hex_dump=False ):
74  """
75  Docks the ligand into the receptor using HEX.
76  @param fn_receptor PDB of the receptor
77  @param fn_ligand PDB of the ligand
78  @param fn_transformations Output fiole with the results, which
79  are the transformations of the ligand to dock it
80  on the receptor.
81  @param fn_docked PDB file of the first solution
82  @param write_hex_dump For debugging. If True, the output of HEX
83  is written to a file with .dump extension
84  """
85  log.info("===> Docking %s into %s with HEX",fn_ligand, fn_receptor)
86  hex_commands = self.text_base % (fn_receptor, fn_ligand,
87  fn_transforms, fn_docked)
88  pid = subprocess.Popen(["hex"],stdin=subprocess.PIPE ,
89  stdout=subprocess.PIPE, stderr=subprocess.PIPE)
90  # send the commands to hex with the PIPE file object for stdin
91  pid.stdin.write(hex_commands)
92  out, err = pid.communicate()
93  if err != "": # there is some error output
94  fn = "hexdock_errors.txt"
95  f = open(fn, "w")
96  f.write(err)
97  f.close()
98  raise IOError("Error running HEX. See file hexdoc_errors.txt")
99 
100  if write_hex_dump:
101  fn = "%s-%s.dump" % (fn_receptor, fn_ligand)
102  f = open(fn, "w")
103  f.write(out)
104  f.close()
105 
106 
107 def parse_hex_transform(row):
108  """
109  Returns an HEX transform with the conventions used by IMP.
110  @param row A list containing a splitted line from the HEX output file
111  NOTE: An HEX rotation is the inverse of the IMP convention
112  """
113  euler = map(float,row[8:11])
114  xyz = map(float, row[5:8])
115  euler = [-euler[0], -euler[1], -euler[2]]
116  R = alg.get_rotation_from_fixed_zyz(*euler)
117  R = R.get_inverse()
118  t = alg.Vector3D(*xyz)
119  return alg.Transformation3D(R, t)
120 
121 def read_hex_transforms(fn_transforms):
122  """
123  Parses a file of hex transforms. It returns the transformations as
124  alg.Transformation3D objects
125  @param fn_transforms
126  """
127  rows = csv_related.read_csv(fn_transforms, delimiter=" ")
128  return [parse_hex_transform(row) for row in rows]
129 
130 def get_xlinks_are_good(coords_rec, coords_lig, distances):
131  """
132  Compare distances between points in coords_rec and coords_lig
133  @param coords_rec A set of Vector3D objects containing coordinates from
134  the receptor
135  @param coords_lig Coordinates from the ligand
136  @param distances The set of maximum distances for each pair of
137  (coord_rec, coord_lig)
138 
139  @return True if all the distances between the pairs of coordinates are
140  below the thresholds contained in the argument "distances""
141  """
142  for vi,vj,d in zip(coords_rec, coords_lig, distances):
143  if alg.get_distance(vi, vj) > d:
144  return False
145  return True
146 
147 def filter_docking_results(h_receptor, h_ligand,
148  list_xlinks, fn_transforms,
149  fn_filtered, max_number=False):
150  """
151  Check if the set of transforms proposed by docking with HEX is
152  compatible with the distances between aminoacids from crosslinking
153 
154  @param h_receptor atom.Hierarchy for the receptor
155  @param h_ligand atom.Hierarchy for the ligand
156  @param list_xlinks - list of tuples with the format
157  (residue receptor, residue ligand, distance)
158  @param list_xlinks - list of Xlink class
159  @param fn_transforms File of transforms as given by HEX
160  @param fn_filtered Output file that will contain only the
161  transformations satisfying the cross-linking
162  restraints
163  """
164  log.info("Filtering results of docking in %s with links:", fn_transforms)
165  coords_rec = []
166  coords_lig = []
167  threshold_distances = []
168  for xl in list_xlinks:
169  log.info("%s", xl.show())
170  coords_rec.append(representation.get_residue_coordinates(
171  h_receptor, xl.first_chain, xl.first_residue))
172  coords_lig.append(representation.get_residue_coordinates(
173  h_ligand, xl.second_chain, xl.second_residue))
174  threshold_distances.append(xl.distance);
175  rows = csv_related.read_csv(fn_transforms, delimiter=" ",
176  max_number=max_number)
177  good = []
178  for r in rows:
179  T = parse_hex_transform(r)
180  transformed = [T.get_transformed(c) for c in coords_lig ]
181  if get_xlinks_are_good(coords_rec, transformed, threshold_distances):
182  good.append(r)
183 
184  # If there are not good transforms because all were filtered, the best
185  # bet is to keep all of them and trust the randomization procedure.
186  f_output = open(fn_filtered, "w")
187  w = csv.writer(f_output, delimiter=" ")
188  if len(good) ==0:
189  log.warning("No docking solution satisfies all the x-linking " \
190  "restraints. Keeping all of them")
191  w.writerows(rows)
192  else:
193  w.writerows(good)
194  f_output.close()
195 
196 
197 def get_internal_transform(Thex, rb_receptor, rb_ligand):
198  """
199  Get the internal transform (Ti) resulting from applying
200  the transformation Thex (from HEX) that docks the ligand.
201  The relationship betwen the reference frame of the ligand docked
202  (Tdock) and the internal transformation (Ti) respect to the receptor
203  (Trec) is Tdock = Trec * Tinternal
204 
205  @param Thex HEX transformation
206  @param rb_receptor Rigid body of the receptor
207  @param rb_ligand Rigid body of the ligand
208  @return The internal transformation
209  """
210  Tlig = rb_ligand.get_reference_frame().get_transformation_to()
211  Trec = rb_receptor.get_reference_frame().get_transformation_to()
212  Tdock = alg.compose(Thex, Tlig)
213  Ti = alg.compose(Trec.get_inverse(), Tdock)
214  return Ti
215 
216 def get_docked_reference_frames(hex_transforms, rb_ligand):
217  """
218  @param hex_transforms A list of HEX transformations
219  @param rb_ligand The rigid body of the ligand
220  @return The reference frames that the rigid body of the ligand has
221  when docked. The docked reference frame has transformation:
222  Tdock = Thex * Tinitial
223  """
224  Trb = rb_ligand.get_reference_frame().get_transformation_to()
225  docked_ref_frames = []
226  for Thex in hex_transforms:
227  Tdock = alg.compose(Thex, Trb)
228  docked_rf = alg.ReferenceFrame3D(Tdock)
229  docked_ref_frames.append(docked_rf)
230  return docked_ref_frames
231 
232 
233 if __name__ == "__main__":
234  parser = IMP.OptionParser(
235  usage="""%prog [options] fn_receptor fn_ligand
236 
237 fn_receptor and fn_ligand are the filenames of the PDB files for the
238 receptor and ligand, respectively.""",
239  description="Functions related to docking with HEX",
240  imp_module=IMP.em2d)
241  parser.add_option("--dock",
242  action="store_true", dest="dock", default=False,
243  help="Dock the ligand into the receptor")
244  parser.add_option("--hex", dest="fn_transforms", default=False,
245  help="File where to write the hex transforms.")
246  parser.add_option("--w", type=int, dest="write", default=False,
247  help="Write a number of docking transformations" )
248  parser.add_option("--int", dest="fn_internal_transforms", default=False,
249  help="Input/Output file of internal transformations")
250  parser.add_option("--log", dest="log", default=None,
251  help="File for logging ")
252 
253  opts, args = parser.parse_args()
254  if len(args) != 2:
255  parser.error("Wrong number of arguments")
256  fn_receptor, fn_ligand = args
257 
258  if(opts.log):
259  logging.basicConfig(filename=opts.log, filemode="w")
260  else:
261  logging.basicConfig(stream=sys.stdout)
262  logging.root.setLevel(logging.DEBUG)
263 
264  sel = atom.ATOMPDBSelector()
265  m = IMP.Model()
266  h_receptor = atom.read_pdb(fn_receptor, m, sel)
267  rb_receptor = atom.create_rigid_body(h_receptor)
268  h_ligand = atom.read_pdb(fn_ligand, m, sel)
269  rb_ligand = atom.create_rigid_body(h_ligand)
270  if opts.dock:
271  check_for_hexdock()
272  if not opts.fn_transforms or not opts.fn_internal_transforms:
273  raise IOError("Docking requires the --int and --hex arguments")
274  hex_docking = HexDocking()
275  hex_docking.dock(fn_receptor, fn_ligand, opts.fn_transforms)
276  # read the HEX file of solutions and get the internal transformations
277  # giving the relative orientation of the ligand respect to the receptor
278  Ts = read_hex_transforms(opts.fn_transforms)
279  Tis = [get_internal_transform(T, rb_receptor, rb_ligand) for T in Ts]
280  io.write_transforms(Tis, opts.fn_internal_transforms)
281  elif opts.write:
282  # To write the positions correctly, the script requires that the
283  # ligand file is the same that was used for the docking
284  Tinternal = io.read_transforms(opts.fn_internal_transforms)
285  max_number = min(opts.write, len(Tinternal))
286  Trec = rb_receptor.get_reference_frame().get_transformation_to()
287  for i in range(max_number):
288  Tdock = alg.compose(Trec, Tinternal[i])
289  ref = alg.ReferenceFrame3D(Tdock)
290  rb_ligand.set_reference_frame(ref)
291  atom.write_pdb(h_ligand,"docked-%03d.pdb" % i)