IMP logo
IMP Reference Guide  2.11.0
The Integrative Modeling Platform
buildxlinks.py
1 """@namespace IMP.EMageFit.buildxlinks
2  Utility functions to handle cross links.
3 """
4 
5 import IMP
6 import IMP.atom as atom
7 import IMP.algebra as alg
8 import IMP.core as core
9 import sys
10 import logging
11 from IMP.EMageFit.argminmax import keymax
12 log = logging.getLogger("buildxlinks")
13 
14 import networkx as nx
15 
16 try:
17  set = set
18 except NameError:
19  from sets import Set as set
20 
21 
22 class Xlink:
23 
24  """
25  Class defining a cross-link
26  """
27 
28  def __init__(self, id1, chain1, residue1,
29  id2, chain2, residue2,
30  distance):
31  """
32  Initialize the class
33  @param[in] id1 Id of the first component of the cross-link
34  @param[in] chain1 Chain of the first component
35  @param[in] residue1 Residue cross-linked in the first component
36  @param[in] id2 Id of the second component of the cross-link
37  @param[in] chain2 chain1 Chain of the second component
38  @param[in] residue2 Residue cross-linked in the second component
39  @param[in] distance Maximum distance
40  """
41  self.first_id = id1
42  self.first_chain = chain1
43  self.first_residue = residue1
44  self.second_id = id2
45  self.second_chain = chain2
46  self.second_residue = residue2
47  self.distance = distance
48 
49  def clone(self):
50  xl = Xlink(self.first_id, self.first_chain, self.first_residue,
51  self.second_id, self.second_chain, self.second_residue,
52  self.distance)
53  return xl
54 
55  def __eq__(self, other):
56  if self.first_id != other.first_id or \
57  self.second_id != other.second_id or \
58  self.first_chain != other.first_chain or \
59  self.second_chain != other.second_chain or \
60  self.first_residue != other.first_residue or \
61  self.second_residue != other.second_residue or \
62  abs(other.distance - self.distance) > 0.01:
63  return False
64  return True
65 
66  def swap(self):
67  """
68  swaps the order of the residues in the restraint
69  """
70  self.first_id, self.second_id = self.second_id, self.first_id
71  self.first_residue, self.second_residue = \
72  self.second_residue, self.first_residue
73  self.first_chain, self.second_chain = \
74  self.second_chain, self.first_chain
75 
76  def show(self):
77  s = "Cross Link: %s %s %d - %s %s %d. Distance %f" % (self.first_id,
78  self.first_chain, self.first_residue,
79  self.second_id, self.second_chain,
80  self.second_residue, self.distance)
81  return s
82 
83  def get_name(self):
84  """
85  Generate a unique name for the restraint.
86  @note: The name cannot start with a number, upsets sqlite
87  """
88  name = "cl_%s_%s%d_%s_%s%d" % (self.first_id, self.first_chain,
89  self.first_residue, self.second_id,
90  self.second_chain, self.second_residue)
91  return name
92 
93 
94 def build_xlinks_graph(xlinks_dict):
95  """
96  Build a set of nodes and edges from the set of crosslinking
97  restraints.
98  @param xlinks_dict a XlinksDict class
99 
100  """
101  subunits = set()
102  edges = []
103  for key in xlinks_dict.keys():
104  subunits.add(key[0])
105  subunits.add(key[1])
106  edge = sorted([key[0], key[1]])
107  if edge not in edges:
108  edges.append(edge)
109  log.debug("Subunits %s", subunits)
110  log.debug("Edges %s", edges)
111  return subunits, edges
112 
113 
114 class XlinksDict(dict):
115 
116  """
117  Description of crosslinking restraints as a python
118  dictionary.
119  The keys are a pair with the ids of the cross-linked subunits.
120  Note: The pairs are considered in alphabetic order
121  """
122 
123  def add(self, xlink):
124  """
125  Add a xlink. It is ensured that the id of the first element is
126  is lower that the second
127  """
128  if xlink.second_id < xlink.first_id:
129  xlink.swap()
130  key = (xlink.first_id, xlink.second_id)
131  if key not in self.keys():
132  self[key] = []
133  self[key].append(xlink)
134 
135  def get_xlinks_for_pair(self, pair_ids):
136  """
137  @param pair_ids Ids fo the subunits that are cross-linked
138  """
139  try:
140  xlinks_list = self[pair_ids]
141  ys = [xl.clone() for xl in xlinks_list]
142 
143  except KeyError:
144  try:
145  xlinks_list = self[(pair_ids[1], pair_ids[0])]
146  ys = [xl.clone() for xl in xlinks_list]
147  # Ensure that the first element of the pair corresponds to the
148  # first id in the xlink
149  for y in ys:
150  y.swap()
151  except KeyError as e:
152  raise e
153  return ys
154 
155  """
156  try:
157  return self[pair_ids]
158  except KeyError:
159  print pair_ids,"NOT FOUND, swapping to",(pair_ids[1], pair_ids[0])
160  # If not found, invert the pair
161  xlinks_list = self[(pair_ids[1], pair_ids[0])]
162  # swap the xlinks, so the first element of the pair corresponds to the
163  # first id in the xlink
164  for xl in xlinks_list:
165  xl.swap()
166  for xl in xlinks_list:
167  print xl.show()
168  return xlinks_list
169  """
170 
171 
172 class DockOrder (object):
173 
174  """
175  Compute the order of the docking experiments. The order is derived
176  from the cross-linking restraints:
177  1) The subunit with the highest number of cross-links with others
178  works as the first receptor (A)
179  2) All the subunits cross-linked with A are docked into A.
180  3) The next receptor (B) is the subunit that has the highest number of
181  cross-links after A.
182  4) All the subunits cross-linked to B are docked into B (except A)
183  5) The procedure is repeated until there are no more cross-links
184  """
185 
186  def __init__(self,):
187  self.G = nx.Graph()
188 
189  def set_components_and_connections(self, subunits, edges):
190  """
191  Instead of setting the xlink restraints, init the graph directly
192  Example of what to pass with the structure in 3sfd
193  G.add_nodes_from(["A", "B", "C", "D"])
194  G.add_edges_from([("B","C"), ("B","A"), ("B","D"), ("C","D")])
195  ("B","C") means that there are cross-links between B and C
196  """
197  self.G.add_nodes_from(subunits)
198  self.G.add_edges_from(edges)
199 
200  def set_xlinks(self, xlinks_dict):
201  """
202  Sets the xlinks used for computing the docking order
203  @param xlinks_dict XlinksDict class
204  """
205  subunits, edges = build_xlinks_graph(xlinks_dict)
206  self.set_components_and_connections(subunits, edges)
207 
208  def get_docking_order(self):
209  """ return the order to dock components from the cross links """
210  docking_pairs = []
211  degs = self.G.degree(self.G.nodes())
212  log.debug("Degrees: %s", degs)
213  sorted_degrees = sorted(((v,k) for (k,v) in dict(degs).items()),
214  reverse=True)
215 
216  receptors_considered = []
217  for degree, node in sorted_degrees:
218  for n in self.G.neighbors(node):
219  if not n in receptors_considered:
220  docking_pairs.append((node, n))
221  receptors_considered.append(node)
222  log.info("The suggested order for the docking pairs is %s",
223  docking_pairs)
224  return docking_pairs
225 
226 
228 
229  """
230  Puts two subunits together using the Xlinkins restraints. The solutions
231  offered by this class are just an initial position of the components
232  to be fed to HEX
233  """
234 
235  def clear_xlinks(self):
236  self.xlinks_list = []
237 
238  def set_xlinks(self, xlinks_list):
239  """
240  Sets the xlinks used for the computation fo the initial rough
241  docking solution
242  @param xlinks_list A list of Xlink classes
243  residue1 belongs to the receptor and residue2 belongs to the ligand
244  """
245  self.xlinks_list = xlinks_list
246 
247  def set_pdbs(self, fn_receptor, fn_ligand):
248  """
249  Set the name of the PDB files of the receptor and the ligand
250  @param fn_receptor
251  @param fn_ligand
252  """
253  sel = atom.ATOMPDBSelector()
254  self.m_receptor = IMP.Model()
255  self.h_receptor = atom.read_pdb(fn_receptor, self.m_receptor, sel)
256  self.m_ligand = IMP.Model()
257  self.h_ligand = atom.read_pdb(fn_ligand, self.m_ligand, sel)
258 
259  def set_hierarchies(self, h_receptor, h_ligand):
260  """
261  Set the hierarchies (atom.Hierarchy objects) for the receptor and
262  the ligand
263  @param h_receptor
264  @param h_ligand
265  """
266  self.h_receptor = h_receptor
267  self.h_ligand = h_ligand
268 
269  def set_rigid_bodies(self, rb_receptor, rb_ligand):
270  """
271  Sets the rigid bodies (core.RigidBody objects) for the receptor and
272  the ligand
273  @param rb_receptor
274  @param rb_ligand
275  """
276  self.rb_receptor = rb_receptor
277  self.rb_ligand = rb_ligand
278 
279  def move_ligand(self):
280  """
281  Movest the ligand close to the receptor based on the xlinks
282  provided by set_xlinks()
283  """
284  log.debug("Moving ligand close to the receptor using the xlinks")
285  n = len(self.xlinks_list)
286  if n == 1:
287  self.move_one_xlink()
288  elif n >= 2:
289  self.move_xlinks()
290 
291  def get_residue_particle(self, h, ch, res):
292  """
293  Get the particle representing a residue in a hierarchy
294  @param h atom.Hierarchy containing the residue
295  @param ch The chain id
296  @param res index of the residue
297  """
298  s = IMP.atom.Selection(h, chain=ch, residue_index=res)
299  return s.get_selected_particles()[0]
300 
301  def get_residue_coordinates(self, h, ch, res):
302  """
303  Get the coordinates for a residue in a molecular hierarchy
304  @param h atom.Hierarchy object
305  @param ch The chain id
306  @param res Residue index
307  """
308  p = self.get_residue_particle(h, ch, res)
309  return core.XYZ(p).get_coordinates()
310 
311  def write_ligand(self, fn):
312  """
313  Write a pdb file the coordinates of the ligand
314  @param fn
315  """
316  atom.write_pdb(self.h_ligand, fn)
317 
318  def move_one_xlink(self):
319  """
320  Put the residues in a random distance between 0 and the maximum
321  cross-linkin distance
322  """
323  xl = self.xlinks_list[0]
324  center = self.get_residue_coordinates(self.h_receptor, xl.first_chain,
325  xl.first_residue)
326  sph = alg.Sphere3D(center, xl.distance)
327  v = alg.get_random_vector_in(sph)
328  ref = self.rb_ligand.get_reference_frame()
329  coords = ref.get_transformation_to().get_translation()
330  R = ref.get_transformation_to().get_rotation()
331  lig = self.get_residue_coordinates(self.h_ligand, xl.second_chain,
332  xl.second_residue)
333  log.debug("Ligand residue before moving %s", lig)
334  displacement = v - lig
335  T = alg.Transformation3D(R, coords + displacement)
336  self.rb_ligand.set_reference_frame(alg.ReferenceFrame3D(T))
337  new_coords = self.get_residue_coordinates(
338  self.h_ligand, xl.second_chain,
339  xl.second_residue)
340  log.debug("ligand after moving %s", new_coords)
341 
342  def move_xlinks(self, ):
343  """
344  Function equivalent to move_one_xlink() for the case where there
345  are more than one cross-link restraints available.
346  Puts the ligand residues as close as possible to the receptor
347  residues
348  """
349  rec_coords = []
350  lig_coords = []
351  for xl in self.xlinks_list:
352  c = self.get_residue_coordinates(self.h_receptor, xl.first_chain,
353  xl.first_residue)
354  rec_coords.append(c)
355  c = self.get_residue_coordinates(self.h_ligand, xl.second_chain,
356  xl.second_residue)
357  lig_coords.append(c)
358  log.debug("Receptor residues before moving %s", rec_coords)
359  log.debug("Ligand residues before moving %s", lig_coords)
360  ref = self.rb_ligand.get_reference_frame()
361  Tr = alg.get_transformation_aligning_first_to_second(lig_coords,
362  rec_coords)
363  T = ref.get_transformation_to()
364  newT = alg.compose(Tr, T)
365  self.rb_ligand.set_reference_frame(alg.ReferenceFrame3D(newT))
366 
367  moved_lig_coords = []
368  for xl in self.xlinks_list:
369  c = self.get_residue_coordinates(self.h_ligand, xl.second_chain,
370  xl.second_residue)
371  moved_lig_coords.append(c)
372  log.debug("Ligand residues after moving %s", moved_lig_coords)
Utility functions to extract min/max from the inputs.
Definition: argminmax.py:1
Class for storing model, its restraints, constraints, and particles.
Definition: Model.h:72
Select all non-alternative ATOM records.
Definition: pdb.h:64
A decorator for a particle with x,y,z coordinates.
Definition: XYZ.h:30
std::ostream & show(Hierarchy h, std::ostream &out=std::cout)
Print the hierarchy using a given decorator to display each node.
Basic functionality that is expected to be used by a wide variety of IMP users.
General purpose algebraic and geometric methods that are expected to be used by a wide variety of IMP...
Functionality for loading, creating, manipulating and scoring atomic structures.
Select hierarchy particles identified by the biological name.
Definition: Selection.h:66