#######################Load your Packages########################
#################################################################
library("WGCNA")
library("colorspace")
library(pbkrtest)
library("rgl")
library("car")
library("plot3D")
library("pca3d")
library("rgl")
library("sva")
library("Biobase")
library("flashClust")
library("Rcpp")
library("ggplot2")
library("FactoMineR") #not available for 3.2.1
library("factoextra")
library("pca3d")
library(scatterplot3d)
library(ggplot2)
library(grid)
library(splitstackshape)
library(scatterplot3d)
library(factoextra)
library(gplots) 
library(RColorBrewer)

#################### A: Load your Inputs ########################
#################################################################
rawData=read.csv("WT_aCSF_v_NMDA_MFI.csv")
head(rawData)
datExpr=as.data.frame((rawData[,-c(1:3,12)]))
row.names(datExpr)=rawData$IP_Probe
head(datExpr)

########### B: Eliminate MFI values less than 100 (noise) #######
#################################################################
aves=rowMeans(datExpr)
datExpr["new.col"]=aves
datlog=subset(datExpr, new.col>100)
nrow(datlog)                                              
head(datlog)
ncol(datlog)->cols
datlog=datlog[,-cols]
exprs = as.matrix(datlog,header=TRUE, sep="\t", row.names=5,as.is=TRUE)
class(exprs)
dim(exprs)
colnames(exprs)
traits=read.csv("traits.csv", row.names =1)
rownames(traits)                     
all(rownames(traits)==colnames(exprs))                            # doublecheck that this is true, if false your inputs files disagree
colnames(exprs)
names(traits)

#################### C: Run COMBAT #############################
################################################################

                                                                  #ensure that row.names matches names(traits), label description can be whatever
metadata=data.frame(labelDescription=c( "experiment", "treatment" ),
                      row.names=c("experiment", "treatment"))
phenoData=new("AnnotatedDataFrame",
               data=traits,varMetadata=metadata)
all(names(traits)==metadata$row.names)                           #double check rownames match
dataS = ExpressionSet(assayData=exprs,phenoData=phenoData)      
pheno = pData(dataS)                                             
edata = exprs(dataS)
batch = pheno$experiment                                         #match to traits
modcombat = model.matrix(~1, data=pheno)
combat_edata = ComBat(dat=edata, batch=batch, mod=modcombat, 
                      par.prior=TRUE, prior.plots=TRUE)
log2data = log(combat_edata,2)                                   # assess plots and decide whether to use para vs. non para (par.prior=FALSE) combat
log2data2 = log(edata,2)                                         # Reference: Johnson, WE, Rabinovic, A, and Li, C (2007). 
tdata = t(log2data)                                              # Adjusting batch effects in microarray expression data using Empirical Bayes methods.
tdata2 = t(log2data2)                                            # Biostatistics 8(1):118-127
plot(hclust(dist(tdata)), main="postCombat")
plot(hclust(dist(tdata2)),main="preCombat")
write.csv(combat_edata,"combated_data.csv")

#################### D: Reformat Combatted data #################
#################################################################

rawData=read.csv("combated_data.csv")
rawData2=as.data.frame(rawData)
head(rawData2)
datExpr=as.data.frame(t(rawData2[,-c(1)]))
datExpr=log(datExpr,2)
names(datExpr)=rawData2$X
head(datExpr)
dim(datExpr)
nGenes=ncol(datExpr)
nSamples=nrow(datExpr)
datTraits<-traits
sampleTree2 = hclust(dist(datExpr))
traitColors = numbers2colors(datTraits, signed = FALSE)
plotDendroAndColors(sampleTree2, traitColors, 
                    groupLabels = names(datTraits),
                    main = "Sample dendrogram and trait heatmap")

################E: PCA on post-combat data ######################
#################################################################
dev.off()
PCA(datExpr)->expPCA
labels=as.factor(c(1,1,1,1,2,2,2,2))
fviz_pca_ind(expPCA, axes = c(1,2),
             addEllipses = TRUE, 
             habillage = labels, 
             ellipse.level=.6)
