IMP  2.0.0
The Integrative Modeling Platform
dock_with_crosslinks.py
1 ## \example kernel/dock_with_crosslinks.py
2 ## This example docks several proteins using excluded volume and crosslinking terms. To set the pdbs and crosslinks to use, edit the data at the start of the python script.
3 ##
4 
5 import IMP
6 import IMP.core
7 import IMP.algebra
8 import IMP.base
9 import IMP.atom
10 import IMP.container
11 import random
12 import sys
13 try:
14  import networkx as nx
15  import networkx.algorithms
16 except:
17  print "Script requires networkx to run"
18  sys.exit()
19 
20 IMP.base.setup_from_argv(sys.argv, "Dock several proteins using excluded volume and crosslinking")
21 
22 # remove internal checks
23 IMP.base.set_check_level(IMP.base.USAGE)
24 pdbs=[IMP.get_example_path('dock_data/chainf.pdb'),
25  IMP.get_example_path('dock_data/chaind.pdb')]
26 xlinks=[
27  (256, "D", 89, "F", 17, 0.1),
28  (259, "D", 84, "F", 17, 0.1),
29  (259, "D", 89, "F", 17, 0.1),
30  (269, "D", 89, "F", 17, 0.1),
31  (128, "F", 275, "D", 17, 0.1),
32  (130, "F", 275, "D", 17, 0.1),
33  (144, "F", 269, "D", 17, 0.1),
34  (150, "F", 256, "D", 17, 0.1),
35  (150, "F", 259, "D", 17, 0.1),
36  (150, "F", 269, "D", 17, 0.1),
37  (150, "F", 179, "D", 17, 0.1),
38  (166, "F", 186, "D", 17, 0.1),
39  (170, "D", 89, "F", 17, 0.1),
40  (179, "D", 84, "F", 17, 0.1)]
41 
42 
43 def read_pdbs(list_pdb_file):
44  """read pdbs from an external list file
45  create a simplified representation"""
46  chains=[]
47  chain_id={}
48  for pdb in pdbs:
49  prot = IMP.atom.read_pdb(pdb, m,
50  IMP.atom.CAlphaPDBSelector())
51  ch=IMP.atom.get_by_type(prot, IMP.atom.CHAIN_TYPE)
52  chains+=ch
53  print ch
54 
55  for c in chains:
56  chain_id[c]=IMP.atom.Chain(c).get_id()
57  return chains,chain_id
58 
59 def setup_radii(chains):
60  "setup the radii for each residue"
61  for c in chains:
62  IMP.atom.create_rigid_body(c)
63  for res in c.get_children():
64  dres=IMP.core.XYZR(res)
65  rt=IMP.atom.Residue(res).get_residue_type()
66  vol=IMP.atom.get_volume_from_residue_type(rt)
67  rg=IMP.algebra.get_ball_radius_from_volume_3d(vol)
68  dres.set_radius(rg)
69 
70 
71 def add_excluded_volume(chains):
72  """add excluded volume score on the coarse grained c-alpha model
73  residues are represented by beads with conveniently defined excluded volume radius"""
74  rs = IMP.RestraintSet('excluded_volume')
75  kappa_=0.1
76  lsc= IMP.container.ListSingletonContainer(m)
77  IMP.atom.get_by_type
78  for c in chains:
79  lsc.add_particles(IMP.atom.get_leaves(c))
80  evr=IMP.core.ExcludedVolumeRestraint(lsc,kappa_)
81  rs.add_restraint(evr)
82  m.add_restraint(rs)
83  restraints['excluded_volume']=rs
84 
85 
86 def setup_move_them_all_MonteCarlo_internal(chains,mc_dx=0.3,mc_dang=0.1,mc_kt=1.0):
87  """setup rigid body monte carlo move set
88  mc_dx is the traslation move
89  mc_dang is the rotation angle"""
90  mc=IMP.core.MonteCarlo(m)
91 
92  mc.set_return_best(False)
93  mc.set_kt(10)
94  for c in chains:
95  p=IMP.atom.get_leaves(c)[0]
96  rb=IMP.core.RigidMember(p).get_rigid_body()
97  mv= IMP.core.RigidBodyMover(rb, mc_dx, mc_dang)
98  mc.add_mover(mv)
99  return mc
100 
101 def setup_move_them_all_MonteCarlo_external(connected_chain_list,rb_ext_list=[],
102  mc_dx=0.3,mc_dang=0.1,mc_kt=1.0,return_best=False):
103  """setup rigid body monte carlo move set
104  mc_dx is the traslation move
105  mc_dang is the rotation angle"""
106 
107  print 'RIGID EXT BODY LIST', [r.get_name() for r in rb_ext_list]
108  for rb in rb_ext_list:
109  IMP.core.RigidBody.teardown_particle(rb)
110 
111  rb_ext_list=[]
112 
113  mc=IMP.core.MonteCarlo(m)
114  mc.set_return_best(return_best)
115  mc.set_kt(mc_kt)
116 
117  print connected_chain_list
118  for c in connected_chain_list:
119  rb_tmp_list=[]
120  for s in c:
121  p=IMP.atom.get_leaves(s)[0]
122  rba=IMP.core.RigidMember(p).get_rigid_body()
123  rb_tmp_list.append(rba)
124 
125  print 'RIGID BODY LIST', [r.get_name( ) for r in rb_tmp_list]
126  p=IMP.Particle(m)
127  p.set_name("root rigid body "+str(chain_id[s]))
128  rb=IMP.core.RigidBody.setup_particle(p,rb_tmp_list)
129  rb_ext_list.append(rb)
130  mv= IMP.core.RigidBodyMover(rb, mc_dx, mc_dang)
131  mc.add_mover(mv)
132  return mc,rb_ext_list
133 
134 def setup_move_them_all_MonteCarlo_internal_2(chains,rb_ext_list=[],
135  mc_dx=0.3,mc_dang=0.1,mc_kt=1.0):
136  """setup rigid body monte carlo move set
137  mc_dx is the traslation move
138  mc_dang is the rotation angle"""
139 
140  print 'RIGID EXT BODY LIST', [r.get_name() for r in rb_ext_list]
141  for rb in rb_ext_list:
142  IMP.core.RigidBody.teardown_particle(rb)
143 
144  rb_ext_list=[]
145 
146  mc=IMP.core.MonteCarlo(m)
147 
148  mc.set_return_best(False)
149  mc.set_kt(10)
150  for c in chains:
151  p=IMP.atom.get_leaves(c)[0]
152  rb=IMP.core.RigidMember(p).get_rigid_body()
153  rb_ext_list.append(rb)
154  mv= IMP.core.RigidBodyMover(rb, mc_dx, mc_dang)
155  mc.add_mover(mv)
156  return mc,rb_ext_list
157 
158 def shuffle_configuration(chains,bounding_box_length=200.0):
159  "shuffle configuration, used to restart the optimization"
160  hbbl=bounding_box_length/2
161  for c in chains:
162  p=IMP.atom.get_leaves(c)[0]
163  rb=IMP.core.RigidMember(p).get_rigid_body()
164  ub = IMP.algebra.Vector3D(-hbbl,-hbbl,-hbbl)
165  lb = IMP.algebra.Vector3D( hbbl, hbbl, hbbl)
166  bb = IMP.algebra.BoundingBox3D(ub, lb)
167  translation = IMP.algebra.get_random_vector_in(bb)
168  rotation = IMP.algebra.get_random_rotation_3d()
169  transformation = IMP.algebra.Transformation3D(rotation, translation)
170  rb.set_reference_frame(IMP.algebra.ReferenceFrame3D(transformation))
171 
172 
173 def add_restraints():
174  """Handle the restraints defined above
175  sintax: resid1 chain1 resid2 chain2 max_distance strength_of_potential
176  example: 226 H 272 N 15.0 10.0"""
177  rs=IMP.RestraintSet('cross_links')
178  restraints={}
179  restraints_map={}
180  for tokens in xlinks:
181  r1=tokens[0]
182  c1=tokens[1]
183  r2=tokens[2]
184  c2=tokens[3]
185  distance=tokens[4]
186  strength=tokens[5]
187  try:
188  s1=IMP.atom.Selection(chains, residue_index=r1, chains=c1)
189  p1=s1.get_selected_particles()[0]
190  except:
191  print "residue %d of chain %s not found" % (r1,c1)
192  continue
193  try:
194  s2=IMP.atom.Selection(chains, residue_index=r2, chains=c2)
195  p2=s2.get_selected_particles()[0]
196  except:
197  print "residue %d of chain %s not found" % (r2,c2)
198  continue
199  d=distance
200  k=strength
201  hub= IMP.core.HarmonicUpperBound(d,k)
202  df= IMP.core.DistancePairScore(hub)
203  dr= IMP.core.PairRestraint(df, (p1, p2))
204  rs.add_restraint(dr)
205  a='cross_link-'+`r1`+'.'+c1+'-'+`r2`+'.'+c2+'.lower'
206  restraints_map
207  restraints[a]=dr
208  restraints_map[a]=(c1,c2)
209 
210 
211  m.add_restraint(rs)
212  restraints['cross_links']=rs
213  return restraints,restraints_map
214 
215 def setup_md(temp=300.0, tau=0.01):
216  md=IMP.atom.MolecularDynamics()
217  md.set_model(m)
218  md.assign_velocities(temp)
219  md.set_time_step(1.0)
220  st = IMP.atom.LangevinThermostatOptimizerState(md.get_simulation_particles(), temp,tau)
221  md.add_optimizer_state(st)
222  return md
223 
224 def init_pdb(prefix='models'):
225  "append models to a pdb file"
226  pdbfile = prefix+'.pdb'
227  flpdb=open(pdbfile,'w')
228  flpdb.close()
229 
230 def write_pdb(prefix='models'):
231  "append zero score models (solutions) to a pdb file"
232  pdbfile = prefix+'.pdb'
233  flpdb=open(pdbfile,'a')
234  tmpfile=".tmp.pdb"
235  IMP.atom.write_pdb(chains,tmpfile)
236  fltmp=open(tmpfile,'r')
237  pdb=fltmp.read()
238  flpdb.write(pdb)
239  flpdb.close()
240  fltmp.close()
241 
242 def chain_pair_list(chains,nop_cutoff=0):
243  pair_list=[]
244  for c in chains:
245  sc1= IMP.container.ListSingletonContainer(m)
246  sc1.add_particles(IMP.atom.get_leaves(c))
247  for s in chains:
248  if chain_id[s]>chain_id[c]:
249  sc2= IMP.container.ListSingletonContainer(m)
250  sc2.add_particles(IMP.atom.get_leaves(s))
251  aaa= IMP.container.CloseBipartitePairContainer(sc1,sc2,0.0)
252  m.update()
253  nop=aaa.get_number_of_particle_pairs()
254  if nop > nop_cutoff :
255  pair_list.append((c,s))
256  else:
257  pair_list.append((c,c))
258  pair_list.append((s,s))
259  return pair_list
260 
261 
262 def chain_pair_list_based_on_restraint(chains,nop_cutoff=0):
263  pair_list=[]
264  id_chain=invert_dict(chain_id)
265  for i in restraints.keys():
266  (chaina_id,chainb_id)=restraints_map[i]
267  chaina=invert_dict(chaina_id)
268  chainb=invert_dict(chainb_id)
269 
270 
271  for c in chains:
272  sc1= IMP.container.ListSingletonContainer(m)
273  sc1.add_particles(IMP.atom.get_leaves(c))
274  for s in chains:
275  if chain_id[s]>chain_id[c]:
276  sc2= IMP.container.ListSingletonContainer(m)
277  sc2.add_particles(IMP.atom.get_leaves(s))
278  aaa= IMP.container.CloseBipartitePairContainer(sc1,sc2,0.0)
279  m.update()
280  nop=aaa.get_number_of_particle_pairs()
281  if nop > nop_cutoff :
282  pair_list.append((c,s))
283  else:
284  pair_list.append((c,c))
285  pair_list.append((s,s))
286  return pair_list
287 
288 
289 def create_graph(chains,nop_cutoff=0):
290  connected=[]
291  while len(connected)<=1:
292  G=nx.Graph()
293  pair_list=chain_pair_list(chains,nop_cutoff)
294  G.add_edges_from(pair_list)
295  connected=nx.algorithms.connected_components(G)
296  G.clear()
297  nop_cutoff+=10
298  return connected,nop_cutoff
299 
300 def invert_dict(d):
301  inv = {}
302  for k, v in d.iteritems():
303  keys = inv.setdefault(v, [])
304  keys.append(k)
305  return inv
306 
307 
308 m = IMP.Model()
309 chain_id={}
310 rb_ext_list=[]
311 restraints={}
312 restraints_map={}
313 
314 
315 
316 init_pdb("models")
317 init_pdb("solutions")
318 init_pdb("best")
319 chains,chain_id=read_pdbs("pdb.list")
320 write_pdb("best")
321 setup_radii(chains)
322 add_excluded_volume(chains)
323 restraints,restraints_map=add_restraints()
324 shuffle_configuration(chains,300.0)
325 
326 
327 
328 delta_r=5.0
329 delta_a=0.2
330 delta_r_best=0.5
331 delta_a_best=0.1
332 temp=1.0
333 
334 
335 mc_int=setup_move_them_all_MonteCarlo_internal(chains,delta_r,delta_a,temp)
336 write_pdb("models")
337 
338 
339 
340 nsteps_int=100
341 nsteps_ext=100
342 
343 cg=IMP.core.ConjugateGradients()
344 cg.set_model(m)
345 
346 number_of_steps = 10000
347 
348 if IMP.base.get_bool_flag("run_quick_test"):
349  number_of_steps = 2
350 
351 for steps in range(number_of_steps):
352 
353  nop_cutoff=10000000
354  connected_chains_list,nop_cutoff=create_graph(chains,nop_cutoff)
355  mc_int,rb_ext_list = setup_move_them_all_MonteCarlo_external(connected_chains_list,rb_ext_list,delta_r,delta_a,temp,return_best=False)
356  mc_int.optimize(nsteps_ext)
357  print ' INTERNAL', m.evaluate(False)
358  write_pdb("models")
359 
360 
361  nop_cutoff=0
362  connected_chains_list,nop_cutoff=create_graph(chains,nop_cutoff)
363  mc_ext,rb_ext_list = setup_move_them_all_MonteCarlo_external(connected_chains_list,rb_ext_list,delta_r,delta_a,temp,return_best=False)
364  mc_ext.optimize(nsteps_ext)
365  print ' EXTERNAL', m.evaluate(False)
366  write_pdb("models")
367 
368  nfs=mc_ext.get_number_of_forward_steps()
369  acceptance=float(nfs)/nsteps_ext
370 
371  if acceptance<0.05:
372  delta_r=delta_r/2
373  if acceptance>0.3:
374  delta_r=delta_r*2
375  print 'ACCEPTANCE 1', acceptance, delta_r, delta_a
376 
377 
378  nop_cutoff=10000000
379  connected_chains_list,nop_cutoff=create_graph(chains,nop_cutoff)
380  mc_int,rb_ext_list = setup_move_them_all_MonteCarlo_external(connected_chains_list,rb_ext_list,delta_r_best,delta_a_best,temp,return_best=True)
381  mc_int.optimize(nsteps_ext)
382  print ' INTERNAL', m.evaluate(False)
383  write_pdb("models")
384 
385 
386  nop_cutoff=0
387  connected_chains_list,nop_cutoff=create_graph(chains,nop_cutoff)
388  mc_ext,rb_ext_list = setup_move_them_all_MonteCarlo_external(connected_chains_list,rb_ext_list,delta_r_best,delta_a_best,temp,return_best=True)
389  mc_ext.optimize(nsteps_ext)
390  print ' EXTERNAL', m.evaluate(False)
391  write_pdb("models")
392 
393  nfs=mc_ext.get_number_of_forward_steps()
394  acceptance=float(nfs)/nsteps_ext
395 
396  if acceptance<0.5:
397  delta_r_best=delta_r_best/2
398 
399  if acceptance>0.8:
400  delta_r_best=delta_r_best*2
401  print 'ACCEPTANCE 2', acceptance, delta_r_best, delta_a_best
402 
403  cg.optimize(10)
404 
405  score=m.evaluate(False)
406  for i in restraints.keys():
407  print i, restraints[i].evaluate(False)
408 
409 
410  if score==0.0:
411  write_pdb("solutions")
412  shuffle_configuration(chains,300.0)
413  exit()
414 
415 
416 
417 
418 print 'ALL DONE'