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