workingDir = "E:/WGCNA/analysis/ES"######working dir should be changed
setwd(workingDir)
getwd()
# Load the WGCNA package
library(WGCNA)


######################################################input data##########################
# The following setting is important, do not omit.
options(stringsAsFactors = FALSE)
#Read in the data set
fpkm = read.table("ES.txt",header=T,comment.char = "",check.names=F)#########file name can be changed
# Take a quick look at what is in the data set
dim(fpkm)
names(fpkm)
datExpr0 = as.data.frame(t(fpkm[,-1]))
names(datExpr0) = fpkm$Transcript_ID; 
rownames(datExpr0) = names(fpkm[,-1])
datExpr0

##check missing value
gsg = goodSamplesGenes(datExpr0, verbose = 3)
gsg$allOK

###############################Sample cluster############################## 
sampleTree = hclust(dist(datExpr0), method = "average")
# Plot the sample tree: Open a graphic output window of size 12 by 9 inches
# The user should change the dimensions if the window is too large or too small.
dev.new()
sizeGrWindow(13,9)
par(cex = 0.6)
par(mar = c(0,4,2,0))
plot(sampleTree, main = "Sample clustering to detect outliers", sub="", xlab="", cex.lab = 1.5,
     cex.axis = 1.5, cex.main = 2)

### Plot a line to show the cut
abline(h = 2000, col = "red")##



### Determine cluster under the line
clust = cutreeStatic(sampleTree, cutHeight = 2000, minSize = 10)
table(clust)


### clust 1 contains the samples we want to keep.
keepSamples = (clust==1)
datExpr0 = datExpr0[keepSamples, ]



##################################
save(datExpr0, file = "ES.RData")






#############################Visualize the scale-free index R2 under different power values#######################################
# Allow multi-threading within WGCNA. At present this call is necessary.
enableWGCNAThreads()

# Choose a set of soft-thresholding powers
powers = c(1:20 )

# Call the network topology analysis function
sft = pickSoftThreshold(datExpr0, powerVector = powers, verbose = 5)
sft
write.table(sft,file = "ES_sft.txt",quote = FALSE,sep = "\t")

# Plot the results:
sizeGrWindow(12, 6)
par(mfrow = c(1,2))
cex1 = 0.9
# Scale-free topology fit index as a function of the soft-thresholding power
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("ES_Scale independence"));
text(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2],
     labels=powers,cex=cex1,col="red");
# this line corresponds to using an R^2 cut-off of h
abline(h=0.9,col="red")
# Mean connectivity as a function of the soft-thresholding power
plot(sft$fitIndices[,1], sft$fitIndices[,5],
     xlab="Soft Threshold (power)",ylab="Mean Connectivity", type="n",
     main = paste("ES_Mean connectivity"))
text(sft$fitIndices[,1], sft$fitIndices[,5], labels=powers, cex=cex1,col="red")