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.pmi1
8 import IMP.pmi1.output
9 from collections import defaultdict
10 import numpy as np
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,quantitative_proteomic_data=None):
16  """Set up a new graphXL object
17  @param model The IMP model
18  @param selection_dict A dictionary containing component names.
19  Keys are labels
20  values are either moleculename or start,stop,moleculename
21  @param cutoff The distance cutoff
22  @param frequency_cutoff The frequency cutoff
23  @param colors A dictionary of colors (HEX code,values) for subunits (keywords)
24  @param fixed A list of subunits that are kept fixed
25  @param pos A dictionary with positions (tuple, values) of subunits (keywords)
26  @param proteomic_edges A list edges to represent proteomic data
27  @param quantitative_proteomic_data A dictionary of edges to represent
28  quantitative proteomic data such as PE Scores, or genetic interactions
29  """
30  import itertools\
31 
32  self.model = model
33  self.selections = selections
34  self.contact_counts={}
35  self.edges=defaultdict(int)
36  for (name1,name2) in itertools.combinations(self.selections.keys(),2):
37  self.edges[tuple(sorted((name1,name2)))]=0
38  self.cutoff = cutoff
39  self.frequency_cutoff=frequency_cutoff
40  self.gcpf = IMP.core.GridClosePairsFinder()
41  self.gcpf.set_distance(self.cutoff)
42  self.names = list(self.selections.keys())
43  self.colors=colors
44  self.fixed=fixed
45  self.pos=pos
46  self.proteomic_edges=proteomic_edges
47  self.quantitative_proteomic_data=quantitative_proteomic_data
48  self.num_rmf=0
49 
50  @property
51  @IMP.deprecated_method("3.0", "Model should be accessed with `.model`.")
52  def mdl(self):
53  return self.model
54 
55  def add_rmf(self,rmf_fn,nframe):
56  """Add selections from an RMF file"""
57  print('reading from RMF file',rmf_fn)
58  rh = RMF.open_rmf_file_read_only(rmf_fn)
59  prots = IMP.rmf.create_hierarchies(rh, self.model)
60  hier = prots[0]
61  IMP.rmf.load_frame(rh, RMF.FrameID(0))
62  ps_per_component=defaultdict(list)
63  if self.num_rmf==0:
64  self.size_per_component=defaultdict(int)
65  self.model.update()
66 
67  #gathers particles for all components
69  all_particles_by_resolution = []
70  for name in part_dict:
71  all_particles_by_resolution += part_dict[name]
72 
73  for component_name in self.selections:
74  for seg in self.selections[component_name]:
75  if type(seg) == str:
76  s = IMP.atom.Selection(hier,molecule=seg)
77  elif type(seg) == tuple:
78  s = IMP.atom.Selection(hier,molecule=seg[2],
79  residue_indexes=range(seg[0], seg[1] + 1))
80  else:
81  raise Exception('could not understand selection tuple '+str(seg))
82  parts = list(set(s.get_selected_particles()) & set(all_particles_by_resolution))
83  ps_per_component[component_name] += IMP.get_indexes(parts)
84  if self.num_rmf==0:
85  self.size_per_component[component_name] += sum(len(IMP.pmi1.tools.get_residue_indexes(p)) for p in parts)
86 
87  for n1,name1 in enumerate(self.names):
88  for name2 in self.names[n1+1:]:
89  ncontacts = len(self.gcpf.get_close_pairs(self.model,
90  ps_per_component[name1],
91  ps_per_component[name2]))
92  if ncontacts>0:
93  self.edges[tuple(sorted((name1,name2)))]+=1.0
94 
95  self.num_rmf+=1
96 
97  def make_plot(self,groups,out_fn,quantitative_proteomic_data=False):
98  '''
99  plot the interaction matrix
100  @param groups is the list of groups of domains, eg,
101  [["protA_1-10","prot1A_11-100"],["protB"]....]
102  it will plot a space between different groups
103  @param out_fn name of the plot file
104  @param quantitative_proteomic_data plot the quantitative proteomic data
105  '''
106  import numpy as np
107  import matplotlib.pyplot as plt
108  from matplotlib import cm
109 
110  ax=plt.gca()
111  ax.set_aspect('equal', 'box')
112  ax.xaxis.set_major_locator(plt.NullLocator())
113  ax.yaxis.set_major_locator(plt.NullLocator())
114 
115  largespace=0.6
116  smallspace=0.5
117  squaredistance=1.0
118  squaresize=0.99
119  domain_xlocations={}
120  domain_ylocations={}
121 
122  xoffset=squaredistance
123  yoffset=squaredistance
124  xlabels=[]
125  ylabels=[]
126  for group in groups:
127  xoffset+=largespace
128  yoffset+=largespace
129  for subgroup in group:
130  xoffset+=smallspace
131  yoffset+=smallspace
132  for domain in subgroup:
133  domain_xlocations[domain]=xoffset
134  domain_ylocations[domain]=yoffset
135  #rect = plt.Rectangle([xoffset- squaresize / 2, yoffset - squaresize / 2], squaresize, squaresize,
136  # facecolor=(1,1,1), edgecolor=(0.1,0.1,0.1))
137 
138  #ax.add_patch(rect)
139  #ax.text(xoffset , yoffset ,domain,horizontalalignment='left',verticalalignment='center',rotation=-45.0)
140  xoffset+=squaredistance
141  yoffset+=squaredistance
142 
143  for edge,count in self.edges.items():
144 
145  if quantitative_proteomic_data:
146  #normalize
147  maxqpd=max(self.quantitative_proteomic_data.values())
148  minqpd=min(self.quantitative_proteomic_data.values())
149  if edge in self.quantitative_proteomic_data:
150  value=self.quantitative_proteomic_data[edge]
151  elif (edge[1],edge[0]) in self.quantitative_proteomic_data:
152  value=self.quantitative_proteomic_data[(edge[1],edge[0])]
153  else:
154  value=0.0
155  print(minqpd,maxqpd)
156  density=(1.0-(value-minqpd)/(maxqpd-minqpd))
157  else:
158  density=(1.0-float(count)/self.num_rmf)
159  color=(density,density,1.0)
160  x=domain_xlocations[edge[0]]
161  y=domain_ylocations[edge[1]]
162  if x>y: xtmp=y; ytmp=x; x=xtmp; y=ytmp
163  rect = plt.Rectangle([x - squaresize / 2, y - squaresize / 2], squaresize, squaresize,
164  facecolor=color, edgecolor='Gray', linewidth=0.1)
165  ax.add_patch(rect)
166  rect = plt.Rectangle([y - squaresize / 2, x - squaresize / 2], squaresize, squaresize,
167  facecolor=color, edgecolor='Gray', linewidth=0.1)
168  ax.add_patch(rect)
169 
170  ax.autoscale_view()
171  plt.savefig(out_fn)
172  plt.show()
173  exit()
174 
175  def make_graph(self,out_fn):
176  edges=[]
177  weights=[]
178  print('num edges',len(self.edges))
179  for edge,count in self.edges.items():
180  # filter if frequency of contacts is greater than frequency_cutoff
181  if float(count)/self.num_rmf>self.frequency_cutoff:
182  print(count,edge)
183  edges.append(edge)
184  weights.append(count)
185  for nw,w in enumerate(weights):
186  weights[nw]=float(weights[nw])/max(weights)
187  IMP.pmi1.output.draw_graph(edges,#node_size=1000,
188  node_size=dict(self.size_per_component),
189  node_color=self.colors,
190  fixed=self.fixed,
191  pos=self.pos,
192  edge_thickness=1, #weights,
193  edge_alpha=0.3,
194  edge_color='gray',
195  out_filename=out_fn,
196  validation_edges=self.proteomic_edges)
197 
198 
199 
200 def draw_component_composition(DegreesOfFreedom, max=1000, draw_pdb_names=False,show=True):
201  """A function to plot the representation on the sequence
202  @param DegreesOfFreedom input a IMP.pmi1.dof.DegreesOfFreedom instance"""
203  import matplotlib as mpl
204  mpl.use('Agg')
205  from matplotlib import pyplot
206  from operator import itemgetter
207  import IMP.pmi1.tools
208 
209  #first build the movers dictionary
210  movers_mols_res=IMP.pmi1.tools.OrderedDict()
211  for mv in DegreesOfFreedom.movers_particles_map:
212  hs=DegreesOfFreedom.movers_particles_map[mv]
213  res=[]
214 
215  for h in hs:
217  res=[IMP.atom.Residue(h).get_index()]
219  res=IMP.atom.Fragment(h).get_residue_indexes()
220  if res:
221  is_molecule=False
222  hp=h
223  while not is_molecule:
224  hp=hp.get_parent()
225  is_molecule=IMP.atom.Molecule.get_is_setup(hp)
226  name=IMP.atom.Molecule(hp).get_name()+"."+str(IMP.atom.Copy(hp).get_copy_index())
227  if not mv in movers_mols_res:
228  movers_mols_res[mv]=IMP.pmi1.tools.OrderedDict()
229  movers_mols_res[mv][name]=IMP.pmi1.tools.Segments(res)
230  else:
231  if not name in movers_mols_res[mv]:
232  movers_mols_res[mv][name]=IMP.pmi1.tools.Segments(res)
233  else:
234  movers_mols_res[mv][name].add(res)
235  # then get all movers by type
236  fb_movers=[]
237  rb_movers=[]
238  srb_movers=[]
239  for mv in movers_mols_res:
240  mvtype=None
241  if type(mv) is IMP.core.RigidBodyMover:
242  rb_movers.append(mv)
243  if type(mv) is IMP.core.BallMover:
244  fb_movers.append(mv)
245  if type(mv) is IMP.pmi1.TransformMover:
246  srb_movers.append(mv)
247 
248  # now remove residue assigned to BallMovers from RigidBodies
249  for mv in fb_movers:
250  for mol in movers_mols_res[mv]:
251  for i in movers_mols_res[mv][mol].get_flatten():
252  for mv_rb in rb_movers:
253  try:
254  movers_mols_res[mv_rb][mol].remove(i)
255  except:
256  continue
257 
258  elements={}
259  c=IMP.pmi1.tools.Colors()
260  colors=c.get_list_distant_colors()
261  colors=["blue", "red", "green", "pink", "cyan", "purple",
262  "magenta", "orange", "grey", "brown", "gold",
263  "khaki", "olive drab", "deep blue sky"]
264  rgb_colors=[(240,163,255),(0,117,220),
265  (153,63,0),(76,0,92),(25,25,25),
266  (0,92,49),(43,206,72),(255,204,153),
267  (128,128,128),(148,255,181),(143,124,0),
268  (157,204,0),(194,0,136),(0,51,128),(255,164,5),
269  (255,168,187),(66,102,0),(255,0,16),(94,241,242),
270  (0,153,143),(224,255,102),(116,10,255),(153,0,0),
271  (255,255,128),(255,255,0),(255,80,5)]
272  colors=[(float(c0)/255,float(c1)/255,float(c2)/255) for (c0,c1,c2) in rgb_colors]
273  mvrb_color={}
274  for mv in movers_mols_res:
275  mvtype=None
276  if type(mv) is IMP.core.RigidBodyMover:
277  mvtype="RB"
278  if not mv 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 not mol in elements: elements[mol]=[]
283  for seg in movers_mols_res[mv][mol].segs:
284  try:
285  elements[mol].append((seg[0],seg[-1]," ","pdb",mvrb_color[mv]))
286  except:
287  continue
288  if type(mv) is IMP.core.BallMover:
289  mvtype="FB"
290  for mol in movers_mols_res[mv]:
291  if not mol in elements: elements[mol]=[]
292  for seg in movers_mols_res[mv][mol].segs:
293  elements[mol].append((seg[0],seg[-1]," ","bead"))
294  if type(mv) is IMP.pmi1.TransformMover:
295  mvtype="SRB"
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],elements[mol][-1][1], " ", "end"))
301 
302  for name in elements:
303  k = name
304  tmplist = sorted(elements[name], key=itemgetter(0))
305  try:
306  endres = tmplist[-1][1]
307  except:
308  print("draw_component_composition: missing information for component %s" % name)
309  return
310  fig = pyplot.figure(figsize=(26.0 * float(endres) / max + 2, 2))
311  ax = fig.add_axes([0.05, 0.475, 0.9, 0.15])
312 
313  # Set the colormap and norm to correspond to the data for which
314  # the colorbar will be used.
315  cmap = mpl.cm.cool
316  norm = mpl.colors.Normalize(vmin=5, vmax=10)
317  bounds = [1]
318  colors = []
319 
320  for n, l in enumerate(tmplist):
321  firstres = l[0]
322  lastres = l[1]
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 l in tmplist:
373  if l[3] == "pdb":
374  npdb += 1
375  mid = 1.0 / endres * float(l[0])
376  # t =ax.text(mid, float(npdb-1)/2.0+1.5, l[2], ha="left", va="center", rotation=0,
377  # size=10)
378  # t=ax.annotate(l[0],2)
379  t = ax.annotate(
380  l[2], xy=(mid, 1), xycoords='axes fraction',
381  xytext=(mid + 0.025, float(npdb - 1) / 2.0 + 1.5), textcoords='axes fraction',
382  arrowprops=dict(arrowstyle="->",
383  connectionstyle="angle,angleA=0,angleB=90,rad=10"),
384  )
385  extra_artists.append(t)
386 
387  # set the title of the bar
388  title = ax.text(-0.005, 0.5, k, ha="right", va="center", rotation=90,
389  size=20)
390 
391  extra_artists.append(title)
392  # changing the xticks labels
393  labels = len(bounds) * [" "]
394  ax.set_xticklabels(labels)
395  mid = 1.0 / endres * float(bounds[0])
396  t = ax.annotate(bounds[0], xy=(mid, 0), xycoords='axes fraction',
397  xytext=(mid - 0.01, -0.5), textcoords='axes fraction',)
398  extra_artists.append(t)
399  offsets = [0, -0.5, -1.0]
400  nclashes = 0
401  for n in range(1, len(bounds)):
402  if bounds[n] == bounds[n - 1]:
403  continue
404  mid = 1.0 / endres * float(bounds[n])
405  if (float(bounds[n]) - float(bounds[n - 1])) / max <= 0.01:
406  nclashes += 1
407  offset = offsets[nclashes % 3]
408  else:
409  nclashes = 0
410  offset = offsets[0]
411  if offset > -0.75:
412  t = ax.annotate(
413  bounds[n], xy=(mid, 0), xycoords='axes fraction',
414  xytext=(mid, -0.5 + offset), textcoords='axes fraction')
415  else:
416  t = ax.annotate(
417  bounds[n], xy=(mid, 0), xycoords='axes fraction',
418  xytext=(mid, -0.5 + offset), textcoords='axes fraction', arrowprops=dict(arrowstyle="-"))
419  extra_artists.append(t)
420 
421  cb2.add_lines(bounds, ["black"] * len(bounds), [1] * len(bounds))
422  # cb2.set_label(k)
423 
424  pyplot.savefig(
425  k + "structure.pdf")
426  #transparent="True")
427  #dpi=150,
428  #bbox_extra_artists=(extra_artists),
429  #bbox_inches='tight')
430  if show:
431  pyplot.show()
static bool get_is_setup(const IMP::ParticleAdaptor &p)
Definition: Molecule.h:35
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)
A class to read RMF files and make a network contact map.
Legacy PMI1 module to represent, score, sample and analyze models.
Modify the transformation of a rigid body.
def get_particles_at_resolution_one
Get particles at res 1, or any beads, based on the name.
Move continuous particle variables by perturbing them within a ball.
def get_residue_indexes
Retrieve the residue indexes for the given particle.
Definition: /tools.py:1064
Tools for clustering and cluster analysis.
Definition: pmi1/Analysis.py:1
def make_plot
plot the interaction matrix
static bool get_is_setup(Model *m, ParticleIndex pi)
Definition: Fragment.h:46
Miscellaneous utilities.
Definition: /tools.py:1
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 deprecated_method
Python decorator to mark a method as deprecated.
Definition: __init__.py:9538
This class stores integers in ordered compact lists eg: [[1,2,3],[6,7,8]] the methods help splitting ...
Definition: /tools.py:1207
Modify the transformation of a rigid body.
void load_frame(RMF::FileConstHandle file, RMF::FrameID frame)
Load the given RMF frame into the state of the linked objects.
def add_rmf
Add selections from an RMF file.
Find all nearby pairs by testing all pairs.
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.
Classes for writing output files and processing them.
Definition: /output.py:1
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
def __init__
Set up a new graphXL object.
Support for the RMF file format for storing hierarchical molecular data and markup.