IMP  2.1.0
The Integrative Modeling Platform
buildxlinks.py
1 """@namespace IMP.em2d.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.em2d.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 
23 
24 class Xlink:
25  """
26  Class defining a cross-link
27  """
28 
29  def __init__(self, id1, chain1, residue1,
30  id2, chain2, residue2,
31  distance ):
32  """
33  Initialize the class
34  @param[in] id1 Id of the first component of the cross-link
35  @param[in] chain1 Chain of the first component
36  @param[in] residue1 Residue cross-linked in the first component
37  @param[in] id2 Id of the second component of the cross-link
38  @param[in] chain2 chain1 Chain of the second component
39  @param[in] residue2 Residue cross-linked in the second component
40  @param[in] distance Maximum distance
41  """
42  self.first_id = id1
43  self.first_chain = chain1
44  self.first_residue = residue1
45  self.second_id = id2
46  self.second_chain = chain2
47  self.second_residue = residue2
48  self.distance = distance
49 
50  def clone(self):
51  xl = Xlink( self.first_id, self.first_chain, self.first_residue,
52  self.second_id, self.second_chain, self.second_residue,
53  self.distance)
54  return xl
55 
56  def __eq__(self, other):
57  if self.first_id != other.first_id or \
58  self.second_id != other.second_id or \
59  self.first_chain != other.first_chain or \
60  self.second_chain != other.second_chain or \
61  self.first_residue != other.first_residue or \
62  self.second_residue != other.second_residue or \
63  abs(other.distance-self.distance) > 0.01:
64  return False
65  return True
66 
67  def swap(self):
68  """
69  swaps the order of the residues in the restraint
70  """
71  self.first_id, self.second_id = self.second_id, self.first_id
72  self.first_residue, self.second_residue = \
73  self.second_residue, self.first_residue
74  self.first_chain, self.second_chain = \
75  self.second_chain, self.first_chain
76 
77 
78  def show(self):
79  s = "Cross Link: %s %s %d - %s %s %d. Distance %f" % (self.first_id,
80  self.first_chain, self.first_residue,
81  self.second_id, self.second_chain,
82  self.second_residue, self.distance )
83  return s
84 
85  def get_name(self):
86  """
87  Generate a unique name for the restraint.
88  @note: The name cannot start with a number, upsets sqlite
89  """
90  name = "cl_%s_%s%d_%s_%s%d" % (self.first_id, self.first_chain,
91  self.first_residue, self.second_id,
92  self.second_chain, self.second_residue)
93  return name
94 
95 
96 
97 def build_xlinks_graph(xlinks_dict):
98  """
99  Build a set of nodes and edges from the set of crosslinking
100  restraints.
101  @param xlinks_dict a XlinksDict class
102 
103  """
104  subunits = set()
105  edges = []
106  for key in xlinks_dict.keys():
107  subunits.add(key[0])
108  subunits.add(key[1])
109  edge = [key[0],key[1]]
110  edge.sort()
111  if edge not in edges:
112  edges.append(edge)
113  log.debug("Subunits %s", subunits)
114  log.debug("Edges %s", edges)
115  return subunits, edges
116 
117 class XlinksDict(dict):
118  """
119  Description of crosslinking restraints as a python
120  dictionary.
121  The keys are a pair with the ids of the cross-linked subunits.
122  Note: The pairs are considered in alphabetic order
123  """
124  def add(self, xlink):
125  """
126  Add a xlink. It is ensured that the id of the first element is
127  is lower that the second
128  """
129  if xlink.second_id < xlink.first_id:
130  xlink.swap()
131  key = (xlink.first_id, xlink.second_id)
132  if key not in self.keys():
133  self[key] = []
134  self[key].append(xlink)
135 
136  def get_xlinks_for_pair(self, pair_ids):
137  """
138  @param pair_ids Ids fo the subunits that are cross-linked
139  """
140  try:
141  xlinks_list = self[pair_ids]
142  ys = [xl.clone() for xl in xlinks_list]
143 
144  except KeyError:
145  try:
146  xlinks_list = self[(pair_ids[1], pair_ids[0])]
147  ys = [xl.clone() for xl in xlinks_list]
148  # Ensure that the first element of the pair corresponds to the
149  # first id in the xlink
150  for y in ys:
151  y.swap()
152  except KeyError,e:
153  raise e
154  return ys
155 
156 
157 
158  """
159  try:
160  return self[pair_ids]
161  except KeyError:
162  print pair_ids,"NOT FOUND, swapping to",(pair_ids[1], pair_ids[0])
163  # If not found, invert the pair
164  xlinks_list = self[(pair_ids[1], pair_ids[0])]
165  # swap the xlinks, so the first element of the pair corresponds to the
166  # first id in the xlink
167  for xl in xlinks_list:
168  xl.swap()
169  for xl in xlinks_list:
170  print xl.show()
171  return xlinks_list
172  """
173 
174 class DockOrder (object):
175  """
176  Compute the order of the docking experiments. The order is derived
177  from the cross-linking restraints:
178  1) The subunit with the highest number of cross-linkns with others
179  works as the first receptor (A)
180  2) All the subunits cross-linked with A are docked into A.
181  3) The next receptor (B) is the subunit that has the highest number of
182  cross-links after A.
183  4) All the subunits cross-linked to B are docked into B (expect A)
184  5) The procedure is repeated until there are no more cross-links
185  """
186 
187  def __init__(self,):
188  self.G = nx.Graph()
189 
190  def set_components_and_connections(self, subunits, edges):
191  """
192  Instead of setting the xlink restraints, init the graph directly
193  Example of what to pass with the structure in 3sfd
194  G.add_nodes_from(["A", "B", "C", "D"])
195  G.add_edges_from([("B","C"), ("B","A"), ("B","D"), ("C","D")])
196  ("B","C") means that there are cross-links between B and C
197  """
198  self.G.add_nodes_from(subunits)
199  self.G.add_edges_from(edges)
200 
201  def set_xlinks(self, xlinks_dict):
202  """
203  Sets the xlinks used for computing the docking order
204  @param xlinks_dict XlinksDict class
205  """
206  subunits, edges = build_xlinks_graph(xlinks_dict)
207  self.set_components_and_connections(subunits, edges)
208 
209  def get_docking_order(self):
210  """ return the order to dock components from the cross links """
211  docking_pairs = []
212  degs = self.G.degree(self.G.nodes())
213  log.debug("Degrees: %s",degs )
214  sorted_degrees = [(v,k) for v,k in zip(degs.values(),degs.keys())]
215  sorted_degrees.sort()
216  sorted_degrees.reverse() # descending order
217 
218  receptors_considered = []
219  for degree, node in sorted_degrees:
220  for n in self.G.neighbors(node):
221  if not n in receptors_considered:
222  docking_pairs.append((node, n))
223  receptors_considered.append(node)
224  log.info("The suggested order for the docking pairs is %s",
225  docking_pairs)
226  return docking_pairs
227 
228 
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  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 intial 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.kernel.Model()
255  self.h_receptor = atom.read_pdb(fn_receptor, self.m_receptor, sel)
256  self.m_ligand = IMP.kernel.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(self.h_ligand, xl.second_chain,
338  xl.second_residue)
339  log.debug("ligand after moving %s", new_coords)
340 
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)
void show(const cv::Mat &m, std::ostream &out=std::cout)
Prints a OpenCV matrix.
Utility functions to extract min/max from the inputs.
Definition: argminmax.py:1
Select all non-alternative ATOM records.
Definition: pdb.h:63
A decorator for a particle with x,y,z coordinates.
Definition: XYZ.h:32
See IMP.core for more information.
See IMP.algebra for more information.
See IMP.atom for more information.
Class for storing model, its restraints, constraints, and particles.