"""
Plotting networks
"""
import numpy as np
import networkx as nx
import matplotlib.pyplot as plt
from matplotlib.patches import (Circle, Rectangle, FancyArrowPatch,
Arc, Polygon, ArrowStyle)
# Edge colors
activ = plt.get_cmap('tab10')(2)
inhib = plt.get_cmap('tab10')(3)
[docs]
def build_pos(inter, method=None):
"""
Compute node layout given an interaction matrix.
"""
G = inter.shape[0]
# Define graph
graph = nx.Graph()
for i in range(G):
graph.add_node(i)
for j in range(0,G):
graph.add_node(j)
if np.abs(inter[i,j] != 0):
graph.add_edge(i, j)
# Compute graph layout
p = np.random.normal(size=(G,2))
if method is None:
p = nx.kamada_kawai_layout(graph, pos=p)
if method == 'graphviz':
p = nx.nx_agraph.graphviz_layout(graph, prog='neato')
# Return node positions
pos = np.zeros((G,2))
for k in range(0,G):
pos[k,:] = p[k]
scale = 1 if method is None else 2/(pos.max()-pos.min())
center = pos.mean(axis=0) if G > 0 else 0
return scale * (pos - center)
#### Plotting functions ####
def node(k, ax, pos, name, scale=1., color=None, fontsize=None, nodesize=1.):
if color is None:
color = 'gray'
if fontsize is None:
fontsize = 8*scale
x, y = pos[k,0] * scale, pos[k,1] * scale
r = nodesize * 0.1 * scale
# Node shape
circle1 = Circle((x, y), radius=r, fill=True, facecolor='white',
edgecolor='lightgray', lw=1*scale, zorder=3, clip_on=False)
ax.add_artist(circle1)
# Node label
ax.text(x, y - 0.007*scale, name, color=color, fontsize=fontsize,
zorder=4, horizontalalignment='center', verticalalignment='center')
def link(k1, k2, ax, pos, weight, bend=0., scale=1., nodesize=1., alpha=None):
# Node coordinates
x1, y1 = pos[k1,0]*scale, pos[k1,1]*scale
x2, y2 = pos[k2,0]*scale, pos[k2,1]*scale
shrink = 9 * scale * nodesize
if alpha is None: alpha = 1
# Case 1: activation
if weight > 0:
style = ArrowStyle('Simple', tail_width=1.1*scale,
head_width=3.5*scale, head_length=5*scale)
arrow = FancyArrowPatch((x1,y1), (x2,y2), arrowstyle=style,
shrinkA=shrink, shrinkB=shrink, fc=activ, lw=0, zorder=0,
connectionstyle='arc3,rad={}'.format(bend),
clip_on=False, alpha=alpha)
ax.add_artist(arrow)
# Case 2: inhibition
if weight < 0:
style = ArrowStyle('Simple', tail_width=1.1*scale,
head_width=0, head_length=1e-9*scale)
arrow = FancyArrowPatch((x1,y1), (x2,y2), arrowstyle=style,
shrinkA=shrink, shrinkB=shrink, fc=inhib, lw=0, zorder=0,
connectionstyle='arc3,rad={}'.format(bend),
clip_on=False, alpha=alpha)
ax.add_artist(arrow)
r = 0.125*nodesize
u = pos[k2] - pos[k1]
u0 = u/np.sqrt(np.sum(u**2))
x0 = (pos[k1,0] + r*u0[0])*scale
y0 = (pos[k1,1] + r*u0[1])*scale
dx = (u[0] - 2*r*u0[0])*scale
dy = (u[1] - 2*r*u0[1])*scale
h_width = 0.07*scale
h_height = 0.015*scale
x1, y1 = x0 + dx, y0 + dy
d = np.sqrt(dx**2 + dy**2)
u = np.array([dx, dy])/d
v = np.array([dy,-dx])/d
x = x1 + 0.2*h_height*u[0] + 0.5*h_width*v[0]
y = y1 + 0.2*h_height*u[1] + 0.5*h_width*v[1]
if v[0] > 0: angle = (180/np.pi)*np.arctan(v[1]/v[0]) - 180
elif v[0] < 0: angle = (180/np.pi)*np.arctan(v[1]/v[0])
else: angle = np.sign(u[0]) * 90
angle += np.tanh(bend) * 90
theta = np.tanh(bend) * (np.pi/2)
x3 = x2 + np.cos(theta)*(x-x2) - np.sin(theta)*(y-y2)
y3 = y2 + np.sin(theta)*(x-x2) + np.cos(theta)*(y-y2)
head = Rectangle((x3,y3), h_width, h_height, clip_on=False,
angle=angle, fc=inhib, zorder=0, alpha=alpha)
ax.add_artist(head)
def link_auto(k, ax, pos, weight, v=None, scale=1., nodesize=1., alpha=None):
# Node coordinates
x0, y0 = pos[k,0]*scale, pos[k,1]*scale
# Orientation
if v is None:
v = np.array([1,0])
v /= np.sqrt(np.sum(v**2))
if v[0] > 0: angle0 = np.arcsin(v[1]) * (180/np.pi)
else: angle0 = 180 - np.arcsin(v[1]) * (180/np.pi)
angle1 = 58 + angle0
angle2 = 300 + angle0
theta = angle0 * (np.pi/180)
x, y = x0 + 0.145*scale, y0
x1 = x0 + np.cos(theta)*(x-x0) - np.sin(theta)*(y-y0)
y1 = y0 + np.sin(theta)*(x-x0) + np.cos(theta)*(y-y0)
# Loop diameter
d = 0.15*scale
# Case 1: activation
if weight > 0:
c = activ
arc = Arc((x1,y1), d, d, angle=180, theta1=angle1, theta2=angle2-30,
lw=1.1*scale, color=c, zorder=0, clip_on=False)
ax.add_artist(arc)
# Activation head
head_width = 0.045*scale
head_length = 0.065*scale
angle = angle2 - 180 - 30
theta = angle * (np.pi/180)
triangle = np.array([
[x1 + d/2 + head_width/2, y1 - head_length/4],
[x1 + d/2, y1 + head_length * 3/4],
[x1 + d/2 - head_width/2, y1 - head_length/4]])
x, y = triangle[:,0].copy(), triangle[:,1].copy()
triangle[:,0] = x1 + np.cos(theta)*(x-x1) - np.sin(theta)*(y-y1)
triangle[:,1] = y1 + np.sin(theta)*(x-x1) + np.cos(theta)*(y-y1)
triangle -= 0.03 * (triangle[0]-triangle[2])
head = Polygon(triangle, fc=activ, zorder=0, clip_on=False)
ax.add_artist(head)
# Case 2: inhibition
if weight < 0:
c = inhib
arc = Arc((x1,y1), d, d, angle=180, theta1=angle1, theta2=angle2-2,
lw=1.1*scale, color=c, zorder=0, clip_on=False)
ax.add_artist(arc)
# Inhibition head
width = 0.06*scale
height = 0.015*scale
angle = angle2 - 180 - 2
theta = angle * (np.pi/180)
x = x1 - width/2 + d/2
y = y1 - height/2
x2 = x1 + np.cos(theta)*(x-x1) - np.sin(theta)*(y-y1)
y2 = y1 + np.sin(theta)*(x-x1) + np.cos(theta)*(y-y1)
head = Rectangle((x2,y2), width, height, angle=angle, fc=c,
zorder=0, clip_on=False)
ax.add_artist(head)
def show_empty_plot_warning():
print('Warning: nothing to show in this plot')
def is_isolated(gene, inter):
"""Test whether a gene is isolated."""
G = inter.shape[0]
others = (np.arange(G) != gene)
targets = np.nonzero(inter[gene,others])[0]
factors = np.nonzero(inter[others,gene])[0]
return (targets.size == 0) and (factors.size == 0)
def is_stimulus_leaf(gene, inter):
"""Test whether a gene is a leaf of the stimulus."""
G = inter.shape[0]
others = (np.arange(G) != gene)
targets = np.nonzero(inter[gene,others])[0]
factors = np.nonzero(inter[others,gene])[0]
if factors.size != 1:
return False
else:
return (targets.size == 0) and (factors[0] == 0)
#### Main function ####
[docs]
def plot_network(inter, pos, width=1., height=1., scale=1., names=None,
vdict=None, tol=None, root=False, axes=None, nodes=None, n0=True,
file=None, verb=False, fontsize=None, vcolor=None, nodesize=1.,
bend=0.14, bend_all=False, alpha=None,
hide_isolated_genes=False, hide_stimulus_leaves=False):
"""Plot a gene regulatory network.
Parameters
----------
hide_isolated_genes : bool
Hide genes that have no neighbors.
hide_stimulus_leaves : bool
Hide genes whose only neighbor is the stimulus.
"""
G, G = inter.shape
w, h = width/2.54, height/2.54 # Centimeters
if names is None: names = [''] + ['{}'.format(k) for k in range(1,G)]
if vcolor is None: vcolor = G * ['#5C5C5C']
if vdict is None: vdict = {}
if axes is None:
fig = plt.figure(figsize=(w,h), dpi=100, frameon=False)
plt.axes((0,0,1,1))
ax = fig.gca()
I, J = inter.nonzero()
ax.axis('off')
else:
ax = axes
fig = plt.gcf()
size = fig.get_size_inches()
w, h = size[0], size[1]
if nodes is None: I, J = inter.nonzero()
else: I, J = nodes, nodes
ax.axis('equal')
ax.axis('off')
ax_pos = ax.get_position()
plt.xlim([-w*ax_pos.width/2, w*ax_pos.width/2])
plt.ylim([-h*ax_pos.height/2, h*ax_pos.height/2])
scale = scale * np.min([ax_pos.width,ax_pos.height])
# Decide which genes to show
v = set(range(G))
e = set(zip(*inter.nonzero()))
if G == 0:
show_empty_plot_warning()
if hide_isolated_genes:
for k in range(1,G):
if is_isolated(k, inter):
v.discard(k)
e.discard((k,k))
if hide_stimulus_leaves:
for k in range(1,G):
if is_stimulus_leaf(k, inter):
v.discard(k)
e.discard((0,k))
# Draw nodes
for k in v:
node(k, ax, pos, names[k], scale, fontsize=fontsize, color=vcolor[k],
nodesize=nodesize)
# Draw links
for k1, k2 in e:
weight = inter[k1,k2]
if k1 != k2:
if (k2, k1) in e:
link(k1, k2, ax, pos, weight, bend, scale, nodesize=nodesize,
alpha=alpha)
else: link(k1, k2, ax, pos, weight, bend*bend_all, scale,
nodesize=nodesize, alpha=alpha)
else:
v = vdict.get(k1)
if v is None:
b, c = 0, 0
for k in range(G):
if (k != k1) and (((k,k1) in e) or ((k1, k) in e)):
b += pos[k]
c += 1
if c == 0:
b, c = 0, 0
for k in np.union1d(I, J):
if (k != k1):
b += pos[k]
c += 1
v = pos[k1] - b/c
else: v = np.array(v)
d = np.sqrt(np.sum(v**2))
if d > 0: v = v/d
else: v = np.array([0,1])
link_auto(k1, ax, pos, weight, v, scale, nodesize=nodesize,
alpha=alpha)
if file is None: file = 'network.pdf'
if axes is None: fig.savefig(file, bbox_inches='tight')