1 from __future__
import print_function
6 import pyRMSD.RMSDCalculator
9 def parse_custom_ranges(ranges_file):
10 with open(ranges_file)
as fh:
13 return d[
'density_custom_ranges']
16 def get_particles_from_superposed(
17 cluster_conform_i, cluster_conform_0, align, ps, trans):
18 def _to_vector3ds(numpy_array):
24 calculator = pyRMSD.RMSDCalculator.RMSDCalculator(
25 "QCP_SERIAL_CALCULATOR",
26 numpy.array([cluster_conform_0, cluster_conform_i]))
28 calculator = pyRMSD.RMSDCalculator.RMSDCalculator(
29 "NOSUP_SERIAL_CALCULATOR",
30 numpy.array([cluster_conform_0, cluster_conform_i]))
32 rmsd, superposed_fit = calculator.pairwise(
33 0, 1, get_superposed_coordinates=
True)
40 _to_vector3ds(superposed_fit[0]), _to_vector3ds(cluster_conform_0))
42 for particle_index
in range(len(superposed_fit[1])):
47 return rmsd, ps, trans
50 def get_particles_from_superposed_amb(
51 cluster_conform_i, cluster_conform_0, align, ps, trans, symm_groups):
53 '''Modified superposed function to work with symmetric copies'''
55 def _to_vector3ds(numpy_array):
62 superposed_final_coords = []
64 for perm
in pyRMSD.symmTools.symm_permutations(symm_groups):
67 new_cluster_conform_i = cluster_conform_i
71 for [particle0, particle1]
in sg:
74 if particle0 > particle1:
75 pyRMSD.symmTools.swap_atoms(
76 new_cluster_conform_i, particle0, particle1)
79 calculator = pyRMSD.RMSDCalculator.RMSDCalculator(
80 "QCP_SERIAL_CALCULATOR",
81 numpy.array([cluster_conform_0, new_cluster_conform_i]))
83 calculator = pyRMSD.RMSDCalculator.RMSDCalculator(
84 "NOSUP_SERIAL_CALCULATOR",
85 numpy.array([cluster_conform_0, new_cluster_conform_i]))
87 rmsd, superposed_fit = calculator.pairwise(
88 0, 1, get_superposed_coordinates=
True)
92 superposed_final_coords = superposed_fit
100 _to_vector3ds(superposed_final_coords[0]),
101 _to_vector3ds(cluster_conform_0))
103 for particle_index
in range(len(superposed_final_coords[1])):
108 superposed_final_coords[1][particle_index]))
110 return min_rmsd, ps, trans
114 """Compute mean density maps from structures.
115 Keeps a dictionary of density maps,
116 keys are in the custom ranges. When you call add_subunits_density, it adds
117 particle coordinates to the existing density maps.
120 def __init__(self, custom_ranges=None, resolution=20.0, voxel=5.0,
123 @param list of particles decorated with mass, radius, and XYZ
124 @param resolution The MRC resolution of the output map
126 @param voxel The voxel size for the output map (lower is slower)
129 self.MRCresolution = resolution
131 self.count_models = 0.0
133 self.bead_names = bead_names
134 self.custom_ranges = custom_ranges
138 self.particle_indices_in_custom_ranges = {}
140 for density_name
in self.custom_ranges:
141 self.particle_indices_in_custom_ranges[density_name] = []
144 for index, beadname
in enumerate(self.bead_names):
145 for density_name
in self.custom_ranges:
147 for domain
in self.custom_ranges[density_name]:
148 if self._is_contained(beadname, domain):
149 self.particle_indices_in_custom_ranges[
150 density_name].append(index)
153 def normalize_density(self):
156 def _create_density_from_particles(self, ps, name,
157 kernel_type=
'GAUSSIAN'):
158 '''Internal function for adding to densities.
159 pass XYZR particles with mass and create a density from them.
160 kernel type options are GAUSSIAN, BINARIZED_SPHERE, and SPHERE.'''
163 dmap.set_was_used(
True)
165 if name
not in self.densities:
166 self.densities[name] = dmap
172 dmap3.set_was_used(
True)
174 dmap3.add(self.densities[name])
175 self.densities[name] = dmap3
177 def _is_contained(self, bead_name, domain):
178 """ domain can be the name of a single protein or a tuple
179 (start_residue,end_residue,protein_name)
180 bead is a string of type moleculeName_startResidue_endResidue
183 (bead_protein, bead_res_start,
184 bead_res_end, bead_copy) = bead_name.split(
"_")
187 if isinstance(domain, tuple):
188 domain_protein = domain[2]
190 domain_protein = domain
192 if "." in domain_protein:
193 spl = domain_protein.split(
".")
194 domain_protein = spl[0]
195 domain_copy = int(spl[1])
197 domain_copy = bead_copy = -1
199 if bead_protein != domain_protein
or int(bead_copy) != domain_copy:
203 if isinstance(domain, tuple):
204 bead_residues = set(range(int(bead_res_start),
205 int(bead_res_end)+1))
206 domain_residues = set(range(int(domain[0]),
208 return not domain_residues.isdisjoint(bead_residues)
213 """Add a frame to the densities.
214 @param ps List of particles decorated with XYZR and Mass.
216 self.count_models += 1.0
218 particles_custom_ranges = {}
219 for density_name
in self.custom_ranges:
220 particles_custom_ranges[density_name] = []
223 for density_name
in self.custom_ranges:
225 in self.particle_indices_in_custom_ranges[density_name]:
226 particles_custom_ranges[density_name].append(
230 for density_name
in self.custom_ranges:
231 self._create_density_from_particles(
232 particles_custom_ranges[density_name], density_name)
234 def get_density_keys(self):
235 return list(self.densities.keys())
238 """Get the current density for some component name"""
239 if name
not in self.densities:
242 return self.densities[name]
244 def write_mrc(self, path=".", file_prefix=""):
245 for density_name
in self.densities:
246 mrc = os.path.join(path, file_prefix +
"_" + density_name +
".mrc")
247 self.densities[density_name].
multiply(1. / self.count_models)
249 self.densities[density_name], mrc,
251 if len(self.densities) == 1:
254 return os.path.join(path, file_prefix +
"_*.mrc")
def get_density
Get the current density for some component name.
Compute mean density maps from structures.
Class for sampling a density map from particles.
DensityMap * multiply(const DensityMap *m1, const DensityMap *m2)
Return a density map for which voxel i contains the result of m1[i]*m2[i].
DensityMap * create_density_map(const IMP::algebra::GridD< 3, S, V, E > &arg)
Create a density map from an arbitrary IMP::algebra::GridD.
def add_subunits_density
Add a frame to the densities.
Basic utilities for handling cryo-electron microscopy 3D density maps.
A decorator for a particle with x,y,z coordinates.
algebra::BoundingBoxD< 3 > get_bounding_box(const DensityMap *m)
Transformation3D get_transformation_aligning_first_to_second(Vector3Ds a, Vector3Ds b)