WGCNA.Rmd

---
title: "WGCNA"
author: "Brian Yandell"
date: "`r format(Sys.time(), '%d %B %Y')`"
output: html_document
params:
  echo: no
  rootdir: ~/founder_diet_study
---

```{r}
rawdir <- file.path(params$rootdir, "RawData")
harmonizeddir <- file.path(params$rootdir, "HarmonizedData")
```

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, warning = FALSE, message = FALSE)
knitr::opts_knit$set(root.dir = harmonizeddir)
```

```{r}
options(stringsAsFactors = FALSE)
getwd()
```

Goal: Have one set of routines in package modulr to create module objects and to read module objects provided by Mark.
Note that `geneTree` object is being changed to `dendro`; need to change already stored objects.
Probably get away for now from value, signal, cellmean, noise construction as too confusing to audience;
keep around for later study, perhaps.
Probably simplify from `listof_...` to modules; later consider generic module object (which may not rely on dendro).

```{r eval = FALSE}
if (!require("BiocManager", quietly = TRUE))
  install.packages("BiocManager")
BiocManager::install("GO.db")
BiocManager::install("impute")
BiocManager::install("preprocessCore")
```

```{r}
library(foundr)
library(tidyverse)
library(modulr)
```

```{r}
datasets <- c("Physio","LivMet","PlaMet","Lipid","LivMet","PlaEnrich","LivRna","LivIso")
pvalmax <- 0.5
```

# Contrasts

```{r}
filename <- file.path(harmonizeddir, "traitContrast.rds")
traitContrast <- list()
for(mod in c(datasets)) {
  dirname <- file.path(harmonizeddir, mod)
  modname <- paste0(mod, "Module")
  cat(modname, "\n")

  traitContrast[[mod]] <- 
    foundr::conditionContrasts(
      dplyr::mutate(
        readRDS(file.path(harmonizeddir, mod, paste0(mod, "Signal.rds"))),
        dataset = mod),
      dplyr::mutate(
        readRDS(file.path(harmonizeddir, mod, paste0(mod, "Stats.rds"))),
        dataset = mod),
      termname = "signal")
}
traitContrast <- dplyr::bind_rows(traitContrast, .id = "dataset")
class(traitContrast) <- c("conditionContrasts", class(traitContrast))
saveRDS(traitContrast, filename)
```

# WGCNA 

Save RDS module information for each dataset.

```{r}
if(!exists("traitContrast")) {
  filename <- file.path(harmonizeddir, "traitContrast.rds")
  traitContrast <- readRDS(filename)
}
```

```{r}
for(mod in datasets) {
  modname <- paste0(mod, "Module")
  cat(modname, "\n")
  filename <- file.path(harmonizeddir, mod, paste0(mod, "Module.rds")) 
  # Create list of WGCNA modules.
  contr <- dplyr::filter(traitContrast, .data$dataset == mod, .data$p.value <= pvalmax)
  # Assign list of WGCNA modules 
  assign(modname, modulr::listof_wgcnaModules(dplyr::select(contr, -p.value)))
    
  # Save listof_wgcnaModules object.
  saveRDS(get(modname), filename)
}
```

Create `traitModule` object as list of dataset module objects.

```{r}
filename <- file.path(harmonizeddir, "traitModule.rds")
traitModule <- list()
for(mod in c(datasets)) {
  dirname <- file.path(harmonizeddir, mod)
  modname <- paste0(mod, "Module")
  traitModule[[mod]] <- readRDS(file.path(dirname, paste0(modname, ".rds")))
}
attr(traitModule, "p.value") <- pvalmax
saveRDS(traitModule, filename)
```

```{r}
summary(traitModule[["LivMet"]], "Both Sexes")
```
```{r}
PlaTraitModule <- traitModule[c("Physio","PlaMet","LivMet","Lipid","Module","Enrich","PlaEnrich","LivEnrich")]
LivTraitModule <- traitModule[c("Physio","LivMet","Lipid","LivRna","LivIso","MixMod")]
```

# Eigentraits

```{r}
if(!exists("traitContrast")) {
  filename <- file.path(harmonizeddir, "traitContrast.rds")
  traitContrast <- readRDS(filename)
}
```

```{r}
if(!exists("traitModule")) {
  filename <- file.path(harmonizeddir, "traitModule.rds")
  traitModule <- readRDS(filename)
}
```

```{r}
filename <- file.path(harmonizeddir, "traitEigen.rds")
traitEigen <- list()
for(mod in c(datasets)) {
  dirname <- file.path(harmonizeddir, mod)
  modname <- paste0(mod, "Eigen")
  cat(modname, "\n")
  traitEigen[[mod]] <- foundr::eigen_contrast(traitModule[[mod]], 
    dplyr::filter(traitContrast, .data$dataset %in% mod))
}
saveRDS(traitEigen, filename)
```

## Specific Examples worked out

```{r}
LivMetEigen <- 
  foundr::eigen_contrast(traitModule$LivMet,
                         traitContrast %>% filter(dataset == "LivMet"))
