IMP logo
IMP Reference Guide  2.15.0
The Integrative Modeling Platform
precision_rmsd.py
1 from __future__ import print_function
2 import os
3 import numpy
4 import IMP
5 import IMP.em
6 import pyRMSD.RMSDCalculator
7 
8 
9 def parse_custom_ranges(ranges_file):
10  with open(ranges_file) as fh:
11  d = {}
12  exec(fh.read(), d)
13  return d['density_custom_ranges']
14 
15 
16 def get_particles_from_superposed(
17  cluster_conform_i, cluster_conform_0, align, ps, trans):
18  def _to_vector3ds(numpy_array):
19  # No need to fit the whole array - we only need 4 non-coplanar points,
20  # so 100 should be plenty
21  return [IMP.algebra.Vector3D(c) for c in numpy_array[:100]]
22 
23  if align:
24  calculator = pyRMSD.RMSDCalculator.RMSDCalculator(
25  "QCP_SERIAL_CALCULATOR",
26  numpy.array([cluster_conform_0, cluster_conform_i]))
27  else:
28  calculator = pyRMSD.RMSDCalculator.RMSDCalculator(
29  "NOSUP_SERIAL_CALCULATOR",
30  numpy.array([cluster_conform_0, cluster_conform_i]))
31 
32  rmsd, superposed_fit = calculator.pairwise(
33  0, 1, get_superposed_coordinates=True)
34  # Get transformation from pyRMSD reference on the first call.
35  # This is somewhat inefficient (since we are essentially repeating
36  # the pyRMSD calculation) but pyRMSD doesn't appear to make its
37  # reference orientation available.
38  if trans is None:
40  _to_vector3ds(superposed_fit[0]), _to_vector3ds(cluster_conform_0))
41 
42  for particle_index in range(len(superposed_fit[1])):
43  # Transform from pyRMSD back to original reference
44  IMP.core.XYZ(ps[particle_index]).set_coordinates(
45  trans * IMP.algebra.Vector3D(superposed_fit[1][particle_index]))
46 
47  return rmsd, ps, trans
48 
49 
50 def get_particles_from_superposed_amb(
51  cluster_conform_i, cluster_conform_0, align, ps, trans, symm_groups):
52 
53  '''Modified superposed function to work with symmetric copies'''
54 
55  def _to_vector3ds(numpy_array):
56  # No need to fit the whole array - we only need 4 non-coplanar points,
57  # so 100 should be plenty
58  return [IMP.algebra.Vector3D(c) for c in numpy_array[:100]]
59 
60  min_rmsd = 10000.0
61 
62  superposed_final_coords = []
63 
64  for perm in pyRMSD.symmTools.symm_permutations(symm_groups):
65  # for each permutation
66 
67  new_cluster_conform_i = cluster_conform_i
68 
69  for sg in perm: # for each symmetric group in perm
70 
71  for [particle0, particle1] in sg:
72  # swap the particles if they are in non-standard order
73  # in this permutation
74  if particle0 > particle1:
75  pyRMSD.symmTools.swap_atoms(
76  new_cluster_conform_i, particle0, particle1)
77 
78  if align:
79  calculator = pyRMSD.RMSDCalculator.RMSDCalculator(
80  "QCP_SERIAL_CALCULATOR",
81  numpy.array([cluster_conform_0, new_cluster_conform_i]))
82  else:
83  calculator = pyRMSD.RMSDCalculator.RMSDCalculator(
84  "NOSUP_SERIAL_CALCULATOR",
85  numpy.array([cluster_conform_0, new_cluster_conform_i]))
86 
87  rmsd, superposed_fit = calculator.pairwise(
88  0, 1, get_superposed_coordinates=True)
89 
90  if rmsd < min_rmsd:
91  min_rmsd = rmsd
92  superposed_final_coords = superposed_fit
93 
94  # Get transformation from pyRMSD reference on the first call.
95  # This is somewhat inefficient (since we are essentially repeating
96  # the pyRMSD calculation) but pyRMSD doesn't appear to make its
97  # reference orientation available.
98  if trans is None:
100  _to_vector3ds(superposed_final_coords[0]),
101  _to_vector3ds(cluster_conform_0))
102 
103  for particle_index in range(len(superposed_final_coords[1])):
104 
105  # Transform from pyRMSD back to original reference
106  IMP.core.XYZ(ps[particle_index]).set_coordinates(
107  trans * IMP.algebra.Vector3D(
108  superposed_final_coords[1][particle_index]))
109 
110  return min_rmsd, ps, trans
111 
112 
113 class GetModelDensity(object):
114  """Compute mean density maps from structures.
115  Keeps a dictionary of density maps,
116  keys are in the custom ranges. When you call add_subunits_density, it adds
117  particle coordinates to the existing density maps.
118  """
119 
120  def __init__(self, custom_ranges=None, resolution=20.0, voxel=5.0,
121  bead_names=None):
122  """Constructor.
123  @param list of particles decorated with mass, radius, and XYZ
124  @param resolution The MRC resolution of the output map
125  (in Angstrom unit)
126  @param voxel The voxel size for the output map (lower is slower)
127  """
128 
129  self.MRCresolution = resolution
130  self.voxel = voxel
131  self.count_models = 0.0
132  self.densities = {}
133  self.bead_names = bead_names
134  self.custom_ranges = custom_ranges
135 
136  # for each custom range get the particle indices that will be
137  # added to the density for that custom range
138  self.particle_indices_in_custom_ranges = {}
139 
140  for density_name in self.custom_ranges:
141  self.particle_indices_in_custom_ranges[density_name] = []
142 
143  # go through each bead, put it in the appropriate custom range(s)
144  for index, beadname in enumerate(self.bead_names):
145  for density_name in self.custom_ranges:
146  # each domain in the list custom_ranges[density_name]
147  for domain in self.custom_ranges[density_name]:
148  if self._is_contained(beadname, domain):
149  self.particle_indices_in_custom_ranges[
150  density_name].append(index)
151  break # already added particle to this custom range
152 
153  def normalize_density(self):
154  pass
155 
156  def _create_density_from_particles(self, ps, name,
157  kernel_type='GAUSSIAN'):
158  '''Internal function for adding to densities.
159  pass XYZR particles with mass and create a density from them.
160  kernel type options are GAUSSIAN, BINARIZED_SPHERE, and SPHERE.'''
161  dmap = IMP.em.SampledDensityMap(ps, self.MRCresolution, self.voxel)
162  dmap.calcRMS()
163  dmap.set_was_used(True)
164 
165  if name not in self.densities:
166  self.densities[name] = dmap
167  else:
168  bbox1 = IMP.em.get_bounding_box(self.densities[name])
169  bbox2 = IMP.em.get_bounding_box(dmap)
170  bbox1 += bbox2
171  dmap3 = IMP.em.create_density_map(bbox1, self.voxel)
172  dmap3.set_was_used(True)
173  dmap3.add(dmap)
174  dmap3.add(self.densities[name])
175  self.densities[name] = dmap3
176 
177  def _is_contained(self, bead_name, domain):
178  """ domain can be the name of a single protein or a tuple
179  (start_residue,end_residue,protein_name)
180  bead is a string of type moleculeName_startResidue_endResidue
181  """
182 
183  (bead_protein, bead_res_start,
184  bead_res_end, bead_copy) = bead_name.split("_")
185 
186  # protein name and copy number check
187  if isinstance(domain, tuple):
188  domain_protein = domain[2]
189  else:
190  domain_protein = domain
191  # A period indicates that we have a copy number
192  if "." in domain_protein:
193  spl = domain_protein.split(".")
194  domain_protein = spl[0]
195  domain_copy = int(spl[1])
196  else:
197  domain_copy = bead_copy = -1
198 
199  if bead_protein != domain_protein or int(bead_copy) != domain_copy:
200  return False
201 
202  # residue range check
203  if isinstance(domain, tuple):
204  bead_residues = set(range(int(bead_res_start),
205  int(bead_res_end)+1))
206  domain_residues = set(range(int(domain[0]),
207  int(domain[1])+1))
208  return not domain_residues.isdisjoint(bead_residues)
209  else:
210  return True
211 
212  def add_subunits_density(self, ps):
213  """Add a frame to the densities.
214  @param ps List of particles decorated with XYZR and Mass.
215  """
216  self.count_models += 1.0
217  # initialize custom list of particles
218  particles_custom_ranges = {}
219  for density_name in self.custom_ranges:
220  particles_custom_ranges[density_name] = []
221 
222  # add each particle to the relevant custom list
223  for density_name in self.custom_ranges:
224  for particle_index \
225  in self.particle_indices_in_custom_ranges[density_name]:
226  particles_custom_ranges[density_name].append(
227  ps[particle_index])
228 
229  # finally, add each custom particle list to the density
230  for density_name in self.custom_ranges:
231  self._create_density_from_particles(
232  particles_custom_ranges[density_name], density_name)
233 
234  def get_density_keys(self):
235  return list(self.densities.keys())
236 
237  def get_density(self, name):
238  """Get the current density for some component name"""
239  if name not in self.densities:
240  return None
241  else:
242  return self.densities[name]
243 
244  def write_mrc(self, path=".", file_prefix=""):
245  for density_name in self.densities:
246  mrc = os.path.join(path, file_prefix + "_" + density_name + ".mrc")
247  self.densities[density_name].multiply(1. / self.count_models)
248  IMP.em.write_map(
249  self.densities[density_name], mrc,
251  if len(self.densities) == 1:
252  return mrc
253  else:
254  return os.path.join(path, file_prefix + "_*.mrc")
def get_density
Get the current density for some component name.
Compute mean density maps from structures.
Class for sampling a density map from particles.
DensityMap * multiply(const DensityMap *m1, const DensityMap *m2)
Return a density map for which voxel i contains the result of m1[i]*m2[i].
DensityMap * create_density_map(const IMP::algebra::GridD< 3, S, V, E > &arg)
Create a density map from an arbitrary IMP::algebra::GridD.
Definition: DensityMap.h:650
def add_subunits_density
Add a frame to the densities.
Basic utilities for handling cryo-electron microscopy 3D density maps.
A decorator for a particle with x,y,z coordinates.
Definition: XYZ.h:30
algebra::BoundingBoxD< 3 > get_bounding_box(const DensityMap *m)
Definition: DensityMap.h:481
VectorD< 3 > Vector3D
Definition: VectorD.h:421
Transformation3D get_transformation_aligning_first_to_second(Vector3Ds a, Vector3Ds b)