IMP logo
IMP Reference Guide  2.5.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)
214  for v, k in zip(degs.values(), degs.keys())])
215  sorted_degrees.reverse() # descending order
216 
217  receptors_considered = []
218  for degree, node in sorted_degrees:
219  for n in self.G.neighbors(node):
220  if not n in receptors_considered:
221  docking_pairs.append((node, n))
222  receptors_considered.append(node)
223  log.info("The suggested order for the docking pairs is %s",
224  docking_pairs)
225  return docking_pairs
226 
227 
229 
230  """
231  Puts two subunits together using the Xlinkins restraints. The solutions
232  offered by this class are just an initial position of the components
233  to be fed to HEX
234  """
235 
236  def clear_xlinks(self):
237  self.xlinks_list = []
238 
239  def set_xlinks(self, xlinks_list):
240  """
241  Sets the xlinks used for the computation fo the intial rough
242  docking solution
243  @param xlinks_list A list of Xlink classes
244  residue1 belongs to the receptor and residue2 belongs to the ligand
245  """
246  self.xlinks_list = xlinks_list
247 
248  def set_pdbs(self, fn_receptor, fn_ligand):
249  """
250  Set the name of the PDB files of the receptor and the ligand
251  @param fn_receptor
252  @param fn_ligand
253  """
254  sel = atom.ATOMPDBSelector()
255  self.m_receptor = IMP.Model()
256  self.h_receptor = atom.read_pdb(fn_receptor, self.m_receptor, sel)
257  self.m_ligand = IMP.Model()
258  self.h_ligand = atom.read_pdb(fn_ligand, self.m_ligand, sel)
259 
260  def set_hierarchies(self, h_receptor, h_ligand):
261  """
262  Set the hierarchies (atom.Hierarchy objects) for the receptor and
263  the ligand
264  @param h_receptor
265  @param h_ligand
266  """
267  self.h_receptor = h_receptor
268  self.h_ligand = h_ligand
269 
270  def set_rigid_bodies(self, rb_receptor, rb_ligand):
271  """
272  Sets the rigid bodies (core.RigidBody objects) for the receptor and
273  the ligand
274  @param rb_receptor
275  @param rb_ligand
276  """
277  self.rb_receptor = rb_receptor
278  self.rb_ligand = rb_ligand
279 
280  def move_ligand(self):
281  """
282  Movest the ligand close to the receptor based on the xlinks
283  provided by set_xlinks()
284  """
285  log.debug("Moving ligand close to the receptor using the xlinks")
286  n = len(self.xlinks_list)
287  if n == 1:
288  self.move_one_xlink()
289  elif n >= 2:
290  self.move_xlinks()
291 
292  def get_residue_particle(self, h, ch, res):
293  """
294  Get the particle representing a residue in a hierarchy
295  @param h atom.Hierarchy containing the residue
296  @param ch The chain id
297  @param res index of the residue
298  """
299  s = IMP.atom.Selection(h, chain=ch, residue_index=res)
300  return s.get_selected_particles()[0]
301 
302  def get_residue_coordinates(self, h, ch, res):
303  """
304  Get the coordinates for a residue in a molecular hierarchy
305  @param h atom.Hierarchy object
306  @param ch The chain id
307  @param res Residue index
308  """
309  p = self.get_residue_particle(h, ch, res)
310  return core.XYZ(p).get_coordinates()
311 
312  def write_ligand(self, fn):
313  """
314  Write a pdb file the coordinates of the ligand
315  @param fn
316  """
317  atom.write_pdb(self.h_ligand, fn)
318 
319  def move_one_xlink(self):
320  """
321  Put the residues in a random distance between 0 and the maximum
322  cross-linkin distance
323  """
324  xl = self.xlinks_list[0]
325  center = self.get_residue_coordinates(self.h_receptor, xl.first_chain,
326  xl.first_residue)
327  sph = alg.Sphere3D(center, xl.distance)
328  v = alg.get_random_vector_in(sph)
329  ref = self.rb_ligand.get_reference_frame()
330  coords = ref.get_transformation_to().get_translation()
331  R = ref.get_transformation_to().get_rotation()
332  lig = self.get_residue_coordinates(self.h_ligand, xl.second_chain,
333  xl.second_residue)
334  log.debug("Ligand residue before moving %s", lig)
335  displacement = v - lig
336  T = alg.Transformation3D(R, coords + displacement)
337  self.rb_ligand.set_reference_frame(alg.ReferenceFrame3D(T))
338  new_coords = self.get_residue_coordinates(
339  self.h_ligand, xl.second_chain,
340  xl.second_residue)
341  log.debug("ligand after moving %s", new_coords)
342 
343  def move_xlinks(self, ):
344  """
345  Function equivalent to move_one_xlink() for the case where there
346  are more than one cross-link restraints available.
347  Puts the ligand residues as close as possible to the receptor
348  residues
349  """
350  rec_coords = []
351  lig_coords = []
352  for xl in self.xlinks_list:
353  c = self.get_residue_coordinates(self.h_receptor, xl.first_chain,
354  xl.first_residue)
355  rec_coords.append(c)
356  c = self.get_residue_coordinates(self.h_ligand, xl.second_chain,
357  xl.second_residue)
358  lig_coords.append(c)
359  log.debug("Receptor residues before moving %s", rec_coords)
360  log.debug("Ligand residues before moving %s", lig_coords)
361  ref = self.rb_ligand.get_reference_frame()
362  Tr = alg.get_transformation_aligning_first_to_second(lig_coords,
363  rec_coords)
364  T = ref.get_transformation_to()
365  newT = alg.compose(Tr, T)
366  self.rb_ligand.set_reference_frame(alg.ReferenceFrame3D(newT))
367 
368  moved_lig_coords = []
369  for xl in self.xlinks_list:
370  c = self.get_residue_coordinates(self.h_ligand, xl.second_chain,
371  xl.second_residue)
372  moved_lig_coords.append(c)
373  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:63
A decorator for a particle with x,y,z coordinates.
Definition: XYZ.h:30
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...
void show(Hierarchy h, std::ostream &out=std::cout)
Print out a molecular hierarchy.
Functionality for loading, creating, manipulating and scoring atomic structures.
Select hierarchy particles identified by the biological name.
Definition: Selection.h:65