博主认为重新写个整体的脚本比较合适,一段一段写容易混淆。层次聚类先将矩阵转置,然后求powers(阈值),根据阈值进行聚类:
library(WGCNA)
options(stringsAsFactors=FALSE)
enableWGCNAThreads()
myfile=read.table("da1.nom",sep="\t",header=TRUE)
mydata=as.data.frame(t(myfile[,-c(1)]))
names(mydata)=myfile$ENST
rownames(mydata)=names(myfile)[-c(1)]
powers = c(c(1:20), seq(from=22, to=40, by=2))
sft = pickSoftThreshold(mydata, powerVector=powers, verbose=5)
sizeGrWindow(10,5)
par(mfrow = c(1,2))
plot(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2],xlab="Soft Threshold(power)",ylab="Scale Free Topology Model Fit, signed R^2",type="n",main=paste("Scale independence"))
text(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2],labels=powers,cex=0.9,col="red")
plot(sf