Redes no Igraph

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 =)

Atributes

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.

Plotting

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")

Analysis

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)

Exporting

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()