library(pheatmap)
library(tidyverse)
library(useful)
library(RColorBrewer)

Load files

load("Expression.UT.RData")
load("2.2.c_WGCNA.Modules.UT.RData")
load("2.2.z_WGCNA.Modules.UT.RData")
load("2.3_ModuleEigengene.KWtest.UT.RData")

load("../Enviroment/Environment.RData")
corner(matEnv)

##          BGA-PC 15min Max BGA-PC 15min Mean BGA-PC 15min Min BGA-PC 1d Max
## 201411Is             0.02              0.02             0.02          0.13
## 201411Ky               NA                NA               NA            NA
## 201411Mi             0.08              0.08             0.08          0.08
## 201502Is            -0.05             -0.05            -0.05         -0.01
## 201502Ky               NA                NA               NA            NA
##          BGA-PC 1d Mean
## 201411Is    0.001770833
## 201411Ky             NA
## 201411Mi    0.010416667
## 201502Is   -0.040208333
## 201502Ky             NA

load("../Enviroment/ZooxanthellaeDensity.RData")
matZooDen = matZooDen %>% column_to_rownames("ID")
head(matZooDen)

##             Value
## 201411Is1 2594918
## 201411Is2 5707933
## 201411Is3 2950626
## 201411Is4 4634314
## 201411Is5 2972293
## 201411Mi1 4322314

Corals

# Eigengenes
t = 0.2
eig.c = modules.c[[as.character(t)]]$newMEs

# Average by biological duplicates
samples = apply(meta[rownames(eig.c), c("Month", "Colony")], 1, paste, collapse="")
eig.c = apply(eig.c, 2, tapply, samples, mean)
colnames(eig.c) = str_to_title(gsub("^ME", "", colnames(eig.c)))
corner(eig.c)

##           Midnightblue     Skyblue        Black      Green Floralwhite
## 201411Is1   0.06063190 -0.09174994  0.015707449 0.10931765 -0.04523056
## 201411Ky1  -0.05223023 -0.26343291 -0.034844846 0.09382846 -0.21129775
## 201411Mi1  -0.29148413 -0.28495066  0.042130294 0.07483644 -0.04149472
## 201502Is1  -0.07221333 -0.09794980 -0.010553743 0.06600492 -0.05705480
## 201502Ky1  -0.04773254 -0.07447161 -0.009168224 0.03987089 -0.08074528

# Correlation with environmental paramters
target = "1h Mean"
sites = gsub("[0-9]$", "", rownames(eig.c))
env = data.frame(matEnv[sites, grep(target, colnames(matEnv))], 
                 Zooxanthellae=matZooDen[rownames(eig.c), 1],
                 check.names = F)
cor.c = psych::corr.test(eig.c, env, use="pairwise.complete.obs", adjust="BH")

# Color
colors = rev(colorRampPalette(brewer.pal(11, "RdBu"))(100))

# Visualization
figmat = cor.c$r
figmat[figmat > 0.8] = 0.8; figmat[figmat < -0.8] = -0.8
figmat = figmat[kruskal_test.c$Module, ]
gaps_row = cumsum(table(kruskal_test.c$Category))[1]
gaps_col = ncol(figmat) - 1
#
text = matrix("", nrow=nrow(figmat), ncol=ncol(figmat))
colnames(text) = colnames(figmat)
text[cor.c$p[rownames(figmat), colnames(figmat)] < 0.05] = "*"
text[cor.c$p[rownames(figmat), colnames(figmat)] < 0.01] = "**"
text[cor.c$p[rownames(figmat), colnames(figmat)] < 1e-5] = "***"
pheatmap(figmat, 
         display_numbers = text, 
         fontsize_number = 13, 
         number_color = "white", 
         breaks = seq(-0.8, 0.8, length.out = 100),
         color = colors,
         cluster_rows = F, cluster_cols = F,
         gaps_row = gaps_row, gaps_col = gaps_col,
         filename = paste0("2.4_ModuleEigengene_Environment", gsub(" ","_",target), ".Coral.UT.pdf"),
         fontsize_row = 13, fontsize_col = 13,
         labels_col = gsub(target, "", colnames(figmat)),
         width = 5, height = 8)

Zooxanthellaes

# Eigengenes
t = 0.2
eig.z = modules.z[[as.character(t)]]$newMEs

# Average by biological duplicates
samples = apply(meta[rownames(eig.z), c("Month", "Colony")], 1, paste, collapse="")
eig.z = apply(eig.z, 2, tapply, samples, mean)
colnames(eig.z) = str_to_title(gsub("^ME", "", colnames(eig.z)))
corner(eig.z)

##            Orangered4   Royalblue        Cyan Lightyellow      Purple
## 201411Is1 -0.01405369 -0.09182565 -0.21909679 -0.29301672 -0.13556501
## 201411Ky1 -0.25109303 -0.13616311 -0.25844657 -0.33651410 -0.35182204
## 201411Mi1  0.01056370 -0.09299885 -0.24866809 -0.26394491 -0.09350877
## 201502Is1 -0.14433250  0.04549765  0.03225102  0.02611006 -0.18768293
## 201502Ky1 -0.07394598  0.05860824  0.02189739  0.03246651 -0.11401280

# Correlation with environmental paramters
target = "1h Mean"
sites = gsub("[0-9]$", "", rownames(eig.z))
env = data.frame(matEnv[sites, grep(target, colnames(matEnv))], 
                 Zooxanthellae=matZooDen[rownames(eig.z), 1],
                 check.names = F)
cor.z = psych::corr.test(eig.z, env, use="pairwise.complete.obs", adjust="BH")

# Color
colors = rev(colorRampPalette(brewer.pal(11, "RdBu"))(100))

# Visualization
figmat = cor.z$r
figmat[figmat > 0.8] = 0.8; figmat[figmat < -0.8] = -0.8
figmat = figmat[kruskal_test.z$Module, ]
gaps_row = cumsum(table(kruskal_test.z$Category))[1]
gaps_col = ncol(figmat) - 1
#
text = matrix("", nrow=nrow(figmat), ncol=ncol(figmat))
colnames(text) = colnames(figmat)
text[cor.z$p[rownames(figmat), colnames(figmat)] < 0.05] = "*"
text[cor.z$p[rownames(figmat), colnames(figmat)] < 0.01] = "**"
text[cor.z$p[rownames(figmat), colnames(figmat)] < 1e-5] = "***"
pheatmap(figmat, 
         display_numbers = text, 
         fontsize_number = 13, 
         number_color = "white", 
         breaks = seq(-0.8, 0.8, length.out = 100),
         color = colors,
         cluster_rows = F, cluster_cols = F,
         gaps_row = gaps_row, gaps_col = gaps_col,
         filename = paste0("2.4_ModuleEigengene_Environment", gsub(" ","_",target), ".Zooxanthellae.UT.pdf"),
         fontsize_row = 13, fontsize_col = 13,
         labels_col = gsub(target, "", colnames(figmat)),
         width = 5, height = 8)