7 import IMP.em2d.imp_general.representation
as representation
8 import IMP.em2d.imp_general.io
as io
15 log = logging.getLogger(
"emagefit_dock")
19 This script can be used for:
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
27 NOTE: The script requires the HEX docking program to perform the docking.
32 def check_for_hexdock():
33 log.debug(
"Checking for HEXDOCK")
35 subprocess.Popen([
"hex"],stdin=subprocess.PIPE ,
36 stdout=subprocess.PIPE, stderr=subprocess.PIPE)
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."
46 class HexDocking(object):
57 max_docking_solutions 100000
58 receptor_range_angle 90
60 docking_receptor_stepsize 10.0
61 docking_ligand_stepsize 10.0
62 docking_alpha_stepsize 6.0
71 def dock(self, fn_receptor, fn_ligand, fn_transforms,
72 fn_docked=
"docked.pdb",
73 write_hex_dump=
False ):
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
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
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)
91 pid.stdin.write(hex_commands)
92 out, err = pid.communicate()
94 fn =
"hexdock_errors.txt"
98 raise IOError(
"Error running HEX. See file hexdoc_errors.txt")
101 fn =
"%s-%s.dump" % (fn_receptor, fn_ligand)
107 def parse_hex_transform(row):
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
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)
118 t = alg.Vector3D(*xyz)
119 return alg.Transformation3D(R, t)
121 def read_hex_transforms(fn_transforms):
123 Parses a file of hex transforms. It returns the transformations as
124 alg.Transformation3D objects
127 rows = csv_related.read_csv(fn_transforms, delimiter=
" ")
128 return [parse_hex_transform(row)
for row
in rows]
130 def get_xlinks_are_good(coords_rec, coords_lig, distances):
132 Compare distances between points in coords_rec and coords_lig
133 @param coords_rec A set of Vector3D objects containing coordinates from
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)
139 @return True if all the distances between the pairs of coordinates are
140 below the thresholds contained in the argument "distances""
142 for vi,vj,d
in zip(coords_rec, coords_lig, distances):
143 if alg.get_distance(vi, vj) > d:
147 def filter_docking_results(h_receptor, h_ligand,
148 list_xlinks, fn_transforms,
149 fn_filtered, max_number=
False):
151 Check if the set of transforms proposed by docking with HEX is
152 compatible with the distances between aminoacids from crosslinking
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
164 log.info(
"Filtering results of docking in %s with links:", fn_transforms)
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)
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):
186 f_output = open(fn_filtered,
"w")
187 w = csv.writer(f_output, delimiter=
" ")
189 log.warning(
"No docking solution satisfies all the x-linking " \
190 "restraints. Keeping all of them")
197 def get_internal_transform(Thex, rb_receptor, rb_ligand):
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
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
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)
216 def get_docked_reference_frames(hex_transforms, rb_ligand):
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
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
233 if __name__ ==
"__main__":
235 usage=
"""%prog [options] fn_receptor fn_ligand
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",
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 ")
253 opts, args = parser.parse_args()
255 parser.error(
"Wrong number of arguments")
256 fn_receptor, fn_ligand = args
259 logging.basicConfig(filename=opts.log, filemode=
"w")
261 logging.basicConfig(stream=sys.stdout)
262 logging.root.setLevel(logging.DEBUG)
264 sel = atom.ATOMPDBSelector()
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)
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)
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)
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)