【R】微生物网络构建简单方法汇总

偶然发现一本书:

里面刚好有微生物网络分析的内容，我就当个搬运工把代码整理了一下（注：书籍和示例数据在公众号PLANTOMIX后台回复“微生物组分析”即可获取下载链接噢）：

library(ggraph)
library(vegan)
library(MCL)
library(tidyverse)
library(qgraph)
library(SpiecEasi)
library(ggnet)

# load OTU table
load('../data/data.OTU.RData')

data.otu = data.otu[,seq(1,90,3)]
colnames(data.otu) = paste('S',1:ncol(data.otu), sep = '')
# 基于相异度的网络（Dissimilarity-Based Network）
diss.mat = vegan::vegdist(t(data.otu),
                           method = 'bray') %>%
        as.matrix()

diss.cutoff = 0.6 # 设置阈值
diss.adj = ifelse(diss.mat <= diss.cutoff, 1, 0)

# 从邻接矩阵中构建网络
diss.net = graph.adjacency(diss.adj,
                           mode = 'undirected',
                           diag = FALSE)

# 基于相关性的网络
cor.matrix = cor(diss.mat,method = 'pearson')
cor.cutoff = 0.3
cor.adj = ifelse(abs(cor.matrix) >= cor.cutoff,1,0) # 相关性矩阵转换成二元邻接矩阵
cor.net = graph.adjacency(cor.adj,
                          mode = 'undirected',
                          diag = FALSE)

build.cor.net =function(MB, method, num_perms, sig_level){

        taxa =dim(MB)[2]
        MB.mat =array(0, dim = c(taxa, taxa, num_perms + 1))
        # Perform permutation Constructing and Analyzing Microbiome Networks in R 251
        MBperm =permatswap(MB, "quasiswap", times = num_perms)
        # Convert to relative abundance
        MB.relative =MB / rowSums(MB)
        MB.mat[,,1]<-as.matrix(cor(MB.relative,method=method))
        for(p in 2:num_perms) {
                MBperm.relative=MBperm$perm[[p-1]] / rowSums(MBperm$perm[[p-1]])
                MB.mat[, , p] =as.matrix(cor(MBperm.relative, 
                                               method=method))
        }
        # Get p-values
        pvals=sapply(1:taxa,
                       function(i) sapply(1:taxa, function(j)
                               sum(MB.mat[i, j, 1]> MB.mat[i, j, 2:num_perms])))
        pvals=pvals / num_perms
        # p-value correction
        pvals_BH=array(p.adjust(pvals, method = "BH"),
                         dim=c(nrow(pvals), ncol(pvals)))
        # Build adjacency matrix
        adj.mat=ifelse(pvals_BH>= (1 - sig_level), 1, 0)
        # Add names to rows & cols
        rownames(adj.mat)=colnames(MB)
        colnames(adj.mat)=colnames(MB)
        # Build and return the network
        graph<-graph.adjacency(adj.mat,
                               mode="undirected",
                               diag=FALSE)
}

cor.net.2 = build.cor.net(as.matrix(data.otu), # 此处矩阵要求输入的是整数，小数会报错
                          method = 'pearson',
                          num_perms = 999,
                          sig_level = 0.01)

# Graphical Model Networks
ebic.cor = cor_auto(data.otu) # extended Bayesian information criterion
ebic.graph = EBICglasso(ebic.cor,
                        ncol(data.otu),
                        0.5)
ebic.qgnet = qgraph(ebic.graph,
                    DoNotPlot = FALSE) # build network
ebic.net = as.igraph(ebic.qgnet,
                     attributes=T) # convert to igraph

# local false discovery rate network
fdr.cor = cor_auto(data.otu)
fdr.graph = FDRnetwork(fdr.cor,
                       cutoff = 0.01,
                       method = 'pval')
fdr.qgnet = qgraph(fdr.graph, DoNotPlot = T)
fdr.net = as.igraph(fdr.qgnet)

# SparCC and SPIEC-EASI Networks
sparcc.matrix = sparcc(data.otu)
sparcc.cutoff = 0.3
sparcc.adj = ifelse(abs(sparcc.matrix$Cor) >= sparcc.cutoff,
                    1,0)
rownames(sparcc.adj) = colnames(data.otu)
colnames(sparcc.adj) = colnames(data.otu)

sparcc.net = graph.adjacency(sparcc.adj,
                             mode = 'undirected',
                             diag = FALSE)

# SPIEC-EASI network
spieceasi.matrix = spiec.easi(as.matrix(data.otu),
                              method = 'glasso',
                              lambda.min.ratio = 1e-2,
                              nlambda = 20,
                              icov.select.params = list(rep.num = 50))
rownames(spieceasi.matrix[["refit"]][["stars"]]) = colnames(data.otu)

spieceasi.net = graph.adjacency(spieceasi.matrix[["refit"]][["stars"]],
                                mode = 'undirected',
                                diag = FALSE)

# hub decetion
net = sparcc.net
net.cn = closeness(net)
net.bn = betweenness(net)
net.pr = page_rank(net)$vector
net.hs = hub_score(net)$vector

net.hs.sort = sort(net.hs, decreasing = T)
net.hs.top5 = head(net.hs.sort, 5)

# cluster decetion
wt = walktrap.community(net)
ml = multilevel.community(net)
membership(wt)
adj = as_adjacency_matrix(net)
mc = mcl(adj, addLoops = TRUE)

compare(membership(wt), membership(ml))
compare(membership(wt), mc$Cluster)
expected.cls =sample(1:5, vcount(net), replace = T) %>%
        as_membership
compare(expected.cls, membership(wt))
plot_dendrogram(wt)

modularity(net, membership(wt))

# Network Simulation
num.nodes = 50
regular.net = k.regular.game(num.nodes, k = 4)
random.net = erdos.renyi.game(num.nodes, p = 0.037)
smallworld.net = sample_smallworld(dim = 1, num.nodes, nei =
                                            2, p = 0.2)
scalefree.net = barabasi.game(num.nodes)

# 网络可视化
net = diss.net
wt = walktrap.community(net)

layout =  c(layout.fruchterman.reingold,
            layout_as_tree,
            layout_nicely,
            layout_on_sphere)

plot.igraph(net,
            vertex.size = 10,
            vertex.label = NA,
            edge.curved = T,
            edge.size = 1,
            layout = layout.fruchterman.reingold)

# 带聚类
ceb = cluster_edge_betweenness(net)
dendPlot(ceb, mode="hclust")
plot(ceb, net)

# 导出网络
#将网络图导出为"graphml"、"gml"格式，方便导入Gephi等软件中使用
#支持的格式有"edgelist", "pajek", "ncol", "lgl"
#"graphml", "dimacs", "gml", "dot", "leda"
write_graph(net, "../figures/g1.graphml", format="graphml")
write_graph(net, "../figures/g1.gephi", format="gml")
write_graph(net, '../figures/net.txt', format = 'edgelist')

大概图形是这样的：

我的 Cytoscape 和 Gephi 都罢工了，无法进行网络图绘制，大家自行尝试一下噢！