rm(list=ls())
setwd("C:/rshuju")
otu <- read.csv("OTUZUIXIN2.csv",row.name = 1) 

a<-rowSums(otu)#计算相对丰度总和并付给a
u <- which(a > 0.0001)
otu1 <- otu[u,]
otu2=otu1
otu1[otu1>0] <- 1#将丰度值大于0的值替换为1，便于计算不同OTU的出现率；
otu3<- otu2[which(rowSums(otu1) >=5), ]#例如只保留在 3 个及以上样本中出现的属
write.csv(otu3,"jieguo2.csv")
library(psych)
# 读取otu-sample矩阵，行为sample，列为otu
#otu <- read.csv(file.choose(), head=T, row.names=1)
otu4 <- read.csv("jieguo2.csv",row.name = 1) 
occor <- corr.test(t(otu4),method="spearman",adjust="fdr",alpha=0.05)#adjust为校准r
occor.r <-  occor$r # 取相关性矩阵R值
occor.p <- occor$p # 取相关性矩阵p值
# p 值校正，这里使用 BH 法校正 p 值
p <- p.adjust(occor.p, method = 'BH') 
# 确定物种间存在相互作用关系的阈值，将相关性R矩阵内不符合的数据转换为0
occor.r[occor.p>0.05|abs(occor.r)<0.8] = 0
diag(occor.r) <- 0
#将occor.r保存为csv文件
write.csv(occor.r,file="xiangguanxijieguo2.csv")
#根据上述筛选的 r 值和 p 值保留数据
z <- occor.r * occor.p
diag(z) <- 0    #将相关矩阵中对角线中的值（代表了自相关）转为 0
head(z)[1:6,1:6]
z
write.csv(z,"huitu2.csv")
##获得网络
library(igraph)
#将邻接矩阵转化为 igraph 网络的邻接列表
#构建含权的无向网络，权重代表了微生物属间丰度的 spearman 相关系数 
g <- graph.adjacency(z, weighted = TRUE, mode = 'undirected')
g
#自相关也可以通过该式去除
g <- simplify(g)
#孤立节点的删除（删除度为 0 的节点）
g <- delete.vertices(g, names(degree(g)[degree(g) == 0]))
#该模式下，边权重代表了相关系数
#由于权重通常为正值，因此最好取个绝对值，相关系数重新复制一列
E(g)$correlation <- E(g)$weight
E(g)$weight <- abs(E(g)$weight)
E(g)$cor[E(g)$correlation>0]<- 1
E(g)$cor[E(g)$correlation<0]<- -1
tax <- read.csv(file.choose(), row.name = 1)
tax <- tax[as.character(V(g)$name), ]
V(g)$Domain <- tax$Domain
V(g)$Phylum <- tax$Phylum
g
plot(g)
write.graph(g, 'xiyou.graphml', format = 'graphml')