IMP  2.2.0
The Integrative Modeling Platform
nup84.py
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
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
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
73  IMP.get_example_path(os.path.join("data", file)), model,
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
83  # make the simplified structure rigid
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):
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 = []
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")
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]),
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]),
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]),
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]),
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]),
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:
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(
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"
A class to store a set of configurations of a model.
A Monte Carlo optimizer.
Definition: MonteCarlo.h:45
void save_frame(RMF::FileHandle file, unsigned int, std::string name="")
Definition: frames.h:42
void add_restraints(RMF::NodeHandle fh, const kernel::Restraints &hs)
See IMP.container for more information.
bool get_bool_flag(std::string name)
std::string create_temporary_file_name(std::string prefix="imp_temp", std::string suffix="")
Modify the transformation of a rigid body.
static bool get_is_setup(const IMP::kernel::ParticleAdaptor &p)
Definition: rigid_bodies.h:157
kernel::Restraint * create_distance_restraint(const Selection &n0, const Selection &n1, double x0, double k, std::string name="Distance%1%")
Modify a set of continuous variables by perturbing them within a ball.
Vector3D get_random_vector_in(const Cylinder3D &c)
Generate a random vector in a cylinder with uniform density.
See IMP.statistics for more information.
Select all non-alternative ATOM records.
Definition: pdb.h:63
Hierarchy create_protein(kernel::Model *m, std::string name, double target_radius, const Ints domain_boundaries)
void add_hierarchy(RMF::FileHandle fh, atom::Hierarchy hs)
void add_bool_flag(std::string name, std::string description)
static Hierarchy setup_particle(kernel::Model *m, kernel::ParticleIndex pi, kernel::ParticleIndexesAdaptor children=kernel::ParticleIndexesAdaptor())
void transform(XYZ a, const algebra::Transformation3D &tr)
Apply a transformation to the particle.
Hierarchies get_by_type(Hierarchy mhd, GetByType t)
kernel::Restraint * create_connectivity_restraint(const Selections &s, double x0, double k, std::string name="Connectivity%1%")
A decorator for a particle with x,y,z coordinates.
Definition: XYZ.h:30
Class to handle individual model particles.
void destroy(Hierarchy d)
Delete the Hierarchy.
See IMP.example for more information.
Hierarchy create_simplified_along_backbone(Chain input, const IntRanges &residue_segments, bool keep_detailed=false)
VectorD< 3 > Vector3D
Definition: VectorD.h:395
IMP::core::RigidBody create_rigid_body(Hierarchy h)
Store info for a chain of a protein.
Definition: Chain.h:21
Strings setup_from_argv(const Strings &argv, std::string description, std::string positional_description, int num_positional)
See IMP.display for more information.
Definition: ChimeraWriter.h:17
See IMP.atom for more information.
void read_pdb(base::TextInput input, int model, Hierarchy h)
Hierarchies get_leaves(const Selection &h)
Apply a list of movers one at a time.
Definition: SerialMover.h:23
Class for storing model, its restraints, constraints, and particles.
Definition: kernel/Model.h:72
A decorator for a particle with x,y,z coordinates and a radius.
Definition: XYZR.h:27