360th Post

pH Female Broodstock Histology Analysis

Load libraries

library(ggplot2)
library(ggpubr)
library(rstatix)

## 
## Attaching package: 'rstatix'

## The following object is masked from 'package:stats':
## 
##     filter

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, union

library(broom)
library(lme4)

## Loading required package: Matrix

## Registered S3 methods overwritten by 'lme4':
##   method                          from
##   cooks.distance.influence.merMod car 
##   influence.merMod                car 
##   dfbeta.influence.merMod         car 
##   dfbetas.influence.merMod        car

library(kableExtra)

## 
## Attaching package: 'kableExtra'

## The following object is masked from 'package:dplyr':
## 
##     group_rows

library(tidyr)

## 
## Attaching package: 'tidyr'

## The following objects are masked from 'package:Matrix':
## 
##     expand, pack, unpack

Read in data

female_data <- read.csv("../../data/histology/Female_Gonad/20200820_Female_Gonad_RESULTS.csv", stringsAsFactors = F,colClasses = c(rep("character",4),rep("numeric",7)))


egg_size <- read.csv("../../data/histology/Female_Gonad/egg_size.csv", stringsAsFactors = F, colClasses = c("character","character","numeric"))

format egg size df

#rename columns
colnames(egg_size) <- c("date", "sample_ID", "egg_size")

#merge egg size data with treatment info
egg_size_m <- merge(female_data[,c("date", "sample_ID","pH")], egg_size)

egg_size_m$date <- gsub("20190123", "72",egg_size_m$date)

egg_size_m$date <- gsub("20190221", "93 + 8 day recovery",egg_size_m$date)

plot egg size

#first check distribution

p <- ggplot(data = egg_size_m,aes(x = egg_size, color = pH, group = sample_ID)) + geom_density() + theme_bw() + labs(x = "density", y =expression(paste(log[2]," egg size (",mu,"m"^2,")", sep = "")))+ theme(plot.margin = unit(c(rep(1,4)), "lines"), axis.text.x = element_text(angle = 45)) + facet_wrap(~date)

#it's left skewed so try log transformation
q <- ggplot(data = egg_size_m,aes(x = log(egg_size,2), color = pH, group = sample_ID)) + geom_density() + theme_bw() + labs(x = "density", y = expression(paste(log[2]," egg size (",mu,"m"^2,")", sep = "")))+ theme(plot.margin = unit(c(rep(1,4)), "lines"), axis.text.x = element_text(angle = 45)) + facet_wrap(~date)

#plot boxplots
r <- ggplot(data = egg_size_m,aes(x = date, y = egg_size, color = pH, group = interaction(pH,date, sample_ID))) + geom_boxplot(outlier.shape = NA) + geom_point(pch = 16, size=0.2,position = position_jitterdodge(jitter.width = 0.05)) + theme_bw() + labs(x = "exposure time (days)", y = expression(paste(log[2]," egg size (",mu,"m"^2,")", sep = "")))+ theme(plot.margin = unit(c(rep(1,4)), "lines"))

#after log transformation
s <- ggplot(data = egg_size_m,aes(x = date, y = log(egg_size,2), color = pH, group = interaction(pH,date, sample_ID))) + geom_boxplot(outlier.shape = NA) + geom_point(pch = 16, size=0.2,position = position_jitterdodge(jitter.width = 0.05)) + theme_bw() + labs(x = "exposure time (days)", y = expression(paste(log[2]," egg size (",mu,"m"^2,")", sep = "")))+ theme(plot.margin = unit(c(rep(1,4)), "lines"))

#plot egg size without transformation
ggpubr::ggarrange(p,r,labels = "AUTO", common.legend = T)

#plot eggs size after transformation
ggpubr::ggarrange(q,s,labels = "AUTO", common.legend = T)

Save plots

jpeg("egg_size_dist.jpg", width = 12, height = 6, units = "in", res = 300)
ggpubr::ggarrange(p,r,labels = "AUTO", common.legend = T)
dev.off()

## quartz_off_screen 
##                 2