dimdesc(expPCA, axes=c(1))->PCAdrivers1                               #gives a correlation between PCs and PiSCES
PCAdrivers1

########################## F: Run CNA ###########################
#################################################################

powers = c(c(1:10), seq(from = 12, to=20, by=2))
sft = pickSoftThreshold(datExpr, 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")

softPower= 6                                                         #make sure to change this number
minModuleSize= 10
adjacency=adjacency(datExpr,power = softPower)
TOM=TOMsimilarity(adjacency)
dissTOM = 1-TOM
geneTree = flashClust(as.dist(dissTOM), method = "average")
dev.off()
plot(geneTree)
dynamicMods = cutreeDynamic(dendro = geneTree, distM = dissTOM, deepSplit = 2, pamRespectsDendro = FALSE, minClusterSize = minModuleSize)
dynamicColors = labels2colors(dynamicMods)
MEList = moduleEigengenes(datExpr, colors = dynamicColors)
MEs = MEList$eigengenes
MEDiss = 1-cor(MEs)
METree = flashClust(as.dist(MEDiss), method = "average")
MEDissThres = 0.25
merge = mergeCloseModules(datExpr, dynamicColors, cutHeight = MEDissThres, verbose = 3)
mergedColors = merge$colors
mergedMEs = merge$newMEs
PPInames = colnames(datExpr)
sizeGrWindow(120,4)
plotDendroAndColors(geneTree, cbind(dynamicColors, mergedColors), c("Dynamic Tree Cut", "Merged dynamic"), cex.dendroLabels = 0.7, dendroLabels = PPInames, addGuide = TRUE, guideHang = 0.05)

moduleColors = mergedColors
colorOrder = c("grey", standardColors(50));
moduleLabels = match(moduleColors, colorOrder)-1;
MEs = mergedMEs;

################################################################
#### G:Check the TOM plot and make sure CNA did a good job #####


plotTOM = dissTOM^3
diag(plotTOM) = NA
TOMplot(plotTOM, geneTree, moduleColors, main = "Network heatmap plot, all genes")          


#pdf(file = paste0("PPI",softPower,minModuleSize,".pdf"),wi = 36, he = 9)
#plotDendroAndColors(geneTree, cbind(dynamicColors, mergedColors), c("Dynamic Tree Cut", "Merged dynamic"), dendroLabels = PPInames, addGuide = TRUE, guideHang = 0.05)
#dev.off()

#################################################################
########### H: Generate module-trait relationships ##############
dev.off() 
moduleTraitCor = cor(MEs, datTraits, use = "p");
moduleTraitPvalue = corPvalueStudent(moduleTraitCor, nSamples);
sizeGrWindow(10,6)
textMatrix = paste(signif(moduleTraitCor, 2), "\n(",
                   signif(moduleTraitPvalue, 1), ")", sep = "");
dim(textMatrix) = dim(moduleTraitCor)
par(mar = c(6, 8.5, 3, 3));
palette<- colorRampPalette(c("blue", "white", "white","red"))(n=299)
labeledHeatmap(Matrix = moduleTraitCor,
               xLabels = names(datTraits),
               yLabels = names(MEs),
               ySymbols = names(MEs),
               colorLabels = FALSE,
               colors = palette,
               textMatrix = textMatrix,
               setStdMargins = FALSE,
               cex.text = 0.5,
               zlim = c(-1,1),
               main = paste("Module-trait relationships"))

#################################################################
###### I: Define variables for output file and heatmaps ######### 
analTrait = "magenta"                                          # define the color you're interested in 
analCondition = "treatment"                                  # and the datTrait condition you're interested in
weight = as.data.frame(datTraits$treatment)                  # and make sure analCondition matches the datTrait column here
names(weight) = "treatment"                                  #...and here
aggregate(datExpr, list(datTraits$treatment), mean) -> means # ... and here
meanNames=c("NMDA", "ACSF")                                  # Finally, name the same variable for the output file, e.g. "Mutant")

###################### and run it ###############################
module = analTrait
rownames(means)=meanNames
means=means[,-c(1)]
foldChange= means[1,]/means[2,]
means=2^means
rownames(foldChange)="FoldChange"
modNames = substring(names(MEs), 3)
geneModuleMembership = as.data.frame(cor(datExpr, 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(datExpr, weight, use = "p"))
GSPvalue = as.data.frame(corPvalueStudent(as.matrix(geneTraitSignificance), nSamples))
names(geneTraitSignificance) = paste("GS.", names(weight), sep="")
names(GSPvalue) = paste("p.GS.", names(weight), sep="")
column = match(module, modNames)
moduleGenes = moduleColors==module
probes = colnames(datExpr)
geneInfo0 = data.frame(probes = probes, moduleColor = moduleColors, t(means), geneTraitSignificance, GSPvalue)
modOrder = order(-abs(cor(MEs, weight, use = "p")))
for (mod in 1:ncol(geneModuleMembership))
{
  oldNames = names(geneInfo0)
  geneInfo0 = data.frame(geneInfo0, geneModuleMembership[, modOrder[mod]],
                         MMPvalue[, modOrder[mod]]);
  names(geneInfo0) = c(oldNames, paste("MM.", modNames[modOrder[mod]], sep=""),
                       paste("p.MM.", modNames[modOrder[mod]], sep=""))
}
moi1="p.MM"
geneorderspacer=paste(moi1,analTrait, sep=".")
geneOrder = order(abs(geneInfo0[,geneorderspacer])) 
geneInfo = geneInfo0[geneOrder, ]


#################################################################
############### J: name the output CNA file and save ############

write.csv(geneInfo, file = "CNA_outputs.csv")

#################################################################
################# K: Make module heatmaps #######################

MOI="black"                                                       #define the Module Of Interest (MOI)

########################run it###################################
my_palette <- colorRampPalette(c("blue","white","red"))(n = 299)
colors = c(seq(-1,-.11, length=145),                                
           seq(-.1,.1,length=10),
           seq(0.11,1,length=145))
datHeatmap=data.frame(t(datExpr),MMPvalue,t(foldChange),geneTraitSignificance)
moiname1=paste(moi1,MOI, sep=".")
moiname2=paste(moi1,MOI, sep="")
mapOrder = order(abs(geneInfo0[,moiname1]))     
datHeatM = datHeatmap[mapOrder, ]
sortedHeat=subset(datHeatM,datHeatM[,moiname2]<.05)
as.matrix(sortedHeat[,c(1:8)])->aa                 
hclust_fun= function(x) hclust(x, method = "ward")

heatmap.2 (aa,
           col= my_palette,
           dendrogram = "both",
           scale = "row", 
           key = T, 
           keysize = 1,
           density.info = "none",
           trace = "none",
           hclustfun = hclust_fun,
           cexCol = .8,
           cexRow = .5,
           Colv = TRUE,
           Rowv = TRUE,
           margin = c(5,15),
           srtCol = 45,
           sepcolor = "black",
           sepwidth= c(.01, .01),
           key.xlab = "Relative Scaled Value",
           key.title = NA,
           main = MOI)

################## L: To create ANCandCNA heatmaps ##############
#################################################################
write.csv(aa, "aa.csv")                                          #manually edit before the next step
read.csv("aa.csv", row.names = 1)->bb
as.matrix(bb)->bb
heatmap.2 (bb,
           col= my_palette,
           dendrogram = "col",
           scale = "row", 
           key = T, 
           keysize = 1,
           density.info = "none",
           trace = "none",
           hclustfun = hclust_fun,
           cexCol = .8,
           cexRow = .5,
           Colv = TRUE,
           Rowv = FALSE,
           margin = c(5,15),
           srtCol = 45,
           sepcolor = "black",
           sepwidth= c(.01, .01),
           key.xlab = "Relative Scaled Value",
           key.title = NA,
           main = "ANCnCNA Hits")

#################################################################
#########################END OF PROGRAM##########################

