enbik.Rmd

---
title: "Informed Design of Experiments?"
subtitle: ""
author: "Martin Modrák"
date: "2018/06/11"
output:
  xaringan::moon_reader:
    lib_dir: libs
    css: ["default", "metropolis-fonts", "slides.css" ]
    nature:
      highlightStyle: github
      highlightLines: true
      countIncrementalSlides: true
---
background-image:url("enbik_img/matrix.jpg")
background-position:50% 50%
class: center bottom inverse

# Simulations!

.copyright[photo by Maurizio Pesce, CC-BY]

---

```{r setup, include=FALSE}
options(htmltools.dir.version = FALSE)
library(tidyverse)
library(knitr)
library(cowplot)
library(scales)
library(DESeq2)
library(here)
library(readr)
source(here::here("simulateDeSeq2.R"))
source(here::here("DeSeq2_helpers.R"))

opts_chunk$set(echo = FALSE)

#Colors to be used: https://coolors.co/f2efea-66d7d1-403d58-dea54b-380a30
```



# Why & What
--

1. Design of experiments

--

  * No. of replicates, comparison groups, ...

--

1. Understanding the methods you use

--

1. Case Studies

--

  * t-test

--

  * DESeq2


---
# Power Analysis

--

* Simulations:

--

  * Easier
  
--

  * Test the whole process
  
--

  * More assumptions

---
background-image: url("enbik_img/dive_in.jpg")
background-position: 50% 0%
class: inverse, center, bottom


.copyright[photo: U.S. government work]


--

# Case Study 1

## Two sample t-test

---

# A Hypothetical Experiment

--

* Cell culture

--

* Does unoptanium increase midichlorian production?

--

* 5 replicates

--

* Analyze with t-test, significant if $p < 0.05$

--

* Simulation assumptions

  * Unoptanium helps ( $+2\mu g$ on average)

--

  * $\mathrm{sd} = 8\mu g$


---


# What do we care about?

--

* Observed effect size

--

* How frequently will we claim significance
  
--

  * a.k.a. power

--

  * But there's more!

--

* Let's simulate 10000 datasets


```{r t_simulate, cache=TRUE}
set.seed(14682456)

result = data.frame(id = 1:10000, effect = NA, p = NA, lower_confidence = NA)
true_effect = 2
sd = 8
sample_size = 5

for(i in 1:10000) {
  baseline = rnorm(sample_size, 0, sd)
  treatment = rnorm(sample_size,  true_effect, sd)
  test_result = t.test(treatment, baseline)
  if(test_result$conf.int[2] < 0) {
    result[i, "lower_confidence"] = test_result$conf.int[2] 
  } else {
    result[i, "lower_confidence"] = test_result$conf.int[1] 
  }
  result[i, "effect"] = mean(treatment) - mean(baseline)
  result[i, "p"] = test_result$p.value
}
```

---
background-image: url("enbik_img/what_could_go_wrong.jpg")
background-position: 50% 0%
background-size: cover

.copyright[photo: U.S. government work]

---

# What We Observe

```{r t_observed_effects, fig.height = 3, dev='svg', warning=FALSE}
hist_scale <- scale_x_continuous(limits = c(-20,20))
large_hist_scale_y <- scale_y_continuous(limits = c(0,1200))
result %>% ggplot(aes(x = effect)) + geom_histogram(bins = 30) + geom_vline(xintercept = true_effect, color = "blue", size = 2) + hist_scale + large_hist_scale_y 
```

---

# Filter for Significance

```{r t_filtered_effects, fig.height = 3, dev='svg', warning=FALSE}
results_sig <- result %>% filter(p < 0.05)  

sig_plot <-  results_sig %>% ggplot(aes(x = effect)) + geom_histogram(bins = 30) + geom_vline(xintercept = true_effect, color = "blue", size = 2)  + hist_scale

sig_plot + large_hist_scale_y

```

--

**Power:**
```{r t_power}
cat("p < 0.05 in",mean(result$p < 0.05),"cases")
```


---

# A Closer Look

```{r t_filtered_zoomed, fig.height = 3, dev='svg', warning=FALSE}
sig_plot <-  results_sig %>% ggplot(aes(x = effect)) + geom_histogram(bins = 30) + geom_vline(xintercept = true_effect, color = "blue", size = 2)  + hist_scale

sig_plot

```

--

**Type S Error** (wrong **S**ign)

--

```{r t_type_s}
results_sig %>% summarize("Type S error" = mean(effect < 0) %>% percent, "95% CI excludes true" = mean(lower_confidence > true_effect) %>% percent) %>% 
  kable(format = "html", format.args = )
```

---
# A Closer Look

```{r t_filtered_zoomed_2, fig.height = 3, dev='svg', warning = FALSE}
sig_plot
```

**Type M Error** (wrong **M**agnitude)

--

