1 """@namespace IMP.EMageFit.imp_general.comparisons
2 Utility functions for comparisons.
20 log = logging.getLogger(
"comparisons")
23 def get_coordinates(hierarchy):
24 xyz = [
core.XYZ(l)
for l
in atom.get_leaves(hierarchy)]
25 coords = [x.get_coordinates()
for x
in xyz]
31 Computes the placement score of an assembly respect to the native_assembly.
32 @param assembly An atom.Molecule object
33 @param native_assembly An atom.Molecule object
34 @param align If True, assembly is aligned to native_assembly before
35 calculating the placement score
39 native_assembly, align)
40 n = 1. * len(distances)
41 return sum(distances) / n, sum(angles) / n
46 Compute the placement score of each of the children of an assembly.
47 The function does not do any time of alignment of the coordinates
48 @param assembly An atom.Molecule object
49 @param native_assembly An atom.Molecule object with the native conformation
50 Obviously the atoms in assembly and native assembly
52 @param align If True, assembly is aligned to native_assembly before
53 calculating the placement score
54 @return The function returns 2 lists. The first list contains the
55 placement distances of the children. The second list contains the
58 model_coords_per_child = [get_coordinates(c)
59 for c
in assembly.get_children()]
60 native_coords_per_child = [get_coordinates(c)
61 for c
in native_assembly.get_children()]
64 nil = [model_coords.extend(x)
for x
in model_coords_per_child]
66 nil = [native_coords.extend(x)
for x
in native_coords_per_child]
67 T = alg.get_transformation_aligning_first_to_second(model_coords,
70 new_model_coords_per_child = []
71 for c
in model_coords_per_child:
72 coords = [T.get_transformed(x)
for x
in c]
73 new_model_coords_per_child.append(coords)
74 model_coords_per_child = new_model_coords_per_child
76 native_coords_per_child, model_coords_per_child)
77 return distances, angles
81 native_components_coords):
83 Computes the placement score for each of the components
84 @param model_components_coords A list with the coordinates for each
86 @param native_components_coords A list with the coordinates for each
87 component in the native assembly
91 for model_coords, native_coords
in zip(
92 model_components_coords, native_components_coords):
95 distances.append(distance)
97 return distances, angles
102 Computes the position error (placement distance) and the orientation
103 error (placement angle) of the coordinates in model_coords respect to
104 the coordinates in native_coords.
105 placement distance - translation between the centroids of the
107 placement angle - Angle in the axis-angle formulation of the rotation
108 aligning the two rigid bodies.
110 native_centroid = alg.get_centroid(native_coords)
111 model_centroid = alg.get_centroid(model_coords)
112 translation_vector = native_centroid - model_centroid
113 distance = translation_vector.get_magnitude()
114 if(len(model_coords) != len(native_coords)):
116 "Mismatch in the number of members %d %d " % (
119 TT = alg.get_transformation_aligning_first_to_second(model_coords,
121 P = alg.get_axis_and_angle(TT.get_rotation())
123 return distance, angle
126 def get_rmsd(hierarchy1, hierarchy2):
127 xyz1 = [
core.XYZ(l)
for l
in atom.get_leaves(hierarchy1)]
128 xyz2 = [
core.XYZ(l)
for l
in atom.get_leaves(hierarchy2)]
129 return atom.get_rmsd(xyz1, xyz2)
132 def get_ccc(native_assembly, assembly, resolution, voxel_size,
133 threshold, write_maps=
False):
135 Threshold - threshold used for the map of the native assembly. Pixels
136 with values above this threshold in the native map are used for the
137 calculation of the cross_correlation_coefficient
140 particles_native = atom.get_leaves(native_assembly)
141 particles_solution = atom.get_leaves(assembly)
142 bb_native = core.get_bounding_box(
core.XYZs(particles_native))
143 bb_solution = core.get_bounding_box(
core.XYZs(particles_solution))
146 bb_union = alg.get_union(bb_native, bb_solution)
148 border = 4 * voxel_size
149 bottom = bb_union.get_corner(0)
150 bottom += alg.Vector3D(-border, -border, -border)
151 top = bb_union.get_corner(1)
152 top += alg.Vector3D(border, border, border)
153 bb_union = alg.BoundingBox3D(bottom, top)
156 header = em.create_density_header(bb_union, voxel_size)
157 header.set_resolution(resolution)
160 map_native.set_particles(particles_native)
161 map_native.resample()
164 map_solution.set_particles(particles_solution)
165 map_solution.resample()
168 em.write_map(map_solution,
"map_solution.mrc", mrw)
169 em.write_map(map_native,
"map_native.mrc", mrw)
171 map_solution.calcRMS()
177 ccc = coarse_cc.cross_correlation_coefficient(map_solution,
178 map_native, threshold)
179 log.debug(
"cross_correlation_coefficient (based on native_map "
180 "treshold %s) %s", threshold, ccc)
186 Measure the DRMS ob the backbone between two assemblies.
187 @param assembly The DRMS is computed for this assembly
188 @param native_assembly The assembly that acts as a reference
191 1) The components of the assembly can be proteins or nucleic acids
192 2) If a protein, the c-alphas are used for calculating the drms
193 3) If a nucleic acid, the backbone of C4' atoms is used
194 4) The chains are treated as rigid bodies to speed the calculation.
196 WARNING: if the function fails with a segmentation fault, one of the
197 possible problems is that IMP reads some HETATM as calphas. Check that
198 the chain does not have heteroatoms.
200 log.debug(
"Measuring DRMS of the backbone")
204 h_chains = atom.get_by_type(assembly, atom.CHAIN_TYPE)
206 atoms = representation.get_backbone(h)
212 res = atom.Residue(hr)
221 backbone.extend(atoms)
222 end_range = begin_range + len(atoms)
223 ranges.append((begin_range, end_range))
224 begin_range = end_range
225 log.debug(
"Ranges %s number of atoms %s", ranges, len(backbone))
226 xyzs = [
core.XYZ(l)
for l
in backbone]
227 native_chains = atom.get_by_type(native_assembly, atom.CHAIN_TYPE)
228 names = [
atom.Chain(ch).get_id()
for ch
in native_chains]
230 for h
in native_chains:
231 native_backbone.extend(representation.get_backbone(h))
232 native_xyzs = [
core.XYZ(l)
for l
in native_backbone]
233 if len(xyzs) != len(native_xyzs):
235 "Cannot compute DRMS for sets of atoms of different size")
236 log.debug(
"Getting rigid bodies rmsd")
237 drms = atom.get_rigid_bodies_drms(xyzs, native_xyzs, ranges)
238 if drms < 0
or math.isnan(drms):
240 "len(xyzs) = %s. len(native_xyzs) = %s",
243 log.debug(
"drms = %s", drms)
244 atom.write_pdb(assembly,
"drms_model_calphas.pdb")
245 atom.write_pdb(native_assembly,
"drms_native_calphas.pdb")
246 raise ValueError(
"There is a problem with the drms. I wrote the pdbs "
247 "for you: drms_model_calphas.pdb drms_native_calphas.pdb")
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.
Utility functions to handle representation.
Utility functions to handle alignments.
def get_drms_for_backbone
Measure the DRMS ob the backbone between two assemblies.
Responsible for performing coarse fitting between two density objects.
Class for sampling a density map from particles.
double get_rmsd(const Vector3DsOrXYZs0 &m1, const Vector3DsOrXYZs1 &m2)
def get_ccc
Threshold - threshold used for the map of the native assembly.
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...
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...
Store info for a chain of a protein.
def get_components_placement_scores
Compute the placement score of each of the children of an assembly.
Functionality for loading, creating, manipulating and scoring atomic structures.