IMP logo
IMP Reference Guide  2.21.0
The Integrative Modeling Platform
precision_rmsd.py
1 """@namespace IMP.sampcon.precision_rmsd
2  Calculation of precision and RMSD."""
3 
4 from __future__ import print_function
5 import os
6 import numpy
7 import IMP
8 import IMP.em
9 import pyRMSD.RMSDCalculator
10 from IMP.sampcon.rmsd_calculation import parse_symm_groups_for_pyrmsd
11 
12 
13 def parse_custom_ranges(ranges_file):
14  if not ranges_file:
15  return []
16  with open(ranges_file) as fh:
17  d = {}
18  exec(fh.read(), d)
19  return d['density_custom_ranges']
20 
21 
22 def get_particles_from_superposed(
23  cluster_conform_i, cluster_conform_0, align, ps, trans,
24  symm_groups=None):
25  def _to_vector3ds(numpy_array):
26  # No need to fit the whole array - we only need 4 non-coplanar points,
27  # so 100 should be plenty
28  return [IMP.algebra.Vector3D(c) for c in numpy_array[:100]]
29 
30  if align:
31  calculator_name = "QCP_SERIAL_CALCULATOR"
32  else:
33  calculator_name = "NOSUP_SERIAL_CALCULATOR"
34 
35  conforms = numpy.array([cluster_conform_0, cluster_conform_i])
36 
37  if symm_groups is None:
38  calculator = pyRMSD.RMSDCalculator.RMSDCalculator(
39  calculator_name,
40  conforms)
41  else:
42  if calculator_name == 'NOSUP_SERIAL_CALCULATOR':
43  # calc_symm_groups are enough without any fitting
44  s1 = parse_symm_groups_for_pyrmsd(symm_groups)
45  calculator = pyRMSD.RMSDCalculator.RMSDCalculator(
46  calculator_name,
47  fittingCoordsets=conforms,
48  calculationCoordsets=conforms,
49  calcSymmetryGroups=s1,
50  fitSymmetryGroups=[])
51  else:
52  calculator = pyRMSD.RMSDCalculator.RMSDCalculator(
53  calculator_name,
54  fittingCoordsets=conforms,
55  calcSymmetryGroups=[],
56  fitSymmetryGroups=symm_groups)
57 
58  check1 = (calculator_name == 'NOSUP_SERIAL_CALCULATOR')
59  check2 = symm_groups is not None
60  if check1 and check2:
61  # superposed calc_coords returned if calc-coords
62  # specified while creating the calculator
63  rmsd, superposed_fit, _calc_fit = calculator.pairwise(
64  0, 1, get_superposed_coordinates=True)
65  else:
66  rmsd, superposed_fit = calculator.pairwise(
67  0, 1, get_superposed_coordinates=True)
68  # Get transformation from pyRMSD reference on the first call.
69  # This is somewhat inefficient (since we are essentially repeating
70  # the pyRMSD calculation) but pyRMSD doesn't appear to make its
71  # reference orientation available.
72  if trans is None:
74  _to_vector3ds(superposed_fit[0]), _to_vector3ds(cluster_conform_0))
75 
76  for particle_index in range(len(superposed_fit[1])):
77  # Transform from pyRMSD back to original reference
78  IMP.core.XYZ(ps[particle_index]).set_coordinates(
79  trans * IMP.algebra.Vector3D(superposed_fit[1][particle_index]))
80 
81  return rmsd, ps, trans
82 
83 
84 class GetModelDensity(object):
85  """Compute mean density maps from structures.
86  Keeps a dictionary of density maps,
87  keys are in the custom ranges. When you call add_subunits_density, it adds
88  particle coordinates to the existing density maps.
89  """
90 
91  def __init__(self, custom_ranges=None, resolution=20.0, voxel=5.0,
92  bead_names=None):
93  """Constructor.
94  @param list of particles decorated with mass, radius, and XYZ
95  @param resolution The MRC resolution of the output map
96  (in Angstrom unit)
97  @param voxel The voxel size for the output map (lower is slower)
98  """
99 
100  self.MRCresolution = resolution
101  self.voxel = voxel
102  self.count_models = 0.0
103  self.densities = {}
104  self.bead_names = bead_names
105  self.custom_ranges = custom_ranges
106 
107  # for each custom range get the particle indices that will be
108  # added to the density for that custom range
109  self.particle_indices_in_custom_ranges = {}
110 
111  for density_name in self.custom_ranges:
112  self.particle_indices_in_custom_ranges[density_name] = []
113 
114  # go through each bead, put it in the appropriate custom range(s)
115  for index, beadname in enumerate(self.bead_names):
116  for density_name in self.custom_ranges:
117  # each domain in the list custom_ranges[density_name]
118  for domain in self.custom_ranges[density_name]:
119  if self._is_contained(beadname, domain):
120  self.particle_indices_in_custom_ranges[
121  density_name].append(index)
122  break # already added particle to this custom range
123 
124  def normalize_density(self):
125  pass
126 
127  def _create_density_from_particles(self, ps, name,
128  kernel_type='GAUSSIAN'):
129  '''Internal function for adding to densities.
130  pass XYZR particles with mass and create a density from them.
131  kernel type options are GAUSSIAN, BINARIZED_SPHERE, and SPHERE.'''
132  dmap = IMP.em.SampledDensityMap(ps, self.MRCresolution, self.voxel)
133  dmap.calcRMS()
134  dmap.set_was_used(True)
135 
136  if name not in self.densities:
137  self.densities[name] = dmap
138  else:
139  bbox1 = IMP.em.get_bounding_box(self.densities[name])
140  bbox2 = IMP.em.get_bounding_box(dmap)
141  bbox1 += bbox2
142  dmap3 = IMP.em.create_density_map(bbox1, self.voxel)
143  dmap3.set_was_used(True)
144  dmap3.add(dmap)
145  dmap3.add(self.densities[name])
146  self.densities[name] = dmap3
147 
148  def _is_contained(self, bead_name, domain):
149  """ domain can be the name of a single protein or a tuple
150  (start_residue,end_residue,protein_name)
151  bead is a string of type moleculeName_startResidue_endResidue
152  """
153 
154  (bead_protein, bead_res_start,
155  bead_res_end, bead_copy) = bead_name.rsplit("_", 3)
156 
157  # protein name and copy number check
158  if isinstance(domain, tuple):
159  domain_protein = domain[2]
160  else:
161  domain_protein = domain
162  # A period indicates that we have a copy number
163  if "." in domain_protein:
164  spl = domain_protein.split(".")
165  domain_protein = spl[0]
166  domain_copy = int(spl[1])
167  else:
168  domain_copy = bead_copy = -1
169 
170  if bead_protein != domain_protein or int(bead_copy) != domain_copy:
171  return False
172 
173  # residue range check
174  if isinstance(domain, tuple):
175  bead_residues = set(range(int(bead_res_start),
176  int(bead_res_end)+1))
177  domain_residues = set(range(int(domain[0]),
178  int(domain[1])+1))
179  return not domain_residues.isdisjoint(bead_residues)
180  else:
181  return True
182 
183  def add_subunits_density(self, ps):
184  """Add a frame to the densities.
185  @param ps List of particles decorated with XYZR and Mass.
186  """
187  self.count_models += 1.0
188  # initialize custom list of particles
189  particles_custom_ranges = {}
190  for density_name in self.custom_ranges:
191  particles_custom_ranges[density_name] = []
192 
193  # add each particle to the relevant custom list
194  for density_name in self.custom_ranges:
195  for particle_index \
196  in self.particle_indices_in_custom_ranges[density_name]:
197  particles_custom_ranges[density_name].append(
198  ps[particle_index])
199 
200  # finally, add each custom particle list to the density
201  for density_name in self.custom_ranges:
202  self._create_density_from_particles(
203  particles_custom_ranges[density_name], density_name)
204 
205  def get_density_keys(self):
206  return list(self.densities.keys())
207 
208  def get_density(self, name):
209  """Get the current density for some component name"""
210  if name not in self.densities:
211  return None
212  else:
213  return self.densities[name]
214 
215  def write_mrc(self, path=".", file_prefix=""):
216  densmap = {}
217  for density_name in self.densities:
218  mrc = os.path.join(path, file_prefix + "_" + density_name + ".mrc")
219  self.densities[density_name].multiply(1. / self.count_models)
220  densmap[density_name] = mrc
221  IMP.em.write_map(
222  self.densities[density_name], mrc,
224  return densmap
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:674
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:505
Utilities to help with RMSD calculation.
VectorD< 3 > Vector3D
Definition: VectorD.h:408
Transformation3D get_transformation_aligning_first_to_second(Vector3Ds a, Vector3Ds b)