metab.Rmd

---
title: "Metabolites"
author: "Brian Yandell"
date: "`r format(Sys.time(), '%d %B %Y')`"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = FALSE, warning = FALSE, message = FALSE)
```

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

# Trait names across metabolites

Strategy is to try to try to match up names when difference is only about capitalization. Then we look at 0-120 pairs for PlaMet to see how they are associated. In particular, consider the model

<center>
```
PlaMet120 ~ PlaMet0 + strain * sex * diet
PlaMet120 ~ PlaMet0 + signal + rest
```
<\center>

If PlaMet0 explains the signal at least, and ideally the rest,
then all the information about factor effects is there at the beginning and there are no GTT adjustments. If PlaMet0 remains significant after adjusting for signal and rest, then there is additional information in PlaMet0 beyond the experimental design.

Coding will involve adjustment to strainstats to allow for covariate.

Also consider what to do with LivMet. See stuff later in this document.

Goal also is to have beginning of document static and rest dynamic.


```{r}
traitSignal <- readRDS("traitSignal.rds") %>%
  filter(dataset %in% c("LivMet", "PlaMet0", "PlaMet120"))
```

```{r}
trnames <- traitSignal %>%
  mutate(Trait = as.character(trait),
         trait = str_to_title(Trait)) %>%
  select(dataset, trait, Trait) %>%
  distinct(dataset, trait, Trait) %>%
   group_by(dataset, trait) %>%
   summarize(Trait = paste(Trait, collapse = "; "),
             trait = trait[1],
             .groups = "drop") %>%
   ungroup() %>%
   select(dataset, trait, Trait) %>%
   arrange(trait)
```


```{r}
traitct <- traitSignal %>%
   mutate(trait = str_to_title(trait)) %>%
   count(dataset, trait) %>%
   pivot_wider(names_from = "dataset", values_from = "n", values_fill = 0) %>%
   mutate(status = ifelse(LivMet > 0 & PlaMet0 > 0 & PlaMet120 > 0,
                          "all", NA),
          status = ifelse(is.na(status) & (PlaMet120 > 0 & PlaMet0 > 0),
                          "PlaMets", status),
          status = ifelse(is.na(status) & (PlaMet120 > 0 & PlaMet0 == 0),
                          "PlaMet120", status),
          status = ifelse(is.na(status) & (PlaMet0 > 0 & PlaMet120 == 0),
                          "PlaMet0", status))
```

```{r}
table(traitct$status)
```

  Cytidine              192
2 N-Acetylneuraminate   192
3 Quinate               192

```{r}
LivMetData <- readRDS("LivMetData.rds") %>%
  filter(tolower(trait) %in% c("cytidine", "quinate", "n-acetylneuraminate"))
```

```{r}
traitall <- 
  traitct %>%
    filter(status %in% c("PlaMets","all")) %>%
    select(status, trait) %>%
    arrange(trait)
```

```{r}
traitpla <- 
  traitct %>%
    filter(status %in% c("PlaMet0","PlaMet120")) %>%
    select(status, trait) %>%
    arrange(trait)
```

```{r}
traitjoin <- function(traitpla, trnames, dataname) {
  trnames <- trnames %>%
    filter(dataset == dataname) %>%
    select(-dataset)
  names(trnames)[match("Trait", names(trnames))] <-dataname
  
  arrange(
    select(
    left_join(
    traitpla,
    trnames,
    by = "trait"),
    status, everything()),
    trait)
}
```

```{r}
traitpla <- 
  traitjoin(
    traitjoin(
      traitjoin(
        traitpla,
        trnames,
        "PlaMet0"),
      trnames,
      "PlaMet120"),
    trnames,
    "LivMet")    
```

Following data frame has 0 rows since trait names forced into Title case (stringr::str_to_title() in R/Met.R).

```{r}
traitall <- 
  traitjoin(
    traitjoin(
      traitjoin(
        traitall,
        trnames,
        "PlaMet0"),
      trnames,
      "PlaMet120"),
    trnames,
    "LivMet") %>%
  filter(!(PlaMet0 == PlaMet120 & LivMet == PlaMet0))
```

```{r}
traitall <- bind_rows(
    traitall,
    traitpla)
```

```{r eval = FALSE}
write_csv(
  traitall,
  "PlaMetNames.csv"
)
```

# Compare two datasets

```{r}
datasets <- c("PlaMet0", "LivMet")
```

```{r}
traits <- 
  (traitct %>%
     filter(.data[[datasets[1]]] > 0,
            .data[[datasets[2]]] > 0))$trait
```

```{r}
traitData <- readRDS("traitData.rds") %>%
  filter(dataset %in% datasets,
         trait %in% traits) %>%
  mutate(trait = factor(trait, traits))
