1 from __future__
import print_function
7 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(cluster_conform_i, cluster_conform_0, align, ps, trans):
17 def _to_vector3ds(numpy_array):
23 calculator = pyRMSD.RMSDCalculator.RMSDCalculator(
"QCP_SERIAL_CALCULATOR", numpy.array([cluster_conform_0, cluster_conform_i]))
25 calculator = pyRMSD.RMSDCalculator.RMSDCalculator(
"NOSUP_SERIAL_CALCULATOR", numpy.array([cluster_conform_0, cluster_conform_i]))
27 rmsd, superposed_fit = calculator.pairwise(0, 1, get_superposed_coordinates =
True)
34 _to_vector3ds(superposed_fit[0]), _to_vector3ds(cluster_conform_0))
36 for particle_index
in range(len(superposed_fit[1])):
40 return rmsd, ps, trans
42 def get_particles_from_superposed_amb(cluster_conform_i, cluster_conform_0, align, ps, trans, symm_groups):
44 '''Modified superposed function to work with symmetric copies'''
46 def _to_vector3ds(numpy_array):
53 superposed_final_coords=[]
55 for perm
in pyRMSD.symmTools.symm_permutations(symm_groups):
57 new_cluster_conform_i = cluster_conform_i
61 for [particle0, particle1]
in sg:
63 if particle0 > particle1:
65 pyRMSD.symmTools.swap_atoms(new_cluster_conform_i, particle0, particle1)
68 calculator = pyRMSD.RMSDCalculator.RMSDCalculator(
"QCP_SERIAL_CALCULATOR", numpy.array([cluster_conform_0, new_cluster_conform_i]))
70 calculator = pyRMSD.RMSDCalculator.RMSDCalculator(
"NOSUP_SERIAL_CALCULATOR", numpy.array([cluster_conform_0, new_cluster_conform_i]))
72 rmsd, superposed_fit = calculator.pairwise(0, 1, get_superposed_coordinates =
True)
76 superposed_final_coords = superposed_fit
85 for particle_index
in range(len(superposed_final_coords[1])):
90 return min_rmsd, ps, trans
93 """Compute mean density maps from structures.
94 Keeps a dictionary of density maps,
95 keys are in the custom ranges. When you call add_subunits_density, it adds
96 particle coordinates to the existing density maps.
99 def __init__(self, custom_ranges=None, resolution=20.0, voxel=5.0, bead_names = None):
101 @param list of particles decorated with mass, radius, and XYZ
102 @param resolution The MRC resolution of the output map (in Angstrom unit)
103 @param voxel The voxel size for the output map (lower is slower)
106 self.MRCresolution = resolution
108 self.count_models = 0.0
110 self.bead_names = bead_names
111 self.custom_ranges=custom_ranges
114 self.particle_indices_in_custom_ranges={}
116 for density_name
in self.custom_ranges:
117 self.particle_indices_in_custom_ranges[density_name]=[]
120 for index,beadname
in enumerate(self.bead_names):
121 for density_name
in self.custom_ranges:
122 for domain
in self.custom_ranges[density_name]:
123 if self._is_contained(beadname,domain):
124 self.particle_indices_in_custom_ranges[density_name].append(index)
128 def normalize_density(self):
131 def _create_density_from_particles(self, ps, name,
132 kernel_type=
'GAUSSIAN'):
133 '''Internal function for adding to densities.
134 pass XYZR particles with mass and create a density from them.
135 kernel type options are GAUSSIAN, BINARIZED_SPHERE, and SPHERE.'''
137 'GAUSSIAN': IMP.em.GAUSSIAN,
138 'BINARIZED_SPHERE': IMP.em.BINARIZED_SPHERE,
139 'SPHERE': IMP.em.SPHERE}
142 dmap.set_was_used(
True)
144 if name
not in self.densities:
145 self.densities[name] = dmap
151 dmap3.set_was_used(
True)
153 dmap3.add(self.densities[name])
154 self.densities[name] = dmap3
156 def _is_contained(self,bead_name,domain):
157 """ domain can be the name of a single protein or a tuple (start_residue,end_residue,protein_name)
158 bead is a string of type moleculeName_startResidue_endResidue
161 (bead_protein, bead_res_start,
162 bead_res_end, bead_copy) = bead_name.split(
"_")
165 if isinstance(domain, tuple):
166 domain_protein = domain[2]
168 domain_protein = domain
170 if "." in domain_protein:
171 spl = domain_protein.split(
".")
172 domain_protein = spl[0]
173 domain_copy = int(spl[1])
175 domain_copy = bead_copy = -1
177 if bead_protein != domain_protein
or int(bead_copy) != domain_copy:
181 if isinstance(domain, tuple):
182 bead_residues = set(range(int(bead_res_start),int(bead_res_end)+1))
183 domain_residues = set(range(int(domain[0]),int(domain[1])+1))
184 return not domain_residues.isdisjoint(bead_residues)
189 """Add a frame to the densities.
190 @param ps List of particles decorated with XYZR and Mass.
192 self.count_models += 1.0
194 particles_custom_ranges={}
195 for density_name
in self.custom_ranges:
196 particles_custom_ranges[density_name]=[]
199 for density_name
in self.custom_ranges:
200 for particle_index
in self.particle_indices_in_custom_ranges[density_name]:
201 particles_custom_ranges[density_name].append(ps[particle_index])
204 for density_name
in self.custom_ranges:
205 self._create_density_from_particles(particles_custom_ranges[density_name],density_name)
207 def get_density_keys(self):
208 return list(self.densities.keys())
211 """Get the current density for some component name"""
212 if name
not in self.densities:
215 return self.densities[name]
217 def write_mrc(self, path=".", file_prefix=""):
218 for density_name
in self.densities:
219 mrc = os.path.join(path, file_prefix +
"_" + density_name +
".mrc")
220 self.densities[density_name].
multiply(1. / self.count_models)
222 self.densities[density_name], mrc,
224 if len(self.densities) == 1:
227 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)