IMP logo
IMP Reference Guide  2.6.0
The Integrative Modeling Platform
reference.py
1 #!/usr/bin/env python
2 
3 from __future__ import print_function
4 import IMP.multifit
5 from IMP import OptionParser
6 
7 __doc__ = "Compare output models to a reference structure."
8 
9 # analyse the ensemble, first we will do the rmsd stuff
10 
11 def get_placement_scores_from_coordinates(model_components_coords,
12  native_components_coords):
13  """
14  Computes the placement score for each of the components
15  @param model_components_coords A list with the coordinates for each
16  component
17  @param native_components_coords A list with the coordinates for each
18  component in the native assembly
19  """
20  distances = []
21  angles = []
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,
25  native_coords)
26  distances.append(distance)
27  angles.append(angle)
28  return distances, angles
29 
30 
31 def get_placement_score_from_coordinates(model_coords, native_coords):
32  """
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
37  coordinates
38  placement angle - Angle in the axis-angle formulation of the rotation
39  aligning the two rigid bodies.
40  """
41  native_centroid = IMP.algebra.get_centroid(native_coords)
42  model_centroid = IMP.algebra.get_centroid(model_coords)
43  translation_vector = native_centroid - model_centroid
44  distance = translation_vector.get_magnitude()
45  if(len(model_coords) != len(native_coords)):
46  raise ValueError(
47  "Mismatch in the number of members %d %d " % (
48  len(model_coords),
49  len(native_coords)))
51  native_coords)
52  P = IMP.algebra.get_axis_and_angle(TT.get_rotation())
53  angle = P.second
54  return distance, angle
55 
56 
57 def get_rmsd(hierarchy1, hierarchy2):
58  xyz1 = [IMP.core.XYZ(l) for l in IMP.atom.get_leaves(hierarchy1)]
59  xyz2 = [IMP.core.XYZ(l) for l in IMP.atom.get_leaves(hierarchy2)]
60  return IMP.atom.get_rmsd(xyz1, xyz2)
61 
62 
63 def get_components_placement_scores(assembly, native_assembly, align=False):
64  """
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
71  calculated.
72  @return The function returns 2 lists. The first list contains the
73  placement distances of the children. The second list contains
74  the placement angles
75  """
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()]
80  if align:
81  model_coords = []
82  nil = [model_coords.extend(x) for x in model_coords_per_child]
83  native_coords = []
84  nil = [native_coords.extend(x) for x in native_coords_per_child]
85  T = alg.get_transformation_aligning_first_to_second(model_coords,
86  native_coords)
87  # get aligned coordinates
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
96 
97 
98 def parse_args():
99  usage = """%prog [options] <asmb.input> <proteomics.input>
100  <mapping.input> <combinations>
101 
102 Compare output models to a reference structure.
103 The reference structure for each subunit is read from the rightmost column
104 of the asmb.input file.
105 """
106  parser = OptionParser(usage)
107  parser.add_option("-m", "--max", type="int", dest="max", default=None,
108  help="maximum number of models to compare")
109  (options, args) = parser.parse_args()
110  if len(args) != 4:
111  parser.error("incorrect number of arguments")
112  return options, args
113 
114 
115 def run(asmb_fn, proteomics_fn, mapping_fn, combs_fn, max_comb):
116  # get rmsd for subunits
117  mdl = IMP.Model()
118  combs = IMP.multifit.read_paths(combs_fn)
119  sd = IMP.multifit.read_settings(asmb_fn)
120  sd.set_was_used(True)
121  prot_data = IMP.multifit.read_proteomics_data(proteomics_fn)
122  mapping_data = IMP.multifit.read_protein_anchors_mapping(prot_data,
123  mapping_fn)
124  ensmb = IMP.multifit.load_ensemble(sd, mdl, mapping_data)
125  ensmb.set_was_used(True)
126  mhs = ensmb.get_molecules()
127  mhs_ref = []
128  for j, mh in enumerate(mhs):
129  mhs_ref.append(
131  sd.get_component_header(
132  j).get_reference_fn(
133  ),
134  mdl))
135  print("number of combinations:", len(combs), max_comb)
136  results = []
137  for i, comb in enumerate(combs[:max_comb]):
138  if i % 500 == 0:
139  print(i)
140  ensmb.load_combination(comb)
141  scores = []
142  for j, mh in enumerate(mhs):
143  mh_ref = mhs_ref[j]
144  coords1 = []
145  for xyz in IMP.core.XYZs(IMP.core.get_leaves(mh)):
146  coords1.append(xyz.get_coordinates())
147  coords2 = []
148  for xyz in IMP.core.XYZs(IMP.core.get_leaves(mh_ref)):
149  coords2.append(xyz.get_coordinates())
150  scores.append(
151  get_placement_score_from_coordinates(coords1, coords2))
152  # scores=get_placement_score_from_coordinates(IMP.core.XYZs(IMP.atom.get_leaves(mh)),
153  # IMP.core.XYZs(IMP.atom.get_leaves(mh_ref)))
154 
155  rmsd = get_rmsd(mhs, mhs_ref)
156  print(i, rmsd, scores)
157  results.append((rmsd, scores))
158  ensmb.unload_combination(comb)
159  return results
160 
161 
162 def main():
163  options, args = parse_args()
164  return run(args[0], args[1], args[2], args[3], options.max)
165 
166 if __name__ == "__main__":
167  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.
Definition: Model.h:72
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.
Definition: XYZ.h:30
ProteomicsData * read_proteomics_data(const char *proteomics_fn)
Proteomics reader.
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)
Hierarchies get_leaves(const Selection &h)