/
growth_curves_analysis.Rmd
154 lines (140 loc) · 5.34 KB
/
growth_curves_analysis.Rmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
---
title: "Analysis: Growth Curves"
author: "Silverman, Kopf, Gordon, Bebout, Som. Morphological and isotopic changes of heterocystous cyanobacteria in response to N2 partial pressure."
date: "`r format(Sys.Date(), '%d %b %Y')`"
output:
html_document:
css: stylesheet.css
fig_caption: yes
number_sections: yes
toc: yes
toc_float: true
toc_depth: 3
code_folding: show
df_print: paged
subtitle: "Source file: [SI GitHub Repository](http://www.github.com/KopfLab/2018_Silverman_et_al/) / [growth_curves_analysis.Rmd](http://www.github.com/KopfLab/2018_Silverman_et_al/blob/master/growth_curves_analysis.Rmd)"
editor_options:
chunk_output_type: inline
---
```{r setup, echo = FALSE, message=FALSE, warning=FALSE}
library(tidyverse) # dplyr, tidyr, ggplot
library(readxl) # reading excel file
library(openxlsx) # writing excel file
library(knitr) # generating reports
library(latex2exp) # for latex plot labels
library(boot) # bootstrap sample means
library(broom) # processing regression data
source(file.path("lib", "functions.R"))
opts_chunk$set(
dev=c("png", "pdf"), dev.args=list(pdf = list(encoding="WinAnsi", useDingbats=FALSE)),
fig.keep="all", fig.path=file.path("figures", "2018_Silverman_et_al-"))
```
# Load Data
```{r}
OD_data <- read_excel(file.path("data", "2018_Silverman_et_al-OD_data.xlsx"))
```
# Growth curves
## Bootstrapped means for OD data
```{r}
n_bootstrap <- 1000
OD_summary <- OD_data %>%
nest(-organism, -pN2, -type, -day) %>%
mutate(
# number of samples
n_samples = map_int(data, nrow),
# bootstrap data to estimate mean and standard error
bootstrap = map(data, function(.x) {
if (nrow(.x) > 1) boot(data = .x$OD750, statistic = function(x, idx) mean(x[idx]), R = n_bootstrap)
else list(t = .x$OD750)
}),
# bootstrapped based mean and standard error
OD750_mean = map_dbl(bootstrap, ~mean(.x$t)),
OD750_mean_se = map_dbl(bootstrap, ~sd(.x$t)),
# error bars upper and lower
OD750_mean_ue = OD750_mean + OD750_mean_se,
OD750_mean_le = OD750_mean - OD750_mean_se,
# pN2 formatted for legend
pN2_formatted = sprintf("%0.1f bar", pN2)
)
```
## Visualization of all data and means
```{r SI_growth_curves, fig.width = 10, fig.height = 7}
# plot OD both on linear and log scale
bind_rows(
# linear scale data
mutate(OD_summary, scale = "OD750"),
# log scale data
mutate(filter(OD_summary, OD750_mean > 0, day > 0, type == "sample"), scale = "log10 (OD750)",
data = map(data, ~mutate(.x, OD750 = log10(OD750))),
OD750_mean = log10(OD750_mean),
OD750_mean_le = log10(OD750_mean_le),
OD750_mean_ue = log10(OD750_mean_ue))
) %>%
ggplot() +
# plot-wide aesthetics
aes(x = day, y = OD750_mean, fill = pN2_formatted, color = pN2_formatted, shape = type) +
# all data (small translucent points)
geom_point(data = function(df) unnest(df, data), map = aes(y = OD750), size = 2, alpha = 0.5) +
# error bars (1 SE) for averages
geom_errorbar(data = function(df) filter(df, !is.na(OD750_mean_se)),
map = aes(ymin = OD750_mean_le, ymax = OD750_mean_ue), width = 0.2, alpha = 0.6) +
# averages and lines connection them
geom_line(size = 1) +
geom_point(size = 4, color = "black") +
# OD cutoff indicator for growth rate calculations
geom_hline(
data = data_frame(scale = c("OD750", "log10 (OD750)"), cutoff = c(0.4, log10(0.4))),
map = aes(yintercept = cutoff, color = NULL, fill = NULL, shape = NULL)
) +
# panels
facet_grid(scale~organism, space = "free_x", scales = "free") +
# scale
scale_x_continuous(breaks = 0:25) +
scale_color_manual(values = cbPalette) +
scale_fill_manual(values = cbPalette) +
scale_shape_manual(values = c(21:26)) +
# labels
labs(x = "day", y = TeX("$OD_{750}$"), color = expression("pN"[2]), fill = expression("pN"[2]), size = 5) +
# theme
theme_figure(grid = FALSE, text_size = 20) +
guides(fill = guide_legend(override.aes = list(shape = 21)))
```
# Growth rates
## Analysis
```{r calculate_growth_rates}
# calculate growth rates
growth_rates <-
OD_data %>%
filter(day > 0, OD750 > 0, type == "sample", pN2 > 0, OD750 > 0.01, OD750 <= 0.4) %>%
nest(-organism, -pN2) %>%
mutate(
# bootstrap slope from regression model
bootstrap = map(data, ~boot(data = .x, R = n_bootstrap,
statistic = function (x, idx) {
# model fit for resampled data
m <- lm(log(OD750/min(OD750)) ~ day, data = x[idx,])
r2 <- glance(m)$r.squared
slope <- tidy(m)$estimate[2]
return(c(slope, r2))
})),
growth_rate.day = map_dbl(bootstrap, ~mean(.x$t[,1])),
growth_rate_se.day = map_dbl(bootstrap, ~sd(.x$t[,1])),
r2 = map_dbl(bootstrap, ~mean(.x$t[,2]))
) %>%
select(-bootstrap, -data)
growth_rates %>% write.xlsx(file.path("data", "2018_Silverman_et_al-growth_rate_data.xlsx"))
growth_rates
```
## Visualization
```{r SI_growth_rates, fig.width = 8, fig.height = 6}
# plot growth rates
growth_rates %>%
ggplot() +
aes(pN2, growth_rate.day, color = organism, shape = organism) +
geom_errorbar(aes(ymin = growth_rate.day - growth_rate_se.day, ymax = growth_rate.day + growth_rate_se.day), width = 0) +
geom_point(size = 5) +
scale_color_manual("Species", values = cbPalette) +
scale_shape_manual("Species", values = c(16, 15)) +
labs(y = "growth rate [1/day]", x = TeX("$pN_{2}$ \\[bar\\]")) +
theme_figure()
```