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: viridisLite
library(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, union

We 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...
## done

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