```

```{r}
summary(LivMetEigen)
```

** Note change on workbench to ggplot_conditionContrasts for ordername.**

```{r}
plot(LivMetEigen, bysex = "Both Sexes")
```

```{r}
plot(LivMetEigen, bysex = "Both Sexes", volcano = TRUE,
                 threshold = c(SD = 0.5))
```

```{r}
plotly::ggplotly(plot(LivMetEigen, bysex = "Both Sexes", volcano = TRUE,
                 threshold = c(SD = 0.5)))
```

Compare eigentrait to traits

```{r}
LivMetModule <- traitModule[["LivMet"]]
LivMetContr <- traitContrast %>% filter(dataset == "LivMet")
```

```{r}
object <- eigen_traits(LivMetModule, "Both Sexes", "B_green",
                       LivMetContr,
                       LivMetEigen)
```

```{r}
object %>% arrange(p.value)
```

```{r}
ggplot(object) + aes(abs(kME), -log10(p.value)) + geom_point()
```

** Plot by p.value as alternative to kME. **
** Assign min p.value to each module. **

```{r}
modules <- (LivMetEigen %>% filter(sex == "Both Sexes") %>% distinct(trait))$trait
names(modules) <- modules
(eigen_p.value <- unlist(purrr::map(modules,
           function(x) min(eigen_traits(LivMetModule, "Both Sexes", x,
                       LivMetContr,
                       LivMetEigen)$p.value, na.rm = TRUE))))

```

```{r}
plot(object, bysex = "Both Sexes")
```

```{r}
summary(LivMetEigen)
```

```{r}
plot(object %>% mutate(kME = abs(kME)), bysex = "Both Sexes", volcano = TRUE,
     threshold = c(kME = 0.8))
```

```{r}
objectp <- object
attr(objectp, "ordername") <- "p.value"
plot(objectp, bysex = "Both Sexes", volcano = TRUE)
```

# Modules on Contrasts

```{r}
LivMetSignal <- readRDS(file.path(harmonizeddir, "traitSignal.rds")) %>%
  filter(dataset == "LivMet")
LivMetStats <- readRDS(file.path(harmonizeddir, "traitStats.rds")) %>%
  filter(dataset == "LivMet")
```

```{r}
sexes <- c("Both Sexes", "Female", "Male", "Sex Contrast")
names(sexes) <- c("F+M", "F", "M", "F-M")
```

```{r}
LivMetContr <- conditionContrasts(LivMetSignal, LivMetStats,
                                  termname = "signal") %>%
  filter(p.value <= 0.05)
```

## Module Topology

```{r}
top <- wgcna_topology(
  LivMetContr %>% 
    select(-p.value) %>%
    filter(sex == "Both Sexes"),
  power = seq(from = 5, to = 30, by = 5))
```

```{r}
ggplot_wgcna_topology(top)
```

```{r}
top <- wgcna_topology(
  LivMetContr %>%
    select(-p.value) %>%
    filter(sex == "Sex Contrast"),
  power = c(5,10,15,20,25,30,35,40))
```

```{r}
ggplot_wgcna_topology(top)
```

```{r}
LivMetMod <- list()
# More than 13 causes crash.
for(size in seq(from = 2, to = 13)) {
  cat(size)
  LivMetMod[[paste0("size=", size)]] <- 
    listof_wgcnaModules(LivMetContr %>% select(-p.value),
                        params = list(power = 30, minSize = size))
}
LivMetMod <- transpose(LivMetMod)
```

```{r}
ggplot_listof_wgcnaModules(LivMetMod[["Both Sexes"]], "size=4")
```

```{r}
LivMetMod <- list()
for(beta in seq(from = 5, to = 30, by = 5)) {
  cat(beta)
  LivMetMod[[paste0("beta=", beta)]] <- 
    listof_wgcnaModules(LivMetContr %>% select(-p.value),
                        params = list(power = beta, minSize = 4))
}
LivMetMod <- transpose(LivMetMod)
```

```{r}
ggplot_listof_wgcnaModules(LivMetMod[["Both Sexes"]], "beta=30")
```

```{r}
LivMetModule <- listof_wgcnaModules(LivMetContr %>% select(-p.value),
                                    params = list(power = 30, minSize = 6))
