IMP  2.2.1
The Integrative Modeling Platform
kernel/nup84.py

Modify solve for a Nup84-like structure using a mix of rigid bodies and coarse grained models using crosslinking data. In addition, show how to visualize restraints and visualize the rejected conformations. Both are useful things to do when trying to figure out why optimization is not converging.

1 ## \example kernel/nup84.py
2 # Modify solve for a Nup84-like structure using a mix of rigid bodies
3 # and coarse grained models using crosslinking data. In
4 # addition, show how to visualize restraints and visualize the
5 # rejected conformations. Both are useful things to do when trying to
6 # figure out why optimization is not converging.
7 
8 import IMP
9 import IMP.atom
10 import IMP.container
11 import IMP.display
12 import IMP.statistics
13 import IMP.example
14 import os
15 import sys
16 
17 # not finished
18 IMP.base.add_bool_flag("run", "Whether to run the program")
19 
20 # parse command line arguments so, eg profiling can be used
21 IMP.base.setup_from_argv(sys.argv, "Nup84 example")
22 
23 if IMP.base.get_bool_flag("run") != "yes":
24  exit(0)
25 
26 # First we define some basic parameters for the modeling effort
27 
28 # the spring constant to use, it doesn't really matter
29 k = 10
30 # the target resolution for the representation
31 resolution = 100
32 # the box to perform everything in, make it flat as it is a 2D structure
33 bb = IMP.algebra.BoundingBox3D(IMP.algebra.Vector3D(-300, -300, -50),
34  IMP.algebra.Vector3D(300, 300, 50))
35 
36 # how many times to try to find a good solution
37 number_of_sampling_attempts = 1
38 number_of_mc_steps = 10000
39 
40 # Create a coarse grained protein with a given name, adding it to universe
41 
42 
43 def add_protein_from_length(model, name, residues, parent, restraints,
44  excluded_volume_particles, optimized_particles):
45  # Create a coarse grained protein with the passed residue information
46  h = IMP.atom.create_protein(model, name, resolution, residues)
47 
48  parent.add_child(h)
49  # Here, each of the constituent particles will be optimized independently
50  leaves = IMP.atom.get_leaves(h)
51  optimized_particles.extend(leaves)
52  excluded_volume_particles.extend(leaves)
53 
54  # Ensure that the various particles of the protein stay connected
55  r = IMP.atom.create_connectivity_restraint([IMP.atom.Selection(c)
56  for c in h.get_children()], k,
57  "connect " + name)
58 
59  if r:
60  # make sure there is an actual restraint
61  restraints.append(r)
62  r.set_maximum_score(k)
63 
64 # Create protein as a rigid body from a pdb file
65 
66 
67 def add_protein_from_pdb(model, name, file, parent, restraints,
68  excluded_volume_particles, optimized_particles):
69  # we should keep the original particles around so they get written
70 
71  # create an atomic representation from the pdb file
72  t = IMP.atom.read_pdb(
73  IMP.get_example_path(os.path.join("data", file)), model,
74  IMP.atom.ATOMPDBSelector())
75  # extract the chain from the file
76  c = IMP.atom.Chain(IMP.atom.get_by_type(t, IMP.atom.CHAIN_TYPE)[0])
77  # there is no reason to use all atoms, just approximate the pdb shape
78  # instead
79  s = IMP.atom.create_simplified_along_backbone(c, resolution / 2.0, True)
80  s.set_name(name)
81  # tear down what is left
82  IMP.atom.destroy(t)
83  # make the simplified structure rigid
84  rb = IMP.atom.create_rigid_body(s)
85  rb.set_coordinates_are_optimized(True)
86  optimized_particles.append(rb)
87  excluded_volume_particles.extend(s.get_children())
88  parent.add_child(s)
89 
90 
91 # Create protein as a rigid body from several pdb file
92 def add_protein_from_pdbs(model, name, files, parent, restraints,
93  excluded_volume_particles, optimized_particles):
94  h = IMP.atom.Hierarchy.setup_particle(IMP.kernel.Particle(model, name))
95  for i, f in enumerate(files):
96  add_protein_from_pdb(model, name + str(i), f, h, restraints,
97  excluded_volume_particles, optimized_particles)
98  r = IMP.atom.create_connectivity_restraint([IMP.atom.Selection(c, hierarchy_types=[IMP.atom.FRAGMENT_TYPE])
99  for c in h.get_children()],
100  k, "connect " + name)
101  if r:
102  restraints.append(r)
103  r.set_maximum_score(k)
104 
105 # Create all the needed representation using the above functions
106 
107 
108 def create_representation(model):
109  restraints = []
110  optimized_particles = []
111  excluded_volume_particles = []
112  universe = IMP.atom.Hierarchy.setup_particle(
113  IMP.kernel.Particle(model, "the universe"))
114 
115  add_protein_from_length(model, "Nup85", 570, universe, restraints,
116  excluded_volume_particles, optimized_particles)
117 
118  # pin the c-terminus
119  ct = IMP.atom.Selection(universe, molecule="Nup85",
120  hierarchy_types=[IMP.atom.FRAGMENT_TYPE],
121  terminus=IMP.atom.Selection.C)
122  d = IMP.core.XYZ(ct.get_selected_particles()[0])
123  d.set_coordinates(IMP.algebra.Vector3D(0, 0, 0))
124  d.set_coordinates_are_optimized(False)
125 
126  add_protein_from_length(model, "Nup84", 460, universe, restraints,
127  excluded_volume_particles, optimized_particles)
128  add_protein_from_length(model, "Nup145C", 442, universe, restraints,
129  excluded_volume_particles, optimized_particles)
130  add_protein_from_length(
131  model, "Nup120", [0, 500, 761], universe, restraints,
132  excluded_volume_particles, optimized_particles)
133  add_protein_from_length(
134  model, "Nup133", [0, 450, 778, 1160], universe, restraints,
135  excluded_volume_particles, optimized_particles)
136  add_protein_from_pdb(model, "Seh1", "seh1.pdb", universe, restraints,
137  excluded_volume_particles, optimized_particles)
138  add_protein_from_pdb(model, "Sec13", "sec13.pdb", universe, restraints,
139  excluded_volume_particles, optimized_particles)
140  return universe, restraints, excluded_volume_particles, optimized_particles
141 
142 
143 def add_distance_restraint(selection0, selection1, name, restraints):
144  r = IMP.atom.create_distance_restraint(selection0, selection1, 0, k, name)
145  r.set_maximum_score(k)
146  restraints.append(r)
147 
148 
149 def encode_data_as_restraints(universe, restraints):
150  s0 = IMP.atom.Selection(
151  hierarchy=universe, hierarchy_types=[IMP.atom.FRAGMENT_TYPE],
152  molecule="Nup145C", residue_indexes=[(0, 423)])
153  s1 = IMP.atom.Selection(
154  hierarchy=universe, hierarchy_types=[IMP.atom.FRAGMENT_TYPE],
155  molecule="Nup84")
156  s2 = IMP.atom.Selection(
157  hierarchy=universe, hierarchy_types=[IMP.atom.FRAGMENT_TYPE],
158  molecule="Sec13")
159  r = IMP.atom.create_connectivity_restraint(
160  [s0, s1, s2], k, "Nup145C Nup84 Sec13")
161  r.set_maximum_score(k)
162  restraints.append(r)
163 
164  add_distance_restraint(
165  IMP.atom.Selection(hierarchy=universe, molecule="Nup145C",
166  residue_indexes=[(0, 423)],
167  hierarchy_types=[
168  IMP.atom.FRAGMENT_TYPE]),
169  IMP.atom.Selection(
170  hierarchy=universe, molecule="Nup85",
171  hierarchy_types=[
172  IMP.atom.FRAGMENT_TYPE]),
173  "Num145C, Nup85", restraints)
174  add_distance_restraint(
175  IMP.atom.Selection(hierarchy=universe, molecule="Nup145C",
176  residue_indexes=[(0, 423)],
177  hierarchy_types=[
178  IMP.atom.FRAGMENT_TYPE]),
179  IMP.atom.Selection(
180  hierarchy=universe, molecule="Nup120",
181  residue_indexes=[(500, 762)],
182  hierarchy_types=[
183  IMP.atom.FRAGMENT_TYPE]),
184  "Nup145C Nup120", restraints)
185  add_distance_restraint(
186  IMP.atom.Selection(hierarchy=universe, molecule="Nup84",
187  hierarchy_types=[
188  IMP.atom.FRAGMENT_TYPE]),
189  IMP.atom.Selection(
190  hierarchy=universe, molecule="Nup133",
191  residue_indexes=[(778, 1160)],
192  hierarchy_types=[
193  IMP.atom.FRAGMENT_TYPE]),
194  "Nup84 Nup133", restraints)
195  add_distance_restraint(
196  IMP.atom.Selection(hierarchy=universe, molecule="Nup85",
197  hierarchy_types=[
198  IMP.atom.FRAGMENT_TYPE]),
199  IMP.atom.Selection(
200  hierarchy=universe, molecule="Seh1",
201  hierarchy_types=[
202  IMP.atom.FRAGMENT_TYPE]),
203  "Nup85 Seh1", restraints)
204  add_distance_restraint(
205  IMP.atom.Selection(hierarchy=universe, molecule="Nup145C",
206  residue_indexes=[(0, 423)],
207  hierarchy_types=[
208  IMP.atom.FRAGMENT_TYPE]),
209  IMP.atom.Selection(
210  hierarchy=universe, molecule="Sec13",
211  hierarchy_types=[
212  IMP.atom.FRAGMENT_TYPE]),
213  "Nup145C Sec13", restraints)
214 
215 
216 # find acceptable conformations of the model
217 def get_configurations(
218  model,
219  restraints,
220  excluded_volume_particles,
221  optimized_particles):
222  #cpc= IMP.container.ClosePairContainer(representation.get_particles(), 0, 10)
223  # evr= IMP.container.PairRestraint(IMP.core.SoftSpherePairScore(k), cpc,
224  # "Excluded Volume")
225  scale = .5
226  mc = IMP.core.MonteCarlo(model)
227  movers = []
228  for p in optimized_particles:
229  if IMP.core.RigidBody.get_is_setup(p):
230  mover = IMP.core.RigidBodyMover(
231  p, IMP.core.XYZR(p).get_radius() * scale,
232  .2 * scale)
233  movers.append(mover)
234  else:
235  mover = IMP.core.BallMover(
236  [p], IMP.core.XYZR(p).get_radius() * scale)
237  movers.append(mover)
238  serial_mover = IMP.core.SerialMover(movers)
239  mc.add_mover(serial_mover)
240  scoring_function = IMP.core.IncrementalScoringFunction(
241  optimized_particles, restraints)
242  scoring_function.add_close_pair_score(IMP.core.SoftSpherePairScore(k), 0.0,
243  excluded_volume_particles)
244 
245  configuration_set = IMP.ConfigurationSet(model)
246  # must write our own sampler as IMP.core.MCCGSampler doesn't handle rigid
247  # bodies
248  for i in range(number_of_sampling_attempts):
249  for p in optimized_particles:
250  IMP.core.XYZ(p).set_coordinates(
251  IMP.algebra.get_random_vector_in(bb))
252  mc.optimize(number_of_mc_steps)
253  if scoring_function.get_had_good_score():
254  configuration_set.save()
255  return configuration_set
256 
257 
258 model = IMP.kernel.Model()
259 universe, restraints, excluded_volume_particles, optimized_particles = create_representation(
260  model)
261 encode_data_as_restraints(universe, restraints)
262 
263 configuration_set = get_configurations(model, restraints,
264  excluded_volume_particles,
265  optimized_particles)
266 
267 print "Found", configuration_set.get_number_of_configurations(), "good configurations"
268 
269 # now lets save them all to a file
270 rmf_file_name = IMP.base.create_temporary_file_name("nup84", ".rmf")
271 rmf = RMF.create_rmf_file(rmf_file_name)
272 
273 # we want to see the scores of the various restraints also
274 IMP.rmf.add_restraints(rmf, restraints)
275 # and the actual structures
276 IMP.rmf.add_hierarchy(rmf, universe)
277 
278 for i in range(0, configuration_set.get_number_of_configurations()):
279  configuration_set.load_configuration(i)
280  # align all the configurations with the first so they display nicely
281  # if we want to be fancy we can account for flips too
282  if i == 0:
283  base_coords = [IMP.core.XYZ(p).get_coordinates()
284  for p in optimized_particles]
285  else:
286  tr = IMP.algebra.get_transform_taking_first_to_second(
287  optimized_particles, base_coords)
288  IMP.core.transform(optimized_particles, tr)
289  # update the restraint scores
290  sf.evaluate(False)
291  IMP.rmf.save_frame(rmf, i)
292 
293 print "You can now open", rmf_file_name, "in chimera"