IMP  2.3.0
The Integrative Modeling Platform
fit_fft.py
1 #!/usr/bin/env python
2 
3 __doc__ = "Fit subunits into a density map with FFT."
4 
5 import math
6 import IMP.multifit
7 import IMP.atom
8 import IMP.em
9 from IMP import OptionParser
10 import os
11 import sys
12 
13 multiproc_exception = None
14 try:
15  from multiprocessing import Pool
16  # Detect whether we are running Windows Python via Wine. Wine does not
17  # currently support some named pipe functions which the multiprocessing
18  # module needs: http://bugs.winehq.org/show_bug.cgi?id=17273
19  if sys.platform == 'win32' and 'WINELOADERNOEXEC' in os.environ:
20  multiproc_exception = "Wine does not currently support multiprocessing"
21 except ImportError as detail:
22  multiproc_exception = str(detail)
23 
24 
25 class Fitter(object):
26 
27  def __init__(
28  self,
29  em_map,
30  spacing,
31  resolution,
32  origin,
33  density_threshold,
34  pdb,
35  fits_fn,
36  angle,
37  num_fits,
38  angles_per_voxel,
39  ref_pdb=''):
40  self.em_map = em_map
41  self.spacing = spacing
42  self.resolution = resolution
43  self.threshold = density_threshold
44  self.originx = origin[0]
45  self.originy = origin[1]
46  self.originz = origin[2]
47  self.pdb = pdb
48  self.fits_fn = fits_fn
49  self.angle = angle
50  self.num_fits = num_fits
51  self.angles_per_voxel = angles_per_voxel
52  self.ref_pdb = ref_pdb
53 
54  def run(self):
55  print "resolution is:", self.resolution
56  dmap = IMP.em.read_map(self.em_map)
57  dmap.get_header().set_resolution(self.resolution)
58  dmap.update_voxel_size(self.spacing)
59  dmap.set_origin(IMP.algebra.Vector3D(self.originx,
60  self.originy,
61  self.originz))
62  dmap.set_was_used(True)
63  dmap.get_header().show()
64  mdl = IMP.kernel.Model()
65  mol2fit = IMP.atom.read_pdb(self.pdb, mdl)
66  mh_xyz = IMP.core.XYZs(IMP.core.get_leaves(mol2fit))
67  rb = IMP.atom.create_rigid_body(mol2fit)
69  ff.set_was_used(True)
70  fits = ff.do_global_fitting(dmap, self.threshold, mol2fit,
71  self.angle / 180.0 * math.pi,
72  self.num_fits, self.spacing, 0.5,
73  True, self.angles_per_voxel)
74  fits.set_was_used(True)
75  final_fits = fits.best_fits_
76  if self.ref_pdb != '':
77  ref_mh = IMP.atom.read_pdb(self.ref_pdb, mdl)
78  ref_mh_xyz = IMP.core.XYZs(IMP.core.get_leaves(ref_mh))
79  cur_low = [1e4, 0]
80  for i, fit in enumerate(final_fits):
81  fit.set_index(i)
82  if self.ref_pdb != '':
83  trans = fit.get_fit_transformation()
84  IMP.atom.transform(mol2fit, trans)
85  rmsd = IMP.atom.get_rmsd(mh_xyz, ref_mh_xyz)
86  if rmsd < cur_low[0]:
87  cur_low[0] = rmsd
88  cur_low[1] = i
89  fit.set_rmsd_to_reference(rmsd)
90  IMP.atom.transform(mol2fit, trans.get_inverse())
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>
102 
103 Fit subunits into a density map with FFT."""
104  parser = OptionParser(usage)
105  parser.add_option("-c", "--cpu", dest="cpus", type="int", default=1,
106  help="number of cpus to use (default 1)")
107  parser.add_option("-a", "--angle", dest="angle", type="float",
108  default=30,
109  help="angle delta (degrees) for FFT rotational "
110  "search (default 30)")
111 
112  parser.add_option("-n", "--num", dest="num", type="int",
113  default=100,
114  help="Number of fits to report"
115  "(default 100)")
116 
117  parser.add_option("-v", "--angle_voxel", dest="angle_voxel", type="int",
118  default=10,
119  help="Number of angles to keep per voxel"
120  "(default 10)")
121 
122  # parser.add_option("-n", "--num", dest="num", type="int",
123  # default=100,
124  # help="Number of fits to report"
125  # "(default 100)")
126 
127  options, args = parser.parse_args()
128  if len(args) != 1:
129  parser.error("incorrect number of arguments")
130  return options, args
131 
132 
133 def run(asmb_fn, options):
134  if multiproc_exception is None and options.cpus > 1:
135  work_units = []
136  asmb_input = IMP.multifit.read_settings(asmb_fn)
137  asmb_input.set_was_used(True)
138  em_map = asmb_input.get_assembly_header().get_dens_fn()
139  resolution = asmb_input.get_assembly_header().get_resolution()
140  spacing = asmb_input.get_assembly_header().get_spacing()
141  origin = asmb_input.get_assembly_header().get_origin()
142  for i in range(asmb_input.get_number_of_component_headers()):
143  fits_fn = asmb_input.get_component_header(i).get_transformations_fn()
144  pdb_fn = asmb_input.get_component_header(i).get_filename()
145  f = Fitter(
146  em_map,
147  spacing,
148  resolution,
149  origin,
150  asmb_input.get_assembly_header().get_threshold(),
151  pdb_fn,
152  fits_fn,
153  options.angle,
154  options.num,
155  options.angle_voxel)
156  if multiproc_exception is None and options.cpus > 1:
157  work_units.append(f)
158  else:
159  if options.cpus > 1:
160  options.cpus = 1
161  print >> sys.stderr, """
162 The Python 'multiprocessing' module (available in Python 2.6 and later) is
163 needed to run on multiple CPUs, and could not be found
164 (Python error: '%s').
165 Running on a single processor.""" % multiproc_exception
166  f.run()
167  if multiproc_exception is None and options.cpus > 1:
168  # No point in spawning more processes than components
169  nproc = min(options.cpus, asmb_input.get_number_of_component_headers())
170  p = Pool(processes=nproc)
171  out = list(p.imap_unordered(do_work, work_units))
172 
173 
174 def main():
175  options, args = parse_args()
176  asmb_input = args[0]
177  run(asmb_input, options)
178 
179 if __name__ == "__main__":
180  main()
double get_resolution(kernel::Model *m, kernel::ParticleIndex pi)
Fit a molecule inside its density by local or global FFT.
double get_rmsd(const core::XYZs &s0, const core::XYZs &s1, const IMP::algebra::Transformation3D &tr_for_second)
SettingsData * read_settings(const char *filename)
GenericHierarchies get_leaves(Hierarchy mhd)
Get all the leaves of the bit of hierarchy.
void transform(Hierarchy h, const algebra::Transformation3D &tr)
Transform a hierarchy. This is aware of rigid bodies.
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.
DensityMap * read_map(std::string filename)
Read a density map from a file and return it.
IMP::kernel::OptionParser OptionParser
void show(Hierarchy h, std::ostream &out=std::cout)
Print out a molecular hierarchy.
VectorD< 3 > Vector3D
Definition: VectorD.h:395
IMP::core::RigidBody create_rigid_body(Hierarchy h)
Functionality for loading, creating, manipulating and scoring atomic structures.
void read_pdb(base::TextInput input, int model, Hierarchy h)
Class for storing model, its restraints, constraints, and particles.
Definition: kernel/Model.h:73