---
title: "Avian Muscle Satellite Cell Counts"
output: 
  html_document 
---


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

```

The data presented are from manual cell counts using a haemocytometer. Primary muscle satellite cells were isolated from Ruffs (*Philomachus pugnax*), Sanderlings (*Calidris alba*), Western Sandpipers (*Calidris mauri*) Swainson's Thrushes (*Catharus ustulatus*), and Yellow-rumped Warblers (*Setophaga coronata*). Initially 20,000 cells were plated into 35 mm wells and counted at the timepoints shown below.


```{r Set-up Data}
## Load packages
library(ggplot2)
library(dplyr)
library(growthrates)


## Standard error function
se_calc = function(x){
  sqrt(var(x)/length(x))
}

#### Import and Organize Data ####
##################################################################.

counts <- read.csv("all_species_counts.csv") %>%
  ## Calculate cells/mL and cell total
  mutate(cellsml = count*(Dilution/squares_counted)*10000 ) %>%
  mutate(celltotal = cellsml*volume_ml)

## Summarise technical replicates and remove unused columns
counts <- counts %>% 
  suppressMessages() %>%
  group_by(t, id, bFGF, sp) %>%
  summarise(total = mean(celltotal) ) %>%
  ## Add cells/cm^2 ---> 35mm dish == 9cm^2
  mutate(cell_area = total/9,
         log_total = log10(total),
         log_cell_area = log10(total/9)) %>%
  as.data.frame()

## Summarise data by species
counts_sum <- counts %>%
  suppressMessages() %>%
  group_by(t, sp, bFGF) %>%
  summarise(mean_total = mean(total), 
            sd_total = sd(total), 
            se_total = se_calc(total),
            
            mean_area = mean(cell_area),
            sd_area = sd(cell_area),
            se_area = se_calc(cell_area),
            
            mean_log_total = mean(log_total),
            sd_log_total = sd(log_total),
            se_log_total = se_calc(log_total),
            
            mean_log_area = mean(log_cell_area),
            sd_log_area = sd(log_cell_area),
            se_log_area = se_calc(log_cell_area)) %>%
  as.data.frame()
```

Calculate doubling times

```{r Bird Growth Rates to g}

## Create a column to combine id, bfgf and sp identifiers
counts <- counts %>%
  mutate(name = paste(id, bFGF, sp, sep = "_"))

## generate a list with each item containing a unique identifier
counts_list <- replicate(length(unique(counts$name)), counts, simplify = FALSE)
## Assign the unique identifier to each list name
id_bfgf_sp <- unique(counts$name)
names(counts_list) <- id_bfgf_sp

## Filter each item in list to contain only data matched by the named identifier
for(i in 1:length(counts_list)){
  counts_list[[i]] = counts_list[[i]] %>%
    filter(name == id_bfgf_sp[i])
}

## Calculate splines for each 
counts_spline <- counts_list
counts_spline_dt <- counts_list

for(i in 1:length(counts_list)){
  counts_spline[[i]] <- fit_easylinear(time = counts_list[[i]]$t, 
                                       y = counts_list[[i]]$total,
                                       h = 4)
  counts_spline_dt[[i]] <- log(2)/coef(counts_spline[[i]])[[3]]
}

## collapse doubling times together
spline_dt <- bind_rows(counts_spline_dt, .id = "names") %>%
  t()

spline_dt <- cbind(spline_dt, rownames(spline_dt)) %>%
  as.data.frame()

## Change column names
colnames(spline_dt) <- c("dt", "id_bfgf_sp")

## Convert dt column to a numeric and create empty columns to fill
spline_dt <- spline_dt %>%
  mutate(dt = as.numeric(dt),
         id = "",
         bfgf = "",
         sp = "")

## Replace each identifier column with real values from id_bfgf_sp
for(i in 1:length(spline_dt[,1])){
  spline_dt$id[i] <- strsplit(spline_dt$id_bfgf_sp[i], "_")[[1]][1]
  
  spline_dt$bfgf[i] <- strsplit(spline_dt$id_bfgf_sp[i], "_")[[1]][2]
  
  spline_dt$sp[i] <- strsplit(spline_dt$id_bfgf_sp[i], "_")[[1]][3]
  
}

```

mumax is the maximum growth rate -> convert this to doubling time for comparison

growth rate k is reciprocal to time (usually hr^-1) 
k = ln(N2/N1)/(t2-t1) ... where N is the concentration of cells, and t is time

$$ dt = T*ln(2)/ln(Xe/Xb) $$ -- dt also referred to as g
Where Xb = the first timepoint, Xe = the second timepoint, and T = the total time between the two timepoints.

dt = ln2/k


```{r summarise plot and test dt}
## Summary table
spline_dt_sum <- spline_dt %>%
  group_by(sp, bfgf) %>%
  summarise(mean_dt = mean(dt), sd_dt = sd(dt), se_dt = se_calc(dt))

## Initial plot of doubling times by species and bFGF treatment
ggplot(spline_dt, aes(x = sp, y = dt, colour = bfgf)) + geom_boxplot()

## Models 
lm1.1 <- aov(dt ~ sp*bfgf, data = spline_dt)
anova(lm1.1)

## Check model assumptions
spline_dt_res1.1 <- spline_dt %>%
  mutate(stres = rstandard(lm1.1) )

