Script FRUITY

##############################################################
####### WEIGHTED (GENE) CORRELATION NETWORK ANALYSIS##########
##############################################################

##WGCNA is available from CRAN but can be badly compiled the script below allows manual loading

install.packages(c("matrixStats", "Hmisc", "splines", "foreach", "doParallel", "fastcluster", "dynamicTreeCut", "survival"))
source("http://bioconductor.org/biocLite.R")
biocLite(c("GO.db", "preprocessCore", "impute"))
source("https://bioconductor.org/biocLite.R")
    biocLite("BiocUpgrade")

#also install flashClust from CRAN

library(impute)
library(WGCNA)
library(flashClust)

# loading data, again in separate files. Select your own method.
d()
#rename working directory
workingDir = "/Users/Natasha.spadafora/OneDrive - Markes International Limited/R software/Peach VOCs Sensor Phytochem 2017"
#C:\R software\Fruitydata2018
#C:\TOF-DS\FRUITY2017\FRUITY WCNA\R tempate Analysis\FRUITY MARKES INT
#C:\R software\Sqrt VOCas trait Fruity2017 by CTM 2019
# set new working directory:
setwd(workingDir)

datCompds0=read.csv("Sqrtalldata2017NoSensorialNogeneNoMDA.csv", header = T)
#datCompds0=read.csv("D:/BI3153/NormAreaSqrt.csv", header = T)# equivalent to YData

dim(datCompds0)
names(datCompds0)

traitData = read.csv("ParamSensorialall.csv")
#traitData = read.csv("D:/BI3153/xdata2.csv");# equivalent to XData

dim(traitData)
names(traitData)

sampleTree = flashClust(dist(datCompds0), 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.
sizeGrWindow(12,10) #this command is redundant in RStudio
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 = 40, col = "red");

# Usually with scent data there is no cut required as there's no outgroup; if there is go to appendix 1

# Convert traits to a color representation: white means low, red means high, grey means missing entry

traitColors = numbers2colors(traitData ,signed = FALSE);
# Plot the sample dendrogram and the colors underneath.
plotDendroAndColors(sampleTree, traitColors,
groupLabels = names(traitData),
main = "Sample dendrogram and trait heatmap") #if reclustering is required the term 'sampletree' needs to be amended to sampletree2 (see appendix I)

