---
title: "Using Growthcurver"
author: "Kathleen Sprouffske"
date: "`r Sys.Date()`"
output: 
  rmarkdown::html_vignette:
    toc: true
vignette: >
  %\VignetteIndexEntry{Using growthcurver}
  %\VignetteEncoding{UTF-8}
  %\VignetteEngine{knitr::rmarkdown}
editor_options: 
  chunk_output_type: inline
---

## Introduction

`Growthcurver` calculates simple metrics to summarize growth curves.

Growth curves are commonly used in a variety of microbial experiments, including experimental evolution. The data are typically obtained by repeatedly measuring the cell density.  Modern microbial growth curves can be conducted in a plate reader and may result in hundreds of absorbance measurements over the course of 24 hours.

In the `Growthcurver` package, we fit growth curve data to a standard form of the logistic equation common in ecology and evolution whose parameters (the growth rate, the initial population size, and the carrying capacity) provide meaningful population-level information with straight-forward biological interpretation. The logistic equation describes the population size $N_t$ at time $t$ using:

$$
\label{nt}
N_t = \frac{K}{1 + \left( \frac{K-N_0}{N_0} \right) e^{-rt}} \
$$

Here, the population size at the beginning of the growth curve is given by $N_0$. The maximum possible population size in a particular environment, or the carrying capacity, is given by $K$. The intrinsic growth rate of the population, $r$, is the growth rate that would occur if there were no restrictions imposed on total population size. `Growthcurver` finds the best values of $K$, $r$, and $N_0$ for the growth curve data. 

## Input data 
We have provided simulated sample data in the simplest format for `Growthcurver`. The sample data that we've provided has one column for the time (in hours), and one column for each well in a 96-well plate. Each row contains the absorbance reading for each well at a given time, and the rows are sorted by time. 

Preparing your data in this format will allow you to easily adapt our sample simple code, as well as the more complicated pipeline code, to your data.

Below, we show the first few rows and columns of the sample data.

```{r, echo = FALSE, eval = TRUE}
# First, load the package and the dataset. 
library(growthcurver)

# Load the sample growth curve data provided with the package 
# The first column is the time in hours, and there is one column 
# for each well in a 96-well plate.
d <- growthdata
knitr::kable(d[1:10, 1:8])
```

