1 """@namespace IMP.pmi.restraints.em
2 Restraints for handling electron microscopy maps.
5 from __future__
import print_function
22 """Fit Gaussian-decorated particles to an EM map
23 (also represented with a set of Gaussians)
24 @note This class wraps an isd::GaussianEMRestraint
29 cutoff_dist_model_model=0.0,
30 cutoff_dist_model_data=0.0,
31 target_mass_scale=1.0,
33 target_radii_scale=3.0,
34 model_radii_scale=1.0,
36 spherical_gaussians=
False,
37 close_pair_container=
None,
39 scale_target_to_mass=
False,
41 target_is_rigid_body=
False,
45 @param densities The Gaussian-decorated particles to be restrained
46 @param target_fn GMM file of the target density map
47 (alternatively, pass the ps)
48 @param target_ps List of Gaussians of the target map
49 (alternatively, pass the filename)
50 @param cutoff_dist_model_model Distance in model-model close
52 @param cutoff_dist_model_data Distance in model-data close pair
53 container. Usually can set to zero because we multiply the
55 @param target_mass_scale Scale up the target densities so that
57 Needed if the GMM you generated was not already scaled.
58 To make it the same as model mass, set scale_to_target_mass=True
59 @param target_mass Sets the mass of the target density to the given value. Default is None. This
60 will override target_mass_scale argument
61 @param target_radii_scale Scale the target density radii -
62 only used for the close pair container.
63 If you keep this at 3.0 or so you don't have to use cutoff dist.
64 @param model_radii_scale Scale the model density radii - only used
65 for the close pair container
66 @param slope Linear term added to help bring model into the density
67 @param spherical_gaussians Set to True for a speed bonus when
68 the model densities are spheres. (This means you don't have
69 to do a matrix multiplication if they rotate.)
70 @param close_pair_container Pass a close pair container for
71 the model if you already have one (e.g. for an excluded
72 volume restraint.) May give a speed bonus.
73 @param backbone_slope Only apply slope to backbone particles -
74 only matters for atomic
75 @param scale_target_to_mass Set True if you would need to scale
76 target to EXACTLY the model mass
77 @param weight The restraint weight
78 @param target_is_rigid_body Set True if you want to put the target density particles
79 into a rigid body that need to be sampled (e.g.,when you need to fit one density
80 against another one). Default is False.
81 @param local Only consider density particles that are within the
82 specified model-density cutoff (experimental)
87 self.sigmaissampled =
False
88 self.sigmamaxtrans = 0.3
93 self.densities = densities
94 self.em_root_hier =
None
97 self.m = self.densities[0].get_model()
99 if scale_target_to_mass:
100 def hierarchy_mass(h):
103 target_mass = sum(hierarchy_mass(h)
for h
in densities)
104 print(
'will set target mass to', target_mass)
105 print(
'will scale target mass by', target_mass_scale)
108 self._set_dataset(target_fn, representation)
114 radius_scale=target_radii_scale,
115 mass_scale=target_mass_scale)
116 elif target_ps != []:
117 self.target_ps = target_ps
119 print(
'Gaussian EM restraint: must provide target density file or properly set up target densities')
124 scale=target_mass/tmass
125 for p
in self.target_ps:
130 for p, state
in representation._protocol_output:
131 p.add_em3d_restraint(state, self.target_ps, self.densities,
136 for h
in self.densities:
138 if model_radii_scale != 1.0:
139 for p
in self.model_ps:
147 if target_is_rigid_body:
155 self.sigmaglobal = IMP.pmi.tools.SetupNuisance(self.m, self.sigmainit,
156 self.sigmamin, self.sigmamax,
157 self.sigmaissampled).get_particle()
160 print(
'target num particles', len(self.target_ps), \
162 for p
in self.target_ps]))
163 print(
'model num particles', len(self.model_ps), \
165 for p
in self.model_ps]))
167 update_model=
not spherical_gaussians
173 self.sigmaglobal.get_particle().
get_index(),
174 cutoff_dist_model_model,
175 cutoff_dist_model_data,
177 update_model, backbone_slope, local)
179 self.gaussianEM_restraint.set_density_filename(target_fn)
181 print(
'done EM setup')
183 self.rs.add_restraint(self.gaussianEM_restraint)
184 self.set_weight(weight)
186 def _set_dataset(self, target_fn, representation):
187 """Set the dataset to point to the input file"""
189 self.dataset = representation.get_file_dataset(target_fn)
192 l = ihm.location.InputFileLocation(target_fn,
193 details=
"Electron microscopy density map, "
194 "represented as a Gaussian Mixture "
196 self.dataset = ihm.dataset.EMDensityDataset(l)
198 m = re.match(
'(.*\.mrc)\..*\.txt$', target_fn)
199 if m
and os.path.exists(m.group(1)):
200 l = ihm.location.InputFileLocation(path=m.group(1),
201 details=
'Original MRC file from which the GMM was derived')
202 self.dataset.parents.append(ihm.dataset.EMDensityDataset(l))
204 def center_target_density_on_model(self):
207 for p
in self.target_ps:
210 target_com += pos * mass
212 target_com /= target_mass
213 print(
'target com', target_com)
216 for p
in self.model_ps:
219 model_com += pos * mass
221 model_com /= model_mass
222 print(
'model com', model_com)
224 v = target_com - model_com
225 print(
'translating with', -v)
227 for p
in self.target_ps:
232 '''Returns the geometric center of the GMM particles
233 @param target = True - returns target map gmm COM
234 @param target = False - returns model gmm COM'''
249 def center_target_density_on_origin(self):
251 print(
'target com', target_com)
253 print(
'model com', model_com)
254 v = target_com - model_com
255 print(
'translating with', -v)
257 for p
in self.target_ps:
261 def center_model_on_target_density(self, input_object):
263 hier = input_object.prot
265 hier = input_object.get_hierarchy()
267 raise Exception(
"Input must be a Representation or topology.State object")
269 print(
'target com', target_com)
271 print(
'model com', model_com)
272 v = target_com - model_com
273 print(
'translating with', v)
286 for rb
in list(rigid_bodies):
289 for p
in list(XYZRs):
292 def center_on_target_density(self):
294 print(
'target com', target_com)
296 print(
'model com', model_com)
297 v = target_com - model_com
298 print(
'translating with', v)
304 for p
in self.model_ps:
311 for rb
in list(rigid_bodies):
314 for p
in list(XYZRs):
320 def set_weight(self,weight):
322 self.rs.set_weight(weight)
324 def set_label(self, label):
327 def add_to_model(self):
330 def get_particles_to_sample(self):
332 if self.sigmaissampled:
333 ps[
"Nuisances_GaussianEMRestraint_sigma_" +
334 self.label] = ([self.sigmaglobal], self.sigmamaxtrans)
336 ps[
"Rigid_Bodies_GaussianEMRestraint"] = (
342 def get_rigid_body(self):
344 raise Exception(
"No rigid body created for GMM particles. Ensure target_is_rigid_body is set to True")
347 def get_density_as_hierarchy(self):
348 if self.em_root_hier
is None:
350 self.em_root_hier.set_name(
"GaussianEMRestraint_density_"+self.label)
351 for p
in self.target_ps:
352 self.em_root_hier.add_child(p)
353 return self.em_root_hier
356 ''' Can add a target GMM to a Hierarchy.
357 For PMI2 a state object may also be passed'''
359 inp.get_hierarchy().add_child(self.get_density_as_hierarchy())
361 inp.add_child(self.get_density_as_hierarchy())
363 raise Exception(
"Can only add a density to a PMI State object or IMP.atom.Hierarchy. You passed a", type(inp))
365 def get_restraint(self):
368 def get_restraint_set(self):
371 def get_output(self):
373 score = self.weight * self.rs.unprotected_evaluate(
None)
374 ccc = self.gaussianEM_restraint.get_cross_correlation_coefficient()
376 output[
"_TotalScore"] = str(score)
377 output[
"GaussianEMRestraint_" +
378 self.label] = str(score)
379 output[
"GaussianEMRestraint_%s_CCC" % self.label] = ccc
380 output[
"GaussianEMRestraint_sigma_" +
381 self.label] = str(self.sigmaglobal.get_scale())
385 return self.weight * self.rs.unprotected_evaluate(
None)
388 '''Writes target GMM file to MRC'''
390 fileout=
"Gaussian_map_" + self.label +
".mrc"
398 """Fit particles to an EM map. This creates a simulate density map and updates them every eval.
399 @note Wraps an em::FitRestraint
410 @param ps The particles to restrain. Currently these must be atomic particles.
411 @param map_fn The EM density map to fit to
412 @param resolution Map resolution
413 @param origin In case you need to tell IMP the correct origin
414 @param voxel_size In case you need to tell IMP the angstroms per pixel
415 @param weight The data weight
416 @param label Extra PMI label
418 print(
'FitRestraint: setup')
420 print(
'\tresolution',resolution)
421 print(
'\tvoxel_size',voxel_size)
422 print(
'\torigin',origin)
423 print(
'\tweight',weight)
426 self.mdl = ps[0].get_model()
432 self.dmap.update_voxel_size(voxel_size)
433 if origin
is not None:
434 if type(origin)==IMP.algebra.Vector3D:
435 self.dmap.set_origin(origin)
436 elif type(origin)==list:
437 self.dmap.set_origin(*origin)
439 print(
'FitRestraint did not recognize format of origin')
443 self.rs.add_restraint(fr)
444 self.set_weight(weight)
446 def set_weight(self,weight):
448 self.rs.set_weight(weight)
450 def set_label(self, label):
453 def add_to_model(self):
456 def get_restraint_set(self):
459 def get_output(self):
461 score = self.weight * self.rs.unprotected_evaluate(
None)
462 output[
"_TotalScore"] = str(score)
463 output[
"EMRestraint_" + self.label] = str(score)
467 return self.weight * self.rs.unprotected_evaluate(
None)
472 class ElectronMicroscopy2D(object):
481 self.m = representation.prot.get_model()
488 resolution=resolution)
491 self.rs.add_restraint(em2d)
493 def set_label(self, label):
496 def add_to_model(self):
499 def get_restraint(self):
502 def set_weight(self,weight):
504 self.rs.set_weigth(self.weight)
506 def get_output(self):
508 score = self.weight*self.rs.unprotected_evaluate(
None)
509 output[
"_TotalScore"] = str(score)
510 output[
"ElectronMicroscopy2D_" + self.label] = str(score)
Fit Gaussian-decorated particles to an EM map (also represented with a set of Gaussians) ...
A member of a rigid body, it has internal (local) coordinates.
static bool get_is_setup(const IMP::ParticleAdaptor &p)
def write_target_gmm_to_mrc
Writes target GMM file to MRC.
Various classes to hold sets of particles.
static XYZR setup_particle(Model *m, ParticleIndex pi)
double get_mass(ResidueType c)
Get the mass from the residue type.
def add_target_density_to_hierarchy
Can add a target GMM to a Hierarchy.
static bool get_is_setup(const IMP::ParticleAdaptor &p)
Object used to hold a set of restraints.
ParticleIndexPairs get_indexes(const ParticlePairsTemp &ps)
The standard decorator for manipulating molecular structures.
Ints get_index(const ParticlesTemp &particles, const Subset &subset, const Subsets &excluded)
void transform(XYZ a, const algebra::Transformation3D &tr)
Apply a transformation to the particle.
Restraint between two Gaussian Mixture Models, "model" and "density".
A decorator for a particle with x,y,z coordinates.
def get_center_of_mass
Returns the geometric center of the GMM particles.
static Copy setup_particle(Model *m, ParticleIndex pi, Int number)
Basic functionality that is expected to be used by a wide variety of IMP users.
General purpose algebraic and geometric methods that are expected to be used by a wide variety of IMP...
Calculate score based on fit to EM map.
The general base class for IMP exceptions.
IMP::core::RigidBody create_rigid_body(Hierarchy h)
Class to handle individual particles of a Model object.
Stores a list of Molecules all with the same State index.
A decorator for a rigid body.
Set up the representation of all proteins and nucleic acid macromolecules.
Functionality for loading, creating, manipulating and scoring atomic structures.
Hierarchies get_leaves(const Selection &h)
Fit particles to an EM map.
Inferential scoring building on methods developed as part of the Inferential Structure Determination ...
A decorator for a particle with x,y,z coordinates and a radius.