# Choose a set of soft-thresholding powers
powers = c(c(1:10), seq(from = 12, to=20, by=2))
# Call the network topology analysis function
sft = pickSoftThreshold(datCompds0, powerVector = powers, verbose = 5)
# Plot the results:
sizeGrWindow(9, 5)
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("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.90,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("Mean connectivity"))
text(sft$fitIndices[,1], sft$fitIndices[,5], labels=powers, cex=cex1,col="red")
# all above is plot only!
#SELECT SOFTPOWER TO GAIN LARGEST SCALE FREE TOPOLOGY FIT AND LOWWEST MEAN CONNECTIVITY TO CONTINUE

# Calculating adjeacency with soft power 4
softPower= 5
adjacency = adjacency(datCompds0, power = softPower)

# Turn adjacency into topological overlap
TOM = TOMsimilarity(adjacency);
dissTOM = 1-TOM

# Call the hierarchical clustering function
geneTree = flashClust(as.dist(dissTOM), method = "average");
# Plot the resulting clustering tree (dendrogram)
sizeGrWindow(8,5)
plot(geneTree, xlab="", sub="", main = "Gene clustering on TOM-based dissimilarity",
labels = FALSE, hang = 0.04)

# We like large modules, so we set the minimum module size relatively high:
minModuleSize =4;
# Module identification using dynamic tree cut:
dynamicMods = cutreeDynamic(dendro = geneTree, distM = dissTOM,
deepSplit =3, pamRespectsDendro = FALSE,
minClusterSize = minModuleSize);
table(dynamicMods)
# Convert numeric lables into colors
dynamicColors = labels2colors(dynamicMods)
table(dynamicColors)

# Plot the dendrogram and colors underneath
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")



# Calculate eigengenes
MEList = moduleEigengenes(datCompds0, colors = dynamicColors)
MEs = MEList$eigengenes
# Calculate dissimilarity of module eigengenes
MEDiss = 1-cor(MEs);
# Cluster module eigengenes
METree = flashClust(as.dist(MEDiss), method = "average");
# Plot the result
sizeGrWindow(7, 6)
plot(METree, main = "Clustering of module eigengenes",
xlab = "", sub = "")

MEDissThres = 0.01
# Plot the cut line into the dendrogram
abline(h=MEDissThres, col = "red")
###### IF MERGER IS WANTED, E.G. COMBINING MODULES GO TO APPENDIX 2

# Define numbers of genes and samples
nGenes = ncol(datCompds0);
nSamples = nrow(datCompds0);
# Recalculate MEs with color labels
#MEs0 = moduleEigengenes(datCompds0, moduleColors)$eigengenes
#MEs = orderMEs(MEs0)
moduleTraitCor = cor(MEs, traitData, use = "p");
moduleTraitPvalue = corPvalueStudent(moduleTraitCor, nSamples);

# START PLOTTING HEAT MAP MODULE-TRAIT RELATION AND SIGNIFICANCE


sizeGrWindow(10,6)
# Will display correlations and their p-values
textMatrix = paste(signif(moduleTraitCor, 2), "\n(",
signif(moduleTraitPvalue, 1), ")", sep = "");
dim(textMatrix) = dim(moduleTraitCor)
par(mar = c(6, 8.5, 3, 3));
# Display the correlation values within a heatmap plot
labeledHeatmap(Matrix = moduleTraitCor,
xLabels = names(traitData),
yLabels = names(MEs),
ySymbols = names(MEs),
colorLabels = FALSE,
colors = greenWhiteRed(50),
textMatrix = textMatrix,
setStdMargins = FALSE,
cex.text = 0.5,
zlim = c(-1,1),
main = paste("Module-trait relationships"))

#END PLOTTING##

#PREPARING OUTPUT TABLE#

# Define variable 'Sex' containing the 'Sex' column of datTrait

S1 = as.data.frame(traitData$S1);
names(S1) = "S1"
S2 = as.data.frame(traitData$S2)
names(S2)="S2"
S3 = as.data.frame(traitData$S3)
names(S3)="S3"
S4 = as.data.frame(traitData$S4)
names(S4)="S4"
S5 = as.data.frame(traitData$S5)
names(S5)="S5"
S6 = as.data.frame(traitData$S6)
names(S6)="S6"
S7 = as.data.frame(traitData$S7)
names(S7)="S7"
S8 = as.data.frame(traitData$S8)
names(S8)="S8"
S9 = as.data.frame(traitData$S9)
names(S9)="S9"
# names (colors) of the modules
modNames = substring(names(MEs), 5)
geneModuleMembership = as.data.frame(cor(datCompds0, MEs, use = "p"));
MMPvalue = as.data.frame(corPvalueStudent(as.matrix(geneModuleMembership), nSamples));
names(geneModuleMembership) = paste("MM", modNames, sep="");
names(MMPvalue) = paste("p.MM", modNames, sep="");
geneTraitSignificance = as.data.frame(cor(datCompds0, S1, use = "p"));
GSPvalue = as.data.frame(corPvalueStudent(as.matrix(geneTraitSignificance), nSamples));
names(geneTraitSignificance) = paste("GS.", names(S1), sep="");
names(GSPvalue) = paste("p.GS.", names(S1), sep="");

geneTraitSignificance1 = as.data.frame(cor(datCompds0, S2, use = "p"));
GSPvalue1 = as.data.frame(corPvalueStudent(as.matrix(geneTraitSignificance1), nSamples));
names(geneTraitSignificance1) = paste("GS.", names(S2), sep="");
names(GSPvalue1) = paste("p.GS.", names(S2), sep="");

geneTraitSignificance2 = as.data.frame(cor(datCompds0,S3, use = "p"));
GSPvalue2 = as.data.frame(corPvalueStudent(as.matrix(geneTraitSignificance2), nSamples));
names(geneTraitSignificance2) = paste("GS.", names(S3), sep="");
names(GSPvalue2) = paste("p.GS.", names(S3), sep="");

geneTraitSignificance3 = as.data.frame(cor(datCompds0,S4, use = "p"));
GSPvalue3 = as.data.frame(corPvalueStudent(as.matrix(geneTraitSignificance3), nSamples));
names(geneTraitSignificance3) = paste("GS.", names(S4), sep="");
names(GSPvalue3) = paste("p.GS.", names(S4), sep="");

geneTraitSignificance4 = as.data.frame(cor(datCompds0,S5, use = "p"));
GSPvalue4 = as.data.frame(corPvalueStudent(as.matrix(geneTraitSignificance4), nSamples));
names(geneTraitSignificance4) = paste("GS.", names(S5), sep="");
names(GSPvalue4) = paste("p.GS.", names(S5), sep="");

geneTraitSignificance5 = as.data.frame(cor(datCompds0,S6, use = "p"));
GSPvalue5 = as.data.frame(corPvalueStudent(as.matrix(geneTraitSignificance5), nSamples));
names(geneTraitSignificance5) = paste("GS.", names(S6), sep="");
names(GSPvalue5) = paste("p.GS.", names(S6), sep="");

geneTraitSignificance6 = as.data.frame(cor(datCompds0,S7, use = "p"));
GSPvalue6 = as.data.frame(corPvalueStudent(as.matrix(geneTraitSignificance6), nSamples));
names(geneTraitSignificance6) = paste("GS.", names(S7), sep="");
names(GSPvalue6) = paste("p.GS.", names(S7), sep="");

geneTraitSignificance7 = as.data.frame(cor(datCompds0,S8, use = "p"));
GSPvalue7 = as.data.frame(corPvalueStudent(as.matrix(geneTraitSignificance7), nSamples));
names(geneTraitSignificance7) = paste("GS.", names(S8), sep="");
names(GSPvalue7) = paste("p.GS.", names(S8), sep="");

geneTraitSignificance8 = as.data.frame(cor(datCompds0,S9, use = "p"));
GSPvalue8 = as.data.frame(corPvalueStudent(as.matrix(geneTraitSignificance8), nSamples));
names(geneTraitSignificance8) = paste("GS.", names(S9), sep="");
names(GSPvalue8) = paste("p.GS.", names(S9), sep="");






WCNASensorial2017=data.frame(moduleColors=dynamicColors,geneModuleMembership,MMPvalue,
geneTraitSignificance,GSPvalue,geneTraitSignificance1,GSPvalue1,geneTraitSignificance2,GSPvalue2,geneTraitSignificance3,GSPvalue3,
geneTraitSignificance4,GSPvalue4,geneTraitSignificance5,GSPvalue5,geneTraitSignificance6,GSPvalue6,geneTraitSignificance7,GSPvalue7,
geneTraitSignificance8,GSPvalue8)
write.csv(WCNASensorial2017, file="WCNA_TraitSensor_1.csv")

#####OUTPUT TABLE CREATED######

# ASSESSMENT CORRELATION MODULE V TRAIT SIGNIFICANCE OF COMPOUNDS WITHIN MODULES
moduleColors = dynamicColors
module = "blue"
column = match(module, modNames);
moduleGenes = moduleColors==module;
sizeGrWindow(7, 7);
par(mfrow = c(1,1));
verboseScatterplot(abs(geneModuleMembership[moduleGenes, column]),
abs(geneTraitSignificance[moduleGenes, 1]),
xlab = paste("Module Membership in", module, "module"),
ylab = "Compound significance for Age",
main = paste("Module membership vs. gene significance\n"),
cex.main = 1.2, cex.lab = 1.2, cex.axis = 1.2, col = module)

names(datCompds0)
names(datCompds0)[moduleColors=="purple"]
names(datCompds0)[moduleColors=="turquoise"]


# Stopped here for project all necessary data available

#Appendix 1: IF CUT: Determine cluster under the line
#clust = cutreeStatic(sampleTree, cutHeight = 15, minSize = 10)
#table(clust)
# clust 1 contains the samples we want to keep.
#keepSamples = (clust==1)
#datExpr = datExpr0[keepSamples, ] #generating new compounds file leaving out cut-offs
#nCompds = ncol(datCompds0)
#nSamples = nrow(datCompds0)
#Re-cluster samples
#sampleTree2 = flashClust(dist(datCompds0), method = "average")

#Appendix 2: IF MERGER
# Call an automatic merging function
#merge = mergeCloseModules(datCompds0, dynamicColors, cutHeight = MEDissThres, verbose = 3)
# The merged module colors
mergedColors = merge$colors;
# Eigengenes of the new merged modules:
#mergedMEs = merge$newMEs;
#sizeGrWindow(10, 7)
#pdf(file = "Plots/geneDendro-3.pdf", wi = 9, he = 6)
#plotDendroAndColors(geneTree, cbind(dynamicColors, mergedColors),
#c("Dynamic Tree Cut", "Merged dynamic"),
#dendroLabels = FALSE, hang = 0.03,
#addGuide = TRUE, guideHang = 0.05)
#dev.off()
# Rename to moduleColors
moduleColors = dynamicColors #mergedColors
# Construct numerical labels corresponding to the colors
#colorOrder = c("grey", standardColors(50));
#moduleLabels = match(moduleColors, colorOrder)-1;


getwd()
#rename working directory
workingDir = "/R software/TranscriptomePeach"
#C:\R software\Fruitydata2018
#C:\TOF-DS\FRUITY2017\FRUITY WCNA\R tempate Analysis\FRUITY MARKES INT
# set new working directory:
setwd(workingDir)

datCompds0=read.csv("SqrtTranscript.csv", header = T)
#datCompds0=read.csv("D:/BI3153/NormAreaSqrt.csv", header = T)# equivalent to YData

dim(datCompds0)
names(datCompds0)

traitData = read.csv("ParamTranscript-rep.csv")
#traitData = read.csv("D:/BI3153/xdata2.csv");# equivalent to XData

dim(traitData)
names(traitData)

sampleTree = flashClust(dist(datCompds0), 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.
sizeGrWindow(8,5) #this command is redundant in RStudio
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 = 350, col = "red");

# Usually with scent data there is no cut required as there's no outgroup; if there is go to appendix 1

# Convert traits to a color representation: white means low, red means high, grey means missing entry

traitColors = numbers2colors(traitData ,signed = FALSE);
# Plot the sample dendrogram and the colors underneath.
plotDendroAndColors(sampleTree, traitColors,
groupLabels = names(traitData),
main = "Sample dendrogram and trait heatmap") #if reclustering is required the term 'sampletree' needs to be amended to sampletree2 (see appendix I)

# Choose a set of soft-thresholding powers
powers = c(c(1:10), seq(from = 12, to=20, by=2))
# Call the network topology analysis function
sft = pickSoftThreshold(datCompds0, powerVector = powers, verbose = 5)
# Plot the results:
sizeGrWindow(9, 5)
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("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.90,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("Mean connectivity"))
text(sft$fitIndices[,1], sft$fitIndices[,5], labels=powers, cex=cex1,col="red")
# all above is plot only!
#SELECT SOFTPOWER TO GAIN LARGEST SCALE FREE TOPOLOGY FIT AND LOWWEST MEAN CONNECTIVITY TO CONTINUE

# Calculating adjeacency with soft power 4
softPower =6
adjacency = adjacency(datCompds0, power = softPower)

# Turn adjacency into topological overlap
TOM = TOMsimilarity(adjacency);
dissTOM = 1-TOM

# Call the hierarchical clustering function
geneTree = flashClust(as.dist(dissTOM), method = "average");
# Plot the resulting clustering tree (dendrogram)
sizeGrWindow(8,5)
plot(geneTree, xlab="", sub="", main = "Gene clustering on TOM-based dissimilarity",
labels = FALSE, hang = 0.04)

# We like large modules, so we set the minimum module size relatively high:
minModuleSize =60;
# Module identification using dynamic tree cut:
dynamicMods = cutreeDynamic(dendro = geneTree, distM = dissTOM,
deepSplit = 3, pamRespectsDendro = FALSE,
minClusterSize = minModuleSize);
table(dynamicMods)
# Convert numeric lables into colors
dynamicColors = labels2colors(dynamicMods)
table(dynamicColors)

# Plot the dendrogram and colors underneath
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")



# Calculate eigengenes
MEList = moduleEigengenes(datCompds0, colors = dynamicColors)
MEs = MEList$eigengenes
# Calculate dissimilarity of module eigengenes
MEDiss = 1-cor(MEs);
# Cluster module eigengenes
METree = flashClust(as.dist(MEDiss), method = "average");
# Plot the result
sizeGrWindow(7, 6)
plot(METree, main = "Clustering of module eigengenes",
xlab = "", sub = "")

MEDissThres = 0.01
# Plot the cut line into the dendrogram
abline(h=MEDissThres, col = "red")
###### IF MERGER IS WANTED, E.G. COMBINING MODULES GO TO APPENDIX 2

# Define numbers of genes and samples
nGenes = ncol(datCompds0);
nSamples = nrow(datCompds0);
# Recalculate MEs with color labels
#MEs0 = moduleEigengenes(datCompds0, moduleColors)$eigengenes
#MEs = orderMEs(MEs0)
moduleTraitCor = cor(MEs, traitData, use = "p");
moduleTraitPvalue = corPvalueStudent(moduleTraitCor, nSamples);

# START PLOTTING HEAT MAP MODULE-TRAIT RELATION AND SIGNIFICANCE


sizeGrWindow(10,6)
# Will display correlations and their p-values
textMatrix = paste(signif(moduleTraitCor, 2), "\n(",
signif(moduleTraitPvalue, 1), ")", sep = "");
dim(textMatrix) = dim(moduleTraitCor)
par(mar = c(6, 8.5, 3, 3));
# Display the correlation values within a heatmap plot
labeledHeatmap(Matrix = moduleTraitCor,
xLabels = names(traitData),
yLabels = names(MEs),
ySymbols = names(MEs),
colorLabels = FALSE,
colors = greenWhiteRed(50),
textMatrix = textMatrix,
setStdMargins = FALSE,
cex.text = 0.5,
zlim = c(-1,1),
main = paste("Module-trait relationships"))

#END PLOTTING##

#PREPARING OUTPUT TABLE#

# Define variable 'Sex' containing the 'Sex' column of datTrait
Day = as.data.frame(traitData$Day);
names(Day) = "Day"
Cultivar = as.data.frame(traitData$Cultivar)
names(Cultivar)="Cultivar"
Temperature = as.data.frame(traitData$Temperature)
names(Temperature)="Temperature"
Sample = as.data.frame(traitData$Sample)
names(Sample)="Sample"
# names (colors) of the modules
modNames = substring(names(MEs), 5)
geneModuleMembership = as.data.frame(cor(datCompds0, MEs, use = "p"));
MMPvalue = as.data.frame(corPvalueStudent(as.matrix(geneModuleMembership), nSamples));
names(geneModuleMembership) = paste("MM", modNames, sep="");
names(MMPvalue) = paste("p.MM", modNames, sep="");
geneTraitSignificance = as.data.frame(cor(datCompds0, Day, use = "p"));
GSPvalue = as.data.frame(corPvalueStudent(as.matrix(geneTraitSignificance), nSamples));
names(geneTraitSignificance) = paste("GS.", names(Day), sep="");
names(GSPvalue) = paste("p.GS.", names(Day), sep="");
geneTraitSignificance1 = as.data.frame(cor(datCompds0, Cultivar, use = "p"));
GSPvalue1 = as.data.frame(corPvalueStudent(as.matrix(geneTraitSignificance1), nSamples));
names(geneTraitSignificance1) = paste("GS.", names(Cultivar), sep="");
names(GSPvalue1) = paste("p.GS.", names(Cultivar), sep="");
geneTraitSignificance2 = as.data.frame(cor(datCompds0, Temperature, use = "p"));
GSPvalue2 = as.data.frame(corPvalueStudent(as.matrix(geneTraitSignificance2), nSamples));
names(geneTraitSignificance2) = paste("GS.", names(Temperature), sep="");
names(GSPvalue2) = paste("p.GS.", names(Temperature), sep="");
geneTraitSignificance3 = as.data.frame(cor(datCompds0, Sample, use = "p"));
GSPvalue3 = as.data.frame(corPvalueStudent(as.matrix(geneTraitSignificance3), nSamples));
names(geneTraitSignificance3) = paste("GS.", names(Sample), sep="");
names(GSPvalue3) = paste("p.GS.", names(Sample), sep="");

WCNASag_BT=data.frame(moduleColors=dynamicColors,geneModuleMembership,MMPvalue,geneTraitSignificance,GSPvalue,geneTraitSignificance1,GSPvalue1,geneTraitSignificance2,GSPvalue2,geneTraitSignificance3,GSPvalue3)
write.csv(WCNASag_BT, file="WCNASag_BTnooverload.csv")

#####OUTPUT TABLE CREATED######

# ASSESSMENT CORRELATION MODULE V TRAIT SIGNIFICANCE OF COMPOUNDS WITHIN MODULES
moduleColors = dynamicColors
module = "blue"
column = match(module, modNames);
moduleGenes = moduleColors==module;
sizeGrWindow(7, 7);
par(mfrow = c(1,1));
verboseScatterplot(abs(geneModuleMembership[moduleGenes, column]),
abs(geneTraitSignificance[moduleGenes, 1]),
xlab = paste("Module Membership in", module, "module"),
ylab = "Compound significance for Age",
main = paste("Module membership vs. gene significance\n"),
cex.main = 1.2, cex.lab = 1.2, cex.axis = 1.2, col = module)

names(datCompds0)
names(datCompds0)[moduleColors=="purple"]
names(datCompds0)[moduleColors=="turquoise"]


# Stopped here for project all necessary data available

#Appendix 1: IF CUT: Determine cluster under the line
#clust = cutreeStatic(sampleTree, cutHeight = 15, minSize = 10)
#table(clust)
# clust 1 contains the samples we want to keep.
#keepSamples = (clust==1)
#datExpr = datExpr0[keepSamples, ] #generating new compounds file leaving out cut-offs
#nCompds = ncol(datCompds0)
#nSamples = nrow(datCompds0)
#Re-cluster samples
#sampleTree2 = flashClust(dist(datCompds0), method = "average")

#Appendix 2: IF MERGER
# Call an automatic merging function
#merge = mergeCloseModules(datCompds0, dynamicColors, cutHeight = MEDissThres, verbose = 3)
# The merged module colors
mergedColors = merge$colors;
# Eigengenes of the new merged modules:
#mergedMEs = merge$newMEs;
#sizeGrWindow(10, 7)
#pdf(file = "Plots/geneDendro-3.pdf", wi = 9, he = 6)
#plotDendroAndColors(geneTree, cbind(dynamicColors, mergedColors),
#c("Dynamic Tree Cut", "Merged dynamic"),
#dendroLabels = FALSE, hang = 0.03,
#addGuide = TRUE, guideHang = 0.05)
#dev.off()
# Rename to moduleColors
moduleColors = dynamicColors #mergedColors
# Construct numerical labels corresponding to the colors
#colorOrder = c("grey", standardColors(50));
#moduleLabels = match(moduleColors, colorOrder)-1;
#MEs = mergedMEs;