Flow Cytometry Analysis
Maggi Brisbin
7/25/2022
Packages
load all necessary packages:
library(flowCore)
library(ggcyto)
library(tidyverse)
library(scales)
library(growthrates)
library(ggpubr)
library(flowStats)
library(kableExtra)
library(magick)
library(patchwork)
"%ni%" <- Negate("%in%")Trial 1
Data
Load Flow Cytometry data:
rnd1fs1<- read.flowSet(path = "20220803_BACTB12_3_1_1200/", pattern = ".fcs",alter.names = T)
rnd1fs2 <- read.flowSet(path = "20220803_BactB12_3_1_1500/", pattern = ".fcs",alter.names = T)
rnd1fs3 <- read.flowSet(path = "20220803_bactB12_3_1_2100/", pattern = ".fcs",alter.names = T)
rnd1fs4 <- read.flowSet(path = "20220804_bactb12_3_2_0900/", pattern = ".fcs",alter.names = T)
rnd1fs5 <- read.flowSet(path = "20220804_bactb12_3_2_1200/", pattern = ".fcs",alter.names = T)Add sample names for each run:
samples <- c( "Rpom_full_cold_1", "Rpom_full_cold_2", "Rpom_full_cold_3","Amac_full_cold_1", "Amac_full_cold_2", "Amac_full_cold_3", "Rpom_noB12_cold_1", "Rpom_noB12_cold_2", "Rpom_noB12_cold_3", "Amac_noB12_cold_1", "Amac_noB12_cold_2","Amac_noB12_cold_3",
"Rpom_full_med_1", "Rpom_full_med_2", "Rpom_full_med_3", "Amac_full_med_1", "Amac_full_med_2", "Amac_full_med_3",
"Rpom_noB12_med_1", "Rpom_noB12_med_2", "Rpom_noB12_med_3", "Amac_noB12_med_1", "Amac_noB12_med_2", "Amac_noB12_med_3",
"Rpom_full_hot_1", "Rpom_full_hot_2", "Rpom_full_hot_3", "Amac_full_hot_1", "Amac_full_hot_2", "Amac_full_hot_3",
"Rpom_noB12_hot_1", "Rpom_noB12_hot_2", "Rpom_noB12_hot_3","Amac_noB12_hot_1", "Amac_noB12_hot_2", "Amac_noB12_hot_3")
pData(rnd1fs1)$well <- paste0(samples, "_1")
pData(rnd1fs2)$well <- paste0(samples, "_2")
pData(rnd1fs3)$well <- paste0(samples, "_3")
pData(rnd1fs4)$well <- paste0(samples, "_4")
pData(rnd1fs5)$well <- paste0(samples, "_5")Create gate for live cells:
g.singlets <- polygonGate(filterId = "Singlets","FL1.A"=c(3e3,1e6,1e6, 3e3),"FL3.A"=c(0,5000, 1e5, 1200)) # define gateCreate log version of gate for plotting and plot log transformed data with log transformed gate:
logg.singlets <- polygonGate(filterId = "Singlets", "FL1.A"=c(log(3e3),log(1e6), log(1e6), log(3e3)),"FL3.A"=c(4, log(5000), log(1e5), log(1200)))
autoplot(transform(rnd1fs1[[1]], `FL1.A` =log(`FL1.A` ), `FL3.A` =log(`FL3.A` ) ), "FL1.A","FL3.A", bins = 400, gate = logg.singlets) +geom_gate(logg.singlets) +geom_stats(type = c("percent", "count"))
Filter flowSet based on gate and then replot: (sanity check)
filteredresult <- flowCore::Subset(rnd1fs1,g.singlets)
autoplot(transform(filteredresult[[1]], `FL1.A` =log(`FL1.A` ), `FL3.A` =log(`FL3.A` ) ), "FL1.A","FL3.A", bins = 400, gate = logg.singlets) +geom_gate(logg.singlets) +geom_stats(type = c("percent"))
plot gate for an entire plate:
autoplot(transform(filteredresult, `FL1.A` =log(`FL1.A` ), `FL3.A` =log(`FL3.A` ) ), "FL1.A","FL3.A", bins = 400, gate = logg.singlets) +geom_gate(logg.singlets) 
Make a forward scatter gate:
rg1 <- rectangleGate("FSC.A"=c(100, 100000), filterId="NonDebris")Check forward scatter gate:
autoplot(rnd1fs1[[5]], x = "FSC.A") +geom_gate(rg1) +geom_stats()
Filter based on forward scatter:
filteredresult <- flowCore::Subset(rnd1fs1,g.singlets)
filteredresult2 <- flowCore::Subset(filteredresult,rg1)plot filtered forward scatter (sanity check):
ggcyto(filteredresult2, aes(x = "FSC.A")) +geom_gate(rg1) +geom_density()
Get the total volume injected into the instrument for each sample at each time point and save as a character vector:
rnd1vols1<- c()
for (i in 1:length(samples)) {
rnd1vols1[i] = rnd1fs1[[i]]@description["$VOL"][[1]]
}
rnd1vols2<- c()
for (i in 1:length(samples)) {
rnd1vols2[i] = rnd1fs2[[i]]@description["$VOL"][[1]]
}
rnd1vols3<- c()
for (i in 1:length(samples)) {
rnd1vols3[i] = rnd1fs3[[i]]@description["$VOL"][[1]]
}
rnd1vols4<- c()
for (i in 1:length(samples)) {
rnd1vols4[i] = rnd1fs4[[i]]@description["$VOL"][[1]]
}
rnd1vols5<- c()
for (i in 1:length(samples)) {
rnd1vols5[i] = rnd1fs5[[i]]@description["$VOL"][[1]]
}Create gatingset from filtered flowset to get event counts for each treatment at each time point. Add a volume column using the volume vectors for each day.
#1
filtered_gs1 <- GatingSet(filteredresult2)
gs_pop_add(filtered_gs1,g.singlets) # add gate to GatingSet
recompute(filtered_gs1)
filtered_ps1 <- gs_pop_get_count_with_meta(filtered_gs1)
filtered_ps1$vol <- rnd1vols1
#2
fr_2 <- flowCore::Subset(rnd1fs2,g.singlets)
fr2_2 <- flowCore::Subset(fr_2,rg1)
filtered_gs2 <- GatingSet(fr2_2)
gs_pop_add(filtered_gs2,g.singlets) # add gate to GatingSet
recompute(filtered_gs2)
filtered_ps2 <- gs_pop_get_count_with_meta(filtered_gs2)
filtered_ps2$vol <- rnd1vols2
#3
fr_3 <- flowCore::Subset(rnd1fs3,g.singlets)
fr2_3 <- flowCore::Subset(fr_3,rg1)
filtered_gs3 <- GatingSet(fr2_3)
gs_pop_add(filtered_gs3, g.singlets) # add gate to GatingSet
recompute(filtered_gs3)
filtered_ps3 <- gs_pop_get_count_with_meta(filtered_gs3)
filtered_ps3$vol <- rnd1vols3
#4
fr_4 <- flowCore::Subset(rnd1fs4,g.singlets)
fr2_4 <- flowCore::Subset(fr_4,rg1)
filtered_gs4 <- GatingSet(fr2_4)
gs_pop_add(filtered_gs4,g.singlets) # add gate to GatingSet
recompute(filtered_gs4)
filtered_ps4 <- gs_pop_get_count_with_meta(filtered_gs4)
filtered_ps4$vol <- rnd1vols4
#5
fr_5 <- flowCore::Subset(rnd1fs5,g.singlets)
fr2_5 <- flowCore::Subset(fr_5,rg1)
filtered_gs5 <- GatingSet(fr2_5)
gs_pop_add(filtered_gs5,g.singlets) # add gate to GatingSet
recompute(filtered_gs5)
filtered_ps5 <- gs_pop_get_count_with_meta(filtered_gs5)
filtered_ps5$vol <- rnd1vols5Bind together count frames and create new column for volume in µl. Divide counts by volume to normalize. Separate names to get columns for grouping and faceting. Add factors for temperature.
results_t1 <- rbind(filtered_ps1, filtered_ps2, filtered_ps3, filtered_ps4, filtered_ps5)
results_t1$vol_ul <- as.numeric(results_t1$vol) / 1000
results_t1$cells_per_ul <- results_t1$Count / results_t1$vol_ul
results_t1 <- results_t1 %>% separate(well,c("Species", "Media", "Temp", "Rep", "Time_point") )
results_t1$Temp <- factor(results_t1$Temp, levels = c("cold", "med", "hot"))Growth Curve Plot
Convert “time points” to hours and plot growth curve:
r2_t1<- results_t1 %>%
mutate(hours = case_when(
Time_point == 1 ~ 3,
Time_point == 2 ~ 6,
Time_point == 3 ~ 12,
Time_point == 4 ~ 24,
Time_point == 5 ~ 27))
write.csv(r2_t1, "r2_t1.csv", row.names = FALSE)r2_t1<- read.csv("r2_t1.csv")
r2_t1$Temp <- factor(r2_t1$Temp, levels = c("cold", "med", "hot"))
r2_t1 %>%
group_by(hours, Media, Species, Temp) %>%
dplyr::summarize(mean_cells_per_ul = mean(cells_per_ul), sd = sd(cells_per_ul)) %>% ggplot(aes(x = hours, y = mean_cells_per_ul, color = Temp, shape = Media)) + facet_grid(Temp~Species) + geom_point(size =2, alpha = .7) + theme_bw() + scale_color_manual(values = c("#2D55CC", "#FAAB33", "#FF4041")) +geom_line(aes(linetype = Media))+
geom_errorbar(aes(ymin=mean_cells_per_ul-sd, ymax=mean_cells_per_ul+sd), width=0) + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
log scale
plot with Log scale on y-axis:
r2_t1 %>%
group_by(hours, Media, Species, Temp) %>%
dplyr::summarize(mean_cells_per_ul = mean(cells_per_ul), sd = sd(cells_per_ul)) %>% ggplot(aes(x = hours, y = mean_cells_per_ul, color = Temp, shape = Media)) + facet_grid(Temp~Species) + geom_point(size =2, alpha = .7) + theme_bw() + scale_color_manual(values = c("#2D55CC", "#FAAB33", "#FF4041")) +
geom_line(aes(linetype = Media))+
geom_errorbar(aes(ymin=mean_cells_per_ul-sd, ymax=mean_cells_per_ul+sd), width=0) + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())+ scale_y_log10()
results_Amac_Hot_Full <- r2_t1 %>% group_by(hours, Media, Species, Temp) %>%
dplyr::summarize(mean_cells_per_ul = mean(cells_per_ul), sd = sd(cells_per_ul)) %>% dplyr::filter(Media == "full") %>% dplyr::filter(Species == "Amac") %>% dplyr::filter(Temp == "hot")
print(results_Amac_Hot_Full)## # A tibble: 5 × 6
## # Groups: hours, Media, Species [5]
## hours Media Species Temp mean_cells_per_ul sd
## <int> <chr> <chr> <fct> <dbl> <dbl>
## 1 3 full Amac hot 10899. 738.
## 2 6 full Amac hot 29318. 2441.
## 3 12 full Amac hot 34562. 3947.
## 4 24 full Amac hot 42676. 8487.
## 5 27 full Amac hot 43231. 8166.results_Amac_Hot_noB12 <- r2_t1 %>% group_by(hours, Media, Species, Temp) %>%
dplyr::summarize(mean_cells_per_ul = mean(cells_per_ul), sd = sd(cells_per_ul)) %>% dplyr::filter(Media == "noB12") %>% dplyr::filter(Species == "Amac") %>% dplyr::filter(Temp == "hot")
print(results_Amac_Hot_noB12)## # A tibble: 5 × 6
## # Groups: hours, Media, Species [5]
## hours Media Species Temp mean_cells_per_ul sd
## <int> <chr> <chr> <fct> <dbl> <dbl>
## 1 3 noB12 Amac hot 3656. 1044.
## 2 6 noB12 Amac hot 18813. 4766.
## 3 12 noB12 Amac hot 18361. 2210.
## 4 24 noB12 Amac hot 17819. 2374.
## 5 27 noB12 Amac hot 18542. 3285.( 43230.92 - 18542.19 ) / 43230.92## [1] 0.5710896Max cell counts (biomass)
Get and plot maximum cell concentrations:
df_mean <- data.frame(r2_t1 %>%
group_by(Media, Species, Temp, Rep) %>%
dplyr::summarize(max_cells_per_ul = max(cells_per_ul)) %>%
dplyr::summarize(mean_max_cells = mean(max_cells_per_ul), sd = sd(max_cells_per_ul)) %>% ungroup())
df_sd <- data.frame(r2_t1 %>%
group_by(Media, Species, Temp, Rep) %>%
dplyr::summarize(max_cells_per_ul = max(cells_per_ul)) %>%
dplyr::summarize(mean_max_cells = mean(max_cells_per_ul), sd = sd(max_cells_per_ul)) %>% ungroup() %>% mutate(ymin = mean_max_cells-sd) %>% mutate(ymax = mean_max_cells+sd))
df_points <- data.frame(r2_t1 %>%
group_by(Media, Species, Temp, Rep) %>%
dplyr::summarize(max_cells_per_ul = max(cells_per_ul)) %>% ungroup())
ggplot() + geom_point(data = df_mean, aes(x = Media, y = mean_max_cells, color = Media), size =10, shape ="-") + theme_bw() + scale_color_manual(values = c( "black", "#6495ed"))+ geom_errorbar(data = df_sd, aes(x = Media, y = mean_max_cells, ymin=ymin, ymax=ymax, color = Media), width=0) + geom_jitter(data =df_points, aes(x =Media, y = max_cells_per_ul, color = Media), shape = 1, width = 0.114, fill = "white") +facet_grid(Species~ Temp, scales = "free_y") + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) + ggtitle("Trial 1")
biomass stats
Amac
amac_biomass <- df_points %>% dplyr::filter(Species == "Amac")
amac_biomass$TempMedia <- paste(amac_biomass$Media, amac_biomass$Temp)Amac biomass difference by media in cold temp treatment:
amac_biomass_cold_noB12 <- amac_biomass %>% dplyr::filter(TempMedia == "noB12 cold")
amac_biomass_cold_full <- amac_biomass %>% dplyr::filter(TempMedia == "full cold")
t.test(amac_biomass_cold_noB12$max_cells_per_ul, amac_biomass_cold_full$max_cells_per_ul)##
## Welch Two Sample t-test
##
## data: amac_biomass_cold_noB12$max_cells_per_ul and amac_biomass_cold_full$max_cells_per_ul
## t = 2.5758, df = 3.9112, p-value = 0.06297
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -482.8366 11505.9798
## sample estimates:
## mean of x mean of y
## 56052.99 50541.42Amac biomass difference by media in mid temp treatment:
amac_biomass_mid_noB12 <- amac_biomass %>% dplyr::filter(TempMedia == "noB12 med")
amac_biomass_mid_full <- amac_biomass %>% dplyr::filter(TempMedia == "full med")
t.test(amac_biomass_mid_noB12$max_cells_per_ul, amac_biomass_mid_full$max_cells_per_ul)##
## Welch Two Sample t-test
##
## data: amac_biomass_mid_noB12$max_cells_per_ul and amac_biomass_mid_full$max_cells_per_ul
## t = 0.40098, df = 2.2303, p-value = 0.7236
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -4135.226 5082.012
## sample estimates:
## mean of x mean of y
## 61630.82 61157.43Amac biomass difference by media in hot temp treatment:
amac_biomass_hot_noB12 <- amac_biomass %>% dplyr::filter(TempMedia == "noB12 hot")
amac_biomass_hot_full <- amac_biomass %>% dplyr::filter(TempMedia == "full hot")
t.test(amac_biomass_hot_noB12$max_cells_per_ul, amac_biomass_hot_full$max_cells_per_ul)##
## Welch Two Sample t-test
##
## data: amac_biomass_hot_noB12$max_cells_per_ul and amac_biomass_hot_full$max_cells_per_ul
## t = -4.6326, df = 2.6477, p-value = 0.02481
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -41081.813 -6092.685
## sample estimates:
## mean of x mean of y
## 19643.67 43230.92Rpom
rpom_biomass <- df_points %>% dplyr::filter(Species == "Rpom")
rpom_biomass$TempMedia <- paste(rpom_biomass$Media, rpom_biomass$Temp)Rpom biomass difference by media in cold temp treatment:
rpom_biomass_cold_noB12 <- rpom_biomass %>% dplyr::filter(TempMedia == "noB12 cold")
rpom_biomass_cold_full <- rpom_biomass %>% dplyr::filter(TempMedia == "full cold")
t.test(rpom_biomass_cold_noB12$max_cells_per_ul, rpom_biomass_cold_full$max_cells_per_ul)##
## Welch Two Sample t-test
##
## data: rpom_biomass_cold_noB12$max_cells_per_ul and rpom_biomass_cold_full$max_cells_per_ul
## t = 1.4125, df = 2.1848, p-value = 0.2834
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -19198.35 40387.16
## sample estimates:
## mean of x mean of y
## 51180.4 40586.0Rpom biomass difference by media in mid temp treatment:
rpom_biomass_mid_noB12 <- rpom_biomass %>% dplyr::filter(TempMedia == "noB12 med")
rpom_biomass_mid_full <- rpom_biomass %>% dplyr::filter(TempMedia == "full med")
t.test(rpom_biomass_mid_noB12$max_cells_per_ul, rpom_biomass_mid_full$max_cells_per_ul)##
## Welch Two Sample t-test
##
## data: rpom_biomass_mid_noB12$max_cells_per_ul and rpom_biomass_mid_full$max_cells_per_ul
## t = 0.63281, df = 3.9539, p-value = 0.5616
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -5558.866 8821.360
## sample estimates:
## mean of x mean of y
## 60360.54 58729.29Rpom biomass difference by media in hot temp treatment:
rpom_biomass_hot_noB12 <- rpom_biomass %>% dplyr::filter(TempMedia == "noB12 hot")
rpom_biomass_hot_full <- rpom_biomass %>% dplyr::filter(TempMedia == "full hot")
t.test(rpom_biomass_hot_noB12$max_cells_per_ul, rpom_biomass_hot_full$max_cells_per_ul)##
## Welch Two Sample t-test
##
## data: rpom_biomass_hot_noB12$max_cells_per_ul and rpom_biomass_hot_full$max_cells_per_ul
## t = -1.4241, df = 2.3808, p-value = 0.2714
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -21914.184 9747.558
## sample estimates:
## mean of x mean of y
## 55386.61 61469.92Cell size
get mean cell sizes:
FSC1<- data.frame(fsApply(filteredresult2,each_col,mean))
FSC2<- data.frame(fsApply(fr2_2,each_col,mean))
FSC3<- data.frame(fsApply(fr2_3,each_col,mean))
FSC4<- data.frame(fsApply(fr2_4,each_col,mean))
FSC5<- data.frame(fsApply(fr2_5,each_col,mean))
FSC1$well <- paste0(samples, "_1")
FSC2$well <- paste0(samples, "_2")
FSC3$well <- paste0(samples, "_3")
FSC4$well <- paste0(samples, "_4")
FSC5$well <- paste0(samples, "_5")
Fscatter_t1<- rbind(FSC1,FSC2,FSC3,FSC4,FSC5) %>% separate(well,c("Species", "Media", "Temp", "Rep", "Time_point"))
fscatt1<- Fscatter_t1 %>%
mutate(hours = case_when(
Time_point == 1 ~ 3,
Time_point == 2 ~ 6,
Time_point == 3 ~ 12,
Time_point == 4 ~ 24,
Time_point == 5 ~ 27))
write.csv(fscatt1, "fscatter_rnd1.csv", row.names = FALSE)plot:
fscatt1<- read.csv("fscatter_rnd1.csv")
fscatt1$Temp <- factor(fscatt1$Temp, levels = c("cold", "med", "hot"))
fscatt1 %>%
group_by(hours, Media, Species, Temp) %>%
dplyr::summarize(mean_size = mean(FSC.A), sd = sd(FSC.A)) %>% ggplot(aes(x = hours, y = mean_size, color = Temp, shape = Media)) + facet_grid(Temp~Species) + geom_point(size =2, alpha = .7) + theme_bw() + scale_color_manual(values = c("#2D55CC", "#FAAB33", "#FF4041")) +geom_line(aes(linetype = Media))+
geom_errorbar(aes(ymin=mean_size-sd, ymax=mean_size+sd), width=0) + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
Trial 2
Data
Data import:
rnd2fs1<- read.flowSet(path = "20220808_bactb12_4_1_0900/", pattern = ".fcs",alter.names = T)
rnd2fs2 <- read.flowSet(path = "20220808_BACTB12_4_1_1000//", pattern = ".fcs",alter.names = T)
rnd2fs3 <- read.flowSet(path = "20220808_bactb12_4_1_1100/", pattern = ".fcs",alter.names = T)
rnd2fs4 <- read.flowSet(path = "20220808_bactb12_4_1_1200/", pattern = ".fcs",alter.names = T)
rnd2fs5 <- read.flowSet(path = "20220808_bactb12_4_1_1300/", pattern = ".fcs",alter.names = T)
rnd2fs6 <- read.flowSet(path = "20220808_bactb12_4_1_1400/", pattern = ".fcs",alter.names = T)
rnd2fs7 <- read.flowSet(path = "20220808_bactb12_4_1_1500/", pattern = ".fcs",alter.names = T)
rnd2fs8 <- read.flowSet(path = "20220808_bactb12_4_1_1600/", pattern = ".fcs",alter.names = T)
rnd2fs9 <- read.flowSet(path = "20220808_bactb12_4_1_1700/", pattern = ".fcs",alter.names = T)
rnd2fs10 <- read.flowSet(path = "20220808_bactb12_4_1_2000/", pattern = ".fcs",alter.names = T)
rnd2fs11 <- read.flowSet(path = "20220808_bactb12_4_1_2100/", pattern = ".fcs",alter.names = T)
rnd2fs12 <- read.flowSet(path = "20220809_bactb12_4_2_0900/", pattern = ".fcs",alter.names = T)
rnd2fs13 <- read.flowSet(path = "20220809_bactb12_4_2_0900_2/", pattern = ".fcs",alter.names = T)Add sample names:
sampleNames3 <- c( "Rpom_full_cold_1", "Rpom_full_cold_2", "Rpom_full_cold_3","Amac_full_cold_1", "Amac_full_cold_2", "Amac_full_cold_3", "Rpom_noB12_cold_1", "Rpom_noB12_cold_2", "Rpom_noB12_cold_3", "Amac_noB12_cold_1", "Amac_noB12_cold_2","Amac_noB12_cold_3",
"Rpom_full_med_1", "Rpom_full_med_2", "Rpom_full_med_3", "Amac_full_med_1", "Amac_full_med_2", "Amac_full_med_3",
"Rpom_noB12_med_1", "Rpom_noB12_med_2", "Rpom_noB12_med_3", "Amac_noB12_med_1", "Amac_noB12_med_2", "Amac_noB12_med_3",
"Rpom_full_hot_1", "Rpom_full_hot_2", "Rpom_full_hot_3", "Amac_full_hot_1", "Amac_full_hot_2", "Amac_full_hot_3",
"Rpom_noB12_hot_1", "Rpom_noB12_hot_2", "Rpom_noB12_hot_3","Amac_noB12_hot_1", "Amac_noB12_hot_2", "Amac_noB12_hot_3")
sampleNames4 <- c( "Rpom_full_cold_1", "Rpom_full_cold_2", "Rpom_full_cold_3","Amac_full_cold_1", "Amac_full_cold_2", "Amac_full_cold_3", "Rpom_noB12_cold_1", "Rpom_noB12_cold_2", "Rpom_noB12_cold_3", "Amac_noB12_cold_1", "Amac_noB12_cold_2","Amac_noB12_cold_3",
"Amac_full_med_1", "Amac_full_med_2", "Amac_full_med_3",
"Amac_noB12_med_1", "Amac_noB12_med_2", "Amac_noB12_med_3",
"Amac_full_hot_1", "Amac_full_hot_2", "Amac_full_hot_3",
"Amac_noB12_hot_1", "Amac_noB12_hot_2", "Amac_noB12_hot_3")
sampleNames5 <- c( "Rpom_full_cold_1", "Rpom_full_cold_2", "Rpom_full_cold_3","Amac_full_cold_1", "Amac_full_cold_2", "Amac_full_cold_3", "Rpom_noB12_cold_1", "Rpom_noB12_cold_2", "Rpom_noB12_cold_3", "Amac_noB12_cold_1", "Amac_noB12_cold_2","Amac_noB12_cold_3")
sampleNames6 <- c( "Amac_full_med_1", "Amac_full_med_2", "Amac_full_med_3",
"Amac_noB12_med_1", "Amac_noB12_med_2", "Amac_noB12_med_3",
"Amac_full_hot_1", "Amac_full_hot_2", "Amac_full_hot_3",
"Amac_noB12_hot_1", "Amac_noB12_hot_2", "Amac_noB12_hot_3")
pData(rnd2fs1)$well <- paste0(sampleNames3, "_1")
pData(rnd2fs2)$well <- paste0(sampleNames3, "_2")
pData(rnd2fs3)$well <- paste0(sampleNames3, "_3")
pData(rnd2fs4)$well <- paste0(sampleNames3, "_4")
pData(rnd2fs5)$well <- paste0(sampleNames3, "_5")
pData(rnd2fs6)$well <- paste0(sampleNames3, "_6")
pData(rnd2fs7)$well <- paste0(sampleNames3, "_7")
pData(rnd2fs8)$well <- paste0(sampleNames3, "_8")
pData(rnd2fs9)$well <- paste0(sampleNames3, "_9")
pData(rnd2fs10)$well <- paste0(sampleNames4, "_10")
pData(rnd2fs11)$well <- paste0(sampleNames4, "_11")
pData(rnd2fs12)$well <- paste0(sampleNames5, "_12")
pData(rnd2fs13)$well <- paste0(sampleNames6, "_13")Double check gates:
autoplot(transform(rnd2fs1[[1]], `FL1.A` =log(`FL1.A` ), `FL3.A` =log(`FL3.A` ) ), "FL1.A","FL3.A", bins = 400, gate = logg.singlets) +geom_gate(logg.singlets) +geom_stats(type = c("percent", "count"))
Filter flowSet based on gate and then replot: (sanity check)
filteredresult <- flowCore::Subset(rnd2fs1, g.singlets)
autoplot(transform(filteredresult[[1]], `FL1.A` =log(`FL1.A` ), `FL3.A` =log(`FL3.A` ) ), "FL1.A","FL3.A", bins = 400, gate = logg.singlets) +geom_gate(logg.singlets) +geom_stats(type = c("percent"))
Forward scatter gate:
autoplot(rnd2fs1[[1]], x = "FSC.A") +geom_gate(rg1) +geom_stats()
Filter flowSet based on gate and then replot: (sanity check)
filteredresult <- flowCore::Subset(rnd2fs1,g.singlets)
filteredresult_2 <- flowCore::Subset(filteredresult, rg1)autoplot(filteredresult_2[[1]], x = "FSC.A") +geom_gate(rg1)
Whole plate:
ggcyto(filteredresult_2, aes(x = "FSC.A")) +geom_gate(rg1) +geom_density()
Get sample volumes for each time point:
rnd2vols1<- c()
for (i in 1:length(sampleNames3)) {
rnd2vols1[i] = rnd2fs1[[i]]@description["$VOL"][[1]]
}
rnd2vols2<- c()
for (i in 1:length(sampleNames3)) {
rnd2vols2[i] = rnd2fs2[[i]]@description["$VOL"][[1]]
}
rnd2vols3<- c()
for (i in 1:length(sampleNames3)) {
rnd2vols3[i] = rnd2fs3[[i]]@description["$VOL"][[1]]
}
rnd2vols4<- c()
for (i in 1:length(sampleNames3)) {
rnd2vols4[i] = rnd2fs4[[i]]@description["$VOL"][[1]]
}
rnd2vols5<- c()
for (i in 1:length(sampleNames3)) {
rnd2vols5[i] = rnd2fs5[[i]]@description["$VOL"][[1]]
}
rnd2vols6<- c()
for (i in 1:length(sampleNames3)) {
rnd2vols6[i] = rnd2fs6[[i]]@description["$VOL"][[1]]
}
rnd2vols7<- c()
for (i in 1:length(sampleNames3)) {
rnd2vols7[i] = rnd2fs7[[i]]@description["$VOL"][[1]]
}
rnd2vols8<- c()
for (i in 1:length(sampleNames3)) {
rnd2vols8[i] = rnd2fs8[[i]]@description["$VOL"][[1]]
}
rnd2vols9<- c()
for (i in 1:length(sampleNames3)) {
rnd2vols9[i] = rnd2fs9[[i]]@description["$VOL"][[1]]
}
rnd2vols10<- c()
for (i in 1:length(sampleNames4)) {
rnd2vols10[i] = rnd2fs10[[i]]@description["$VOL"][[1]]
}
rnd2vols11<- c()
for (i in 1:length(sampleNames4)) {
rnd2vols11[i] = rnd2fs11[[i]]@description["$VOL"][[1]]
}
rnd2vols12<- c()
for (i in 1:length(sampleNames5)) {
rnd2vols12[i] = rnd2fs12[[i]]@description["$VOL"][[1]]
}
rnd2vols13<- c()
for (i in 1:length(sampleNames6)) {
rnd2vols13[i] = rnd2fs13[[i]]@description["$VOL"][[1]]
}create gating sets for each timepoint:
filtered_gs1 <- GatingSet(filteredresult_2)
gs_pop_add(filtered_gs1,g.singlets) # add gate to GatingSet
recompute(filtered_gs1)
filtered_ps1 <- gs_pop_get_count_with_meta(filtered_gs1)
filtered_ps1$vol <- rnd2vols1
fr_2 <- flowCore::Subset(rnd2fs2,g.singlets)
fr2_2 <- flowCore::Subset(fr_2,rg1)
filtered_gs2 <- GatingSet(fr2_2)
gs_pop_add(filtered_gs2,g.singlets) # add gate to GatingSet
recompute(filtered_gs2)
filtered_ps2 <- gs_pop_get_count_with_meta(filtered_gs2)
filtered_ps2$vol <- rnd2vols2
#3
fr_3 <- flowCore::Subset(rnd2fs3,g.singlets)
fr2_3 <- flowCore::Subset(fr_3,rg1)
filtered_gs3 <- GatingSet(fr2_3)
gs_pop_add(filtered_gs3, g.singlets) # add gate to GatingSet
recompute(filtered_gs3)
filtered_ps3 <- gs_pop_get_count_with_meta(filtered_gs3)
filtered_ps3$vol <- rnd2vols3
#4
fr_4 <- flowCore::Subset(rnd2fs4,g.singlets)
fr2_4 <- flowCore::Subset(fr_4,rg1)
filtered_gs4 <- GatingSet(fr2_4)
gs_pop_add(filtered_gs4,g.singlets) # add gate to GatingSet
recompute(filtered_gs4)
filtered_ps4 <- gs_pop_get_count_with_meta(filtered_gs4)
filtered_ps4$vol <- rnd2vols4
#5
fr_5 <- flowCore::Subset(rnd2fs5,g.singlets)
fr2_5 <- flowCore::Subset(fr_5,rg1)
filtered_gs5 <- GatingSet(fr2_5)
gs_pop_add(filtered_gs5,g.singlets) # add gate to GatingSet
recompute(filtered_gs5)
filtered_ps5 <- gs_pop_get_count_with_meta(filtered_gs5)
filtered_ps5$vol <- rnd2vols5
#6
fr_6 <- flowCore::Subset(rnd2fs6,g.singlets)
fr2_6 <- flowCore::Subset(fr_6,rg1)
filtered_gs6 <- GatingSet(fr2_6)
gs_pop_add(filtered_gs6,g.singlets) # add gate to GatingSet
recompute(filtered_gs6)
filtered_ps6 <- gs_pop_get_count_with_meta(filtered_gs6)
filtered_ps6$vol <- rnd2vols6
#7
fr_7 <- flowCore::Subset(rnd2fs7,g.singlets)
fr2_7 <- flowCore::Subset(fr_7,rg1)
filtered_gs7 <- GatingSet(fr2_7)
gs_pop_add(filtered_gs7,g.singlets) # add gate to GatingSet
recompute(filtered_gs7)
filtered_ps7 <- gs_pop_get_count_with_meta(filtered_gs7)
filtered_ps7$vol <- rnd2vols7
#8
fr_8 <- flowCore::Subset(rnd2fs8,g.singlets)
fr2_8 <- flowCore::Subset(fr_8,rg1)
filtered_gs8 <- GatingSet(fr2_8)
gs_pop_add(filtered_gs8,g.singlets) # add gate to GatingSet
recompute(filtered_gs8)
filtered_ps8 <- gs_pop_get_count_with_meta(filtered_gs8)
filtered_ps8$vol <- rnd2vols8
#9
fr_9 <- flowCore::Subset(rnd2fs9,g.singlets)
fr2_9 <- flowCore::Subset(fr_9,rg1)
filtered_gs9 <- GatingSet(fr2_9)
gs_pop_add(filtered_gs9,g.singlets) # add gate to GatingSet
recompute(filtered_gs9)
filtered_ps9 <- gs_pop_get_count_with_meta(filtered_gs9)
filtered_ps9$vol <- rnd2vols9
#10
fr_10 <- flowCore::Subset(rnd2fs10,g.singlets)
fr2_10 <- flowCore::Subset(fr_10,rg1)
filtered_gs10 <- GatingSet(fr2_10)
gs_pop_add(filtered_gs10,g.singlets) # add gate to GatingSet
recompute(filtered_gs10)
filtered_ps10 <- gs_pop_get_count_with_meta(filtered_gs10)
filtered_ps10$vol <- rnd2vols10
#11
fr_11 <- flowCore::Subset(rnd2fs11,g.singlets)
fr2_11 <- flowCore::Subset(fr_11,rg1)
filtered_gs11 <- GatingSet(fr2_11)
gs_pop_add(filtered_gs11,g.singlets) # add gate to GatingSet
recompute(filtered_gs11)
filtered_ps11 <- gs_pop_get_count_with_meta(filtered_gs11)
filtered_ps11$vol <- rnd2vols11
#12
fr_12 <- flowCore::Subset(rnd2fs12,g.singlets)
fr2_12 <- flowCore::Subset(fr_12,rg1)
filtered_gs12 <- GatingSet(fr2_12)
gs_pop_add(filtered_gs12,g.singlets) # add gate to GatingSet
recompute(filtered_gs12)
filtered_ps12 <- gs_pop_get_count_with_meta(filtered_gs12)
filtered_ps12$vol <- rnd2vols12
#13
fr_13 <- flowCore::Subset(rnd2fs13,g.singlets)
fr2_13 <- flowCore::Subset(fr_13,rg1)
filtered_gs13 <- GatingSet(fr2_13)
gs_pop_add(filtered_gs13,g.singlets) # add gate to GatingSet
recompute(filtered_gs13)
filtered_ps13 <- gs_pop_get_count_with_meta(filtered_gs13)
filtered_ps13$vol <- rnd2vols13Growth Curve Plot
Put together time points into single dataframe, compute concentrations, split sample names for grouping and faceting:
results <- rbind(filtered_ps1, filtered_ps2, filtered_ps3, filtered_ps4, filtered_ps5, filtered_ps6, filtered_ps7, filtered_ps8, filtered_ps9, filtered_ps10, filtered_ps11, filtered_ps12, filtered_ps13)
results$vol_ul <- as.numeric(results$vol) / 1000
results$cells_per_ul <- results$Count / results$vol_ul
results <- results %>% separate(well,c("Species", "Media", "Temp", "Rep", "Time_point") )
results$Temp <- factor(results$Temp, levels = c("cold", "med", "hot"))convert time points to hours and plot!
r2<- results %>%
mutate(hours = case_when(
Time_point == 1 ~ 0,
Time_point == 2 ~ 1,
Time_point == 3 ~ 2,
Time_point == 4 ~ 3,
Time_point == 5 ~ 4,
Time_point == 6 ~ 5,
Time_point == 7 ~ 6,
Time_point == 8 ~ 7,
Time_point == 9 ~ 8,
Time_point == 10 ~ 9,
Time_point == 11 ~ 12,
Time_point == 12 ~ 24,
Time_point == 13 ~ 24
))
write.csv(r2, "r2_t2.csv", row.names = FALSE)r2<- read.csv("r2_t2.csv")
r2$Temp <- factor(r2$Temp, levels = c("cold", "med", "hot"))
r2 %>% group_by(hours, Media, Species, Temp) %>%
dplyr::summarize(mean_cells_per_ul = mean(cells_per_ul), sd = sd(cells_per_ul)) %>% ggplot(aes(x = hours, y = mean_cells_per_ul, color = Temp, shape = Media)) + facet_grid(Temp~Species) + geom_point(size =2, alpha = .7) + theme_bw() + scale_color_manual(values = c("#2D55CC", "#FAAB33", "#FF4041")) + geom_line(aes(linetype = Media))+
geom_errorbar(aes(ymin=mean_cells_per_ul-sd, ymax=mean_cells_per_ul+sd), width=0) + scale_y_continuous(labels = function(x) format(x, scientific = TRUE))
log scale
r2 %>%
group_by(hours, Media, Species, Temp) %>%
dplyr::summarize(mean_cells_per_ul = mean(cells_per_ul), sd = sd(cells_per_ul)) %>% dplyr::filter(!(Temp == "cold" &hours ==2)) %>% dplyr::filter(!(Temp == "hot" & hours ==1)) %>% ggplot(aes(x = hours, y = mean_cells_per_ul, color = Temp, shape = Media)) + facet_grid(Temp~Species) + geom_point(size =2, alpha = .7) + theme_light() + scale_color_manual(values = c("#2D55CC", "#FAAB33", "#FF4041")) +
geom_line(aes(linetype = Media))+
geom_errorbar(aes(ymin=mean_cells_per_ul-sd, ymax=mean_cells_per_ul+sd), width=0) + scale_y_log10() + theme_bw() + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) + xlim(0, 26)
results_Amac_Hot_Full <- r2 %>% group_by(hours, Media, Species, Temp) %>%
dplyr::summarize(mean_cells_per_ul = mean(cells_per_ul), sd = sd(cells_per_ul)) %>% dplyr::filter(Media == "full") %>% dplyr::filter(Species == "Amac") %>% dplyr::filter(Temp == "hot")
print(results_Amac_Hot_Full)## # A tibble: 12 × 6
## # Groups: hours, Media, Species [12]
## hours Media Species Temp mean_cells_per_ul sd
## <int> <chr> <chr> <fct> <dbl> <dbl>
## 1 0 full Amac hot 823. 24.5
## 2 1 full Amac hot 2117. 193.
## 3 2 full Amac hot 1117. 23.7
## 4 3 full Amac hot 10059. 218.
## 5 4 full Amac hot 23305. 1428.
## 6 5 full Amac hot 30140. 2021.
## 7 6 full Amac hot 34347. 2073.
## 8 7 full Amac hot 36114. 2482.
## 9 8 full Amac hot 39178. 3673.
## 10 9 full Amac hot 40039. 3503.
## 11 12 full Amac hot 39376. 3680.
## 12 24 full Amac hot 35815. 2528.results_Amac_Hot_noB12 <- r2 %>% group_by(hours, Media, Species, Temp) %>%
dplyr::summarize(mean_cells_per_ul = mean(cells_per_ul), sd = sd(cells_per_ul)) %>% dplyr::filter(Media == "noB12") %>% dplyr::filter(Species == "Amac") %>% dplyr::filter(Temp == "hot")
print(results_Amac_Hot_noB12)## # A tibble: 12 × 6
## # Groups: hours, Media, Species [12]
## hours Media Species Temp mean_cells_per_ul sd
## <int> <chr> <chr> <fct> <dbl> <dbl>
## 1 0 noB12 Amac hot 775. 82.9
## 2 1 noB12 Amac hot 1776. 130.
## 3 2 noB12 Amac hot 1331. 31.2
## 4 3 noB12 Amac hot 6201. 1777.
## 5 4 noB12 Amac hot 14712. 2561.
## 6 5 noB12 Amac hot 23356. 1129.
## 7 6 noB12 Amac hot 27518. 1888.
## 8 7 noB12 Amac hot 29538. 1803.
## 9 8 noB12 Amac hot 29211. 1075.
## 10 9 noB12 Amac hot 31299. 970.
## 11 12 noB12 Amac hot 30923. 1707.
## 12 24 noB12 Amac hot 27738. 1147.biomass reduction
(40039 - 31299) / 40039## [1] 0.218287222%
Max cell counts (biomass)
df_mean <- data.frame(r2 %>%
group_by(Media, Species, Temp, Rep) %>%
dplyr::summarize(max_cells_per_ul = max(cells_per_ul)) %>%
dplyr::summarize(mean_max_cells = mean(max_cells_per_ul), sd = sd(max_cells_per_ul)) %>% ungroup())
df_sd <- data.frame(r2 %>%
group_by(Media, Species, Temp, Rep) %>%
dplyr::summarize(max_cells_per_ul = max(cells_per_ul)) %>%
dplyr::summarize(mean_max_cells = mean(max_cells_per_ul), sd = sd(max_cells_per_ul)) %>% ungroup() %>% mutate(ymin = mean_max_cells-sd) %>% mutate(ymax = mean_max_cells+sd))
df_points <- data.frame(r2 %>%
group_by(Media, Species, Temp, Rep) %>%
dplyr::summarize(max_cells_per_ul = max(cells_per_ul)) %>% ungroup())
ggplot() + geom_point(data = df_mean, aes(x = Media, y = mean_max_cells, color = Media), size =10, shape ="-") + theme_bw() + scale_color_manual(values = c( "black", "#6495ed"))+ geom_errorbar(data = df_sd, aes(x = Media, y = mean_max_cells, ymin=ymin, ymax=ymax, color = Media), width=0) + geom_jitter(data =df_points, aes(x =Media, y = max_cells_per_ul, color = Media), shape = 1, width = 0.114, fill = "white") +facet_grid(Species~ Temp, scales = "free_y") + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) + ggtitle("Trial 2")
biomass stats
Amac
amac_biomass <- df_points %>% dplyr::filter(Species == "Amac")
amac_biomass$TempMedia <- paste(amac_biomass$Media, amac_biomass$Temp)Amac biomass difference by media in cold temp treatment:
amac_biomass_cold_noB12 <- amac_biomass %>% dplyr::filter(TempMedia == "noB12 cold")
amac_biomass_cold_full <- amac_biomass %>% dplyr::filter(TempMedia == "full cold")
t.test(amac_biomass_cold_noB12$max_cells_per_ul, amac_biomass_cold_full$max_cells_per_ul)##
## Welch Two Sample t-test
##
## data: amac_biomass_cold_noB12$max_cells_per_ul and amac_biomass_cold_full$max_cells_per_ul
## t = 4.0004, df = 2.5674, p-value = 0.0372
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## 1096.749 16737.342
## sample estimates:
## mean of x mean of y
## 47626.44 38709.40Amac biomass difference by media in med temp treatment:
amac_biomass_mid_noB12 <- amac_biomass %>% dplyr::filter(TempMedia == "noB12 med")
amac_biomass_mid_full <- amac_biomass %>% dplyr::filter(TempMedia == "full med")
t.test(amac_biomass_mid_noB12$max_cells_per_ul, amac_biomass_mid_full$max_cells_per_ul)##
## Welch Two Sample t-test
##
## data: amac_biomass_mid_noB12$max_cells_per_ul and amac_biomass_mid_full$max_cells_per_ul
## t = -0.85157, df = 2.1333, p-value = 0.4793
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -12433.705 8117.794
## sample estimates:
## mean of x mean of y
## 60282.25 62440.21Amac biomass difference by media in hot temp treatment:
amac_biomass_hot_noB12 <- amac_biomass %>% dplyr::filter(TempMedia == "noB12 hot")
amac_biomass_hot_full <- amac_biomass %>% dplyr::filter(TempMedia == "full hot")
t.test(amac_biomass_hot_noB12$max_cells_per_ul, amac_biomass_hot_full$max_cells_per_ul)##
## Welch Two Sample t-test
##
## data: amac_biomass_hot_noB12$max_cells_per_ul and amac_biomass_hot_full$max_cells_per_ul
## t = -4.5772, df = 2.4798, p-value = 0.02923
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -17482.352 -2099.897
## sample estimates:
## mean of x mean of y
## 31449.15 41240.28Rpom
rpom_biomass <- df_points %>% dplyr::filter(Species == "Rpom")
rpom_biomass$TempMedia <- paste(rpom_biomass$Media, rpom_biomass$Temp)Rpom biomass difference by media in cold temp treatment:
rpom_biomass_cold_noB12 <- rpom_biomass %>% dplyr::filter(TempMedia == "noB12 cold")
rpom_biomass_cold_full <- rpom_biomass %>% dplyr::filter(TempMedia == "full cold")
t.test(rpom_biomass_cold_noB12$max_cells_per_ul, rpom_biomass_cold_full$max_cells_per_ul)##
## Welch Two Sample t-test
##
## data: rpom_biomass_cold_noB12$max_cells_per_ul and rpom_biomass_cold_full$max_cells_per_ul
## t = 1.6876, df = 3.2092, p-value = 0.1841
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -1147.687 3953.320
## sample estimates:
## mean of x mean of y
## 28700.86 27298.05Amac biomass difference by media in mid temp treatment:
rpom_biomass_mid_noB12 <- rpom_biomass %>% dplyr::filter(TempMedia == "noB12 med")
rpom_biomass_mid_full <- rpom_biomass %>% dplyr::filter(TempMedia == "full med")
t.test(rpom_biomass_mid_noB12$max_cells_per_ul, rpom_biomass_mid_full$max_cells_per_ul)##
## Welch Two Sample t-test
##
## data: rpom_biomass_mid_noB12$max_cells_per_ul and rpom_biomass_mid_full$max_cells_per_ul
## t = -1.9538, df = 2.1821, p-value = 0.1792
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -19691.69 6715.77
## sample estimates:
## mean of x mean of y
## 25413.64 31901.60Amac biomass difference by media in hot temp treatment:
rpom_biomass_hot_noB12 <- rpom_biomass %>% dplyr::filter(TempMedia == "noB12 hot")
rpom_biomass_hot_full <- rpom_biomass %>% dplyr::filter(TempMedia == "full hot")
t.test(rpom_biomass_hot_noB12$max_cells_per_ul, rpom_biomass_hot_full$max_cells_per_ul)##
## Welch Two Sample t-test
##
## data: rpom_biomass_hot_noB12$max_cells_per_ul and rpom_biomass_hot_full$max_cells_per_ul
## t = -0.67715, df = 2.193, p-value = 0.5627
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -14173.54 10033.59
## sample estimates:
## mean of x mean of y
## 39253.28 41323.26maximum growth rates
All line fits:
all_line_fits <- all_easylinear(cells_per_ul ~ hours|Species+Temp+Rep+Media,
data = (r2 %>% dplyr::filter(!(Temp == "cold" & hours ==2)) %>% dplyr::filter(!(Temp == "hot" & hours ==1))), h=3, quota = 0.95)
par(mfrow = c(6, 6))
par(mar = c(2.5, 4, 2, 1))
plot(all_line_fits, log = "y")
Rpom:
rpom_line_fits <- all_easylinear(cells_per_ul ~ hours|Species+Temp+Rep+Media,
data = (r2 %>% dplyr::filter(!(Temp == "cold" & hours ==2)) %>% dplyr::filter(!(Temp == "hot" & hours ==1)) %>% dplyr::filter(Species == "Rpom")), h=3, quota = 0.95)
#pdf("rpom_line_fits.pdf")
par(mfrow = c(6, 3))
par(mar = c(2.5, 4, 2, 1))
plot(rpom_line_fits, log = "y")
#dev.off()kable((data.frame(results(rpom_line_fits)) %>% select(Species, Temp, Rep, Media, mumax, r2) %>% arrange(Temp)), format = "html", digits = 3) %>%
kable_styling("striped") #%>% save_kable("rpom_line_fits.png", density = 10000)| Species | Temp | Rep | Media | mumax | r2 | |
|---|---|---|---|---|---|---|
| Rpom:cold:1:full | Rpom | cold | 1 | full | 0.487 | 0.916 |
| Rpom:cold:2:full | Rpom | cold | 2 | full | 0.500 | 0.861 |
| Rpom:cold:3:full | Rpom | cold | 3 | full | 0.524 | 0.872 |
| Rpom:cold:1:noB12 | Rpom | cold | 1 | noB12 | 0.605 | 0.936 |
| Rpom:cold:2:noB12 | Rpom | cold | 2 | noB12 | 0.636 | 0.857 |
| Rpom:cold:3:noB12 | Rpom | cold | 3 | noB12 | 0.626 | 0.921 |
| Rpom:hot:1:full | Rpom | hot | 1 | full | 1.744 | 0.923 |
| Rpom:hot:2:full | Rpom | hot | 2 | full | 1.740 | 0.921 |
| Rpom:hot:3:full | Rpom | hot | 3 | full | 1.750 | 0.931 |
| Rpom:hot:1:noB12 | Rpom | hot | 1 | noB12 | 1.677 | 0.927 |
| Rpom:hot:2:noB12 | Rpom | hot | 2 | noB12 | 1.622 | 0.950 |
| Rpom:hot:3:noB12 | Rpom | hot | 3 | noB12 | 1.654 | 0.947 |
| Rpom:med:1:full | Rpom | med | 1 | full | 0.763 | 1.000 |
| Rpom:med:2:full | Rpom | med | 2 | full | 0.796 | 0.998 |
| Rpom:med:3:full | Rpom | med | 3 | full | 0.815 | 0.998 |
| Rpom:med:1:noB12 | Rpom | med | 1 | noB12 | 0.819 | 0.993 |
| Rpom:med:2:noB12 | Rpom | med | 2 | noB12 | 0.785 | 0.999 |
| Rpom:med:3:noB12 | Rpom | med | 3 | noB12 | 0.810 | 0.997 |
Amac:
amac_line_fits <- all_easylinear(cells_per_ul ~ hours|Species+Temp+Rep+Media,
data = (r2 %>% dplyr::filter(!(Temp == "cold" & hours ==2)) %>% dplyr::filter(!(Temp == "hot" & hours ==1)) %>% dplyr::filter(Species == "Amac")), h=3, quota = 0.95)
#pdf("amac_line_fits.pdf")
par(mfrow = c(6, 3))
par(mar = c(2.5, 4, 2, 1))
plot(amac_line_fits, log = "y")
#dev.off()kable((data.frame(results(amac_line_fits)) %>% select(Species, Temp, Rep, Media, mumax, r2) %>% arrange(Temp)), format = "html", digits = 3) %>%
kable_styling("striped") #%>% save_kable("amac_line_fits.png", density = 10000)| Species | Temp | Rep | Media | mumax | r2 | |
|---|---|---|---|---|---|---|
| Amac:cold:1:full | Amac | cold | 1 | full | 0.716 | 0.905 |
| Amac:cold:2:full | Amac | cold | 2 | full | 0.682 | 0.881 |
| Amac:cold:3:full | Amac | cold | 3 | full | 0.679 | 0.842 |
| Amac:cold:1:noB12 | Amac | cold | 1 | noB12 | 0.543 | 0.917 |
| Amac:cold:2:noB12 | Amac | cold | 2 | noB12 | 0.539 | 0.931 |
| Amac:cold:3:noB12 | Amac | cold | 3 | noB12 | 0.597 | 0.961 |
| Amac:hot:1:full | Amac | hot | 1 | full | 1.503 | 0.935 |
| Amac:hot:2:full | Amac | hot | 2 | full | 1.512 | 0.932 |
| Amac:hot:3:full | Amac | hot | 3 | full | 1.541 | 0.945 |
| Amac:hot:1:noB12 | Amac | hot | 1 | noB12 | 1.204 | 0.970 |
| Amac:hot:2:noB12 | Amac | hot | 2 | noB12 | 1.103 | 0.998 |
| Amac:hot:3:noB12 | Amac | hot | 3 | noB12 | 1.282 | 0.952 |
| Amac:med:1:full | Amac | med | 1 | full | 0.912 | 0.998 |
| Amac:med:2:full | Amac | med | 2 | full | 0.926 | 0.998 |
| Amac:med:3:full | Amac | med | 3 | full | 0.895 | 0.992 |
| Amac:med:1:noB12 | Amac | med | 1 | noB12 | 0.824 | 0.973 |
| Amac:med:2:noB12 | Amac | med | 2 | noB12 | 0.734 | 0.975 |
| Amac:med:3:noB12 | Amac | med | 3 | noB12 | 0.741 | 0.974 |
Plot mean growth rates with raw data as points:
many_mus <- rbind(results(amac_line_fits), results(rpom_line_fits))
many_mus$Temp <- factor(many_mus$Temp, levels = c("cold", "med", "hot"))
df_mean <- data.frame(many_mus %>% dplyr::filter(Species == "Amac") %>% group_by(Species, Temp, Media) %>%
dplyr::summarize(mean_mumax = mean(mumax), sd = sd(mumax)) %>% ungroup())
df_sd <- data.frame(many_mus %>% dplyr::filter(Species == "Amac") %>% group_by(Species, Temp, Media) %>%
dplyr::summarize(mean_mumax = mean(mumax), sd = sd(mumax)) %>% ungroup()) %>% mutate(ymin = mean_mumax-sd) %>% mutate(ymax = mean_mumax+sd)
df_points <- data.frame(many_mus %>% dplyr::filter(Species == "Amac"))
amac_mumax<- ggplot() + geom_point(data = df_mean, aes(x = Media, y = mean_mumax, color = Media), size =10, shape ="-") + theme_bw() + scale_color_manual(values = c( "black", "#6495ed"))+ geom_errorbar(data = df_sd, aes(x = Media, y = mean_mumax, ymin=ymin, ymax=ymax, color = Media), width=0) + geom_jitter(data =df_points, aes(x =Media, y = mumax, color = Media), shape = 1, width = 0.114, fill = "white") +facet_grid(~Temp) + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) + ggtitle("A. macleodii") +ylim(0.5, 1.7)df_mean <- data.frame(many_mus %>% dplyr::filter(Species == "Rpom") %>% group_by(Species, Temp, Media) %>%
dplyr::summarize(mean_mumax = mean(mumax), sd = sd(mumax)) %>% ungroup())
df_sd <- data.frame(many_mus %>% dplyr::filter(Species == "Rpom") %>% group_by(Species, Temp, Media) %>%
dplyr::summarize(mean_mumax = mean(mumax), sd = sd(mumax)) %>% ungroup()) %>% mutate(ymin = mean_mumax-sd) %>% mutate(ymax = mean_mumax+sd)
df_points <- data.frame(many_mus %>% dplyr::filter(Species == "Rpom"))
rpom_mumax <- ggplot() + geom_point(data = df_mean, aes(x = Media, y = mean_mumax, color = Media), size =10, shape ="-") + theme_bw() + scale_color_manual(values = c( "black", "#6495ed"))+ geom_errorbar(data = df_sd, aes(x = Media, y = mean_mumax, ymin=ymin, ymax=ymax, color = Media), width=0) + geom_jitter(data =df_points, aes(x =Media, y = mumax, color = Media), shape = 1, width = 0.114, fill = "white") +facet_grid(~Temp) + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) + ggtitle("R. pomeroyi") +ylim(0.4, 1.9)(amac_mumax / rpom_mumax) + plot_layout(guides = "collect")
Growth rate comparison stats
Amac
amac_mu <- many_mus %>% dplyr::filter(Species == "Amac")
amac_mu$TempMedia <- paste(amac_mu$Media, amac_mu$Temp)Amac growth rate difference by media in cold temp treatment:
#cold
amac_mu_cold_noB12 <- amac_mu %>% dplyr::filter(TempMedia == "noB12 cold")
amac_mu_cold_full <- amac_mu %>% dplyr::filter(TempMedia == "full cold")
t.test(amac_mu_cold_noB12$mumax, amac_mu_cold_full$mumax)##
## Welch Two Sample t-test
##
## data: amac_mu_cold_noB12$mumax and amac_mu_cold_full$mumax
## t = -5.9585, df = 3.4317, p-value = 0.006391
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -0.19870978 -0.06658427
## sample estimates:
## mean of x mean of y
## 0.5596798 0.6923268Amac growth rate difference by media in mid temp treatment:
#warm
amac_mu_warm_noB12 <- amac_mu %>% dplyr::filter(TempMedia == "noB12 med")
amac_mu_warm_full <- amac_mu %>% dplyr::filter(TempMedia == "full med")
t.test(amac_mu_warm_noB12$mumax, amac_mu_warm_full$mumax)##
## Welch Two Sample t-test
##
## data: amac_mu_warm_noB12$mumax and amac_mu_warm_full$mumax
## t = -4.7573, df = 2.387, p-value = 0.02901
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -0.25717828 -0.03219524
## sample estimates:
## mean of x mean of y
## 0.7665456 0.9112323Amac growth rate difference by media in hot temp treatment:
#warm
amac_mu_hot_noB12 <- amac_mu %>% dplyr::filter(TempMedia == "noB12 hot")
amac_mu_hot_full <- amac_mu %>% dplyr::filter(TempMedia == "full hot")
t.test(amac_mu_hot_noB12$mumax, amac_mu_hot_full$mumax)##
## Welch Two Sample t-test
##
## data: amac_mu_hot_noB12$mumax and amac_mu_hot_full$mumax
## t = -6.0668, df = 2.1915, p-value = 0.02089
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -0.5326034 -0.1117969
## sample estimates:
## mean of x mean of y
## 1.196487 1.518687####Rpom
Rpom growth rate difference by media in hot temp treatment:
rpom_mu <- many_mus %>% dplyr::filter(Species == "Rpom")
rpom_mu$TempMedia <- paste(rpom_mu$Temp, rpom_mu$Media)
rpom_mu_hot_noB12 <- rpom_mu %>% dplyr::filter(TempMedia == "hot noB12")
rpom_mu_hot_full <- rpom_mu %>% dplyr::filter(TempMedia == "hot full")
t.test(rpom_mu_hot_noB12$mumax, rpom_mu_hot_full$mumax)##
## Welch Two Sample t-test
##
## data: rpom_mu_hot_noB12$mumax and rpom_mu_hot_full$mumax
## t = -5.7599, df = 2.1483, p-value = 0.02442
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -0.15862142 -0.02802538
## sample estimates:
## mean of x mean of y
## 1.651267 1.744591Rpom growth rate difference by media in mid temp treatment:
rpom_mu_mid_noB12 <- rpom_mu %>% dplyr::filter(TempMedia == "med noB12")
rpom_mu_mid_full <- rpom_mu %>% dplyr::filter(TempMedia == "med full")
t.test(rpom_mu_mid_noB12$mumax, rpom_mu_mid_full$mumax)##
## Welch Two Sample t-test
##
## data: rpom_mu_mid_noB12$mumax and rpom_mu_mid_full$mumax
## t = 0.73964, df = 3.4638, p-value = 0.5065
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -0.04074797 0.06796380
## sample estimates:
## mean of x mean of y
## 0.8047011 0.7910931Rpom growth rate difference by media in cool temp treatment:
rpom_mu_cool_noB12 <- rpom_mu %>% dplyr::filter(TempMedia == "cold noB12")
rpom_mu_cool_full <- rpom_mu %>% dplyr::filter(TempMedia == "cold full")
t.test(rpom_mu_cool_noB12$mumax, rpom_mu_cool_full$mumax)##
## Welch Two Sample t-test
##
## data: rpom_mu_cool_noB12$mumax and rpom_mu_cool_full$mumax
## t = 8.5117, df = 3.8698, p-value = 0.001208
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## 0.07939223 0.15780072
## sample estimates:
## mean of x mean of y
## 0.6223185 0.5037220cell size
Mean
FSC1<- data.frame(fsApply(filteredresult_2,each_col,mean))
FSC2<- data.frame(fsApply(fr2_2,each_col,mean))
FSC3<- data.frame(fsApply(fr2_3,each_col,mean))
FSC4<- data.frame(fsApply(fr2_4,each_col,mean))
FSC5<- data.frame(fsApply(fr2_5,each_col,mean))
FSC6<- data.frame(fsApply(fr2_6,each_col,mean))
FSC7<- data.frame(fsApply(fr2_7,each_col,mean))
FSC8<- data.frame(fsApply(fr2_8,each_col,mean))
FSC9<- data.frame(fsApply(fr2_9,each_col,mean))
FSC10<- data.frame(fsApply(fr2_10,each_col,mean))
FSC11<- data.frame(fsApply(fr2_11,each_col,mean))
FSC12<- data.frame(fsApply(fr2_12,each_col,mean))
FSC13<- data.frame(fsApply(fr2_13,each_col,mean))
FSC1$well <- paste0(sampleNames3, "_1")
FSC2$well <- paste0(sampleNames3, "_2")
FSC3$well <- paste0(sampleNames3, "_3")
FSC4$well <- paste0(sampleNames3, "_4")
FSC5$well <- paste0(sampleNames3, "_5")
FSC6$well <- paste0(sampleNames3, "_6")
FSC7$well <- paste0(sampleNames3, "_7")
FSC8$well <- paste0(sampleNames3, "_8")
FSC9$well <- paste0(sampleNames3, "_9")
FSC10$well <- paste0(sampleNames4, "_10")
FSC11$well <- paste0(sampleNames4, "_11")
FSC12$well <- paste0(sampleNames5, "_12")
FSC13$well <- paste0(sampleNames6, "_13")
Fscatter<- rbind(FSC1,FSC2,FSC3,FSC4,FSC5, FSC6, FSC7, FSC8, FSC9, FSC10, FSC11, FSC12, FSC13) %>% separate(well,c("Species", "Media", "Temp", "Rep", "Time_point"))
Fscatter$Temp <- factor(Fscatter$Temp, levels = c("cold", "med", "hot"))
fscat2 <- Fscatter %>%
mutate(hours = case_when(
Time_point == 1 ~ 0,
Time_point == 2 ~ 1,
Time_point == 3 ~ 2,
Time_point == 4 ~ 3,
Time_point == 5 ~ 4,
Time_point == 6 ~ 5,
Time_point == 7 ~ 6,
Time_point == 8 ~ 7,
Time_point == 9 ~ 8,
Time_point == 10 ~ 9,
Time_point == 11 ~ 12,
Time_point == 12 ~ 24,
Time_point == 13 ~ 24
))
write.csv(fscat2, "fscatter_rnd2.csv", row.names = FALSE)fscat2<- read.csv("fscatter_rnd2.csv")
fscat2$Temp <- factor(fscat2$Temp, levels = c("cold", "med", "hot"))
fscat2 %>% dplyr::filter(!(Temp == "cold" &hours ==2)) %>% dplyr::filter(!(Temp == "hot" & hours ==1)) %>% ggplot(aes(x = hours, y = FSC.A, color = Temp, shape = Media)) + facet_grid(Temp~Species) + geom_point(size =2, alpha = .7) + theme_light() + scale_color_manual(values = c("#2D55CC", "#FAAB33", "#FF4041")) +
theme(
strip.text.x = element_text(
size = 12, color = "white", face = "bold.italic"
)) + geom_line(aes(linetype = Media))
fscat2 %>% dplyr::filter(!(Temp == "cold" &hours ==2)) %>% dplyr::filter(!(Temp == "hot" & hours ==2)) %>% group_by(hours, Media, Species, Temp) %>%
dplyr::summarize(mean_size = mean(FSC.A), sd = sd(FSC.A)) %>% ggplot(aes(x = hours, y = mean_size, color = Temp, shape = Media)) + facet_grid(Temp~Species) + geom_point(size =2, alpha = .7) + theme_bw() + scale_color_manual(values = c("#2D55CC", "#FAAB33", "#FF4041")) +geom_line(aes(linetype = Media))+
geom_errorbar(aes(ymin=mean_size-sd, ymax=mean_size+sd), width=0) + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
fscatt1$trial <- "1"
fscat2$trial <- "2"
fscatt1_2 <- rbind(fscatt1, fscat2) %>% filter(hours == 24) %>% filter(Temp == "hot") %>% filter(Species == "Amac")
df_mean <- data.frame(fscatt1_2 %>%
group_by(Media, trial) %>%
dplyr::summarize(meanFSc = mean(FSC.A), sd = sd(FSC.A)) %>% ungroup())
df_sd <- data.frame(fscatt1_2 %>%
group_by(Media, trial) %>%
dplyr::summarize(meanFSc = mean(FSC.A), sd = sd(FSC.A)) %>% ungroup() %>% mutate(ymin = meanFSc-sd) %>% mutate(ymax =meanFSc+sd))
df_points <- fscatt1_2
ggplot() + geom_point(data = df_mean, aes(x = Media, y = meanFSc, color = Media), size =10, shape ="-") + theme_bw() + scale_color_manual(values = c( "black", "#6495ed"))+ geom_errorbar(data = df_sd, aes(x = Media, y = meanFSc, ymin=ymin, ymax=ymax, color = Media), width=0) + geom_jitter(data =df_points, aes(x =Media, y = FSC.A, color = Media), shape = 1, width = 0.114, fill = "white") +facet_grid(~trial) + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) + ggtitle("Amac Cell Size at 24 hours")
Cell size stats
Trial 1
amac_cellsize_full <- df_points %>% dplyr::filter(trial == 1) %>% dplyr::filter(Media == "full")
amac_cellsize_noB12 <- df_points %>% dplyr::filter(trial == 1) %>% dplyr::filter(Media == "noB12")
t.test(amac_cellsize_full$FSC.A, amac_cellsize_noB12$FSC.A)##
## Welch Two Sample t-test
##
## data: amac_cellsize_full$FSC.A and amac_cellsize_noB12$FSC.A
## t = -5.1522, df = 3.1303, p-value = 0.01276
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -16713.191 -4134.111
## sample estimates:
## mean of x mean of y
## 24233.84 34657.49Trial 2
amac_cellsize_full <- df_points %>% dplyr::filter(trial == 2) %>% dplyr::filter(Media == "full")
amac_cellsize_noB12 <- df_points %>% dplyr::filter(trial == 2) %>% dplyr::filter(Media == "noB12")
t.test(amac_cellsize_full$FSC.A, amac_cellsize_noB12$FSC.A)##
## Welch Two Sample t-test
##
## data: amac_cellsize_full$FSC.A and amac_cellsize_noB12$FSC.A
## t = -8.1546, df = 3.5762, p-value = 0.001959
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -9108.574 -4316.035
## sample estimates:
## mean of x mean of y
## 24639.44 31351.75B12 Measurements
Standard Curve
stdcurve <- read.csv("cobalamin_standard_curve.csv")calcurve <-ggplot(stdcurve, aes(x=cynanocobalamin_ug_per_L, y=Total_Fragment_Area)) +geom_smooth(method = "lm") +geom_point(shape = 21) +theme_bw()+ theme(panel.grid.minor = element_blank(), panel.grid.major = element_blank()) + ylab("Total fragment area") + xlab("Cyanocobalamin (µg L-1)")
calcurve
zoom<- ggplot(stdcurve, aes(x=cynanocobalamin_ug_per_L, y=Total_Fragment_Area)) +geom_smooth(method = "lm") +geom_point(shape = 21)+ ylim(c(0,1e+09)) + xlim(c(0,40)) +theme_bw()+ theme(panel.grid.minor = element_blank(), panel.grid.major = element_blank())
zoom
dbl_zoom<- ggplot(stdcurve, aes(x=cynanocobalamin_ug_per_L, y=Total_Fragment_Area)) +geom_smooth(method = "lm") +geom_point(shape = 21)+ ylim(c(0,6e+07)) + xlim(c(0,4)) +theme_bw()+ theme(panel.grid.minor = element_blank(), panel.grid.major = element_blank()) +ylab("") +xlab("")
dbl_zoom
linear<- lm(Total_Fragment_Area ~ cynanocobalamin_ug_per_L, data=stdcurve)
summary(linear)##
## Call:
## lm(formula = Total_Fragment_Area ~ cynanocobalamin_ug_per_L,
## data = stdcurve)
##
## Residuals:
## Min 1Q Median 3Q Max
## -173050484 -55431536 33687621 45363866 149568464
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -44436300 26055431 -1.705 0.112
## cynanocobalamin_ug_per_L 19612170 144924 135.327 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 87720000 on 13 degrees of freedom
## Multiple R-squared: 0.9993, Adjusted R-squared: 0.9992
## F-statistic: 1.831e+04 on 1 and 13 DF, p-value: < 2.2e-16with_inset <- calcurve+
annotation_custom(grob=ggplotGrob(dbl_zoom),
ymin = 4.5e+09, ymax=8e+09, xmin=-Inf, xmax=200)
with_inset
ggsave(“calcurve.pdf”, width = 6, height = 4)
final cell counts
To get number of cells at time of sacrifice to which to normalize b12 mas spec measurements
write.csv(filtered_ps9, "vols.csv", row.names = FALSE)vols<- read.csv("vols.csv")
vols$vol_ml <- as.numeric(vols$vol) /1000 /1000
vols$cells_per_ml <- vols$Count / vols$vol_ml
vols <- vols %>% separate(well,c("Species", "Media", "Temp", "Rep", "Time_point") )
vols$Temp <- factor(vols$Temp, levels = c("cold", "med", "hot"))r4<- vols %>%
mutate(hours = case_when(
Time_point == 10 ~ 9)) %>% filter(Species == "Rpom") %>% filter(Media == "noB12") %>% filter(Temp != "cold")
r4$trial <- "2"
r4$Tube <- c(19,20,21,31,32,33)write.csv(filtered_ps12, "cold_final.csv", row.names = FALSE)results <- read.csv("cold_final.csv")
results$vol_ml <- as.numeric(results$vol) /1000 /1000 #B12 data is mg/ml so getting cells /ml instedd of cells/ml
results$cells_per_ml <- results$Count / results$vol_ml
results <- results %>% separate(well,c("Species", "Media", "Temp", "Rep", "Time_point") )
results$Temp <- factor(results$Temp, levels = c("cold", "med", "hot"))r5<- results %>%
mutate(hours = case_when(
Time_point == 13 ~ 24)) %>% filter(Species == "Rpom") %>% filter(Media == "noB12") %>% filter(Temp == "cold")
r5$trial <- "2"
r5$Tube <- c(7,8,9)
df1 <- rbind(data.frame(r4 %>% group_by(hours, Media, Species, Temp) %>%
dplyr::summarize(mean_cells_per_ml = mean(cells_per_ml), sd = sd(cells_per_ml)) %>% ungroup()), data.frame(r5 %>% group_by(hours, Media, Species, Temp) %>%
dplyr::summarize(mean_cells_per_ml = mean(cells_per_ml), sd = sd(cells_per_ml)) %>% ungroup()))
df2 <- rbind(data.frame(r4 %>% group_by(hours, Media, Species, Temp) %>%
dplyr::summarize(mean_cells_per_ml = mean(cells_per_ml), sd = sd(cells_per_ml)) %>% ungroup() %>% mutate(ymin = mean_cells_per_ml-sd) %>% mutate(ymax = mean_cells_per_ml+sd)), data.frame(r5 %>% group_by(hours, Media, Species, Temp) %>%
dplyr::summarize(mean_cells_per_ml = mean(cells_per_ml), sd = sd(cells_per_ml)) %>% ungroup() %>% mutate(ymin = mean_cells_per_ml-sd) %>% mutate(ymax = mean_cells_per_ml+sd)))
df3 <- rbind(r4, r5)Check that final cell counts look right:
ggplot() + geom_point(data = df1, aes(x = Temp, y = mean_cells_per_ml, color = Temp), size =10, shape ="-") + theme_light() + scale_color_manual(values = c("#2D55CC", "#FAAB33", "#FF4041")) +
theme(
strip.text.x = element_text(
size = 12, color = "white", face = "bold.italic"
)) +
geom_errorbar(data = df2, aes(x = Temp, y = mean_cells_per_ml, ymin=ymin, ymax=ymax, color = Temp), width=0) + geom_jitter(data =df3, aes(x = Temp, y = cells_per_ml, color = Temp), shape = 1, width = 0.115)
print(df3 %>% select(Species, Media, Temp, Rep, Tube, cells_per_ml))## Species Media Temp Rep Tube cells_per_ml
## 1 Rpom noB12 med 1 19 24407279
## 2 Rpom noB12 med 2 20 26743498
## 3 Rpom noB12 med 3 21 23723542
## 4 Rpom noB12 hot 1 31 38518822
## 5 Rpom noB12 hot 2 32 38677891
## 6 Rpom noB12 hot 3 33 38907743
## 7 Rpom noB12 cold 1 7 27917922
## 8 Rpom noB12 cold 2 8 29341732
## 9 Rpom noB12 cold 3 9 28842938These cell count values were added to the spreadsheet with B12 results.
plot B12 concentrations normalized to cell counts
cob<- read.csv("cobalamin_supernatant_pellets_t2.csv")
str(cob)## 'data.frame': 18 obs. of 9 variables:
## $ Sample_name : chr "7P_8_9 " "8P_8_9 " "9P_8_8 " "19P_8_8 " ...
## $ conc_mgL_in_final_extract: int 20 3 0 6 57 36 1 100 0 0 ...
## $ S_P : chr "P" "P" "P" "P" ...
## $ Tube : int 7 8 9 19 21 31 32 33 19 20 ...
## $ Temp : chr "cool" "cool" "cool" "warm" ...
## $ Volume_mL : int 30 30 30 30 30 30 30 30 25 25 ...
## $ cells.perml : int 29341732 28842938 27917922 23723542 26743498 38518822 38677891 38907743 23723542 24407279 ...
## $ Hours_from_Innoc : int 24 24 24 9 9 9 9 9 9 9 ...
## $ Treatment : chr "noB12" "noB12" "noB12" "noB12" ...normalize to mls and # of cells / ml
cob$conc_mg_ml_in_extract <- cob$conc_mgL_in_final_extract / 1000
cob$mg_in_sample <- cob$conc_mg_ml_in_extract * cob$Volume_mL
cob$cells_per_sample <- cob$Volume_mL * cob$cells.perml
cob$ug_per_cell <- cob$mg_in_sample / cob$cells_per_sample * 1000cob$Temp <- factor(cob$Temp, levels = c("cool", "warm", "hot"))
cob %>% ggplot(aes(x=S_P, y=ug_per_cell)) + geom_boxplot(aes(color = Temp)) + scale_color_manual(values = c("#2D55CC", "#FAAB33", "#FF4041")) + theme_light() +
theme(
strip.text.x = element_text(
size = 12, color = "white", face = "bold.italic"
)) +
stat_compare_means(label.y = c(3e-06), method = "t.test", group.by = "S_P")
convenient because ggplot separates treatments, but boxplots are not appropriate for n=3 …
df4 <- data.frame(cob %>% group_by(S_P, Temp) %>%
dplyr::summarize(mean_ug_per_cell = mean(ug_per_cell), sd = sd(ug_per_cell)) %>% ungroup())
df5 <- data.frame(cob %>% group_by(S_P, Temp) %>%
dplyr::summarize(mean_ug_per_cell = mean(ug_per_cell), sd = sd(ug_per_cell)) %>% ungroup()) %>% mutate(ymin = mean_ug_per_cell-sd) %>% mutate(ymax = mean_ug_per_cell+sd)
ggplot() + geom_point(data = df4, aes(x = Temp, y = mean_ug_per_cell, color = Temp), size =10, shape ="-") + theme_bw() + scale_color_manual(values = c("#2D55CC", "#FAAB33", "#FF4041")) + geom_errorbar(data = df5, aes(x = Temp, y = mean_ug_per_cell, ymin=ymin, ymax=ymax, color = Temp), width=0) + geom_jitter(data =cob, aes(x =Temp, y = ug_per_cell, color = Temp), shape = 1, width = 0.114, fill = "white") +facet_grid(~S_P) +ylim( -1.124685e-07, 3e-06) + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
ggplot() + geom_point(data = df4, aes(x = Temp, y = mean_ug_per_cell, color = Temp), size =12, shape ="-") + theme_light() + scale_color_manual(values = c("#2D55CC", "#FAAB33", "#FF4041")) + geom_errorbar(data = df5, aes(x = Temp, y = mean_ug_per_cell, ymin=ymin, ymax=ymax, color = Temp), width=0) + geom_dotplot(data =cob, aes(x =Temp, y = ug_per_cell, color = Temp), shape = 1, dotsize=0.5, binaxis='y', stackdir = 'center', fill = "white" ) +facet_grid(~S_P) +ylim( -1.124685e-07, 3e-06) + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
B12 comparison stats
Two-Way ANOVA - Temperature and Supernatant/Pellet
cob$TempPellet <- paste(cob$Temp, cob$S_P)
anova_1 <- aov(ug_per_cell~Temp+S_P, data = cob)
summary(anova_1)## Df Sum Sq Mean Sq F value Pr(>F)
## Temp 2 6.760e-13 3.378e-13 0.750 0.4904
## S_P 1 2.949e-12 2.949e-12 6.547 0.0227 *
## Residuals 14 6.306e-12 4.504e-13
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Cool v. Warm:
cob_coolP <- cob %>% dplyr::filter(Temp == "cool") %>% dplyr::filter(S_P == "P")
cob_warmP <- cob%>% dplyr::filter(Temp == "warm") %>% dplyr::filter(S_P == "P")
t.test(cob_coolP$ug_per_cell,cob_warmP$ug_per_cell)##
## Welch Two Sample t-test
##
## data: cob_coolP$ug_per_cell and cob_warmP$ug_per_cell
## t = -1.2336, df = 2.5412, p-value = 0.3192
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -2.903646e-06 1.400859e-06
## sample estimates:
## mean of x mean of y
## 2.618782e-07 1.013272e-06cool v. hot:
cob_hotP <- cob %>% dplyr::filter(Temp == "hot") %>% dplyr::filter(S_P == "P")
t.test(cob_coolP$ug_per_cell, cob_hotP$ug_per_cell)##
## Welch Two Sample t-test
##
## data: cob_coolP$ug_per_cell and cob_hotP$ug_per_cell
## t = -1.1822, df = 2.3223, p-value = 0.344
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -3.839299e-06 2.009292e-06
## sample estimates:
## mean of x mean of y
## 2.618782e-07 1.176882e-06Supernatant v. Pellet:
cob_S <- cob%>% dplyr::filter(S_P == "S")
cob_P <- cob %>% dplyr::filter(S_P == "P")
t.test(cob_P$ug_per_cell, cob_S$ug_per_cell)##
## Welch Two Sample t-test
##
## data: cob_P$ug_per_cell and cob_S$ug_per_cell
## t = 2.5997, df = 8.0044, p-value = 0.03162
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## 9.151903e-08 1.527504e-06
## sample estimates:
## mean of x mean of y
## 8.173438e-07 7.832202e-09Session Info
sessionInfo()## R version 4.2.1 (2022-06-23)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Big Sur ... 10.16
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRlapack.dylib
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] patchwork_1.1.2 magick_2.7.4 kableExtra_1.3.4
## [4] flowStats_4.8.2 ggpubr_0.6.0 growthrates_0.8.4
## [7] deSolve_1.35 lattice_0.21-8 scales_1.2.1
## [10] lubridate_1.9.2 forcats_1.0.0 stringr_1.5.0
## [13] dplyr_1.1.2 purrr_1.0.1 readr_2.1.4
## [16] tidyr_1.3.0 tibble_3.2.1 tidyverse_2.0.0
## [19] ggcyto_1.24.1 flowWorkspace_4.8.0 ncdfFlow_2.42.1
## [22] BH_1.81.0-1 RcppArmadillo_0.12.4.0.0 ggplot2_3.4.2
## [25] flowCore_2.8.0
##
## loaded via a namespace (and not attached):
## [1] minqa_1.2.5 colorspace_2.1-0 ggsignif_0.6.4
## [4] deldir_1.0-9 mclust_6.0.0 cytolib_2.8.0
## [7] corpcor_1.6.10 base64enc_0.1-3 rstudioapi_0.14
## [10] farver_2.1.1 hexbin_1.28.3 IDPmisc_1.1.20
## [13] fansi_1.0.4 mvtnorm_1.2-2 xml2_1.3.4
## [16] codetools_0.2-19 splines_4.2.1 mnormt_2.1.1
## [19] cachem_1.0.8 rootSolve_1.8.2.3 robustbase_0.95-1
## [22] knitr_1.43 jsonlite_1.8.5 broom_1.0.5
## [25] cluster_2.1.4 png_0.1-8 graph_1.74.0
## [28] rrcov_1.7-3 compiler_4.2.1 httr_1.4.6
## [31] backports_1.4.1 rainbow_3.7 Matrix_1.5-4.1
## [34] fastmap_1.1.1 cli_3.6.1 htmltools_0.5.5
## [37] tools_4.2.1 coda_0.19-4 gtable_0.3.3
## [40] glue_1.6.2 Rcpp_1.0.10 carData_3.0-5
## [43] Biobase_2.56.0 jquerylib_0.1.4 vctrs_0.6.2
## [46] nlme_3.1-162 svglite_2.1.1 xfun_0.39
## [49] rvest_1.0.3 timechange_0.2.0 lifecycle_1.0.3
## [52] rstatix_0.7.2 XML_3.99-0.14 DEoptimR_1.0-14
## [55] zlibbioc_1.42.0 MASS_7.3-60 RProtoBufLib_2.8.0
## [58] hms_1.1.3 parallel_4.2.1 RColorBrewer_1.1-3
## [61] yaml_2.3.7 curl_5.0.1 aws.signature_0.6.0
## [64] gridExtra_2.3 sass_0.4.6 latticeExtra_0.6-30
## [67] stringi_1.7.12 highr_0.10 S4Vectors_0.34.0
## [70] pcaPP_2.0-3 BiocGenerics_0.42.0 hdrcde_3.4
## [73] FME_1.3.6.2 flowViz_1.60.2 systemfonts_1.0.4
## [76] rlang_1.1.1 pkgconfig_2.0.3 matrixStats_1.0.0
## [79] bitops_1.0-7 pracma_2.4.2 evaluate_0.21
## [82] fda_6.1.4 labeling_0.4.2 ks_1.14.0
## [85] tidyselect_1.2.0 plyr_1.8.8 magrittr_2.0.3
## [88] R6_2.5.1 generics_0.1.3 DBI_1.1.3
## [91] mgcv_1.8-42 pillar_1.9.0 withr_2.5.0
## [94] abind_1.4-5 RCurl_1.98-1.12 car_3.1-2
## [97] interp_1.1-4 KernSmooth_2.23-21 utf8_1.2.3
## [100] tzdb_0.4.0 rmarkdown_2.22 fds_1.8
## [103] aws.s3_0.3.21 jpeg_0.1-10 grid_4.2.1
## [106] minpack.lm_1.2-3 data.table_1.14.8 Rgraphviz_2.40.0
## [109] webshot_0.5.4 digest_0.6.31 RcppParallel_5.1.7
## [112] stats4_4.2.1 munsell_0.5.0 viridisLite_0.4.2
## [115] bslib_0.5.0