## CD-1 PFOA/GenX gestational exposure and latent health outcomes study
## Week 6 necropsy data
## Body weight, liver weight, relative liver weight
## HAC 01/22/2020 (w/ QC data)

## Load packages (needs to be done every time R is opened)
library(dplyr)
library(lme4)
library(ggplot2)
library(multcomp)
library(tidyverse)

# Set your working directory (replace this example with your file folder destination)
setwd("/Users/harliecope/Desktop/GenX:PFOA Study/GenX PFOA Mouse Study/GenXPFOA Study")

# Read in your files
x<- read_csv("~/Desktop/GenX:PFOA Study/GenX PFOA Mouse Study/QC Data Files/csv for R/QCWeek6NecropsyClean.csv")
x$ID <-as.factor(x$ID)
x$Dam.ID<-as.factor(x$Dam.ID)
x$Group<-as.factor(x$Group)
x$Diet<-as.factor(x$Diet)

str(x)

w6necrop<-x


names(w6necrop)

w6necrop.summary.stats<-w6necrop%>%
  group_by(Group,Sex,Diet)%>%
  summarise(
    bw.mean=mean(BW, na.rm=T),
    bw.sd=sd(BW, na.rm=T),
    liv.mean=mean(LW, na.rm=T),
    liv.sd=sd(LW, na.rm=T),
    rel.liv.mean=mean(RelLivWt, na.rm=T),
    rel.liv.sd=sd(RelLivWt, na.rm=T))
w6necrop.summary.stats


write.csv(w6necrop.summary.stats, "Week6Necropsy Summary Statistics.csv",row.names = F)


### Statistical Analyses

## Statistical analysis
## Males

w6m<-subset(w6necrop, Sex=="M")
w6f<-subset(w6necrop, Sex=="F")

w6mhf<-subset(w6m, Diet=="H")
w6mlf<-subset(w6m, Diet=="L")
w6fhf<-subset(w6f, Diet=="H")
w6flf<-subset(w6f, Diet=="L")

w6fveh<-subset(w6f, Group=="Control")
w6mveh<-subset(w6m, Group=="Control")

## Male, HFD
summary(glht(aov(BW~Group,data=w6mhf), linfct=mcp(Group="Dunnett")))
## No differences

summary(glht(aov(LW~Group,data=w6mhf), linfct=mcp(Group="Dunnett")))
## No differences

summary(glht(aov(RelLivWt~Group,data=w6mhf), linfct=mcp(Group="Dunnett")))
## No differences


## Male, LFD
summary(glht(aov(BW~Group,data=w6mlf), linfct=mcp(Group="Dunnett")))
## Difference in GenX 2.0 mg/kg compared to control (p=0.0359)

summary(glht(aov(LW~Group,data=w6mlf), linfct=mcp(Group="Dunnett")))
## Difference between control and 1.0 mg/kg GenX (p=0.00447) and 2.0 mg/kg GenX (0.00998)

summary(glht(aov(RelLivWt~Group,data=w6mlf), linfct=mcp(Group="Dunnett")))
## Difference between 1.0 mg/kg GenX vs Control (p = 0.00237), 1.0 mg/kg PFOA and control (p=0.03580)


## Female, HFD
summary(glht(aov(BW~Group,data=w6fhf), linfct=mcp(Group="Dunnett")))
## No differences (Vehicle vs 0.2 mg/kg GenX p=0.0956, vehicle vs 1.0 mg/kg GenX p=0.0976)

summary(glht(aov(LW~Group,data=w6fhf), linfct=mcp(Group="Dunnett")))
## No differences

summary(glht(aov(RelLivWt~Group,data=w6fhf), linfct=mcp(Group="Dunnett")))
## No differences


## Female, LFD
summary(glht(aov(BW~Group,data=w6flf), linfct=mcp(Group="Dunnett")))
## No differences

summary(glht(aov(LW~Group,data=w6flf), linfct=mcp(Group="Dunnett")))
## No differences

summary(glht(aov(RelLivWt~Group,data=w6flf), linfct=mcp(Group="Dunnett")))
## No differences

## compare LFD and HFD
summary(glht(aov(BW~Diet, data=w6fveh)))
summary(glht(aov(BW~Diet, data=w6mveh)))
summary(glht(aov(LW~Diet, data=w6fveh)))
summary(glht(aov(LW~Diet, data=w6mveh)))
summary(glht(aov(RelLivWt~Diet, data=w6fveh)))
summary(glht(aov(RelLivWt~Diet, data=w6mveh)))