Here we go:

Prepare the data

WT_VSD <- read_delim("~/biocond/WT_VSD.csv", delim = ";", escape_double = FALSE, trim_ws = TRUE) # load the expression matrix
datExpr<-t(data.frame(WT_VSD, row.names = 1)) # convert the expression matrix into a data frame and transpose

Calculation of network adjacencies and Topological Overlap Matrix (TOM)

adj = adjacency(datExpr, power = 14, type="signed") # build adjacency matrix signed
distTOM = TOMdist(adj, TOMType = "signed") # calculate dissimilarity TOM signed

Clustering and module identification

geneTree = hclust(as.dist(distTOM), method = "average") # Calculate the gene tree dendrogram
minModuleSize = 100 # Choose the minimum module size
dynamicMods = cutreeDynamic(dendro = geneTree, distM = distTOM,deepSplit = 3, pamRespectsDendro = FALSE,minClusterSize = minModuleSize) # Module identification using dynamic tree cut
Colors = labels2colors(dynamicMods) # Change module labels to colors
plotDendroAndColors(geneTree, Colors, "Dynamic Tree Cut", dendroLabels = FALSE, hang = 0.03, addGuide = TRUE,  guideHang = 0.05, main = "Gene dendrogram and module colors")  # Plots the gene dendrogram

Merge modules sharing similar expression profiles

MEs = moduleEigengenes(datExpr, Colors)$eigengenes # Calculate the module eigengenes
MEDiss = 1-cor(MEs) # Calculate dissimilarity of module eigengenes
METree = hclust(as.dist(MEDiss), method = "average") # Cluster module eigengenes
# Plot the MEs dendrogram and selecte a dissimilarity threshold of 0.25 (0.75 similarity)
MEDissThres = 0.25
plot(METree, main = "Clustering of module eigengenes", xlab = "", sub = "")
abline(h=MEDissThres, col = "red")

merge = mergeCloseModules(datExpr, Colors, cutHeight = MEDissThres, verbose = 3) # Call an automatic merging function

# The merged MEs and module colors
mergedMEs = merge$newMEs
mergedColors = merge$colors

# Plots the gene dendrogram with merged and non-merged modules
plotDendroAndColors(geneTree, cbind(Colors, mergedColors),  c("Unmerged", "Merged"), dendroLabels = FALSE, hang = 0.03, addGuide = TRUE,  guideHang = 0.05, main = "Gene dendrogram and module colors")

Plot an heatmap for each module (mergedColors)

library(gplots)
col = colorpanel(300, 'purple','black','yellow')
colorsA1 = names(table(mergedColors))
pdf("Heatmap.pdf",height=9,width=9)
for (c in 1:length(colorsA1)){
    mergedColors == colorsA1[c]
    heatmap.2(t(datExpr[, mergedColors==colorsA1[c]]), scale = "row", col=col, density.info ="none", trace="none", cexCol=0.7, cexRow=0.7, margin=c(19,11), main = colorsA1[c], Colv = FALSE)
}; 
dev.off()

If you check the pdf file with the heatmaps of each module you will definitely see that almost every module is strongly associated with your experimental traits. Here is some example:

1. In the brown module all genes are expressed at very high level in N3P and WTP samples but downregulated in SM samples
2. In the royalblue all genes are expressed at very high level in SM samples but downregulated in N3P and WTP samples
3. In the cyan all genes are up-regulated in CA1 samples
4. In the darkred all genes are up-regulated in DG samples

Go ahead with the module-trait relationship analysis and remember that the final modules and the module eigengenes are in mergedColors and mergedMEs, respectively

Thank you for the datExpr

Also, when I checked soft Thresholds for "signed" network, it is now different than the default, the power 12 or 14 seems to be adequate

This is ok. The soft Thresholds for signed networks almost double.