Networkx basics, some code and materials adapted from Hiroki Sayama's awesome textbook. http://textbooks.opensuny.org/introduction-to-the-modeling-and-analysis-of-complex-systems/
from __future__ import division
import networkx as nx
import matplotlib.pyplot as plt
import numpy as np
%matplotlib inline
plt.rcParams['figure.figsize'] = (20.0, 16.0)
Lots of stuff imported there, let's review: networkx is the tool we are going to use to analyze networks, matplolib.pyplot is for graphing stuff and graphical handling. numpy is math, matplotlib inline lets data show up in the notebook, and the next command makes our figures bigger so we can see them more easily.
First we need to talk about what a 'network' is. A network is any list of relationships between objects. It could be synaptic connections, friendships, facebook likes, interactions, roads, casual public intereactions. It can be directed (ie we know which object acts on the other), undirected (ie objects interact mutually), quantitative (we know how strongly objects interact in a quantitative sense), or binary (we jsut know that objects interact).
Lets start by making a hypothetical network g, and adding some nodes to it. A node is any object that interacts (or doesn't interact) in the network.
g=nx.Graph() #graph object
g.add_nodes_from (['hippocampus','piriform','entorhinal','perirhinal','frontal cortex'])
Now that we made some nodes, let's add some connections between those nodes, these will be called 'edges'
g.add_edges_from([('hippocampus','entorhinal'),('frontal cortex','hippocampus'),('piriform','entorhinal')])
Now that we've added some objects to our network let's print those details. Networkx gives us a few ways to do that.
print(g.nodes())
print(g.edges())
Today we won't be working with it, but in theory nodes and edges can have any number of traits. This is useful if you want to store lots of inforaation about the interaction, this comes in handy for filtering for analysis or visualization. Below I print all traits associated with each node, you'll notice it's empty.
g.node
Now we'll draw the network I made, don't worry about the syntax quite yet, but you'll notice it's pretty uninteresting
nx.draw(g,with_labels=True)
There is A LOT you can do with networks, today we'll focus on some basic visualization and quantiifcation of traits. In ther interest of space I'm going to grab one of the prebuilt graphs they have in NetworkX. In this case we'll be using the 'davis southern woman' network. They represent observed attendance at 14 social events by 18 Southern women. I'll declare this graph as G
G=nx.davis_southern_women_graph()
Let's draw it!
nx.draw(G)
print'basic'
There are lots of ways to plot this data, the default is called a spring network, imagine the nodes are springs, lots of nodes with other nodes in their neighborhood are more closely pulled together, nodes with few connections dangle far away. Next we'll draw nodes with random positions.
nx.draw_random(G)
print 'random positions'
We can also place all the nodes on the outside and see their connections, this is called a circular graph.
nx.draw_circular(G)
print 'circular graph'
The next class is called spectral format, it's complicated what it means
nx.draw_spectral(G)
print 'spectral graph'
net we are going to add the name of the nodes to the nodes, we do this with the tag with_labels=True, we are also defining the positions explicitly with the spring_layout command.
positions = nx.spring_layout(G)
nx.draw(G,positions, with_labels=True)
Next we will have each node be equivelent in size to the number connections it has, we'll discuss the degree command i ndepth soon, but for now, know that it means 'number of nodes that connect with this one'
node_sizes=[G.degree(i)*50 for i in G.nodes()]
nx.draw(G,positions,node_size=node_sizes)
That's a fun toy network but let's play with a REAL network, I've given you the allen brain atlas matrix, it's a list of all regions and their connective strengths in a mouse brain let's import it and get it plotted. We are going to import it with np.genfromtext, and we'll import it as a string so we can get the names associated with each row.
x = np.genfromtxt('ad_matrix.csv', delimiter=',', dtype='string'
)
Now let's grab the bit with numbers and convert it to a matrix! it's just that easy, this is called an 'adjacency matrix' where columns and rows meet at connection points. We'll also grab a variable called names for later
values=x[1:,1:].astype('float')
new_g=nx.from_numpy_matrix(values)
names=x[0][1:]
idx_name_dict=dict(zip(new_g.nodes(),names))
let's plot it! with the same size trick as before
node_sizes=[new_g.degree(i)**1.2 for i in new_g.nodes()]
positions = nx.spring_layout(new_g)
nx.draw(new_g,positions,node_size=node_sizes)
Well that's a bit of a mess, in reality, netowrkx is actually not very good for this kind of plotting, gephi is better. though networkx can help us clean the graph up for gephi. for the moment we won't worry about that, we'll work here. Lets add names to our nodes, to do this we make a paired list: with a range the length of our nodes and the names file, we'll then make this a dictionary, don't get too bogged down with the specifics of this, jsut know that it works for this conversion. We'll then use relabel and make a new matrix cleanedup
renamer=dict(zip(range(len(new_g.node)),names))
cleanedup=nx.relabel_nodes(new_g,renamer)
Run .node to confirm it worked.
cleanedup.node
and replot it with names, let's try the shell layout this time!
node_sizes=[cleanedup.degree(i)**1.2 for i in cleanedup.nodes()]
positions = nx.shell_layout(cleanedup)
nx.draw(cleanedup,positions,node_size=node_sizes,with_labels=True)
Still a total mess, let's see why by checking number of nodes, number of edges, and density (% of total connections)
cleanedup.number_of_edges()
cleanedup.number_of_nodes()
nx.density(cleanedup)
And we find 60% of the possible connections occur, let's shave this down! Luckily this network is quantitative, so let's say we only care about the top 5% most powerful connections. First go back to our original numpy graph. find the 5% cut off value and replot! numpy has a nice command percentile. Then we'll filter by only connections stronger than that and redo all the stuff to make a new grarph: powerful connections. Then we'll plot it
values=x[1:,1:].astype('float')
p=np.percentile(values,95)
top_5_percent=(values>p)*1
top_5_g=nx.from_numpy_matrix(top_5_percent)
top_5_g_with_proper_names=nx.relabel_nodes(top_5_g,renamer)
#now plot
node_sizes=[top_5_g_with_proper_names.degree(i)**1.2 for i in top_5_g_with_proper_names.nodes()]
positions = nx.spring_layout(top_5_g_with_proper_names)
nx.draw(top_5_g_with_proper_names,positions,node_size=node_sizes,with_labels=True)
This is way more usable, for fun try other percentile values, here we'll try 1
values=x[1:,1:].astype('float')
p=np.percentile(values,99)
top_1_percent=(values>p)*1
top_1_g=nx.from_numpy_matrix(top_1_percent)
top_1_g_with_proper_names=nx.relabel_nodes(top_1_g,renamer)
#now plot
node_sizes=[top_1_g_with_proper_names.degree(i)**1.2 for i in top_1_g_with_proper_names.nodes()]
positions = nx.spring_layout(top_1_g_with_proper_names)
nx.draw(top_1_g_with_proper_names,positions,node_size=node_sizes,with_labels=True)
Now we are getting cool. Of course I'm cheating here and using this as an undirected graph, for the moment I'm going to cheat and keep using that. For further analysis I'm going to use the 5% graph it's nicely connected without being painfully over connected. Depenging on your personal mathematical sophistication, most of the commands I'm going to use next can be adpated to both quantitative and directed graphs. Next we are going to evaluate network 'topology'/
print('density')
print(nx.density(top_5_g_with_proper_names))
next we'll look at connected components, this is a realatively meaningless analysis in the brain as all regions should be connected, but lets see what is stongly connected
ccs=nx.connected_components(top_5_g_with_proper_names)
ccsl=[i for i in ccs]
print ccsl
The two regions disconnected are subzones of the striatum (for the SF) and a subzone of the cerebellum. It's not surprising. Let's run an iterative shaving analysis:
for perc in range(95,100):
values=x[1:,1:].astype('float')
p=np.percentile(values,perc)
shaved=(values>p)*1
shaved=nx.from_numpy_matrix(shaved)
shaved=nx.relabel_nodes(shaved,renamer)
ccs=nx.connected_components(shaved)
ccsl=[i for i in ccs]
print(np.float(nx.density(shaved)))
print(ccsl)
Now this is a crazy feature of the brain, we can cut 99% of the connections and get down to 2% connection density and still have only a few regions knocked off. This is a common feature of preferentially attached networks, these are networks where very few nodes draw in a lot of the connections. Next lets see the average length between any two paths, this shoudl be farmiliar to you if you've ever played 6 degrees of Kevin bacon, in essence it's the average 'bacon number' of any two nodes in our network. First I have to remove the disconnected nodes!
top_5_g_with_proper_names.remove_node('SF')
top_5_g_with_proper_names.remove_node('PYR')
nx.average_shortest_path_length(top_5_g_with_proper_names)
There are lots of other metrics covered in Sayama including eccentricity, diameter and radius, rather than reproducing them here, I'll just reccomend you check out page 379.
Two cool metrics are 'perphiery' which is points that are furthest from all others, and 'center' whcih are points closest to all others.
print(nx.periphery(top_5_g_with_proper_names))
print(nx.center(top_5_g_with_proper_names))
The peripheral structures are a bunch of striatal, cerebellar, and spinal nuclei (shocker) and the central is msotly thalamic, midbrain, and isocortex.
There are lots of measures you can use to capture this more purely, by looking at neighbors and neighbors of neighbors. This is all pretty subtle and complicated check out sayama and other referencs for more explanation, so I'm only gonna look at one here: pagerank. Pagerank is basically, how many of your neighbros are well connected (ie who is popular amongst the popular kids)
page_ranks=nx.pagerank(top_5_g_with_proper_names)
pr_values=[page_ranks[i] for i in page_ranks]
#the next two lines are super dirty, I don't fully understand them and you don't ahve to either,I just looked up 'sort dictionary
#by values in python and using what showed up in stack exchange, sometimes you gotta do that
import operator
sorted_x = sorted(page_ranks.items(), key=operator.itemgetter(1))
sorted_x
The answers are actually a bit surprising, networks are weird! reticular nucleus, caudate and periductal greay (probably an artifact, PERI is the first one I actually believe) remember these are thresholded and bicirectional, in a proper nonthresholded analysis this probably won't be true
Last metric for today, clustering: clustering is measuring neighborhoods, a high degree of clustering means a suepr connected neighborhood, low clustering means more spread out, basically it's measure of how connected your neighbors are to eachother.
print(nx.clustering(top_5_g_with_proper_names))
print(nx.average_clustering(top_5_g_with_proper_names))
Plotting degree distribution
plt.plot(np.histogram(top_5_g_with_proper_names.degree().values(),bins=20)[1][1:]
,np.histogram(top_5_g_with_proper_names.degree().values(),bins=20)[0])
There is lots of other things and if you are going to write a real paper, I'd reccomend looking up assorativity and modularity cause it's cool, but we don't really have the time to discuss it here. Now go and do this with Swanson's network!
node_colors=[page_ranks[i] for i in top_5_g_with_proper_names.nodes()]
clust=nx.clustering(top_5_g_with_proper_names)
pos=nx.spring_layout(top_5_g_with_proper_names,iterations=500)
nx.draw(top_5_g_with_proper_names,pos,node_color=node_colors,cmap=plt.cm.Blues)
