1 from __future__
import print_function
9 from collections
import defaultdict
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.
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
38 self.selections = selections
39 self.contact_counts = {}
40 self.edges = defaultdict(int)
41 for (name1, name2)
in itertools.combinations(self.selections.keys(),
43 self.edges[tuple(sorted((name1, name2)))] = 0
45 self.frequency_cutoff = frequency_cutoff
47 self.gcpf.set_distance(self.cutoff)
48 self.names = list(self.selections.keys())
52 self.proteomic_edges = proteomic_edges
53 self.quantitative_proteomic_data = quantitative_proteomic_data
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)
63 ps_per_component = defaultdict(list)
65 self.size_per_component = defaultdict(int)
70 all_particles_by_resolution = []
71 for name
in part_dict:
72 all_particles_by_resolution += part_dict[name]
74 for component_name
in self.selections:
75 for seg
in self.selections[component_name]:
76 if isinstance(seg, str):
78 elif isinstance(seg, tuple):
80 hier, molecule=seg[2],
81 residue_indexes=range(seg[0], seg[1] + 1))
83 raise Exception(
'could not understand selection tuple '
85 parts = list(set(s.get_selected_particles())
86 & set(all_particles_by_resolution))
89 self.size_per_component[component_name] += \
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]))
99 self.edges[tuple(sorted((name1, name2)))] += 1.0
103 def make_plot(self, groups, out_fn, quantitative_proteomic_data=False):
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
112 import matplotlib.pyplot
as plt
115 ax.set_aspect(
'equal',
'box')
116 ax.xaxis.set_major_locator(plt.NullLocator())
117 ax.yaxis.set_major_locator(plt.NullLocator())
123 domain_xlocations = {}
124 domain_ylocations = {}
126 xoffset = squaredistance
127 yoffset = squaredistance
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
140 for edge, count
in self.edges.items():
142 if quantitative_proteomic_data:
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],
153 print(minqpd, maxqpd)
154 density = (1.0-(value-minqpd)/(maxqpd-minqpd))
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]]
165 rect = plt.Rectangle([x - squaresize / 2, y - squaresize / 2],
166 squaresize, squaresize,
167 facecolor=color, edgecolor=
'Gray',
170 rect = plt.Rectangle([y - squaresize / 2, x - squaresize / 2],
171 squaresize, squaresize,
172 facecolor=color, edgecolor=
'Gray',
181 def make_graph(self, out_fn):
184 print(
'num edges', len(self.edges))
185 for edge, count
in self.edges.items():
187 if float(count)/self.num_rmf > self.frequency_cutoff:
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,
202 validation_edges=self.proteomic_edges)
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
211 from matplotlib
import pyplot
212 from operator
import itemgetter
216 movers_mols_res = IMP.pmi.tools.OrderedDict()
217 for mv
in DegreesOfFreedom.movers_particles_map:
218 hs = DegreesOfFreedom.movers_particles_map[mv]
229 while not is_molecule:
234 if mv
not in movers_mols_res:
235 movers_mols_res[mv] = IMP.pmi.tools.OrderedDict()
238 if name
not in movers_mols_res[mv]:
241 movers_mols_res[mv][name].add(res)
246 for mv
in movers_mols_res:
252 srb_movers.append(mv)
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:
260 movers_mols_res[mv_rb][mol].remove(i)
261 except (ValueError, KeyError):
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]
276 for mv
in movers_mols_res:
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:
284 for seg
in movers_mols_res[mv][mol].segs:
286 elements[mol].append((seg[0], seg[-1],
" ",
"pdb",
288 except (KeyError, IndexError):
291 for mol
in movers_mols_res[mv]:
292 if mol
not in elements:
294 for seg
in movers_mols_res[mv][mol].segs:
295 elements[mol].append((seg[0], seg[-1],
" ",
"bead"))
299 elements[mol].sort(key=
lambda tup: tup[0])
300 elements[mol].append((elements[mol][-1][1],
301 elements[mol][-1][1],
" ",
"end"))
303 for name
in elements:
305 tmplist = sorted(elements[name], key=itemgetter(0))
307 endres = tmplist[-1][1]
309 print(
"draw_component_composition: missing information for "
310 "component %s" % name)
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])
318 norm = mpl.colors.Normalize(vmin=5, vmax=10)
322 for n, l
in enumerate(tmplist):
324 if bounds[-1] != l[0]:
325 colors.append(
"white")
328 colors.append(
"#99CCFF")
330 colors.append(
"white")
332 colors.append(
"#33CCCC")
334 bounds.append(l[1] + 1)
339 colors.append(
"white")
341 colors.append(
"#33CCCC")
343 bounds.append(l[1] + 1)
345 if bounds[-1] - 1 == l[0]:
349 colors.append(
"white")
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')
358 norm = mpl.colors.BoundaryNorm(bounds, cmap.N)
359 cb2 = mpl.colorbar.ColorbarBase(ax, cmap=cmap,
365 spacing=
'proportional',
366 orientation=
'horizontal')
375 mid = 1.0 / endres * float(tl[0])
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',
382 connectionstyle=
"angle,angleA=0,angleB=90,rad=10"))
383 extra_artists.append(t)
386 title = ax.text(-0.005, 0.5, k, ha=
"right", va=
"center", rotation=90,
389 extra_artists.append(title)
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]
399 for n
in range(1, len(bounds)):
400 if bounds[n] == bounds[n - 1]:
402 mid = 1.0 / endres * float(bounds[n])
403 if (float(bounds[n]) - float(bounds[n - 1])) / max <= 0.01:
405 offset = offsets[nclashes % 3]
410 t = ax.annotate(bounds[n], xy=(mid, 0),
411 xycoords=
'axes fraction',
412 xytext=(mid, -0.5 + offset),
413 textcoords=
'axes fraction')
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)
422 cb2.add_lines(bounds, [
"black"] * len(bounds), [1] * len(bounds))
425 pyplot.savefig(k +
"structure.pdf")
static bool get_is_setup(const IMP::ParticleAdaptor &p)
A class to read RMF files and make a network contact map.
static bool get_is_setup(const IMP::ParticleAdaptor &p)
A decorator to associate a particle with a part of a protein/DNA/RNA.
atom::Hierarchies create_hierarchies(RMF::FileConstHandle fh, Model *m)
Modify the transformation of a rigid body.
def make_plot
plot the interaction matrix
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)
A decorator for keeping track of copies of a molecule.
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.
Tools for clustering and cluster analysis.
Find all nearby pairs by testing all pairs.
Classes for writing output files and processing them.
A decorator for a residue.
The general base class for IMP exceptions.
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.
Select hierarchy particles identified by the biological name.
Support for the RMF file format for storing hierarchical molecular data and markup.