```

**There are kludges below because using outdated conditionContrasts for now.**

```{r}
LivMetEigen <- eigen_contrast(LivMetModule, LivMetContr)
```

```{r}
summary(LivMetEigen)
```

** Note change on workbench to ggplot_conditionContrasts for ordername.**

```{r}
plot(LivMetEigen, bysex = "Both Sexes")
```

```{r}
plot(LivMetEigen, bysex = "Both Sexes", volcano = TRUE,
                 threshold = c(SD = 0.5))
```

```{r}
plotly::ggplotly(plot(LivMetEigen, bysex = "Both Sexes", volcano = TRUE,
                 threshold = c(SD = 0.5)))
```

Compare eigentrait to traits

```{r}
object <- eigen_traits(LivMetModule, "Both Sexes", "green",
                       LivMetContr,
                       LivMetEigen)
```

```{r}
plot(object, bysex = "Both Sexes")
```

```{r}
summary(LivMetEigen)
```

```{r}
plot(object %>% mutate(kME = abs(kME)), bysex = "Both Sexes", volcano = TRUE,
     threshold = c(kME = 0.8))
```

## To do

Can conditionContrasts be modified to have multiple p.value's?

- download 4 datasets for Mark to play with
- run on LivRna dataset
- create data frame with
  + p.value (signal, diet:strain, diet:strain:sex)
  + modules (all 4)
- correlations among modules

# Contrasts over Time: Enrich

```{r}
EnrichSignal <- readRDS(file.path(harmonizeddir, "traitSignal.rds")) %>%
  filter(dataset == "Enrich")
EnrichStats <- readRDS(file.path(harmonizeddir, "traitStats.rds")) %>%
  filter(dataset == "Enrich")
```

```{r}
EnrichContr <- conditionContrasts(EnrichSignal, EnrichStats,
                                  termname = "signal")
EnrichContrSignal <- foundr:::contrast2signal(EnrichContr)
```

```{r}
(timetrait_selection <- timetraits(EnrichContrSignal, "minute"))
```

```{r}
EnrichTimes <- traitTimes(
        EnrichContrSignal, EnrichContrSignal, EnrichStats,
        timetrait_selection[1], "minute", "cellmean")
```

```{r}
ggplot_traitTimes(EnrichTimes$traits)
```

```{r}
ggplot_traitTimes(EnrichTimes$traits, facet_strain = TRUE)
```

Basic idea:

- clone shinyContrastPanel to shinyContrast
- redesign shinyContrastPanel to call
  + shinyContrast: stuff as is
  + shinyTimesPanel: use stuff above
  + shinyModuleStuff: use stuff above on sex-based modules
  
Other important idea:

- examine WGCNA params, esp beta
- set up graphic using WGCNAcontrasts with one level deeper
  + dataset, sex, beta
- determine if beta (or other parameters matter)

# Mark's WGCNA Modules

All metabolites: liver, plasma at 0min, and plasma at 120min.
Rankz transformation to normalize.
Unsigned modules
beta=12
modSize = 4
cutHeight = 0.25

This is the basis for a read WGCN routine. can use 

```{r}
markHarmony <- function(mod, moddir, modRdata, modparams) {
  paramcode <- modparams
  paramcode$signType <- ifelse(paramcode$signType == "unsigned", "U", "S")
  paramcode <- paste(paramcode, collapse = "")

  modname <- paste0(mod, "Module")
  if(file.exists(
    filename <- file.path(harmonizeddir, mod,
                          paste0(modname, "_", paramcode, ".rds")))) {
    assign(modname, readRDS(filename))
  } else {
    assign(
      modname,
      modulr::load_wgcna(
        moddir,
        modRdata,
        params = modparams))
    if(!dir.exists(file.path(harmonizeddir, mod)))
      dir.create(file.path(harmonizeddir, mod))
    saveRDS(get(modname), filename)
  }
}
```


```{r}
mod <- "Metab"
moddir <- file.path(rawdir, "WGCNA_modules",
                    "WGCNA_modules_all_metabolites")
modRdata <- "WGCNA_objects_ms10.Rdata"
modparams <- list(
        power = 12, 
        signType = "unsigned",
        minSize = 4)
