IMP logo
IMP Reference Guide  2.6.0
The Integrative Modeling Platform
refine_fft.py
1 #!/usr/bin/env python
2 
3 from __future__ import print_function
4 import math
5 import IMP.multifit
6 import IMP.atom
7 import IMP.em
8 from IMP import OptionParser
9 import os
10 import sys
11 
12 __doc__ = "Refine fitting subunits into a density map with FFT."
13 
14 class Fitter(object):
15 
16  def __init__(
17  self,
18  em_map,
19  spacing,
20  resolution,
21  origin,
22  density_threshold,
23  pdb,
24  fits_fn,
25  angle,
26  num_fits,
27  angles_per_voxel,
28  max_trans,
29  max_angle,
30  ref_pdb=''):
31  self.em_map = em_map
32  self.spacing = spacing
33  self.resolution = resolution
34  self.threshold = density_threshold
35  self.originx = origin[0]
36  self.originy = origin[1]
37  self.originz = origin[2]
38  self.pdb = pdb
39  self.fits_fn = fits_fn
40  self.angle = angle
41  self.num_fits = num_fits
42  self.angles_per_voxel = angles_per_voxel
43  self.max_trans = max_trans
44  self.max_angle = max_angle
45  self.ref_pdb = ref_pdb
46 
47  # TODO - update function
48  def run_local_fitting(self, mol2fit, rb, initial_transformation):
49  print("resolution is:", self.resolution)
50  dmap = IMP.em.read_map(self.em_map)
51  dmap.get_header().set_resolution(self.resolution)
52  dmap.update_voxel_size(self.spacing)
53  dmap.set_origin(IMP.algebra.Vector3D(self.originx,
54  self.originy,
55  self.originz))
56  dmap.set_was_used(True)
57  dmap.get_header().show()
58  mh_xyz = IMP.core.XYZs(IMP.core.get_leaves(mol2fit))
60  ff.set_was_used(True)
61  #
62  do_cluster_fits = True
63  max_clustering_translation = 3
64  max_clustering_rotation = 5
65  num_fits_to_report = 100
66  #
67  fits = ff.do_local_fitting(dmap, self.threshold, mol2fit,
68  self.angle / 180.0 * math.pi,
69  self.max_angle / 180.0 *
70  math.pi, self.max_trans, num_fits_to_report,
71  do_cluster_fits, self.angles_per_voxel,
72  max_clustering_translation, max_clustering_rotation)
73  fits.set_was_used(True)
74  final_fits = fits.best_fits_
75  if self.ref_pdb != '':
76  ref_mh = IMP.atom.read_pdb(self.ref_pdb, mdl)
77  ref_mh_xyz = IMP.core.XYZs(IMP.core.get_leaves(ref_mh))
78  cur_low = [1e4, 0]
79  for i, fit in enumerate(final_fits):
80  fit.set_index(i)
81  if self.ref_pdb != '':
82  trans = fit.get_fit_transformation()
83  IMP.atom.transform(mol2fit, trans)
84  rmsd = IMP.atom.get_rmsd(mh_xyz, ref_mh_xyz)
85  if rmsd < cur_low[0]:
86  cur_low[0] = rmsd
87  cur_low[1] = i
88  fit.set_rmsd_to_reference(rmsd)
89  IMP.atom.transform(mol2fit, trans.get_inverse())
90  fit.set_fit_transformation(trans * initial_transformation)
91  if self.ref_pdb != '':
92  print('from all fits, lowest rmsd to ref:', cur_low)
93  IMP.multifit.write_fitting_solutions(self.fits_fn, final_fits)
94 
95 
96 def do_work(f):
97  f.run()
98 
99 
100 def parse_args():
101  usage = """%prog [options] <assembly input> <refined assembly input> <proteomics.input> <mapping.input> <combinations file> <combination index>
102 
103 Fit subunits locally around a combination solution with FFT."""
104  parser = OptionParser(usage)
105  parser.add_option("-a", "--angle", dest="angle", type="float",
106  default=5,
107  help="angle delta (degrees) for FFT rotational "
108  "search (default 5)")
109 
110  parser.add_option("-n", "--num", dest="num", type="int",
111  default=100,
112  help="Number of fits to report"
113  "(default 100)")
114 
115  parser.add_option("-v", "--angle_voxel", dest="angle_voxel", type="int",
116  default=10,
117  help="Number of angles to keep per voxel"
118  "(default 10)")
119 
120  parser.add_option("-t", "--max_trans", dest="max_trans", type="float",
121  default=10.,
122  help="maximum translational search in A"
123  "(default 10)")
124 
125  parser.add_option("-m", "--max_angle", dest="max_angle", type="float",
126  default=30.,
127  help="maximum angular search in degrees"
128  "(default 50)")
129 
130  options, args = parser.parse_args()
131  if len(args) != 6:
132  parser.error("incorrect number of arguments")
133  return options, args
134 
135 
136 def run(
137  asmb_fn,
138  asmb_refined_fn,
139  proteomics_fn,
140  mapping_fn,
141  combs_fn,
142  comb_ind,
143  options):
144  # get rmsd for subunits
145  mdl1 = IMP.Model()
146  mdl2 = IMP.Model()
147  combs = IMP.multifit.read_paths(combs_fn)
148  asmb_input = IMP.multifit.read_settings(asmb_fn)
149  asmb_input.set_was_used(True)
150  asmb_refined_input = IMP.multifit.read_settings(asmb_refined_fn)
151  asmb_refined_input.set_was_used(True)
152  prot_data = IMP.multifit.read_proteomics_data(proteomics_fn)
153  mapping_data = IMP.multifit.read_protein_anchors_mapping(prot_data,
154  mapping_fn)
155  ensmb = IMP.multifit.load_ensemble(asmb_input, mdl1, mapping_data)
156  ensmb.set_was_used(True)
157  mhs = ensmb.get_molecules()
158 
159  ensmb_ref = IMP.multifit.load_ensemble(asmb_input, mdl2, mapping_data)
160  ensmb_ref.set_was_used(True)
161  mhs_ref = ensmb_ref.get_molecules()
162 
163  ensmb.load_combination(combs[comb_ind])
164 
165  em_map = asmb_input.get_assembly_header().get_dens_fn()
166  resolution = asmb_input.get_assembly_header().get_resolution()
167  spacing = asmb_input.get_assembly_header().get_spacing()
168  origin = asmb_input.get_assembly_header().get_origin()
169 
170  rbs_ref = ensmb_ref.get_rigid_bodies()
171  rbs = ensmb.get_rigid_bodies()
172 
173  for i, mh in enumerate(mhs):
174  fits_fn = asmb_refined_input.get_component_header(
175  i).get_transformations_fn()
176 
177  # todo - get the initial transformation
178  rb_ref = rbs_ref[i]
179  rb = rbs[i]
180 
182  rb_ref.get_reference_frame(),
183  rb.get_reference_frame())
184 
185  pdb_fn = asmb_input.get_component_header(i).get_filename()
186 
187  f = Fitter(
188  em_map, spacing, resolution, origin, asmb_input.get_assembly_header(
189  ).get_threshold(
190  ), pdb_fn, fits_fn, options.angle, options.num, options.angle_voxel,
191  options.max_trans, options.max_angle)
192  f.run_local_fitting(mh, rb, initial_transformation)
193 
194 
195 def main():
196  options, args = parse_args()
197  asmb_input = args[0]
198  asmb_refined_input = args[1]
199  proteomics_fn = args[2]
200  mapping_fn = args[3]
201  combinations_fn = args[4]
202  combination_ind = int(args[5])
203  run(asmb_input, asmb_refined_input, proteomics_fn,
204  mapping_fn, combinations_fn, combination_ind, options)
205 
206 if __name__ == "__main__":
207  main()
Fit a molecule inside its density by local or global FFT.
SettingsData * read_settings(const char *filename)
GenericHierarchies get_leaves(Hierarchy mhd)
Get all the leaves of the bit of hierarchy.
Transformation3D get_transformation_from_first_to_second(const ReferenceFrame3D &a, const ReferenceFrame3D &b)
ProteinsAnchorsSamplingSpace read_protein_anchors_mapping(multifit::ProteomicsData *prots, const std::string &anchors_prot_map_fn, int max_paths=INT_MAX)
void read_pdb(TextInput input, int model, Hierarchy h)
Class for storing model, its restraints, constraints, and particles.
Definition: Model.h:72
double get_rmsd(const Selection &s0, const Selection &s1)
void transform(Hierarchy h, const algebra::Transformation3D &tr)
Transform a hierarchy. This is aware of rigid bodies.
def __init__
input a list of particles, the slope and theta of the sigmoid potential theta is the cutoff distance ...
Ensemble * load_ensemble(multifit::SettingsData *sd, Model *mdl, const ProteinsAnchorsSamplingSpace &mapping_data)
Fitting atomic structures into a cryo-electron microscopy density map.
Basic utilities for handling cryo-electron microscopy 3D density maps.
void write_fitting_solutions(const char *fitting_fn, const FittingSolutionRecords &fit_sols, int num_sols=-1)
Write fitting solutions to a file.
ProteomicsData * read_proteomics_data(const char *proteomics_fn)
Proteomics reader.
IntsList read_paths(const char *txt_filename, int max_paths=INT_MAX)
Read paths.
VectorD< 3 > Vector3D
Definition: VectorD.h:395
double get_resolution(Model *m, ParticleIndex pi)
Estimate the resolution of the hierarchy as used by Representation.
Functionality for loading, creating, manipulating and scoring atomic structures.