Identificaiton of Time- & Colony-dependent modulesΒΆ
library(pheatmap)
library(RColorBrewer)
library(tidyverse)
library(openxlsx)Load files
t = 0.2
load("Expression.UT.RData")
load("2.1.c_WGCNA_ThresholdingPower.UT.RData")
load("2.1.z_WGCNA_ThresholdingPower.UT.RData")
load("2.2.c_WGCNA.Modules.UT.RData")
load("2.2.z_WGCNA.Modules.UT.RData")
load("../Annotation/coral/Annotation.Coral.UT.RData")
load("../Annotation/zooxanthellae/rnaSPAdes.1Mreads.UT/Annotation.Zooxanthellae.UT.RData")Corals
Proportion of annotated genes in the modules
modules = modules.c[[as.character(t)]]$colors
annot = data.frame(Gene=colnames(mat.c), Module=modules, stringsAsFactors=F) %>%
mutate(Swissprot = Gene %in% gene2GO.c$Gene,
KAAS = Gene %in% gene2KO.c$Gene) %>%
group_by(Module) %>%
summarize(Size = length(Swissprot),
Annotated.Swissprot = sum(Swissprot) / length(Swissprot) * 100,
Annotated.KAAS = sum(KAAS) / length(KAAS) * 100) %>%
mutate(Module = Hmisc::capitalize(Module))## `summarise()` ungrouping output (override with `.groups` argument)head(annot)## # A tibble: 6 x 4
## Module Size Annotated.Swissprot Annotated.KAAS
## <chr> <int> <dbl> <dbl>
## 1 Bisque4 55 52.7 25.5
## 2 Black 754 66.2 43.2
## 3 Blue 1735 62.4 42.4
## 4 Brown4 58 79.3 55.2
## 5 Cyan 268 62.3 42.2
## 6 Darkgrey 173 52.6 27.7Classification of modules
rename = function(mat){
x2 = rownames(mat)
x2[rownames(mat) == "201602_Is1"] = "201602_Mi1"
x2[rownames(mat) == "201602_Mi1"] = "201602_Is1"
rownames(mat) = x2
mat = mat[sort(rownames(mat)), ]
mat
}eigen = modules.c[[as.character(t)]]$newMEs
eigen = rename(eigen)
# Separate module by expression pattern
kruskal_test.c = matrix(nrow=0, ncol=5)
colnames(kruskal_test.c) = c("Module", "pvalue.colony", "pvalue.time", "pvalue.site", "Category")
for(i in colnames(eigen)){
tmp = data.frame(value=eigen[,i], meta[rownames(eigen),])
pvalues = c(kruskal.test(value~factor(Colony), tmp)$p.value,
kruskal.test(value~factor(Month), tmp)$p.value,
kruskal.test(value~factor(Site), tmp)$p.value)
min_p = min(pvalues)
if(min_p >= 0.01) category = ""
else category = c("Colony", "Time", "Site")[pvalues == min_p]
kruskal_test.c = rbind(kruskal_test.c, c(i, pvalues, category))
}
kruskal_test.c = data.frame(kruskal_test.c, stringsAsFactors=F)
for(i in 2:4) kruskal_test.c[,i] = as.numeric(kruskal_test.c[,i])
kruskal_test.c = kruskal_test.c %>% mutate(Module=str_to_title(gsub("ME", "", Module)))
write.xlsx(kruskal_test.c, file="2.3_ModuleEigengene.Coral.Kruskal-Wallis.UT.xlsx")
kruskal_test.c## Module pvalue.colony pvalue.time pvalue.site Category
## 1 Midnightblue 5.573652e-08 5.418771e-02 1.636627e-03 Colony
## 2 Skyblue 2.241974e-03 7.709554e-07 1.293630e-03 Time
## 3 Black 3.015374e-01 3.046807e-06 1.857722e-01 Time
## 4 Green 1.383875e-01 1.430471e-06 2.267926e-01 Time
## 5 Floralwhite 6.707788e-03 1.860948e-04 1.157967e-02 Time
## 6 Brown4 5.739024e-01 4.936269e-09 1.352950e-01 Time
## 7 Cyan 7.398517e-03 1.591865e-08 1.652944e-03 Time
## 8 Darkorange2 3.091725e-01 5.073069e-09 1.690469e-01 Time
## 9 Paleturquoise 4.508452e-01 1.444891e-06 9.216888e-01 Time
## 10 Skyblue3 8.894206e-02 8.360159e-09 2.147567e-01 Time
## 11 Darkturquoise 3.270258e-08 9.549869e-01 1.991831e-01 Colony
## 12 Mediumpurple3 1.901059e-09 3.026667e-01 6.500297e-02 Colony
## 13 Bisque4 6.575650e-02 2.120217e-03 3.839335e-02 Time
## 14 Lightgreen 7.194623e-07 1.019709e-01 6.692426e-04 Colony
## 15 Darkmagenta 3.854757e-05 1.406723e-04 2.261562e-04 Colony
## 16 Steelblue 1.685718e-08 6.913375e-01 3.315714e-06 Colony
## 17 Grey60 5.224896e-10 4.316311e-01 8.772746e-04 Colony
## 18 Lightcyan 1.904185e-09 1.277668e-01 6.042484e-01 Colony
## 19 White 3.444674e-03 1.976662e-04 1.128058e-02 Time
## 20 Yellowgreen 3.815667e-05 7.167705e-02 2.435243e-01 Colony
## 21 Lightyellow 5.466979e-01 3.941902e-08 7.066833e-01 Time
## 22 Lightcyan1 3.320162e-01 4.276202e-08 1.796819e-01 Time
## 23 Turquoise 2.169349e-02 6.771278e-06 5.998884e-02 Time
## 24 Greenyellow 4.904319e-04 1.966914e-07 4.106867e-03 Time
## 25 Tan 5.171021e-03 2.762848e-08 1.673173e-02 Time
## 26 Darkred 1.964268e-08 5.912912e-01 5.740096e-01 Colony
## 27 Saddlebrown 2.780873e-08 2.851287e-01 3.952965e-02 Colony
## 28 Darkorange 1.510346e-07 9.215046e-01 1.527360e-06 Colony
## 29 Darkgrey 2.415698e-05 1.066430e-03 5.217150e-01 Colony
## 30 Ivory 3.852822e-09 5.030189e-01 4.804959e-01 Colony
## 31 Purple 8.427431e-05 1.771514e-01 2.111140e-01 Colony
## 32 Royalblue 8.005544e-08 4.068414e-01 3.793982e-01 Colony
## 33 Orange 5.828074e-01 3.136213e-08 1.982500e-01 Time
## 34 Orangered4 6.109914e-03 1.466537e-05 9.828483e-01 Time
## 35 Blue 5.676562e-01 2.866131e-07 2.752429e-01 Time
## 36 Sienna3 2.504945e-01 2.559881e-07 4.319452e-01 TimeVisualization
# Average by biological duplicates
samples = apply(meta[rownames(eigen), c("Month", "Colony")], 1, paste, collapse="_")
figmat = apply(eigen, 2, tapply, samples, mean)
colnames(figmat) = str_to_title(gsub("^ME", "", colnames(figmat)))# Time-dependent modules
figmat.Time = t(figmat[, kruskal_test.c$Module[kruskal_test.c$Category == "Time"]])
modules = rownames(figmat.Time)
figmat.Time = figmat.Time %>% as.data.frame() %>% rownames_to_column("Module") %>%
gather(key="Sample", value="value", -Module) %>%
mutate(Module=factor(Module, levels=modules), Time=meta[Sample,"Month"], Colony=meta[Sample,"Colony"]) %>%
mutate(value=ifelse(value > 0.3, 0.3, ifelse(value < -0.3, -0.3, value)))
figs = list(); size = c(); times = unique(figmat.Time$Time)
for(i in times){
size = c(size, length(unique(figmat.Time$Colony[figmat.Time$Time == i])))
figs[[i]] = ggplot(subset(figmat.Time, Time==i)) +
geom_tile(aes(y=Module, x=Colony, fill=value)) +
scale_fill_gradientn(colors=rev(brewer.pal(9, "RdBu")), limits=c(-0.3, 0.3)) +
theme_minimal() +
theme(legend.position="none", plot.margin = unit(c(0, 0, 0, 0), "cm"),
axis.text.x=element_text(angle=90, vjust=0.5, hjust=1)) +
xlab(NULL) + ylab(NULL)
if(i != times[1]) figs[[i]] = figs[[i]] + theme(axis.text.y=element_blank())
}
size[1] = size[1] + 8
cowplot::plot_grid(plotlist = figs, nrow=1, align="h", rel_widths=size) +
ggsave("2.3_ModuleEigengene.Coral.Time.UT.pdf", width=6, height=2.8)
legend = cowplot::get_legend(figs[[1]] + theme(legend.position="right", legend.title=element_blank()))
cowplot::plot_grid(legend) + ggsave("2.3_ModuleEigengene.Coral.Time.UT.Legend.pdf", width=1, height=2.8)
figs = list()
annot.Time = subset(annot, Module %in% modules)
annot.Time$Module = factor(annot.Time$Module, levels=modules)
figs[["size"]] = ggplot(annot.Time) +
geom_bar(aes(y=Size, x=Module, fill=I(tolower(Module))), stat="identity", color="grey60") +
theme_minimal() + coord_flip() + xlab(NULL) + ylab("Module size") +
scale_y_log10(limits = c(1, 5000)) +
theme(axis.title.x=element_text(size=8))
size = c(18)
figs[["annot1"]] = ggplot(annot.Time) +
geom_bar(aes(y=Annotated.Swissprot, x=Module, fill=I(tolower(Module))), stat="identity", color="grey60") +
theme_minimal() + coord_flip() + xlab(NULL) + ylab("Swissprot [%]") +
scale_y_continuous(limits=c(0, 80)) +
theme(axis.text.y = element_blank(), axis.title.x=element_text(size=8))
size = c(size, 10)
figs[["annot2"]] = ggplot(annot.Time) +
geom_bar(aes(y=Annotated.KAAS, x=Module, fill=I(tolower(Module))), stat="identity", color="grey60") +
theme_minimal() + coord_flip() + xlab(NULL) + ylab("KAAS [%]") +
scale_y_continuous(limits=c(0, 80)) +
theme(axis.text.y = element_blank(), axis.title.x=element_text(size=8))
size = c(size, 10)
cowplot::plot_grid(plotlist = figs, nrow=1, align="h", rel_widths=size) +
ggsave("2.3_ModuleInfo.Coral.Time.UT.pdf", width=4, height=2.8)
figmat.Colony = t(figmat[, kruskal_test.c$Module[kruskal_test.c$Category == "Colony"]])
modules = rownames(figmat.Colony)
figmat.Colony = figmat.Colony %>% as.data.frame() %>% rownames_to_column("Module") %>%
gather(key="Sample", value="value", -Module) %>%
mutate(Module=factor(Module, levels=modules), Time=meta[Sample,"Month"], Colony=meta[Sample,"Colony"]) %>%
mutate(value=ifelse(value > 0.3, 0.3, ifelse(value < -0.3, -0.3, value)))
figs = list(); size = c(); colonies = sort(unique(figmat.Colony$Colony))
for(i in colonies){
size = c(size, length(unique(figmat.Colony$Sample[figmat.Colony$Colony == i])))
figs[[i]] = ggplot(subset(figmat.Colony, Colony==i)) +
geom_tile(aes(y=Module, x=Time, fill=value)) +
scale_fill_gradientn(colors=rev(brewer.pal(9, "RdBu")), limits=c(-0.3, 0.3)) +
theme_minimal() +
theme(legend.position="none", plot.margin = unit(c(0, 0, 0, 0), "cm"),
axis.text.x=element_text(angle=90, vjust=0.5, hjust=1)) +
xlab(NULL) + ylab(NULL)
if(i != colonies[1]) figs[[i]] = figs[[i]] + theme(axis.text.y=element_blank())
}
size[1] = size[1] + 8
cowplot::plot_grid(plotlist = figs, nrow=1, align="h", rel_widths=size) +
ggsave("2.3_ModuleEigengene.Coral.Colony.UT.pdf", width=6, height=2.5)
figs = list()
# Proportion of annotated genes
annot.Colony = subset(annot, Module %in% modules)
annot.Colony$Module = factor(annot.Colony$Module, levels=modules)
figs[["size"]] = ggplot(annot.Colony) +
geom_bar(aes(y=Size, x=Module, fill=I(tolower(Module))), stat="identity", color="grey60") +
theme_minimal() + coord_flip() + xlab(NULL) + ylab("Module size") +
scale_y_log10(limits = c(1, 5000)) +
theme(axis.title.x=element_text(size=8))
size = c(18)
figs[["annot1"]] = ggplot(annot.Colony) +
geom_bar(aes(y=Annotated.Swissprot, x=Module, fill=I(tolower(Module))), stat="identity", color="grey60") +
theme_minimal() + coord_flip() + xlab(NULL) + ylab("Swissprot [%]") +
scale_y_continuous(limits=c(0, 80)) +
theme(axis.text.y = element_blank(), axis.title.x=element_text(size=8))
size = c(size, 10)
figs[["annot2"]] = ggplot(annot.Colony) +
geom_bar(aes(y=Annotated.KAAS, x=Module, fill=I(tolower(Module))), stat="identity", color="grey60") +
theme_minimal() + coord_flip() + xlab(NULL) + ylab("KAAS [%]") +
scale_y_continuous(limits=c(0, 80)) +
theme(axis.text.y = element_blank(), axis.title.x=element_text(size=8))
size = c(size, 10)
cowplot::plot_grid(plotlist = figs, nrow=1, align="h", rel_widths=size) +
ggsave("2.3_ModuleInfo.Coral.Colony.UT.pdf", width=4, height=2.4)
Zooxanthellae
Proportion of annotated genes in the modules
modules = modules.z[[as.character(t)]]$colors
annot = data.frame(Gene=colnames(mat.z), Module=modules, stringsAsFactors=F) %>%
mutate(Swissprot = Gene %in% gene2GO.z$Gene,
KAAS = Gene %in% gene2KO.z$Gene) %>%
group_by(Module) %>%
summarize(Size = length(Swissprot),
Annotated.Swissprot = sum(Swissprot) / length(Swissprot) * 100,
Annotated.KAAS = sum(KAAS) / length(KAAS) * 100) %>%
mutate(Module = Hmisc::capitalize(Module))## `summarise()` ungrouping output (override with `.groups` argument)head(annot)## # A tibble: 6 x 4
## Module Size Annotated.Swissprot Annotated.KAAS
## <chr> <int> <dbl> <dbl>
## 1 Bisque4 78 34.6 24.4
## 2 Brown4 297 38.7 30.0
## 3 Cyan 3986 36.9 26.2
## 4 Darkorange 239 45.6 34.3
## 5 Darkslateblue 73 45.2 26.0
## 6 Darkturquoise 318 30.5 20.8Classification of modules
eigen = modules.z[[as.character(t)]]$newMEs
eigen = eigen[,colnames(eigen) != "MEgrey"]
eigen = rename(eigen)
# Separate module by expression pattern
kruskal_test.z = matrix(nrow=0, ncol=5)
colnames(kruskal_test.z) = c("Module", "pvalue.colony", "pvalue.time", "pvalue.site", "Category")
for(i in colnames(eigen)){
tmp = data.frame(value=eigen[,i], meta[rownames(eigen),])
pvalues = c(kruskal.test(value~factor(Colony), tmp)$p.value,
kruskal.test(value~factor(Month), tmp)$p.value,
kruskal.test(value~factor(Site), tmp)$p.value)
min_p = min(pvalues)
if(min_p >= 0.01) category = ""
else category = c("Colony", "Time", "Site")[pvalues == min_p]
kruskal_test.z = rbind(kruskal_test.z, c(i, pvalues, category))
}
kruskal_test.z = data.frame(kruskal_test.z, stringsAsFactors=F)
for(i in 2:4) kruskal_test.z[,i] = as.numeric(kruskal_test.z[,i])
kruskal_test.z = kruskal_test.z %>% mutate(Module=str_to_title(gsub("ME", "", Module)))
write.xlsx(kruskal_test.z, file="2.3_ModuleEigengene.Zooxanthellae.Kruskal-Wallis.UT.xlsx")
kruskal_test.z## Module pvalue.colony pvalue.time pvalue.site Category
## 1 Orangered4 2.702082e-07 5.362137e-01 0.0060108325 Colony
## 2 Royalblue 9.402718e-05 1.628554e-02 0.0123262922 Colony
## 3 Cyan 2.465296e-01 2.217970e-09 0.1324792347 Time
## 4 Lightyellow 5.904067e-03 1.253073e-05 0.0090867092 Time
## 5 Purple 4.092058e-04 4.127219e-05 0.0033310430 Time
## 6 Skyblue3 1.543838e-01 7.051268e-09 0.0206123002 Time
## 7 Darkslateblue 9.133840e-01 3.441319e-08 0.6806134178 Time
## 8 Green 8.756114e-01 3.277916e-09 0.9343579004 Time
## 9 Maroon 1.601848e-01 1.925670e-05 0.1399016678 Time
## 10 Paleturquoise 2.193426e-06 1.958502e-01 0.0359871931 Colony
## 11 Darkorange 1.896955e-06 5.122280e-01 0.0795840679 Colony
## 12 Mediumpurple3 2.102014e-05 6.094637e-03 0.7187224756 Colony
## 13 Thistle1 3.069444e-01 1.962841e-06 0.1674380654 Time
## 14 Red 4.443135e-02 7.592496e-08 0.0345419446 Time
## 15 Brown4 2.721276e-03 3.104265e-09 0.0050099499 Time
## 16 Darkturquoise 6.527110e-01 4.866589e-09 0.3600186559 Time
## 17 Navajowhite2 1.120833e-01 7.532803e-08 0.0386683561 Time
## 18 Salmon4 4.414551e-01 6.929513e-08 0.2086075622 Time
## 19 Lightsteelblue1 9.133408e-01 3.386661e-08 0.7724143679 Time
## 20 Orange 6.742590e-02 3.150048e-08 0.2384770636 Time
## 21 Bisque4 7.589080e-01 3.123259e-07 0.6590444568 Time
## 22 Floralwhite 6.370270e-01 1.094796e-06 0.4490368368 Time
## 23 Lavenderblush3 6.434733e-02 4.574367e-06 0.0421053036 Time
## 24 Skyblue 9.768642e-01 2.401600e-09 0.7410346509 Time
## 25 Thistle2 6.366325e-01 2.254847e-08 0.2054583124 Time
## 26 Plum2 8.117549e-04 5.829597e-08 0.0005323902 Time
## 27 Steelblue 4.163195e-01 1.804720e-08 0.0583196918 TimeVisualization
# Average by biological duplicates
samples = apply(meta[rownames(eigen), c("Month", "Colony")], 1, paste, collapse="_")
figmat = apply(eigen, 2, tapply, samples, mean)
colnames(figmat) = str_to_title(gsub("^ME", "", colnames(figmat)))# Time-dependent modules
figmat.Time = t(figmat[, kruskal_test.z$Module[kruskal_test.z$Category == "Time"]])
modules = rownames(figmat.Time)
figmat.Time = figmat.Time %>% as.data.frame() %>% rownames_to_column("Module") %>%
gather(key="Sample", value="value", -Module) %>%
mutate(Module=factor(Module, levels=modules), Time=meta[Sample,"Month"], Colony=meta[Sample,"Colony"]) %>%
mutate(value=ifelse(value > 0.3, 0.3, ifelse(value < -0.3, -0.3, value)))
figs = list(); size = c(); times = unique(figmat.Time$Time)
for(i in times){
size = c(size, length(unique(figmat.Time$Colony[figmat.Time$Time == i])))
figs[[i]] = ggplot(subset(figmat.Time, Time==i)) +
geom_tile(aes(y=Module, x=Colony, fill=value)) +
scale_fill_gradientn(colors=rev(brewer.pal(9, "RdBu")), limits=c(-0.3, 0.3)) +
theme_minimal() +
theme(legend.position="none", plot.margin = unit(c(0, 0, 0, 0), "cm"),
axis.text.x=element_text(angle=90, vjust=0.5, hjust=1)) +
xlab(NULL) + ylab(NULL)
if(i != times[1]) figs[[i]] = figs[[i]] + theme(axis.text.y=element_blank())
}
size[1] = size[1] + 9
cowplot::plot_grid(plotlist = figs, nrow=1, align="h", rel_widths=size) +
ggsave("2.3_ModuleEigengene.Zooxanthellae.Time.UT.pdf", width=6, height=2.8)
# Proportion of annotated genes
figs = list()
annot.Time = subset(annot, Module %in% modules)
annot.Time$Module = factor(annot.Time$Module, levels=modules)
figs[["size"]] = ggplot(annot.Time) +
geom_bar(aes(y=Size, x=Module, fill=I(tolower(Module))), stat="identity", color="grey60") +
theme_minimal() + coord_flip() + xlab(NULL) + ylab("Module size") +
scale_y_log10(limits=c(1, 40000), breaks=c(1, 100, 10000)) +
theme(axis.title.x=element_text(size=8))
size = c(18)
figs[["annot1"]] = ggplot(annot.Time) +
geom_bar(aes(y=Annotated.Swissprot, x=Module, fill=I(tolower(Module))), stat="identity", color="grey60") +
theme_minimal() + coord_flip() + xlab(NULL) + ylab("Swissprot [%]") +
scale_y_continuous(limits=c(0, 90), breaks=seq(0, 80, by=20)) +
theme(axis.text.y = element_blank(), axis.title.x=element_text(size=8))
size = c(size, 10)
figs[["annot2"]] = ggplot(annot.Time) +
geom_bar(aes(y=Annotated.KAAS, x=Module, fill=I(tolower(Module))), stat="identity", color="grey60") +
theme_minimal() + coord_flip() + xlab(NULL) + ylab("KAAS [%]") +
scale_y_continuous(limits=c(0, 80)) +
theme(axis.text.y = element_blank(), axis.title.x=element_text(size=8))
size = c(size, 10)
cowplot::plot_grid(plotlist = figs, nrow=1, align="h", rel_widths=size) +
ggsave("2.3_ModuleInfo.Zooxanthellae.Time.UT.pdf", width=4, height=2.8)
figmat.Colony = t(figmat[, kruskal_test.z$Module[kruskal_test.z$Category == "Colony"]])
modules = rownames(figmat.Colony)
figmat.Colony = figmat.Colony %>% as.data.frame() %>% rownames_to_column("Module") %>%
gather(key="Sample", value="value", -Module) %>%
mutate(Module=factor(Module, levels=modules), Time=meta[Sample,"Month"], Colony=meta[Sample,"Colony"]) %>%
mutate(value=ifelse(value > 0.3, 0.3, ifelse(value < -0.3, -0.3, value)))
figs = list(); size = c(); colonies = sort(unique(figmat.Colony$Colony))
for(i in colonies){
size = c(size, length(unique(figmat.Colony$Sample[figmat.Colony$Colony == i])))
figs[[i]] = ggplot(subset(figmat.Colony, Colony==i)) +
geom_tile(aes(y=Module, x=Time, fill=value)) +
scale_fill_gradientn(colors=rev(brewer.pal(9, "RdBu")), limits=c(-0.3, 0.3)) +
theme_minimal() +
theme(legend.position="none", plot.margin = unit(c(0, 0, 0, 0), "cm"),
axis.text.x=element_text(angle=90, vjust=0.5, hjust=1)) +
xlab(NULL) + ylab(NULL)
if(i != colonies[1]) figs[[i]] = figs[[i]] + theme(axis.text.y=element_blank())
}
size[1] = size[1] + 8
cowplot::plot_grid(plotlist = figs, nrow=1, align="h", rel_widths=size) +
ggsave("2.3_ModuleEigengene.Zooxanthellae.Colony.UT.pdf", width=6, height=1.1)
figs = list()
# Proportion of annotated genes
annot.Colony = subset(annot, Module %in% modules)
annot.Colony$Module = factor(annot.Colony$Module, levels=modules)
figs[["size"]] = ggplot(annot.Colony) +
geom_bar(aes(y=Size, x=Module, fill=I(tolower(Module))), stat="identity", color="grey60") +
theme_minimal() + coord_flip() + xlab(NULL) + ylab("Module size") +
scale_y_log10(limits=c(1, 40000), breaks=c(1, 100, 10000)) +
theme(axis.title.x=element_text(size=8))
size = c(18)
figs[["annot1"]] = ggplot(annot.Colony) +
geom_bar(aes(y=Annotated.Swissprot, x=Module, fill=I(tolower(Module))), stat="identity", color="grey60") +
theme_minimal() + coord_flip() + xlab(NULL) + ylab("Swissprot [%]") +
scale_y_continuous(limits=c(0, 80)) +
theme(axis.text.y = element_blank(), axis.title.x=element_text(size=8))
size = c(size, 10)
figs[["annot2"]] = ggplot(annot.Colony) +
geom_bar(aes(y=Annotated.KAAS, x=Module, fill=I(tolower(Module))), stat="identity", color="grey60") +
theme_minimal() + coord_flip() + xlab(NULL) + ylab("KAAS [%]") +
scale_y_continuous(limits=c(0, 80)) +
theme(axis.text.y = element_blank(), axis.title.x=element_text(size=8))
size = c(size, 10)
cowplot::plot_grid(plotlist = figs, nrow=1, align="h", rel_widths=size) +
ggsave("2.3_ModuleInfo.Zooxanthellae.Colony.UT.pdf", width=4, height=1.2)
Save results of KW test
save(kruskal_test.c, kruskal_test.z, file="2.3_ModuleEigengene.KWtest.UT.RData")Color bar
colorBar <-function(colors=NULL,
min,
max=-min,
nticks=11,
ticks=seq(min, max, len=nticks),
title="") {
plot(c(min,max), c(0,10), type="n", bty="n", xaxt="n", xlab=title, yaxt="n", ylab="")
axis(1, ticks, las=1, )
scale <- (length(colors)-1)/(max-min)
for (i in 1:(length(colors)-1)) {
x = (i-1)/scale + min
rect(x, 0, x+1/scale, 10, col=colors[i], border=NA)
}
}
colors = rev(colorRampPalette(brewer.pal(9, "RdBu"))(100))
pdf("Legend1.pdf", 3.8, 2.)
colorBar(colors, min=-0.3, max=0.3, nticks=5, title="Module eigengene")
dev.off()## quartz_off_screen
## 2pdf("Legend2.pdf", 3.8, 2.)
colors = rev(colorRampPalette(brewer.pal(11, "PuOr"))(100))
colorBar(colors, min=-0.8, max=0.8, nticks=5, title="Pearson correlation")
dev.off()## quartz_off_screen
## 2