## Mean residuals approximate 0
ggplot(spline_dt_res1.1, aes(x = sp, y = stres, fill = bfgf), xlab = "Species", ylab = "Standardized Residuals") + geom_boxplot() + theme_classic() + geom_hline(yintercept = 0)

## qqplot to check normality of residuals
qplot(data = spline_dt_res1.1, sample = stres, geom = "qq", xlab = "Theoretical Quantile", ylab = "Sample Quantile") + geom_abline(intercept = 0, slope = 1)
```

Log-transform doubling times to fix model fit.

```{r}
## Log transform doubling time
spline_dt <- spline_dt %>%
  mutate(log_dt = log(dt) )

## Models 
lm1.2 <- aov(dt ~ sp*bfgf, data = spline_dt)
anova(lm1.2)

## Check model assumptions
spline_dt_res1.2 <- spline_dt %>%
  mutate(stres = rstandard(lm1.2) )

## Mean residuals approximate 0
ggplot(spline_dt_res1.2, aes(x = sp, y = stres, fill = bfgf), xlab = "Species", ylab = "Standardized Residuals") + geom_boxplot() + theme_classic() + geom_hline(yintercept = 0)

## qqplot to check normality of residuals
qplot(data = spline_dt_res1.2, sample = stres, geom = "qq", xlab = "Theoretical Quantile", ylab = "Sample Quantile") + geom_abline(intercept = 0, slope = 1)
```


This model isn't valid because of an extreme outlier. SWTH grown without bFGF did not reach exponential growth. 


Remove outlier and reproduce anova.

```{r}

no_swth_spline_dt <- spline_dt %>%
  filter(id_bfgf_sp != "A_N_SWTH", id_bfgf_sp != "B_N_SWTH")

## No swth
lm2.1 <- aov(dt ~ sp*bfgf, data = no_swth_spline_dt)
summary(lm2.1)

## Check model assumptions
spline_dt_res2.1 <- no_swth_spline_dt %>%
  mutate(stres = rstandard(lm2.1) )

## Mean residuals approximate 0
ggplot(spline_dt_res2.1, aes(x = sp, y = stres, fill = bfgf)) + geom_boxplot() + 
  theme_classic() + geom_hline(yintercept = 0)

## qqplot to check normality of residuals 
qplot(data = spline_dt_res2.1, sample = stres, geom = "qq", xlab = "Theoretical Quantile", ylab = "Sample Quantile") + geom_abline(intercept = 0, slope = 1)

```

Interaction is not significant here, remove it from anova

```{r}
lm2.2 <- aov(dt ~ sp + bfgf, data = no_swth_spline_dt)
summary(lm2.2)

## Check model assumptions
spline_dt_res2.2 <- no_swth_spline_dt %>%
  mutate(stres = rstandard(lm2.2) )

## Mean residuals approximate 0
ggplot(spline_dt_res2.2, aes(x = sp, y = stres, fill = bfgf)) + geom_boxplot() + 
  theme_classic() + geom_hline(yintercept = 0)

## qqplot to check normality of residuals 
qplot(data = spline_dt_res2.2, sample = stres, geom = "qq", xlab = "Theoretical Quantile", ylab = "Sample Quantile") + geom_abline(intercept = 0, slope = 1)

```

Using Tukey's honest significant differences for multiple comparisons

```{r}
lm2.2_tukey <- TukeyHSD(lm2.2)

lm2.2_tukey

spline_dt_sum
```

Produce final plot of doubling times

```{r Pub plots}

## Generate label list for showing similarities/differences among species in boxplots

## Grouping label including only one label per species to conform to the plot - this produces missing values warning in the final plot but successfully produces the differences labels. 
label_list <- c(
  c("A", rep(NA, 9)), # RUFF
  c("B", rep(NA, 7)), # WESA
  c("A", rep(NA, 7)), # SAND
  c("B", rep(NA, 7)), # YRWA
  c("B", rep(NA, 3))  # SWTH
)

## Position for label
label_pos_list <- c(
  c(46, rep(NA, 9)), # RUFF
  c(46, rep(NA, 7)), # WESA
  c(46, rep(NA, 7)), # SAND
  c(46, rep(NA, 7)), # YRWA
  c(46, rep(NA, 3))  # SWTH
)

## Add labels to dataframe
spline_plot_dt <- no_swth_spline_dt %>%
  cbind(label_list, label_pos_list)

## Colours to include in plot
bwcolors = c("gray50", "gray15")


noswth_plot <- ggplot(spline_plot_dt, aes(x = sp, y = dt, fill = bfgf)) +
  geom_boxplot() + guides(fill = FALSE) +
  theme_classic() + 
  theme(text = element_text(size = 20)) +
  labs(x = "Species", y = "Doubling Time (hrs)") + 
  scale_fill_manual(values = bwcolors) + 
  scale_y_continuous(breaks = seq(5,45,5), limits = c(5,46)) + 
  geom_text(data = spline_plot_dt, aes(x = sp, y = label_pos_list, label = label_list), size = 7)
  
noswth_plot

## Save plot
#ggsave("Doubling_Time.png", plot = noswth_plot, device = "png", width = 12, height = 9.3, dpi = 300)
```