```{r t_type_m}
results_sig %>% filter(effect > 0) %>% summarize("Mean exaggeration" = mean(effect) / true_effect, "Min. exaggeration" = min(effect) / true_effect)  %>% kable(format = "html", digits = 1)

```


---
background-image: url("enbik_img/kangaroo.jpg")
background-position: 50% 0%
background-size: 60%
class: center, bottom

# Significance is Not a Savior!

---

# Impact on the Literature

--

* Published effects are exaggerated

--

  * Exaggeration depends on amount of noise
  
--

  * Negligible in high-powered studies
  
--

* If a results looks too good given the noise
--
 it probably is.

---
background-image: url("enbik_img/challenge.jpg")
background-position: 50% 0%
class: center, bottom

.copyright[photo by Llann Wé, CC-BY]

--

# Case Study 2

## Differential Expression (DESeq2)

---

# Less Hypothetical Experiment

--

* Differential expression upon unoptanium stress

--

* Control, treatment, 3 replicates each

--

* 1000 genes

--

* We use DESeq2 to test for effect = $|log_2(fc)| > 1$

---

# Simulating DESeq2

--

* Where do the read counts come from?

--

  * From a previous experiment
  
--

* How to set $log_2(fc)$ ?

--

  * 80% genes have $log_2(fc) = 0$
  
--

  * 0, 2, 4 and 6 for the other 20%
  
--

* 100 simulations each


---
# Some results

```{r DE_simulate, warning=FALSE, message=FALSE, results="hide", cache = TRUE}
set.seed(79314522)
inputs <- data.frame(num_replicates = c(3,3,3,3), effect_size = c(0,2,4,6), lfcThreshold = c(1,1, 1, 1))
results_base <- deSeqMultiTest(inputs, num_simulations = 100)
```


```{r DE_results_1}
avg_func <- function(x) { mean(x, na.rm  = TRUE)}

default_presentation_transform <- function(x) {
  x %>% ungroup() %>% 
    dplyr::rename(log_fc = effect_size) %>%
    select(-num_replicates, -lfcThreshold) %>% 
    kable(format = "html",digits = 1)
}

show_DeSeq_results <- function(results_base) {
  results_base %>% group_by(num_replicates, effect_size, lfcThreshold) %>% 
    summarise(
              "True Pos." = avg_func(TP_), 
              "False Pos." = avg_func(FP_), 
              "Type S error" = avg_func(S_error_),
              #FDR = avg_func(FP_ / (FP_ + TP_)),
              #"Median exaggeration" = median(mean_true_eff_shrunk, na.rm = TRUE),
              "Mean exaggeration" = avg_func(mean_true_eff / effect_size),
              "Mean shrunk exaggeration" = avg_func(mean_true_eff_shrunk / effect_size) 
    ) %>%     
    default_presentation_transform()
}

show_DeSeq_results(results_base)


```

We tested for $|log_2(fc)| > 1$

---

# Replicating DeSeq2 results

--

* Exact experiment replication (3 replicates each)

--

* Replicated = significant in both

---

# Replication results

```{r DE_replication_simulate, warning=FALSE, message=FALSE, results = "hide", cache = TRUE}
set.seed(324588)
inputs_repl <- data.frame(effect_size = c(2,4,6), lfcThreshold = c(1, 1, 1))
results_repl <- deseq_replication_multi(inputs_repl, num_simulations = 100)
```

```{r DE_replication_show}
results_repl %>% group_by(num_replicates, effect_size, lfcThreshold) %>%
  summarise("Significant 1st experiment" = mean(significant_), 
    Replicated =  mean(replicated),
    #"Smaller effect" = mean(smaller_eff, na.rm = TRUE) / 1000,
    "Smaller effect - significant" =  sum(smaller_eff_significant, na.rm = TRUE) / sum(significant_, na.rm = TRUE)
  ) %>%
  default_presentation_transform()

```

---

# DESeq2 Summary

--

* DE experiments have low power

--

* DESeq2 rocks!

--

* DESeq2 avoids false positives at all costs
--
 -> high false negatives

---
class:inverse

# Take Home 

--

* Worry about Type S & M errors

--

* Simulate experiments before investing money

--

* Simulate to understand published research

--

* Code available at https://github.com/cas-bioinf/statistical-simulations

--

.thanks[
Thanks for your attention!
]


---

# What about 6 replicates?

```{r DE_simulate_2, warning=FALSE, message=FALSE, results="hide", cache = TRUE}
set.seed(578195214)
inputs2 <- data.frame(num_replicates = c(6,6,6,6), effect_size = c(0,2,4,6), lfcThreshold = c(1, 1, 1, 1))
results_base2 <- deSeqMultiTest(inputs2, num_simulations = 100)
```

```{r DE_results_2}
show_DeSeq_results(results_base2)
```

We tested for $|log_2(fc)| > 1$