FounderDietStudy_cor.Rmd

---
title: "Founder Diet Study Correlations"
author: "Brian Yandell"
date: "`r format(Sys.time(), '%d %B %Y')`"
output: html_document
params:
  dataset:
    label: "Measurement Set"
    value: plasma
    input: select
    choices: [physio, plasma, liver]
  traits:
    label: "Number of traits"
    value: 3
    input: slider
    min: 2
    max: 6
    
---

```{r setup, include=FALSE, echo = FALSE}
knitr::opts_chunk$set(echo = FALSE, message = FALSE, warning = FALSE,
                      comment = "", fig.height = 7, fig.width = 7)
```

```{r}
library(tidyverse)
library(readxl)
library(readr)
library(foundr) # github/byandell/foundr
```

# Compare raw and adjusted correlations

```{r}
plasmaData <- readRDS("traitData.rds") %>%
  filter(datatype == "plasma") %>%
  rename(dataset = "datatype")
```

```{r}
plasmaSignal <- readRDS("plasmaSignal.rds") %>%
  mutate(dataset = "plasma")
```

```{r}
plasmaDataSignal <- 
  left_join(
    plasmaData %>% select(-dataset),
    plasmaSignal,
    by = c("trait", "strain", "sex"))
```

```{r}
corval <- 
  bind_rows(
    map(
    split(plasmaDataSignal, plasmaDataSignal$trait),
    function(x) {
      terms <- c("signal","mean")
      as.data.frame(cor(x[,"value"], x[,terms], use = "pair"))
    }),
    .id = "trait")
apply(corval[,-1], 2, summary)
```

Pattern below makes sense because `mean` captures most of `value`, so their correlations range up toward 1. The `signal` has no clear pattern of correlation with `value`. However, `signal` is that part of the `value` consists of the terms `strain:condition` and `strain:condition:sex`, the important part of the signal.

```{r}
GGally::ggpairs(corval[,-1])
```


Spearman rank correlation of residuals after removing `strain*sex` and `condition*sex` effects.
That is, looking at correlation due to `strain*condition` and `strain*condition*sex` effects only.

```{r}
if(file.exists(filerds <- "corsraw.rds")) {
  corsraw <- readRDS(filerds)
} else {
  corsraw <- abscorcalc(plasmaData)
  saveRDS(corsraw, filerds)  
}
if(file.exists(filerds <- "cors.rds")) {
  cors <- readRDS(filerds)
} else {
  cors <- abscorcalc(plasmaSignal %>% select(-mean) %>% rename(value = "signal"))
  saveRDS(cors, filerds)  
}
```

```{r eval=FALSE}
traitStats <- readRDS("traitStats.rds")
tmp <- (traitStats %>%
          filter(datatype == "plasma",
                 term == "signal"))$p.value
pvalr <- matrix(tmp, byrow = TRUE, nrow = length(tmp), ncol = length(tmp))
pvalr <- pvalr[upper.tri(pvalr)]
pvalc <- matrix(tmp, byrow = FALSE, nrow = length(tmp), ncol = length(tmp))
pvalc <- pvalc[upper.tri(pvalc)]
```

Each point is a pair of traits showing their raw and adjusted (signal only) correlations. Coloring is by the most significant `p_signal`
for the two traits being correlated. This shows essentially no correlation between  until the raw are fairly large (>0.75). 

Find large correlations for `raw` or `cor_adj` and relatively large `logpmax` for `p_signal`. Add these to plots of correlations (black circles) and use for pair
plots below.

```{r}
plasmaStats <- readRDS("traitStats.rds") %>%
  filter(datatype == "plasma") %>%
  rename(dataset = "datatype")
```

```{r}
dat <- cor_compare(plasmaStats, corsraw, cors)
```

```{r}
ggplot(dat %>% rename(cor_raw = "cor1", cor_adj = "cor2")) +
  aes(cor_raw, cor_adj) +
  geom_hex()
```


```{r}
ggplot(dat %>%
         rename(cor_raw = "cor1", cor_adj = "cor2") %>%
         filter(cor_raw > .75, cor_adj > .75)) +
  aes(cor_raw, cor_adj) +
  geom_hex()
```

# Extreme Comparisons Identified by Correlation and P-value

```{r}
(datsub <- cor_extreme(plasmaStats, dat))
```