jpeg("log_egg_size_dist.jpg", width = 12, height = 6, units = "in", res = 300)
ggpubr::ggarrange(q,s,labels = "AUTO", common.legend = T)
dev.off()

## quartz_off_screen 
##                 2

Egg size stats

tidy(TukeyHSD(aov(log(egg_size,2) ~ pH * date, data = egg_size_m))) %>%
  kbl(caption = "2-way ANOVA: log2 egg size ~ pH * time") %>%
  kable_classic(full_width = F, html_font = "Cambria")

2-way ANOVA: log2 egg size ~ pH * time
term	contrast	estimate	conf.low	conf.high	adj.p.value
pH	low-ambient	-0.2422256	-0.3107394	-0.1737119	0.0000000
date	93 + 8 day recovery-72	0.4998445	0.4350260	0.5646630	0.0000000
pH:date	low:72-ambient:72	-0.5010485	-0.6330825	-0.3690145	0.0000000
pH:date	ambient:93 + 8 day recovery-ambient:72	0.3190786	0.2147047	0.4234524	0.0000000
pH:date	low:93 + 8 day recovery-ambient:72	0.3643507	0.2344607	0.4942408	0.0000000
pH:date	ambient:93 + 8 day recovery-low:72	0.8201271	0.6948747	0.9453795	0.0000000
pH:date	low:93 + 8 day recovery-low:72	0.8653993	0.7182066	1.0125920	0.0000000
pH:date	low:93 + 8 day recovery-ambient:93 + 8 day recovery	0.0452722	-0.0777181	0.1682624	0.7799909

Save 2-way ANOVA results

write.table(tidy(TukeyHSD(aov(log(egg_size,2) ~ pH * date, data = egg_size_m))),"2wayAOV_Tuk_egg_size.txt",sep = "\t", quote = F, row.names = F)

Tissue proportions analysis

#calculate proportions

female_data$prop_egg <- (female_data$egg.area..based.on.ind.eggs. / female_data$total.area..um.) *100

female_data$prop_no_tissue <- (female_data$no.tissue.area / female_data$total.area..um.) *100

female_data$prop_cnx_tissue_indve <- (female_data$connective.tissue.area..egg.based.on.indiv./ female_data$total.area..um.)*100

#create column for relative # of eggs
female_data$number.of.eggs.per.um2 <- (female_data$number.of.eggs/female_data$total.area..um.)

female_data$pH_date <- paste(female_data$date, female_data$pH, sep = "_")

#convert female data to long format

female_data_STACKED <- tidyr::gather(female_data,"metric", "um", 5:(ncol(female_data)-1))

female_data_STACKED$um <- as.numeric(female_data_STACKED$um)




#reformat data and calculate mean and sd proportions
female_data_STACKED_props <- female_data_STACKED[grep("prop", female_data_STACKED$metric),]%>%group_by(date, pH,metric) %>% summarise(mean = mean(um), sd =sd(um)) %>% mutate(y_pos = cumsum(mean))

## `summarise()` regrouping output by 'date', 'pH' (override with `.groups` argument)

#convert metric to ordered factor
female_data_STACKED_props$metric <- factor(female_data_STACKED_props$metric, levels = rev(c("prop_cnx_tissue_indve","prop_egg", "prop_no_tissue")))


#chi sq to compare 20190123 amb vs. low
ct <- data.frame(female_data_STACKED_props[grep("20190123", female_data_STACKED_props$date),c("pH","mean","metric")]%>% pivot_wider(names_from = metric, values_from = mean))

rownames(ct) <- ct$pH

ct <- ct[,-1]

ct_a <- chisq.test(ct)
ct_a <- tidy(ct_a)
ct_a$comparison <- "ambient:72-low:72"


#chi sq to compare 20190123 amb vs. 20190221 amb
ct <- data.frame(female_data_STACKED_props[grep("ambient", female_data_STACKED_props$pH),c("date","mean","metric")]%>% pivot_wider(names_from = metric, values_from = mean))

rownames(ct) <- ct$pH

ct <- ct[,-1]

ct_b <- chisq.test(ct)
ct_b <- tidy(ct_b)
ct_b$comparison <- "ambient:72-ambient:93 + 8 day recovery"