traitSignal <- readRDS("traitSignal.rds") %>%
  filter(dataset %in% datasets,
         trait %in% traits) %>%
  mutate(trait = factor(trait, traits))
traitStats <- readRDS("traitStats.rds") %>%
  filter(dataset %in% datasets,
         trait %in% traits) %>%
  mutate(trait = factor(trait, traits))
```

```{r}
traitJoin <-
  left_join(
    traitData,
    traitSignal,
    by = c("dataset", "strain", "sex", "condition", "trait"))

```

# Correlation of values

```{r}
cortrait <- function(object, datasets) {
  dplyr::ungroup(
    dplyr::summarize(
      dplyr::group_by(
        tidyr::pivot_wider(
          object,
          names_from = "dataset", values_from = "value"),
        trait),
      value = stats::cor(
        .data[[datasets[1]]],
        .data[[datasets[2]]],
        use = "pair"),
      .groups = "drop"))
}
```

```{r}
cors <- bind_cols(
  # Correlation of values
  traitData %>%
    cortrait(datasets),
  
  # Correlation of cellmean
  traitSignal %>%
    select(-signal) %>%
    rename(value = "cellmean") %>%
    cortrait(datasets) %>%
    select(-trait),
  
  # Correlation of signal
  traitSignal %>%
    select(-cellmean) %>%
    rename(value = "signal") %>%
    cortrait(datasets) %>%
    select(-trait),
  
  # Correlation of rest
  traitSignal %>%
    mutate(value = cellmean - signal) %>%
    select(-signal, -cellmean) %>%
    cortrait(datasets) %>%
    select(-trait),

  # Correlation of individual
  traitJoin %>%
    mutate(value = value - cellmean) %>%
    select(-signal, -cellmean) %>%
    cortrait(datasets) %>%
    select(-trait))
```

```{r}
ggplot(
  cors %>%
    pivot_longer(cellmean:noise,
                 names_to = "response", values_to = "cors") %>%
    mutate(response = factor(response, c("cellmean","signal","ind_signal","noise")))) +
  aes(individual, cors) +
  facet_wrap(~ response) +
  geom_point(shape = 1) +
  geom_smooth(method = "lm", se = FALSE, formula = "y~x", col = "red") +
  geom_abline(slope = 1, intercept = 0, col = "blue") +
  ylab("extracted correlation") +
  xlab("raw correlation")
```


```{r}
ggplot(cors) +
  aes(ind_signal, signal) +
  geom_point(shape = 1) +
  geom_smooth(method = "lm", se = FALSE, formula = "y~x", col = "red") +
  geom_abline(slope = 1, intercept = 0, col = "blue")
```


```{r}
ggplot(cors) +
  aes(ind_signal, noise) +
  geom_point(shape = 1) +
  geom_smooth(method = "lm", se = FALSE, formula = "y~x", col = "red") +
  geom_abline(slope = 1, intercept = 0, col = "blue")
```

# Correlation of Stats

```{r}
(stats <- traitStats %>%
   mutate(p.value = -log10(p.value)) %>%
   pivot_longer(SD:p.value, names_to = "stats", values_to = "value") %>%
   pivot_wider(names_from = "dataset", values_from = "value"))
```

```{r}
CB_colors <- RColorBrewer::brewer.pal(n = 3, name = "Dark2")
c("#1B9E77", "#D95F02", "#7570B3")
```

```{r}
plot_stats <- function(stats, terms = unique(stats$term)) {
  p1 <- ggplot(
    stats %>%
      filter(stats == "p.value",
             term %in% terms) %>%
      mutate(term = factor(term, terms))) +
    aes(PlaMet, LivMet) +
    facet_grid(stats ~ term, scale = "free") +
    geom_point(shape = 1, col = CB_colors[3]) +
    geom_smooth(method = "lm", se = FALSE, formula = "y~x", col = CB_colors[2]) +
    geom_abline(slope = 1, intercept = 0, col = CB_colors[1]) +
    xlab("")
  p2 <- ggplot(
    stats %>%
      filter(stats == "SD",
             term %in% terms) %>%
      mutate(term = factor(term, terms))) +
    aes(PlaMet, LivMet) +
    facet_grid(stats ~ term, scale = "free") +
    geom_point(shape = 1, col = CB_colors[3]) +
    geom_smooth(method = "lm", se = FALSE, formula = "y~x", col = CB_colors[2]) +
    geom_abline(slope = 1, intercept = 0, col = CB_colors[1])
cowplot::plot_grid(p1, p2, nrow = 2)
}
```


```{r}
plot_stats(stats, c("cellmean", "signal", "rest"))
```

```{r}
plot_stats(stats, c("signal", "strain_condition", "strain_sex_condition"))
```

```{r}
plot_stats(stats, c("strain", "sex", "strain_sex"))
```

```{r}
plot_stats(stats, c("condition", "sex_condition"))
```