```{r}
traits <- c(datsub$traitx, datsub$traity) 
sampleData <- plasmaData %>% filter(trait %in% traits)
for(i in seq_along(traits)) {
  sampleData <- mutate(sampleData, trait = ifelse(trait == traits[i], LETTERS[i], trait))
}
sampleData <- sampleData %>% mutate(condition = ifelse(condition == "HC_LF", "X", "Y"))
sampleData$dataset <- "sample"
anim <- sort(unique(sampleData$animal))
shuf <- sample(anim, length(anim))
names(shuf) <- anim
sampleData$animal <- shuf[sampleData$animal]
sampleData <- sampleData %>%
  group_by(trait) %>%
  mutate(value = (value - mean(value, na.rm = TRUE)) / sd(value, na.rm = TRUE)) %>%
  ungroup()
write_csv(sampleData, "sampleData.csv")
save(sampleData, file = "sampleData.RData")

ggplot(sampleData) + aes(strain, value) + geom_point() + facet_grid(sex+condition ~ trait)
```

```{r}
sampleSignal <- partition(sampleData)
```

## Scatterplots

- C on A: mostly signal
  + mean cor = 0.0107; signal cor = 0.90
  + -10*log(p.value) = 1.97, 2.81
- D on B: negligible signal
  + mean cor = 0.8110; signal cor = 0.28
  + -10*log(p.value) = 0.32, 2.73

#### Plots of Means

```{r}
out <- traitSolos(sampleData, sampleSignal,
                  response = "mean")
plot(out)
#ggplot_extreme(plasmaSignal, "mean", datsub[1,], condition = "sex + condition")
```

```{r}
out2 <- traitPairs(
  out,
  traitnames = attr(out, "traitnames"),
  pair = c(
    paste(attr(out, "traitnames")[1:2], collapse = " ON "),
    paste(attr(out, "traitnames")[3:4], collapse = " ON ")))
plot(out2)
```

```{r}
out <- traitSolos(sampleData, sampleSignal,
                  response = "signal")
plot(out)
#ggplot_extreme(plasmaSignal, "mean", datsub[1,], condition = "sex + condition")
```

```{r}
out2 <- traitPairs(
  out,
  traitnames = attr(out, "traitnames"),
  pair = c(
    paste(attr(out, "traitnames")[1:2], collapse = " ON "),
    paste(attr(out, "traitnames")[3:4], collapse = " ON ")))
plot(out2)
```

# Correlation Comparisons by P-value

```{r}
ggplot(
  dat %>%
    rename(cor_raw = "cor1", cor_adj = "cor2") %>%
    filter(cor_raw > .75, cor_adj > .75) %>%
    mutate(logpmax = pmax(logpvalr, logpvalc)) %>%
    arrange(logpmax)) +
  aes(cor_raw, cor_adj, color = logpmax) +
  geom_point() +
  geom_smooth(col = "black") + 
  scale_colour_gradientn(colours=rainbow(4))
```


```{r}
ggplot(
  dat %>%
    rename(cor_raw = "cor1", cor_adj = "cor2") %>%
    mutate(logpmax = pmax(logpvalr, logpvalc)) %>%
    arrange(logpmax)) +
  aes(cor_raw, cor_adj, color = logpmax) +
  geom_point() +
  geom_hline(yintercept = 0.5, col = "darkgray") +
  geom_smooth(col = "black") + 
  scale_colour_gradientn(colours=rainbow(4)) +
  geom_point(data = datsub, 
    aes(x = cor1, y = cor2), 
    shape = 21, fill = NA, size = 5, col = "black", stroke = 2)
```

Same plot but coloring by the least significant `p_signal` among a pair of traits.

```{r}
ggplot(dat %>%
    rename(cor_raw = "cor1", cor_adj = "cor2") %>%
    mutate(logpmin = pmin(logpvalr, logpvalc)) %>%
    arrange(logpmin)) +
  aes(cor_raw, cor_adj, color = logpmin) +
  geom_point() +
  geom_hline(yintercept = 0.5, col = "darkgray") +
  geom_smooth(col = "black") + 
  scale_colour_gradientn(colours=rainbow(4)) +
  geom_point(data = datsub, 
    aes(x = cor1, y = cor2), 
    shape = 21, fill = NA, size = 5, col = "black", stroke = 2)
```