Load example data:
library(microbiome)
data(dietswap)
pseq <- dietswap
# Keep only the prevalent taxa to speed up examples
pseq <- core(pseq, detection = 5^2, prevalence = 80/100)
pseq <- subset_samples(pseq, nationality == "AFR" & group == "DI" & bmi_group == "lean")
See the phyloseq tutorial for additional network visualization tools.
The widely reported compositionality bias in similarity measures can be fixed with SpiecEasi or SparCC; the implementations are available via the SpiecEasi package. Note that the execution is slow.
# Pick the OTU table
library(phyloseq)
otu <- abundances(pseq)
# SPIEC-EASI network reconstruction
# In practice, use more repetitions
library(devtools)
#install_github("zdk123/SpiecEasi")
library(SpiecEasi) #install_github("zdk123/SpiecEasi")
net <- spiec.easi(t(otu), method='mb', lambda.min.ratio=1e-2, nlambda=5, icov.select.params=list(rep.num=1))
## Create graph object
n <- net$refit
#colnames(n) <- rownames(n) <- rownames(otu)
# Network format
library(network)
#netw <- network(as.matrix(n), directed = FALSE)
# igraph format
library(igraph)
# ig <- graph.adjacency(n, mode='undirected', add.rownames = TRUE)
# Network layout
# coord <- layout.fruchterman.reingold(ig)
## set size of vertex to log2 mean abundance
# vsize <- log2(rowMeans(otu))
# Visualize the network
# print(plot(ig, layout = coord, vertex.size = vsize, vertex.label = names(vsize)))
Investigate degree distribution with the following:
#dd <- degree.distribution(ig)
#plot(0:(length(dd)-1), dd, ylim = c(0,.35), type = 'b', ylab = "Frequency", xlab = "Degree", main = "Degree Distributions")
Visualize the network with ggnet2:
library(GGally)
#library(ggnet)
library(network)
library(sna)
library(ggplot2)
library(intergraph) # ggnet2 works also with igraph with this
phyla <- map_levels(rownames(otu),
from = "Genus", to = "Phylum",
tax_table(pseq))
#netw %v% "Phylum" <- phyla
#p <- ggnet2(netw, color = "Phylum", label = TRUE, label.size = 2)
print(p)