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