MetabModule <- markHarmony(
  mod,
  moddir,
  modRdata,
  modparams)
```

The following does not work yet, since these data have "A129_nn" naming,
different from Metab. Need to use `annot` to study. Wait for another day.

```{r eval=FALSE}
mod <- "LivGene"
moddir <- file.path(rawdir, "WGCNA_modules",
                    "WGCNA_modules_liverRNA",
                    "Gene\ level\ quantitation\ from\ Colin")
modparams <- list(
        power = 10, 
        signType = "unsigned",
        minSize = 4)
LivGeneModule <- markHarmony(
  mod,
  moddir,
  "WGCNA_objects_ms10_log2_unsigned_beta10.Rdata",
  modparams)

modparams <- list(
        power = 12, 
        signType = "signed",
        minSize = 4)
LivGeneModule <- markHarmony(
  mod,
  moddir,
  "WGCNA_objects_ms10_log2_unsigned_beta10.Rdata",
  modparams)
```

```{r}
getwd()
filename <- "traitModule.rds"
traitModule <- readRDS(filename)
```

```{r}
traitModule$Metab <- MetabModule
```

```{r}
saveRDS(traitModule, file.path("Module", "traitModule.rds"))
```

```{r}
knitr::knit_exit()
```

# Module Harmony

Use code below to create `moduleData`, `moduleStats`, `moduleSignal`.
This does not work for contrasts as cannot (yet) return to full data on 192 mice.

```{r}
moduleData <- 
  dplyr::mutate(
    foundr::moduleHarmony("Module",, traitModule, response = "Both Sexes"),
    dataset = "Module")
```

```{r}
moduleStats <- foundr::strainstats(moduleData)
```

```{r}
moduleSignal <- foundr::partition(moduleData)
```

```{r}
(modstats <- dplyr::arrange(
  tidyr::pivot_wider(
  dplyr::mutate(
    dplyr::select(
      dplyr::filter(
        moduleStats,
        !(trait %in% c("grey")),
        term %in% c("signal","cellmean","rest")),
      -SD),
    p.value = signif(-log10(p.value), 3)),
  names_from = "term", values_from = "p.value"),
dplyr::desc(signal)))
```

```{r}
(traits <-
   dplyr::filter(
     tidyr::unite(
       modstats,
       datatraits,
       dataset, trait,
       sep = ": "),
     signal > 1)$datatraits)
```

```{r}
moduleSolo <- foundr::traitSolos(
  moduleData, moduleSignal, response = "cellmean",
  traitnames = traits[1:4])
```

Cell mean for top LivMet and PlaMet module

```{r}
foundr::ggplot_traitSolos(moduleSolo)
```

```{r}
moduleSolo <- foundr::traitSolos(
  moduleData, moduleSignal, response = "signal",
  traitnames = traits[1:4])
```

Signal for top LivMet and PlaMet module

```{r}
foundr::ggplot_traitSolos(moduleSolo)
```

```{r}
traitData <- readRDS(file.path(harmonizeddir, normdir, "traitData.rds"))
traitSignal <- readRDS(file.path(harmonizeddir, normdir, "traitSignal.rds"))
```

### LivMet greenyellow module

```{r}
(GYLivMet <- 
  dplyr::arrange(
    dplyr::filter(
      traitModule$LivMet$value$modules,
      module == "greenyellow"),
    dplyr::desc(abs(kME))))
```

```{r}
GYSolo <- foundr::traitSolos(
  traitData, traitSignal, response = "cellmean",
  traitnames = paste("LivMet", GYLivMet$trait[1], sep = ": "))
```

```{r}
foundr::ggplot_traitSolos(GYSolo)
```

```{r}
foundr::ggplot_traitSolos(GYSolo, facet_strain = TRUE)
```

### LivMet black module

```{r}
(BlLivMet <- 
  dplyr::arrange(
    dplyr::filter(
      traitModule$LivMet$value$modules,
      module == "black"),
    dplyr::desc(abs(kME))))
```

```{r}
BlSolo <- foundr::traitSolos(
  traitData, traitSignal, response = "cellmean",
  traitnames = paste("LivMet", BlLivMet$trait[2], sep = ": "))
```

```{r}
foundr::ggplot_traitSolos(GYSolo)
```

```{r}
BlSolo <- foundr::traitSolos(
  traitData, traitSignal, response = "signal",
  traitnames = paste("LivMet", BlLivMet$trait[2], sep = ": "))
