用到的包
library(WGCNA)
library(igraph)
library(reshape2)
读表,后面的 %>%t()是对读的表格进行转置
OTU <- read.table("D:/aSbreviflora/data/N.txt", header=T, row.names=1) %>%t()
occor <- corAndPvalue(x=OTU,use = "pairwise.complete.obs", alternative="two.sided",method="spearman")
#提取相关矩阵的r、p值;
cor.cor<-occor$cor
cor.p<- occor$p
cor.cor[cor.p>0.001|abs(cor.cor)<0.9] = 0#p和cor的大小根据需要的以及最后算出来的节点和边数量进行设定,p最小是0.001
cor.cor2<-abs(cor.cor)
A<-data.frame(cor.cor2)
#用于计算Pi Zi值的时候用A,如果正常导出,用cor.cor
write.csv (A,paste0('D:/aSbreviflora/data/corr.csv'),row.names = T)#导出corre表格
r_value <- cbind( melt(cor.cor), melt(cor.p)$value ) #cbind()合并两个子集
r_value <- subset(r_value, r_value[,3]!=0)#subset 去除数据中符合某个条件的值
r_value <- r_value[r_value$Var1 != r_value$Var2, ]#删除第一列和第二列相同的行
#删除第一列与第二列互为相关的行
r_value$pair <- apply(r_value[, c("Var1", "Var2")], 1, function(x) paste(sort(x), collapse = "_"))#找到第一列与第二列相关的行,将其定义为pair列
r_value <- r_value[!duplicated(r_value$pair), ]#删除第一列与第二列相关,但是第二列和第一列也相关的行
r_value$pair <- NULL#删除多于的列
r_value <- na.omit(r_value)#去除NA
#对r表格增补绝对值、正负型等信息
abs <- abs(r_value$value)
linktype <- r_value$value
linktype[linktype>0]=1
linktype[linktype<0]=-1
r_value <- cbind(r_value, linktype)#并子集
names(r_value) <- c("Source","Target","r_value","p_value", "linktype")#names()重命名r、p表头
p_value <- melt(cor.p)
names(p_value) <- c("Source","Target","p_value")
导出边信息
write.csv(r_value,'D:/aSbreviflora/data/edge.csv')