3 from __future__
import print_function
5 from IMP
import ArgumentParser
7 __doc__ =
"Compare output models to a reference structure."
11 def get_placement_scores_from_coordinates(model_components_coords,
12 native_components_coords):
14 Computes the placement score for each of the components
15 @param model_components_coords A list with the coordinates for each
17 @param native_components_coords A list with the coordinates for each
18 component in the native assembly
22 for model_coords, native_coords
in zip(
23 model_components_coords, native_components_coords):
24 distance, angle = get_placement_score_from_coordinates(model_coords,
26 distances.append(distance)
28 return distances, angles
31 def get_placement_score_from_coordinates(model_coords, native_coords):
33 Computes the position error (placement distance) and the orientation
34 error (placement angle) of the coordinates in model_coords respect to
35 the coordinates in native_coords.
36 placement distance - translation between the centroids of the
38 placement angle - Angle in the axis-angle formulation of the rotation
39 aligning the two rigid bodies.
43 translation_vector = native_centroid - model_centroid
44 distance = translation_vector.get_magnitude()
45 if(len(model_coords) != len(native_coords)):
47 "Mismatch in the number of members %d %d " % (
54 return distance, angle
57 def get_rmsd(hierarchy1, hierarchy2):
63 def get_components_placement_scores(assembly, native_assembly, align=False):
65 Compute the placement score of each of the children of an assembly.
66 @param assembly An atom.Molecule object
67 @param native_assembly An atom.Molecule object with the native
68 conformation. Obviously the atoms in assembly and
69 native_assembly must be the same.
70 @param align if True, the coordinates are aligned before the score is
72 @return The function returns 2 lists. The first list contains the
73 placement distances of the children. The second list contains
76 model_coords_per_child = [get_coordinates(c)
77 for c
in assembly.get_children()]
78 native_coords_per_child = [get_coordinates(c)
79 for c
in native_assembly.get_children()]
82 nil = [model_coords.extend(x)
for x
in model_coords_per_child]
84 nil = [native_coords.extend(x)
for x
in native_coords_per_child]
85 T = alg.get_transformation_aligning_first_to_second(model_coords,
88 new_model_coords_per_child = []
89 for c
in model_coords_per_child:
90 coords = [T.get_transformed(x)
for x
in c]
91 new_model_coords_per_child.append(coords)
92 model_coords_per_child = new_model_coords_per_child
93 distances, angles = get_placement_scores_from_coordinates(
94 native_coords_per_child, model_coords_per_child)
95 return distances, angles
100 Compare output models to a reference structure.
101 The reference structure for each subunit is read from the rightmost column
102 of the asmb.input file.
104 p = ArgumentParser(description=desc)
105 p.add_argument(
"-m",
"--max", type=int, dest=
"max", default=
None,
106 help=
"maximum number of models to compare")
107 p.add_argument(
"assembly_file", help=
"assembly file name")
108 p.add_argument(
"proteomics_file", help=
"proteomics file name")
109 p.add_argument(
"mapping_file", help=
"mapping file name")
110 p.add_argument(
"combinations_file", help=
"combinations file name")
111 return p.parse_args()
114 def run(asmb_fn, proteomics_fn, mapping_fn, combs_fn, max_comb):
119 sd.set_was_used(
True)
124 ensmb.set_was_used(
True)
125 mhs = ensmb.get_molecules()
127 for j, mh
in enumerate(mhs):
130 sd.get_component_header(
134 print(
"number of combinations:", len(combs), max_comb)
136 for i, comb
in enumerate(combs[:max_comb]):
139 ensmb.load_combination(comb)
141 for j, mh
in enumerate(mhs):
145 coords1.append(xyz.get_coordinates())
148 coords2.append(xyz.get_coordinates())
150 get_placement_score_from_coordinates(coords1, coords2))
155 print(i, rmsd, scores)
156 results.append((rmsd, scores))
157 ensmb.unload_combination(comb)
163 return run(args.assembly_file, args.proteomics_file, args.mapping_file,
164 args.combinations_file, args.max)
166 if __name__ ==
"__main__":
SettingsData * read_settings(const char *filename)
GenericHierarchies get_leaves(Hierarchy mhd)
Get all the leaves of the bit of hierarchy.
void read_pdb(TextInput input, int model, Hierarchy h)
ProteinsAnchorsSamplingSpace read_protein_anchors_mapping(multifit::ProteomicsData *prots, const std::string &anchors_prot_map_fn, int max_paths=INT_MAX)
Class for storing model, its restraints, constraints, and particles.
double get_rmsd(const Selection &s0, const Selection &s1)
Ensemble * load_ensemble(multifit::SettingsData *sd, Model *mdl, const ProteinsAnchorsSamplingSpace &mapping_data)
double get_rmsd(const Vector3DsOrXYZs0 &m1, const Vector3DsOrXYZs1 &m2)
Fitting atomic structures into a cryo-electron microscopy density map.
std::pair< Vector3D, double > get_axis_and_angle(const Rotation3D &rot)
Decompose a Rotation3D object into a rotation around an axis.
A decorator for a particle with x,y,z coordinates.
ProteomicsData * read_proteomics_data(const char *proteomics_fn)
Proteomics reader.
Vector3D get_centroid(const Vector3Ds &ps)
Return the centroid of a set of vectors.
IntsList read_paths(const char *txt_filename, int max_paths=INT_MAX)
Read paths.
Transformation3D get_transformation_aligning_first_to_second(Vector3Ds a, Vector3Ds b)
Hierarchies get_leaves(const Selection &h)