library(tidyverse)
library(lubridate)
library(dplyr)
library(tidyr)
library(ggplot2)
library("Rmisc")
library("plyr")
library("parallel")
library(ggpubr)
library(ggsignif)
rm(list = ls())
t5<-Sys.time()
#load raw data files
total.data<- read.csv(file = paste0(getwd(),'/total.data.csv'),header = TRUE)
genotypecsv.name<- list.files(paste0( getwd(),"/Input/Fly genotype/"))
cbPalette <- c("#999999", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7") #8 colors
#cbPalette <- c("#666666", "#999999", "#CC79A7")# 3 groups, mutation
#cbPalette <- c( "#999999", "#D55E00")# 2 groups, RNAi
#cbPalette <- c( "#999999", "#009E73")# 2 groups, R20D02
#cbPalette <- c( "#999999", "#0072B2")# 2 groups, R23E10
#cbPalette <- c("#999999","#E69F00","#56B4E9", "#CC79A7" )# 4 groups, pdfr5304;Dh31r10974, w c c e
#cbPalette <- c("#999999", "#999999","#009E73","#999999","#0072B2") #5 groups, RNAi, R20D02, R23E10
#cbPalette <- c("#999999","#E69F00","#999999", "#CC79A7", "#999999", "#F0E442")
#cbPalette <- c("#666666", "#999999", "#0072B2")# 3 groups, mutation
genotypecsv_path <- paste0(getwd(),'/Input/Fly genotype/',list.files(paste0( getwd(),"/Input/Fly genotype/")))
#path of genotype###.csv, it was not matter the formal of the name.
levelscsv_path <- paste0(getwd(),'/Input/Levels/',list.files(paste0( getwd(),"/Input/Levels/")))
f_genotype <- read.csv(file = genotypecsv_path,header = TRUE) #read csv for genotype
#If there are workspace.RData saved in "Input" folders
rdata_path <- paste0(getwd(),'/Input/Workspace/')
if (length(list.files(rdata_path) )== 1 ) {
load(list.files(rdata_path,full.names = TRUE))
}else {
f_levels <- read.csv(file= levelscsv_path,header = F) %>% .$V1 #read csv for genotype levels
print("Genotype arrange as")
print(f_levels)
f_levels
#experiment period (days)
period_exp <-2
#days after experiment
after_exp <- 0
#setting
total.data$e.time <- ymd_hms(total.data$e.time, tz="japan")
F_ID <- c(formatC(1:32 , width = 2, flag = "0"))
file_name <- list.files(paste0( getwd(),"/Input/Analysis data/"))
files_number <- length(file_name)
print("file list : ")
file_name
print('"............................."')
print('"............................."')
#``
#sleep judge function
sleep.judge <- function(x){
z <- x-5
y <- 5
sleep.c <- 0
if(x<=6){
z <- 1
}
if(x>=nr.e.data-5){
y <- nr.e.data-x
}
for (i in c(z:x)) {
if(sum(exp.monitor.data[c(i:(i+y)),'activity.logic'])==0){
sleep.c <- 1
}
}
sleep.c
}
#n_days rep function
rep_d <- function(x){
rep(x,length=1440)
}
rep_hr <- function(x){
rep(x,length=60)
}
# For loop. Load monitor data and calculate the sleep mins with activity.logic.
for (i in 1:files_number) {
monitor.path <- paste0(getwd(),"/Input/Analysis data/",file_name[i])
#show loading file name
print(paste("loading file ===> ",file_name[i]))
raw.monitor.data <- read.delim ( monitor.path, header=FALSE)%>%
.[c(-4,-5,-6,-7,-8,-9)] #delete the useless column
names(raw.monitor.data) <- c('NO.','date','time','light',F_ID ) #rename the title
raw.monitor.data$date <- gsub(" ","",raw.monitor.data$date)#delete all of the space in the date
raw.monitor.data$date <- dmy(raw.monitor.data$date) #transform date frame, need package "lubridate"
#raw.monitor.data$time <- as.integer(raw.monitor.data$time)
#raw.monitor.data$time <- as.data.frame(raw.monitor.data$time)
raw.monitor.data$time <- format( strptime(raw.monitor.data$time, format= "%H:%M:%S"),"%H:%M:%S")
e_date <- max(raw.monitor.data$date)-after_exp
s_date <- e_date- period_exp
s_times <- filter(raw.monitor.data, `date` == s_date&light==1)%>%.[,'time']%>%min(.)%>%paste(s_date, .)%>% ymd_hms(.,tz="japan")
e_times <- strftime(s_times,"%H:%M:%S")%>%paste(e_date,.)%>%ymd_hms(.,tz="japan")-60
raw.monitor.data <- mutate(raw.monitor.data, e.time =ymd_hms(paste(raw.monitor.data$date, raw.monitor.data$time, sep='_'),tz="japan"))%>%
filter(., e.time>=s_times&e.time<=e_times)%>%
mutate(n.mins= rep(1:1440,length= nrow(.) ) )%>%
mutate(n.hrs= ((n.mins-1)%/%60)+1)%>%
mutate(n.days= unlist(lapply(1:period_exp,rep_d)))
exp.monitor.data <- tidyr::gather(
raw.monitor.data,
key = F_ID,
value = activity,
-date,
-e.time,
-time,
-NO.,
-light,
-n.mins,
-n.hrs,
-n.days
)
n.monitor <- file_name[i]
exp.monitor.data <- add_column(exp.monitor.data,monitor = substr(n.monitor, 1, nchar(n.monitor)-4),.before = 1) # put the name of monitor files into DATA "monitor"
exp.monitor.data$monitor <- as.character(exp.monitor.data$monitor) # the data must be character for joining the 'monitor' of 'genotype###.csv'
exp.monitor.data$F_ID <- as.integer(exp.monitor.data$F_ID)
exp.monitor.data <- mutate(exp.monitor.data,activity.logic = if_else(activity == 0, 0, 1 ))%>%
arrange(monitor,F_ID,e.time)%>%
mutate(sleep = 0)
exp.monitor.data <- left_join(exp.monitor.data,f_genotype,by = c("monitor"="monitor","F_ID"="F_ID"))
nr.e.data <- nrow(exp.monitor.data)
t1<-Sys.time()
print("calculating... ")
exp.monitor.data[1:nr.e.data, 'sleep'] <- sapply(1:nr.e.data, sleep.judge)
#从5开始，算式中带-5, 无法为负值。计算sleep logic值，0无睡，1睡眠
#calculate must start from row 5, because it will -5 during the calculation, and the result cant below 0. Calculate the sleep logic value, 0 means no sleep, 1 means flies sleep during this min.
t2<-Sys.time()
rt<- t2-t1
print(paste("Calculation took ",round(rt,1)," seconds"))
total.data <- rbind(total.data,exp.monitor.data)
print('"............................."')
ifelse(i<files_number,print('"............next............."'),print('"...........finish............"'))
print('"............................."')
}
total.data$sleep <- as.integer(total.data$sleep)
total.data$activity <- as.integer(total.data$activity)
total.data$genotype <- as.character(total.data$genotype)
total.data <- mutate(total.data, F_NO= paste0(monitor,"_", F_ID))%>%mutate(exp_period= ifelse(light==1, "light","dark")) %>%
select( monitor,F_ID,F_NO,n.days,n.hrs,n.mins, light,exp_period, genotype, activity, activity.logic, sleep ) %>% filter(genotype != "")
n_genotype <- f_genotype%>%group_by(genotype)%>%dplyr::summarise(n_f= n())
n_genotype$genotype <- as.character(n_genotype$genotype)
sleep_bd_judge <- function(x){
#x <- number of genotype
bd_geno <- sleep_bd_list[x] %>% as.character()
bd_geno_df <- total.data %>%
filter(genotype== bd_geno)
list_F_ID <- bd_geno_df$F_NO %>% unlist() %>% unique()
n_F_ID <- length(list_F_ID)
for (i in 1:n_F_ID) {
# each flies of each genotypes
bd_fid <- list_F_ID[i]
df1 <- bd_geno_df %>% filter(F_NO == bd_fid)
sbd <- 0
for (j in (1:nrow(df1))) {
# calculate each data of sleep
k1 <- df1[j,'sleep']
if (k1 ==1) {
sbd <- sbd+1
}else{
k2 <- c(k2,df1[j,"genotype"])
k3 <- c(k3,sbd)
k4 <- c(k4,df1[j,"light"])
k5 <- c(k5,df1[j,'F_NO'])
sbd <- 0
}
}
k2 <- c(k2,df1[nrow(df1),"genotype"])
k3 <- c(k3,sbd)
k4 <- c(k4,df1[nrow(df1),"light"])
k5 <- c(k5,df1[nrow(df1),'F_NO'])
}
k2 <- unlist(k2)
k3 <- unlist(k3)
k4 <- unlist(k4)
k5 <- unlist(k5)
rbind(k2,k3,k4,k5) %>% as.data.frame()
}
k2 <- 'a' %>% as.character()
k3 <- 0 %>% as.integer()
k4 <- 0 %>% as.character()
k5 <- 0 %>% as.character()
sleep_bd_list <- total.data$genotype %>% unlist()%>% unique()
sleep_bd_n <- count(total.data$genotype) %>%nrow()
raw_bd<- sapply(1:sleep_bd_n, sleep_bd_judge) %>% as.data.frame() %>% t() #ERROR, K1 ==1. means data not found. there are something wrong before data input
print(paste("Calculation took ",round(rt,1)," seconds"))
k6 <- as.character(raw_bd[,1])
k7 <- as.integer(raw_bd[,2])
k8 <- as.character(raw_bd[,3])
k9 <- as.character(raw_bd[,4])
bd_frame <- data.frame(genotype= k6, sleep_bout_duration=k7,light=k8,F_NO=k9)
bd_frame <- filter(bd_frame, genotype!=0&sleep_bout_duration>=5) %>% mutate(exp_period= ifelse(light==1,'light','dark'))
df_dj <- bd_frame %>% group_by(F_NO) %>%
dplyr::summarise(
dead_judge_logic = ifelse(sleep_bout_duration > 720,1,0)
) %>% group_by(F_NO) %>%
dplyr::summarise(
sum_dj_logic= sum(dead_judge_logic)
) %>% mutate(dead_judge= ifelse(sum_dj_logic>0,"D","L")) %>%
select(F_NO,dead_judge)
bd_frame <- left_join(bd_frame,df_dj,by= 'F_NO')
total.data <- left_join(total.data,df_dj,by= 'F_NO')
n_genotype0<- total.data %>% select(genotype, F_NO,dead_judge) %>% unique()%>% filter(genotype!=""&dead_judge=="L") %>%
group_by(genotype,dead_judge) %>%
dplyr::summarise(n_f=n()) %>% select(genotype,n_f)
#change the levels of light and dark, for _facet
bd_frame$exp_period <- factor(bd_frame$exp_period, levels = c("light","dark"))
total.data$exp_period <- factor(total.data$exp_period, levels = c("light","dark"))
#change the arrange of different genotype, as the levles.csv
bd_frame$genotype <- factor(bd_frame$genotype, levels = f_levels)
total.data$genotype <- factor(total.data$genotype, levels = f_levels)
count_bd <- bd_frame %>% filter(dead_judge=='L') %>%
left_join(., n_genotype0,by= 'genotype') %>%
group_by(sleep_bout_duration,genotype,exp_period,n_f,dead_judge) %>%
dplyr::summarise(
n=n()/n_f
)
F_NO_list <- total.data$F_NO %>% unique()
}
## [1] "Genotype arrange as"
## [1] "R23E10/+" "R23E10>Dh31R-RNAi" "tim/+"
## [4] "tim>Dh31R-RNAi" "pdf/+" "pdf>v37763"
## [7] "AstA/+" "AstA>v37763"
## [1] "file list : "
## [1] "\".............................\""
## [1] "\".............................\""
## [1] "loading file ===> 220328_1.txt"
## [1] "calculating... "
## [1] "Calculation took 19 seconds"
## [1] "\".............................\""
## [1] "\"............next.............\""
## [1] "\".............................\""
## [1] "loading file ===> 220404_5.txt"
## [1] "calculating... "
## [1] "Calculation took 19 seconds"
## [1] "\".............................\""
## [1] "\"............next.............\""
## [1] "\".............................\""
## [1] "loading file ===> 220404_6.txt"
## [1] "calculating... "
## [1] "Calculation took 15.5 seconds"
## [1] "\".............................\""
## [1] "\"...........finish............\""
## [1] "\".............................\""
## [1] "Calculation took 15.5 seconds"
#raw data
head(raw.monitor.data)
#analysised data
head(total.data)
ld <- total.data %>% select(genotype, F_NO,dead_judge) %>% unique()%>% filter(genotype!="") %>%
group_by(genotype,dead_judge) %>%
dplyr::summarise(n=n())
mean_hr_sleep <- total.data%>%filter(dead_judge=='L') %>%
group_by(genotype,n.hrs,n.days)%>%
dplyr::summarise(n_f= n()/60,
sleep_sum_hr=sum(sleep),
sleep_per_hr=(sleep_sum_hr/n_f))
mean_hr_sleep_2 <- total.data%>%filter(dead_judge=='L') %>%
group_by(genotype,F_NO,n.hrs)%>%
dplyr::summarise(n_f= n()/60,
sleep_sum_hr=sum(sleep),
sleep_per_hr=(sleep_sum_hr/n_f),
sleep_per_hr_2=mean(sleep)*60,
)
mean_hr_sleep_error <- mean_hr_sleep_2 %>% group_by(genotype,n.hrs)%>%
dplyr::summarise(sleep_per_hr_total=mean(sleep_per_hr),
sleep_per_hr_sd= sd(sleep_per_hr),
sleep_per_hr_n = n(),
sleep_per_hr_sem= sleep_per_hr_sd/sqrt(sleep_per_hr_n)
)
g_ld <- ggplot(data= ld,mapping=aes(x=genotype,y=n,fill=dead_judge))+
geom_bar(stat="identity",width=0.5,position='fill')+
scale_fill_manual(values=c('#999999','#D55E00'))+
geom_text(stat='identity',aes(label=n), color="white", size=3.5,position=position_fill(0.5))+
labs(x="genotype",y = "N of flies",title = "Dead flies")+
theme_classic()+
theme(axis.title=element_text(size=14, lineheight=.9, family="Times"),
axis.text.x = element_text(size=12,angle = 45,hjust=1),
plot.title=element_text(size=14, lineheight=.9, family="Times"))
g_ld
print(ld)
## # A tibble: 4 × 3
## # Groups: genotype [2]
## genotype dead_judge n
## <fct> <chr> <int>
## 1 tim/+ D 1
## 2 tim/+ L 56
## 3 tim>Dh31R-RNAi D 2
## 4 tim>Dh31R-RNAi L 37
ld_s<- spread(ld, key= dead_judge, value = n) %>% as.tibble()
g_total_sleep_curve_df <- total.data%>% filter(dead_judge=='L') %>%
group_by(genotype,n.hrs)%>%
dplyr::summarise(n_f= n()/60,
sleep_sum_hr=sum(sleep),
sleep_per_hr=(sleep_sum_hr/n_f))
g_total_sleep_curve <-
ggplot(g_total_sleep_curve_df,aes(n.hrs,sleep_per_hr,color=genotype, fill= genotype))+
annotate("rect",xmin=-1,xmax=12,ymin=-1,ymax=0,fill="snow1",color="black",show.legend = TRUE) +
annotate("rect",xmin=12,xmax=25,ymin=-1,ymax=0,alpha=.9,fill="gray11",color="black")+
scale_x_continuous(breaks = scales::breaks_width(6))+
theme_classic()+
theme(legend.position = c(0.72, 0.2),
legend.key.height = unit(0.5, "cm"),
legend.key.width = unit(0.8, "cm")
)+
theme(legend.title = element_blank(),legend.text = element_text(size=6))+
geom_line(show.legend = TRUE)+
geom_vline(aes(xintercept= 12),colour="black", linetype=2)+
theme(axis.title=element_text(size=14, lineheight=.9, family="Times"),
axis.text.x = element_text(size=12),
plot.title=element_text(size=14, lineheight=.9, family="Times"))+
labs(x="ZT",y = "Sleep [per hr]",title = "Sleep curve of total")+
geom_point(show.legend = FALSE)+
scale_colour_manual(values=cbPalette)
g_total_sleep_curve_sem <-
ggplot(mean_hr_sleep_error,aes(n.hrs,sleep_per_hr_total,color=genotype, fill= genotype))+
annotate("rect",xmin=-1,xmax=12,ymin=-1,ymax=0,fill="snow1",color="black",show.legend = TRUE) +
annotate("rect",xmin=12,xmax=25,ymin=-1,ymax=0,alpha=.9,fill="gray11",color="black")+
scale_x_continuous(breaks = scales::breaks_width(6))+
theme_classic()+
theme(legend.position = c(0.72, 0.2),
legend.key.height = unit(0.5, "cm"),
legend.key.width = unit(0.8, "cm")
)+
theme(legend.title = element_blank(),legend.text = element_text(size=6))+
geom_line(show.legend = TRUE)+
geom_vline(aes(xintercept= 12),colour="black", linetype=2)+
theme(axis.title=element_text(size=14, lineheight=.9, family="Times"),
axis.text.x = element_text(size=12),
plot.title=element_text(size=14, lineheight=.9, family="Times"))+
labs(x="ZT",y = "Sleep [per hr]",title = "Sleep curve of total")+
geom_point(show.legend = FALSE)+
geom_errorbar(aes(ymin=sleep_per_hr_total-sleep_per_hr_sem, ymax=sleep_per_hr_total+sleep_per_hr_sem), width=.2) +
scale_colour_manual(values=cbPalette)
g_total_sleep_curve_sd <-
ggplot(mean_hr_sleep_error,aes(n.hrs,sleep_per_hr_total,color=genotype, fill= genotype))+
annotate("rect",xmin=-1,xmax=12,ymin=-1,ymax=0,fill="snow1",color="black",show.legend = TRUE) +
annotate("rect",xmin=12,xmax=25,ymin=-1,ymax=0,alpha=.9,fill="gray11",color="black")+
scale_x_continuous(breaks = scales::breaks_width(6))+
theme_classic()+
theme(legend.title = element_blank(),legend.text = element_text(size=6))+
geom_line(show.legend = TRUE)+
geom_vline(aes(xintercept= 12),colour="black", linetype=2)+
theme(axis.title=element_text(size=14, lineheight=.9, family="Times"),
axis.text.x = element_text(size=12),
plot.title=element_text(size=14, lineheight=.9, family="Times"))+
labs(x="ZT",y = "Sleep [per hr]",title = "Sleep curve of total")+
geom_point(show.legend = FALSE)+
geom_errorbar(aes(ymin=sleep_per_hr_total-sleep_per_hr_sd, ymax=sleep_per_hr_total+sleep_per_hr_sd), width=.2) +
scale_colour_manual(values=cbPalette)
#sleep curve in one day
g_total_sleep_curve
g_total_sleep_curve_sem
g_total_sleep_curve_sd
#sleep curve each day
g_sleep_curve_days<- mean_hr_sleep%>%
ggplot(aes(n.hrs,sleep_per_hr,color=genotype, fill= genotype))+
annotate("rect",xmin=1,xmax=12,ymin=-1,ymax=0,fill="snow1",color="black") +
annotate("rect",xmin=12,xmax=24,ymin=-1,ymax=0,alpha=.9,fill="gray11",color="black")+
scale_x_continuous(breaks = scales::breaks_width(3))+
geom_line()+
geom_point()+
geom_vline(aes(xintercept= 12),colour="black", linetype=2)+
theme_classic()+
theme(axis.title=element_text(size=14, lineheight=.9, family="Times"),
axis.text.x = element_text(size=12),
plot.title=element_text(size=14, lineheight=.9, family="Times"))+
labs(x="hrs",y = "Sleep [per hr]",title = "Sleep curve of each day")+
facet_grid(n.days~.)+
scale_colour_manual(values=cbPalette)
g_sleep_curve_days
#以10位单位，填充1到总行数
g_sleep_curve_days_vline1 <- data.frame(ZT=seq(24, (period_exp)*24, by=24) )
g_sleep_curve_days_v<- mean_hr_sleep%>%mutate(ZT = (n.days-1)*24+n.hrs) %>%
ggplot(aes(ZT,sleep_per_hr,color=genotype, fill= genotype))+
annotate("rect",xmin=g_sleep_curve_days_vline1$ZT-24,xmax=g_sleep_curve_days_vline1$ZT-12,ymin=-1,ymax=0,fill="snow1",color="black") +
annotate("rect",xmin=g_sleep_curve_days_vline1$ZT-12,xmax=g_sleep_curve_days_vline1$ZT,ymin=-1,ymax=0,alpha=.9,fill="gray11",color="black")+
geom_line()+
geom_point()+
geom_vline(aes(xintercept= ZT) ,data= g_sleep_curve_days_vline1, colour="grey" , linetype=4,size=.5 ,alpha=.9)+
theme_classic()+
xlim (0,period_exp*24)+
scale_x_continuous(breaks = scales::breaks_width(12))+
theme(axis.title=element_text(size=14, lineheight=.9, family="Times"),
axis.text.x = element_text(size=12),
plot.title=element_text(size=14, lineheight=.9, family="Times"))+
labs(x="ZT",y = "Sleep [per hr]",title = "Sleep curve of each day")+
scale_colour_manual(values=cbPalette)
## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.
g_sleep_curve_days_v
g_total_sleep_bar<- total.data%>%filter(dead_judge=='L') %>% group_by(genotype,monitor,F_ID)%>%dplyr::summarise(total_sleep= sum(sleep)/period_exp)%>%
ggplot(aes(genotype, total_sleep,color=genotype,fill = genotype))+
theme_classic()+
theme(axis.title=element_text(size=14, lineheight=.9, family="Times"),
axis.text.x = element_text(size=12,angle = 45,hjust=1),
plot.title=element_text(size=14, lineheight=.9, family="Times"))+
stat_summary(fun = "mean", fun.min = function(x)mean(x), fun.max = function(x)mean(x), geom = "bar",colour="grey", width = 0.5,alpha= 0.8,show.legend = FALSE)+
stat_summary(fun.min = function(x)mean(x) - sd(x), fun.max = function(x)mean(x) + sd(x), geom = "errorbar", colour = "black", width = 0.2) +
geom_jitter(color= "black",alpha=.5,size=1.5,width = .05,show.legend = FALSE)+
labs(x="",y = "sleep [mins/day]",title = "Total sleep")+
scale_fill_manual(values=cbPalette)
#total sleep
#total.data%>% group_by(genotype,monitor,F_ID)%>%dplyr::summarise(total_sleep= sum(sleep)/period_exp)
avo_total_sleep <- total.data%>%filter(dead_judge=='L') %>%
group_by(genotype,F_NO)%>%dplyr::summarise(total_sleep=sum(sleep)/period_exp)%>%select(.,genotype,total_sleep)%>%
ungroup()%>%aov(total_sleep~genotype,data=.)%>%TukeyHSD()#One-way ANOVA
g_total_sleep_bar
avo_total_sleep
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = total_sleep ~ genotype, data = .)
##
## $genotype
## diff lwr upr p adj
## tim>Dh31R-RNAi-tim/+ 18.1945 -24.84574 61.23474 0.4032741
sum_day_sleep<- total.data%>%filter(dead_judge=='L') %>%
group_by(genotype,monitor,F_ID,exp_period)%>%dplyr::summarise(total_sleep= sum(sleep)/period_exp)
#total sleep of light and dark period
g_total_sleep_period_bar<- sum_day_sleep%>%
ggplot(aes(genotype, total_sleep, fill=genotype))+
stat_summary(fun = "mean", fun.min = function(x)mean(x), fun.max = function(x)mean(x), geom = "bar",colour="grey", width = 0.5,alpha= 0.9,show.legend = FALSE)+
stat_summary(fun.min = function(x)mean(x) - sd(x), fun.max = function(x)mean(x) + sd(x), geom = "errorbar", colour = "black", width = 0.2) +
geom_jitter(color= "black",alpha=.5,size=1.5,width = .1,show.legend = FALSE)+
labs(x="",y = "sleep [mins/day]",title = "Total sleep by period")+
theme_classic()+
theme(axis.title=element_text(size=14, lineheight=.9, family="Times"),
axis.text.x = element_text(size=12,angle = 45,hjust=1),
plot.title=element_text(size=14, lineheight=.9, family="Times"))+
facet_grid(.~exp_period) +
scale_fill_manual(values=cbPalette)
#One-way ANOVA, exp_period = light
avo_total_sleep_period_light <- total.data%>%filter(dead_judge=='L') %>%
filter(exp_period == 'light') %>%
group_by(genotype,exp_period,F_NO)%>%dplyr::summarise(total_sleep=sum(sleep)/period_exp)%>%
select(.,genotype,total_sleep)%>%
ungroup()%>%aov(total_sleep~genotype,data=.,)%>%TukeyHSD()
#One-way ANOVA, exp_period = dark
avo_total_sleep_period_dark <- total.data%>%filter(dead_judge=='L') %>%
filter(exp_period == 'dark') %>%
group_by(genotype,exp_period,F_NO)%>%dplyr::summarise(total_sleep=sum(sleep)/period_exp)%>%
select(.,genotype,total_sleep)%>%
ungroup()%>%aov(total_sleep~genotype,data=.)%>%TukeyHSD()
g_total_sleep_period_bar
print('One-way ANOVA, total_sleep, exp_period = light')
## [1] "One-way ANOVA, total_sleep, exp_period = light"
avo_total_sleep_period_light
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = total_sleep ~ genotype, data = .)
##
## $genotype
## diff lwr upr p adj
## tim>Dh31R-RNAi-tim/+ 45.35907 17.69367 73.02447 0.001584
print('One-way ANOVA, total_sleep, exp_period = dark')
## [1] "One-way ANOVA, total_sleep, exp_period = dark"
avo_total_sleep_period_dark
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = total_sleep ~ genotype, data = .)
##
## $genotype
## diff lwr upr p adj
## tim>Dh31R-RNAi-tim/+ -27.16458 -50.26529 -4.063864 0.0216978
core_exp_data <- total.data %>% filter(dead_judge=="L") %>% filter(n.hrs>=4&n.hrs<=20) %>% filter(n.hrs<=8|n.hrs>=16)
core_exp_data_df <- core_exp_data %>% group_by(genotype,monitor,F_ID,exp_period)%>%dplyr::summarise(total_sleep= sum(sleep)/period_exp)
#core sleep of total
g_core_sleep_bar<- core_exp_data%>%filter(dead_judge=='L') %>% group_by(genotype,monitor,F_ID)%>%dplyr::summarise(total_sleep= sum(sleep)/period_exp)%>%
ggplot(aes(genotype, total_sleep,color=genotype,fill = genotype))+
theme_classic()+
theme(axis.title=element_text(size=14, lineheight=.9, family="Times"),
axis.text.x = element_text(size=12,angle = 45,hjust=1),
plot.title=element_text(size=14, lineheight=.9, family="Times"))+
stat_summary(fun = "mean", fun.min = function(x)mean(x), fun.max = function(x)mean(x), geom = "bar",colour="grey", width = 0.5,alpha= 0.8,show.legend = FALSE)+
stat_summary(fun.min = function(x)mean(x) - sd(x), fun.max = function(x)mean(x) + sd(x), geom = "errorbar", colour = "black", width = 0.2) +
geom_jitter(color= "black",alpha=.5,size=1.5,width = .05,show.legend = FALSE)+
labs(x="",y = "sleep [mins per day]",title = "Core time sleep (4-9hrs, 16-21hrs, total)")+
scale_fill_manual(values=cbPalette)
avo_core_sleep <- total.data%>%filter(dead_judge=='L') %>%
group_by(genotype,F_NO)%>%dplyr::summarise(total_sleep=sum(sleep)/period_exp)%>%select(.,genotype,total_sleep)%>%
ungroup()%>%aov(total_sleep~genotype,data=.)%>%TukeyHSD()#One-way ANOVA
#core sleep of light and dark period
g_core_sleep_period_bar<- core_exp_data%>%filter(dead_judge=='L') %>% group_by(genotype,monitor,exp_period, F_ID)%>%
dplyr::summarise(total_sleep= sum(sleep)/period_exp)%>%
ggplot(aes(genotype, total_sleep, fill=genotype))+
stat_summary(fun = "mean", fun.min = function(x)mean(x), fun.max = function(x)mean(x), geom = "bar",colour="grey", width = 0.5,alpha= 0.9,show.legend = FALSE)+
stat_summary(fun.min = function(x)mean(x) - sd(x), fun.max = function(x)mean(x) + sd(x), geom = "errorbar", colour = "black", width = 0.2) +
geom_jitter(color= "black",alpha=.5,size=1.5,width = .1,show.legend = FALSE)+
labs(x="",y = "sleep [mins per day]",title = "Core time sleep (4-9hrs, 16-21hrs, by period)")+
theme_classic()+
theme(axis.title=element_text(size=14, lineheight=.9, family="Times"),
axis.text.x = element_text(size=12,angle = 45,hjust=1),
plot.title=element_text(size=14, lineheight=.9, family="Times"))+
facet_grid(.~exp_period) +
scale_fill_manual(values=cbPalette)
#One-way ANOVA, exp_period = light
avo_core_sleep_period_light <- core_exp_data%>%filter(dead_judge=='L') %>%
filter(exp_period == 'light') %>%
group_by(genotype,exp_period,F_NO)%>%dplyr::summarise(total_sleep=sum(sleep)/period_exp)%>%
select(.,genotype,total_sleep)%>%
ungroup()%>%aov(total_sleep~genotype,data=.,)%>%TukeyHSD()
#One-way ANOVA, exp_period = dark
avo_core_sleep_period_dark <- core_exp_data%>%filter(dead_judge=='L') %>%
filter(exp_period == 'dark') %>%
group_by(genotype,exp_period,F_NO)%>%dplyr::summarise(total_sleep=sum(sleep)/period_exp)%>%
select(.,genotype,total_sleep)%>%
ungroup()%>%aov(total_sleep~genotype,data=.)%>%TukeyHSD()
g_core_sleep_bar
print('One-way ANOVA, core_sleep_total')
## [1] "One-way ANOVA, core_sleep_total"
avo_core_sleep
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = total_sleep ~ genotype, data = .)
##
## $genotype
## diff lwr upr p adj
## tim>Dh31R-RNAi-tim/+ 18.1945 -24.84574 61.23474 0.4032741
g_core_sleep_period_bar
print('One-way ANOVA, core_sleep, exp_period = dark')
## [1] "One-way ANOVA, core_sleep, exp_period = dark"
avo_core_sleep_period_dark
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = total_sleep ~ genotype, data = .)
##
## $genotype
## diff lwr upr p adj
## tim>Dh31R-RNAi-tim/+ -23.08398 -32.15483 -14.01312 2.2e-06
print('One-way ANOVA, core_sleep, exp_period = light')
## [1] "One-way ANOVA, core_sleep, exp_period = light"
avo_core_sleep_period_light
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = total_sleep ~ genotype, data = .)
##
## $genotype
## diff lwr upr p adj
## tim>Dh31R-RNAi-tim/+ 5.290058 -2.919813 13.49993 0.2038235
dat_vlines <- count_bd %>%
group_by(genotype,exp_period) %>%
dplyr::summarise(
bd_mean= mean(sleep_bout_duration),
F_NO=n_f,
n_sbd=n(),
n_sbd_perday=(n()/F_NO)/period_exp
) %>% unique()
sbd_lim <- max(count_bd$sleep_bout_duration)*.2
n_lim <- max(count_bd$n)*.3
g_sbd_point_v <- count_bd %>%
ggplot(.,aes(n,sleep_bout_duration))+
geom_point(aes(color= genotype ),alpha=.5,size=2,show.legend = FALSE)+
labs(x="n of sleep bout duration [each fly]",y = "Sleep bout duration [mins]",title = "Sleep bout duration ") +
coord_cartesian(ylim = c(0, sbd_lim),xlim = c(0, n_lim))+
theme_classic()+
theme(axis.title=element_text(size=14, lineheight=.9, family="Times"),
axis.text.x = element_text(size=12),
plot.title=element_text(size=14, lineheight=.9, family="Times"))+
geom_hline(aes(yintercept= bd_mean),data= dat_vlines,colour="#990000", linetype="dashed")+
geom_text(aes(x=n_lim*0.6,y=bd_mean+7, label= round(bd_mean,1)),data = dat_vlines, color="#990000", size=3)+
geom_text(aes(x=n_lim*0.6,y=.9*sbd_lim,hjust = 0,label= paste0("Flies number=",round(F_NO,1))),data = dat_vlines, color="black", size=3)+
geom_text(aes(x=n_lim*0.6,y=.82*sbd_lim,hjust = 0, label= paste0("Sbd number=",round(n_sbd,1))),data = dat_vlines, color="black", size=3)+
geom_text(aes(x=n_lim*0.6,y=.74*sbd_lim,hjust = 0, label= paste0("Sbd n per fly=",round(n_sbd_perday,1))),data = dat_vlines, color="black", size=3)+
facet_grid(exp_period~genotype)+
scale_colour_manual(values=cbPalette)
g_sbd_point_v_total<- count_bd %>%
ggplot(aes(n,sleep_bout_duration))+
geom_point(aes(color= genotype ),alpha=.5,size=2,show.legend = FALSE)+
geom_hline(aes(yintercept= bd_mean),data= dat_vlines,colour="#990000", linetype="dashed")+
labs(x="n of sleep bout duration [each fly]",y = "Sleep bout duration [mins]",title = "Sleep bout duration ") +
theme_classic()+
theme(axis.title=element_text(size=14, lineheight=.9, family="Times"),
axis.text.x = element_text(size=12),
plot.title=element_text(size=14, lineheight=.9, family="Times"))+
geom_text(aes(x=2.5,y=bd_mean+25, label= round(bd_mean,1)),data = dat_vlines, color="#990000", size=3.2)+
geom_text(aes(x=.6*max(count_bd$n),y=.9*max(count_bd$sleep_bout_duration),hjust = 0,label= paste0("Flies number=",round(F_NO,1))),data = dat_vlines, color="black", size=3)+
geom_text(aes(x=.6*max(count_bd$n),y=.82*max(count_bd$sleep_bout_duration),hjust = 0, label= paste0("Sbd number=",round(n_sbd,1))),data = dat_vlines, color="black", size=3)+
geom_text(aes(x=.6*max(count_bd$n),y=.74*max(count_bd$sleep_bout_duration),hjust = 0, label= paste0("Sbd n per fly=",round(n_sbd_perday,1))),data = dat_vlines, color="black", size=3)+
facet_grid(exp_period~genotype)+
scale_colour_manual(values=cbPalette)
g_sbd_point_h <- count_bd%>%
ggplot(aes(sleep_bout_duration,n))+
geom_point(aes( color= genotype ),alpha=.5,size=2,show.legend = FALSE)+
geom_vline(aes(xintercept= bd_mean),data= dat_vlines,colour="#990000", linetype="dashed")+
labs(y="n of sleep bout duration [each fly]",x = "Sleep bout duration [mins]",title = "Sleep bout duration ") +
theme_classic()+
theme(axis.title=element_text(size=14, lineheight=.9, family="Times"),
axis.text.x = element_text(size=12),
plot.title=element_text(size=14, lineheight=.9, family="Times"))+
coord_cartesian(xlim = c(0, sbd_lim),ylim = c(0, n_lim))+
geom_text(aes(y=.82*n_lim,x=bd_mean+4, label= round(bd_mean,1)),data = dat_vlines, color="#990000", size=3)+
geom_text(aes(x=sbd_lim*0.6,y=.9*n_lim,hjust = 0,label= paste0("Flies number=",round(F_NO,1))),data = dat_vlines, color="black", size=3)+
geom_text(aes(x=sbd_lim*0.6,y=.82*n_lim,hjust = 0, label= paste0("Sbd number=",round(n_sbd,1))),data = dat_vlines, color="black", size=3)+
geom_text(aes(x=sbd_lim*0.6,y=.74*n_lim,hjust = 0, label= paste0("Sbd n per fly=",round(n_sbd_perday,1))),data = dat_vlines, color="black", size=3)+
facet_grid(genotype~exp_period)+
scale_colour_manual(values=cbPalette)
g_sbd_point_h_total <- count_bd%>%
ggplot(aes(sleep_bout_duration,n))+
geom_point(aes( color= genotype ),alpha=.5,size=2,show.legend = FALSE)+
geom_vline(aes(xintercept= bd_mean),data= dat_vlines,colour="#990000", linetype="dashed")+
labs(y="n of sleep bout duration [each fly]",x = "Sleep bout duration [mins]",title = "Sleep bout duration ") +
theme_classic()+
theme(axis.title=element_text(size=14, lineheight=.9, family="Times"),
axis.text.x = element_text(size=12),
plot.title=element_text(size=14, lineheight=.9, family="Times"))+
geom_text(aes(y=4,x=bd_mean+25, label= round(bd_mean,1)),data = dat_vlines, color="#990000", size=3)+
geom_text(aes(x=.6*max(count_bd$sleep_bout_duration),y=.9*max(count_bd$n),hjust = 0,label= paste0("Flies number=",round(F_NO,1))),data = dat_vlines, color="black", size=3)+
geom_text(aes(x=.6*max(count_bd$sleep_bout_duration),y=.82*max(count_bd$n),hjust = 0, label= paste0("Sbd number=",round(n_sbd,1))),data = dat_vlines, color="black", size=3)+
geom_text(aes(x=.6*max(count_bd$sleep_bout_duration),y=.74*max(count_bd$n),hjust = 0, label= paste0("Sbd n per fly=",round(n_sbd_perday,1))),data = dat_vlines, color="black", size=3)+
facet_grid(genotype~exp_period)+
scale_colour_manual(values=cbPalette)
g_sbd_point_v
g_sbd_point_h
dat_vlines1 <- count_bd %>%
group_by(genotype) %>%
dplyr::summarise(
bd_mean= mean(sleep_bout_duration)
)
max_sbd<- max(count_bd$sleep_bout_duration)
g_sbd_total <- count_bd %>%
ggplot(aes(x=genotype,y= sleep_bout_duration,color=genotype))+
geom_jitter(alpha=.3,size=1.1,show.legend = FALSE)+
stat_summary(fun = "mean", fun.min = function(x)mean(x), fun.max = function(x)mean(x), geom = "crossbar", colour = "red", width = 0.5)+
stat_summary(fun.min = function(x)mean(x) - sd(x), fun.max = function(x)mean(x) + sd(x), geom = "errorbar", colour = "black", width = 0.2) +
labs(x="",y = "Sleep bout duration [mins]",title = "Total sleep bout duration ") +
theme_classic()+
theme(axis.title=element_text(size=14, lineheight=.9, family="Times"),
axis.text.x = element_text(size=12,angle = 45,hjust=1),
plot.title=element_text(size=14, lineheight=.9, family="Times"))+
geom_text(aes(x=genotype,y=bd_mean+20, label= round(bd_mean,1)),nudge_x=.2,data = dat_vlines1, color="#990000", size=3)+
coord_cartesian(ylim = c(-0.1*max_sbd, 0.5* max_sbd))+
scale_colour_manual(values=cbPalette)
g_sbd_total_period <- count_bd %>%
ggplot(aes(x=genotype,y= sleep_bout_duration,color=genotype))+
geom_jitter(alpha=.3,size=1.8,show.legend = FALSE)+
stat_summary(fun = "mean", fun.min = function(x)mean(x), fun.max = function(x)mean(x), geom = "crossbar", colour = "red", width = 0.5)+
stat_summary(fun.min = function(x)mean(x) - sd(x), fun.max = function(x)mean(x) + sd(x), geom = "errorbar", colour = "black", width = 0.2) +
labs(x="",y = "Sleep bout duration [mins]",title = "Total sleep bout duration ") +
theme_classic()+
theme(axis.title=element_text(size=14, lineheight=.9, family="Times"),
axis.text.x = element_text(size=12,angle = 45,hjust=1),
plot.title=element_text(size=14, lineheight=.9, family="Times"))+
geom_text(aes(x=genotype,y=bd_mean+20, label= round(bd_mean,1)),nudge_x=.2,data = dat_vlines, color="#990000", size=3)+
coord_cartesian(ylim = c(-0.1*max_sbd, 0.5* max_sbd))+
facet_grid(.~exp_period)+
scale_colour_manual(values=cbPalette)
#One-way ANOVA, total sleep bout duration
avo_sbd_total <- bd_frame%>%filter(dead_judge=='L') %>%
aov(sleep_bout_duration~genotype,data=.,)%>%TukeyHSD()
avo_sbd_total_light <- bd_frame %>% filter(dead_judge=='L'&exp_period=="light") %>%
aov(sleep_bout_duration~genotype,data=.,)%>%TukeyHSD()
avo_sbd_total_dark <- bd_frame %>% filter(dead_judge=='L'&exp_period=="dark") %>%
aov(sleep_bout_duration~genotype,data=.,)%>%TukeyHSD()
#m <- lm(log(Weight) ~ Drought.Treatment*Plant.Species, larval.data)
#anova(m) #Two-way ANOVA
g_sbd_total
print('One-way ANOVA, total_sbd')
## [1] "One-way ANOVA, total_sbd"
avo_sbd_total
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = sleep_bout_duration ~ genotype, data = .)
##
## $genotype
## diff lwr upr p adj
## tim>Dh31R-RNAi-tim/+ -8.329877 -12.58087 -4.07888 0.000124
g_sbd_total_period
print('One-way ANOVA, total_sbd, exp_period = light')
## [1] "One-way ANOVA, total_sbd, exp_period = light"
avo_sbd_total_light
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = sleep_bout_duration ~ genotype, data = .)
##
## $genotype
## diff lwr upr p adj
## tim>Dh31R-RNAi-tim/+ 3.55214 -2.837576 9.941856 0.2757652
print('One-way ANOVA, total_sbd, exp_period = dark')
## [1] "One-way ANOVA, total_sbd, exp_period = dark"
avo_sbd_total_dark
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = sleep_bout_duration ~ genotype, data = .)
##
## $genotype
## diff lwr upr p adj
## tim>Dh31R-RNAi-tim/+ -23.81914 -29.24275 -18.39552 0
g_sb_total_df <- bd_frame %>% filter(dead_judge=='L') %>%
group_by(genotype,F_NO) %>%
dplyr::summarise(
n_perday= n()/period_exp)
g_sb_total_period_df <- bd_frame %>% filter(dead_judge=='L') %>%
group_by(genotype,F_NO,exp_period) %>%
dplyr::summarise(
n_perday= n()/period_exp)
g_sb_total_bar <- g_sb_total_df %>%
ggplot(aes(genotype, n_perday,color=genotype,fill = genotype))+
theme_classic()+
theme(axis.title=element_text(size=14, lineheight=.9, family="Times"),
axis.text.x = element_text(size=12,angle = 45,hjust=1),
plot.title=element_text(size=14, lineheight=.9, family="Times"))+
stat_summary(fun = "mean", fun.min = function(x)mean(x), fun.max = function(x)mean(x), geom = "bar",colour="grey", width = 0.5,alpha= 0.8,show.legend = FALSE)+
stat_summary(fun.min = function(x)mean(x) - sd(x), fun.max = function(x)mean(x) + sd(x), geom = "errorbar", colour = "black", width = 0.2) +
geom_jitter(color= "black",alpha=.5,size=1.5,width = .05,show.legend = FALSE)+
labs(x="",y = "Amount of sleep bout",title = "Sleep bout per day")+
scale_fill_manual(values=cbPalette)
avo_sb_total <- g_sb_total_df %>% select(.,genotype,n_perday)%>%
ungroup()%>%aov(n_perday~genotype,data=.)%>%TukeyHSD()#One-way ANOVA
g_sb_total_period_bar <- g_sb_total_period_df %>%
ggplot(aes(genotype, n_perday,color=genotype,fill = genotype))+
theme_classic()+
theme(axis.title=element_text(size=14, lineheight=.9, family="Times"),
axis.text.x = element_text(size=12,angle = 45,hjust=1),
plot.title=element_text(size=14, lineheight=.9, family="Times"))+
stat_summary(fun = "mean", fun.min = function(x)mean(x), fun.max = function(x)mean(x), geom = "bar",colour="grey", width = 0.5,alpha= 0.8,show.legend = FALSE)+
stat_summary(fun.min = function(x)mean(x) - sd(x), fun.max = function(x)mean(x) + sd(x), geom = "errorbar", colour = "black", width = 0.2) +
geom_jitter(color= "black",alpha=.5,size=1.5,width = .05,show.legend = FALSE)+
labs(x="",y = "Amount of sleep bout",title = "Sleep bout per day")+
facet_grid(.~exp_period)+
scale_fill_manual(values=cbPalette)
avo_sb_total <- g_sb_total_df %>% select(.,genotype,n_perday)%>%
ungroup()%>%aov(n_perday~genotype,data=.)%>%TukeyHSD()#One-way ANOVA
avo_sb_total_period_light <- g_sb_total_period_df %>%filter(exp_period=="light") %>% select(.,genotype,n_perday)%>%
ungroup()%>%aov(n_perday~genotype,data=.)%>%TukeyHSD()#One-way ANOVA
avo_sb_total_period_dark <- g_sb_total_period_df %>%filter(exp_period=="dark") %>% select(.,genotype,n_perday)%>%
ungroup()%>%aov(n_perday~genotype,data=.)%>%TukeyHSD()#One-way ANOVA
g_sb_total_bar
print("avo_sb_total")
## [1] "avo_sb_total"
avo_sb_total
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = n_perday ~ genotype, data = .)
##
## $genotype
## diff lwr upr p adj
## tim>Dh31R-RNAi-tim/+ 4.872587 2.050716 7.694457 0.0009093
g_sb_total_period_bar
print("avo_sb_total_period_light")
## [1] "avo_sb_total_period_light"
avo_sb_total_period_light
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = n_perday ~ genotype, data = .)
##
## $genotype
## diff lwr upr p adj
## tim>Dh31R-RNAi-tim/+ 0.001206564 -2.086502 2.088915 0.9990865
print("avo_sb_total_period_dark")
## [1] "avo_sb_total_period_dark"
avo_sb_total_period_dark
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = n_perday ~ genotype, data = .)
##
## $genotype
## diff lwr upr p adj
## tim>Dh31R-RNAi-tim/+ 4.87138 3.487628 6.255133 0
g_sbd_sb_total_df <- bd_frame %>% filter(dead_judge=="L") %>% left_join(., n_genotype0,genotype=genotype) %>%
group_by(genotype, F_NO) %>%
dplyr::summarise(
n_bd_day=n()/period_exp,
mean_sbd=mean(sleep_bout_duration)
)
g_sbd_sb_period_df <- bd_frame %>% filter(dead_judge=="L") %>% left_join(., n_genotype0,genotype=genotype) %>%
group_by(genotype, F_NO, exp_period) %>%
dplyr::summarise(
n_bd_day=n()/period_exp,
mean_sbd=mean(sleep_bout_duration)
)
g_sbd_sb_total_0<- g_sbd_sb_total_df %>%
ggplot(aes(n_bd_day,mean_sbd,color=genotype))+
theme_classic()+
theme(axis.title=element_text(size=14, lineheight=.9, family="Times"),
axis.text.x = element_text(size=12,angle = 45,hjust=1),
plot.title=element_text(size=14, lineheight=.9, family="Times"))+
geom_point()+
labs(x="sleep bout",y = "sleep bout duration",title = "Sleep bout and sleep bout duration (per day)")+
geom_smooth(alpha= .2)+
scale_colour_manual(values=cbPalette)
g_sbd_sb_total_1<- g_sbd_sb_total_df %>%
ggplot(aes(1/n_bd_day,mean_sbd,color=genotype))+
theme_classic()+
theme(axis.title=element_text(size=14, lineheight=.9, family="Times"),
axis.text.x = element_text(size=12,angle = 45,hjust=1),
plot.title=element_text(size=14, lineheight=.9, family="Times"))+
geom_point()+
labs(x="1/sleep bout",y = "sleep bout duration",title = "Sleep bout and sleep bout duration (per day)")+
geom_smooth(alpha= .2,method="lm")+
scale_colour_manual(values=cbPalette)
g_sbd_sb_total_0
g_sbd_sb_total_1
n_sum <- 10 #以10位单位
seq.total.data <- total.data %>% filter(genotype!=""&dead_judge=="L")
grp <- seq(1, nrow(seq.total.data), by=n_sum) #以10位单位，填充1到总行数
seq_sleep_sum <- function(x){ #读取数据表中SLEEP数据，以n_sum为单位计算sleep的和。**计算每10分钟睡眠
sum(seq.total.data[x:(x+(n_sum-1)),"sleep"])
}
seq_sleep_sum_c<- sapply(c(grp), seq_sleep_sum) #循环计算
seq_total.data <- seq.total.data[grp,] %>%
mutate(
seq_sleep_sum= seq_sleep_sum_c ,
NO=rep(1:(1440*period_exp/n_sum), length= nrow(.))) %>%
select( NO,F_NO,exp_period, genotype, dead_judge,seq_sleep_sum )
#dat_vlines2 <- seq_total.data %>%filter(dead_judge=="D") %>% #死亡数据的标识横线 **未使用
#select(genotype,F_NO) %>% mutate(NO_m=720*period_exp) %>% unique()
dat_vlines3_c <- seq(0,144*period_exp,by=72) #半日标识线，全时间段每720分钟一根
dat_vlines3 <- data.frame(NO= dat_vlines3_c)
dat_vlines3_c1 <- seq(0,144*period_exp-1,by=144) #nighttime -> daytime 切换时标识线，全时间段每1440分钟一根
dat_vlines4 <- data.frame(NO= dat_vlines3_c1)
g_sbd_tile<- seq_total.data %>% filter(genotype!=""&dead_judge=="L") %>%
ggplot() +
geom_tile(aes(x=NO,y=F_NO,fill = seq_sleep_sum))+
scale_fill_gradient(low ="steelblue" ,high = "white")+
facet_grid(genotype~., scales = "free", space = "free_y")+
labs(x="",y = "",title = "Activity tile ") +
theme(panel.grid =element_blank())+
theme(axis.title=element_text(size=16, lineheight=.9, family="Times"),
axis.text.x = element_blank(),
plot.title=element_text(size=16, lineheight=.9, family="Times"),
legend.title = element_text( size=10, face="bold"))+
#geom_hline(aes(yintercept= F_NO),data= dat_vlines2,colour="black", linetype=1.2,size=2.37,alpha=.5)+
geom_vline(aes(xintercept= NO) ,data= dat_vlines3 , colour="red" , linetype=4,size=.5 ,alpha=.9)+
annotate("rect",xmin=dat_vlines4$NO,xmax=dat_vlines4$NO+72,ymin=-1,ymax=0,fill="snow1",color="black") +
annotate("rect",xmin=dat_vlines4$NO+72,xmax=dat_vlines4$NO+144,ymin=-1,ymax=0,
alpha=.9,fill="gray11",color="black")
plot(g_sbd_tile)
#sleep bout duration tile
dat_vlines5_c <- seq(0,24,by=12) #半日标识线，全时间段每720分钟一根
dat_vlines5 <- data.frame(NO= dat_vlines5_c)
dat_vlines5_c1 <- seq(0,24,by=24) #nighttime -> daytime 切换时标识线，全时间段每1440分钟一根
dat_vlines6 <- data.frame(NO= dat_vlines5_c1)
g_sleepbout_tile<- total.data %>% filter(genotype!=""&dead_judge=="L", n.days== period_exp) %>%
ggplot() +
geom_tile(aes(x=n.mins/60,y=F_NO,fill = sleep))+
scale_fill_gradient(low ="white" ,high = "darkorange3")+
facet_grid(genotype~., scales = "free", space = "free_y")+
labs(x="",y = "",title = "Activity tile ") +
theme(panel.grid =element_blank())+
scale_x_continuous(breaks = scales::breaks_width(6))+
theme(axis.title=element_text(size=16, lineheight=.9, family="Times"),
axis.text.y = element_blank(),
plot.title=element_text(size=16, lineheight=.9, family="Times"),
legend.title = element_text( size=10, face="bold"))+
#geom_hline(aes(yintercept= F_NO),data= dat_vlines2,colour="black", linetype=1.2,size=2.37,alpha=.5)+
geom_vline(aes(xintercept= NO) ,data= dat_vlines5 , colour="red" , linetype=4,size=.5 ,alpha=.9)+
annotate("rect",xmin=0,xmax=12,ymin=-1,ymax=0,fill="snow1",color="black") +
labs(x="Zeitgeber time(h)",y = "Fly NO.",title = "Sleep bout tile in 24hr")+
annotate("rect",xmin=12,xmax=24,ymin=-1,ymax=0,
alpha=.9,fill="gray11",color="black")
g_sleepbout_tile
g_waking_activities_df <- total.data%>%filter(dead_judge=='L', activity.logic=='1') %>% group_by(genotype,F_NO,n.days)%>%
dplyr::summarise(waking_activities_d= sum(activity)/n()) %>%
group_by(genotype,F_NO)%>%dplyr::summarise(waking_activities= mean(waking_activities_d))
g_waking_activities_period_df <- total.data%>%filter(dead_judge=='L', activity.logic=='1') %>% group_by(genotype,F_NO,n.days,exp_period)%>%
dplyr::summarise(waking_activities_d= sum(activity)/n()) %>%
group_by(genotype,F_NO,exp_period)%>%dplyr::summarise(waking_activities= mean(waking_activities_d))
g_waking_activities_bar<- ggplot(g_waking_activities_df,aes(genotype, waking_activities,color=genotype,fill = genotype))+
theme_classic()+
theme(axis.title=element_text(size=14, lineheight=.9, family="Times"),
axis.text.x = element_text(size=12,angle = 45,hjust=1),
plot.title=element_text(size=14, lineheight=.9, family="Times"))+
stat_summary(fun = "mean", fun.min = function(x)mean(x), fun.max = function(x)mean(x), geom = "bar",colour="grey", width = 0.5,alpha= 0.8,show.legend = FALSE)+
stat_summary(fun.min = function(x)mean(x) - sd(x), fun.max = function(x)mean(x) + sd(x), geom = "errorbar", colour = "black", width = 0.2) +
geom_jitter(color= "black",alpha=.5,size=1.5,width = .05,show.legend = FALSE)+
labs(x="",y = "Waking activities [amount/mins]",title = "Waking activities")+
scale_fill_manual(values=cbPalette)
g_waking_activities_period_bar<- ggplot(g_waking_activities_period_df,aes(genotype, waking_activities,color=genotype,fill = genotype))+
theme_classic()+
theme(axis.title=element_text(size=14, lineheight=.9, family="Times"),
axis.text.x = element_text(size=12,angle = 45,hjust=1),
plot.title=element_text(size=14, lineheight=.9, family="Times"))+
stat_summary(fun = "mean", fun.min = function(x)mean(x), fun.max = function(x)mean(x), geom = "bar",colour="grey", width = 0.5,alpha= 0.8,show.legend = FALSE)+
stat_summary(fun.min = function(x)mean(x) - sd(x), fun.max = function(x)mean(x) + sd(x), geom = "errorbar", colour = "black", width = 0.2) +
geom_jitter(color= "black",alpha=.5,size=1.5,width = .05,show.legend = FALSE)+
labs(x="",y = "Waking activities [amount/mins]",title = "Waking activities")+
facet_grid(.~exp_period)+
scale_fill_manual(values=cbPalette)
#One-way ANOVA, total
avo_waking_activities <- g_waking_activities_df %>%
ungroup()%>%aov(waking_activities~genotype,data=.)%>%TukeyHSD()
#One-way ANOVA, exp_period = light
avo_waking_activities_period_light <- g_waking_activities_period_df %>% filter(exp_period == 'light') %>%
ungroup()%>%aov(waking_activities~genotype,data=.,)%>%TukeyHSD()
#One-way ANOVA, exp_period = dark
avo_waking_activities_period_dark <- g_waking_activities_period_df %>% filter(exp_period == 'dark') %>%
ungroup()%>%aov(waking_activities~genotype,data=.,)%>%TukeyHSD()
g_waking_activities_bar
g_waking_activities_period_bar
print('One-way ANOVA, Waking activities, total')
## [1] "One-way ANOVA, Waking activities, total"
avo_waking_activities
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = waking_activities ~ genotype, data = .)
##
## $genotype
## diff lwr upr p adj
## tim>Dh31R-RNAi-tim/+ 0.2720328 0.08250437 0.4615613 0.0053906
print('One-way ANOVA, Waking activities, exp_period = light')
## [1] "One-way ANOVA, Waking activities, exp_period = light"
avo_waking_activities_period_light
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = waking_activities ~ genotype, data = .)
##
## $genotype
## diff lwr upr p adj
## tim>Dh31R-RNAi-tim/+ -0.1616276 -0.3672061 0.0439509 0.1218282
print('One-way ANOVA, Waking activities, exp_period = dark')
## [1] "One-way ANOVA, Waking activities, exp_period = dark"
avo_waking_activities_period_dark
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = waking_activities ~ genotype, data = .)
##
## $genotype
## diff lwr upr p adj
## tim>Dh31R-RNAi-tim/+ 0.4762083 0.2742729 0.6781438 9.8e-06
# core time waking activities
#
g_c_waking_activities_df <- total.data%>%filter(dead_judge=='L', activity.logic=='1',n.hrs>=4&n.hrs<=20,n.hrs<=8|n.hrs>=16) %>% group_by(genotype,F_NO,n.days)%>%
dplyr::summarise(waking_activities_d= sum(activity)/n()) %>%
group_by(genotype,F_NO)%>%dplyr::summarise(waking_activities= mean(waking_activities_d))
g_c_waking_activities_period_df <- total.data%>%filter(dead_judge=='L', activity.logic=='1',n.hrs>=4&n.hrs<=20,n.hrs<=8|n.hrs>=16) %>% group_by(genotype,F_NO,n.days,exp_period)%>%
dplyr::summarise(waking_activities_d= sum(activity)/n()) %>%
group_by(genotype,F_NO,exp_period)%>%dplyr::summarise(waking_activities= mean(waking_activities_d))
g_c_waking_activities_bar<- ggplot(g_c_waking_activities_df,aes(genotype, waking_activities,color=genotype,fill = genotype))+
theme_classic()+
theme(axis.title=element_text(size=14, lineheight=.9, family="Times"),
axis.text.x = element_text(size=12,angle = 45,hjust=1),
plot.title=element_text(size=14, lineheight=.9, family="Times"))+
stat_summary(fun = "mean", fun.min = function(x)mean(x), fun.max = function(x)mean(x), geom = "bar",colour="grey", width = 0.5,alpha= 0.8,show.legend = FALSE)+
stat_summary(fun.min = function(x)mean(x) - sd(x), fun.max = function(x)mean(x) + sd(x), geom = "errorbar", colour = "black", width = 0.2) +
geom_jitter(color= "black",alpha=.5,size=1.5,width = .05,show.legend = FALSE)+
labs(x="",y = "Waking activities [amount/mins]",title = "Core time Waking activities")+
scale_fill_manual(values=cbPalette)
g_c_waking_activities_period_bar<- ggplot(g_c_waking_activities_period_df,aes(genotype, waking_activities,color=genotype,fill = genotype))+
theme_classic()+
theme(axis.title=element_text(size=14, lineheight=.9, family="Times"),
axis.text.x = element_text(size=12,angle = 45,hjust=1),
plot.title=element_text(size=14, lineheight=.9, family="Times"))+
stat_summary(fun = "mean", fun.min = function(x)mean(x), fun.max = function(x)mean(x), geom = "bar",colour="grey", width = 0.5,alpha= 0.8,show.legend = FALSE)+
stat_summary(fun.min = function(x)mean(x) - sd(x), fun.max = function(x)mean(x) + sd(x), geom = "errorbar", colour = "black", width = 0.2) +
geom_jitter(color= "black",alpha=.5,size=1.5,width = .05,show.legend = FALSE)+
labs(x="",y = "Waking activities [amount/mins]",title = "Core time Waking activities")+
facet_grid(.~exp_period)+
scale_fill_manual(values=cbPalette)
#One-way ANOVA, total
avo_c_waking_activities <- g_c_waking_activities_df %>%
ungroup()%>%aov(waking_activities~genotype,data=.)%>%TukeyHSD()
#One-way ANOVA, exp_period = light
avo_c_waking_activities_period_light <- g_c_waking_activities_period_df %>% filter(exp_period == 'light') %>%
ungroup()%>%aov(waking_activities~genotype,data=.,)%>%TukeyHSD()
#One-way ANOVA, exp_period = dark
avo_c_waking_activities_period_dark <- g_c_waking_activities_period_df %>% filter(exp_period == 'dark') %>%
ungroup()%>%aov(waking_activities~genotype,data=.,)%>%TukeyHSD()
g_c_waking_activities_bar
g_c_waking_activities_period_bar
print('One-way ANOVA, Waking activities, total')
## [1] "One-way ANOVA, Waking activities, total"
avo_c_waking_activities
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = waking_activities ~ genotype, data = .)
##
## $genotype
## diff lwr upr p adj
## tim>Dh31R-RNAi-tim/+ 0.3771843 0.08153697 0.6728317 0.0129796
print('One-way ANOVA, Waking activities, exp_period = light')
## [1] "One-way ANOVA, Waking activities, exp_period = light"
avo_c_waking_activities_period_light
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = waking_activities ~ genotype, data = .)
##
## $genotype
## diff lwr upr p adj
## tim>Dh31R-RNAi-tim/+ -0.6435497 -1.215618 -0.07148178 0.0279233
print('One-way ANOVA, Waking activities, exp_period = dark')
## [1] "One-way ANOVA, Waking activities, exp_period = dark"
avo_c_waking_activities_period_dark
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = waking_activities ~ genotype, data = .)
##
## $genotype
## diff lwr upr p adj
## tim>Dh31R-RNAi-tim/+ 0.6914302 0.4027443 0.9801161 7.3e-06
# wacking activity of each hr
g_waking_activity_curve_df <- total.data%>% filter(dead_judge=='L') %>%
group_by(genotype,n.hrs)%>%
dplyr::summarise(n_f= n()/60,
activity_sum_hr=sum(activity),
activity_per_hr=(activity_sum_hr/n_f),
sleep_sum_hr=sum(sleep),
sleep_per_hr=(sleep_sum_hr/n_f),
activity_time=(60-sleep_per_hr),
waking_activity_per_hr=(activity_per_hr/(activity_time)))
g_waking_activity_curve <-
ggplot(g_waking_activity_curve_df,aes(n.hrs,waking_activity_per_hr,color=genotype, fill= genotype))+
annotate("rect",xmin=-1,xmax=12,ymin=-.1,ymax=0,fill="snow1",color="black",show.legend = TRUE) +
annotate("rect",xmin=12,xmax=25,ymin=-.1,ymax=0,alpha=.9,fill="gray11",color="black")+
scale_x_continuous(breaks = scales::breaks_width(6))+
theme_classic()+
theme(legend.title = element_blank(),legend.text = element_text(size=6))+
geom_line(show.legend = TRUE)+
geom_vline(aes(xintercept= 12),colour="black", linetype=2)+
theme(axis.title=element_text(size=14, lineheight=.9, family="Times"),
axis.text.x = element_text(size=12),
plot.title=element_text(size=14, lineheight=.9, family="Times"))+
labs(x="ZT",y = "Waking activity [per hr]",title = "Waking Activity curve of total")+
geom_point(show.legend = FALSE)+
scale_colour_manual(values=cbPalette)
g_activity_curve <-
ggplot(g_waking_activity_curve_df,aes(n.hrs,activity_per_hr,color=genotype, fill= genotype))+
annotate("rect",xmin=-1,xmax=12,ymin=-10,ymax=0,fill="snow1",color="black",show.legend = TRUE) +
annotate("rect",xmin=12,xmax=25,ymin=-10,ymax=0,alpha=.9,fill="gray11",color="black")+
scale_x_continuous(breaks = scales::breaks_width(6))+
theme_classic()+
theme(legend.title = element_blank(),legend.text = element_text(size=6))+
geom_line(show.legend = TRUE)+
geom_vline(aes(xintercept= 12),colour="black", linetype=2)+
theme(axis.title=element_text(size=14, lineheight=.9, family="Times"),
axis.text.x = element_text(size=12),
plot.title=element_text(size=14, lineheight=.9, family="Times"))+
labs(x="ZT",y = "Activity [per hr]",title = "Activity curve of total")+
geom_point(show.legend = FALSE)+
scale_colour_manual(values=cbPalette)
#sleep curve in one day
g_activity_time_curve <-
ggplot(g_waking_activity_curve_df,aes(n.hrs,activity_time,color=genotype, fill= genotype))+
annotate("rect",xmin=-1,xmax=12,ymin=-10,ymax=0,fill="snow1",color="black",show.legend = TRUE) +
annotate("rect",xmin=12,xmax=25,ymin=-10,ymax=0,alpha=.9,fill="gray11",color="black")+
scale_x_continuous(breaks = scales::breaks_width(6))+
theme_classic()+
theme(legend.title = element_blank(),legend.text = element_text(size=6))+
geom_line(show.legend = TRUE)+
geom_vline(aes(xintercept= 12),colour="black", linetype=2)+
theme(axis.title=element_text(size=14, lineheight=.9, family="Times"),
axis.text.x = element_text(size=12),
plot.title=element_text(size=14, lineheight=.9, family="Times"))+
labs(x="ZT",y = "Activity time [per hr]",title = "Activity time curve of total")+
geom_point(show.legend = FALSE)+
scale_colour_manual(values=cbPalette)
g_waking_activity_curve
g_activity_curve
g_activity_time_curve
g_sbd_bar_df <- bd_frame%>%filter(dead_judge=='L') %>% group_by(genotype,F_NO)%>%
dplyr::summarise(sbd_mean= mean(sleep_bout_duration))
## `summarise()` has grouped output by 'genotype'. You can override using the
## `.groups` argument.
g_sbd_period_bar_df <- bd_frame%>%filter(dead_judge=='L') %>% group_by(genotype,F_NO,exp_period)%>%
dplyr::summarise(sbd_mean= mean(sleep_bout_duration))
## `summarise()` has grouped output by 'genotype', 'F_NO'. You can override using
## the `.groups` argument.
g_sbd_bar <- ggplot(g_sbd_bar_df,aes(genotype, sbd_mean,color=genotype,fill = genotype))+
theme_classic()+
theme(axis.title=element_text(size=14, lineheight=.9, family="Times"),
axis.text.x = element_text(size=12,angle = 45,hjust=1),
plot.title=element_text(size=14, lineheight=.9, family="Times"))+
stat_summary(fun = "mean", fun.min = function(x)mean(x), fun.max = function(x)mean(x), geom = "bar",colour="grey", width = 0.5,alpha= 0.8,show.legend = FALSE)+
stat_summary(fun.min = function(x)mean(x) - sd(x), fun.max = function(x)mean(x) + sd(x), geom = "errorbar", colour = "black", width = 0.2) +
geom_jitter(color= "black",alpha=.5,size=1.5,width = .05,show.legend = FALSE)+
labs(x="",y = "Sleep bout duration [mins]",title = "Sleep bout duration Bar")+
scale_fill_manual(values=cbPalette)
g_sbd_period_bar <- ggplot(g_sbd_period_bar_df,aes(genotype, sbd_mean,color=genotype,fill = genotype))+
theme_classic()+
theme(axis.title=element_text(size=14, lineheight=.9, family="Times"),
axis.text.x = element_text(size=12,angle = 45,hjust=1),
plot.title=element_text(size=14, lineheight=.9, family="Times"))+
stat_summary(fun = "mean", fun.min = function(x)mean(x), fun.max = function(x)mean(x), geom = "bar",colour="grey", width = 0.5,alpha= 0.8,show.legend = FALSE)+
stat_summary(fun.min = function(x)mean(x) - sd(x), fun.max = function(x)mean(x) + sd(x), geom = "errorbar", colour = "black", width = 0.2) +
geom_jitter(color= "black",alpha=.5,size=1.5,width = .05,show.legend = FALSE)+
labs(x="",y = "Sleep bout duration [mins]",title = "Sleep bout duration Bar")+
facet_grid(.~exp_period)+
scale_fill_manual(values=cbPalette)
#One-way ANOVA, total
avo_sbd_bar <- g_sbd_bar_df %>%
ungroup()%>%aov(sbd_mean~genotype,data=.)%>%TukeyHSD()
#One-way ANOVA, exp_period = light
avo_sbd_period_bar_light <- g_sbd_period_bar_df %>% filter(exp_period == 'light') %>%
ungroup()%>%aov(sbd_mean~genotype,data=.,)%>%TukeyHSD()
#One-way ANOVA, exp_period = dark
avo_sbd_period_bar_dark <- g_sbd_period_bar_df %>% filter(exp_period == 'dark') %>%
ungroup()%>%aov(sbd_mean~genotype,data=.,)%>%TukeyHSD()
g_sbd_bar
g_sbd_period_bar
print('One-way ANOVA, Sleep bout duration, total')
## [1] "One-way ANOVA, Sleep bout duration, total"
avo_sbd_bar
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = sbd_mean ~ genotype, data = .)
##
## $genotype
## diff lwr upr p adj
## tim>Dh31R-RNAi-tim/+ -12.09186 -20.41343 -3.770284 0.0048673
print('One-way ANOVA, Sleep bout duration, exp_period = light')
## [1] "One-way ANOVA, Sleep bout duration, exp_period = light"
avo_sbd_period_bar_light
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = sbd_mean ~ genotype, data = .)
##
## $genotype
## diff lwr upr p adj
## tim>Dh31R-RNAi-tim/+ 0.9577368 -11.7432 13.65867 0.8812647
print('One-way ANOVA, Sleep bout duration, exp_period = dark')
## [1] "One-way ANOVA, Sleep bout duration, exp_period = dark"
avo_sbd_period_bar_dark
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = sbd_mean ~ genotype, data = .)
##
## $genotype
## diff lwr upr p adj
## tim>Dh31R-RNAi-tim/+ -31.4664 -41.81904 -21.11377 0
#for merge analysis
g_total_sleep_curve_df1 <- total.data%>% filter(dead_judge=='L') %>%
group_by(genotype ,F_NO,light)%>%
dplyr::summarise(sleep_sum_hr=sum(sleep)/(60*period_exp))
## `summarise()` has grouped output by 'genotype', 'F_NO'. You can override using
## the `.groups` argument.
g_total_sleep_curve_df2 <- total.data%>% filter(dead_judge=='L') %>%
group_by(genotype ,F_NO,light)%>%
dplyr::summarise(sleep_sum_d=sum(sleep)/(period_exp))
## `summarise()` has grouped output by 'genotype', 'F_NO'. You can override using
## the `.groups` argument.
g1 <- ggarrange(ggarrange(g_total_sleep_curve, ncol=2, nrow=1,legend = "right"),
g_sleep_curve_days_v,ncol=1, nrow=2,common.legend = T, legend = "none"
)
g1
g2 <- ggarrange(g_total_sleep_bar,g_sb_total_bar,g_sbd_bar,g_waking_activities_bar,
g_total_sleep_period_bar,g_sb_total_period_bar,g_sbd_period_bar,g_waking_activities_period_bar,
ncol=4, nrow=2,common.legend = T, legend = "bottom")
g2
g3 <- ggarrange(
ggarrange(g_total_sleep_curve, g_total_sleep_bar,g_sb_total_bar,g_sbd_bar, ncol=4, nrow=1,legend = "right"),
ncol=1, nrow=2,common.legend = T, legend = "none"
)
g3
#output path
getwd()
## [1] "/Users/leelv/Documents/Programing/R/DAM analysis R v1.9"
now_time <- Sys.time() %>% as.POSIXct()
today_date <- paste0(year(now_time),month(now_time),day(now_time),hour(now_time),minute(now_time))
output_path <- paste0(getwd(),"/Output/",today_date)
dir.create(output_path)
#total data frame output
write.csv(total.data, paste0(output_path,"/","total_data.csv"), col.names=TRUE,row.names = FALSE)
#dead and live flies message output.
ggsave(file = paste0(output_path,"/","1_Survived_Flies.png"), plot = g_ld, dpi = 200, width = 7, height = 5)
write.csv(ld_s, paste0(output_path,"/","Survived_Flies.csv"), col.names=TRUE,row.names = FALSE)
#sleep curve output
ggsave(file = paste0(output_path,"/","2_total_sleep_curve.png"), plot = g_total_sleep_curve_sem, dpi = 500, width = 7, height = 5)
ggsave(file = paste0(output_path,"/","3_Sleep_curve_days.png"), plot = g_sleep_curve_days, dpi = 400, width = 7, height = 7)
write.csv(g_total_sleep_curve_df, paste0(output_path,"/","sleep_total.csv"), col.names=TRUE,row.names = FALSE)
write.csv(g_total_sleep_curve_df1, paste0(output_path,"/","sleep_dn_total.csv"), col.names=TRUE,row.names = FALSE)
write.csv(g_total_sleep_curve_df2, paste0(output_path,"/","sleep_dn_total_d.csv"), col.names=TRUE,row.names = FALSE)
write.csv(mean_hr_sleep, paste0(output_path,"/","sleep_days.csv"), col.names=TRUE,row.names = FALSE)
write.csv(mean_hr_sleep_error, paste0(output_path,"/","sleep_days_sd.csv"), col.names=TRUE,row.names = FALSE)
#sleep bar output
ggsave(file = paste0(output_path,"/","4_total_sleep_period_bar.png"), plot = g_total_sleep_period_bar, dpi = 200, width = 4, height = 5)
ggsave(file = paste0(output_path,"/","5_total_sleep_bar.png"), plot = g_total_sleep_bar, dpi = 200, width = 2.5, height = 5)
#sleep bout duration output
ggsave(file = paste0(output_path,"/","6_sbd_point_h.png"), plot = g_sbd_point_h, dpi = 200, width = 10, height = 6)
ggsave(file = paste0(output_path,"/","7_sbd_point_h_total.png"), plot = g_sbd_point_h_total, dpi = 200, width = 10, height = 6)
ggsave(file = paste0(output_path,"/","8_sbd_point_v.png"), plot = g_sbd_point_v, dpi = 200, width = 10, height = 6)
ggsave(file = paste0(output_path,"/","9_sbd_point_v_total.png"), plot = g_sbd_point_v_total, dpi = 200, width = 10, height = 6)
ggsave(file = paste0(output_path,"/","10_sbd_total.png"), plot = g_sbd_total, dpi = 200, width = 5, height = 5)
ggsave(file = paste0(output_path,"/","11_sbd_total_period.png"), plot = g_sbd_total_period, dpi = 200, width = 5, height = 5)
write.csv(bd_frame, paste0(output_path,"/","sbd_data.csv"), col.names=TRUE,row.names = FALSE)
#sleep bout output
ggsave(file = paste0(output_path,"/","12_sleep_bout_bar.png"), plot = g_sb_total_bar, dpi = 200, width = 2.5, height = 5)
ggsave(file = paste0(output_path,"/","13_sleep_bout_period_bar.png"), plot = g_sb_total_period_bar, dpi = 200, width = 4, height = 5)
write.csv(g_sb_total_df, paste0(output_path,"/","sleep_bout_bar_data.csv"), col.names=TRUE,row.names = FALSE)
write.csv(g_sb_total_period_df, paste0(output_path,"/","sleep_bout_period_bar_data.csv"), col.names=TRUE,row.names = FALSE)
#core sleep bar output
ggsave(file = paste0(output_path,"/","14_core_sleep_bar.png"), plot = g_core_sleep_bar, dpi = 200, width = 5, height = 5)
ggsave(file = paste0(output_path,"/","15_core_sleep_period_bar.png"), plot = g_core_sleep_period_bar, dpi = 200, width = 5, height = 5)
write.csv(core_exp_data_df, paste0(output_path,"/","core_exp_data.csv"), col.names=TRUE,row.names = FALSE)
#sleep tile output
ggsave(file = paste0(output_path,"/","16_sbd_tile.png"), plot = g_sbd_tile, dpi = 200, width = 10, height = 8)
write.csv(seq_total.data, paste0(output_path,"/","seq_total_data.csv"), col.names=TRUE,row.names = FALSE)
#sleep bout and sleep bout duration output
ggsave(file = paste0(output_path,"/","17_sbd_sb_total_0.png"), plot = g_sbd_sb_total_0, dpi = 200, width = 10, height = 7.5)
ggsave(file = paste0(output_path,"/","18_sbd_sb_total_1.png"), plot = g_sbd_sb_total_1, dpi = 200, width = 10, height = 7.5)
#waking activity
ggsave(file = paste0(output_path,"/","19_waking_activities_bar.png"), plot = g_waking_activities_bar, dpi = 500, width = 2.5, height = 5)
ggsave(file = paste0(output_path,"/","20_waking_activities_period_bar.png"), plot = g_waking_activities_period_bar, dpi = 200, width = 4, height = 5)
write.csv(g_waking_activities_df, paste0(output_path,"/","waking_activities_total.csv"), col.names=TRUE,row.names = FALSE)
write.csv(g_waking_activities_period_df, paste0(output_path,"/","waking_activities_period.csv"), col.names=TRUE,row.names = FALSE)
#sleep bout duration bar
ggsave(file = paste0(output_path,"/","21_sbd_bar_bar.png"), plot = g_sbd_bar, dpi = 200, width = 2.5, height = 5)
ggsave(file = paste0(output_path,"/","22_sbd_period_bar.png"), plot = g_sbd_period_bar, dpi = 200, width = 4, height = 5)
ggsave(file = paste0(output_path,"/","23_sleepbout_tile.png"), plot = g_sleepbout_tile, dpi = 200, width = 10, height = 8)
write.csv(g_sbd_bar_df, paste0(output_path,"/","sbd_bar.csv"), col.names=TRUE,row.names = FALSE)
write.csv(g_sbd_period_bar_df, paste0(output_path,"/","sbd_period_bar.csv"), col.names=TRUE,row.names = FALSE)
#ANOVA analysis
write.csv(avo_total_sleep[[1]], paste0(output_path,"/","oneway_ANOVA_1_sleep.CSV"), col.names=TRUE,row.names = TRUE)
write.csv(avo_total_sleep_period_dark[[1]], paste0(output_path,"/","oneway_ANOVA_2_sleep_dark.CSV"), col.names=TRUE,row.names = TRUE)
write.csv(avo_total_sleep_period_light[[1]], paste0(output_path,"/","oneway_ANOVA_3_sleep_light.CSV"), col.names=TRUE,row.names = TRUE)
write.csv(avo_sb_total[[1]], paste0(output_path,"/","oneway_ANOVA_4_sleep_bout.CSV"), col.names=TRUE,row.names = TRUE)
write.csv(avo_sb_total_period_dark[[1]], paste0(output_path,"/","oneway_ANOVA_5_sleep_bout_dark.CSV"), col.names=TRUE,row.names = TRUE)
write.csv(avo_sb_total_period_light[[1]], paste0(output_path,"/","oneway_ANOVA_6_sleep_bout_light.CSV"), col.names=TRUE,row.names = TRUE)
write.csv(avo_sbd_total[[1]], paste0(output_path,"/","oneway_ANOVA_7_sleep_bout_duration.CSV"), col.names=TRUE,row.names = TRUE)
write.csv(avo_sbd_total_dark[[1]], paste0(output_path,"/","oneway_ANOVA_8_sleep_bout_duration_dark.CSV"), col.names=TRUE,row.names = TRUE)
write.csv(avo_sbd_total_light[[1]], paste0(output_path,"/","oneway_ANOVA_9_sleep_bout_duration_light.CSV"), col.names=TRUE,row.names = TRUE)
write.csv(avo_waking_activities[[1]], paste0(output_path,"/","oneway_ANOVA_10_waking_activities.CSV"), col.names=TRUE,row.names = TRUE)
write.csv(avo_waking_activities_period_light[[1]], paste0(output_path,"/","oneway_ANOVA_11_waking_activities_period_light.CSV"), col.names=TRUE,row.names = TRUE)
write.csv(avo_waking_activities_period_dark[[1]], paste0(output_path,"/","oneway_ANOVA_12_waking_activities_period_dark.CSV"), col.names=TRUE,row.names = TRUE)
write.csv(avo_sbd_bar[[1]], paste0(output_path,"/","oneway_ANOVA_13_sbd_bar.CSV"), col.names=TRUE,row.names = TRUE)
write.csv(avo_sbd_period_bar_light[[1]], paste0(output_path,"/","oneway_ANOVA_14_sbd_period_bar_light.CSV"), col.names=TRUE,row.names = TRUE)
write.csv(avo_sbd_period_bar_dark[[1]], paste0(output_path,"/","oneway_ANOVA_15_sbd_period_period_dark.CSV"), col.names=TRUE,row.names = TRUE)
#***core time calculation
write.csv(avo_core_sleep[[1]], paste0(output_path,"/","oneway_ANOVA_16_core_sleep.CSV"), col.names=TRUE,row.names = TRUE)
write.csv(avo_core_sleep_period_dark[[1]], paste0(output_path,"/","oneway_ANOVA_17_core_sleep_dark.CSV"), col.names=TRUE,row.names = TRUE)
write.csv(avo_core_sleep_period_light[[1]], paste0(output_path,"/","oneway_ANOVA_18_core_sleep_light.CSV"), col.names=TRUE,row.names = TRUE)
# merge graph
ggsave(file = paste0(output_path,"/","g1.png"), plot = g1, dpi = 200, width = 16, height = 9)
ggsave(file = paste0(output_path,"/","g2.png"), plot = g2, dpi = 200, width = 12, height = 9)
ggsave(file = paste0(output_path,"/","g3.png"), plot = g3, dpi = 400, width = 20, height = 10)
#***core time waking activity
ggsave(file = paste0(output_path,"/","23_coretime_waking_activities_bar.png"), plot = g_c_waking_activities_bar, dpi = 500, width = 2.5, height = 5)
ggsave(file = paste0(output_path,"/","24_coretime_waking_activities_period_bar.png"), plot = g_c_waking_activities_period_bar, dpi = 200, width = 4, height = 5)
write.csv(avo_c_waking_activities[[1]], paste0(output_path,"/","oneway_ANOVA_19_c_waking_activities.CSV"), col.names=TRUE,row.names = TRUE)
write.csv(avo_c_waking_activities_period_light[[1]], paste0(output_path,"/","oneway_ANOVA_20_c_waking_activities_period_light.CSV"), col.names=TRUE,row.names = TRUE)
write.csv(avo_c_waking_activities_period_dark[[1]], paste0(output_path,"/","oneway_ANOVA_21_c_waking_activities_period_dark.CSV"), col.names=TRUE,row.names = TRUE)
# waking activity & activity curve
ggsave(file = paste0(output_path,"/","25_total_waking_activity_curve.png"), plot = g_waking_activity_curve, dpi = 500, width = 4.5, height = 3.2)
ggsave(file = paste0(output_path,"/","26_total_activity_curve.png"), plot = g_activity_curve, dpi = 500, width = 4.5, height = 3.2)
output_path
## [1] "/Users/leelv/Documents/Programing/R/DAM analysis R v1.9/Output/20224142320"
#file copy
dir.create(paste0(output_path,"/Input/"))
dir.create(paste0(output_path,"/Input/Workspace/"))
save.image(file= paste0(output_path,"/Input/Workspace/Workspace_",today_date,".RData"))
dir.create(paste0(output_path,"/Input/Analysis data/"))
for (i in 1:files_number) {
file.copy(from = paste0(getwd(),"/Input/Analysis data/",file_name[i]),
to = paste0(output_path,"/Input/Analysis data/",file_name[i]))
}
#copy genotype.csv
dir.create(paste0(output_path,"/Input/Fly genotype"))
file.copy(from = genotypecsv_path,
to = paste0(output_path,"/Input/Fly genotype/genotype_",today_date,".csv"))
## [1] TRUE
#copy levels.csv
dir.create(paste0(output_path,"/Input/Levels"))
file.copy(from = levelscsv_path,
to = paste0(output_path,"/Input/Levels/Levels_",today_date,".csv"))
## [1] TRUE
t6<-Sys.time()
difftime(t6,t5)
## Time difference of 3.520257 mins