Your data will probably come from the plate reader as an Excel spreadsheet. Convert the data in your spreadsheet so that the wells and time are the column headers and the data are the rows as shown above (Excel's Paste Special | Transpose checkbox is useful for this). Save your growth curve data file as a tab-separated `txt` or `csv` file, and then read that file into R.  

```{r, eval = FALSE}
# Replace the next line with the location and name of your input data file.
file_name <- "the/path/to/my/data/myfilename.txt"
d <- read.table(file_name, header = TRUE, sep = "\t", stringsAsFactors = FALSE)
```

If you plan to analyze whole plates at a time, then the names of the columns of your input data must adhere to the following requirements
* The column containing the time must be named `time`  
* If you are doing background correction using measurements from a blank well, that well must be named `blank` (see the **Background correction** section for more details)
* The remaining columns must have a unique well name that will be eventually be identified as the sample name

In this example, time is reported in units of hours. This means that all metrics involving time for these data are reported by `Growthcurver` in hours (e.g., $r$ is in hours$^{-1}$). If your data are in a different time unit than you would like, the following code may help you adjust the data in your time column.

```{r, eval = FALSE}
# Convert the "time" column from hours to minutes
d$time <- d$time * 60

# Convert the "time" column from minutes to seconds
d$time <- d$time * 60

# Convert the "time" column from seconds to hours
d$time <- d$time / 60 / 60
```


### Background correction


`Growthcurver` provides two methods for doing a background correction to remove the absorbance signal due to the media. The default method (called `min`) finds the minimum value from each well, and subtracts it from all timepoints (for that well only). This method works well if you do not expect your background readings to change over time. 

You can also select a second method (called `blank`), which is useful when the background readings change over time. This may happen, for example, when the media precipitates over the course of the growth curve. In this case, you provide the background correction data in a column called `blank`, which contains a series of measurements in a single media-only well measured at the same timepoints that the experimental measurements were made.

Finally, you can do your own background correction to the data before you call `Growthcurver`. In this case, you specify that you do not need `Growthcurver` to do any background correction for you (called `none`).  I have provided two examples of custom background correction in the example code for  **Customize growth curves for a plate**.

## A simple first example
In this example, we will use `Growthcurver` to summarize a single growth curve in a single well. 

```{r, eval = TRUE}
# First, load the package. 
library(growthcurver)

# Load the sample growth curve data provided in the Growthcurver package.
# The first column is the time in hours, and there is one column 
# for each well in a 96-well plate.
d <- growthdata

# Now, we'll use Growthcurver to summarize the growth curve data using the 
# simple background correction method (minimum value correction). This is the 
# default method, so we don't need to specify it in the command.
# This returns an object of type "gcfit" that holds information about
# the best parameters, the model fit, and additional metrics summarizing
# the growth curve.
gc_fit <- SummarizeGrowth(d$time, d$A1)

# It is easy to get the most useful metrics from a gcfit object, just type:
gc_fit

# And it is easy to plot the raw data and the best fit logistic curve
plot(gc_fit)
```

## Output metrics

The object returned from `SummarizeGrowth` contains the metrics, and you can view them easily (as you did above), or you can manipulate them with R commands (as we will do now). 

```{R, eval = FALSE}
# The gcfit object returned from SummarizeGrowth also contains further metrics 
# summarizing the growth curve data.
gc_fit$vals

# look at the structure of the gc_fit object
str(gc_fit)
```

You'll see that there are three main parts of a `gcfit` object: the `vals`, the   `model`, and the `data`.  `vals` contains the summarized metrics for your growth curver, `model` contains the details of the model fit (for advanced users), and `data` contains the input data to `SummarizeGrowth`.  For most purposes, `vals` is the most useful. Let's see what else we can access in the `vals`.

```{R, eval = TRUE} 
# To see all the available metrics 
str(gc_fit$vals)

# To access a single metric (for example the growth rate r)
gc_fit$vals$r

```

The most useful values are `k`, `n0`, and `r`, which are the values of the parameters for the logistic equation that best fit the data. The fitting algorithm provides a measure of uncertainty for each, which is available (for `n`) in the `n_p` and `n_se` values, for example. The values `sigma` and `df` are both determined during the nonlinear regression fit. `Df` is the degrees of freedom and `sigma` is a measure of the goodnesss of fit of the parameters of the logistic equation for the data; it is the residual standard error from the nonlinear regression model. Smaller `sigma` values indicate a better fit of the logistic curve to the data than larger values. 

`t_mid` is the time at which the population density reaches $\frac{1}{2}K$ (which occurs at the inflection point), `t_gen` is the fastest possible generation time (also called the doubling time), `auc_l` is the area under the logistic curve obtained by taking the integral of the logistic equation, and `auc_e` is the empirical area under the curve which is obtained by summing up the area under the experimental curve from the measurements in the input data. If you decide to use `auc_l` or `auc_e`, make sure that you specify the parameter `t_trim` so that these metrics are comparable across samples or plates that were grown for different lengths of time.

The `note` value provides additional information about problems with fitting the logistic curve to your data. No common problems were identified if it is empty.

## Get growth curves for a plate

One often measures growth curves in a plate reader for many wells at the same time. Following is some sample R code that uses `Growthcurver` to summarize the growth curve data for a whole plate. 

```{R, eval = TRUE} 
# First, load the package and the sample dataset. 
library(growthcurver)
d <- growthdata
```

```{R, eval = FALSE}
# To analyze your data from Excel, you should read your data into the variable
# called d. To do so, replace the next line with the name and location of 
# your input data file.
file_name <- "the/path/to/my/data/myfilename.txt"
d <- read.table(file_name, header = TRUE, sep = "\t", stringsAsFactors = FALSE)

# Make sure that you have a column called "time" (and a column called "blank" 
# if you are using "blanks" for your background correction). See the
# "Input Data" data section of the Vignette if you need help with this.
```
```{R, eval = TRUE}
# Now, we'll use Growthcurver to summarize the growth curve data for the entire
# plate using the default background correction method ("min").
gc_out <- SummarizeGrowthByPlate(d)
```
```{R, eval = FALSE}
# If you would like to use the "blank" background correction, then call
# Growthcurver as follows
gc_out <- SummarizeGrowthByPlate(d, bg_correct = "blank")

# If you would like to generate plots for all of the growth curves in your
# plate, then call Growthcurver as follows. You can change the name of 
# the output file "gc_plots.pdf" to something that makes sense for you.
gc_out <- SummarizeGrowthByPlate(d, plot_fit = TRUE, 
                                 plot_file = "gc_plots.pdf")

# The summary information for each well is listed as a row in the output
# data frame called gc_out.

# We can look at the first few rows in the output using the head command.
head(gc_out)
```
```{R, eval = TRUE, echo = FALSE}
knitr::kable(gc_out[1:5, ])
```
```{R, eval = FALSE}
# Or, you can save the entire data table to a tab-separated file that can be
# imported into Excel.
output_file_name <- "the/path/to/my/data/myfilename.txt"
write.table(gc_out, file = output_file_name, 
            quote = FALSE, sep = "\t", row.names = FALSE)
```

See the section **Output metrics** for more details on interpreting parameters.

## Customize growth curves for a plate

Advanced users may want more control over background correction or plotting of the curves. You can use the following code to customize summarizing your growth curves. To adapt this code to your data, just ensure that your data are in the same format as the example data (discussed previously in the **Input data** section). Most importantly, this code assumes that the time information is stored in a column named `time`.

Here, we will first create a data frame in which to store the output data. We will then loop through the columns in the experimental growth curve data, call `SummarizeGrowth` for each, and store the metrics for each column in the output data frame.

```{r, message = FALSE, fig.width = 7}
# As in the simple example, load the package and the data. 
library(growthcurver)
d <- growthdata

# Let's create an output data frame to store the results in. 
# We'll create it so that it is the right size (it's faster this way!), 
# but leave it empty.
num_analyses <- length(names(d)) - 1
d_gc <- data.frame(sample = character(num_analyses),
                   k = numeric(num_analyses),
                   n0  = numeric(num_analyses),
                   r = numeric(num_analyses),
                   t_mid = numeric(num_analyses),
                   t_gen = numeric(num_analyses),
                   auc_l = numeric(num_analyses),
                   auc_e = numeric(num_analyses),
                   sigma = numeric(num_analyses),
                   stringsAsFactors = FALSE)

# Truncate or trim the input data to observations occuring in the first 20 hours.
# Remember that the times in these sample data are reported in hours. To use  
# minutes (or to trim at a different time), change the next line of code. 
# For example, if you still would like to trim at 20 hours, but your time data 
# are reported in minutes use: trim_at_time <- 20 * 60
trim_at_time <- 20   

# Now, loop through all of the columns in the data frame. For each column,
# run Growthcurver, save the most useful metrics in the output data frame,
# and make a plot of all the growth curve data and their best fits.

# First, create a plot for each of the wells in the 96-well plate.
# Uncomment the next line to save the plots from your 96-well plate to a 
# pdf file in the working directory.
# pdf("growthcurver.pdf", height = 8.5, width = 11)
par(mfcol = c(8,12))
par(mar = c(0.25,0.25,0.25,0.25))
y_lim_max <- max(d[,setdiff(names(d), "time")]) - min(d[,setdiff(names(d), "time")])

n <- 1    # keeps track of the current row in the output data frame
for (col_name in names(d)) {
  
  # Don't process the column called "time". 
  # It contains time and not absorbance data.
  if (col_name != "time") {

    # Create a temporary data frame that contains just the time and current col
    d_loop <- d[, c("time", col_name)]
    
    # Do the background correction.
    # Background correction option 1: subtract the minimum value in a column
    #                                 from all measurements in that column
        min_value <- min(d_loop[, col_name])
    d_loop[, col_name] <- d_loop[, col_name] - min_value
    # Background correction option 2: subtract the mean value of blank wells
    #                                 over the course the experiment
    #                                 (Replace B2, D8, G11 with the column
    #                                  names of your media-only wells)
    #d$blank <- apply(d[, c("B2", "D8", "G11")], 1, mean)
    #d$A1 <- d$A1 - d$blank
    
    # Now, call Growthcurver to calculate the metrics using SummarizeGrowth
    gc_fit <- SummarizeGrowth(data_t = d_loop[, "time"], 
                              data_n = d_loop[, col_name],
                              t_trim = trim_at_time,
                              bg_correct = "none")
    
    # Now, add the metrics from this column to the next row (n) in the 
    # output data frame, and increment the row counter (n)
    d_gc$sample[n] <- col_name
    d_gc[n, 2:9] <- c(gc_fit$vals$k,
                      gc_fit$vals$n0,
                      gc_fit$vals$r,
                      gc_fit$vals$t_mid,
                      gc_fit$vals$t_gen,
                      gc_fit$vals$auc_l,
                      gc_fit$vals$auc_e,
                      gc_fit$vals$sigma)
    n <- n + 1
    
    # Finally, plot the raw data and the fitted curve
    # Here, I'll just print some of the data points to keep the file size smaller
    n_obs <- length(gc_fit$data$t)
    idx_to_plot <- 1:20 / 20 * n_obs
    plot(gc_fit$data$t[idx_to_plot], gc_fit$data$N[idx_to_plot], 
         pch = 20, 
         xlim = c(0, trim_at_time), 
         ylim = c(0, y_lim_max),
         cex = 0.6, xaxt = "n", yaxt = "n")
     text(x = trim_at_time / 4, y = y_lim_max, labels = col_name, pos = 1)
     lines(gc_fit$data$t, predict(gc_fit$model), col = "red")
  }
}
# Uncomment the next line to save the plots from your 96-well plate to a file
# dev.off()
```

After running the above code, the summary metrics are available for each well (each column in the input data corresponds to a single well). Each column of absorbance data is summarized, and the summary is a row in the output data frame that we created.

```{R, eval = FALSE}
# Look at the first few rows (samples) of data in the output data frame. 
# (I'm only showing the first 4 rows of results, but you may want to see more. 
#  You can either look at everything using the command "d_gc", or adjust the
#  number of rows displayed by changing the 4 to something else, 
#  e.g., "d_gc[1:15,]").
d_gc[1:4, ]
```
```{R, eval = TRUE, message = FALSE, echo = FALSE}
library(dplyr)
d_gc[1:4, ] %>% 
    mutate(k = round(k, digits = 5),
         n0 = round(n0, digits = 5), 
         r = round(r, digits = 5),
         t_mid = round(t_mid, digits = 5),
         t_gen = round(t_gen, digits = 5),
         auc_l = round(auc_l, digits = 5),
         auc_e = round(auc_e, digits = 5), 
         sigma = round(sigma, digits = 5))
  
```


## Quality control and best practices

Sometimes `Growthcurver` doesn't find the best fit for your data. This can happen especially in the following cases:

* your population didn't reach stationary phase during  data collection,
* you don't have enough measurements during the initial or final parts of the growth curve (lag phase or stationary phase),
* you didn't do a background correction, or your background correction method doesn't bring the minimum absorbance reading to near zero,
* you are trying to run `Growthcurver` on an empty well.

To see if `Growthcurver` has encountered any of these problems that can lead to a poor fit, we recommend doing these three quality control steps after you have run `Growthcurver`. Each step is covered in more detail below.

 * Plot the curves and your data
 * Check if `Growthcurver` returned any fit notes
 * Identify outliers

### Plot the curves and your data
Firstly, we strongly recommend plotting all the curves and checking them manually. See the sections **A simple first example** and **Get growth curves for a plate** for code demonstrating how to easily plot your data and `Growthcurver`'s fit.

### Check for poor fit notes

`Growthcurver` returns a `note` if it finds a potential problem with the fit of the logistic curve to your data. This may happen when it cannot fit the logistic curve to your data, or if it finds evidence of a questionable fit. For example, `Growthcurver` returns a `note` when the carrying capacity $K$ is greater than the initial population size $N_0$, or when the inflection point `t_mid` is found to be negative (both things should not happen in a well-behaved growth curve!)

You can examine the `note` after fitting the growthcurves for an individual sample, or for an entire plate.

```{R, eval = FALSE, message = FALSE}
# Check if Growthcurver provided any notes in a plate of growthcurves returned 
# from SummarizeGrowthByPlate
gc_out %>% filter(note != "") 

# Check if Growthcurver provided any notes in a single growthcurve returned 
# from SummarizeGrowth
gc_fit$vals$note
```

We simulated this dataset, so it doesn't have any data noisy enough to cause problems with `Growthcurver`. Therefore, there are no `notes` returned for these samples.

### Identify outliers 

You should look for outliers that have unusually large `sigma` values. Each `sigma` value is the residual sum of squares from the fit of the logistic curve to the data, so larger values mean poorer fits. To simplify looking at the `sigma` values, I use the package `dplyr` for data wrangling and exploration.


```{R, eval = TRUE, message = FALSE}
# Load dplyr and the sample output data
library(dplyr)
gc_out <- as_data_frame(gc_out)

# Plot a histogram of the sigma values in order to check for outliers
hist(gc_out$sigma, main = "Histogram of sigma values", xlab = "sigma")

```
``` {R, eval = FALSE, message = FALSE}
# Show the top 5 samples with the largest sigma value 
# (with the worst model fit to the growth curve data)
gc_out %>% top_n(5, sigma) %>% arrange(desc(sigma))
```
```{R, eval = TRUE, echo = FALSE, message = FALSE}
gc_out %>%  
  mutate(k = round(k, digits = 5),
         n0 = round(n0, digits = 5), 
         r = round(r, digits = 5),
         t_mid = round(t_mid, digits = 5),
         t_gen = round(t_gen, digits = 5),
         auc_l = round(auc_l, digits = 5),
         auc_e = round(auc_e, digits = 5), 
         sigma = round(sigma, digits = 5)) %>%
  top_n(5, sigma) %>% arrange(desc(sigma))
```

We simulated this dataset, so it is not very noisy. Therefore, it is not surprising that there aren't any extreme `sigma` outliers in this case.

One additional method for identifying outlier parameters is to conduct a principal components analysis on the samples, which projects high dimensional data on lower dimensional space. Plotting the samples on the first two principal components (PC1 and PC2) can identify natural clusters within the growth curve data as well as outliers.

```{R, eval = TRUE, message = FALSE}
# Load dplyr, ggplot2, and the sample data
library(dplyr)
library(ggplot2)
pca_gc_out <- as_data_frame(gc_out) 

# Prepare the gc_out data for the PCA
rownames(pca_gc_out) <- pca_gc_out$sample

# Do the PCA
pca.res <- prcomp(pca_gc_out %>% select(k:sigma), center=TRUE, scale=TRUE)

# Plot the results
as_data_frame(list(PC1=pca.res$x[,1],
                   PC2=pca.res$x[,2],
                   samples = pca_gc_out$sample)) %>% 
  ggplot(aes(x=PC1,y=PC2, label=samples)) + 
  geom_text(size = 3)
```


Most of the samples in this plot cluster together, except B8 which should be examined more closely to determine if it had a poor fit.
