source("https://bioconductor.org/biocLite.R")
biocLite(c("GO.db","preprocessCore","impute"))
install.packages(c("matrixStats","Hmisc","splines","foreach","doParallel","reshape","fastcluster","dynamicTreeCut","survival"))
install.packages("WGCNA")
setwd("D:/TCGA_coexpression/WGCNA")
library(WGCNA)
allowWGCNAThreads() #查看核心数
ALLOW_WGCNA_THREADS=4 #设置为4核心运行
rt <- read.table(file="merge.txt",sep="\t",row.names=1,header=T,check.names=F,quote="!")
datExpr <- t(rt)
power1 <- c(seq(1,10,by=1),seq(12,30,by=2))
RpowerTable <- pickSoftThreshold(datExpr,powerVector=power1)[[2]]
#上面的输出结果
# Power SFT.R.sq slope truncated.R.sq mean.k. median.k. max.k.
#1 1 0.115 -0.643 0.860 639.000 6.13e+02 1360.0
#2 2 0.527 -1.160 0.816 209.000 1.73e+02 653.0
#3 3 0.693 -1.030 0.781 96.300 6.48e+01 359.0
#4 4 0.945 -0.919 0.954 54.000 2.95e+01 236.0
#5 5 0.929 -1.080 0.963 34.200 1.44e+01 195.0
#6 6 0.945 -1.160 0.986 23.400 7.41e+00 167.0
#7 7 0.951 -1.210 0.989 16.900 4.02e+00 147.0
#8 8 0.967 -1.250 0.999 12.800 2.29e+00 130.0
#9 9 0.975 -1.270 0.998 9.930 1.34e+00 117.0
#10 10 0.979 -1.280 0.998 7.920 8.25e-01 106.0
#11 12 0.982 -1.280 0.997 5.340 3.62e-01 88.9
#12 14 0.981 -1.270 0.991 3.820 1.75e-01 76.0
#13 16 0.983 -1.270 0.992 2.860 8.68e-02 66.0
#14 18 0.982 -1.260 0.989 2.210 4.28e-02 57.9
#15 20 0.981 -1.250 0.988 1.750 2.24e-02 51.4
#16 22 0.961 -1.250 0.966 1.420 1.21e-02 46.0
#17 24 0.961 -1.240 0.961 1.170 6.49e-03 41.4
#18 26 0.977 -1.230 0.982 0.981 3.44e-03 37.5
#19 28 0.981 -1.240 0.987 0.832 1.87e-03 34.2
#20 30 0.968 -1.240 0.970 0.713 1.01e-03 31.3
cex1 <- 0.7
pdf(file="softThresholding.pdf")
par(mfrow=c(1,2))
plot(RpowerTable[,1],-sign(RpowerTable[,3])*RpowerTable[,2],xlab="Soft Threshold (power)",ylab="Scale Free Topology Model Fit,signed R^2",type="n")
text(RpowerTable[,1],-sign(RpowerTable[,3])*RpowerTable[,2],labels=power1,cex=cex1,col="red")
abline(h=0.85,col="red")
plot(RpowerTable[,1],RpowerTable[,5],xlab="Soft Threshold (power)",ylab="Mean Connectivity",type="n")
text(RpowerTable[,1],RpowerTable[,5],labels=power1,cex=cex1,col="red")
dev.off()
beta1 <- 9
ADJ <- adjacency(datExpr,power=beta1)
visual <- exportNetworkToCytoscape(ADJ,edgeFile="edge.txt",nodeFile="node.txt",threshold=0.8)
disableWGCNAThreads()
