IMP logo
IMP Reference Guide  2.20.0
The Integrative Modeling Platform
plotting/topology.py
1 from __future__ import print_function
2 import IMP
3 import RMF
4 import IMP.atom
5 import IMP.rmf
6 import IMP.pmi
7 import IMP.pmi.analysis
8 import IMP.pmi.output
9 from collections import defaultdict
10 
11 
12 class TopologyPlot(object):
13  """A class to read RMF files and make a network contact map"""
14  def __init__(self, model, selections, cutoff, frequency_cutoff,
15  colors=None, fixed=None, pos=None, proteomic_edges=None,
16  quantitative_proteomic_data=None):
17  """Set up a new graphXL object
18  @param model The IMP model
19  @param selections A dictionary containing component names.
20  Keys are labels
21  values are either moleculename or
22  start,stop,moleculename
23  @param cutoff The distance cutoff
24  @param frequency_cutoff The frequency cutoff
25  @param colors A dictionary of colors (HEX code,values)
26  for subunits (keywords)
27  @param fixed A list of subunits that are kept fixed
28  @param pos A dictionary with positions (tuple, values)
29  of subunits (keywords)
30  @param proteomic_edges A list edges to represent proteomic data
31  @param quantitative_proteomic_data A dictionary of edges to represent
32  quantitative proteomic data such as PE Scores,
33  or genetic interactions
34  """
35  import itertools
36 
37  self.model = model
38  self.selections = selections
39  self.contact_counts = {}
40  self.edges = defaultdict(int)
41  for (name1, name2) in itertools.combinations(self.selections.keys(),
42  2):
43  self.edges[tuple(sorted((name1, name2)))] = 0
44  self.cutoff = cutoff
45  self.frequency_cutoff = frequency_cutoff
46  self.gcpf = IMP.core.GridClosePairsFinder()
47  self.gcpf.set_distance(self.cutoff)
48  self.names = list(self.selections.keys())
49  self.colors = colors
50  self.fixed = fixed
51  self.pos = pos
52  self.proteomic_edges = proteomic_edges
53  self.quantitative_proteomic_data = quantitative_proteomic_data
54  self.num_rmf = 0
55 
56  def add_rmf(self, rmf_fn, nframe):
57  """Add selections from an RMF file"""
58  print('reading from RMF file', rmf_fn)
59  rh = RMF.open_rmf_file_read_only(rmf_fn)
60  prots = IMP.rmf.create_hierarchies(rh, self.model)
61  hier = prots[0]
62  IMP.rmf.load_frame(rh, RMF.FrameID(0))
63  ps_per_component = defaultdict(list)
64  if self.num_rmf == 0:
65  self.size_per_component = defaultdict(int)
66  self.model.update()
67 
68  # gathers particles for all components
70  all_particles_by_resolution = []
71  for name in part_dict:
72  all_particles_by_resolution += part_dict[name]
73 
74  for component_name in self.selections:
75  for seg in self.selections[component_name]:
76  if isinstance(seg, str):
77  s = IMP.atom.Selection(hier, molecule=seg)
78  elif isinstance(seg, tuple):
80  hier, molecule=seg[2],
81  residue_indexes=range(seg[0], seg[1] + 1))
82  else:
83  raise Exception('could not understand selection tuple '
84  + str(seg))
85  parts = list(set(s.get_selected_particles())
86  & set(all_particles_by_resolution))
87  ps_per_component[component_name] += IMP.get_indexes(parts)
88  if self.num_rmf == 0:
89  self.size_per_component[component_name] += \
91  for p in parts)
92 
93  for n1, name1 in enumerate(self.names):
94  for name2 in self.names[n1+1:]:
95  ncontacts = len(self.gcpf.get_close_pairs(
96  self.model, ps_per_component[name1],
97  ps_per_component[name2]))
98  if ncontacts > 0:
99  self.edges[tuple(sorted((name1, name2)))] += 1.0
100 
101  self.num_rmf += 1
102 
103  def make_plot(self, groups, out_fn, quantitative_proteomic_data=False):
104  '''
105  plot the interaction matrix
106  @param groups is the list of groups of domains, eg,
107  [["protA_1-10","prot1A_11-100"],["protB"]....]
108  it will plot a space between different groups
109  @param out_fn name of the plot file
110  @param quantitative_proteomic_data plot the quantitative proteomic data
111  '''
112  import matplotlib.pyplot as plt
113 
114  ax = plt.gca()
115  ax.set_aspect('equal', 'box')
116  ax.xaxis.set_major_locator(plt.NullLocator())
117  ax.yaxis.set_major_locator(plt.NullLocator())
118 
119  largespace = 0.6
120  smallspace = 0.5
121  squaredistance = 1.0
122  squaresize = 0.99
123  domain_xlocations = {}
124  domain_ylocations = {}
125 
126  xoffset = squaredistance
127  yoffset = squaredistance
128  for group in groups:
129  xoffset += largespace
130  yoffset += largespace
131  for subgroup in group:
132  xoffset += smallspace
133  yoffset += smallspace
134  for domain in subgroup:
135  domain_xlocations[domain] = xoffset
136  domain_ylocations[domain] = yoffset
137  xoffset += squaredistance
138  yoffset += squaredistance
139 
140  for edge, count in self.edges.items():
141 
142  if quantitative_proteomic_data:
143  # normalize
144  maxqpd = max(self.quantitative_proteomic_data.values())
145  minqpd = min(self.quantitative_proteomic_data.values())
146  if edge in self.quantitative_proteomic_data:
147  value = self.quantitative_proteomic_data[edge]
148  elif (edge[1], edge[0]) in self.quantitative_proteomic_data:
149  value = self.quantitative_proteomic_data[(edge[1],
150  edge[0])]
151  else:
152  value = 0.0
153  print(minqpd, maxqpd)
154  density = (1.0-(value-minqpd)/(maxqpd-minqpd))
155  else:
156  density = (1.0-float(count)/self.num_rmf)
157  color = (density, density, 1.0)
158  x = domain_xlocations[edge[0]]
159  y = domain_ylocations[edge[1]]
160  if x > y:
161  xtmp = y
162  ytmp = x
163  x = xtmp
164  y = ytmp
165  rect = plt.Rectangle([x - squaresize / 2, y - squaresize / 2],
166  squaresize, squaresize,
167  facecolor=color, edgecolor='Gray',
168  linewidth=0.1)
169  ax.add_patch(rect)
170  rect = plt.Rectangle([y - squaresize / 2, x - squaresize / 2],
171  squaresize, squaresize,
172  facecolor=color, edgecolor='Gray',
173  linewidth=0.1)
174  ax.add_patch(rect)
175 
176  ax.autoscale_view()
177  plt.savefig(out_fn)
178  plt.show()
179  exit()
180 
181  def make_graph(self, out_fn):
182  edges = []
183  weights = []
184  print('num edges', len(self.edges))
185  for edge, count in self.edges.items():
186  # filter if frequency of contacts is greater than frequency_cutoff
187  if float(count)/self.num_rmf > self.frequency_cutoff:
188  print(count, edge)
189  edges.append(edge)
190  weights.append(count)
191  for nw, w in enumerate(weights):
192  weights[nw] = float(weights[nw])/max(weights)
193  IMP.pmi.output.draw_graph(edges,
194  node_size=dict(self.size_per_component),
195  node_color=self.colors,
196  fixed=self.fixed,
197  pos=self.pos,
198  edge_thickness=1, # weights,
199  edge_alpha=0.3,
200  edge_color='gray',
201  out_filename=out_fn,
202  validation_edges=self.proteomic_edges)
203 
204 
205 def draw_component_composition(DegreesOfFreedom, max=1000,
206  draw_pdb_names=False, show=True):
207  """A function to plot the representation on the sequence
208  @param DegreesOfFreedom input a IMP.pmi.dof.DegreesOfFreedom instance"""
209  import matplotlib as mpl
210  mpl.use('Agg')
211  from matplotlib import pyplot
212  from operator import itemgetter
213  import IMP.pmi.tools
214 
215  # first build the movers dictionary
216  movers_mols_res = IMP.pmi.tools.OrderedDict()
217  for mv in DegreesOfFreedom.movers_particles_map:
218  hs = DegreesOfFreedom.movers_particles_map[mv]
219  res = []
220 
221  for h in hs:
223  res = [IMP.atom.Residue(h).get_index()]
225  res = list(IMP.atom.Fragment(h).get_residue_indexes())
226  if res:
227  is_molecule = False
228  hp = h
229  while not is_molecule:
230  hp = hp.get_parent()
231  is_molecule = IMP.atom.Molecule.get_is_setup(hp)
232  name = (IMP.atom.Molecule(hp).get_name() + "."
233  + str(IMP.atom.Copy(hp).get_copy_index()))
234  if mv not in movers_mols_res:
235  movers_mols_res[mv] = IMP.pmi.tools.OrderedDict()
236  movers_mols_res[mv][name] = IMP.pmi.tools.Segments(res)
237  else:
238  if name not in movers_mols_res[mv]:
239  movers_mols_res[mv][name] = IMP.pmi.tools.Segments(res)
240  else:
241  movers_mols_res[mv][name].add(res)
242  # then get all movers by type
243  fb_movers = []
244  rb_movers = []
245  srb_movers = []
246  for mv in movers_mols_res:
247  if type(mv) is IMP.core.RigidBodyMover:
248  rb_movers.append(mv)
249  if type(mv) is IMP.core.BallMover:
250  fb_movers.append(mv)
251  if type(mv) is IMP.pmi.TransformMover:
252  srb_movers.append(mv)
253 
254  # now remove residue assigned to BallMovers from RigidBodies
255  for mv in fb_movers:
256  for mol in movers_mols_res[mv]:
257  for i in movers_mols_res[mv][mol].get_flatten():
258  for mv_rb in rb_movers:
259  try:
260  movers_mols_res[mv_rb][mol].remove(i)
261  except (ValueError, KeyError):
262  continue
263 
264  elements = {}
265  rgb_colors = [(240, 163, 255), (0, 117, 220),
266  (153, 63, 0), (76, 0, 92), (25, 25, 25),
267  (0, 92, 49), (43, 206, 72), (255, 204, 153),
268  (128, 128, 128), (148, 255, 181), (143, 124, 0),
269  (157, 204, 0), (194, 0, 136), (0, 51, 128), (255, 164, 5),
270  (255, 168, 187), (66, 102, 0), (255, 0, 16), (94, 241, 242),
271  (0, 153, 143), (224, 255, 102), (116, 10, 255), (153, 0, 0),
272  (255, 255, 128), (255, 255, 0), (255, 80, 5)]
273  colors = [(float(c0)/255, float(c1)/255, float(c2)/255)
274  for (c0, c1, c2) in rgb_colors]
275  mvrb_color = {}
276  for mv in movers_mols_res:
277  if type(mv) is IMP.core.RigidBodyMover:
278  if mv not in mvrb_color:
279  color = colors[len(mvrb_color.keys())]
280  mvrb_color[mv] = color
281  for mol in movers_mols_res[mv]:
282  if mol not in elements:
283  elements[mol] = []
284  for seg in movers_mols_res[mv][mol].segs:
285  try:
286  elements[mol].append((seg[0], seg[-1], " ", "pdb",
287  mvrb_color[mv]))
288  except (KeyError, IndexError):
289  continue
290  if type(mv) is IMP.core.BallMover:
291  for mol in movers_mols_res[mv]:
292  if mol not in elements:
293  elements[mol] = []
294  for seg in movers_mols_res[mv][mol].segs:
295  elements[mol].append((seg[0], seg[-1], " ", "bead"))
296 
297  # sort everything
298  for mol in elements:
299  elements[mol].sort(key=lambda tup: tup[0])
300  elements[mol].append((elements[mol][-1][1],
301  elements[mol][-1][1], " ", "end"))
302 
303  for name in elements:
304  k = name
305  tmplist = sorted(elements[name], key=itemgetter(0))
306  try:
307  endres = tmplist[-1][1]
308  except IndexError:
309  print("draw_component_composition: missing information for "
310  "component %s" % name)
311  return
312  fig = pyplot.figure(figsize=(26.0 * float(endres) / max + 2, 2))
313  ax = fig.add_axes([0.05, 0.475, 0.9, 0.15])
314 
315  # Set the colormap and norm to correspond to the data for which
316  # the colorbar will be used.
317  cmap = mpl.cm.cool
318  norm = mpl.colors.Normalize(vmin=5, vmax=10)
319  bounds = [1]
320  colors = []
321 
322  for n, l in enumerate(tmplist):
323  if l[3] != "end":
324  if bounds[-1] != l[0]:
325  colors.append("white")
326  bounds.append(l[0])
327  if l[3] == "pdb":
328  colors.append("#99CCFF")
329  if l[3] == "bead":
330  colors.append("white")
331  if l[3] == "helix":
332  colors.append("#33CCCC")
333  if l[3] != "end":
334  bounds.append(l[1] + 1)
335  else:
336  if l[3] == "pdb":
337  colors.append(l[4])
338  if l[3] == "bead":
339  colors.append("white")
340  if l[3] == "helix":
341  colors.append("#33CCCC")
342  if l[3] != "end":
343  bounds.append(l[1] + 1)
344  else:
345  if bounds[-1] - 1 == l[0]:
346  bounds.pop()
347  bounds.append(l[0])
348  else:
349  colors.append("white")
350  bounds.append(l[0])
351 
352  bounds.append(bounds[-1])
353  colors.append("white")
354  cmap = mpl.colors.ListedColormap(colors)
355  cmap.set_over('0.25')
356  cmap.set_under('0.75')
357 
358  norm = mpl.colors.BoundaryNorm(bounds, cmap.N)
359  cb2 = mpl.colorbar.ColorbarBase(ax, cmap=cmap,
360  norm=norm,
361  # to use 'extend', you must
362  # specify two extra boundaries:
363  boundaries=bounds,
364  ticks=bounds, # optional
365  spacing='proportional',
366  orientation='horizontal')
367 
368  extra_artists = []
369  npdb = 0
370 
371  if draw_pdb_names:
372  for tl in tmplist:
373  if tl[3] == "pdb":
374  npdb += 1
375  mid = 1.0 / endres * float(tl[0])
376  t = ax.annotate(
377  tl[2], xy=(mid, 1), xycoords='axes fraction',
378  xytext=(mid + 0.025, float(npdb - 1) / 2.0 + 1.5),
379  textcoords='axes fraction',
380  arrowprops=dict(
381  arrowstyle="->",
382  connectionstyle="angle,angleA=0,angleB=90,rad=10"))
383  extra_artists.append(t)
384 
385  # set the title of the bar
386  title = ax.text(-0.005, 0.5, k, ha="right", va="center", rotation=90,
387  size=20)
388 
389  extra_artists.append(title)
390  # changing the xticks labels
391  labels = len(bounds) * [" "]
392  ax.set_xticklabels(labels)
393  mid = 1.0 / endres * float(bounds[0])
394  t = ax.annotate(bounds[0], xy=(mid, 0), xycoords='axes fraction',
395  xytext=(mid - 0.01, -0.5), textcoords='axes fraction',)
396  extra_artists.append(t)
397  offsets = [0, -0.5, -1.0]
398  nclashes = 0
399  for n in range(1, len(bounds)):
400  if bounds[n] == bounds[n - 1]:
401  continue
402  mid = 1.0 / endres * float(bounds[n])
403  if (float(bounds[n]) - float(bounds[n - 1])) / max <= 0.01:
404  nclashes += 1
405  offset = offsets[nclashes % 3]
406  else:
407  nclashes = 0
408  offset = offsets[0]
409  if offset > -0.75:
410  t = ax.annotate(bounds[n], xy=(mid, 0),
411  xycoords='axes fraction',
412  xytext=(mid, -0.5 + offset),
413  textcoords='axes fraction')
414  else:
415  t = ax.annotate(bounds[n], xy=(mid, 0),
416  xycoords='axes fraction',
417  xytext=(mid, -0.5 + offset),
418  textcoords='axes fraction',
419  arrowprops=dict(arrowstyle="-"))
420  extra_artists.append(t)
421 
422  cb2.add_lines(bounds, ["black"] * len(bounds), [1] * len(bounds))
423  # cb2.set_label(k)
424 
425  pyplot.savefig(k + "structure.pdf")
426  if show:
427  pyplot.show()
static bool get_is_setup(const IMP::ParticleAdaptor &p)
Definition: Molecule.h:35
A class to read RMF files and make a network contact map.
static bool get_is_setup(const IMP::ParticleAdaptor &p)
Definition: Residue.h:158
A decorator to associate a particle with a part of a protein/DNA/RNA.
Definition: Fragment.h:20
atom::Hierarchies create_hierarchies(RMF::FileConstHandle fh, Model *m)
Miscellaneous utilities.
Definition: tools.py:1
Modify the transformation of a rigid body.
def make_plot
plot the interaction matrix
This class stores integers in ordered compact lists eg: [[1,2,3],[6,7,8]] the methods help splitting ...
Definition: tools.py:623
Move continuous particle variables by perturbing them within a ball.
def add_rmf
Add selections from an RMF file.
static bool get_is_setup(Model *m, ParticleIndex pi)
Definition: Fragment.h:46
A decorator for keeping track of copies of a molecule.
Definition: Copy.h:28
ParticleIndexPairs get_indexes(const ParticlePairsTemp &ps)
Get the indexes from a list of particle pairs.
Ints get_index(const ParticlesTemp &particles, const Subset &subset, const Subsets &excluded)
def __init__
Set up a new graphXL object.
void load_frame(RMF::FileConstHandle file, RMF::FrameID frame)
Load the given RMF frame into the state of the linked objects.
def get_particles_at_resolution_one
Get particles at res 1, or any beads, based on the name.
Modify the transformation of a rigid body.
Tools for clustering and cluster analysis.
Definition: pmi/Analysis.py:1
Find all nearby pairs by testing all pairs.
Classes for writing output files and processing them.
Definition: output.py:1
A decorator for a residue.
Definition: Residue.h:137
The general base class for IMP exceptions.
Definition: exception.h:48
int get_copy_index(Hierarchy h)
Walk up the hierarchy to find the current copy index.
Python classes to represent, score, sample and analyze models.
Functionality for loading, creating, manipulating and scoring atomic structures.
A decorator for a molecule.
Definition: Molecule.h:24
Select hierarchy particles identified by the biological name.
Definition: Selection.h:70
Support for the RMF file format for storing hierarchical molecular data and markup.
def get_residue_indexes
Retrieve the residue indexes for the given particle.
Definition: tools.py:512