IMP logo
IMP Reference Guide  2.20.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 ArgumentParser
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  desc = """
102 Fit subunits locally around a combination solution with FFT."""
103  p = ArgumentParser(description=desc)
104  p.add_argument("-a", "--angle", dest="angle", type=float, default=5,
105  help="angle delta (degrees) for FFT rotational "
106  "search (default 5)")
107 
108  p.add_argument("-n", "--num", dest="num", type=int, default=100,
109  help="Number of fits to report (default 100)")
110 
111  p.add_argument("-v", "--angle_voxel", dest="angle_voxel", type=int,
112  default=10,
113  help="Number of angles to keep per voxel (default 10)")
114 
115  p.add_argument("-t", "--max_trans", dest="max_trans", type=float,
116  default=10.,
117  help="maximum translational search in A (default 10)")
118 
119  p.add_argument("-m", "--max_angle", dest="max_angle", type=float,
120  default=30.,
121  help="maximum angular search in degrees (default 50)")
122  p.add_argument("assembly_file", help="assembly file name")
123  p.add_argument("ref_assembly_file", help="refined assembly file name")
124  p.add_argument("proteomics_file", help="proteomics file name")
125  p.add_argument("mapping_file", help="mapping file name")
126  p.add_argument("combinations_file", help="combinations file name")
127  p.add_argument("combination_index", type=int,
128  help="number of the combination to read from the "
129  "combinations file")
130  return p.parse_args()
131 
132 
133 def run(
134  asmb_fn,
135  asmb_refined_fn,
136  proteomics_fn,
137  mapping_fn,
138  combs_fn,
139  comb_ind,
140  options):
141  # get rmsd for subunits
142  mdl1 = IMP.Model()
143  mdl2 = IMP.Model()
144  combs = IMP.multifit.read_paths(combs_fn)
145  asmb_input = IMP.multifit.read_settings(asmb_fn)
146  asmb_input.set_was_used(True)
147  asmb_refined_input = IMP.multifit.read_settings(asmb_refined_fn)
148  asmb_refined_input.set_was_used(True)
149  prot_data = IMP.multifit.read_proteomics_data(proteomics_fn)
150  mapping_data = IMP.multifit.read_protein_anchors_mapping(prot_data,
151  mapping_fn)
152  ensmb = IMP.multifit.load_ensemble(asmb_input, mdl1, mapping_data)
153  ensmb.set_was_used(True)
154  mhs = ensmb.get_molecules()
155 
156  ensmb_ref = IMP.multifit.load_ensemble(asmb_input, mdl2, mapping_data)
157  ensmb_ref.set_was_used(True)
158  mhs_ref = ensmb_ref.get_molecules()
159 
160  ensmb.load_combination(combs[comb_ind])
161 
162  em_map = asmb_input.get_assembly_header().get_dens_fn()
163  resolution = asmb_input.get_assembly_header().get_resolution()
164  spacing = asmb_input.get_assembly_header().get_spacing()
165  origin = asmb_input.get_assembly_header().get_origin()
166 
167  rbs_ref = ensmb_ref.get_rigid_bodies()
168  rbs = ensmb.get_rigid_bodies()
169 
170  for i, mh in enumerate(mhs):
171  fits_fn = asmb_refined_input.get_component_header(
172  i).get_transformations_fn()
173 
174  # todo - get the initial transformation
175  rb_ref = rbs_ref[i]
176  rb = rbs[i]
177 
179  rb_ref.get_reference_frame(),
180  rb.get_reference_frame())
181 
182  pdb_fn = asmb_input.get_component_header(i).get_filename()
183 
184  f = Fitter(
185  em_map, spacing, resolution, origin, asmb_input.get_assembly_header(
186  ).get_threshold(
187  ), pdb_fn, fits_fn, options.angle, options.num, options.angle_voxel,
188  options.max_trans, options.max_angle)
189  f.run_local_fitting(mh, rb, initial_transformation)
190 
191 
192 def main():
193  args = parse_args()
194  run(args.assembly_file, args.ref_assembly_file, args.proteomics_file,
195  args.mapping_file, args.combinations_file, args.combination_index,
196  args)
197 
198 if __name__ == "__main__":
199  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:86
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.
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.
def __init__
input a list of particles, the slope and theta of the sigmoid potential theta is the cutoff distance ...
IntsList read_paths(const char *txt_filename, int max_paths=INT_MAX)
Read paths.
VectorD< 3 > Vector3D
Definition: VectorD.h:425
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.