################################################################### ## Data analysis placental thyroid hormones HPLC (Stapleton Lab) ## Set working directory library(dplyr) library(ggplot2) #Old working directory #setwd("C:/Users/Public/Documents/My Documents/CD1 PFOA GenX/Stapleton HPLC thyroid") setwd("C:/Users/blakebe/Documents/CD1 PFOA GenX/Stapleton HPLC thyroid") x<-read.csv("TH data R_imputed_vals.csv") x$T3.T4<-x$T3/x$T4 x$rT3.T4<-x$rT3/x$T4 x$group<-factor(x$group, levels=c("Vehicle Control","PFOA 1 mg/kg","PFOA 5 mg/kg", "GenX 2 mg/kg","GenX 10 mg/kg")) x$group <- relevel(x$group, ref = "Vehicle Control") library(Hmisc) Hmisc::describe(x) names(x) sex.summary.stats<-x%>%group_by(group, sex)%>%summarise_at(vars(wt.mg:rT3.T4), funs(mean, sd), na.rm=T) sex.summary.stats summary.stats<-x%>%group_by(group)%>%summarise_at(vars(wt.mg:rT3.T4), funs(mean, sd), na.rm=T) summary.stats #write.csv(sex.summary.stats,"Placental thyroid hormone summary stats by sex and group.csv",row.names = F) #write.csv(summary.stats,"Placental thyroid hormone summary stats by group.csv",row.names = F) library(lme4) library(multcomp) xm<-subset(x, sex=="M") xf<-subset(x, sex=="F") m1<-lm(T4~group*sex, data=x) m2<-lm(T4~group+sex, data=x) anova(m1,m2) ## No sex interaction m3<-lm(T4~group, data=x) anova(m2,m3) ## No overall effect of sex summary(glht(aov(T4~sex+group, data=x), linfct=mcp(group="Tukey", sex="Tukey"))) ## Significant effect of 10 mg/kg GenX even with sex covariate summary(glht(lm(T4~group, data=x), linfct=mcp(group="Tukey"))) ## Significant effect of 10 mg/kg GenX without sex covariate summary(glht(lm(T4~group, data=xm))) summary(glht(lm(T4~group, data=xf))) ## No significant effects with post-hoc correction for either sex analyzed separately m1<-lm(T3~group*sex, data=x) m2<-lm(T3~group+sex, data=x) anova(m1,m2) ## No sex interaction m3<-lm(T3~group, data=x) anova(m2,m3) ## No overall effect of sex summary(glht(aov(T3~sex+group, data=x), linfct=mcp(group="Tukey", sex="Tukey"))) ## No effects post-hoc summary(glht(lm(T3~group, data=x), linfct=mcp(group="Tukey"))) ## No effects post-hoc summary(glht(lm(T3~group, data=xm))) ## Errors because too many males with imputed vals (can't reliably estimate effects) summary(glht(lm(T3~group, data=xf))) ## No significant effects with post-hoc correction m1<-lm(rT3~group*sex, data=x) m2<-lm(rT3~group+sex, data=x) anova(m1,m2) ## No sex interaction m3<-lm(rT3~group, data=x) anova(m2,m3) ## No overall effect of sex (marginal) summary(glht(aov(rT3~sex+group, data=x), linfct=mcp(group="Tukey", sex="Tukey"))) ## No effects post-hoc summary(glht(lm(rT3~group, data=x), linfct=mcp(group="Tukey"))) ## No effects post-hoc summary(glht(lm(rT3~group, data=xm))) ## No significant effects with post-hoc (males in 1 mg/kg PFOA marginally lower than controls) summary(glht(lm(rT3~group, data=xf))) ## No significant effects with post-hoc correction (females in 2 mg/kg GenX marginally higher than controls) m1<-lm(T3.T4~group*sex, data=x) m2<-lm(T3.T4~group+sex, data=x) anova(m1,m2) ## No sex interaction m3<-lm(T3.T4~group, data=x) anova(m2,m3) ## No overall effect of sex (marginal) summary(glht(aov(T3.T4~sex+group, data=x), linfct=mcp(group="Tukey", sex="Tukey"))) ## No significant effects post-hoc (10 mg/kg GenX marginal) summary(glht(lm(T3.T4~group, data=x), linfct=mcp(group="Tukey"))) ## No significant effects post-hoc (10 mg/kg GenX marginal) summary(glht(lm(T3.T4~group, data=xm))) ## No significant effects summary(glht(lm(T3.T4~group, data=xf))) ## No significant effects (10 mg/kg GenX females marginal) m1<-lm(rT3.T4~group*sex, data=x) m2<-lm(rT3.T4~group+sex, data=x) anova(m1,m2) ## No sex interaction m3<-lm(rT3.T4~group, data=x) anova(m2,m3) ## No overall effect of sex summary(glht(aov(rT3.T4~sex+group, data=x), linfct=mcp(group="Tukey", sex="Tukey"))) ## No significant effects post-hoc summary(glht(lm(rT3.T4~group, data=x), linfct=mcp(group="Tukey"))) ## No significant effects post-hoc summary(glht(lm(rT3.T4~group, data=xm))) ## No significant effects summary(glht(lm(rT3.T4~group, data=xf))) ## No significant effects x$wt.mg_mean<-x$wt.mg p1<-ggplot(sex.summary.stats, aes(x=group, y=wt.mg_mean, group=group))+ theme_bw()+ theme(legend.position="none", axis.text=element_text(size=12), axis.text.x=element_text(angle =-45, hjust=0), axis.title=element_text(size=14))+ geom_jitter(data=x, aes(color=group),size=2,stat="identity",width=0.1)+ geom_point(data=sex.summary.stats, aes(x=group, y=wt.mg_mean),pch=3, size=4)+ geom_errorbar(data=sex.summary.stats,aes(ymin=wt.mg_mean-wt.mg_sd, ymax=wt.mg_mean+wt.mg_sd),size=0.6,width=0.3)+ #scale_color_manual(values=c("dodgerblue","springgreen3","darkorange","firebrick1","darkgoldenrod1"))+ facet_grid(.~sex)+ ggtitle("Whole placenta weight")+ ylab("Weight (mg)\n")+ xlab("Group") p1 x$T4_mean<-x$T4 p2<-ggplot(summary.stats, aes(x=group, y=T4_mean, group=group))+ theme_bw()+ theme(#legend.position="none", axis.text=element_text(size=12), axis.text.x=element_text(angle =-45, hjust=0), axis.title=element_text(size=14))+ geom_jitter(data=x, aes(x=group, y=T4_mean, color=group, shape=sex),size=2,stat="identity",width=0.2)+ geom_point(data=summary.stats, aes(x=group, y=T4_mean),pch=3, size=4)+ geom_errorbar(data=summary.stats, aes(x=group, ymin=T4_mean-T4_sd, ymax=T4_mean+T4_sd),size=0.6,width=0.3)+ scale_color_manual(values=c("gray","darkgoldenrod1","darkgoldenrod","deepskyblue2","deepskyblue4"))+ #scale_color_manual(values=c("dodgerblue","springgreen3","darkorange","firebrick1","darkgoldenrod1"))+ #facet_grid(.~sex)+ ggtitle("T4")+ ylab("ng/g placenta \n")+ xlab("Group") p2 x$T3_mean<-x$T3 p1<-ggplot(sex.summary.stats, aes(x=group, y=T3_mean, group=group))+ theme_bw()+ theme(legend.position="none", axis.text=element_text(size=12), axis.text.x=element_text(angle =-45, hjust=0), axis.title=element_text(size=14))+ geom_jitter(data=x, aes(color=group),size=2,stat="identity",width=0.1)+ geom_point(data=sex.summary.stats, aes(x=group, y=T3_mean),pch=3, size=4)+ geom_errorbar(data=sex.summary.stats,aes(ymin=T3_mean-T3_sd, ymax=T3_mean+T3_sd),size=0.6,width=0.3)+ #scale_color_manual(values=c("dodgerblue","springgreen3","darkorange","firebrick1","darkgoldenrod1"))+ facet_grid(.~sex)+ ggtitle("T3")+ ylab("ng/g placenta \n")+ xlab("Group") p1 x$rT3_mean<-x$rT3 p1<-ggplot(sex.summary.stats, aes(x=group, y=rT3_mean, group=group))+ theme_bw()+ theme(legend.position="none", axis.text=element_text(size=12), axis.text.x=element_text(angle =-45, hjust=0), axis.title=element_text(size=14))+ geom_jitter(data=x, aes(color=group),size=2,stat="identity",width=0.1)+ geom_point(data=sex.summary.stats, aes(x=group, y=rT3_mean),pch=3, size=4)+ geom_errorbar(data=sex.summary.stats,aes(ymin=rT3_mean-rT3_sd, ymax=rT3_mean+rT3_sd),size=0.6,width=0.3)+ #scale_color_manual(values=c("dodgerblue","springgreen3","darkorange","firebrick1","darkgoldenrod1"))+ facet_grid(.~sex)+ ggtitle("rT3")+ ylab("ng/g placenta \n")+ xlab("Group") p1 ############################################## ## Testing using gather function to facet plots (bc i'm lazy) library(dplyr) library(ggplot2) library(tidyr) library(ggsci) setwd("C:/Users/blakebe/Documents/CD1 PFOA GenX/Stapleton HPLC thyroid") x<-read.csv("TH data R_imputed_vals.csv") x$T3.T4<-x$T3/x$T4 x$rT3.T4<-x$rT3/x$T4 x$group<-factor(x$group, levels=c("Vehicle Control","PFOA 1 mg/kg","PFOA 5 mg/kg", "GenX 2 mg/kg","GenX 10 mg/kg")) x$group <- relevel(x$group, ref = "Vehicle Control") summary.stats<-x%>%group_by(group)%>%summarise_at(vars(rT3:rT3.T4), funs(mean, sd), na.rm=T) summary.stats summary.stats<-as.data.frame(summary.stats) x.gather<-gather(x,"hormone","ng.g", 5:9) summary.gather<-gather(summary.stats, "hormone","ng.g_mean", 2:6) head(summary.gather) summary.gather.test<-summary.gather %>% separate(hormone, c("hormone", NA),"_") head(summary.gather.test) summary.gather.mean<-summary.gather.test[, -c(2:6)] head(summary.gather.mean) summary.gather2<-gather(summary.stats, "hormone","ng.g_sd", 7:11) summary.gather2.test<-summary.gather2 %>% separate(hormone, c("hormone", NA),"_") head(summary.gather2.test) summary.gather.sd<-summary.gather2.test[, -c(2:6)] head(summary.gather.sd) summary.gathered<-merge(summary.gather.mean, summary.gather.sd, by=c("group","hormone")) head(summary.gathered) str(summary.gathered) summary.gathered$hormone<-as.factor(summary.gathered$hormone) str(summary.gathered) x.gather$ng.g_mean<-x.gather$ng.g str(x.gather) x.gather$hormone<-as.factor(x.gather$hormone) summary.gathered$hormone<-factor(summary.gathered$hormone, levels=c("T4","T3","rT3","T3.T4","rT3.T4")) t4.only<-subset(summary.gathered, hormone=="T4") x.t4<-subset(x.gather, hormone=="T4") setwd("C:/Users/blakebe/Documents/CD1 PFOA GenX/R plots") tiff("Placental T4.tiff", units="in", width=4, height=4, res=300) p2<-ggplot(t4.only, aes(x=group, y=ng.g_mean, group=group))+ theme_bw()+ theme(legend.position="none", axis.text=element_text(size=12), axis.text.x=element_text(angle =-45, hjust=0), plot.margin = margin(20, 40, 20, 20), axis.title=element_text(size=14))+ geom_jitter(data=x.t4, aes(x=group, y=ng.g_mean, color=group, shape=sex),size=2,stat="identity",width=0.2)+ geom_point(data=t4.only, aes(x=group, y=ng.g_mean),pch=3, size=4)+ geom_errorbar(data=t4.only, aes(x=group, ymin=ng.g_mean-ng.g_sd, ymax=ng.g_mean+ng.g_sd),size=0.6,width=0.3)+ scale_color_manual(values=c("gray","darkgoldenrod1","darkgoldenrod","deepskyblue2","deepskyblue4"))+ #scale_color_manual(values=c("dodgerblue","springgreen3","darkorange","firebrick1","darkgoldenrod1"))+ #facet_wrap(.~hormone, scales="free_y",nrow = 2)+ #ggtitle("Thyroid hormone analysis")+ ylab("ng/g placenta \n")+ ylim(0,10)+ xlab("Group") p2 dev.off()