library(readxl) library(tidyverse) library(agricolae) library(fitdistrplus) library(DHARMa) library(lme4) library(car) library(emmeans) library(gamlss) #####Setting things up #setting contrasts options(contrasts = c("contr.sum", "contr.poly")) #importing data and omiting blank cells resulting from removing outliers spd <- read_excel("Data/updated stats.xlsx") spd$TRT <- as.factor(spd$TRT) spd$Day <- as.factor(spd$Day) #####KDM6B #looking at distribution plotdist(spd$NKdm6b, histo = TRUE, demp = TRUE) descdist(spd$NKdm6b, discrete = FALSE, boot = 500) #this passes model assumptions after removing extreme outlier on low end kdmlog <- glm(log(NKdm6b)~ TRT * Day, family = gaussian, data = spd) summary(kdmlog) #looking at residuals plot(fitted(kdmlog), residuals(kdmlog), xlab = "Fitted Values", ylab = "Residuals") abline(h = 0, lty = 2) lines(smooth.spline(fitted(kdmlog), residuals(kdmlog))) #checking assumptions - passed with log trans and 3 extreme outliers removed. RANDOMLY DIFFERENT THAN EARLIER TODAY kdmlog_simres <- simulateResiduals(kdmlog) plot(kdmlog_simres) drop1(kdmlog, test = "Chisq") #final results Anova(kdmlog, type = "III") #Analysis of Deviance Table (Type III tests) # #Response: log(NKdm6b) #LR Chisq Df Pr(>Chisq) #TRT 118.938 4 < 2.2e-16 *** # Day 87.800 3 < 2.2e-16 *** # TRT:Day 66.789 12 1.271e-09 *** # --- # Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 #creating emmeans kdmemm <- emmeans(kdmlog, ~TRT:Day, type = "response") kdmemm #creating contrast strings DHW.kdm = c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0) DHWF.kdm = c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0) FPTC.kdm = c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0) LHW.kdm = c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0) MPTC.kdm = c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1) #contrasts with last sampling point with fdr adjustment contrast(kdmemm, method = list("DHW - DHWF" = DHW.kdm - DHWF.kdm, "DHW - FPTC" = DHW.kdm - FPTC.kdm, "DHW - LHW" = DHW.kdm - LHW.kdm, "DHW - MPTC" = DHW.kdm - MPTC.kdm, "DHWF - FPTC" = DHWF.kdm - FPTC.kdm, "DHWF - LHW" = DHWF.kdm - LHW.kdm, "DHWF - MPTC" = DHWF.kdm - MPTC.kdm, "FPTC - LHW" = FPTC.kdm - LHW.kdm, "FPTC - MPTC" = FPTC.kdm - MPTC.kdm, "LHW - MPTC" = LHW.kdm - MPTC.kdm), adjust = "fdr") #contrast ratio SE df null t.ratio p.value #DHW / DHWF 4.9279 1.7675 91 1 4.447 0.0001 #DHW / FPTC 3.5829 1.3478 91 1 3.392 0.0021 #DHW / LHW 2.1683 0.8157 91 1 2.057 0.0531 #DHW / MPTC 0.4408 0.1581 91 1 -2.284 0.0412 #DHWF / FPTC 0.7271 0.2735 91 1 -0.847 0.3991 #DHWF / LHW 0.4400 0.1655 91 1 -2.182 0.0452 #DHWF / MPTC 0.0895 0.0321 91 1 -6.730 <.0001 #FPTC / LHW 0.6052 0.2378 91 1 -1.278 0.2271 #FPTC / MPTC 0.1230 0.0463 91 1 -5.570 <.0001 #LHW / MPTC 0.2033 0.0765 91 1 -4.235 0.0001 # #P value adjustment: fdr method for 10 tests #Tests are performed on the log scale #graphing emms y_kdm <- expression(paste("Normalized ", italic("Kdm6b"), "(+IR) expression ± 95% CI")) kdmgraph<-summary(kdmemm, infer = TRUE, type = "response") kdmgraph kdmgraph %>% ggplot(aes(x = Day, y = response, color = TRT, group = TRT)) + geom_line(stat = "identity", position = position_dodge(.2), linewidth = 1) + geom_errorbar(aes(ymin = lower.CL, ymax = upper.CL), width = .2, position = position_dodge(.2)) + labs(x = "Sampling Day", y = y_kdm, colour = "Treatment") + scale_color_manual(values = c("blue", "orange", "light blue", "dark blue", "red"), limits = c("MPTC", "FPTC", "LHW", "DHW", "DHW+F")) + theme_classic() + theme(axis.title = element_text(size = 13), axis.text = element_text(size = 12)) #zooming in on last sampling point kdmgraph %>% ggplot(aes(x = Day, y = response, color = TRT, group = TRT)) + geom_line(stat = "identity", position = position_dodge(.2), linewidth = 1) + geom_errorbar(aes(ymin = lower.CL, ymax = upper.CL), width = .2, position = position_dodge(.2)) + labs(x = "Sampling Day", y = y_kdm, colour = "Treatment") + scale_color_manual(values = c("blue", "orange", "light blue", "dark blue", "red"), limits = c("MPTC", "FPTC", "LHW", "DHW", "DHW+F")) + coord_cartesian(ylim = c(0, 0.0002)) + theme_classic() #####Dmrt1 #looking at distributions plotdist(spd$NDmrt1, histo = TRUE, demp = TRUE) descdist(spd$NDmrt1, discrete = FALSE, boot = 500) #running glm with lognormal distribution dmrlog <- glm(log(NDmrt1)~ TRT * Day, family = gaussian, data = spd) summary(dmrlog) #looking at residuals plot(fitted(dmrlog), residuals(dmrlog), xlab = "Fitted Values", ylab = "Residuals") abline(h = 0, lty = 2) lines(smooth.spline(fitted(dmrlog), residuals(dmrlog))) #checking assumptions - passed with lognormal distribution dmrlog_simres <- simulateResiduals(dmrlog) plot(dmrlog_simres) drop1(dmrlog, test = "Chisq") #final results Anova(dmrlog, type = "III") #Analysis of Deviance Table (Type III tests) # #Response: log(NDmrt1) #R Chisq Df Pr(>Chisq) #TRT 113.363 4 < 2e-16 *** # Day 5.692 3 0.12760 #TRT:Day 22.636 12 0.03098 * # --- # Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 #creating emmeans dmremm <- emmeans(dmrlog, ~ TRT:Day, type = "response") dmremm #creating contrast strings DHW.dmr = c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0) DHWF.dmr = c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0) FPTC.dmr = c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0) LHW.dmr = c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0) MPTC.dmr = c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1) #contrasts with fdr adjustment contrast(dmremm, method = list("DHW - DHWF" = DHW.dmr - DHWF.dmr, "DHW - FPTC" = DHW.dmr - FPTC.dmr, "DHW - LHW" = DHW.dmr - LHW.dmr, "DHW - MPTC" = DHW.dmr - MPTC.dmr, "DHWF - FPTC" = DHWF.dmr - FPTC.dmr, "DHWF - LHW" = DHWF.dmr - LHW.dmr, "DHWF - MPTC" = DHWF.dmr - MPTC.dmr, "FPTC - LHW" = FPTC.dmr - LHW.dmr, "FPTC - MPTC" = FPTC.dmr - MPTC.dmr, "LHW - MPTC" = LHW.dmr - MPTC.dmr), adjust = "fdr") #contrast ratio SE df null t.ratio p.value #DHW / DHWF 1.5014 0.8373 94 1 0.729 0.4679 #DHW / FPTC 2.2787 1.2707 94 1 1.477 0.2043 #DHW / LHW 0.2487 0.1387 94 1 -2.496 0.0239 #DHW / MPTC 0.1491 0.0831 94 1 -3.413 0.0024 #DHWF / FPTC 1.5177 0.8463 94 1 0.748 0.4679 #DHWF / LHW 0.1656 0.0924 94 1 -3.224 0.0035 #DHWF / MPTC 0.0993 0.0554 94 1 -4.142 0.0004 #FPTC / LHW 0.1091 0.0609 94 1 -3.973 0.0005 #FPTC / MPTC 0.0654 0.0365 94 1 -4.890 <.0001 #LHW / MPTC 0.5996 0.3343 94 1 -0.917 0.4516 # #P value adjustment: fdr method for 10 tests #Tests are performed on the log scale #graphing emmeans y_dmr <- expression(paste("Normalized ", italic("Dmrt1"), " expression ± 95% CI")) dmrgraph<-summary(dmremm, infer = TRUE, type = "response") dmrgraph dmrgraph %>% ggplot(aes(x = Day, y = response, color = TRT, group = TRT)) + geom_line(stat = "identity", position = position_dodge(.2), linewidth = 1) + geom_errorbar(aes(ymin = lower.CL, ymax = upper.CL), width = .2, position = position_dodge(.2)) + labs(x = "Sampling Day", y = y_dmr, colour = "Treatment") + scale_color_manual(values = c("blue", "orange", "light blue", "dark blue", "red"), limits = c("MPTC", "FPTC", "LHW", "DHW", "DHW+F")) + theme_classic() + theme(axis.title = element_text(size = 13), axis.text = element_text(size = 12)) #####Sox9 #looking at distribution plotdist(spd$NSox9, histo = TRUE, demp = TRUE) descdist(spd$NSox9, discrete = FALSE, boot = 500) #running glm soxsqrt <- glm(sqrt(NSox9)~ TRT * Day, family = "gaussian", data = spd) summary(soxsqrt) #looking at residuals plot(fitted(soxsqrt), residuals(soxsqrt), xlab = "Fitted Values", ylab = "Residuals") abline(h = 0, lty = 2) lines(smooth.spline(fitted(soxsqrt), residuals(soxsqrt))) #checking model assumptions - passed using sqrt transformation soxsqrt_simres <- simulateResiduals(soxsqrt) plot(soxsqrt_simres) drop1(soxsqrt, test = "Chisq") #final results Anova(soxsqrt, type = "III") #creating emmeans soxemm <- emmeans(soxsqrt, pairwise ~ Day, type = "response", adjust = "fdr") soxemm #nned to create a different one with TRT included to graph appropriately soxemmgraph <- emmeans(soxsqrt, ~TRT:Day, type = "response") #graphing emmeans y_sox <- expression(paste("Normalized ", italic("Sox9"), " expression ± 95% CI")) soxgraph<-summary(soxemmgraph, infer = TRUE, type = "response") soxgraph soxgraph %>% ggplot(aes(x = Day, y = response, color = TRT, group = TRT)) + geom_line(stat = "identity", position = position_dodge(.2), linewidth = 1) + geom_errorbar(aes(ymin = lower.CL, ymax = upper.CL), width = .2, position = position_dodge(.2)) + labs(x = "Sampling Day", y = y_sox, colour = "Treatment") + scale_color_manual(values = c("blue", "orange", "light blue", "dark blue", "red"), limits = c("MPTC", "FPTC", "LHW", "DHW", "DHW+F")) + theme_classic() #####Foxl2 #looking at distribution plotdist(spd$NFoxl2, histo = TRUE, demp = TRUE) descdist(spd$NFoxl2, discrete = FALSE, boot = 500) #running glm foxmod <- glm(log(NFoxl2)~ TRT * Day, family = gaussian, data = spd) summary(foxmod) plot(fitted(foxmod), residuals(foxmod), xlab = "Fitted Values", ylab = "Residuals") abline(h = 0, lty = 2) lines(smooth.spline(fitted(foxmod), residuals(foxmod))) foxmod_simres <- simulateResiduals(foxmod) plot(foxmod_simres) drop1(foxmod, test = "Chisq") Anova(foxmod, type = "III") foxemm <- emmeans(foxmod, ~TRT:Day, type = "response") foxemm DHWfox.em = c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0) DHWFfox.em = c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0) FPTCfox.em = c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0) LHWfox.em = c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0) MPTCfox.em = c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1) contrast(foxemm, method = list("DHW - DHWF" = DHWfox.em - DHWFfox.em, "DHW - FPTC" = DHWfox.em - FPTCfox.em, "DHW - LHW" = DHWfox.em - LHWfox.em, "DHW - MPTC" = DHWfox.em - MPTCfox.em, "DHWF - FPTC" = DHWFfox.em - FPTCfox.em, "DHWF - LHW" = DHWFfox.em - LHWfox.em, "DHWF - MPTC" = DHWFfox.em - MPTCfox.em, "FPTC - LHW" = FPTCfox.em - LHWfox.em, "FPTC - MPTC" = FPTCfox.em - MPTCfox.em, "LHW - MPTC" = LHWfox.em - MPTCfox.em), adjust = "fdr") y_fox <- expression(paste("Normalized ", italic("Foxl2"), " expression ± 95% CI")) foxgraph<-summary(foxemm, infer = TRUE, type = "response") foxgraph foxgraph %>% ggplot(aes(x = Day, y = response, color = TRT, group = TRT)) + geom_line(stat = "identity", position = position_dodge(.2), linewidth = 1) + geom_errorbar(aes(ymin = lower.CL, ymax = upper.CL), width = .2, position = position_dodge(.2)) + labs(x = "Sampling Day", y = y_fox, colour = "Treatment") + scale_color_manual(values = c("blue", "orange", "light blue", "dark blue", "red"), limits = c("MPTC", "FPTC", "LHW", "DHW", "DHW+F")) + theme_classic() + theme(axis.title = element_text(size = 13), axis.text = element_text(size = 12)) #####Cyp19A1 #looking at distribution plotdist(spd$NCyp19A1, histo = TRUE, demp = TRUE) descdist(spd$NCyp19A1, discrete = FALSE, boot = 500) #running glm cypsqrt <- glm(sqrt(NCyp19A1)~ TRT * Day, family = "gaussian", data = spd) summary(cypsqrt) #looking at residuals plot(fitted(cypsqrt), residuals(cypsqrt), xlab = "Fitted Values", ylab = "Residuals") abline(h = 0, lty = 2) lines(smooth.spline(fitted(cypsqrt), residuals(cypsqrt))) #checking assumptions - passed with sqrt transformation cypsqrt_simres <- simulateResiduals(cypsqrt) plot(cypsqrt_simres) drop1(cypsqrt, test = "Chisq") #final results Anova(cypsqrt, type = "III") #creating emmeans cypemm <- emmeans(cypsqrt, ~TRT:Day, type = "response") cypemm #creating contrast strings DHW.cyp = c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0) DHWF.cyp = c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0) FPTC.cyp = c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0) LHW.cyp = c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0) MPTC.cyp = c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1) #constrasts with fdr adjustment contrast(cypemm, method = list("DHW - DHWF" = DHW.cyp - DHWF.cyp, "DHW - FPCT" = DHW.cyp - FPTC.cyp, "DHW - LHW" = DHW.cyp - LHW.cyp, "DHW - MPTC" = DHW.cyp - MPTC.cyp, "DHWF - FPTC" = DHWF.cyp - FPTC.cyp, "DHWF - LHW" = DHWF.cyp - LHW.cyp, "DHWF - MPTC" = DHWF.cyp - MPTC.cyp, "FPTC - LHW" = FPTC.cyp - LHW.cyp, "FPTC - MPTC" = FPTC.cyp - MPTC.cyp, "LHW - MPTC" = LHW.cyp - MPTC.cyp), adjust = "fdr") #graphing emmeans y_cyp <- expression(paste("Normalized ", italic("Cyp19A1"), " expression ± 95% CI")) cypgraph<-summary(cypemm, infer = TRUE, type = "response") cypgraph cypgraph %>% ggplot(aes(x = Day, y = response, color = TRT, group = TRT)) + geom_line(stat = "identity", position = position_dodge(.2), linewidth = 1) + geom_errorbar(aes(ymin = lower.CL, ymax = upper.CL), width = .2, position = position_dodge(.2)) + labs(x = "Sampling Day", y = y_cyp, colour = "Treatment") + scale_color_manual(values = c("blue", "orange", "light blue", "dark blue", "red"), limits = c("MPTC", "FPTC", "LHW", "DHW", "DHW+F")) + theme_classic() + theme(axis.title = element_text(size = 13), axis.text = element_text(size = 12)) #brining in one treatment at a time for SICB 2023 cypgraph %>% filter(TRT == "FPTC") %>% ggplot(aes(x = Day, y = response, color = TRT, group = TRT)) + geom_line(stat = "identity", position = position_dodge(.2), linewidth = 1) + geom_errorbar(aes(ymin = lower.CL, ymax = upper.CL), width = .2, position = position_dodge(.2)) + labs(x = "Sampling Day", y = y_cyp, colour = "Treatment") + scale_color_manual(values = c("blue", "orange", "light blue", "dark blue", "red"), limits = c("MPTC", "FPTC", "LHW", "DHW", "DHW+F")) + theme_classic() + theme(axis.title = element_text(size = 13), axis.text = element_text(size = 12)) + coord_cartesian(ylim = c(0, 0.008)) + theme(legend.position = "none") cypgraph %>% filter(TRT %in% c("FPTC", "MPTC")) %>% ggplot(aes(x = Day, y = response, color = TRT, group = TRT)) + geom_line(stat = "identity", position = position_dodge(.2), linewidth = 1) + geom_errorbar(aes(ymin = lower.CL, ymax = upper.CL), width = .2, position = position_dodge(.2)) + labs(x = "Sampling Day", y = y_cyp, colour = "Treatment") + scale_color_manual(values = c("blue", "orange", "light blue", "dark blue", "red"), limits = c("MPTC", "FPTC", "LHW", "DHW", "DHW+F")) + theme_classic() + theme(axis.title = element_text(size = 13), axis.text = element_text(size = 12)) + coord_cartesian(ylim = c(0, 0.008)) + theme(legend.position = "none") cypgraph %>% filter(TRT %in% c("FPTC", "MPTC", "LHW")) %>% ggplot(aes(x = Day, y = response, color = TRT, group = TRT)) + geom_line(stat = "identity", position = position_dodge(.2), linewidth = 1) + geom_errorbar(aes(ymin = lower.CL, ymax = upper.CL), width = .2, position = position_dodge(.2)) + labs(x = "Sampling Day", y = y_cyp, colour = "Treatment") + scale_color_manual(values = c("blue", "orange", "light blue", "dark blue", "red"), limits = c("MPTC", "FPTC", "LHW", "DHW", "DHW+F")) + theme_classic() + theme(axis.title = element_text(size = 13), axis.text = element_text(size = 12)) + coord_cartesian(ylim = c(0, 0.008)) + theme(legend.position = "none") cypgraph %>% filter(TRT %in% c("FPTC", "MPTC", "LHW", "DHW")) %>% ggplot(aes(x = Day, y = response, color = TRT, group = TRT)) + geom_line(stat = "identity", position = position_dodge(.2), linewidth = 1) + geom_errorbar(aes(ymin = lower.CL, ymax = upper.CL), width = .2, position = position_dodge(.2)) + labs(x = "Sampling Day", y = y_cyp, colour = "Treatment") + scale_color_manual(values = c("blue", "orange", "light blue", "dark blue", "red"), limits = c("MPTC", "FPTC", "LHW", "DHW", "DHW+F")) + theme_classic() + theme(axis.title = element_text(size = 13), axis.text = element_text(size = 12)) + coord_cartesian(ylim = c(0, 0.008)) + theme(legend.position = "none") cypgraph %>% filter(TRT %in% c("FPTC", "MPTC", "LHW", "DHW", "DHW+F")) %>% ggplot(aes(x = Day, y = response, color = TRT, group = TRT)) + geom_line(stat = "identity", position = position_dodge(.2), linewidth = 1) + geom_errorbar(aes(ymin = lower.CL, ymax = upper.CL), width = .2, position = position_dodge(.2)) + labs(x = "Sampling Day", y = y_cyp, colour = "Treatment") + scale_color_manual(values = c("blue", "orange", "light blue", "dark blue", "red"), limits = c("MPTC", "FPTC", "LHW", "DHW", "DHW+F")) + theme_classic() + theme(axis.title = element_text(size = 13), axis.text = element_text(size = 12)) + coord_cartesian(ylim = c(0, 0.008)) + theme(legend.position = "none")