Disclaimer. I am not an expert on Igraph, I just transformed a lab workshop into this tutorial hoping it would be useful for other people outside the lab. Here I assume you know a little about ecological networks and (very little) about matrices. And obviously, some R.
First of all, load the igraph package (!):
library(igraph)
The most annoying thing about Igraph is that it only accepts Igraph objects. It is possible to create networks from scratch or to import them. Below we create a network that is undirected, with four nodes, each represented by a number. The arrows indicate interactions between nodes: 1->2, 2->3, 3->1, 4->2.
g1 <- graph(edges=c(1,2, 2,3, 3, 1, 4,2), n=4, directed=F)
class(g1) #checking the class
## [1] "igraph"
plot(g1) # an ugly plot to get started
There are several ways to create an Igraph object:
graph.full()
graph.star()
graph.ring()
graph.empty()
When you work with ecological networks it is often the case that you want to import interaction matrices, in which rows always represent species from one group (e.g., plants) and columns represent species from the other group (e.g., pollinators). We will use the network from Bezerra et al. (2009), available in the bipartite package:
library(bipartite)
pp<-bezerra2009 #plant-pollinator interaction matrix
Interaction networks can be either binary or quantitative (like bezerra2009), incidence or adjacency (square). Igraph has functions to deal with all of these, for example:
igraph::graph_from_adj_list()
igraph::graph_from_adjacency_matrix()
igraph::graph_from_data_frame()
igraph::graph_from_edgelist()
igraph::graph_from_incidence_matrix()
In our case, we have a quantitative incidence matrix, so:
ppIg<-graph_from_incidence_matrix(pp)
ppIg
## IGRAPH UN-B 26 71 --
## + attr: type (v/l), name (v/c)
## + edges (vertex names):
## [1] Diplopterys.pubipetala--Centris.aenea
## [2] Diplopterys.pubipetala--Centris.fuscata
## [3] Diplopterys.pubipetala--Centris.caxiensis
## [4] Diplopterys.pubipetala--Centris.tarsata
## [5] Diplopterys.pubipetala--Centris.flavifrons
## [6] Diplopterys.pubipetala--Centris.trigonoides
## [7] Diplopterys.pubipetala--Epicharis.sp2
## [8] Diplopterys.pubipetala--Apis.mellifera
## + ... omitted several edges
Igraph objects have several properties: the first line (IGRAPH UN-B 26 71) indicates the type of object we just created: the first slot (“U”) tell us if the network is directed (“D”) or undirected (“U”); the second whether the nodes (species) have names (“N”) or not; the third slot indicated by a “-” tell us weather the network is weighted or not; and finally the B=bipartite tell us what type of network (bipartite or unipartite), that is, if there are different types of nodes in the network. The subsequent numbers indicate the number of nodes/vertices/species (IGRAPH UN-B 26 71) and the total number of edges/interactions interactions (IGRAPH UN-B 26 71). BUT, we have an weighted network, so we have to tell that to Igraph when importing the interaction matrix.
ppIg<-graph_from_incidence_matrix(pp, weighted = TRUE)
ppIg
## IGRAPH UNWB 26 71 --
## + attr: type (v/l), name (v/c), weight (e/n)
## + edges (vertex names):
## [1] Diplopterys.pubipetala--Centris.aenea
## [2] Diplopterys.pubipetala--Centris.fuscata
## [3] Diplopterys.pubipetala--Centris.caxiensis
## [4] Diplopterys.pubipetala--Centris.tarsata
## [5] Diplopterys.pubipetala--Centris.flavifrons
## [6] Diplopterys.pubipetala--Centris.trigonoides
## [7] Diplopterys.pubipetala--Epicharis.sp2
## [8] Diplopterys.pubipetala--Apis.mellifera
## + ... omitted several edges
Now we are all set =)
Igraph allows us to visualize and assign attributes to all vertices (species) and edges (interactions):
V(ppIg) #checking species' names
## + 26/26 vertices, named:
## [1] Diplopterys.pubipetala Byrsonima.gardnerana
## [3] Banisteriopsis.muricata Heteropterys.sp1
## [5] Heteropterys.sp2 Dicella.bracteosa
## [7] Carolus.chasei Stigmaphyllon.paralias
## [9] Banisteriopsis.stellaris Banisteriopsis.schizoptera
## [11] Stigmaphyllon.auriculatum Stigmaphyllon.ciliatum
## [13] Janusia.anisandra Centris.aenea
## [15] Centris.fuscata Centris.caxiensis
## [17] Centris.tarsata Centris.flavifrons
## [19] Centris.trigonoides Centris.obsoleta
## + ... omitted several vertices
E(ppIg) #checking interactions
## + 71/71 edges (vertex names):
## [1] Diplopterys.pubipetala--Centris.aenea
## [2] Diplopterys.pubipetala--Centris.fuscata
## [3] Diplopterys.pubipetala--Centris.caxiensis
## [4] Diplopterys.pubipetala--Centris.tarsata
## [5] Diplopterys.pubipetala--Centris.flavifrons
## [6] Diplopterys.pubipetala--Centris.trigonoides
## [7] Diplopterys.pubipetala--Epicharis.sp2
## [8] Diplopterys.pubipetala--Apis.mellifera
## [9] Diplopterys.pubipetala--Centris.sp3
## [10] Byrsonima.gardnerana --Centris.aenea
## + ... omitted several edges
Another way to check species’ attributes is:
vertex_attr(ppIg) #checking species' attributes
## $type
## [1] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [12] FALSE FALSE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [23] TRUE TRUE TRUE TRUE
##
## $name
## [1] "Diplopterys.pubipetala" "Byrsonima.gardnerana"
## [3] "Banisteriopsis.muricata" "Heteropterys.sp1"
## [5] "Heteropterys.sp2" "Dicella.bracteosa"
## [7] "Carolus.chasei" "Stigmaphyllon.paralias"
## [9] "Banisteriopsis.stellaris" "Banisteriopsis.schizoptera"
## [11] "Stigmaphyllon.auriculatum" "Stigmaphyllon.ciliatum"
## [13] "Janusia.anisandra" "Centris.aenea"
## [15] "Centris.fuscata" "Centris.caxiensis"
## [17] "Centris.tarsata" "Centris.flavifrons"
## [19] "Centris.trigonoides" "Centris.obsoleta"
## [21] "Epicharis.sp2" "Apis.mellifera"
## [23] "Centris.sp3" "Centris.sp1"
## [25] "Xylocopa.sp" "Xylocopa.grisescens"
We can see that in addition to the species’ names Igraph also imported a “$type” attribute, which happens in bipartite graphs, when lines are “one type” and columns are “another type” of node. In our case, the vertices indicated by “FALSE” are the plants, and the vertices indicated “TRUE” are the animals. Likewise, we can check the interactions’ attributes:
edge_attr(ppIg) #checking interactions' attributes
## $weight
## [1] 1368 1364 740 460 416 256 328 364 368 924 320 2108 464 284
## [15] 28 396 468 108 140 272 652 912 364 44 368 164 84 764
## [29] 680 528 308 404 300 28 76 740 656 528 332 324 116 116
## [43] 556 512 356 132 524 604 452 432 200 504 816 292 300 244
## [57] 116 228 224 124 120 240 164 68 196 268 196 164 188 244
## [71] 96
We can see that there is a list of $weights, which are the interactions weights present in the original quantitative matrix. We can add any type of attribute. For example, let us color the plants as green vertices and the animals as orange vertices:
vertex_attr(ppIg)$color<-rep("#FF9C0B", length(V(ppIg))) #coloring all the nodes orange
#Selecting everyone that is "TRUE" in $type to change it to green (take a look at the help of the grep function in case you don't know it).
vertex_attr(ppIg)$color[grep(pattern = "FALSE", vertex_attr(ppIg)$type)]<-"#1C912C"
vertex_attr(ppIg)$color # checking the color vector
## [1] "#1C912C" "#1C912C" "#1C912C" "#1C912C" "#1C912C" "#1C912C" "#1C912C"
## [8] "#1C912C" "#1C912C" "#1C912C" "#1C912C" "#1C912C" "#1C912C" "#FF9C0B"
## [15] "#FF9C0B" "#FF9C0B" "#FF9C0B" "#FF9C0B" "#FF9C0B" "#FF9C0B" "#FF9C0B"
## [22] "#FF9C0B" "#FF9C0B" "#FF9C0B" "#FF9C0B" "#FF9C0B"
Likewise, you can add any other attributes you want, such as size or shape for example.
There are several ways to plot a networks using Igraph (see here a list of possible arguments). The default is very simple (and ugly!):
plot(ppIg)
Sometimes Igraph is kind of smart (!), and it already recognized the color vector we created. Another way is to include it as an argument:
plot(ppIg, vertex.color=vertex_attr(ppIg)$cor)
But it still ugly. How about taking out the names?
plot(ppIg, vertex.color=vertex_attr(ppIg)$cor,
vertex.label=NA)
A little better, We can also change the vertices sizes according to the species’ degree:
plot(ppIg, vertex.color=vertex_attr(ppIg)$cor,vertex.label=NA,
vertex.size=igraph::degree(ppIg)) #We have to use "igraph::degree" to indicate to R that degree should be calculated using the Igraph package.
It is too small now. Let’s increase it a little:
plot(ppIg, vertex.color=vertex_attr(ppIg)$cor,vertex.label=NA,
vertex.size=2*igraph::degree(ppIg))
We can also change the edges’ thickness:
plot(ppIg, vertex.color=vertex_attr(ppIg)$cor,vertex.label=NA,
vertex.size=2*igraph::degree(ppIg), edge.width=2)
Or put it proportional do the weight:
plot(ppIg, vertex.color=vertex_attr(ppIg)$cor,vertex.label=NA,
vertex.size=2*igraph::degree(ppIg),
edge.width=(edge_attr(ppIg)$weight))
Whoops! The weights are not proportional so we ruined it…
edge_attr(ppIg)$weight
## [1] 1368 1364 740 460 416 256 328 364 368 924 320 2108 464 284
## [15] 28 396 468 108 140 272 652 912 364 44 368 164 84 764
## [29] 680 528 308 404 300 28 76 740 656 528 332 324 116 116
## [43] 556 512 356 132 524 604 452 432 200 504 816 292 300 244
## [57] 116 228 224 124 120 240 164 68 196 268 196 164 188 244
## [71] 96
Adjusting:
plot(ppIg, vertex.color=vertex_attr(ppIg)$cor,vertex.label=NA,
vertex.size=2*igraph::degree(ppIg),
edge.width=3*(edge_attr(ppIg)$weight)/1000)
To make the figure a little better I like to curve the lines a little, and because our matrix has edges that are super weak, let us also darken the edges’ color a little:
plot(ppIg, vertex.color=vertex_attr(ppIg)$cor,vertex.label=NA,
vertex.size=2*igraph::degree(ppIg),
edge.width=3*(edge_attr(ppIg)$weight)/1000,
edge.color="grey50", #changing the edge color (!)
edge.curved=0.3) #curving the edges (!)
And finally, we can also change the graph layout according to some algorithms already implemented on Igraph, such as:
?layout_with_dh
?layout_with_fr
?layout_with_kk
?layout_with_sugiyama
To adjust the figure you just have to assign a new layout:
l<-layout_with_dh(ppIg)
plot(ppIg, vertex.color=vertex_attr(ppIg)$cor,vertex.label=NA,
vertex.size=2*igraph::degree(ppIg),
edge.width=3*(edge_attr(ppIg)$weight)/1000,
edge.color="grey50",
edge.curved=0.3,
layout=l)
Another cool layouts:
par(mfrow=c(2,2), oma=c(0,0,0,0), mar=c(1,1,1,1))#To plot four plots
plot(ppIg, vertex.color=vertex_attr(ppIg)$cor,vertex.label=NA,
vertex.size=2*igraph::degree(ppIg),
edge.width=3*(edge_attr(ppIg)$weight)/1000,
edge.color="grey50",
edge.curved=0.3,
layout=layout_in_circle, main="layout_in_circle")
plot(ppIg, vertex.color=vertex_attr(ppIg)$cor,vertex.label=NA,
vertex.size=2*igraph::degree(ppIg),
edge.width=3*(edge_attr(ppIg)$weight)/1000,
edge.color="grey50",
#edge.curved=0.3,
layout=layout_as_bipartite, main="layout_as_bipartite")
plot(ppIg, vertex.color=vertex_attr(ppIg)$cor,vertex.label=NA,
vertex.size=2*igraph::degree(ppIg),
edge.width=3*(edge_attr(ppIg)$weight)/1000,
edge.color="grey50",
#edge.curved=0.3,
layout=layout_as_tree, main="layout_as_tree")
plot(ppIg, vertex.color=vertex_attr(ppIg)$cor,vertex.label=NA,
vertex.size=2*igraph::degree(ppIg),
edge.width=3*(edge_attr(ppIg)$weight)/1000,
edge.color="grey50",
#edge.curved=0.3,
layout=layout_as_star, main="layout_as_star")
Another cool thing of Igraph is that there are several built in functions to analyse graph properties. For example, you can check for the presence of modules (or communities, as they are called outside of the ecological literature):
cluster_spinglass()
cluster_optimal()
cluster_walktrap()
Let us calculate modularity using the cluster_spinglass function, that uses simulated annealing for optimization:
modulos<-cluster_spinglass(ppIg)
str(modulos) #evaluating o output
## List of 6
## $ membership : chr [1:5] "Stigmaphyllon.paralias" "Banisteriopsis.stellaris" "Banisteriopsis.schizoptera" "Janusia.anisandra" ...
## $ csize : chr [1:4] "Heteropterys.sp2" "Carolus.chasei" "Stigmaphyllon.ciliatum" "Centris.aenea"
## $ modularity : chr [1:3] "Diplopterys.pubipetala" "Apis.mellifera" "Centris.sp3"
## $ temperature: chr [1:4] "Byrsonima.gardnerana" "Heteropterys.sp1" "Centris.caxiensis" "Centris.tarsata"
## $ algorithm : chr [1:3] "Dicella.bracteosa" "Stigmaphyllon.auriculatum" "Centris.flavifrons"
## $ vcount : chr [1:7] "Banisteriopsis.muricata" "Centris.trigonoides" "Centris.obsoleta" "Epicharis.sp2" ...
## - attr(*, "class")= chr "communities"
modulos$membership #checking to which module each species belongs to
## [1] 3 4 6 4 2 5 2 1 1 1 5 2 1 2 1 4 4 5 6 6 6 3 3 6 6 6
(Please note that the modules have both plants and pollinators; it is beyond the scope of this tutorial to talk about modularity and the different ways to calculate it).
We can then plot the result emphasizing the groups:
plot(ppIg, vertex.color=vertex_attr(ppIg)$cor,vertex.label=NA,
vertex.size=2*igraph::degree(ppIg),
edge.width=3*(edge_attr(ppIg)$weight)/1000,
edge.color="grey50",
edge.curved=0.3, #layout=l)
mark.groups= modulos, mark.border=NA)
Another cool tool available in Igraph is to calculate the shortest path between two species (Xylocopa.sp and Centris.sp3 for example):
short.path<-shortest_paths(ppIg, from = "Xylocopa.sp", to = "Centris.sp3", output = "both")
short.path
## $vpath
## $vpath[[1]]
## + 5/26 vertices, named:
## [1] Xylocopa.sp Heteropterys.sp1 Epicharis.sp2
## [4] Diplopterys.pubipetala Centris.sp3
##
##
## $epath
## $epath[[1]]
## + 4/71 edges (vertex names):
## [1] Heteropterys.sp1 --Xylocopa.sp
## [2] Heteropterys.sp1 --Epicharis.sp2
## [3] Diplopterys.pubipetala--Epicharis.sp2
## [4] Diplopterys.pubipetala--Centris.sp3
##
##
## $predecessors
## NULL
##
## $inbound_edges
## NULL
If we want to emphasize the path between them, we have to create a new color vector for interactions (edges) the same way we did before for the species:
ecol <- rep("gray80", ecount(ppIg)) #creating a color vector for all the edges with the color "grey50"; the ecount fucntion tells you how many edges the network has
ecol[unlist(short.path$epath)] <- "blue" # coloring in blue only the path in which we are interested, that we calculed using the shortest_paths function.
plot(ppIg, vertex.color=vertex_attr(ppIg)$cor,vertex.label=NA,
vertex.size=2*igraph::degree(ppIg),
edge.width=2,
edge.color=ecol,
layout=l)
And finally, Igraph is kind of annoying to import stuff into, but it is very easy to export stuff from. You just export it as a normal plot:
pdf("FiguraLinda.pdf")
plot(ppIg, vertex.color=vertex_attr(ppIg)$cor,vertex.label=NA,
vertex.size=2*igraph::degree(ppIg),
edge.width=3*(edge_attr(ppIg)$weight)/1000,
edge.color="grey50",
edge.curved=0.3, #layout=l)
mark.groups= modulos, mark.border=NA)
dev.off()