IMP logo
IMP Reference Guide  2.6.1
The Integrative Modeling Platform
system_tools.py
1 from __future__ import print_function, division
2 import IMP
3 import IMP.atom
5 import IMP.pmi
6 import IMP.pmi.tools
7 from collections import defaultdict
8 from math import pi
9 import os
10 
11 def resnums2str(res):
12  """Take iterable of TempResidues and return compatified string"""
13  if len(res)==0:
14  return ''
15  idxs = [r.get_index() for r in res]
16  idxs.sort()
17  all_ranges=[]
18  cur_range=[idxs[0],idxs[0]]
19  for idx in idxs[1:]:
20  if idx!=cur_range[1]+1:
21  all_ranges.append(cur_range)
22  cur_range=[idx,idx]
23  cur_range[1]=idx
24  all_ranges.append(cur_range)
25  ret = ''
26  for nr,r in enumerate(all_ranges):
27  ret+='%i-%i'%(r[0],r[1])
28  if nr<len(all_ranges)-1:
29  ret+=', '
30  return ret
31 
32 def get_structure(mdl,pdb_fn,chain_id,res_range=None,offset=0,model_num=None,ca_only=False):
33  """read a structure from a PDB file and return a list of residues
34  @param mdl The IMP model
35  @param pdb_fn The file to read
36  @param chain_id Chain ID to read
37  @param res_range Add only a specific set of residues. None gets you all
38  @param offset Apply an offset to the residue indexes of the PDB file
39  @param model_num Read multi-model PDB and return that model
40  @param ca_only Read only CA atoms (by default, all non-waters are read)
41  """
43  if ca_only:
45  if model_num is None:
46  mh = IMP.atom.read_pdb(pdb_fn,mdl,
48 
49  else:
50  mhs = IMP.atom.read_multimodel_pdb(pdb_fn,mdl,sel)
51  if model_num>=len(mhs):
52  raise Exception("you requested model num "+str(model_num)+\
53  " but the PDB file only contains "+str(len(mhs))+" models")
54  mh = IMP.atom.Selection(mhs[model_num],chain=chain_id,with_representation=True)
55 
56  if res_range==[] or res_range is None:
57  sel = IMP.atom.Selection(mh,chain=chain_id,atom_type=IMP.atom.AtomType('CA'))
58  else:
59  start = res_range[0]
60  end = res_range[1]
61  if end==-1:
62  end = IMP.atom.Residue(mh.get_children()[0].get_children()[-1]).get_index()
63  sel = IMP.atom.Selection(mh,chain=chain_id,residue_indexes=range(start,end+1),
64  atom_type=IMP.atom.AtomType('CA'))
65  ret=[]
66 
67  # return final and apply offset
68  for p in sel.get_selected_particles():
69  res = IMP.atom.Residue(IMP.atom.Atom(p).get_parent())
70  res.set_index(res.get_index() + offset)
71  ret.append(res)
72  return ret
73 
74 def build_bead(mdl,residues,input_coord=None):
75  """Generates a single bead"""
76 
77  ds_frag = (residues[0].get_index(), residues[-1].get_index())
78  prt = IMP.Particle(mdl)
80  ptem = IMP.core.XYZR(prt)
81  mass = IMP.atom.get_mass_from_number_of_residues(len(residues))
82 
83  if ds_frag[0] == ds_frag[-1]:
84  rt = residues[0].get_residue_type()
85  h = IMP.atom.Residue.setup_particle(prt, rt, ds_frag[0])
86  h.set_name('%i_bead' % (ds_frag[0]))
87  prt.set_name('%i_bead' % (ds_frag[0]))
88  try:
90  except IMP.ValueException:
92  IMP.atom.ResidueType("ALA"))
94  ptem.set_radius(radius)
95  else:
97  h.set_name('%i-%i_bead' % (ds_frag[0], ds_frag[-1]))
98  prt.set_name('%i-%i_bead' % (ds_frag[0], ds_frag[-1]))
99  h.set_residue_indexes(range(ds_frag[0], ds_frag[-1] + 1))
100  volume = IMP.atom.get_volume_from_mass(mass)
101  radius = 0.8 * (3.0 / 4.0 / pi * volume) ** (1.0 / 3.0)
102  ptem.set_radius(radius)
103 
105  try:
106  if tuple(input_coord) is not None:
107  ptem.set_coordinates(input_coord)
108  except TypeError:
109  pass
110  return h
111 
112 def build_necklace(mdl,residues, resolution, input_coord=None):
113  """Generates a string of beads with given length"""
114  out_hiers = []
115  for chunk in list(IMP.pmi.tools.list_chunks_iterator(residues, resolution)):
116  out_hiers.append(build_bead(mdl,chunk, input_coord=input_coord))
117  return out_hiers
118 
119 def build_ca_centers(mdl,residues):
120  """Create a bead on the CA position with coarsened size and mass"""
121  out_hiers = []
122  for tempres in residues:
123  residue = tempres.get_hierarchy()
124  rp1 = IMP.Particle(mdl)
125  rp1.set_name("Residue_%i"%residue.get_index())
126  rt = residue.get_residue_type()
127  this_res = IMP.atom.Residue.setup_particle(rp1,residue)
128  try:
130  except IMP.ValueException:
132  IMP.atom.ResidueType("ALA"))
133  try:
134  mass = IMP.atom.get_mass(rt)
135  except:
137  calpha = IMP.atom.Selection(residue,atom_type=IMP.atom.AT_CA). \
138  get_selected_particles()[0]
140  shape = IMP.algebra.Sphere3D(IMP.core.XYZ(calpha).get_coordinates(),radius)
143  out_hiers.append(this_res)
144  return out_hiers
145 
146 def setup_bead_as_gaussian(mh):
147  """Setup bead as spherical gaussian, using radius as variance"""
148  p = mh.get_particle()
149  center = IMP.core.XYZ(p).get_coordinates()
150  rad = IMP.core.XYZR(p).get_radius()
151  mass = IMP.atom.Mass(p).get_mass()
155 
156 
157 def show_representation(node):
158  print(node)
160  repr = IMP.atom.Representation(node)
161  resolutions = repr.get_resolutions()
162  for r in resolutions:
163  print('---- resolution %i ----' %r)
164  IMP.atom.show_molecular_hierarchy(repr.get_representation(r))
165  return True
166  else:
167  return False
168 
169 def recursive_show_representations(root):
170  pass
171 
172 def build_representation(parent,rep,coord_finder):
173  """Create requested representation.
174  For beads, identifies continuous segments and sets up as Representation.
175  If any volume-based representations (e.g.,densities) are requested,
176  will instead create a single Representation node.
177  All reps are added as children of the passed parent.
178  @param parent The Molecule to which we'll add add representations
179  @param rep What to build. An instance of pmi::topology::_Representation
180  @param coord_finder A _FindCloseStructure object to help localize beads
181  """
182  ret = []
183  atomic_res = 0
184  ca_res = 1
185  mdl = parent.get_model()
186  if rep.color is not None:
187  if type(rep.color) is float:
188  color = IMP.display.get_rgb_color(rep.color)
189  elif hasattr(rep.color,'__iter__') and len(rep.color)==3:
190  color = IMP.display.Color(*rep.color)
191  elif type(rep.color) is IMP.display.Color:
192  color = rep.color
193  else:
194  raise Exception("Color must be float or (r,g,b) tuple")
195  else:
196  color = None
197 
198  # first get the primary representation (currently, the smallest bead size)
199  # eventually we won't require beads to be present at all
200  primary_resolution = min(rep.bead_resolutions)
201 
202  # if collective densities, will return single node with everything
203  # below we sample or read the GMMs and add them as representation
204  single_node = False
205  if rep.density_residues_per_component!=None:
206  single_node = True
207  num_components = len(rep.residues)//rep.density_residues_per_component+1
208  rep_dict = defaultdict(list)
209  root_representation = IMP.atom.Representation.setup_particle(IMP.Particle(mdl),
210  primary_resolution)
211  parent.add_child(root_representation)
212  res_nums = [r.get_index() for r in rep.residues]
213  density_frag = IMP.atom.Fragment.setup_particle(IMP.Particle(mdl),res_nums)
214  density_frag.get_particle().set_name("Densities %i"%rep.density_residues_per_component)
215  density_ps = []
216 
217  if os.path.exists(rep.density_prefix+'.txt') and not rep.density_force_compute:
218  IMP.isd.gmm_tools.decorate_gmm_from_text(rep.density_prefix+'.txt',
219  density_ps,
220  mdl)
221  if len(density_ps)!=num_components or not os.path.exists(rep.density_prefix+'.txt') or rep.density_force_compute:
222  density_ps = []
223  fit_coords = []
224  for r in rep.residues:
225  fit_coords += [IMP.core.XYZ(p).get_coordinates() for p in IMP.core.get_leaves(r.hier)]
227  num_components,
228  mdl,
229  density_ps)
230  IMP.isd.gmm_tools.write_gmm_to_text(density_ps,rep.density_prefix+'.txt')
231  if rep.density_voxel_size>0.0:
232  IMP.isd.gmm_tools.write_gmm_to_map(density_ps,rep.density_prefix+'.mrc',rep.density_voxel_size)
233 
234  for d in density_ps:
235  density_frag.add_child(d)
236  root_representation.add_representation(density_frag,
237  IMP.atom.DENSITIES,
238  rep.density_residues_per_component)
239 
240  # get continuous segments from residues
241  segments = []
242  rsort = sorted(list(rep.residues),key=lambda r:r.get_index())
243  prev_idx = rsort[0].get_index()-1
244  prev_structure = rsort[0].get_has_structure()
245  cur_seg = []
246  for nr,r in enumerate(rsort):
247  if r.get_index()!=prev_idx+1 or r.get_has_structure()!=prev_structure or \
248  r.get_index() in rep.bead_extra_breaks:
249  segments.append(cur_seg)
250  cur_seg = []
251  cur_seg.append(r)
252  prev_idx = r.get_index()
253  prev_structure = r.get_has_structure()
254  if cur_seg!=[]:
255  segments.append(cur_seg)
256 
257  # for each segment, merge into beads
258  name_all = 'frags:'
259  name_count = 0
260  for frag_res in segments:
261  res_nums = [r.get_index() for r in frag_res]
262  rrange = "%i-%i"%(res_nums[0],res_nums[-1])
263  name = "Frag_"+rrange
264  if name_count<3:
265  name_all +=rrange+','
266  elif name_count==3:
267  name_all +='...'
268  name_count+=1
269  segp = IMP.Particle(mdl,name)
270  this_segment = IMP.atom.Fragment.setup_particle(segp,res_nums)
271  if not single_node:
272  this_representation = IMP.atom.Representation.setup_particle(segp,primary_resolution)
273  parent.add_child(this_representation)
274  for resolution in rep.bead_resolutions:
275  fp = IMP.Particle(mdl)
276  this_resolution = IMP.atom.Fragment.setup_particle(fp,res_nums)
277  this_resolution.set_name("%s: Res %i"%(name,resolution))
278  if frag_res[0].get_has_structure():
279  # if structured, merge particles as needed
280  if resolution==atomic_res:
281  for residue in frag_res:
282  this_resolution.add_child(residue.get_hierarchy())
283  elif resolution==ca_res and rep.bead_ca_centers:
284  beads = build_ca_centers(mdl,frag_res)
285  for bead in beads:
286  this_resolution.add_child(bead)
287  else:
289  for residue in frag_res:
290  tempc.add_child(IMP.atom.create_clone(residue.hier))
291  beads = IMP.atom.create_simplified_along_backbone(tempc,resolution)
292  for bead in beads.get_children():
293  this_resolution.add_child(bead)
294  del tempc
295  del beads
296  else:
297  # if unstructured, create necklace
298  input_coord = coord_finder.find_nearest_coord(min(r.get_index() for r in frag_res))
299  if input_coord is None:
300  input_coord = rep.bead_default_coord
301  beads = build_necklace(mdl,
302  frag_res,
303  resolution,
304  input_coord)
305  for bead in beads:
306  this_resolution.add_child(bead)
307 
308  # if requested, color all resolutions the same
309  if color:
310  for lv in IMP.core.get_leaves(this_resolution):
312 
313  # finally decide where to put this resolution
314  # if volumetric, collect resolutions from different segments together
315  if single_node:
316  rep_dict[resolution]+=this_resolution.get_children()
317  else:
318  if resolution==primary_resolution:
319  this_representation.add_child(this_resolution)
320  else:
321  this_representation.add_representation(this_resolution,
322  IMP.atom.BALLS,
323  resolution)
324  # if individual beads to be setup as Gaussians:
325  if rep.setup_particles_as_densities:
326  for p in IMP.core.get_leaves(this_resolution):
327  setup_bead_as_gaussian(p)
328  this_resolution.set_name(this_resolution.get_name()+' Densities %i'%resolution)
329  this_representation.add_representation(this_resolution,
330  IMP.atom.DENSITIES,
331  resolution)
332 
333  if single_node:
334  root_representation.set_name(name_all.strip(',')+": Base")
335  d = root_representation.get_representations(IMP.atom.DENSITIES)
336  d[0].set_name('%s: '%name_all + d[0].get_name())
337  for resolution in rep.bead_resolutions:
338  this_resolution = IMP.atom.Fragment.setup_particle(
339  IMP.Particle(mdl),
340  [r.get_index() for r in rep.residues])
341  this_resolution.set_name("%s: Res %i"%(name_all,resolution))
342  for hier in rep_dict[resolution]:
343  this_resolution.add_child(hier)
344  if resolution==primary_resolution:
345  root_representation.add_child(this_resolution)
346  else:
347  root_representation.add_representation(this_resolution,
348  IMP.atom.BALLS,
349  resolution)
def list_chunks_iterator
Yield successive length-sized chunks from a list.
Definition: tools.py:1115
Tools for handling Gaussian Mixture Models.
Definition: gmm_tools.py:1
Add mass to a particle.
Definition: Mass.h:23
double get_volume_from_residue_type(ResidueType rt)
Return an estimate for the volume of a given residue.
Simple 3D transformation class.
Represent an RGB color.
Definition: Color.h:24
static Gaussian setup_particle(Model *m, ParticleIndex pi)
Definition: Gaussian.h:48
void show_molecular_hierarchy(Hierarchy h)
Print out the molecular hierarchy.
static Fragment setup_particle(Model *m, ParticleIndex pi)
Definition: Fragment.h:63
double get_mass(const Selection &s)
Get the total mass of a hierarchy, in Daltons.
static XYZR setup_particle(Model *m, ParticleIndex pi)
Definition: XYZR.h:48
Select atoms which are selected by both selectors.
Definition: pdb.h:326
double get_mass(ResidueType c)
Get the mass from the residue type.
Color get_rgb_color(double f)
Return the color for f from the RGB color map.
double get_mass_from_number_of_residues(unsigned int num_aa)
Estimate the mass of a protein from the number of amino acids.
Miscellaneous utilities.
Definition: tools.py:1
double get_ball_radius_from_volume_3d(double volume)
Return the radius of a sphere with a given volume.
Definition: Sphere3D.h:35
The type of an atom.
static Residue setup_particle(Model *m, ParticleIndex pi, ResidueType t, int index, int insertion_code)
Definition: Residue.h:157
static Representation setup_particle(Model *m, ParticleIndex pi)
GenericHierarchies get_leaves(Hierarchy mhd)
Get all the leaves of the bit of hierarchy.
void read_pdb(TextInput input, int model, Hierarchy h)
A reference frame in 3D.
A Gaussian distribution in 3D.
Definition: Gaussian3D.h:24
def fit_gmm_to_points
fit a GMM to some points.
Definition: gmm_tools.py:224
A decorator for a representation.
double get_volume_from_mass(double m, ProteinDensityReference ref=ALBER)
Estimate the volume of a protein from its mass.
Ints get_index(const ParticlesTemp &particles, const Subset &subset, const Subsets &excluded)
A decorator for a particle representing an atom.
Definition: atom/Atom.h:234
static Mass setup_particle(Model *m, ParticleIndex pi, Float mass)
Definition: Mass.h:44
The type for a residue.
PDBSelector * get_default_pdb_selector()
Definition: pdb.h:429
A decorator for a particle with x,y,z coordinates.
Definition: XYZ.h:30
static Colored setup_particle(Model *m, ParticleIndex pi, Color color)
Definition: Colored.h:62
def write_gmm_to_map
write density map from GMM.
Definition: gmm_tools.py:80
A decorator for a residue.
Definition: Residue.h:134
static bool get_is_setup(const IMP::ParticleAdaptor &p)
Hierarchies read_multimodel_pdb(TextInput input, Model *model, PDBSelector *selector=get_default_pdb_selector())
Hierarchy create_simplified_along_backbone(Chain input, const IntRanges &residue_segments, bool keep_detailed=false)
Rotation3D get_identity_rotation_3d()
Return a rotation that does not do anything.
Definition: Rotation3D.h:236
Class to handle individual model particles.
Definition: Particle.h:37
Select all CA ATOM records.
Definition: pdb.h:77
Python classes to represent, score, sample and analyze models.
def write_gmm_to_text
write a list of gaussians to text.
Definition: gmm_tools.py:62
Functionality for loading, creating, manipulating and scoring atomic structures.
static Chain setup_particle(Model *m, ParticleIndex pi, std::string id)
Definition: Chain.h:41
An exception for an invalid value being passed to IMP.
Definition: exception.h:137
Select hierarchy particles identified by the biological name.
Definition: Selection.h:66
Select all ATOM and HETATM records with the given chain ids.
Definition: pdb.h:189
def decorate_gmm_from_text
read the output from write_gmm_to_text, decorate as Gaussian and Mass
Definition: gmm_tools.py:22
A decorator for a particle with x,y,z coordinates and a radius.
Definition: XYZR.h:27