#chi sq to compare 20190123 low vs. 20190221 low
ct <- data.frame(female_data_STACKED_props[grep("low", female_data_STACKED_props$pH),c("date","mean","metric")]%>% pivot_wider(names_from = metric, values_from = mean))

rownames(ct) <- ct$pH

ct <- ct[,-1]

ct_c <- chisq.test(ct)
ct_c <- tidy(ct_c)
ct_c$comparison <- "low:72-low:93 + 8 day recovery"



#chi sq to compare 20190221 amb vs. 20190221 low
ct <- data.frame(female_data_STACKED_props[grep("20190221", female_data_STACKED_props$date),c("pH","mean","metric")]%>% pivot_wider(names_from = metric, values_from = mean))

rownames(ct) <- ct$pH

ct <- ct[,-1]

ct_d <- chisq.test(ct)
ct_d <- tidy(ct_d)
ct_d$comparison <- "ambient:93 + 8 day recovery-low:93 + 8 day recovery"


#chi sq to compare 20190123 amb vs. 20190221 low
ct <- data.frame(female_data_STACKED_props[which(female_data_STACKED_props$date=="20190123" & female_data_STACKED_props$pH =="ambient" | female_data_STACKED_props$date=="20190221" & female_data_STACKED_props$pH =="low"),c("pH","mean","metric")]%>% pivot_wider(names_from = metric, values_from = mean))

rownames(ct) <- ct$pH

ct <- ct[,-1]

ct_e <- chisq.test(ct)
ct_e <- tidy(ct_e)
ct_e$comparison <- "ambient:72-low:93 + 8 day recovery"


#chi sq to compare 20190221 amb vs. 20190123 low
ct <- data.frame(female_data_STACKED_props[which(female_data_STACKED_props$date=="20190221" & female_data_STACKED_props$pH =="ambient" | female_data_STACKED_props$date=="20190123" & female_data_STACKED_props$pH =="low"),c("pH","mean","metric")]%>% pivot_wider(names_from = metric, values_from = mean))

rownames(ct) <- ct$pH

ct <- ct[,-1]

ct_f <- chisq.test(ct)
ct_f <- tidy(ct_f)
ct_f$comparison <- "ambient:93 + 8 day recovery-low:72"


ct <- rbind(ct_a, ct_b, ct_c, ct_d, ct_e, ct_f)

#perform FDR correction for multiple comparisons

ct$adj.p.value <- p.adjust(ct$p.value)

#reorder columns

ct <- ct[,c("comparison", "statistic", "p.value", "adj.p.value")]


ct %>%
  kbl(caption = "Chi square test: tissue proportions ~ pH") %>%
  kable_classic(full_width = F, html_font = "Cambria")

Chi square test: tissue proportions ~ pH
comparison	statistic	p.value	adj.p.value
ambient:72-low:72	2.9365981	0.2303169	0.9212677
ambient:72-ambient:93 + 8 day recovery	2.8983530	0.2347635	0.9212677
low:72-low:93 + 8 day recovery	5.8072777	0.0548234	0.3289402
ambient:93 + 8 day recovery-low:93 + 8 day recovery	5.6465427	0.0594113	0.3289402
ambient:72-low:93 + 8 day recovery	0.5131967	0.7736789	1.0000000
ambient:93 + 8 day recovery-low:72	0.7966350	0.6714488	1.0000000

pd <- position_dodge2(width = 0.2)

ggplot(female_data_STACKED_props, aes(pH,mean, fill = pH,alpha = metric)) + 
geom_bar(stat = "identity") + 
geom_errorbar(aes(ymax = y_pos + sd, ymin=y_pos - sd), stat = "identity", width = 0.1, alpha = 0.7, position = pd) + facet_wrap(~date)+ scale_alpha_manual(values=c(seq(0.3,1, length.out = 3))) + theme_bw() + scale_y_continuous(breaks = seq(from = 0, to = 100, by = 10)) + ylab("mean proportion tissue area") + theme(axis.title.x = element_blank(), axis.ticks.x = element_blank(), axis.text.x = element_blank())

