IMP  2.2.1
The Integrative Modeling Platform
reference.py
1 #!/usr/bin/env python
2 
3 __doc__ = "Compare output models to a reference structure."
4 
5 # analyse the ensemble, first we will do the rmsd stuff
6 import IMP.multifit
7 from IMP import OptionParser
8 
9 
10 def get_placement_scores_from_coordinates(model_components_coords,
11  native_components_coords):
12  """
13  Computes the placement score for each of the components
14  @param model_components_coords A list with the coordinates for each
15  component
16  @param native_components_coords A list with the coordinates for each
17  component in the native assembly
18  """
19  distances = []
20  angles = []
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,
24  native_coords)
25  distances.append(distance)
26  angles.append(angle)
27  return distances, angles
28 
29 
30 def get_placement_score_from_coordinates(model_coords, native_coords):
31  """
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
36  coordinates
37  placement angle - Angle in the axis-angle formulation of the rotation
38  aligning the two rigid bodies.
39  """
40  native_centroid = IMP.algebra.get_centroid(native_coords)
41  model_centroid = IMP.algebra.get_centroid(model_coords)
42  translation_vector = native_centroid - model_centroid
43  distance = translation_vector.get_magnitude()
44  if(len(model_coords) != len(native_coords)):
45  raise ValueError(
46  "Mismatch in the number of members %d %d " % (
47  len(model_coords),
48  len(native_coords)))
50  native_coords)
51  P = IMP.algebra.get_axis_and_angle(TT.get_rotation())
52  angle = P.second
53  return distance, angle
54 
55 
56 def get_rmsd(hierarchy1, hierarchy2):
57  xyz1 = [IMP.core.XYZ(l) for l in IMP.atom.get_leaves(hierarchy1)]
58  xyz2 = [IMP.core.XYZ(l) for l in IMP.atom.get_leaves(hierarchy2)]
59  return IMP.atom.get_rmsd(xyz1, xyz2)
60 
61 
62 def get_components_placement_scores(assembly, native_assembly, align=False):
63  """
64  Compute the placement score of each of the children of an assembly.
65  @param assembly An atom.Molecule object
66  @param native_assembly An atom.Molecule object with the native
67  conformation. Obviously the atoms in assembly and
68  native_assembly must be the same.
69  @param align if True, the coordinates are aligned before the score is
70  calculated.
71  @return The function returns 2 lists. The first list contains the
72  placement distances of the children. The second list contains
73  the placement angles
74  """
75  model_coords_per_child = [get_coordinates(c)
76  for c in assembly.get_children()]
77  native_coords_per_child = [get_coordinates(c)
78  for c in native_assembly.get_children()]
79  if align:
80  model_coords = []
81  nil = [model_coords.extend(x) for x in model_coords_per_child]
82  native_coords = []
83  nil = [native_coords.extend(x) for x in native_coords_per_child]
84  T = alg.get_transformation_aligning_first_to_second(model_coords,
85  native_coords)
86  # get aligned coordinates
87  new_model_coords_per_child = []
88  for c in model_coords_per_child:
89  coords = [T.get_transformed(x) for x in c]
90  new_model_coords_per_child.append(coords)
91  model_coords_per_child = new_model_coords_per_child
92  distances, angles = get_placement_scores_from_coordinates(
93  native_coords_per_child, model_coords_per_child)
94  return distances, angles
95 
96 
97 def parse_args():
98  usage = """%prog [options] <asmb.input> <proteomics.input>
99  <mapping.input> <combinations>
100 
101 Compare output models to a reference structure.
102 The reference structure for each subunit is read from the rightmost column
103 of the asmb.input file.
104 """
105  parser = OptionParser(usage)
106  parser.add_option("-m", "--max", type="int", dest="max", default=None,
107  help="maximum number of models to compare")
108  (options, args) = parser.parse_args()
109  if len(args) != 4:
110  parser.error("incorrect number of arguments")
111  return options, args
112 
113 
114 def run(asmb_fn, proteomics_fn, mapping_fn, combs_fn, max_comb):
115  # get rmsd for subunits
116  mdl = IMP.kernel.Model()
117  combs = IMP.multifit.read_paths(combs_fn)
118  sd = IMP.multifit.read_settings(asmb_fn)
119  sd.set_was_used(True)
120  prot_data = IMP.multifit.read_proteomics_data(proteomics_fn)
121  mapping_data = IMP.multifit.read_protein_anchors_mapping(prot_data,
122  mapping_fn)
123  ensmb = IMP.multifit.load_ensemble(sd, mdl, mapping_data)
124  ensmb.set_was_used(True)
125  mhs = ensmb.get_molecules()
126  mhs_ref = []
127  for j, mh in enumerate(mhs):
128  mhs_ref.append(
130  sd.get_component_header(
131  j).get_reference_fn(
132  ),
133  mdl))
134  print "number of combinations:", len(combs), max_comb
135  results = []
136  for i, comb in enumerate(combs[:max_comb]):
137  if i % 500 == 0:
138  print i
139  ensmb.load_combination(comb)
140  scores = []
141  for j, mh in enumerate(mhs):
142  mh_ref = mhs_ref[j]
143  coords1 = []
144  for xyz in IMP.core.XYZs(IMP.core.get_leaves(mh)):
145  coords1.append(xyz.get_coordinates())
146  coords2 = []
147  for xyz in IMP.core.XYZs(IMP.core.get_leaves(mh_ref)):
148  coords2.append(xyz.get_coordinates())
149  scores.append(
150  get_placement_score_from_coordinates(coords1, coords2))
151  # scores=get_placement_score_from_coordinates(IMP.core.XYZs(IMP.atom.get_leaves(mh)),
152  # IMP.core.XYZs(IMP.atom.get_leaves(mh_ref)))
153 
154  rmsd = get_rmsd(mhs, mhs_ref)
155  print i, rmsd, scores
156  results.append((rmsd, scores))
157  ensmb.unload_combination(comb)
158  return results
159 
160 
161 def main():
162  options, args = parse_args()
163  return run(args[0], args[1], args[2], args[3], options.max)
164 
165 if __name__ == "__main__":
166  main()
double get_rmsd(const core::XYZs &s0, const core::XYZs &s1, const IMP::algebra::Transformation3D &tr_for_second)
SettingsData * read_settings(const char *filename)
GenericHierarchies get_leaves(Hierarchy mhd)
Get all the leaves of the bit of hierarchy.
ProteinsAnchorsSamplingSpace read_protein_anchors_mapping(multifit::ProteomicsData *prots, const std::string &anchors_prot_map_fn, int max_paths=INT_MAX)
Ensemble * load_ensemble(multifit::SettingsData *sd, Model *mdl, const ProteinsAnchorsSamplingSpace &mapping_data)
double get_rmsd(const Vector3DsOrXYZs0 &m1, const Vector3DsOrXYZs1 &m2)
See IMP.multifit for more information.
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.
Definition: XYZ.h:30
ProteomicsData * read_proteomics_data(const char *proteomics_fn)
Proteomics reader.
IMP::kernel::OptionParser OptionParser
Vector3D get_centroid(const Vector3Ds &ps)
Returns the centroid of a set of vectors.
Definition: Vector3D.h:56
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)
void read_pdb(base::TextInput input, int model, Hierarchy h)
Hierarchies get_leaves(const Selection &h)
Class for storing model, its restraints, constraints, and particles.
Definition: kernel/Model.h:72