IMP logo
IMP Reference Guide  2.21.0
The Integrative Modeling Platform
score.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__ = "Score each of a set of combinations."
8 
9 # analyse the ensemble, first we will do the rmsd stuff
10 
11 
12 def get_color_map():
13  colors = {}
14  colors["Rpt1"] = [0.78, 0.78, 0.73]
15  colors["Rpt2"] = [0.78, 0.66, 0.58]
16  colors["Rpt3"] = [0.77, 0.43, 0.5]
17  colors["Rpt4"] = [0.76, 0.29, 0.67]
18  colors["Rpt5"] = [0.51, 0.14, 0.75]
19  colors["Rpt6"] = [0.0, 0., 0.75]
20  colors["Rpn1"] = [0.34, 0.36, 0.27]
21  colors["Rpn2"] = [0.42, 0.43, 0.36]
22  colors["Rpn3"] = [0.49, 0.5, 0.44]
23  colors["Rpn5"] = [0.56, 0.57, 0.51]
24  colors["Rpn6"] = [0.64, 0.64, 0.59]
25  colors["Rpn7"] = [0.71, 0.71, 0.66]
26  colors["Rpn8"] = [0.78, 0.78, 0.74]
27  colors["Rpn9"] = [1, 0, 0]
28  colors["Rpn10"] = [0, 1, 0]
29  colors["Rpn11"] = [0, 0, 1]
30  colors["Rpn12"] = [0.5, 0.2, 0.4]
31  colors["a1"] = [0.78, 0.78, 0.73]
32  colors["a2"] = [0.78, 0.66, 0.58]
33  colors["a3"] = [0.77, 0.43, 0.5]
34  colors["a4"] = [0.76, 0.29, 0.67]
35  colors["a5"] = [0.51, 0.14, 0.75]
36  colors["a6"] = [0.0, 0., 0.75]
37  colors["a7"] = [0.34, 0.36, 0.27]
38  colors["a8"] = [0.42, 0.43, 0.36]
39  colors["a9"] = [0.49, 0.5, 0.44]
40  colors["a10"] = [0.56, 0.57, 0.51]
41 
42  colors["a11"] = [0.78, 0.78, 0.73]
43  colors["a12"] = [0.78, 0.66, 0.58]
44  colors["a13"] = [0.77, 0.43, 0.5]
45  colors["a14"] = [0.76, 0.29, 0.67]
46  colors["a15"] = [0.51, 0.14, 0.75]
47  colors["a16"] = [0.0, 0., 0.75]
48  colors["a17"] = [0.34, 0.36, 0.27]
49  colors["a18"] = [0.42, 0.43, 0.36]
50  colors["a19"] = [0.49, 0.5, 0.44]
51  colors["a20"] = [0.56, 0.57, 0.51]
52 
53  colors["a21"] = [0.78, 0.78, 0.73]
54  colors["a22"] = [0.78, 0.66, 0.58]
55  colors["a23"] = [0.77, 0.43, 0.5]
56  colors["a24"] = [0.76, 0.29, 0.67]
57  colors["a25"] = [0.51, 0.14, 0.75]
58  colors["a26"] = [0.0, 0., 0.75]
59  colors["a27"] = [0.34, 0.36, 0.27]
60  colors["a28"] = [0.42, 0.43, 0.36]
61  colors["a29"] = [0.49, 0.5, 0.44]
62  colors["a30"] = [0.56, 0.57, 0.51]
63  return colors
64 
65 
66 def decompose(dmap, mhs):
67  full_sampled_map = IMP.em.SampledDensityMap(dmap.get_header())
68  all_ps = []
69  for mh in mhs:
70  all_ps += IMP.core.get_leaves(mh)
71  full_sampled_map.set_particles(all_ps)
72  full_sampled_map.resample()
73  full_sampled_map.calcRMS()
74  upper = (
75  dmap.get_number_of_voxels(
76  ) * dmap.get_header(
77  ).dmean * full_sampled_map.get_header(
78  ).dmean) / len(
79  mhs)
80  lower = dmap.get_number_of_voxels(
81  ) * dmap.get_header(
82  ).rms * full_sampled_map.get_header(
83  ).rms
84  norm_factors = [upper, lower]
85  print("===============my norm factors:", upper, lower)
86  return norm_factors
87 
88 
89 def score_each_protein(dmap, mhs, sd):
90  norm_factors = decompose(dmap, mhs)
91  scores = []
92  # mdl = mhs[0].get_model()
93  for i in range(len(mhs)):
94  leaves = IMP.core.get_leaves(mhs[i])
95  rb = IMP.core.RigidMember(leaves[0]).get_rigid_body()
96  mh_dmap = IMP.em.SampledDensityMap(dmap.get_header())
97  mh_dmap.set_particles(leaves)
98  mh_dmap.resample()
99  mh_dmap.calcRMS()
101  sd.get_component_header(i).get_transformations_fn())
102  mh_scores = []
103  for fit in fits[:15]:
104  IMP.core.transform(rb, fit.get_fit_transformation())
105  mh_dmap.resample()
106  mh_scores.append(
108  dmap,
109  mh_dmap,
110  0.,
111  False,
112  norm_factors))
113  IMP.core.transform(rb, fit.get_fit_transformation().get_inverse())
114  scores.append(mh_scores)
115  print("=====mol", i, mh_scores)
116  return scores
117 
118 
119 def usage():
120  desc = "Score each of a set of combinations."
121  p = ArgumentParser(description=desc)
122  p.add_argument("-m", "--max", dest="max", type=int, default=999999999,
123  help="maximum number of fits considered")
124  p.add_argument("assembly_file", help="assembly file name")
125  p.add_argument("proteomics_file", help="proteomics file name")
126  p.add_argument("mapping_file", help="mapping file name")
127  p.add_argument("param_file", help="parameter file name")
128  p.add_argument("combinations_file", help="combinations file name")
129  p.add_argument("scores_file", help="output scores file name")
130  return p.parse_args()
131 
132 
133 def run(asmb_fn, proteomics_fn, mapping_fn, params_fn, combs_fn,
134  scored_comb_output_fn, max_comb):
135  asmb = IMP.multifit.read_settings(asmb_fn)
136  dmap = IMP.em.read_map(asmb.get_assembly_header().get_dens_fn())
137  dmap.get_header().set_resolution(
138  asmb.get_assembly_header().get_resolution())
139  dmap.update_voxel_size(asmb.get_assembly_header().get_spacing())
140  dmap.set_origin(asmb.get_assembly_header().get_origin())
141  threshold = asmb.get_assembly_header().get_threshold()
142  combs = IMP.multifit.read_paths(combs_fn)
143  # get rmsd for subunits
144  # colors = get_color_map()
145  # names = list(colors.keys())
146  print(params_fn)
147  alignment_params = IMP.multifit.AlignmentParams(params_fn)
148  alignment_params.show()
149 
150  IMP.set_log_level(IMP.TERSE)
151  print("=========", combs_fn)
152  combs = IMP.multifit.read_paths(combs_fn)
153  print("=========1")
154  # sd=IMP.multifit.read_settings(asmb_fn)
155  print("=========2")
156  prot_data = IMP.multifit.read_proteomics_data(proteomics_fn)
157  print("=========3")
159  prot_data, mapping_fn)
160  print("=========4")
161  _ = mapping_data.get_anchors()
162  print("=========5")
163  ensmb = IMP.multifit.Ensemble(asmb, mapping_data)
164  print("=========6")
165  # load all proteomics restraints
167  mapping_data, asmb, alignment_params)
168  align.set_fast_scoring(False)
169  print("align")
170  mdl = align.get_model()
171  mhs = align.get_molecules()
172  align.add_states_and_filters()
173  # rbs = align.get_rigid_bodies()
174  print(IMP.core.RigidMember(
175  IMP.core.get_leaves(mhs[0])[0]).get_rigid_body())
176  align.set_density_map(dmap, threshold)
177  for i, mh in enumerate(mhs):
178  ensmb.add_component_and_fits(
180  asmb.get_component_header(i).get_transformations_fn()))
181  '''
182  try:
183  rgb = colors[mh.get_name()]
184  except:
185  rgb = colors[names[i]]
186  color = IMP.display.Color(rgb[0], rgb[1], rgb[2])
187  for p in IMP.core.get_leaves(mh):
188  g= IMP.display.XYZRGeometry(p)
189  g.set_color(color)
190  gs.append(g)
191  '''
192  all_leaves = []
193  for mh in mhs:
194  all_leaves += IMP.core.XYZs(IMP.core.get_leaves(mh))
195 
196  align.add_all_restraints()
197  print("====1")
198  rs = align.get_restraint_set().get_restraints()
199  print("Get number of restraints:", len(rs))
200  for r in rs:
201  rr = IMP.RestraintSet.get_from(r)
202  for i in range(rr.get_number_of_restraints()):
203  print(rr.get_restraint(i).get_name())
204  output = open(scored_comb_output_fn, "w")
205  # load ref structure
206  ref_mhs = []
207  all_ref_leaves = []
208  for i in range(asmb.get_number_of_component_headers()):
209  c = asmb.get_component_header(i)
210  fn = c.get_reference_fn()
211  if fn:
212  ref_mhs.append(IMP.atom.read_pdb(fn, mdl))
213  all_ref_leaves += IMP.core.get_leaves(ref_mhs[-1])
214  for r in rs:
215  rr = IMP.RestraintSet.get_from(r)
216  for i in range(rr.get_number_of_restraints()):
217  output.write(rr.get_restraint(i).get_name() + "|")
218  output.write("\n")
219  # add fit restraint
220  fitr = IMP.em.FitRestraint(all_leaves, dmap)
221  sf = IMP.core.RestraintsScoringFunction(rs + [fitr])
222  print("Number of combinations:", len(combs[:max_comb]))
223 
224  print("native score")
225  num_violated = 0
226  for r in rs:
227  rr = IMP.RestraintSet.get_from(r)
228  for j in range(rr.get_number_of_restraints()):
229  print(rr.get_restraint(j).get_name(), rr.evaluate(False))
230 
231  prev_name = ''
232  for i, comb in enumerate(combs[:max_comb]):
233  print("Scoring combination:", comb)
234  ensmb.load_combination(comb)
235  num_violated = 0
236  for r in rs:
237  rr = IMP.RestraintSet.get_from(r)
238  for j in range(rr.get_number_of_restraints()):
239  print(rr.get_restraint(j).get_name())
240  rscore = rr.evaluate(False)
241  if rscore > 5:
242  num_violated = num_violated + 1
243  IMP.atom.write_pdb(mhs, "model.%d.pdb" % (i))
244  print(str(all_leaves[0]) + " :: " + str(all_leaves[-1]))
245  score = sf.evaluate(False)
246  num_violated = 0
247  msg = "COMB" + str(i) + "|"
248  for r in rs:
249  rr = IMP.RestraintSet.get_from(r)
250  for j in range(rr.get_number_of_restraints()):
251  current_name = rr.get_restraint(j).get_name()
252  if current_name != prev_name:
253  msg += ' ' + current_name + ' '
254  prev_name = current_name
255  rscore = rr.get_restraint(j).evaluate(False)
256  msg += str(rscore) + "|"
257  if rscore > 5:
258  num_violated = num_violated + 1
259  # msg+="|"+str(score)+"|"+str(num_violated)+"|\n"
260  msg += "|" + str(
261  score) + "|" + str(
262  num_violated) + "||||" + str(
263  fitr.evaluate(False)) + "||:"
264  if all_ref_leaves:
265  msg += str(IMP.atom.get_rmsd(IMP.core.XYZs(all_leaves),
266  IMP.core.XYZs(all_ref_leaves)))
267  output.write(msg + "\n")
268  print(msg)
269  ensmb.unload_combination(comb)
270  output.close()
271 
272 
273 def main():
274  args = usage()
275  run(args.assembly_file, args.proteomics_file, args.mapping_file,
276  args.param_file, args.combinations_file, args.scores_file, args.max)
277 
278 
279 if __name__ == "__main__":
280  main()
An ensemble of fitting solutions.
double get_coarse_cc_coefficient(const DensityMap *grid1, const DensityMap *grid2, double grid2_voxel_data_threshold, bool allow_padding=false, FloatPair norm_factors=FloatPair(0., 0.))
Calculates the cross correlation coefficient between two maps.
void write_pdb(const Selection &mhd, TextOutput out, unsigned int model=1)
def evaluate
Evaluate the score of the restraint.
Create a scoring function on a list of restraints.
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)
Align proteomics graph to EM density map.
Class for sampling a density map from particles.
double get_rmsd(const Selection &s0, const Selection &s1)
void transform(XYZ a, const algebra::Transformation3D &tr)
Apply a transformation to the particle.
Fitting atomic structures into a cryo-electron microscopy density map.
ProteomicsData * read_proteomics_data(const char *proteomics_fn)
Proteomics reader.
void set_log_level(LogLevel l)
Set the current global log level.
IntsList read_paths(const char *txt_filename, int max_paths=INT_MAX)
Read paths.
Calculate score based on fit to EM map.
Definition: FitRestraint.h:39
FittingSolutionRecords read_fitting_solutions(const char *fitting_fn)
Fitting solutions reader.
double get_resolution(Model *m, ParticleIndex pi)
Estimate the resolution of the hierarchy as used by Representation.