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