IMP logo
IMP Reference Guide  2.20.0
The Integrative Modeling Platform
comparisons.py
1 """@namespace IMP.EMageFit.imp_general.comparisons
2  Utility functions for comparisons.
3 """
4 
5 
6 import IMP
7 import IMP.algebra
8 import IMP.core
9 import IMP.em
10 import IMP.atom
11 import IMP.EMageFit.imp_general.representation as representation
12 
13 import math
14 import logging
15 log = logging.getLogger("comparisons")
16 
17 
18 def get_coordinates(hierarchy):
19  xyz = [IMP.core.XYZ(leaf) for leaf in IMP.atom.get_leaves(hierarchy)]
20  coords = [x.get_coordinates() for x in xyz]
21  return coords
22 
23 
24 def get_assembly_placement_score(assembly, native_assembly, align=False):
25  """
26  Computes the placement score of an assembly respect to the
27  native_assembly.
28  @param assembly An IMP.atom.Molecule object
29  @param native_assembly An IMP.atom.Molecule object
30  @param align If True, assembly is aligned to native_assembly before
31  calculating the placement score
32  """
33 
34  distances, angles = get_components_placement_scores(assembly,
35  native_assembly, align)
36  n = 1. * len(distances)
37  return sum(distances) / n, sum(angles) / n
38 
39 
40 def get_components_placement_scores(assembly, native_assembly, align=False):
41  """
42  Compute the placement score of each of the children of an assembly.
43  The function does not do any time of alignment of the coordinates
44  @param assembly An IMP.atom.Molecule object
45  @param native_assembly An IMP.atom.Molecule object with the native
46  conformation. Obviously the atoms in assembly
47  and native assembly must be the same
48  @param align If True, assembly is aligned to native_assembly before
49  calculating the placement score
50  @return The function returns 2 lists. The first list contains the
51  placement distances of the children. The second list contains
52  the placement angles
53  """
54  model_coords_per_child = [get_coordinates(c)
55  for c in assembly.get_children()]
56  native_coords_per_child = [get_coordinates(c)
57  for c in native_assembly.get_children()]
58  if align:
59  model_coords = []
60  for x in model_coords_per_child:
61  model_coords.extend(x)
62  native_coords = []
63  for x in native_coords_per_child:
64  native_coords.extend(x)
66  model_coords, native_coords)
67  # get aligned coordinates
68  new_model_coords_per_child = []
69  for c in model_coords_per_child:
70  coords = [T.get_transformed(x) for x in c]
71  new_model_coords_per_child.append(coords)
72  model_coords_per_child = new_model_coords_per_child
73  distances, angles = get_placement_scores_from_coordinates(
74  native_coords_per_child, model_coords_per_child)
75  return distances, angles
76 
77 
78 def get_placement_scores_from_coordinates(model_components_coords,
79  native_components_coords):
80  """
81  Computes the placement score for each of the components
82  @param model_components_coords A list with the coordinates for each
83  component
84  @param native_components_coords A list with the coordinates for each
85  component in the native assembly
86  """
87  distances = []
88  angles = []
89  for model_coords, native_coords in zip(
90  model_components_coords, native_components_coords):
91  distance, angle = get_placement_score_from_coordinates(model_coords,
92  native_coords)
93  distances.append(distance)
94  angles.append(angle)
95  return distances, angles
96 
97 
98 def get_placement_score_from_coordinates(model_coords, native_coords):
99  """
100  Computes the position error (placement distance) and the orientation
101  error (placement angle) of the coordinates in model_coords respect to
102  the coordinates in native_coords.
103  placement distance - translation between the centroids of the
104  coordinates
105  placement angle - Angle in the axis-angle formulation of the rotation
106  aligning the two rigid bodies.
107  """
108  native_centroid = IMP.algebra.get_centroid(native_coords)
109  model_centroid = IMP.algebra.get_centroid(model_coords)
110  translation_vector = native_centroid - model_centroid
111  distance = translation_vector.get_magnitude()
112  if len(model_coords) != len(native_coords):
113  raise ValueError(
114  "Mismatch in the number of members %d %d "
115  % (len(model_coords), len(native_coords)))
117  native_coords)
118  P = IMP.algebra.get_axis_and_angle(TT.get_rotation())
119  angle = P.second
120  return distance, angle
121 
122 
123 def get_rmsd(hierarchy1, hierarchy2):
124  xyz1 = [IMP.core.XYZ(leaf) for leaf in IMP.atom.get_leaves(hierarchy1)]
125  xyz2 = [IMP.core.XYZ(leaf) for leaf in IMP.atom.get_leaves(hierarchy2)]
126  return IMP.atom.get_rmsd(xyz1, xyz2)
127 
128 
129 def get_ccc(native_assembly, assembly, resolution, voxel_size,
130  threshold, write_maps=False):
131  """
132  Threshold - threshold used for the map of the native assembly. Pixels
133  with values above this threshold in the native map are used for the
134  calculation of the cross_correlation_coefficient
135  """
136  particles_native = IMP.atom.get_leaves(native_assembly)
137  particles_solution = IMP.atom.get_leaves(assembly)
138  bb_native = IMP.core.get_bounding_box(IMP.core.XYZs(particles_native))
139  bb_solution = IMP.core.get_bounding_box(IMP.core.XYZs(particles_solution))
140  # bounding box enclosing both the particles of the native assembly
141  # and the particles of the model
142  bb_union = IMP.algebra.get_union(bb_native, bb_solution)
143  # add border of 4 voxels
144  border = 4 * voxel_size
145  bottom = bb_union.get_corner(0)
146  bottom += IMP.algebra.Vector3D(-border, -border, -border)
147  top = bb_union.get_corner(1)
148  top += IMP.algebra.Vector3D(border, border, border)
149  bb_union = IMP.algebra.BoundingBox3D(bottom, top)
150 
151  mrw = IMP.em.MRCReaderWriter()
152  header = IMP.em.create_density_header(bb_union, voxel_size)
153  header.set_resolution(resolution)
154 
155  map_native = IMP.em.SampledDensityMap(header)
156  map_native.set_particles(particles_native)
157  map_native.resample()
158 
159  map_solution = IMP.em.SampledDensityMap(header)
160  map_solution.set_particles(particles_solution)
161  map_solution.resample()
162 
163  if write_maps:
164  IMP.em.write_map(map_solution, "map_solution.mrc", mrw)
165  IMP.em.write_map(map_native, "map_native.mrc", mrw)
166  map_native.calcRMS()
167  map_solution.calcRMS()
168  # base the calculation of the cross_correlation coefficient on the
169  # threshold for the native map, because the threshold for the map of
170  # the model changes with each model
171  threshold = 0.25 # threshold AFTER normalization using calcRMS()
172  ccc = IMP.em.get_coarse_cc_coefficient(map_solution,
173  map_native, threshold)
174  log.debug("cross_correlation_coefficient (based on native_map "
175  "threshold %s) %s", threshold, ccc)
176  return ccc
177 
178 
179 def get_drms_for_backbone(assembly, native_assembly):
180  """
181  Measure the DRMS ob the backbone between two assemblies.
182  @param assembly The DRMS is computed for this assembly
183  @param native_assembly The assembly that acts as a reference
184 
185  Notes:
186  1) The components of the assembly can be proteins or nucleic acids
187  2) If a protein, the c-alphas are used for calculating the drms
188  3) If a nucleic acid, the backbone of C4' atoms is used
189  4) The chains are treated as rigid bodies to speed the calculation.
190 
191  WARNING: if the function fails with a segmentation fault, one of the
192  possible problems is that IMP reads some HETATM as calphas. Check that
193  the chain does not have heteroatoms.
194  """
195  log.debug("Measuring DRMS of the backbone")
196  begin_range = 0
197  ranges = []
198  backbone = []
199  h_chains = IMP.atom.get_by_type(assembly, IMP.atom.CHAIN_TYPE)
200  for h in h_chains:
201  atoms = representation.get_backbone(h)
202  """"
203  for a in atoms:
204  print "atom ===> ",
205  at = IMP.atom.Atom(a)
206  hr = at.get_parent()
207  res = IMP.atom.Residue(hr)
208  ch = IMP.atom.Chain(h)
209  ch.show()
210  print " - ",
211  res.show()
212  print " - ",
213  at.show()
214  print ""
215  """
216  backbone.extend(atoms)
217  end_range = begin_range + len(atoms)
218  ranges.append((begin_range, end_range))
219  begin_range = end_range
220  log.debug("Ranges %s number of atoms %s", ranges, len(backbone))
221  xyzs = [IMP.core.XYZ(leaf) for leaf in backbone]
222  native_chains = IMP.atom.get_by_type(native_assembly, IMP.atom.CHAIN_TYPE)
223  native_backbone = []
224  for h in native_chains:
225  native_backbone.extend(representation.get_backbone(h))
226  native_xyzs = [IMP.core.XYZ(leaf) for leaf in native_backbone]
227  if len(xyzs) != len(native_xyzs):
228  raise ValueError(
229  "Cannot compute DRMS for sets of atoms of different size")
230  log.debug("Getting rigid bodies rmsd")
231  drms = IMP.atom.get_rigid_bodies_drms(xyzs, native_xyzs, ranges)
232  if drms < 0 or math.isnan(drms): # or drms > 100:
233  log.debug(
234  "len(xyzs) = %s. len(native_xyzs) = %s",
235  len(xyzs),
236  len(native_xyzs))
237  log.debug("drms = %s", drms)
238  IMP.atom.write_pdb(assembly, "drms_model_calphas.pdb")
239  IMP.atom.write_pdb(native_assembly, "drms_native_calphas.pdb")
240  raise ValueError("There is a problem with the drms. I wrote the pdbs "
241  "for you: drms_model_calphas.pdb "
242  "drms_native_calphas.pdb")
243  return drms
double get_rigid_bodies_drms(const Vector3DsOrXYZs0 &m1, const Vector3DsOrXYZs1 &m2, const IMP::IntRanges &ranges)
DRMS between two sets of rigid bodies.
def get_assembly_placement_score
Computes the placement score of an assembly respect to the native_assembly.
Definition: comparisons.py:24
def get_placement_scores_from_coordinates
Computes the placement score for each of the components.
Definition: comparisons.py:78
double get_coarse_cc_coefficient(const DensityMap *grid1, const DensityMap *grid2, double grid2_voxel_data_threshold, bool allow_padding=false, FloatPair norm_factors=FloatPair(0., 0.))
Calculates the cross correlation coefficient between two maps.
Utility functions to handle representation.
BoundingBoxD< D > get_union(BoundingBoxD< D > a, const BoundingBoxD< D > &b)
Return the union bounding box.
Definition: BoundingBoxD.h:261
void write_pdb(const Selection &mhd, TextOutput out, unsigned int model=1)
def get_drms_for_backbone
Measure the DRMS ob the backbone between two assemblies.
Definition: comparisons.py:179
Class for sampling a density map from particles.
double get_rmsd(const Selection &s0, const Selection &s1)
DensityHeader create_density_header(const algebra::BoundingBoxD< 3 > &bb, float spacing)
Create a header from a bounding box in 3D.
algebra::BoundingBoxD< 3 > get_bounding_box(const XYZRs &ps)
Get the bounding box.
double get_rmsd(const Vector3DsOrXYZs0 &m1, const Vector3DsOrXYZs1 &m2)
def get_ccc
Threshold - threshold used for the map of the native assembly.
Definition: comparisons.py:129
std::pair< Vector3D, double > get_axis_and_angle(const Rotation3D &rot)
Decompose a Rotation3D object into a rotation around an axis.
Basic utilities for handling cryo-electron microscopy 3D density maps.
A decorator for a particle with x,y,z coordinates.
Definition: XYZ.h:30
def get_placement_score_from_coordinates
Computes the position error (placement distance) and the orientation error (placement angle) of the c...
Definition: comparisons.py:98
Vector3D get_centroid(const Vector3Ds &ps)
Return the centroid of a set of vectors.
Definition: Vector3D.h:68
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...
VectorD< 3 > Vector3D
Definition: VectorD.h:425
def get_components_placement_scores
Compute the placement score of each of the children of an assembly.
Definition: comparisons.py:40
Transformation3D get_transformation_aligning_first_to_second(Vector3Ds a, Vector3Ds b)
Functionality for loading, creating, manipulating and scoring atomic structures.
Hierarchies get_leaves(const Selection &h)