IMP  2.0.1
The Integrative Modeling Platform
buildxlinks.py
1 
2 import IMP
3 import IMP.atom as atom
4 import IMP.algebra as alg
5 import IMP.core as core
6 import sys
7 import logging
8 from IMP.em2d.argminmax import keymax
9 log = logging.getLogger("buildxlinks")
10 
11 import networkx as nx
12 
13 try:
14  set = set
15 except NameError:
16  from sets import Set as set
17 
18 
19 
20 
21 class Xlink:
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 id1 Id of the first component of the cross-link
32  @param chain1 Chain of the first component
33  @param residue1 Residue cross-linked in the first component
34  @param id2 Id of the second component of the cross-link
35  @param chain2 chain1 Chain of the second component
36  @param residue2 Residue cross-linked in the second component
37  @param 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  @param
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 
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 
95 def build_xlinks_graph(xlinks_dict):
96  """
97  Build a set of nodes and edges from the set of crosslinking
98  restraints.
99  @param xlinks_dict a XlinksDict class
100 
101  """
102  subunits = set()
103  edges = []
104  for key in xlinks_dict.keys():
105  subunits.add(key[0])
106  subunits.add(key[1])
107  edge = [key[0],key[1]]
108  edge.sort()
109  if edge not in edges:
110  edges.append(edge)
111  log.debug("Subunits %s", subunits)
112  log.debug("Edges %s", edges)
113  return subunits, edges
114 
115 class XlinksDict(dict):
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  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  @param
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,e:
152  raise e
153  return ys
154 
155 
156 
157  """
158  try:
159  return self[pair_ids]
160  except KeyError:
161  print pair_ids,"NOT FOUND, swapping to",(pair_ids[1], pair_ids[0])
162  # If not found, invert the pair
163  xlinks_list = self[(pair_ids[1], pair_ids[0])]
164  # swap the xlinks, so the first element of the pair corresponds to the
165  # first id in the xlink
166  for xl in xlinks_list:
167  xl.swap()
168  for xl in xlinks_list:
169  print xl.show()
170  return xlinks_list
171  """
172 
173 class DockOrder (object):
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-linkns 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 (expect 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 = [(v,k) for v,k in zip(degs.values(),degs.keys())]
214  sorted_degrees.sort()
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  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  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 fo the intial 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(self.h_ligand, xl.second_chain,
337  xl.second_residue)
338  log.debug("ligand after moving %s", new_coords)
339 
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)