1 """@namespace IMP.EMageFit.imp_general.comparisons
2 Utility functions for comparisons.
19 log = logging.getLogger(
"comparisons")
22 def get_coordinates(hierarchy):
24 coords = [x.get_coordinates()
for x
in xyz]
30 Computes the placement score of an assembly respect to the native_assembly.
31 @param assembly An IMP.atom.Molecule object
32 @param native_assembly An IMP.atom.Molecule object
33 @param align If True, assembly is aligned to native_assembly before
34 calculating the placement score
38 native_assembly, align)
39 n = 1. * len(distances)
40 return sum(distances) / n, sum(angles) / n
45 Compute the placement score of each of the children of an assembly.
46 The function does not do any time of alignment of the coordinates
47 @param assembly An IMP.atom.Molecule object
48 @param native_assembly An IMP.atom.Molecule object with the native conformation
49 Obviously the atoms in assembly and native assembly
51 @param align If True, assembly is aligned to native_assembly before
52 calculating the placement score
53 @return The function returns 2 lists. The first list contains the
54 placement distances of the children. The second list contains the
57 model_coords_per_child = [get_coordinates(c)
58 for c
in assembly.get_children()]
59 native_coords_per_child = [get_coordinates(c)
60 for c
in native_assembly.get_children()]
63 nil = [model_coords.extend(x)
for x
in model_coords_per_child]
65 nil = [native_coords.extend(x)
for x
in native_coords_per_child]
67 model_coords, native_coords)
69 new_model_coords_per_child = []
70 for c
in model_coords_per_child:
71 coords = [T.get_transformed(x)
for x
in c]
72 new_model_coords_per_child.append(coords)
73 model_coords_per_child = new_model_coords_per_child
75 native_coords_per_child, model_coords_per_child)
76 return distances, angles
80 native_components_coords):
82 Computes the placement score for each of the components
83 @param model_components_coords A list with the coordinates for each
85 @param native_components_coords A list with the coordinates for each
86 component in the native assembly
90 for model_coords, native_coords
in zip(
91 model_components_coords, native_components_coords):
94 distances.append(distance)
96 return distances, angles
101 Computes the position error (placement distance) and the orientation
102 error (placement angle) of the coordinates in model_coords respect to
103 the coordinates in native_coords.
104 placement distance - translation between the centroids of the
106 placement angle - Angle in the axis-angle formulation of the rotation
107 aligning the two rigid bodies.
111 translation_vector = native_centroid - model_centroid
112 distance = translation_vector.get_magnitude()
113 if(len(model_coords) != len(native_coords)):
115 "Mismatch in the number of members %d %d " % (
122 return distance, angle
125 def get_rmsd(hierarchy1, hierarchy2):
131 def get_ccc(native_assembly, assembly, resolution, voxel_size,
132 threshold, write_maps=
False):
134 Threshold - threshold used for the map of the native assembly. Pixels
135 with values above this threshold in the native map are used for the
136 calculation of the cross_correlation_coefficient
147 border = 4 * voxel_size
148 bottom = bb_union.get_corner(0)
150 top = bb_union.get_corner(1)
156 header.set_resolution(resolution)
159 map_native.set_particles(particles_native)
160 map_native.resample()
163 map_solution.set_particles(particles_solution)
164 map_solution.resample()
167 IMP.em.write_map(map_solution,
"map_solution.mrc", mrw)
168 IMP.em.write_map(map_native,
"map_native.mrc", mrw)
170 map_solution.calcRMS()
176 ccc = coarse_cc.cross_correlation_coefficient(map_solution,
177 map_native, threshold)
178 log.debug(
"cross_correlation_coefficient (based on native_map "
179 "treshold %s) %s", threshold, ccc)
185 Measure the DRMS ob the backbone between two assemblies.
186 @param assembly The DRMS is computed for this assembly
187 @param native_assembly The assembly that acts as a reference
190 1) The components of the assembly can be proteins or nucleic acids
191 2) If a protein, the c-alphas are used for calculating the drms
192 3) If a nucleic acid, the backbone of C4' atoms is used
193 4) The chains are treated as rigid bodies to speed the calculation.
195 WARNING: if the function fails with a segmentation fault, one of the
196 possible problems is that IMP reads some HETATM as calphas. Check that
197 the chain does not have heteroatoms.
199 log.debug(
"Measuring DRMS of the backbone")
203 h_chains = IMP.atom.get_by_type(assembly, IMP.atom.CHAIN_TYPE)
205 atoms = representation.get_backbone(h)
209 at = IMP.atom.Atom(a)
211 res = IMP.atom.Residue(hr)
212 ch = IMP.atom.Chain(h)
220 backbone.extend(atoms)
221 end_range = begin_range + len(atoms)
222 ranges.append((begin_range, end_range))
223 begin_range = end_range
224 log.debug(
"Ranges %s number of atoms %s", ranges, len(backbone))
226 native_chains = IMP.atom.get_by_type(native_assembly, IMP.atom.CHAIN_TYPE)
229 for h
in native_chains:
230 native_backbone.extend(representation.get_backbone(h))
231 native_xyzs = [
IMP.core.XYZ(l)
for l
in native_backbone]
232 if len(xyzs) != len(native_xyzs):
234 "Cannot compute DRMS for sets of atoms of different size")
235 log.debug(
"Getting rigid bodies rmsd")
237 if drms < 0
or math.isnan(drms):
239 "len(xyzs) = %s. len(native_xyzs) = %s",
242 log.debug(
"drms = %s", drms)
245 raise ValueError(
"There is a problem with the drms. I wrote the pdbs "
246 "for you: drms_model_calphas.pdb 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.
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.
Responsible for performing coarse fitting between two density objects.
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...
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.
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)