library(tidyr) library(dplyr) library(ggplot2) getwd() ## Set working directory ## Old working directory #setwd("C:/Users/Public/Documents/My Documents/CD1 PFOA GenX/R CSV files") setwd("C:/Users/blakebe/Documents/CD1 PFOA GenX/R CSV files") ## Read in data file containing animal IDs, treatment groups, etc x<-read.csv("Combined necropsy data E11.5 and E17.5_updated.csv") ## Set working directory ## Old working directory #setwd("C:/Users/Public/Documents/My Documents/CD1 PFOA GenX/HPLC sample lists") setwd("C:/Users/blakebe/Documents/CD1 PFOA GenX/HPLC sample lists") ## Read in data file containing amniotic fluid and maternal serum HPLC values y<-read.csv("UPLC QQQ CD-1 PFOA GenX serum AF gestational exposure study data BEB 05242019.csv") ## Check structure of data str(x) str(y) names(y) test<-spread(y, sample, ng.ml) names(test) af<-na.omit(test[c(1:5)]) ser<-na.omit(test[c(1:4,6)]) test.2<-left_join(ser, af, by=c("dam.id","timepoint","group")) names(test.2)[names(test.2)=="Amniotic fluid"]<-"AF" #write.csv(test.2, "Serum and AF dosimetry.csv", row.names = F) names(test.2) fluid.ug<-test.2 names(fluid.ug) fluid.ug$af.ug<-fluid.ug$AF/1000 fluid.ug$ser.ug<-fluid.ug$Serum/1000 fluid.ug$group<-factor(fluid.ug$group, levels=c("PFOA 1 mg/kg","PFOA 5 mg/kg","GenX 2 mg/kg","GenX 10 mg/kg")) levels(fluid.ug$group) summary.stats<-fluid.ug%>% dplyr::group_by(timepoint, group) %>% dplyr::summarise(#nobs.ser = nobs(Serum, na.rm = TRUE), #nobs.af = nobs(AF, na.rm = TRUE), mean.af = mean(af.ug, na.rm = TRUE), sd.af = sd(af.ug, na.rm = TRUE), mean.ser = mean(ser.ug, na.rm = TRUE), sd.ser = sd(ser.ug, na.rm = TRUE), mean.ser.af = mean(ser.ug/af.ug, na.rm = TRUE), sd.ser.af = sd(ser.ug/af.ug, na.rm = TRUE), mean.af.ser = mean(af.ug/ser.ug, na.rm = TRUE), sd.af.ser = sd(af.ug/ser.ug, na.rm = TRUE)) summary.stats ## Plotting serum fluid.ug$mean.ser<-fluid.ug$ser.ug p1<-ggplot(summary.stats, aes(x=group, y=mean.ser, 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), plot.margin=unit(c(2,2,2,2),"cm"))+ geom_jitter(data=fluid.ug, aes(color=group),size=2,stat="identity",width=0.1)+ geom_point(data=summary.stats, aes(x=group, y=mean.ser),pch=3, size=4)+ geom_errorbar(data=summary.stats,aes(ymin=mean.ser-sd.ser, ymax=mean.ser+sd.ser),size=0.6,width=0.3)+ scale_color_manual(values=c("darkgoldenrod1","darkgoldenrod","deepskyblue2","deepskyblue4"))+ facet_grid(.~timepoint)+ ggtitle("Maternal serum")+ ylab("PFAS concentration (ug/mL)\n")+ xlab("\nGroup") p1 tiff("Maternal serum PFAS.tiff", units="in", width=7, height=6, res=300) # insert ggplot code p1 dev.off() fluid.ug$mean.af<-fluid.ug$af.ug p2<-ggplot(summary.stats, aes(x=group, y=mean.af, 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), plot.margin=unit(c(2,2,2,2),"cm"))+ geom_jitter(data=fluid.ug, aes(color=group),size=2,stat="identity",width=0.1)+ geom_point(data=summary.stats, aes(x=group, y=mean.af),pch=3, size=4)+ geom_errorbar(data=summary.stats,aes(ymin=mean.af-sd.af, ymax=mean.af+sd.af),size=0.6,width=0.3)+ scale_color_manual(values=c("darkgoldenrod1","darkgoldenrod","deepskyblue2","deepskyblue4"))+ #facet_grid(.~timepoint)+ ggtitle("Amniotic fluid")+ ylim(0,15)+ ylab("PFAS concentration (ug/mL)\n")+ xlab("\nGroup") p2 tiff("Amniotic fluid PFAS.tiff", units="in", width=6, height=6, res=300) # insert ggplot code p2 dev.off() y1<-merge(y, x, by=c("dam.id","timepoint")) names(y1)[names(y1)=="group.y"]<-"color.group" names(y1)[names(y1)=="group.x"]<-"group" names(y1) ser<-subset(y1, sample="Serum") af<-subset(y1, sample="Amniotic fluid") names(ser)[names(ser) == "ng.ml"] <- "ng.ml.ser" names(af)[names(af) == "ng.ml"] <- "ng.ml.af" af2 <- af[c(1,6)] ser2 <- ser[c(1,6)] ## Read in data file containing maternal liver and whole embryo HPLC data z<-read.csv("UPLC QQQ CD-1 PFOA GenX Liver Embryo gestational exposure study data BEB 05282019.csv") names(z) z$group<-factor(z$group, levels=c("PFOA 1 mg/kg","PFOA 5 mg/kg","GenX 2 mg/kg","GenX 10 mg/kg")) levels(z$group) z$ug.g<-z$ng.g/1000 summary.stats2<-z%>% group_by(group,timepoint,type) %>% summarise(mean.ug.g=mean(ug.g), sd.ug.g=sd(ug.g)) summary.stats2 liver.p<-subset(z, type=="Dam liver") liver.m<-subset(summary.stats2, type=="Dam liver") liver.p$mean.ug.g<-liver.p$ug.g p3<-ggplot(liver.m, aes(x=group, y=mean.ug.g, 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), plot.margin=unit(c(2,2,2,2),"cm"))+ geom_jitter(data=liver.p, aes(color=group),size=2,stat="identity",width=0.1)+ geom_point(data=liver.m, aes(x=group, y=mean.ug.g),pch=3, size=4)+ geom_errorbar(data=liver.m,aes(ymin=mean.ug.g-sd.ug.g, ymax=mean.ug.g+sd.ug.g),size=0.6,width=0.3)+ scale_color_manual(values=c("darkgoldenrod1","darkgoldenrod","deepskyblue2","deepskyblue4"))+ facet_grid(.~timepoint)+ ggtitle("Maternal liver")+ ylab("PFAS concentration (ug/g)\n")+ xlab("\nGroup") p3 tiff("Maternal liver PFAS.tiff", units="in", width=7, height=6, res=300) # insert ggplot code p3 dev.off() emb.p<-subset(z, !type=="Dam liver") emb.m<-subset(summary.stats2, !type=="Dam liver") emb.p$mean.ug.g<-emb.p$ug.g emb.p11<-subset(emb.p, timepoint=="E11.5") emb.m11<-subset(emb.m, timepoint=="E11.5") p4<-ggplot(emb.m11, aes(x=group, y=mean.ug.g, 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), plot.margin=unit(c(2,2,2,2),"cm"))+ geom_jitter(data=emb.p11, aes(color=group),size=2,stat="identity",width=0.1)+ geom_point(data=emb.m11, aes(x=group, y=mean.ug.g),pch=3, size=4)+ geom_errorbar(data=emb.m11, aes(ymin=mean.ug.g-sd.ug.g, ymax=mean.ug.g+sd.ug.g),size=0.6,width=0.3)+ scale_color_manual(values=c("darkgoldenrod1","darkgoldenrod","deepskyblue2","deepskyblue4"))+ #facet_grid(.~timepoint)+ ggtitle("E11.5 Whole embryo")+ ylim(0,5)+ ylab("PFAS concentration (ug/g)\n")+ xlab("\nGroup") p4 tiff("E11.5 whole embryo PFAS.tiff", units="in", width=6, height=6, res=300) # insert ggplot code p4 dev.off() emb.p17<-subset(emb.p, timepoint=="E17.5") emb.m17<-subset(emb.m, timepoint=="E17.5") p5<-ggplot(emb.m17, aes(x=group, y=mean.ug.g, 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), plot.margin=unit(c(2,2,2,2),"cm"))+ geom_jitter(data=emb.p17, aes(color=group),size=2,stat="identity",width=0.1)+ geom_point(data=emb.m17, aes(x=group, y=mean.ug.g),pch=3, size=4)+ geom_errorbar(data=emb.m17, aes(ymin=mean.ug.g-sd.ug.g, ymax=mean.ug.g+sd.ug.g),size=0.6,width=0.3)+ scale_color_manual(values=c("darkgoldenrod1","darkgoldenrod","deepskyblue2","deepskyblue4"))+ facet_grid(.~type)+ ggtitle("Whole embryo")+ ylim(0,25)+ ylab("PFAS concentration (ug/g)\n")+ xlab("\nGroup") p5 tiff("E17.5 Whole embryo PFAS.tiff", units="in", width=7, height=6, res=300) # insert ggplot code p5 dev.off() e11.5<-subset(z, timepoint=="E11.5") e17.5<-subset(z, timepoint=="E17.5") e11.5.liv<-subset(e11.5, type=="Dam liver") e11.5.emb<-subset(e11.5, type=="Embryo") names(e11.5.liv)[names(e11.5.liv) == "ng.g"] <- "ng.g.liv" names(e11.5.emb)[names(e11.5.emb) == "ng.g"] <- "ng.g.emb" liver2 <- e11.5.liv[c(1,3,15)] embryo2 <- e11.5.emb[c(1,3,15)] e11.5_spread<-merge(liver2, embryo2, by=c("group","dam.id")) e11.5_spread$emb.liv.rat<-e11.5_spread$ng.g.emb/e11.5_spread$ng.g.liv e11.5_spread$liv.emb.rat<-e11.5_spread$ng.g.liv/e11.5_spread$ng.g.emb #write.csv(e11.5_spread, "Embryo Liver internal dosimetry E11.5.csv", row.names=F) summary.stats3<-e11.5_spread%>% group_by(group) %>% summarise(mean.liv.emb.ratio=mean(liv.emb.rat), sd.liv.emb.ratio=sd(liv.emb.rat), mean.emb.liv.ratio=mean(emb.liv.rat), sd.emb.liv.ratio=sd(emb.liv.rat)) summary.stats3 e17.5.liv<-subset(e17.5, type=="Dam liver") e17.5.embf<-subset(e17.5, type=="F Fetus") e17.5.embm<-subset(e17.5, type=="M Fetus") names(e17.5.liv)[names(e17.5.liv) == "ng.g"] <- "ng.g.liv" names(e17.5.embf)[names(e17.5.embf) == "ng.g"] <- "ng.g.embf" names(e17.5.embm)[names(e17.5.embm) == "ng.g"] <- "ng.g.embm" liver2 <- e17.5.liv[c(1,3,15)] embryoF2 <- e17.5.embf[c(1,3,15)] embryoM2 <- e17.5.embm[c(1,3,15)] e17.5_embryo<-merge(embryoF2, embryoM2, by=c("group","dam.id"), all.x = T) e17.5_spread<-merge(liver2, e17.5_embryo, by=c("group","dam.id")) #write.csv(e17.5_spread, "Embryo Liver internal dosimetry E17.5.csv", row.names=F) e17.5_spread$liv.embf.rat<-e17.5_spread$ng.g.embf/e17.5_spread$ng.g.liv e17.5_spread$embf.liv.rat<-e17.5_spread$ng.g.liv/e17.5_spread$ng.g.embf e17.5_spread$liv.embm.rat<-e17.5_spread$ng.g.embm/e17.5_spread$ng.g.liv e17.5_spread$embm.liv.rat<-e17.5_spread$ng.g.liv/e17.5_spread$ng.g.embm summary.stats5<-e17.5_spread%>% group_by(group) %>% summarise(mean.liv.embf.ratio=mean(liv.embf.rat, na.rm = TRUE), sd.liv.embf.ratio=sd(liv.embf.rat, na.rm = TRUE), mean.embf.liv.ratio=mean(embf.liv.rat, na.rm = TRUE), sd.embf.liv.ratio=sd(embf.liv.rat, na.rm = TRUE), mean.liv.embm.ratio=mean(liv.embm.rat, na.rm = TRUE), sd.liv.embm.ratio=sd(liv.embm.rat, na.rm = TRUE), mean.embm.liv.ratio=mean(embm.liv.rat, na.rm = TRUE), sd.embm.liv.ratio=sd(embm.liv.rat, na.rm = TRUE)) summary.stats5 e17.5<-subset(z, timepoint=="E17.5") e17.5.test<-e17.5%>% filter(!type=="Dam liver")%>% group_by(group, dam.id) %>% summarise(ng.g.emb.mean=mean(ng.g)) e17.5.test e17.5.liv<-subset(e17.5, type=="Dam liver") names(e17.5.liv)[names(e17.5.liv) == "ng.g"] <- "ng.g.liv" e17.5_liv.emb<-merge(e17.5.liv, e17.5.test, by=c("group","dam.id")) e17.5_liv.emb$liv.emb.rat<-e17.5_liv.emb$ng.g.liv/e17.5_liv.emb$ng.g.emb.mean e17.5_liv.emb$emb.liv.rat<-e17.5_liv.emb$ng.g.emb.mean/e17.5_liv.emb$ng.g.liv summary.stats6<-e17.5_liv.emb%>% group_by(group) %>% summarise(mean.liv.emb.ratio=mean(liv.emb.rat, na.rm = TRUE), sd.liv.emb.ratio=sd(liv.emb.rat, na.rm = TRUE), mean.emb.liv.ratio=mean(emb.liv.rat, na.rm = TRUE), sd.emb.liv.ratio=sd(emb.liv.rat, na.rm = TRUE), mean.emb.ng.g=mean(ng.g.emb.mean, na.rm = TRUE), sd.emb.ng.g=sd(ng.g.emb.mean, na.rm = TRUE)) summary.stats6