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