```

```{r}
foundr::ggplot_traitSolos(GYSolo)
```

## Plots

```{r}
modband <- foundr::module_band(traitModule$LivMet, "signal")
```

```{r}
foundr::ggplot_module_band(modband)
```

```{r warning = FALSE, message = FALSE}
foundr::ggplot_listof_wgcnaModules(traitModule$LivMet)
```
Standard plot.

```{r eval=FALSE}
plot(traitModule$LivMet$value)
```

# Module Comparisons

```{r}
mods <- foundr::module_kMEs(traitModule$LivMet)
```

```{r}
subset(mods, "signal", "cellmean", "blue")
```

```{r}
dplyr::arrange(
    dplyr::count(
      mods,
      value_col),
    dplyr::desc(n))
```


```{r}
ggplot2::autoplot(mods, "value", "cellmean")
```

```{r}
plotly::ggplotly(ggplot2::autoplot(mods, "value", "cellmean", abskME = TRUE,
                  facetmodules = "blue"))
```

# Eigentraits

```{r}
eigens <- foundr::eigen_cor(traitModule$LivMet)
```

```{r}
subset(eigens, "value", "signal")
```

```{r}
ggplot2::autoplot(eigens, "value", "cellmean")
```

The sign of correlation of modules could provide adjustment to direction in kME plots above. However, using the absolute value seems to show the strong relationship (or not).

# Compare Datasets

```{r}
dmods <- list(LivMet = traitModule$LivMet,
              PlaMet0 = traitModule$PlaMet0)
response <- "value"
```

```{r}
valueMet <- foundr::join_dataset_modules(dmods, response)
```

```{r}
metmods <- foundr::module_kMEs(valueMet)
```

```{r}
plotly::ggplotly(
  ggplot2::autoplot(metmods, "LivMet", "PlaMet0"))
```

```{r}
traitData <- 
  dplyr::filter(
      readRDS(file.path(harmonizeddir, normdir, "traitData.rds")),
      dataset %in% c("LivMet","PlaMet0"))
traitSignal <-
    dplyr::filter(
      readRDS(file.path(harmonizeddir, normdir, "traitSignal.rds")),
      dataset %in% c("LivMet","PlaMet0"))
```

```{r}
soloData <- 
  foundr::traitSolos(
    dplyr::filter(
      traitData,
      trait == trait[1]),
    dplyr::filter(
      traitSignal,
      trait == trait[1]))
```

```{r}
DT::datatable(soloData)
```

```{r}
foundr::ggplot_traitPairs(
  foundr::traitPairs(soloData))
```

```{r}
eigenMet <- foundr::eigen_cor(valueMet)
```

```{r}
ggplot2::autoplot(eigenMet, "LivMet", "PlaMet0")
```

# WGCNA from Mark

```{r}
load(file.path(moddir, "WGCNA_objects_ms10.Rdata"))
```

```{r}
names(merge)
```

```{r}
traitModule <- readRDS(file.path(harmonizeddir, "Module", "traitModule.rds"))
```

```{r}
names(traitModule$PlaMet0$value)
```

```{r}
head(traitModule$PlaMet0$value$ID)
```

```{r}
str(merge$dendro)
```

```{r}
str(traitModule$PlaMet0$value$geneTree)
```

```{r}
str(merge$newMEs)
```

```{r}
str(traitModule$PlaMet0$value$eigen)
```

```{r}
str(merge$colors)
```


```{r}
str(traitModule$PlaMet0$value$modules)
```

```{r}
tmp <- modulr::module_factors(kMEs, merge$colors)
str(tmp)
```

# Temp Fix of Legacy Module Objects

This was done and should not be repeated as files have moved to Obsolete.

```{r eval = FALSE}
datasets <- c("LivMet", "PlaMet0", "PlaMet120", "LivRna")
for(dataset in datasets) {
  iModule <- readRDS(file.path(dataset, paste0(dataset, "Module.rds")))
  paramcode <- attr(iModule, "params")
  for(response in names(iModule)) {
    names(iModule[[response]])[match("geneTree", names(iModule[[response]]), nomatch = 0)] <- "dendro"
    iModule[[response]]$params <- paramcode
  }
  paramcode <- paramcode[c("power", "signType", "minSize")]
  paramcode$signType <- ifelse(paramcode$signType == "unsigned", "U", "S")
  paramcode <- paste(paramcode, collapse = "")

  filename <- file.path(dataset, paste0(dataset, "Module_", paramcode, ".rds"))
  saveRDS(iModule, filename)
}
```