IMP  2.0.0
The Integrative Modeling Platform
kernel/nup84.py

Modify the Nup84 CG example by replacing a couple of the protein representations with higher resolution representations generated from PDB files. 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 the \ref modules/kernel/examples/nup84_cg.py "Nup84 CG"
3 ## example by replacing a couple of the protein representations with
4 ## higher resolution representations generated from PDB files. In
5 ## addition, show how to visualize restraints and visualize the
6 ## rejected conformations. Both are useful things to do when trying to
7 ## figure out why optimization is not converging.
8 
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 def add_protein_from_length(model, name, residues, parent, restraints,
43  excluded_volume_particles, optimized_particles):
44  ## Create a coarse grained protein with the passed residue information
45  h = IMP.atom.create_protein(model, name, resolution, residues)
46 
47  parent.add_child(h)
48  # Here, each of the constituent particles will be optimized independently
49  leaves = IMP.atom.get_leaves(h)
50  optimized_particles.extend(leaves)
51  excluded_volume_particles.extend(leaves)
52 
53  ## Ensure that the various particles of the protein stay connected
54  r=IMP.atom.create_connectivity_restraint([IMP.atom.Selection(c)\
55  for c in h.get_children()], k,
56  "connect "+name)
57 
58  if r:
59  ## make sure there is an actual restraint
60  restraints.append(r)
61  r.set_maximum_score(k)
62 
63 ## Create protein as a rigid body from a pdb file
64 def add_protein_from_pdb(model, name, file, parent, restraints,
65  excluded_volume_particles, optimized_particles):
66  # we should keep the original particles around so they get written
67 
68  # create an atomic representation from the pdb file
69  t=IMP.atom.read_pdb( IMP.get_example_path(os.path.join("data",file)), model,
70  IMP.atom.ATOMPDBSelector() )
71  # extract the chain from the file
72  c=IMP.atom.Chain(IMP.atom.get_by_type(t, IMP.atom.CHAIN_TYPE)[0])
73  # there is no reason to use all atoms, just approximate the pdb shape instead
74  s=IMP.atom.create_simplified_along_backbone(c, resolution/2.0, True)
75  s.set_name(name)
76  # tear down what is left
77  IMP.atom.destroy(t)
78  # make the simplified structure rigid
79  rb=IMP.atom.create_rigid_body(s)
80  rb.set_coordinates_are_optimized(True)
81  optimized_particles.append(rb)
82  excluded_volume_particles.extend(s.get_children())
83  parent.add_child(s)
84 
85 
86 ## Create protein as a rigid body from several pdb file
87 def add_protein_from_pdbs(model, name, files, parent, restraints,
88  excluded_volume_particles, optimized_particles):
89  h= IMP.atom.Hierarchy.setup_particle(IMP.Particle(model, name))
90  for i, f in enumerate(files):
91  add_protein_from_pdb(model, name+str(i), f, h, restraints,
92  excluded_volume_particles, optimized_particles)
93  r=IMP.atom.create_connectivity_restraint([IMP.atom.Selection(c, hierarchy_types=[IMP.atom.FRAGMENT_TYPE])
94  for c in h.get_children()],
95  k, "connect "+name)
96  if r:
97  restraints.append(r)
98  r.set_maximum_score(k)
99 
100 ## Create all the needed representation using the above functions
101 def create_representation(model):
102  restraints = []
103  optimized_particles = []
104  excluded_volume_particles = []
105  universe=IMP.atom.Hierarchy.setup_particle(IMP.Particle(model, "the universe"))
106 
107  add_protein_from_length(model, "Nup85", 570, universe, restraints,
108  excluded_volume_particles, optimized_particles)
109 
110  # pin the c-terminus
111  ct = IMP.atom.Selection(universe, molecule = "Nup85",
112  hierarchy_types = [IMP.atom.FRAGMENT_TYPE],
113  terminus = IMP.atom.Selection.C)
114  d = IMP.core.XYZ(ct.get_selected_particles()[0])
115  d.set_coordinates(IMP.algebra.Vector3D(0,0,0))
116  d.set_coordinates_are_optimized(False)
117 
118  add_protein_from_length(model, "Nup84", 460, universe, restraints,
119  excluded_volume_particles, optimized_particles)
120  add_protein_from_length(model, "Nup145C", 442, universe, restraints,
121  excluded_volume_particles, optimized_particles)
122  add_protein_from_length(model, "Nup120", [0, 500, 761], universe, restraints,
123  excluded_volume_particles, optimized_particles)
124  add_protein_from_length(model, "Nup133", [0, 450, 778, 1160], universe, restraints,
125  excluded_volume_particles, optimized_particles)
126  add_protein_from_pdb(model, "Seh1", "seh1.pdb", universe, restraints,
127  excluded_volume_particles, optimized_particles)
128  add_protein_from_pdb(model, "Sec13", "sec13.pdb", universe, restraints,
129  excluded_volume_particles, optimized_particles)
130  return universe, restraints, excluded_volume_particles, optimized_particles
131 
132 def add_distance_restraint(selection0, selection1, name, restraints):
133  r=IMP.atom.create_distance_restraint(selection0,selection1, 0, k, name)
134  r.set_maximum_score(k)
135  restraints.append(r)
136 
137 
138 def encode_data_as_restraints(universe, restraints):
139  s0=IMP.atom.Selection(hierarchy=universe, hierarchy_types=[IMP.atom.FRAGMENT_TYPE],
140  molecule="Nup145C", residue_indexes=[(0,423)])
141  s1=IMP.atom.Selection(hierarchy=universe, hierarchy_types=[IMP.atom.FRAGMENT_TYPE],
142  molecule="Nup84")
143  s2=IMP.atom.Selection(hierarchy=universe, hierarchy_types=[IMP.atom.FRAGMENT_TYPE],
144  molecule="Sec13")
145  r=IMP.atom.create_connectivity_restraint([s0,s1,s2], k, "Nup145C Nup84 Sec13")
146  r.set_maximum_score(k)
147  restraints.append(r)
148 
149  add_distance_restraint(IMP.atom.Selection(hierarchy=universe, molecule="Nup145C",
150  residue_indexes=[(0,423)],
151  hierarchy_types=[IMP.atom.FRAGMENT_TYPE]),
152  IMP.atom.Selection(hierarchy=universe, molecule="Nup85",
153  hierarchy_types=[IMP.atom.FRAGMENT_TYPE]),
154  "Num145C, Nup85", restraints)
155  add_distance_restraint(IMP.atom.Selection(hierarchy=universe, molecule="Nup145C",
156  residue_indexes=[(0,423)],
157  hierarchy_types=[IMP.atom.FRAGMENT_TYPE]),
158  IMP.atom.Selection(hierarchy=universe, molecule="Nup120",
159  residue_indexes= [(500, 762)],
160  hierarchy_types=[IMP.atom.FRAGMENT_TYPE]),
161  "Nup145C Nup120", restraints)
162  add_distance_restraint(IMP.atom.Selection(hierarchy=universe, molecule="Nup84",
163  hierarchy_types=[IMP.atom.FRAGMENT_TYPE]),
164  IMP.atom.Selection(hierarchy=universe, molecule="Nup133",
165  residue_indexes=[(778, 1160)],
166  hierarchy_types=[IMP.atom.FRAGMENT_TYPE]),
167  "Nup84 Nup133", restraints)
168  add_distance_restraint(IMP.atom.Selection(hierarchy=universe, molecule="Nup85",
169  hierarchy_types=[IMP.atom.FRAGMENT_TYPE]),
170  IMP.atom.Selection(hierarchy=universe, molecule="Seh1",
171  hierarchy_types=[IMP.atom.FRAGMENT_TYPE]),
172  "Nup85 Seh1", restraints)
173  add_distance_restraint(IMP.atom.Selection(hierarchy=universe, molecule="Nup145C",
174  residue_indexes=[(0,423)],
175  hierarchy_types=[IMP.atom.FRAGMENT_TYPE]),
176  IMP.atom.Selection(hierarchy=universe, molecule="Sec13",
177  hierarchy_types=[IMP.atom.FRAGMENT_TYPE]),
178  "Nup145C Sec13", restraints)
179 
180 
181 
182 # find acceptable conformations of the model
183 def get_configurations(model, restraints, excluded_volume_particles, optimized_particles):
184  #cpc= IMP.container.ClosePairContainer(representation.get_particles(), 0, 10)
185  #evr= IMP.container.PairRestraint(IMP.core.SoftSpherePairScore(k), cpc,
186  # "Excluded Volume")
187  scale=.5
188  mc= IMP.core.MonteCarlo(model)
189  movers=[]
190  for p in optimized_particles:
191  if IMP.core.RigidBody.particle_is_instance(p):
192  mover= IMP.core.RigidBodyMover(p, IMP.core.XYZR(p).get_radius()*scale,
193  .2*scale)
194  movers.append(mover)
195  else:
196  mover= IMP.core.BallMover([p], IMP.core.XYZR(p).get_radius()*scale)
197  movers.append(mover)
198  serial_mover= IMP.core.SerialMover(movers)
199  mc.add_mover(serial_mover)
200  scoring_function= IMP.core.IncrementalScoringFunction(optimized_particles, restraints)
201  scoring_function.add_close_pair_score(IMP.core.SoftSpherePairScore(k), 0.0,
202  excluded_volume_particles)
203 
204  configuration_set= IMP.ConfigurationSet(model)
205  # must write our own sampler as IMP.core.MCCGSampler doesn't handle rigid bodies
206  for i in range(number_of_sampling_attempts):
207  for p in optimized_particles:
208  IMP.core.XYZ(p).set_coordinates(IMP.algebra.get_random_vector_in(bb))
209  mc.optimize(number_of_mc_steps)
210  if scoring_function.get_had_good_score():
211  configuration_set.save()
212  return configuration_set
213 
214 
215 model= IMP.Model()
216 universe, restraints, excluded_volume_particles, optimized_particles= create_representation(model)
217 encode_data_as_restraints(universe, restraints)
218 
219 configuration_set= get_configurations(model, restraints,
220  excluded_volume_particles,
221  optimized_particles)
222 
223 print "Found", configuration_set.get_number_of_configurations(), "good configurations"
224 
225 # now lets save them all to a file
226 rmf_file_name=IMP.base.create_temporary_file_name("nup84", ".rmf")
227 rmf= RMF.create_rmf_file(rmf_file_name)
228 
229 # we want to see the scores of the various restraints also
230 IMP.rmf.add_restraints(rmf, restraints)
231 # and the actual structures
232 IMP.rmf.add_hierarchy(rmf, universe)
233 
234 for i in range(0, configuration_set.get_number_of_configurations()):
235  configuration_set.load_configuration(i)
236  # align all the configurations with the first so they display nicely
237  # if we want to be fancy we can account for flips too
238  if i==0:
239  base_coords=[IMP.core.XYZ(p).get_coordinates() for p in optimized_particles]
240  else:
241  tr= IMP.algebra.get_transform_taking_first_to_second(optimized_particles, base_coords)
242  IMP.core.transform(optimized_particles, tr)
243  # update the restraint scores
244  sf.evaluate(False)
245  IMP.rmf.save_frame(rmf, i)
246 
247 print "You can now open", rmf_file_name, "in chimera"