NetworkScience_igraphcluster
Joy Buongiorno
May 3rd, 2020
This is part 2 of our introduction to igraph. Last time, we introduced graph vizualization by using the American gut microbiome dataset included with the SpiecEasi package SpiecEasi documentation. Today, we will use those data again to recreate the networks a second time using the code from our last lesson, so a lot of the first section will be review. This time, however, we will build on this network a little. We will begin by calculating some graph statistics and then cluster vertices using a few options that are available. Finally, we will see how those clusters compare with what we may already know about our data (taxonomy).
Setup (Making the Network)
First we need a network to visualize.
You should already have installed SpiecEasi and igraph in previous lessons:
library(SpiecEasi)
library(igraph)##
## Attaching package: 'igraph'## The following object is masked from 'package:SpiecEasi':
##
## make_graph## The following objects are masked from 'package:stats':
##
## decompose, spectrum## The following object is masked from 'package:base':
##
## union#install.packages("viridis")
library(viridis)## Loading required package: viridisLitelibrary(ggplot2)
#install.packages("dplyr")
library(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:igraph':
##
## as_data_frame, groups, union## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionWe build a network using SpiecEasi.
Note: Later on, our names get lost and we are back to numbers. Not important for the sake of this tutorial.
## Generate network, see https://github.com/zdk123/SpiecEasi for a nice walkthrough of this process with this example
#Load data
data(amgut1.filt)
#Build network w/ spieceasi - It's gonna take a few minutes (sorry!)
se.gl.amgut <- spiec.easi(amgut1.filt, method='glasso', lambda.min.ratio=1e-2,
nlambda=20, pulsar.params=list(rep.num=50))## Applying data transformations...## Selecting model with pulsar using stars...## Fitting final estimate with glasso...## doneWe use the getRefit() function to extract the adjacency matrix from the spieceasi output. This is a square matrix with species on rows and columns and a one if two species are connected (“adjacent”) in the network and a zero otherwise.
#We want weights!
se.cor <- cov2cor(as.matrix(getOptCov(se.gl.amgut)))
weighted.adj.mat <- se.cor*getRefit(se.gl.amgut)grph <- adj2igraph(weighted.adj.mat)You’ll remember what our circular graph looked like:
plot(grph,vertex.size=1,
vertex.label=NA,
edge.width=1,
layout=layout.circle(grph))
We will clean this up as wel did last lesson, however, we aren’t going to make it shine like last time. For the purposes of clusering, we’re going to keep it simple. We’ll change our layout to the Fruchterman-Reingold layout so we can see clusters better. hhttps://igraph.org/r/doc/layout_with_fr.html
#Remove edges with very low weight
weight_threshold <- 0.01
grph <- delete.edges(grph,which(abs(E(grph)$weight)<weight_threshold))#Remove negative edges
grph.pos <- delete.edges(grph,which(E(grph)$weight<0))
plot(grph.pos,
vertex.label=NA,
edge.color="black",
layout=layout_with_fr(grph.pos))
#Remove unconnected vertices
grph.pos <- delete.vertices(grph.pos,which(degree(grph.pos)<1))
plot(grph.pos,
vertex.label=NA,
edge.color="black",
layout=layout_with_fr(grph.pos))
Now, let’s calculate some statistics about our graph. These statistics are often used to identify “hubs” in networks or assess how “clumpy” the network is. Identification of keystones species and other biologically-interesting hypotheses can be drawn in part from network statistics.
Berry and Widder (https://doi.org/10.3389/fmicb.2014.00219) and references therein explain one approach to using properties of co-occurrence networks to draw inferences about microbial community roles and functioning.
Let’s obtain the graph’s degree distribution.
dd.grph.pos <- degree.distribution(grph.pos)
plot(0:(length(dd.grph.pos)-1), dd.grph.pos, type='b',
ylab="Frequency", xlab="Degree", main="Degree Distributions")
Degree is the number of edges that a node (vertex) has. We see that most nodes are connected to few other nodes when consideringly only positive interactions above our weight threshold.
grph.pos_deg<-degree(grph.pos, v=V(grph.pos), mode="all")
fine = 500 # this will adjust the resolving power.
#this gives you the colors you want for every point
graphCol = viridis(fine)[as.numeric(cut(grph.pos_deg,breaks = fine))]
# now plot
plot(grph.pos, vertex.color=graphCol,
edge.color="black",
vertex.label=NA,
layout=layout_with_fr(grph.pos))
We can vizually confirm the degree distribution by rendering node color on the basis of degree value. Our viridis palette runs from deep purple (low degree) to yellow (high degree). Our graph has many deep purple nodes, indicating that most of our nodes have low degree.
Betweeness, or betweeness centrality, is based on shortest paths. In a network, to travel from one node to another, you must pass through one or several intermediate nodes. For a weighted network such as ours, betweenness centrality is caluculated by summing the weights of edges that represent the shortest paths between nodes.
grph.pos_bw<-betweenness(grph.pos, directed=F)
#this gives you the colors you want for every point
graphCol = viridis(fine)[as.numeric(cut(grph.pos_bw,breaks = fine))]
# now plot
plot(grph.pos, vertex.color=graphCol,
vertex.label=NA,
edge.color="black",
layout=layout_with_fr(grph.pos))
Our graph now has the node color rendered to reflect betweeneess centrality. In our graph, nodes with few or no paths running through them are deep purple, while nodes at the centers of denser areas with many paths running through them are bluer or yellow.
Another useful statistic is called transistivity, or the clustering coefficient. This is a measure of the probability that adjacent nodes of a particular node are connected to each other. In other words, it is the measure of how nodes tend to cluster together.
grph.pos_tran<-transitivity(grph.pos, type="local")
grph.pos_tran## [1] 0.6666667 1.0000000 0.0000000 0.8888889 1.0000000 0.7500000 1.0000000
## [8] NaN 1.0000000 0.6666667 0.8666667 0.4000000 1.0000000 0.6666667
## [15] 0.3333333 NaN 1.0000000 NaN 0.1666667 0.8571429 0.4000000
## [22] 0.6666667 1.0000000 0.0000000 0.9333333 NaN 0.9333333 1.0000000
## [29] 0.5000000 0.5909091 NaN NaN 0.0000000 NaN NaN
## [36] 1.0000000 0.6666667 0.8571429 NaN 0.1000000 NaN 0.9333333
## [43] NaN NaN 0.6666667 NaN 0.4000000 NaN 1.0000000
## [50] 0.6666667 0.9333333 NaN NaN NaN 0.5824176 1.0000000
## [57] 1.0000000 NaN 0.6666667 0.8888889 NaN 0.0000000 0.0000000
## [64] 1.0000000 0.8727273 0.7878788 0.3333333 NaN NaN 1.0000000
## [71] 0.0000000 0.4000000 0.9333333 NaN 0.6483516 1.0000000 NaN
## [78] 1.0000000 NaN 0.3333333 1.0000000 NaN#this gives you the colors you want for every point
graphCol = viridis(fine)[as.numeric(cut(grph.pos_tran,breaks = fine))]
# now plot
plot(grph.pos, vertex.color=graphCol,
vertex.label=NA,
edge.color="black",
layout=layout_with_fr(grph.pos))
Local transitivity cannot be calculated for many of our nodes (producing NaN in our results). This is because they do not meet the criteron of having two adjacent nodes. We can also calculate this statistic for the entire graph:
grph.pos_tran_gl<-transitivity(grph.pos, type="global")
grph.pos_tran_gl## [1] 0.7376033The value we get for transitivity for the entire graph is 0.73. This is the ratio of the triangles and the connected triples in the graph.
Another graph-wide statistic that is a nice segue into modularity analysis is assortativity. The assortativity coefficient measures the level of homophyly of the graph, or how likely things of the same origin or type are to cluster together. This can be based on some vertex labeling or values assigned to vertices. We will calculate categorical (or nominal assortivity) based on labeling using made-up family taxonomic assigments for each species.
The assortativity coefficient is positive if similar vertices (based on some external property) tend to connect to each other, and negative otherwise. If the coefficient is high, that means that connected vertices tend to have the same labels or similar assigned values. (https://igraph.org/r/doc/assortativity.html)
fams<-read.csv("gut_families.csv")
assortativity.nominal(grph.pos, as.integer(fams$Family), directed=F)## [1] -0.01096805vertex_attr(grph.pos, index = V(grph.pos)) In “gut_families.csv”, I made the names column by examining the output of ’vertex_attr(grph.pos, index = V(grph.pos)) to see which rows survived our culling process.
We have a negative value, meaning that our fake taxonomy doesn’t play a role in node clusetring.
We can begin to identify modules/clusters/cliques in our graph by testing different clustering functions.
1st up: Greedy. This function implements the fast greedy modularity optimization algorithm for finding community structure, see A Clauset, MEJ Newman, C Moore: Finding community structure in very large networks, http://www.arxiv.org/abs/cond-mat/0408187 for the details.
grph.pos.greedy <- cluster_fast_greedy(grph.pos, weights=E(grph.pos)$weight)
modularity(grph.pos.greedy)## [1] 0.568487sizes(grph.pos.greedy)## Community sizes
## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19
## 7 4 13 2 7 2 2 2 3 2 2 2 3 5 4 2 3 7 10Modularity tells us how good the division is, or how separatedthe different vertex types are from each other. Our value for “greedy” clustering is 0.5.
With this clustering method, there are 19 clusters found, each containing 2 through 13 nodes. Let’s plot this information on our graph.
colourCount = length(unique(grph.pos.greedy$membership)) # this will adjust the resolving power.
cluster_col = rainbow(colourCount)[as.numeric(cut(grph.pos.greedy$membership,breaks = colourCount))]
# now plot
plot(grph.pos, vertex.color=cluster_col,
vertex.label=NA,
edge.color="black",
layout=layout_with_fr(grph.pos))
The nodes in the graph are rendered by cluser membership.
Let’s try another clusering method, called Louvain. This is an unsupervised two-step method for community detection. There’s a nice explanation here: https://towardsdatascience.com/louvain-algorithm-93fde589f58c
grph.pos.louvain <- cluster_louvain(grph.pos, weights=E(grph.pos)$weight)
modularity(grph.pos.louvain)## [1] 0.5868466sizes(grph.pos.louvain)## Community sizes
## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19
## 7 5 3 3 7 7 2 2 2 2 10 2 4 2 3 13 2 4 2The modularity score is about the same. We have the same number of clusters as before.
colourCount = length(unique(grph.pos.louvain$membership)) # this will adjust the resolving power.
cluster_col = rainbow(colourCount)[as.numeric(cut(grph.pos.louvain$membership,breaks = colourCount))]
# now plot
plot(grph.pos, vertex.color=cluster_col,
vertex.label=NA,
edge.color="black",
layout=layout_with_fr(grph.pos))
The results are pretty intuitive. Even without sophistocated algorithms, we could see there are several unconnected components of the network after we remove negative associations and low weight edges. What if we were to use the pre-culled, positive matrix?
grph_whole <- adj2igraph(weighted.adj.mat)
grph_whole<-delete.edges(grph_whole,which(E(grph_whole)$weight<0))
grph.whole.louvain <- cluster_louvain(grph_whole, weights=E(grph_whole)$weight)
modularity(grph.whole.louvain)## [1] 0.5868951sizes(grph.whole.louvain)## Community sizes
## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
## 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 7 1 1 5 1 3 3 1 1
## 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50
## 1 7 1 7 2 1 2 1 2 2 1 1 1 1 11 1 2 4 1 1 1 1 1 1 1
## 51 52 53 54 55 56 57 58 59 60 61 62 63
## 1 2 1 3 1 13 1 1 1 1 2 4 2This gives us 63 clusters!! Geez
colourCount = length(unique(grph.whole.louvain$membership)) # this will adjust the resolving power.
cluster_col = rainbow(colourCount)[as.numeric(cut(grph.whole.louvain$membership,breaks = colourCount))]
# now plot
plot(grph_whole, vertex.color=cluster_col,
vertex.label=NA,
edge.color="black",
layout=layout_with_fr(grph_whole))
It looks like the clustering algorithm identified each lonley node as it’s own cluster. This is why it may be more informative to remove these isolated nodes (as we did in the beginning) if you are seeking to cluster/identify communities your networks.
Going back to our ‘grph.pos’ graph…. Community sizes analysis shows that many of our clusters contain only two nodes. Let’s ignore those.
V(grph.pos)$cluster=grph.pos.louvain$membership
vertex_attr(grph.pos, index = V(grph.pos))## $name
## [1] 1 2 4 5 11 12 13 14 15 16 19 23 24 25 27 28 29
## [18] 31 32 34 35 38 41 43 44 45 46 47 48 50 51 52 54 56
## [35] 57 58 59 60 61 62 63 65 66 67 68 69 71 72 74 75 76
## [52] 77 78 79 82 86 87 88 89 90 92 93 95 99 103 104 105 107
## [69] 109 110 112 114 116 117 119 120 122 123 124 125 126 127
##
## $cluster
## [1] 5 16 13 11 15 6 2 13 16 11 6 18 4 5 16 3 1 9 16 6 5 1 2
## [24] 3 11 2 11 4 1 6 16 7 16 12 17 5 5 6 8 16 14 6 5 7 16 10
## [47] 2 8 1 16 11 9 10 19 6 4 15 1 2 11 12 13 16 18 11 11 16 14 3
## [70] 15 16 1 11 13 11 16 17 18 19 5 18 1ids <- which(sizes(grph.pos.louvain)<=2)
grph.pos.main.communities <- delete_vertices(grph.pos,which(V(grph.pos)$cluster %in% ids))Let’s now see what types of bacteria make up our main clusters.
nodes <- V(grph.pos.main.communities)$name
nodes## [1] 1 2 4 5 11 12 13 14 15 16 19 23 24 25 27 28 29
## [18] 32 34 35 38 41 43 44 45 46 47 48 50 51 54 58 59 60
## [35] 62 65 66 68 71 74 75 76 82 86 87 88 89 90 93 95 99
## [52] 103 104 105 109 110 112 114 116 117 119 120 123 125 126 127#Our species names are numbers for the sake of this exercise. We have 66 species.
cluster_id <- V(grph.pos.main.communities)$cluster
nodes<-as.data.frame(cbind(nodes, cluster_id))
colnames(nodes)<-c("Species","Louvain Cluster")
nodes<-left_join(nodes, fams, by="Species")
nodes## Species Louvain Cluster Family
## 1 1 5 Enterobacteriaceae
## 2 2 16 Lactobacillaceae
## 3 4 13 Streptococcaceae
## 4 5 11 Staphylococcaceae
## 5 11 15 Enterococcaceae
## 6 12 6 Lachnospiraceae
## 7 13 2 Eryipelotrichaceae
## 8 14 13 Ruminococcaceae
## 9 15 16 Other
## 10 16 11 Enterobacteriaceae
## 11 19 6 Lactobacillaceae
## 12 23 18 Streptococcaceae
## 13 24 4 Staphylococcaceae
## 14 25 5 Enterococcaceae
## 15 27 16 Lachnospiraceae
## 16 28 3 Eryipelotrichaceae
## 17 29 1 Ruminococcaceae
## 18 32 16 Enterobacteriaceae
## 19 34 6 Lactobacillaceae
## 20 35 5 Streptococcaceae
## 21 38 1 Staphylococcaceae
## 22 41 2 Enterococcaceae
## 23 43 3 Lachnospiraceae
## 24 44 11 Eryipelotrichaceae
## 25 45 2 Ruminococcaceae
## 26 46 11 Other
## 27 47 4 Enterobacteriaceae
## 28 48 1 Lactobacillaceae
## 29 50 6 Streptococcaceae
## 30 51 16 Staphylococcaceae
## 31 54 16 Lachnospiraceae
## 32 58 5 Other
## 33 59 5 Enterobacteriaceae
## 34 60 6 Lactobacillaceae
## 35 62 16 Staphylococcaceae
## 36 65 6 Lachnospiraceae
## 37 66 5 Eryipelotrichaceae
## 38 68 16 Other
## 39 71 2 Lactobacillaceae
## 40 74 1 Staphylococcaceae
## 41 75 16 Enterococcaceae
## 42 76 11 Lachnospiraceae
## 43 82 6 Enterobacteriaceae
## 44 86 4 Lactobacillaceae
## 45 87 15 Streptococcaceae
## 46 88 1 Staphylococcaceae
## 47 89 2 Enterococcaceae
## 48 90 11 Lachnospiraceae
## 49 93 13 Ruminococcaceae
## 50 95 16 Other
## 51 99 18 Enterobacteriaceae
## 52 103 11 Lactobacillaceae
## 53 104 11 Streptococcaceae
## 54 105 16 Staphylococcaceae
## 55 109 3 Lachnospiraceae
## 56 110 15 Eryipelotrichaceae
## 57 112 16 Ruminococcaceae
## 58 114 1 Other
## 59 116 11 Enterobacteriaceae
## 60 117 13 Lactobacillaceae
## 61 119 11 Streptococcaceae
## 62 120 16 Staphylococcaceae
## 63 123 18 Lachnospiraceae
## 64 125 5 Ruminococcaceae
## 65 126 18 Other
## 66 127 1 EnterobacteriaceaeFamily_breakdown<-table(nodes$Family,nodes$`Louvain Cluster`)
Family_breakdown<-as.data.frame(Family_breakdown)
ggplot(Family_breakdown) +
geom_bar(aes(x = Var2, y = Freq, fill = Var1), stat = 'identity', width = 0.5) +
labs(x = "Louvain cluster",
y = "Count") +
guides(fill=guide_legend(title="Louvain cluster")) +
scale_fill_manual(values = c("#999999", "#E69F00", "#56B4E9", "#009E73",
"#F0E442", "#0072B2", "#D55E00", "#CC79A7", "#661100", "#44AA99")) +
theme_bw() +
theme(panel.grid.major = element_line(size = 0.2),
panel.grid.minor = element_line(size = 0),
axis.text = element_text(size=10),
axis.title = element_text(size=12),
axis.text.x = element_text(angle=90, hjust=1))
We see that cluster 3 is exclusively composed of two families (Eryipelotrichaceae and Lachnospiraceae), while other clusters, such as 5, are composed of 6 distinct families.
This has been fun! There are many more clustering functions avaiable in igraph. I suggested playing around with your options and taking a look at how each changes your modularity estimate. Than’s all for now!