#### figure1
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.8,col="red")
plot(sft$fitIndices[,1], sft$fitIndices[,5],xlab="Soft Threshold (power)",ylab="Mean Connectivity", type="n",main = paste("Mean connectivity"))
text(sft$fitIndices[,1], sft$fitIndices[,5], labels=powers, cex=0.8,col="red")

### figure2  
plotDendroAndColors(geneTree, cbind(module.colours, mergedColors),
                    c("Dynamic Tree Cut", "Merged dynamic"),
                    dendroLabels = FALSE, hang = 0.03,
                    addGuide = TRUE, guideHang = 0.05)
                    
### figure3 
par(mar = c(2, 6, 2, 5))
labeledHeatmap(Matrix = moduleTraitCor4,
               xLabelsAngle = 0,
               cex.lab = 0.7,
               xLabels = colnames(sample),
               yLabels = names(MET),
               ySymbols = names(MET),
               colorLabels = FALSE,
               colors = blueWhiteRed(100),
               #textMatrix = textMatrix,
               setStdMargins = FALSE,
               cex.text = 0.5,
               zlim = c(-0.2,0.2),
               yColorWidth=0.1,
               xColorWidth = 0.05,
               main = paste("Module-trait relationships"))
               
###figure4 
obs<-paste("sam",1:108,sep="")
sample<-as.data.frame(diag(x=1,nrow = 108))
rownames(sample)<-obs
colnames(sample)<-1:108

MEList = moduleEigengenes(WGCNA_matrix, colors = module.colours)
MEs = MEList$eigengenes
MET = orderMEs(MEs)

#-------------merge
# Call an automatic merging function
merge = mergeCloseModules(WGCNA_matrix, module.colours, cutHeight = MEDissThres, verbose = 3)
# The merged module colors
mergedColors = merge$colors
# Eigengenes of the new merged modules:
mergedMEs = merge$newMEs
MET = orderMEs(mergedMEs)
#-------------merge

moduleTraitCor = cor(MET, sample, use = "p")
moduleTraitCor1<-as.matrix(apply(moduleTraitCor[,1:20],1,mean))
moduleTraitCor2<-as.matrix(apply(moduleTraitCor[,21:70],1,mean))
moduleTraitCor3<-as.matrix(apply(moduleTraitCor[,71:108],1,mean))
moduleTraitCor4<-data.frame(cbind(moduleTraitCor1,moduleTraitCor2,moduleTraitCor3))

par(mar = c(2, 9, 2, 1),family = "sans")   
labeledHeatmap(Matrix = moduleTraitCor4,
               xLabelsAngle = 0,
               cex.lab = 0.9,
               xLabels = colnames(sample),
               yLabels = names(MET),
               ySymbols = names(MET),
               colorLabels = FALSE,
               colors = blueWhiteRed(100),
               #textMatrix = textMatrix,
               setStdMargins = FALSE,
               cex.text = 0.5,
               zlim = c(-0.2,0.2),
               yColorWidth=0.08,
               xColorWidth = 0.05,
               main = paste("Module-trait relationships"))

#### figure5 
datExpr<-WGCNA_matrix[,]
datExpr_tree<-hclust(dist(datExpr), method = "average")
par(mar = c(0,5,2,0))
plot(datExpr_tree, main = "Sample clustering", sub="", xlab="", cex.lab = 2, 
     cex.axis = 1, cex.main = 1,cex.lab=1)
sample_colors <- numbers2colors(as.numeric(modules_5$group),colors = c("green","blue","red"),signed = FALSE)
par(mar = c(1,4,3,1),cex=0.8)
plotDendroAndColors(datExpr_tree, sample_colors,
                    #groupLabels = modules_5$group,
                    dendroLabels = FALSE,
                    cex.dendroLabels = 0.8,
                    marAll = c(1, 4, 3, 1),
                    cex.rowText = 0.01,
                    main = "Sample dendrogram and trait heatmap")   
                    
#figure6
#load rf_18.RData
print(rf_19_1000)
plot(rf_19_1000)   

#figure7
varImpPlot(rf_19_1000)
num=50
sig<-as.matrix(importance(rf_19_1000,type = 2))
sig2<-as.matrix(sig[order(sig,decreasing = T),])
try= rep(NA,num)
try[1:num]<-sub("gene","",rownames(sig2)[1:num])

name2<-as.matrix(colnames(modules_7_1000)[as.numeric(try)]) 

sig<-as.matrix(importance(rf_19_1000,type = 1))
sig2<-as.matrix(sig[order(sig,decreasing = T),])
try= rep(NA,num)
try[1:num]<-sub("gene","",rownames(sig2)[1:num])

name1<-as.matrix(colnames(modules_7_1000)[as.numeric(try)]) 


#figure8 
s2<-cbind(rownames(s),s)
s3<-s[,colnames(modules_7_1000)]
s3<-s3[colnames(modules_7_1000),]
s3<-cbind(rownames(s3),s3)

write.table(s3, "C:/Users/Administrator/Desktop/paper/figure/s3.csv", sep=",", row.names=FALSE)
write.table(names(MET), "C:/Users/Administrator/Desktop/paper/figure/MET.csv", sep=",", row.names=FALSE)



##table 2
modules_5<-WGCNA_matrix[,] #5000¸ö×Ü»ùÒò
colnames(modules_5)<-paste("gene",c(1:ncol(modules_5)),sep="")
group = rep(NA,108)
modules_5<-as.data.frame(cbind(modules_5,group)) 
modules_5$group[1:20]<-"c"
modules_5$group[21:70]<-"MMA" 
modules_5$group[71:108]<-"SA"
modules_5$group<-as.factor(modules_5$group)

rf <- randomForest(group ~ ., ntree=1000,data = modules_5,importance=TRUE,proximity=TRUE,
                   method = "repeatedcv",number = 10, repeats = 10,
                   mtry=1)
print(rf)              
#load rf_18.RData
print(rf_19_1000)
plot(rf_19_1000)