IMP  2.0.0
The Integrative Modeling Platform
nup84.py
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"