IMP logo
IMP Reference Guide  2.20.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 ArgumentParser
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  desc = """
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.
103 """
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()
112 
113 
114 def run(asmb_fn, proteomics_fn, mapping_fn, combs_fn, max_comb):
115  # get rmsd for subunits
116  mdl = IMP.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  args = parse_args()
163  return run(args.assembly_file, args.proteomics_file, args.mapping_file,
164  args.combinations_file, args.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:86
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)
Return the centroid of a set of vectors.
Definition: Vector3D.h:68
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)