Write out Chi sq. test results

write.table(ct,"chisq_tissue_prop.txt",sep = "\t", quote = F, row.names = F)

Generate stacked bars showing average tissue proportions

#export plot with significance denoted
#plot stacked bar of tissue areas



#change labels for plotting
female_data_STACKED_props_labs <- female_data_STACKED_props

female_data_STACKED_props_labs$date <- gsub("20190123", "72",female_data_STACKED_props_labs$date)

female_data_STACKED_props_labs$date <- gsub("20190221", "93 + 8 day recovery",female_data_STACKED_props_labs$date)

female_data_STACKED_props_labs$metric <- gsub("prop_no_tissue","follicle area", female_data_STACKED_props_labs$metric)

female_data_STACKED_props_labs$metric <- gsub("prop_cnx_tissue_indve","connective tissue area", female_data_STACKED_props_labs$metric)

female_data_STACKED_props_labs$metric <- gsub("prop_egg","egg area", female_data_STACKED_props_labs$metric)

#convert metric to ordered factor
female_data_STACKED_props_labs <- female_data_STACKED_props_labs[order(female_data_STACKED_props_labs$metric),]

female_data_STACKED_props_labs$metric <- factor(female_data_STACKED_props_labs$metric, levels = rev(c(rep("connective tissue area",1), rep("egg area",1), rep("follicle area",1))))

pd <- position_dodge2(width = 0.2)

ggplot(female_data_STACKED_props_labs, aes(pH,mean, fill = pH,alpha = metric)) + 
geom_bar(stat = "identity") + 
geom_errorbar(aes(ymax = y_pos + sd, ymin=y_pos - sd), stat = "identity", width = 0.1, alpha = 0.7, position = pd) + facet_wrap(~date)+ scale_alpha_manual(values=c(seq(0.3,1, length.out = 3))) + theme_bw() + scale_y_continuous(breaks = seq(from = 0, to = 100, by = 10)) + ylab("mean proportion of total tissue area (%)") + theme(axis.title.x = element_blank(), axis.ticks.x = element_blank(), axis.text.x = element_blank())

jpeg("tissue_proportions_stacked_bar.jpg", width = 7, height = 6, units = "in", res = 300)
ggplot(female_data_STACKED_props_labs, aes(pH,mean, fill = pH,alpha = metric)) + 
geom_bar(stat = "identity") + 
geom_errorbar(aes(ymax = y_pos + sd, ymin=y_pos - sd), stat = "identity", width = 0.1, alpha = 0.7, position = pd) + facet_wrap(~date)+ scale_alpha_manual(values=c(seq(0.3,1, length.out = 3))) + theme_bw() + scale_y_continuous(breaks = seq(from = 0, to = 100, by = 10)) + ylab("mean proportion of total tissue area (%)") + theme(axis.title.x = element_blank(), axis.ticks.x = element_blank(), axis.text.x = element_blank()) 
dev.off()

## quartz_off_screen 
##                 2

filter data for comparing

female_data_STACKED_filt <- female_data_STACKED[grep("prop|per|mean", female_data_STACKED$metric),]

Convert character to factors

female_data_STACKED_filt$pH <- as.factor(female_data_STACKED_filt$pH)

female_data_STACKED_filt$date <- as.factor(female_data_STACKED_filt$date)

female_data_STACKED_filt$metric <- as.factor(female_data_STACKED_filt$metric)

Plot each metric x treatment group

#plot january samples
#jpeg("20190123_boxplots.jpg", width = 6, height = 6, units = "in", res = 300)
a <- ggplot(data = female_data_STACKED_filt[which(female_data_STACKED_filt$date == "20190123"),],aes(x = pH, y = um, color = pH)) + geom_boxplot(outlier.shape = NA) + facet_wrap(~metric, scale = "free") + geom_jitter(alpha=0.8,shape =16,size = 2, position= position_jitter(0.1)) + theme_bw() + ggtitle("20190123_samples; ambient pH (n = 5) x low pH (n = 3)")
#dev.off()

