---
title: "WGCNA in female APOE3 animals"
Date: April 2020
author: Amin Haghani
output: html_document
---

```{r packages, message=FALSE}
library(WGCNA)
```


I used the transformed data by Voom (log2 expression) for WGCNA
```{r data preparation}
datExpr <- t(v$E)
dim(datExpr)
```

```{r Convert trait character variables to dicotomous variables (0 vs 1)}
datTraits <- as.factor(x$samples$Treatment)
levels(datTraits) <- c(0,1)
datTraits
```

```{r picking soft threshold}
powers = c(c(1:10), seq(from = 10, to=160, by=2))
sft = pickSoftThreshold(datExpr, powerVector = powers, verbose = 5, networkType = "signed")
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")
sft$fitIndices

```

```{r Calculate adjacency, Topological Overlap Matrix, and Dissimilarity Matrix}
softPower = 96

adjacency = adjacency(datExpr, power = softPower, type="signed");
TOM = TOMsimilarity(adjacency, TOMType = "signed");
dissTOM = 1-TOM
```

```{r Cluster using HC}
geneTree = hclust(as.dist(dissTOM), method = "average")
```

```{r Plot Dendrogram}
sizeGrWindow(12,9)
plot(geneTree, xlab="", sub="", main = "Gene clustering on TOM-based dissimilarity",
     labels = FALSE, hang = 0.04)
```
```{r Generate Modules Based on Dendrogram}
minModuleSize = 30;
dynamicMods = cutreeDynamic(dendro = geneTree, distM = dissTOM,
                            deepSplit = 2, pamRespectsDendro = FALSE,
                            minClusterSize = minModuleSize)
```

```{r Convert numbers to color labels}
dynamicColors = labels2colors(dynamicMods)
table(dynamicColors)
```

```{r Relate dynamic modules to Traits}
nGenes=ncol(datExpr)
nSamples=nrow(datExpr)
MEsDynamic=moduleEigengenes(datExpr,dynamicColors)$eigengenes
MEsDynamic=orderMEs(MEsDynamic)

plotEigengeneNetworks(MEsDynamic, "", marDendro = c(0,4,1,2), marHeatmap = c(3,4,1,2))

moduleFreqDynamic <- as.data.frame(table(dynamicColors))
```


```{r Plot the dendrogram with modules depicted 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")
```


```{r}
# Test for module overlap and merge collinear modules (correlation of r>=0.85)
MEList = moduleEigengenes(datExpr, 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.15
abline(h=MEDissThres, col = "red")
merge = mergeCloseModules(datExpr, dynamicColors, cutHeight = MEDissThres, verbose = 3)
mergedColors = merge$colors;
mergedMEs = merge$newMEs;


#dev.off()
sizeGrWindow(12, 9)
#pdf(file = "Plots/geneDendro-3.pdf", wi = 9, he = 6)
plotDendroAndColors(geneTree, cbind(dynamicColors, mergedColors),
                    c("Dynamic colors","Merged modules"),
                    dendroLabels = FALSE, hang = 0.03,
                    addGuide = TRUE, guideHang = 0.05)
#dev.off()
```

```{r plot with nPM}
nPM <- datTraits
nPM <- as.numeric(cor(datExpr,nPM,use="p"))
nPM.Color=numbers2colors(nPM,signed=T)

plotDendroAndColors(geneTree, cbind(dynamicColors, mergedColors, nPM.Color),
                    c("Dynamic colors","Merged modules", "nPM correlation"),
                    dendroLabels = FALSE, hang = 0.03,
                    addGuide = TRUE, guideHang = 0.05)
```


```{r Rename using merged modules}
moduleColors = mergedColors
colorOrder = c("grey", standardColors(50));
moduleLabels = match(moduleColors, colorOrder)-1;
```

```{r calculate intamodule connectivity}
colorh1 <- moduleColors
datME=moduleEigengenes(datExpr,colorh1, softPower = 92)$eigengenes

#Calculating KME
datKME=signedKME(datExpr, datME, outputColumnName="MM.")

datKME$modules <- colorh1


# Display the first few rows of the data frame

datKME <- cbind(genes, datKME)
head(datKME)
#write.csv(datKME, file = "./Results/Female, E3, nPM/Modules KME.csv")
```


```{r Relate to Traits}
table(moduleColors)
nGenes=ncol(datExpr)
nSamples=nrow(datExpr)
MEs0=moduleEigengenes(datExpr,moduleColors)$eigengenes
MEs=orderMEs(MEs0)

plotEigengeneNetworks(MEs, "", marDendro = c(0,4,1,2), marHeatmap = c(3,4,1,2))


moduleFreq <- as.data.frame(table(moduleColors))
```



```{r Relate traits to MEs}
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));


association <- data.frame(correlation = moduleTraitCor, pval = moduleTraitPvalue)
association$qval <- qvalue(association$pval)$qvalues

association$moduleColors <- substring(rownames(association), 3)

association <- merge(association, moduleFreq, by="moduleColors", all=TRUE)


#write.csv(association, "./Results/Female, E3, nPM/module association.csv")
```






















