1 """@namespace IMP.EMageFit.imp_general.comparisons
2 Utility functions for comparisons.
15 log = logging.getLogger(
"comparisons")
18 def get_coordinates(hierarchy):
20 coords = [x.get_coordinates()
for x
in xyz]
26 Computes the placement score of an assembly respect to the
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
35 native_assembly, align)
36 n = 1. * len(distances)
37 return sum(distances) / n, sum(angles) / n
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
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()]
60 for x
in model_coords_per_child:
61 model_coords.extend(x)
63 for x
in native_coords_per_child:
64 native_coords.extend(x)
66 model_coords, native_coords)
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
74 native_coords_per_child, model_coords_per_child)
75 return distances, angles
79 native_components_coords):
81 Computes the placement score for each of the components
82 @param model_components_coords A list with the coordinates for each
84 @param native_components_coords A list with the coordinates for each
85 component in the native assembly
89 for model_coords, native_coords
in zip(
90 model_components_coords, native_components_coords):
93 distances.append(distance)
95 return distances, angles
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
105 placement angle - Angle in the axis-angle formulation of the rotation
106 aligning the two rigid bodies.
110 translation_vector = native_centroid - model_centroid
111 distance = translation_vector.get_magnitude()
112 if len(model_coords) != len(native_coords):
114 "Mismatch in the number of members %d %d "
115 % (len(model_coords), len(native_coords)))
120 return distance, angle
123 def get_rmsd(hierarchy1, hierarchy2):
129 def get_ccc(native_assembly, assembly, resolution, voxel_size,
130 threshold, write_maps=
False):
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
144 border = 4 * voxel_size
145 bottom = bb_union.get_corner(0)
147 top = bb_union.get_corner(1)
153 header.set_resolution(resolution)
156 map_native.set_particles(particles_native)
157 map_native.resample()
160 map_solution.set_particles(particles_solution)
161 map_solution.resample()
164 IMP.em.write_map(map_solution,
"map_solution.mrc", mrw)
165 IMP.em.write_map(map_native,
"map_native.mrc", mrw)
167 map_solution.calcRMS()
173 map_native, threshold)
174 log.debug(
"cross_correlation_coefficient (based on native_map "
175 "threshold %s) %s", threshold, ccc)
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
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.
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.
195 log.debug(
"Measuring DRMS of the backbone")
199 h_chains = IMP.atom.get_by_type(assembly, IMP.atom.CHAIN_TYPE)
201 atoms = representation.get_backbone(h)
205 at = IMP.atom.Atom(a)
207 res = IMP.atom.Residue(hr)
208 ch = IMP.atom.Chain(h)
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))
222 native_chains = IMP.atom.get_by_type(native_assembly, IMP.atom.CHAIN_TYPE)
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):
229 "Cannot compute DRMS for sets of atoms of different size")
230 log.debug(
"Getting rigid bodies rmsd")
232 if drms < 0
or math.isnan(drms):
234 "len(xyzs) = %s. len(native_xyzs) = %s",
237 log.debug(
"drms = %s", drms)
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")
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.
def get_placement_scores_from_coordinates
Computes the placement score for each of the components.
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.
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.
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.
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.
def get_placement_score_from_coordinates
Computes the position error (placement distance) and the orientation error (placement angle) of the c...
Vector3D get_centroid(const Vector3Ds &ps)
Return the centroid of a set of vectors.
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...
def get_components_placement_scores
Compute the placement score of each of the children of an assembly.
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)