加载工具包

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