library(WGCNA)  
options(stringsAsFactors = FALSE)   
acc<-read.csv("data_of_TPM_for_analysis.csv",sep="\t")
sample<-read.csv("acc_clinical_data.csv",sep = "\t")
acc$id<-row.names(brca)
both<-merge(sample,acc,by.x = "id",by.y = "id")
both<-na.omit(both)
sample<-both[,1:8]
acc<-both[,c(1,9:c(dim(both)[2]))]
row.names(acc)<-acc$id
acc<-acc[,-1]
row.names(sample)<-sample$id
sample<-sample[,-1]
write.table(acc,file="acc_tpm.csv",sep = "\t",quote = F)
write.table(sample,file="acc_clinical.csv",sep = "\t",quote = F)

femData = acc 
datExpr0 = femData  
sampleTree = hclust(dist(datExpr0), method = "average")  

pdf("Sample_clustering_to_detect_outliers.pdf")
par(cex = 0.6);par(mar = c(0,4,2,0))
plot(sampleTree,main = "Sample clustering to detect outliers",
     sub = "",xlab = "",cex.lab = 1,cex.axis = 1,cex.main = 1)
abline(h = 1500000,col = "red") #save image as Sample cluster 1st
clust = cutreeStatic(sampleTree,cutHeight = 1500000,minSize = 10)
table(clust)
keepSamples = (clust == 1)
nGenes = ncol(acc)
nSamples = ncol(acc)
traitColors = numbers2colors(sample,signed = F)
plotDendroAndColors(sampleTree, traitColors,
                    groupLabels = names(sample),
                    main = "Sample dendrogram and trait heatmap")
dev.off()
####################################

powers = c(c(1:10), seq(from = 12, to=20, by=2)) 
sft = pickSoftThreshold(datExpr0, powerVector = powers, verbose = 5)  

    sizeGrWindow(9, 5)  
    par(mfrow = c(1,2))
    cex1 = 0.9  
    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=cex1,col="red")  
    abline(h=0.90,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=cex1,col="red") 
######################################################################
library(WGCNA)
options(stringsAsFactors = FALSE)
femData = read.csv("acc_tpm.csv",sep="\t")
datExpr0 = femData
gsg = goodSamplesGenes(datExpr0, verbose = 3)
gsg$allOK

softPower = 6 #depend on power analysis above
adjacency = adjacency(datExpr0, power = softPower)
TOM = TOMsimilarity(adjacency)
dissTOM = 1-TOM
geneTree = hclust(as.dist(dissTOM), method = "average")
sizeGrWindow(12,9)
plot(geneTree, xlab="", sub="", main = "Gene clustering on TOM-based dissimilarity",
     labels = FALSE, hang = 0.04)
minModuleSize = 30
dynamicMods = cutreeDynamic(dendro = geneTree, distM = dissTOM,
                            deepSplit = 2, pamRespectsDendro = FALSE,
                            minClusterSize = minModuleSize)
table(dynamicMods)
dynamicColors = labels2colors(dynamicMods)
table(dynamicColors)
sizeGrWindow(8,6)
plotDendroAndColors(geneTree, dynamicColors, "Dynamic Tree Cut",
                    dendroLabels = FALSE, hang = 0.03,
                    addGuide = TRUE, guideHang = 0.05,
                    main = "Gene dendrogram and module colors")
MEList = moduleEigengenes(datExpr0, colors = dynamicColors)
MEs = MEList$eigengenes
MEDiss = 1-cor(MEs)
METree = hclust(as.dist(MEDiss), method = "average")
sizeGrWindow(7, 6)
plot(METree, main = "Clustering of module eigengenes",
     xlab = "", sub = "")
MEDissThres = 0.25
abline(h=MEDissThres, col = "red")
merge = mergeCloseModules(datExpr0, dynamicColors, cutHeight = MEDissThres, verbose = 3)
mergedColors = merge$colors
mergedMEs = merge$newMEs
sizeGrWindow(12, 9)
plotDendroAndColors(geneTree, cbind(dynamicColors, mergedColors),
                    c("Dynamic Tree Cut", "Merged dynamic"),
                    dendroLabels = FALSE, hang = 0.03,
                    addGuide = TRUE, guideHang = 0.05)
moduleColors = mergedColors
colorOrder = c("grey", standardColors(50))
moduleLabels = match(moduleColors, colorOrder)-1
MEs = mergedMEs
save(MEs, moduleLabels, moduleColors, geneTree, file = "networkConstruction-stepByStep.RData")
########################################################################################
library(WGCNA)
options(stringsAsFactors = FALSE)
femData = read.csv("acc_tpm.csv",sep="\t")
datExpr0 = femData
lnames = load(file = "networkConstruction-stepByStep.RData")
lnames
nGenes = ncol(datExpr0)
nSamples = nrow(datExpr0)

dissTOM = 1-TOMsimilarityFromExpr(datExpr0, power = 4)

nSelect = 500
set.seed(10)
select = sample(nGenes, size = nSelect)
selectTOM = dissTOM[select, select]
selectTree = hclust(as.dist(selectTOM), method = "average")
selectColors = moduleColors[select]
sizeGrWindow(9,9)
plotDiss = selectTOM^7
diag(plotDiss) = NA
TOMplot(plotDiss, selectTree, selectColors, main = "Network heatmap plot, all genes")

###############################################################
library(WGCNA)
options(stringsAsFactors = FALSE)
femData = read.csv("acc_tpm.csv",sep="\t")
datExpr0 = femData
lnames = load(file = "networkConstruction-stepByStep.RData")

TOM = TOMsimilarityFromExpr(datExpr0, power = 6)
annot = read.csv(file = "GeneAnnotation.csv",sep="\t")
module = "blue"
probes = names(datExpr0)
inModule = (moduleColors==module)
modProbes = probes[inModule]
modTOM = TOM[inModule, inModule]
dimnames(modTOM) = list(modProbes, modProbes)
vis = exportNetworkToVisANT(modTOM,
file = paste("VisANTInput-", module, ".txt", sep=""),
weighted = TRUE,
threshold = 0,
probeToGene = data.frame(annot$Gene, annot$gene_symbol) )
inModule = is.finite(match(moduleColors, module))
modProbes = probes[inModule]
modGenes = annot$gene_symbol[match(modProbes, annot$Gene)]
modTOM = TOM[inModule, inModule]
dimnames(modTOM) = list(modProbes, modProbes)
cyt = exportNetworkToCytoscape(modTOM,
edgeFile = paste("CytoscapeInput-edges-", paste(module, collapse="-"), ".txt", sep=""),
nodeFile = paste("CytoscapeInput-nodes-", paste(module, collapse="-"), ".txt", sep=""),
weighted = TRUE,
threshold = 0.02,
nodeNames = modProbes,
altNodeNames = modGenes,
nodeAttr = moduleColors[inModule])

###############################################################