#Title: Seasonal variation in DNA damage is not related to time active
#Authors: Emma S. Rosen, Lily Mikolajczak, Ursula K. Beattie, L. Michael Romero
#Corresponding author: Ursula K. Beattie, ursula.beattie@tufts.edu


#change to your own directory
setwd("Documents/Grad School/2022 Summer/Emma Thesis/paper/PeerJ submission/")
#nothing else to change

#Load libraries
library(ggplot2)
library(dplyr) #this is for summary stats
library(plotrix) #this is for doing standard error
library(car) #levene test
library(multcomp) #pairwise comparisons
library(lme4) #lmer
library(tidyr)
library(patchwork)
library(ggprism) #for minor tick marks
library(sasLM)



#### Diel Rhythm Experiment ####
dielData=read.table("diel rhythm data.csv", sep=",", header=T)
dielData$Time = factor(dielData$Time, levels=c("midnight", "4am", "8am", "noon", "4pm", "8pm"))


dielAverages = dielData %>% 
  group_by(Time24) %>% 
  summarize_at(vars(TailMoment:UricAcid), list(mean=mean, se=std.error), na.rm=T)

#comets
comets = ggplot(data=dielAverages, aes(x=Time24, y=TailMoment_mean)) +
  geom_rect(xmin=-Inf, xmax=8, ymin=-Inf, ymax=Inf, fill="gray90") +
  geom_rect(xmin=20, xmax=Inf, ymin=-Inf, ymax=Inf, fill="gray90") +
  geom_point() + geom_line() + 
  geom_errorbar(aes(ymin = TailMoment_mean-TailMoment_se, 
                    ymax=TailMoment_mean+TailMoment_se), width=0.25) +
  xlab("Time of Day") + ylab("DNA Damage\n(Tail Moment)") +
  scale_x_continuous(breaks = c(0, 4, 8, 12, 16, 20), 
                     labels = c("00:00", "04:00", "08:00", "12:00", "16:00", "20:00"),
                     guide = guide_prism_minor()) +
  scale_y_continuous(limits = c(20,46), guide = guide_prism_minor()) +
  theme_classic()



leveneTest(TailMoment ~ factor(Time24), data=dielData) #over .05 so were good
model = lmer(TailMoment ~ Time24 + (1|Bird), data=dielData)
Anova(model, type = "III", test.statistic = "F") #time of day is not sig



#uric acid
uric = ggplot(data=dielAverages, aes(x=Time24, y=UricAcid_mean)) +
  geom_rect(xmin=-Inf, xmax=8, ymin=-Inf, ymax=Inf, fill="gray90") +
  geom_rect(xmin=20, xmax=Inf, ymin=-Inf, ymax=Inf, fill="gray90") +
  geom_point() + geom_line() + 
  geom_errorbar(aes(ymin = UricAcid_mean-UricAcid_se, 
                    ymax=UricAcid_mean+UricAcid_se), width=0.25) +
  geom_text(data=dielAverages, aes(x=0, y=12, label="a")) +
  geom_text(data=dielAverages, aes(x=4, y=8, label="b")) +
  geom_text(data=dielAverages, aes(x=8, y=12, label="abc")) +
  geom_text(data=dielAverages, aes(x=12, y=17, label="d")) +
  geom_text(data=dielAverages, aes(x=16, y=15, label="cd")) +
  geom_text(data=dielAverages, aes(x=20, y=12, label="a")) +
  xlab("Time of Day") + ylab("Uric Acid (mg/dL)") +
  scale_x_continuous(breaks = c(0, 4, 8, 12, 16, 20), 
                     labels = c("00:00", "04:00", "08:00", "12:00", "16:00", "20:00"),
                     guide = guide_prism_minor()) +
  scale_y_continuous(limits = c(0,20), guide = guide_prism_minor()) +
  theme_classic()


