lnames=load(file="ES-networkConstruction-stepByStep.RData")
ESMEs=MEs
ESmoduleLabels=moduleLabels
ESmoduleColors=moduleColors
ESgeneTree=geneTree

lnames=load(file="RS-networkConstruction-stepByStep.RData")
RSMEs=MEs
RSSmoduleLabels=moduleLabels
RSmoduleColors=moduleColors
RSgeneTree=geneTree

lnames = load(file = "ES.RData")
datExprT=datExpr0
lnames = load(file = "RS.RData")
datExprN=datExpr0

setLabels=c("ES","RS")
colorES=ESmoduleColors
colorRS=RSmoduleColors

#module level comparison
multiExpr = list(ES=list(data=datExprT),RS=list(data=datExprN))
multiColor = list(ES=colorES,RS=colorRS)
system.time({mp=modulePreservation(multiExpr,multiColor,referenceNetworks=1,nPermutation=100,
                                   randomSeed=1,quickCor=0,verbose=3)})

ref = 1
test = 2
statsObs = cbind(mp$quality$observed[[ref]][[test]][, -1], mp$preservation$observed[[ref]][[test]][, -1])
statsZ = cbind(mp$quality$Z[[ref]][[test]][, -1], mp$preservation$Z[[ref]][[test]][, -1])
print( cbind(statsObs[, c("medianRank.pres", "medianRank.qual")],signif(statsZ[, c("Zsummary.pres", "Zsummary.qual")], 2)) )
statResult <- cbind(statsObs[, c("medianRank.pres", "medianRank.qual")],signif(statsZ[, c("Zsummary.pres", "Zsummary.qual")], 2))
write.csv(statResult,"statResult.ES-RS.csv")

# Module labels and module sizes are also contained in the results
modColors = rownames(mp$preservation$observed[[ref]][[test]])
moduleSizes = mp$preservation$Z[[ref]][[test]][, 1]
# leave grey and gold modules out
plotMods = !(modColors %in% c("grey", "gold"))
# Text labels for points
text = modColors[plotMods]
# Auxiliary convenience variable
plotData = cbind(mp$preservation$observed[[ref]][[test]][, 2], mp$preservation$Z[[ref]][[test]][, 2])
# Main titles for the plot
mains = c("Preservation Median rank", "Preservation Zsummary")
# Start the plot
sizeGrWindow(35,18)
par(mfrow = c(1,2))
par(mar = c(4,2.5,2.5,0.5))
for (p in 1:2)
{
  min = min(plotData[, p], na.rm = TRUE);
  max = max(plotData[, p], na.rm = TRUE);
  # Adjust ploting ranges appropriately
  if (p==2)
  {
    if (min > -max/10) min = -max/10
    ylim = c(min - 0.1 * (max-min), max + 0.1 * (max-min))
  } else
    ylim = c(max + 0.1 * (max-min), min - 0.1 * (max-min))
  plot(moduleSizes[plotMods], plotData[plotMods, p], col = 1, bg = modColors[plotMods], pch = 21,
       main = mains[p],
       cex = 2.2,
       ylab = " ", xlab = "Module size", log = "x",
       cex.lab = 1, cex.axis = 1, cex.main =1.5)
  labelPoints(moduleSizes[plotMods], plotData[plotMods, p], text, cex = 0.7, offs = 0.06,protectEdges=TRUE);
  # For Zsummary, add threshold lines
  if (p==2) {
    abline(h=0)
    abline(h=2, col = "red", lty = 2)
    abline(h=10, col = "red", lty = 2)
  } }