#plot feb samples
#jpeg("20190221_boxplots.jpg", width = 6, height = 6, units = "in", res = 300)
b <- ggplot(data = female_data_STACKED_filt[which(female_data_STACKED_filt$date == "20190221"),],aes(x = pH, y = um, color = pH)) + geom_boxplot(outlier.shape = NA) + facet_wrap(~metric, scale = "free") + geom_jitter(alpha=0.8,shape =16,size = 2, position= position_jitter(0.1)) + theme_bw() + ggtitle("20190221_samples; ambient pH (n = 8) x low pH (n = 3)")
#dev.off()

jpeg("all_samples_boxplots.jpg", width = 12, height = 6, units = "in", res = 300)
ggpubr::ggarrange(a,b,labels = "AUTO", common.legend = T)
dev.off()

#plot all samples together
#jpeg("dates_combined_boxplots.jpg", width = 6, height = 6, units = "in", res = 300)
#ggplot(data = female_data_STACKED_filt,aes(x = pH, y = um)) +geom_boxplot(aes(color = pH), outlier.shape = NA) + facet_wrap(~metric, scale = "free") + geom_jitter(aes(color = pH),alpha = 0.8,shape =16, size = 2,position= position_jitter(0.2)) + theme_bw() + ggtitle("all samples; ambient pH (n = 13) x low pH (n = 6)")
#dev.off()

clarify pH treatment and dates

female_data_STACKED_filt$pH <- gsub("low", "low (6.8)", female_data_STACKED_filt$pH)

female_data_STACKED_filt$pH <- gsub("ambient", "ambient (7.8)", female_data_STACKED_filt$pH)

female_data_STACKED_filt$date <- gsub("20190123", "72",female_data_STACKED_filt$date)

female_data_STACKED_filt$date <- gsub("20190221", "93 + 8 day recovery",female_data_STACKED_filt$date)

plot egg size only

c <- ggplot(data = female_data_STACKED_filt[which(female_data_STACKED_filt$metric == "mean.egg.size"),],aes(x = date, y = um, color = pH, group = interaction(pH,date))) + geom_boxplot(outlier.shape = NA) + geom_point(pch = 16, position = position_jitterdodge()) + theme_bw() + labs(x = "exposure time (days)", y = expression(paste("mean egg size (",mu,"m"^2,")", sep = "")))+ theme(plot.margin = unit(c(rep(1,4)), "lines"),, axis.text.x = element_text(angle = 45))

d <- ggplot(data = female_data_STACKED_filt[which(female_data_STACKED_filt$metric == "number.of.eggs.per.um2"),],aes(x = date, y = um, color = pH, group = interaction(pH,date))) + geom_boxplot(outlier.shape = NA) + geom_point(pch = 16, position = position_jitterdodge()) + theme_bw() + labs(x = "exposure time (days)", y = expression(paste("number of eggs per ",mu,"m"^2," tissue", sep = "")))+ theme(plot.margin = unit(c(rep(1,4)), "lines"),, axis.text.x = element_text(angle = 45))

e <- ggplot(data = female_data_STACKED_filt[which(female_data_STACKED_filt$metric == "prop_cnx_tissue_indve"),],aes(x = date, y = um, color = pH, group = interaction(pH,date))) + geom_boxplot(outlier.shape = NA) + geom_point(pch = 16, position = position_jitterdodge()) + theme_bw() + labs(x = "exposure time (days)", y = "connective tissue\nproportion (%)")+ theme(plot.margin = unit(c(rep(1,4)), "lines"),, axis.text.x = element_text(angle = 45))

f <- ggplot(data = female_data_STACKED_filt[which(female_data_STACKED_filt$metric == "prop_no_tissue"),],aes(x = date, y = um, color = pH, group = interaction(pH,date))) + geom_boxplot(outlier.shape = NA) + geom_point(pch = 16, position = position_jitterdodge()) + theme_bw() + labs(x = "exposure time (days)", y = "follicle area proportion (%)") + theme(plot.margin = unit(c(rep(1,4)), "lines"), axis.text.x = element_text(angle = 45))


