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