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