IMP logo
IMP Reference Guide  2.13.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 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  def add_rmf(self,rmf_fn,nframe):
51  """Add selections from an RMF file"""
52  print('reading from RMF file',rmf_fn)
53  rh = RMF.open_rmf_file_read_only(rmf_fn)
54  prots = IMP.rmf.create_hierarchies(rh, self.model)
55  hier = prots[0]
56  IMP.rmf.load_frame(rh, RMF.FrameID(0))
57  ps_per_component=defaultdict(list)
58  if self.num_rmf==0:
59  self.size_per_component=defaultdict(int)
60  self.model.update()
61 
62  #gathers particles for all components
64  all_particles_by_resolution = []
65  for name in part_dict:
66  all_particles_by_resolution += part_dict[name]
67 
68  for component_name in self.selections:
69  for seg in self.selections[component_name]:
70  if type(seg) == str:
71  s = IMP.atom.Selection(hier,molecule=seg)
72  elif type(seg) == tuple:
73  s = IMP.atom.Selection(hier,molecule=seg[2],
74  residue_indexes=range(seg[0], seg[1] + 1))
75  else:
76  raise Exception('could not understand selection tuple '+str(seg))
77  parts = list(set(s.get_selected_particles()) & set(all_particles_by_resolution))
78  ps_per_component[component_name] += IMP.get_indexes(parts)
79  if self.num_rmf==0:
80  self.size_per_component[component_name] += sum(len(IMP.pmi.tools.get_residue_indexes(p)) for p in parts)
81 
82  for n1,name1 in enumerate(self.names):
83  for name2 in self.names[n1+1:]:
84  ncontacts = len(self.gcpf.get_close_pairs(self.model,
85  ps_per_component[name1],
86  ps_per_component[name2]))
87  if ncontacts>0:
88  self.edges[tuple(sorted((name1,name2)))]+=1.0
89 
90  self.num_rmf+=1
91 
92  def make_plot(self,groups,out_fn,quantitative_proteomic_data=False):
93  '''
94  plot the interaction matrix
95  @param groups is the list of groups of domains, eg,
96  [["protA_1-10","prot1A_11-100"],["protB"]....]
97  it will plot a space between different groups
98  @param out_fn name of the plot file
99  @param quantitative_proteomic_data plot the quantitative proteomic data
100  '''
101  import numpy as np
102  import matplotlib.pyplot as plt
103  from matplotlib import cm
104 
105  ax=plt.gca()
106  ax.set_aspect('equal', 'box')
107  ax.xaxis.set_major_locator(plt.NullLocator())
108  ax.yaxis.set_major_locator(plt.NullLocator())
109 
110  largespace=0.6
111  smallspace=0.5
112  squaredistance=1.0
113  squaresize=0.99
114  domain_xlocations={}
115  domain_ylocations={}
116 
117  xoffset=squaredistance
118  yoffset=squaredistance
119  xlabels=[]
120  ylabels=[]
121  for group in groups:
122  xoffset+=largespace
123  yoffset+=largespace
124  for subgroup in group:
125  xoffset+=smallspace
126  yoffset+=smallspace
127  for domain in subgroup:
128  domain_xlocations[domain]=xoffset
129  domain_ylocations[domain]=yoffset
130  #rect = plt.Rectangle([xoffset- squaresize / 2, yoffset - squaresize / 2], squaresize, squaresize,
131  # facecolor=(1,1,1), edgecolor=(0.1,0.1,0.1))
132 
133  #ax.add_patch(rect)
134  #ax.text(xoffset , yoffset ,domain,horizontalalignment='left',verticalalignment='center',rotation=-45.0)
135  xoffset+=squaredistance
136  yoffset+=squaredistance
137 
138  for edge,count in self.edges.items():
139 
140  if quantitative_proteomic_data:
141  #normalize
142  maxqpd=max(self.quantitative_proteomic_data.values())
143  minqpd=min(self.quantitative_proteomic_data.values())
144  if edge in self.quantitative_proteomic_data:
145  value=self.quantitative_proteomic_data[edge]
146  elif (edge[1],edge[0]) in self.quantitative_proteomic_data:
147  value=self.quantitative_proteomic_data[(edge[1],edge[0])]
148  else:
149  value=0.0
150  print(minqpd,maxqpd)
151  density=(1.0-(value-minqpd)/(maxqpd-minqpd))
152  else:
153  density=(1.0-float(count)/self.num_rmf)
154  color=(density,density,1.0)
155  x=domain_xlocations[edge[0]]
156  y=domain_ylocations[edge[1]]
157  if x>y: xtmp=y; ytmp=x; x=xtmp; y=ytmp
158  rect = plt.Rectangle([x - squaresize / 2, y - squaresize / 2], squaresize, squaresize,
159  facecolor=color, edgecolor='Gray', linewidth=0.1)
160  ax.add_patch(rect)
161  rect = plt.Rectangle([y - squaresize / 2, x - squaresize / 2], squaresize, squaresize,
162  facecolor=color, edgecolor='Gray', linewidth=0.1)
163  ax.add_patch(rect)
164 
165  ax.autoscale_view()
166  plt.savefig(out_fn)
167  plt.show()
168  exit()
169 
170  def make_graph(self,out_fn):
171  edges=[]
172  weights=[]
173  print('num edges',len(self.edges))
174  for edge,count in self.edges.items():
175  # filter if frequency of contacts is greater than frequency_cutoff
176  if float(count)/self.num_rmf>self.frequency_cutoff:
177  print(count,edge)
178  edges.append(edge)
179  weights.append(count)
180  for nw,w in enumerate(weights):
181  weights[nw]=float(weights[nw])/max(weights)
182  IMP.pmi.output.draw_graph(edges,#node_size=1000,
183  node_size=dict(self.size_per_component),
184  node_color=self.colors,
185  fixed=self.fixed,
186  pos=self.pos,
187  edge_thickness=1, #weights,
188  edge_alpha=0.3,
189  edge_color='gray',
190  out_filename=out_fn,
191  validation_edges=self.proteomic_edges)
192 
193 
194 
195 def draw_component_composition(DegreesOfFreedom, max=1000, draw_pdb_names=False,show=True):
196  """A function to plot the representation on the sequence
197  @param DegreesOfFreedom input a IMP.pmi.dof.DegreesOfFreedom instance"""
198  import matplotlib as mpl
199  mpl.use('Agg')
200  from matplotlib import pyplot
201  from operator import itemgetter
202  import IMP.pmi.tools
203 
204  #first build the movers dictionary
205  movers_mols_res=IMP.pmi.tools.OrderedDict()
206  for mv in DegreesOfFreedom.movers_particles_map:
207  hs=DegreesOfFreedom.movers_particles_map[mv]
208  res=[]
209 
210  for h in hs:
212  res=[IMP.atom.Residue(h).get_index()]
214  res=IMP.atom.Fragment(h).get_residue_indexes()
215  if res:
216  is_molecule=False
217  hp=h
218  while not is_molecule:
219  hp=hp.get_parent()
220  is_molecule=IMP.atom.Molecule.get_is_setup(hp)
221  name=IMP.atom.Molecule(hp).get_name()+"."+str(IMP.atom.Copy(hp).get_copy_index())
222  if not mv in movers_mols_res:
223  movers_mols_res[mv]=IMP.pmi.tools.OrderedDict()
224  movers_mols_res[mv][name]=IMP.pmi.tools.Segments(res)
225  else:
226  if not name in movers_mols_res[mv]:
227  movers_mols_res[mv][name]=IMP.pmi.tools.Segments(res)
228  else:
229  movers_mols_res[mv][name].add(res)
230  # then get all movers by type
231  fb_movers=[]
232  rb_movers=[]
233  srb_movers=[]
234  for mv in movers_mols_res:
235  mvtype=None
236  if type(mv) is IMP.core.RigidBodyMover:
237  rb_movers.append(mv)
238  if type(mv) is IMP.core.BallMover:
239  fb_movers.append(mv)
240  if type(mv) is IMP.pmi.TransformMover:
241  srb_movers.append(mv)
242 
243  # now remove residue assigned to BallMovers from RigidBodies
244  for mv in fb_movers:
245  for mol in movers_mols_res[mv]:
246  for i in movers_mols_res[mv][mol].get_flatten():
247  for mv_rb in rb_movers:
248  try:
249  movers_mols_res[mv_rb][mol].remove(i)
250  except:
251  continue
252 
253  elements={}
254  rgb_colors=[(240,163,255),(0,117,220),
255  (153,63,0),(76,0,92),(25,25,25),
256  (0,92,49),(43,206,72),(255,204,153),
257  (128,128,128),(148,255,181),(143,124,0),
258  (157,204,0),(194,0,136),(0,51,128),(255,164,5),
259  (255,168,187),(66,102,0),(255,0,16),(94,241,242),
260  (0,153,143),(224,255,102),(116,10,255),(153,0,0),
261  (255,255,128),(255,255,0),(255,80,5)]
262  colors=[(float(c0)/255,float(c1)/255,float(c2)/255) for (c0,c1,c2) in rgb_colors]
263  mvrb_color={}
264  for mv in movers_mols_res:
265  mvtype=None
266  if type(mv) is IMP.core.RigidBodyMover:
267  mvtype="RB"
268  if not mv in mvrb_color:
269  color=colors[len(mvrb_color.keys())]
270  mvrb_color[mv]=color
271  for mol in movers_mols_res[mv]:
272  if not mol in elements: elements[mol]=[]
273  for seg in movers_mols_res[mv][mol].segs:
274  try:
275  elements[mol].append((seg[0],seg[-1]," ","pdb",mvrb_color[mv]))
276  except:
277  continue
278  if type(mv) is IMP.core.BallMover:
279  mvtype="FB"
280  for mol in movers_mols_res[mv]:
281  if not mol in elements: elements[mol]=[]
282  for seg in movers_mols_res[mv][mol].segs:
283  elements[mol].append((seg[0],seg[-1]," ","bead"))
284  if type(mv) is IMP.pmi.TransformMover:
285  mvtype="SRB"
286 
287  # sort everything
288  for mol in elements:
289  elements[mol].sort(key=lambda tup: tup[0])
290  elements[mol].append((elements[mol][-1][1],elements[mol][-1][1], " ", "end"))
291 
292  for name in elements:
293  k = name
294  tmplist = sorted(elements[name], key=itemgetter(0))
295  try:
296  endres = tmplist[-1][1]
297  except:
298  print("draw_component_composition: missing information for component %s" % name)
299  return
300  fig = pyplot.figure(figsize=(26.0 * float(endres) / max + 2, 2))
301  ax = fig.add_axes([0.05, 0.475, 0.9, 0.15])
302 
303  # Set the colormap and norm to correspond to the data for which
304  # the colorbar will be used.
305  cmap = mpl.cm.cool
306  norm = mpl.colors.Normalize(vmin=5, vmax=10)
307  bounds = [1]
308  colors = []
309 
310  for n, l in enumerate(tmplist):
311  firstres = l[0]
312  lastres = l[1]
313  if l[3] != "end":
314  if bounds[-1] != l[0]:
315  colors.append("white")
316  bounds.append(l[0])
317  if l[3] == "pdb":
318  colors.append("#99CCFF")
319  if l[3] == "bead":
320  colors.append("white")
321  if l[3] == "helix":
322  colors.append("#33CCCC")
323  if l[3] != "end":
324  bounds.append(l[1] + 1)
325  else:
326  if l[3] == "pdb":
327  colors.append(l[4])
328  if l[3] == "bead":
329  colors.append("white")
330  if l[3] == "helix":
331  colors.append("#33CCCC")
332  if l[3] != "end":
333  bounds.append(l[1] + 1)
334  else:
335  if bounds[-1] - 1 == l[0]:
336  bounds.pop()
337  bounds.append(l[0])
338  else:
339  colors.append("white")
340  bounds.append(l[0])
341 
342  bounds.append(bounds[-1])
343  colors.append("white")
344  cmap = mpl.colors.ListedColormap(colors)
345  cmap.set_over('0.25')
346  cmap.set_under('0.75')
347 
348  norm = mpl.colors.BoundaryNorm(bounds, cmap.N)
349  cb2 = mpl.colorbar.ColorbarBase(ax, cmap=cmap,
350  norm=norm,
351  # to use 'extend', you must
352  # specify two extra boundaries:
353  boundaries=bounds,
354  ticks=bounds, # optional
355  spacing='proportional',
356  orientation='horizontal')
357 
358  extra_artists = []
359  npdb = 0
360 
361  if draw_pdb_names:
362  for l in tmplist:
363  if l[3] == "pdb":
364  npdb += 1
365  mid = 1.0 / endres * float(l[0])
366  # t =ax.text(mid, float(npdb-1)/2.0+1.5, l[2], ha="left", va="center", rotation=0,
367  # size=10)
368  # t=ax.annotate(l[0],2)
369  t = ax.annotate(
370  l[2], xy=(mid, 1), xycoords='axes fraction',
371  xytext=(mid + 0.025, float(npdb - 1) / 2.0 + 1.5), textcoords='axes fraction',
372  arrowprops=dict(arrowstyle="->",
373  connectionstyle="angle,angleA=0,angleB=90,rad=10"),
374  )
375  extra_artists.append(t)
376 
377  # set the title of the bar
378  title = ax.text(-0.005, 0.5, k, ha="right", va="center", rotation=90,
379  size=20)
380 
381  extra_artists.append(title)
382  # changing the xticks labels
383  labels = len(bounds) * [" "]
384  ax.set_xticklabels(labels)
385  mid = 1.0 / endres * float(bounds[0])
386  t = ax.annotate(bounds[0], xy=(mid, 0), xycoords='axes fraction',
387  xytext=(mid - 0.01, -0.5), textcoords='axes fraction',)
388  extra_artists.append(t)
389  offsets = [0, -0.5, -1.0]
390  nclashes = 0
391  for n in range(1, len(bounds)):
392  if bounds[n] == bounds[n - 1]:
393  continue
394  mid = 1.0 / endres * float(bounds[n])
395  if (float(bounds[n]) - float(bounds[n - 1])) / max <= 0.01:
396  nclashes += 1
397  offset = offsets[nclashes % 3]
398  else:
399  nclashes = 0
400  offset = offsets[0]
401  if offset > -0.75:
402  t = ax.annotate(
403  bounds[n], xy=(mid, 0), xycoords='axes fraction',
404  xytext=(mid, -0.5 + offset), textcoords='axes fraction')
405  else:
406  t = ax.annotate(
407  bounds[n], xy=(mid, 0), xycoords='axes fraction',
408  xytext=(mid, -0.5 + offset), textcoords='axes fraction', arrowprops=dict(arrowstyle="-"))
409  extra_artists.append(t)
410 
411  cb2.add_lines(bounds, ["black"] * len(bounds), [1] * len(bounds))
412  # cb2.set_label(k)
413 
414  pyplot.savefig(
415  k + "structure.pdf")
416  #transparent="True")
417  #dpi=150,
418  #bbox_extra_artists=(extra_artists),
419  #bbox_inches='tight')
420  if show:
421  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:156
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:665
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)
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:135
The general base class for IMP exceptions.
Definition: exception.h:49
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:66
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:553