加载工具包
library(scales)
library(RColorBrewer)
library(igraph)
library(dplyr)
library(ggrepel)
library(reshape2)
加载数据
load('../00_load_data/.RData')
rm(physeq, rarefy)
otutab = data.frame(t(t(otutab)/colSums(otutab)))
source('zipi.R')
定义一些颜色
col_g <- "#C1C1C1"
cols <- c("#DEB99B","#5ECC6D","#5DAFD9","#F16E1D","#6E4821","#A4B423",
"#DC95D8","#326530","#F0027F","#E6AB02","#F96C72")
show_col(cols,ncol = 1)
筛选优势物种
top_otutab = arrange(otutab, desc(rowMeans(otutab)))[1:round(1*nrow(otutab)),] %>%
.[rowSums(.!=0)>=0.1*ncol(.),]
topTax = c('Alphaproteobacteria','Gammaproteobacteria','Actinobacteriota','Cyanobacteria',
'Bacteroidota','Acidobacteriota','Chloroflexi','Planctomycetota')
构建函数计算相关性和共现性网络
network <- function(x){
occor<-WGCNA::corAndPvalue(t(x),method = 'spearman')
mtadj<-multtest::mt.rawp2adjp(unlist(occor$p),proc='BH')
adpcor<-mtadj$adjp[order(mtadj$index),2]
occor.p<-matrix(adpcor,dim(t(x)/colSums(x))[2])
## R value
occor.r<-occor$cor
diag(occor.r) <- 0
occor.r[occor.p>0.01|abs(occor.r)<0.6] = 0
occor.r[is.na(occor.r)]=0
g <- graph.adjacency(occor.r, weighted = TRUE, mode = 'undirected')
# 删除自相关
g <- simplify(g)
# 删除孤立节点
g <- delete.vertices(g, which(degree(g)==0) )
return(g)
}
network_stat = function(x) {
result = data.frame(
node.number = length(V(x)), # number of nodes
edges.number = length(E(x)), # number of edges
average.degree = length(E(x))/length(V(x)), # average degree
Clusting.coeff = transitivity(x), # clustering coefficient
aver.path.len = average.path.length(x), # average.path.length
density = graph.density(x)# graph density
)
return(result)
}
edge_nodes = function(x){
edges = data.frame(get.edgelist(x))
names(edges) = c('Source','Target')
edges$level= ifelse(get.edge.attribute(x)[[1]]>0,'Positive','Negative')
edges$weight = abs(get.edge.attribute(x)[[1]])
write.csv(edges, file = 'edges.csv', row.names = F)
E(x)$correlation <- E(x)$weight
E(x)$weight <- abs(E(x)$weight)
E(x)$width <- abs(E(x)$weight)
set.seed(007)
V(x)$modularity <- membership(cluster_fast_greedy(x))
V(x)$label <- V(x)$name
V(x)$label <- NA
size = data.frame(proportion = rowMeans(top_otutab))
V(x)$size <- size[V(x)$name,'proportion']*100
V(x)$taxonomy = taxa[V(x)$name,'mixed']
zp = ZiPi(x,modules=V(x)$modularity)
nodes = data.frame(id =V(x)$name, size = V(x)$size) %>%
merge(.,taxa, by.x = 'id', by.y = 'otuid')
nodes = merge(nodes, zp, by.x = 'id',by.y='names')
nodes$mixed = ifelse(nodes$mixed %in% topTax, as.character(nodes$mixed), 'Others')
write.csv(nodes, file ='nodes.csv', row.names = F)
}
运行函数，生成结果
net = network(top_otutab)
net.stat = network_stat(net)
net.attr = edge_nodes(net)
输出结果，在gephi中进行精加工
Network.jpg