leveneTest(UricAcid ~ factor(Time24), data=dielData) #over .05 so were good
model = lmer(UricAcid ~ Time + (1|Bird), data=dielData)
Anova(model, type = "III", test.statistic = "F") #time of day is significant
summary(glht(model, mcp(Time="Tukey"))) #tukey posthoc test
cld(summary(glht(model, mcp(Time="Tukey")))) #to obtain significance letters

layout="
A
B"
#png("figure 2.png", width=5, height=4, units="in", res=300)
comets + theme(axis.title.x=element_blank(),
               axis.text.x=element_blank()) +
  uric + plot_layout(design=layout)  +
  plot_annotation(tag_levels = "A")
#dev.off()




#### Resonance Experiment ####
##### Tissue Data #####
cometData = read.table("tissue comet data.csv", sep=",", header=T)
cometData$Tissue = factor(cometData$Tissue, levels=c("B", "AF", "Hi", "Hy", "L")) #reorder tissues for graph
cometData$Total_Light = factor(cometData$Total_Light)

#take averages
cometAverages = cometData %>%
  group_by(Total_Light,Type,Tissue) %>%
  summarize_at(vars(AvgTailMoment:Initial.Comets), list(mean=mean, se=std.error), na.rm=T)
cometAverages$LightCycle = c(rep("9:15",5), rep("4.5:7.5",5), rep("12:12", 5), rep("6:6", 5))
cometAverages$LightCycle = factor(cometAverages$LightCycle, levels=c("12:12", "6:6", "9:15", "4.5:7.5"))

#graph
sig_markings = read.table("tissue_sig markings.csv", sep=",", header=T)
sig_markings$Tissue = factor(sig_markings$Tissue, levels=c("B", "AF", "Hi", "Hy", "L")) #reorder tissues for graph

labels = c(AF="Abdominal Fat", B="Blood", Hi="Hippocampus", Hy="Hypothalamus", L="Liver")
#png("resubmission/figure 3.png", width=6, height=4, units="in", res=300)
ggplot(data=cometAverages, aes(x=LightCycle, y=AvgTailMoment_mean)) + 
  geom_col(data=cometAverages, aes(x=LightCycle, y=Initial.Comets_mean, linetype="dashed"),
           color="grey50", fill="white") +
  geom_errorbar(aes(ymin=cometAverages$Initial.Comets_mean-
                      cometAverages$Initial.Comets_se, 
                    ymax=cometAverages$Initial.Comets_mean+
                      cometAverages$Initial.Comets_se), 
                width=0.2, linetype="solid", color="grey50")+
  geom_col(color="black", aes(fill=Type)) + 
  geom_errorbar(aes(ymin=cometAverages$AvgTailMoment_mean-
                      cometAverages$AvgTailMoment_se, 
                    ymax=cometAverages$AvgTailMoment_mean+
                      cometAverages$AvgTailMoment_se), 
                width=0.2, linetype="solid") +
  geom_text(data=sig_markings, aes(x=Text_x, y=Text_y), label="*", size=6) +
  geom_segment(data=sig_markings, aes(x=Line_x, xend=Line_xend, y=Line_y, yend=Line_yend)) +
  ylab("DNA Damage (Tail Moment)") + 
  scale_x_discrete(name="Light Cycle") +
  scale_y_continuous(limits=c(0,60), expand=c(0,0), guide = guide_prism_minor()) +
  scale_fill_manual(name = "Light Cycle Type", values=c("white", "grey80")) +
  scale_linetype_manual(name = "Sample Time", values=c("dashed"),
                        labels=c("Pre-resonance\n(12:12)")) +
  theme_classic() + theme(text = element_text(size = 10))+
  facet_wrap(~Tissue, scales="free", labeller=labeller(Tissue=labels)) +
  theme(axis.text.x = element_text(angle = 30, hjust=.9))
#dev.off()

