Avian Muscle Satellite Cell Counts
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.
## Load packages
library(ggplot2)
library(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(growthrates)## Loading required package: lattice## Loading required package: deSolve## 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()` regrouping output by 't', 'id', 'bFGF' (override with `.groups` argument)## 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()## `summarise()` regrouping output by 't', 'sp' (override with `.groups` argument)Calculate doubling times
## 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
## 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))## `summarise()` regrouping output by 'sp' (override with `.groups` argument)## 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)## Analysis of Variance Table
##
## Response: dt
## Df Sum Sq Mean Sq F value Pr(>F)
## sp 4 645715 161429 41.836 7.342e-12 ***
## bfgf 1 76653 76653 19.865 0.0001074 ***
## sp:bfgf 4 614562 153641 39.818 1.364e-11 ***
## Residuals 30 115758 3859
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 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.
## 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)## Analysis of Variance Table
##
## Response: dt
## Df Sum Sq Mean Sq F value Pr(>F)
## sp 4 645715 161429 41.836 7.342e-12 ***
## bfgf 1 76653 76653 19.865 0.0001074 ***
## sp:bfgf 4 614562 153641 39.818 1.364e-11 ***
## Residuals 30 115758 3859
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## 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.
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)## Df Sum Sq Mean Sq F value Pr(>F)
## sp 4 2548.6 637.1 20.634 3.86e-08 ***
## bfgf 1 218.4 218.4 7.072 0.0126 *
## sp:bfgf 3 133.6 44.5 1.442 0.2508
## Residuals 29 895.5 30.9
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 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
lm2.2 <- aov(dt ~ sp + bfgf, data = no_swth_spline_dt)
summary(lm2.2)## Df Sum Sq Mean Sq F value Pr(>F)
## sp 4 2548.6 637.1 19.813 2.72e-08 ***
## bfgf 1 218.4 218.4 6.791 0.0138 *
## Residuals 32 1029.1 32.2
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## 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
lm2.2_tukey <- TukeyHSD(lm2.2)
lm2.2_tukey## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = dt ~ sp + bfgf, data = no_swth_spline_dt)
##
## $sp
## diff lwr upr p adj
## SAND-RUFF 1.420447 -7.040850 9.881744 0.9882029
## SWTH-RUFF 16.741527 4.049581 29.433473 0.0050122
## WESA-RUFF 15.082128 7.754429 22.409826 0.0000120
## YRWA-RUFF 18.704958 11.377260 26.032657 0.0000002
## SWTH-SAND 15.321080 1.942595 28.699566 0.0183520
## WESA-SAND 13.661681 5.200384 22.122978 0.0004738
## YRWA-SAND 17.284512 8.823214 25.745809 0.0000136
## WESA-SWTH -1.659399 -14.351345 11.032547 0.9954484
## YRWA-SWTH 1.963431 -10.728514 14.655377 0.9913371
## YRWA-WESA 3.622831 -3.704868 10.950529 0.6144218
##
## $bfgf
## diff lwr upr p adj
## Y-N -4.679661 -8.432527 -0.9267952 0.0161412spline_dt_sum## # A tibble: 10 x 5
## # Groups: sp [5]
## sp bfgf mean_dt sd_dt se_dt
## <chr> <chr> <dbl> <dbl> <dbl>
## 1 RUFF N 9.27 1.31 0.585
## 2 RUFF Y 7.88 1.18 0.528
## 3 SAND N 11.4 2.03 1.17
## 4 SAND Y 8.59 0.902 0.521
## 5 SWTH N 856. 339. 240.
## 6 SWTH Y 25.3 2.29 1.62
## 7 WESA N 29.1 8.90 3.98
## 8 WESA Y 18.2 2.79 1.25
## 9 YRWA N 29.1 9.94 4.45
## 10 YRWA Y 25.4 5.59 2.50Produce final plot of doubling times
## 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## Warning: Removed 33 rows containing missing values (geom_text).
## Save plot
#ggsave("Doubling_Time.png", plot = noswth_plot, device = "png", width = 12, height = 9.3, dpi = 300)