3 __doc__ =
"Compare output models to a reference structure."
7 from optparse
import OptionParser
10 def get_placement_scores_from_coordinates(model_components_coords,
11 native_components_coords):
13 Computes the placement score for each of the components
14 @param model_components_coords A list with the coordinates for each
16 @param native_components_coords A list with the coordinates for each
17 component in the native assembly
21 for model_coords, native_coords
in zip(
22 model_components_coords,native_components_coords):
23 distance, angle = get_placement_score_from_coordinates(model_coords,
25 distances.append(distance)
27 return distances, angles
30 def get_placement_score_from_coordinates(model_coords, native_coords):
32 Computes the position error (placement distance) and the orientation
33 error (placement angle) of the coordinates in model_coords respect to
34 the coordinates in native_coords.
35 placement distance - translation between the centroids of the
37 placement angle - Angle in the axis-angle formulation of the rotation
38 aligning the two rigid bodies.
42 translation_vector = native_centroid - model_centroid
43 distance = translation_vector.get_magnitude()
44 if(len(model_coords) != len(native_coords) ):
46 "Mismatch in the number of members %d %d " % (
53 return distance, angle
56 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 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)
for c
in assembly.get_children()]
77 native_coords_per_child = [get_coordinates(c)
for c
in native_assembly.get_children()]
80 nil = [model_coords.extend(x)
for x
in model_coords_per_child]
82 nil = [native_coords.extend(x)
for x
in native_coords_per_child]
83 T = alg.get_transformation_aligning_first_to_second(model_coords,
86 new_model_coords_per_child = []
87 for c
in model_coords_per_child:
88 coords = [T.get_transformed(x)
for x
in c]
89 new_model_coords_per_child.append(coords)
90 model_coords_per_child = new_model_coords_per_child
91 distances, angles = get_placement_scores_from_coordinates(
92 native_coords_per_child, model_coords_per_child)
93 return distances, angles
97 usage =
"""%prog [options] <asmb.input> <proteomics.input>
98 <mapping.input> <combinations>
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.
105 parser.add_option(
"-m",
"--max", type=
"int", dest=
"max", default=
None,
106 help=
"maximum number of models to compare")
107 (options, args) = parser.parse_args()
109 parser.error(
"incorrect number of arguments")
112 def run(asmb_fn,proteomics_fn,mapping_fn,combs_fn,max_comb):
117 sd.set_was_used(
True)
122 ensmb.set_was_used(
True)
123 mhs=ensmb.get_molecules()
125 for j,mh
in enumerate(mhs):
126 mhs_ref.append(
IMP.atom.read_pdb(sd.get_component_header(j).get_reference_fn(),mdl))
127 print "number of combinations:",len(combs),max_comb
129 for i,comb
in enumerate(combs[:max_comb]):
132 ensmb.load_combination(comb)
134 for j,mh
in enumerate(mhs):
138 coords1.append(xyz.get_coordinates())
141 coords2.append(xyz.get_coordinates())
142 scores.append(get_placement_score_from_coordinates(coords1,coords2))
148 results.append((rmsd, scores))
149 ensmb.unload_combination(comb)
153 options,args = parse_args()
154 return run(args[0],args[1],args[2],args[3],options.max)
156 if __name__ ==
"__main__":