#stats
#full model - slice on total light
SLICE(AvgTailMoment ~ Tissue * Total_Light * Type, cometData, "Type", "Total_Light")
#the model is significant in both 9hrs light and 12hrs light -> split to further break down differences

    #separate into groups of total light - now slice on tissue to get effect of light cycle type in each tissue
    #9 hours light 
    SLICE(AvgTailMoment ~ Tissue * Type, cometData[cometData$Total_Light==9,], "Type", "Tissue")
    #resonance group has higher DNA damage than natural in hypothalamus and liver
    
    #12 hrs light
    SLICE(AvgTailMoment ~ Tissue * Type, cometData[cometData$Total_Light==12,], "Type", "Tissue")
    #resonance group has lower DNA damage than natural in blood, abdominal fat, and hippocampus
       
        
##### Resonance - testes weights #####
testesData = read.table("testes data.csv", sep=",", header=T)
testesData$LightCycle = factor(testesData$LightCycle, levels=c("6:6", "4.5:7.5"))
testesAvg = testesData %>%
  group_by(Total_light,LightCycle) %>%
  summarize_at(vars(Weight:Weight), list(mean=mean, se=std.error), na.rm=T)

#png("resubmission/figure 5.png", width=3, height=3, units="in", res=300)
ggplot(data=testesAvg, aes(x=LightCycle, y=mean)) +
  geom_bar(stat="identity", color="black", fill="grey") +
  geom_point(data=testesData, aes(x=LightCycle, y=Weight), position=position_jitter(width=0.1), shape=1) +
  geom_errorbar(aes(ymin=mean-se, ymax=mean+se), 
                width=0.2, linetype="solid", position=position_dodge(width=.9)) +
  geom_segment(x=1, y=625, xend=2, yend=625) +
  geom_text(x=1.5, y=640, label="*", size=8) +
  ylab("Testes Weight (mg)") +
  xlab("Light Cycle") +
  scale_y_continuous(expand=c(0,0), limits=c(0,700), guide = guide_prism_minor()) +
  theme_classic() 
#dev.off()

t.test(Weight~LightCycle, data=testesData)



        
##### Resonance - before and after #####
beforeAfterData=read.table("initial v final data.csv", sep=",", header=T)
beforeAfterData$Time = factor(beforeAfterData$Time, levels=c("Pre-resonance", "Post-resonance"))
beforeAfterData$LightCycle = factor(beforeAfterData$LightCycle, levels=c("6:6", "4.5:7.5"))
beforeAfterAvg = beforeAfterData %>%
   group_by(LightCycle,Time) %>%
   summarize_at(vars(Weight:Comets), list(mean=mean, se=std.error), na.rm=T)

#weight
weight = ggplot(data=beforeAfterAvg, aes(x=Time, y=Weight_mean, group=LightCycle)) +
  geom_line(aes(linetype=LightCycle), position=position_dodge(width=0.15)) +
  geom_point(position=position_dodge(width=0.15), size=1) +
  geom_errorbar(aes(ymin=Weight_mean-Weight_se, 
                    ymax=Weight_mean+Weight_se), width=0.13, position=position_dodge(width=0.15)) +
  scale_linetype_manual(labels = c(expression(12:12 %->% 6:6), expression(12:12 %->% 4.5:7.5)), 
                        values = c("solid", "twodash"),
                        name = "Light Cycle") +
  scale_x_discrete(labels=c("Pre-\nresonance", "Post-\nresonance"), expand=c(.15,.15)) +
  scale_y_continuous(limits=c(20,27), name = "Weight (g)", guide = guide_prism_minor()) +
  theme_classic() +
  theme(legend.text.align = 0) +
  theme(axis.title.x = element_blank())

leveneTest(Weight ~ Time * LightCycle, data=beforeAfterData) #over .05 so were good
model = lmer(Weight ~ Time * LightCycle + (1|Bird), data=beforeAfterData)
Anova(model, type = "III", test.statistic = "F") #not sig