jpeg("all_samples_boxplots_newlabs.jpg", width = 10, height = 3, units = "in", res = 300)

ggpubr::ggarrange(c,d,e,f, labels = "AUTO", vjust = 1,common.legend = T, nrow = 1,ncol = 4,legend = "bottom", font.label = c(size = 18), align = "hv")
dev.off()

Test for treatment and time effect

## egg size test---------------------------------
#egg size normality with shapiro test
female_data_STACKED_filt %>% tidyr::pivot_wider(names_from = metric,values_from = um ) %>% shapiro_test(mean.egg.size)
#  variable      statistic     p
#  <chr>             <dbl> <dbl>
#1 mean.egg.size     0.951 0.407

#levene test
female_data_STACKED_filt %>% tidyr::pivot_wider(names_from = metric,values_from = um ) %>% levene_test(mean.egg.size~pH)
# A tibble: 1 x 4
#    df1   df2 statistic     p
#  <int> <int>     <dbl> <dbl>
#1     1    17  0.000532 0.982



#try a two way ANOVA for effects of pH and time
egg_size_pHxTIME <- aov(um~date*pH, data = female_data_STACKED_filt[which(female_data_STACKED_filt$metric=="mean.egg.size"),])

summary(egg_size_pHxTIME)
#            Df Sum Sq Mean Sq F value Pr(>F)  
#date         1  31.49  31.491   4.941  0.042 *
#pH           1   6.07   6.074   0.953  0.344  
#date:pH      1   4.97   4.974   0.780  0.391  
#Residuals   15  95.60   6.373

#try a 1way anova for time effect
egg_size_TIME  <- aov(um~date, data = female_data_STACKED_filt[which(female_data_STACKED_filt$metric=="mean.egg.size"),])

summary(egg_size_TIME)

#          Df Sum Sq Mean Sq F value Pr(>F)  
#date         1  31.49  31.491    5.02 0.0387 *
#Residuals   17 106.65   6.273                 
#---
#Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


## Date has an effect

#Does pH have an effect?
egg_size_pH  <- aov(um~pH, data = female_data_STACKED_filt[which(female_data_STACKED_filt$metric=="mean.egg.size"),])

summary(egg_size_pH)
#            Df Sum Sq Mean Sq F value Pr(>F)
#pH           1   9.36   9.361   1.236  0.282
#Residuals   17 128.78   7.575     


#is there an effect of treatment at a specific time?
egg_size_jan_pH  <- aov(um~pH, data = female_data_STACKED_filt[which(female_data_STACKED_filt$metric=="mean.egg.size" & female_data_STACKED_filt$date == "20190123"),])

#yes there is 
summary(egg_size_jan_pH)
#            Df Sum Sq Mean Sq F value Pr(>F)  
#pH           1 10.963  10.963   8.825 0.0249 *
#Residuals    6  7.453   1.242                 
#---
#Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

#BUT NOT in the Feb group
#summary(egg_size_feb_pH)
#            Df Sum Sq Mean Sq F value Pr(>F)
#pH           1   0.08   0.085   0.009  0.928
#Residuals    9  88.14   9.794

Run wilcox test to see if follicle area is significantly different

wt <- female_data_STACKED_filt %>% group_by(metric) %>% wilcox_test(um ~ pH )


glm_res <- female_data_STACKED_filt %>% group_by(metric) %>% do(glm_models = glm(um ~ pH*date, data = .))

full.lmer <-  female_data_STACKED_filt %>% group_by(metric) %>% do(lmer_models = lmer(um ~ pH*date + (1|tank), data = ., REML = F))


red.lmer <-  female_data_STACKED_filt %>% group_by(metric) %>% do(lmer_models = lmer(um ~ 1 + (1|tank), data = ., REML = F))

for(i in 1:length(full.lmer$metric)){
  a <- anova(red.lmer$lmer_models[[i]], full.lmer$lmer_models[[i]])
  print(a)
}


aov_res_summ <- broom::tidy(aov_res)

wt <- wt %>% add_xy_position(x = "pH")

bxp <- ggboxplot(female_data_STACKED_filt, x = "pH", y = "um", color = "pH",facet.by = "metric", scales = "free", add = "jitter") 