#cort
cort = ggplot(data=beforeAfterAvg, aes(x=Time, y=Cort_mean, group=LightCycle)) +
  geom_point(position=position_dodge(width=0.15), size=1) +  
  geom_line(aes(linetype=LightCycle), position=position_dodge(width=0.15)) +
  geom_errorbar(aes(ymin=Cort_mean-Cort_se, 
                    ymax=Cort_mean+Cort_se),width=0.13, position=position_dodge(width=0.15)) +
  geom_segment(aes(x=1, y=24.5, xend=2, yend=24.5)) +
  geom_segment(aes(x=1, y=23.5, xend=1, yend=24.5)) +
  geom_segment(aes(x=2, y=23.5, xend=2, yend=24.5)) +
  geom_text(x=1.5, y=25, label="*", size=6) +
  scale_linetype_manual(labels = c(expression(12:12 %->% 6:6), expression(12:12 %->% 4.5:7.5)), 
                        values = c("solid", "twodash"),
                        name = "Light Cycle") +
  scale_x_discrete(labels=c("Pre-\nresonance", "Post-\nresonance"), expand=c(.15,.15)) +
  scale_y_continuous(breaks=c(0,5,10,15,20,25), limits=c(0,26), 
                     name="Baseline Corticosterone (ng/mL)", guide = guide_prism_minor()) +
  theme_classic() +
  theme(legend.text.align = 0) +
  theme(axis.title.x = element_blank())
        
leveneTest(Cort ~ Time * LightCycle, data=beforeAfterData) #over .05 so were good
model = lmer(Cort ~ Time * LightCycle + (1|Bird), data=beforeAfterData)
Anova(model, type = "III", test.statistic = "F") #time is sig



#comets
comets = ggplot(data=beforeAfterAvg, aes(x=Time, y=Comets_mean, group=LightCycle)) +
  geom_point(position=position_dodge(width=0.15), size=1) +   
  geom_line(aes(linetype=LightCycle), position=position_dodge(width=0.15)) +
  geom_errorbar(aes(ymin=Comets_mean-Comets_se, 
                    ymax=Comets_mean+Comets_se),width=0.13, position=position_dodge(width=0.15)) +
  geom_segment(aes(x=1, y=40.5, xend=2, yend=40.5)) +
  geom_segment(aes(x=1, y=39, xend=1, yend=40.5)) +
  geom_segment(aes(x=2, y=39, xend=2, yend=40.5)) +
  geom_text(x=1.5, y=42, label="*", size=6) +
  scale_linetype_manual(labels = c(expression(12:12 %->% 6:6), expression(12:12 %->% 4.5:7.5)), 
                        values = c("solid", "twodash"),
                        name = "Light Cycle") +
  scale_x_discrete(labels=c("Pre-\nresonance", "Post-\nresonance"), expand=c(.15,.15)) +
  scale_y_continuous(name="DNA Damage (Tail Moment)", limits=c(0,44), guide = guide_prism_minor())+
  theme_classic() +
  theme(legend.text.align = 0) +
  theme(axis.title.x = element_blank())

leveneTest(Comets ~ Time * LightCycle, data=beforeAfterData) #over .05 so were good
model = lmer(Comets ~ Time * LightCycle + (1|Bird), data=beforeAfterData)
Anova(model, type = "III", test.statistic = "F") #time is sig


#png("figure 4.png", width=6, height=3, units="in", res=300)
weight + theme(legend.position = c(0.45, 0.21), text=element_text(size=10)) +
  theme(legend.key.size = unit(.5, 'cm'), #change legend key size
        legend.key.height = unit(.5, 'cm'), #change legend key height
        legend.key.width = unit(.5, 'cm'), #change legend key width
        legend.title = element_text(size=7.5), #change legend title font size
        legend.text = element_text(size=6.5)) + #change legend text font size 
  cort + theme(legend.position = "none", text=element_text(size=10)) + 
  comets+ theme(legend.position = "none", text=element_text(size=10)) +  
  plot_annotation(tag_levels = 'A')
#dev.off()