jpeg("dates_combined_boxplots_pvals.jpg", width = 6, height = 6, units = "in", res = 300)
bxp + stat_pvalue_manual(wt,bracket.nudge.y = -2,label = "{p}") + scale_y_continuous(expand = expansion(mult = c(0.05, 0.1)))
dev.off()


wwt <- female_data_STACKED_filt %>% group_by(date,metric) %>% do(broom::tidy(wilcox.test(um ~ pH,data = . )))



#run stats on all data regardless of date
wt <- female_data_STACKED_filt %>% group_by(metric) %>% do(broom::tidy(wilcox.test(um ~ pH,data = . )))


#run stats on all data regardless of date
Tt <- female_data_STACKED_filt %>% group_by(metric) %>% do(broom::tidy(t.test(um ~ pH,data = . )))

aovt <- female_data_STACKED_filt %>% group_by(metric) %>% do(broom::tidy(aov(um ~ pH,data = . )))

glmt <- female_data_STACKED_filt %>% group_by(metric) %>% do(broom::tidy(glm(um ~ pH,data = . )))

Plot percent follicle area for each tank

ggplot(data = female_data, aes(x = tank,y = perc_follicle_area, group = as.factor(tank), fill = pH)) + geom_violin(trim = FALSE) + geom_boxplot(width = 0.2) + geom_jitter(shape =16, position= position_jitter(0.1)) + ylab("follicle area (%)") + theme_bw()

anova to see if there is a tank or a pH effect

model <- aov(perc_follicle_area ~ pH * tank, data = female_data)
summary(model)

Plot percent egg area for each treatment group

ggplot(data = female_data, aes(x = pH,y = perc_egg_area, group = as.factor(pH), fill = pH)) + geom_violin(trim = FALSE) + geom_boxplot(width = 0.2) + geom_jitter(shape =16, position= position_jitter(0.1)) + ylab("egg area (%)") + theme_bw()

Run wilcox test to see if egg area is significantly different

pH6.8 <- subset(female_data, pH == "6.8", perc_egg_area, drop = TRUE)
pHamb <- subset(female_data, pH == "amb", perc_egg_area, drop = TRUE)
wt <- wilcox.test(pH6.8, pHamb)
print(wt$p.value)

Plot percent egg area for each tank

ggplot(data = female_data, aes(x = tank,y = perc_egg_area, group = as.factor(tank), fill = pH)) + geom_violin(trim = FALSE) + geom_boxplot(width = 0.2) + geom_jitter(shape =16, position= position_jitter(0.1)) + ylab("egg area (%)") + theme_bw()

anova to see if there is a tank or a pH effect

model <- aov(perc_egg_area ~ pH * tank, data = female_data)
summary(model)

Plot egg:follicle ratio for each treatment group

ggplot(data = female_data, aes(x = pH,y = follicle_egg_ratio, group = as.factor(pH), fill = pH)) + geom_violin(trim = FALSE) + geom_boxplot(width = 0.2) + geom_jitter(shape =16, position= position_jitter(0.1)) + ylab("egg:follicle ratio") + theme_bw()

Run wilcox test to see if egg:follicle ratio is significantly different

pH6.8 <- subset(female_data, pH == "6.8", follicle_egg_ratio, drop = TRUE)
pHamb <- subset(female_data, pH == "amb", follicle_egg_ratio, drop = TRUE)
wt <- wilcox.test(pH6.8, pHamb)
print(wt$p.value)

Plot egg:follicle ratio for each tank

ggplot(data = female_data, aes(x = tank,y = follicle_egg_ratio, group = as.factor(tank), fill = pH)) + geom_violin(trim = FALSE) + geom_boxplot(width = 0.2) + geom_jitter(shape =16, position= position_jitter(0.1)) + ylab("egg:follicle ratio") + theme_bw()

anova to see if there is a tank or a pH effect

model <- aov(follicle_egg_ratio ~ pH * tank, data = female_data)
summary(model)

360th Post

pH Female Broodstock Histology Analysis

Shelly Trigg

10/07/2020