setwd("~/Desktop")
library(ggplot2)
library(broom)
mut1=read.table(text="weight length feeding_amount feeding_duration
1.381 1.093 1.386 1.273
1.313 0.805 1.427 1.307
1.2 0.965 1.231 1.186
1.277 1.101 1.282 1.499
1.285 1.005 1.36 0.923
1.263 0.815 1.332 1.138
1.34 0.937 1.22 0.932
1.428 1.072 1.392 1.832
1.031 1.075 1.202 1.192
1.333 1.105 1.29 1.091
1.126 1.197 1.209 1.234
1.325 0.942 1.204 1.287
1.226 1.161 1.373 1.21
1.396 1.054 1.174 1.199
1.125 0.988 1.464 1.05
1.287 1.106 1.328 1.694
1.205 1.098 1.331 0.934
1.137 1.069 1.523 1.932
1.079 1.083 1.101 1.166
1.245 1.021 1.18 1.171
1.362 0.982 1.359 1.076
1.207 0.911 1.46 0.939
1.116 1.015 1.142 1.177
1.214 1.047 1.397 1.312
1.173 1.093 1.221 1.141
1.176 1.085 1.344 1.443
1.209 1.302 1.387 0.932
1.165 1.225 1.09 0.968
1.19 1.137 1.545 1.152
1.189 1.097 1.396 1.157",header=T,sep="\t")
mut2=read.table(text="weight length feeding_amount feeding_duration
1.093 1.031 1.392 0.729
1.046 1.108 0.99 0.893
1.103 0.925 0.489 0.765
0.901 0.975 0.698 1.189
0.904 0.82 0.892 1.129
0.893 1.008 1.051 1.038
1.06 0.811 0.949 0.793
1.039 0.9 0.633 1.468
1.076 0.851 0.946 1.4
1.124 0.799 1.135 0.834
1.105 0.893 0.902 1.323
1.164 1.07 1.096 1.221
1.304 1.037 1.082 0.712
0.736 1.007 1.088 1.233
0.955 0.925 0.802 1.32
1.076 1.086 0.693 1.118
1.118 1.115 0.783 1.3
0.97 0.787 0.884 0.984
1.031 0.852 0.925 1.247
1.28 1.089 0.665 0.707
0.936 0.943 0.778 0.722
0.926 1.114 0.779 1.028
1.106 1.004 1.1 0.932
1.144 1.011 1.113 1.355
1.005 1.035 0.819 0.532
1.08 1.049 1.087 0.939
0.919 1.152 0.882 0.679
0.916 0.849 1.238 0.764
0.936 0.844 0.984 0.764
1.015 0.952 1.121 1.343",header=T,sep="\t")
mut3=read.table(text="weight length feeding_amount feeding_duration
0.959 1.023 0.980 1.291
1.139 0.937 0.724 0.769
1.212 1.022 1.041 1.026
1.187 0.910 0.852 0.618
1.059 1.040 0.851 1.530
1.074 0.947 1.017 0.888
0.861 0.961 0.924 0.937
0.926 0.909 0.892 1.900
0.960 1.034 0.925 1.349
0.881 0.831 1.008 1.168
0.900 1.081 0.952 0.929
1.029 1.023 0.765 0.925
1.120 0.844 0.871 1.202
0.913 0.903 1.091 1.186
0.866 0.982 1.193 1.409
1.007 1.011 0.747 1.674
1.061 1.036 1.097 1.001
1.065 1.064 1.210 0.784
1.054 0.931 1.085 1.372
1.001 0.850 1.303 1.027
0.995 0.960 0.783 1.304
0.884 1.024 0.933 0.930
0.915 0.854 1.019 1.114
0.846 1.143 0.923 1.983
0.993 1.029 0.736 1.192
0.942 1.098 1.099 1.200
1.090 1.027 1.109 1.208
1.009 0.958 1.264 1.727
1.088 0.860 1.261 0.686
0.915 1.046 0.834 0.944",header=T,sep="\t")
mut4=read.table(text="weight length feeding_amount feeding_duration
1.303 1.035 1.626 2.349
1.38 0.978 1.574 1.109
1.347 1.014 1.536 1.393
1.422 1.12 1.646 1.183
1.237 1.015 1.52 0.976
1.347 0.965 1.73 1.197
1.351 0.989 1.553 1.939
1.44 0.925 1.821 0.54
1.52 1.051 1.591 0.997
1.115 1.166 1.6 1.106
1.432 0.935 1.493 0.883
1.287 1.005 1.545 2.594
1.093 1.108 1.666 1.743
1.308 1.11 1.514 1.938
1.299 0.931 1.399 2.809
1.337 1.177 1.562 0.957
1.378 1.149 1.509 1.002
1.124 1.208 1.935 1.106
1.354 1.047 1.386 0.54
1.203 0.811 1.228 1.293
1.401 1.119 1.457 1.119
1.238 1.039 1.536 0.796
1.351 1.101 1.753 1.273
1.345 0.972 1.327 0.966
1.473 1.191 1.505 1.015
1.253 0.876 1.34 1.32
1.304 1.328 1.288 0.564
1.261 0.993 1.501 2.302
1.304 1.106 1.627 1.017
1.202 1.148 1.506 1.31",header=T,sep="\t")
wt=1
#weight
sample_IDs=c("mutant1","mutant2","mutant3","mutant4")
metric=rep("weight",4)
var="weight"
m1=mut1$weight
m2=mut2$weight
m3=mut3$weight
m4=mut4$weight
mutant1=data.frame(group=rep("mutant1",30),value=m1)
mutant2=data.frame(group=rep("mutant2",30),value=m2)
mutant3=data.frame(group=rep("mutant3",30),value=m3)
mutant4=data.frame(group=rep("mutant4",30),value=m4)
df=rbind(mutant1,mutant2,mutant3,mutant4)
hist(df[,2], main=var, xlab="Ratio")
dev.copy(png,"Q2_weight_hist.png")
## quartz_off_screen
## 3
dev.off()
## quartz_off_screen
## 2
ggplot(df,aes(x=group,y=value,color=group)) + scale_x_discrete(labels=c("1","2","3","4"))+labs(title=var,x=NULL,y="Ratio",color="BDNF mutant") + scale_color_manual(labels=c("1","2","3","4"),values=c("red","blue","orange","green")) +geom_point()+ geom_point(stat="summary",color="black",size=4, shape=18)+geom_errorbar(stat="summary",width=.5)+theme(text=element_text(size=24,face="bold"))
## No summary function supplied, defaulting to `mean_se()`
## No summary function supplied, defaulting to `mean_se()`
dev.copy(png,"Q2_weight.png")
## quartz_off_screen
## 3
dev.off()
## quartz_off_screen
## 2
shapiro.test(df$value)
##
## Shapiro-Wilk normality test
##
## data: df$value
## W = 0.98032, p-value = 0.0763
m1=tidy(shapiro.test(mutant1$value))
m2=tidy(shapiro.test(mutant2$value))
m3=tidy(shapiro.test(mutant3$value))
m4=tidy(shapiro.test(mutant4$value))
pval_shapiro.test=rbind(m1$p.value,m2$p.value,m3$p.value,m4$p.value)
m1=tidy(t.test(mutant1$value,mu=wt))
m2=tidy(t.test(mutant2$value,mu=wt))
m3=tidy(t.test(mutant3$value,mu=wt))
m4=tidy(t.test(mutant4$value,mu=wt))
pval_t.test=rbind(m1$p.value,m2$p.value,m3$p.value,m4$p.value)
m1=tidy(wilcox.test(mutant1$value,mu=wt,exact=FALSE))
m2=tidy(wilcox.test(mutant2$value,mu=wt, exact=FALSE))
m3=tidy(wilcox.test(mutant3$value,mu=wt, exact=FALSE))
m4=tidy(wilcox.test(mutant4$value,mu=wt, exact=FALSE))
pval_w.test=rbind(m1$p.value,m2$p.value,m3$p.value,m4$p.value)
number.tests=rep(4,4)
padj.t_bonferroni=pval_t.test*number.tests
padj.w_bonferroni=pval_w.test*number.tests
weight_pvalues=cbind(sample_IDs,metric,pval_shapiro.test,pval_t.test,padj.t_bonferroni,pval_w.test,padj.w_bonferroni)
colnames(weight_pvalues)=c("sample_ID","metric","pval_shapiro.test","pval_t.test","padj.t_bonferroni","pval_w.test","padj.w_bonferroni")
#length
sample_IDs=c("mutant1","mutant2","mutant3","mutant4")
metric=rep("length",4)
var="length"
m1=mut1$length
m2=mut2$length
m3=mut3$length
m4=mut4$length
mutant1=data.frame(group=rep("mutant1",30),value=m1)
mutant2=data.frame(group=rep("mutant2",30),value=m2)
mutant3=data.frame(group=rep("mutant3",30),value=m3)
mutant4=data.frame(group=rep("mutant4",30),value=m4)
df=rbind(mutant1,mutant2,mutant3,mutant4)
hist(df[,2], main=var, xlab="Ratio")
dev.copy(png,"Q2_length_hist.png")
## quartz_off_screen
## 3
dev.off()
## quartz_off_screen
## 2
ggplot(df,aes(x=group,y=value,color=group)) + scale_x_discrete(labels=c("1","2","3","4"))+labs(title=var,x=NULL,y="Ratio",color="BDNF mutant") + scale_color_manual(labels=c("1","2","3","4"),values=c("red","blue","orange","green")) +geom_point()+ geom_point(stat="summary",color="black",size=4, shape=18)+geom_errorbar(stat="summary",width=.5)+theme(text=element_text(size=24,face="bold"))
## No summary function supplied, defaulting to `mean_se()`
## No summary function supplied, defaulting to `mean_se()`
dev.copy(png,"Q2_length.png")
## quartz_off_screen
## 3
dev.off()
## quartz_off_screen
## 2
shapiro.test(df$value)
##
## Shapiro-Wilk normality test
##
## data: df$value
## W = 0.98514, p-value = 0.2116
m1=tidy(shapiro.test(mutant1$value))
m2=tidy(shapiro.test(mutant2$value))
m3=tidy(shapiro.test(mutant3$value))
m4=tidy(shapiro.test(mutant4$value))
pval_shapiro.test=rbind(m1$p.value,m2$p.value,m3$p.value,m4$p.value)
m1=tidy(t.test(mutant1$value,mu=wt))
m2=tidy(t.test(mutant2$value,mu=wt))
m3=tidy(t.test(mutant3$value,mu=wt))
m4=tidy(t.test(mutant4$value,mu=wt))
pval_t.test=rbind(m1$p.value,m2$p.value,m3$p.value,m4$p.value)
m1=tidy(wilcox.test(mutant1$value,mu=wt,exact=FALSE))
m2=tidy(wilcox.test(mutant2$value,mu=wt, exact=FALSE))
m3=tidy(wilcox.test(mutant3$value,mu=wt, exact=FALSE))
m4=tidy(wilcox.test(mutant4$value,mu=wt, exact=FALSE))
pval_w.test=rbind(m1$p.value,m2$p.value,m3$p.value,m4$p.value)
number.tests=rep(4,4)
padj.t_bonferroni=pval_t.test*number.tests
padj.w_bonferroni=pval_w.test*number.tests
length_pvalues=cbind(sample_IDs,metric,pval_shapiro.test,pval_t.test,padj.t_bonferroni,pval_w.test,padj.w_bonferroni)
colnames(length_pvalues)=c("sample_ID","metric","pval_shapiro.test","pval_t.test","padj.t_bonferroni","pval_w.test","padj.w_bonferroni")
#feeding_amount
sample_IDs=c("mutant1","mutant2","mutant3","mutant4")
metric=rep("feeding_amount",4)
var="feeding amount"
m1=mut1$feeding_amount
m2=mut2$feeding_amount
m3=mut3$feeding_amount
m4=mut4$feeding_amount
mutant1=data.frame(group=rep("mutant1",30),value=m1)
mutant2=data.frame(group=rep("mutant2",30),value=m2)
mutant3=data.frame(group=rep("mutant3",30),value=m3)
mutant4=data.frame(group=rep("mutant4",30),value=m4)
df=rbind(mutant1,mutant2,mutant3,mutant4)
hist(df[,2], main=var, xlab="Ratio")
dev.copy(png,"Q2_feeding_amount_hist.png")
## quartz_off_screen
## 3
dev.off()
## quartz_off_screen
## 2
ggplot(df,aes(x=group,y=value,color=group)) + scale_x_discrete(labels=c("1","2","3","4"))+labs(title=var,x=NULL,y="Ratio",color="BDNF mutant") + scale_color_manual(labels=c("1","2","3","4"),values=c("red","blue","orange","green")) +geom_point()+ geom_point(stat="summary",color="black",size=4, shape=18)+geom_errorbar(stat="summary",width=.5)+theme(text=element_text(size=24,face="bold"))
## No summary function supplied, defaulting to `mean_se()`
## No summary function supplied, defaulting to `mean_se()`
dev.copy(png,"Q2_feeding_amount.png")
## quartz_off_screen
## 3
dev.off()
## quartz_off_screen
## 2
shapiro.test(df$value)
##
## Shapiro-Wilk normality test
##
## data: df$value
## W = 0.98745, p-value = 0.3364
m1=tidy(shapiro.test(mutant1$value))
m2=tidy(shapiro.test(mutant2$value))
m3=tidy(shapiro.test(mutant3$value))
m4=tidy(shapiro.test(mutant4$value))
pval_shapiro.test=rbind(m1$p.value,m2$p.value,m3$p.value,m4$p.value)
m1=tidy(t.test(mutant1$value,mu=wt))
m2=tidy(t.test(mutant2$value,mu=wt))
m3=tidy(t.test(mutant3$value,mu=wt))
m4=tidy(t.test(mutant4$value,mu=wt))
pval_t.test=rbind(m1$p.value,m2$p.value,m3$p.value,m4$p.value)
m1=tidy(wilcox.test(mutant1$value,mu=wt,exact=FALSE))
m2=tidy(wilcox.test(mutant2$value,mu=wt, exact=FALSE))
m3=tidy(wilcox.test(mutant3$value,mu=wt, exact=FALSE))
m4=tidy(wilcox.test(mutant4$value,mu=wt, exact=FALSE))
pval_w.test=rbind(m1$p.value,m2$p.value,m3$p.value,m4$p.value)
number.tests=rep(4,4)
padj.t_bonferroni=pval_t.test*number.tests
padj.w_bonferroni=pval_w.test*number.tests
feeding_amount_pvalues=cbind(sample_IDs,metric,pval_shapiro.test,pval_t.test,padj.t_bonferroni,pval_w.test,padj.w_bonferroni)
colnames(feeding_amount_pvalues)=c("sample_ID","metric","pval_shapiro.test","pval_t.test","padj.t_bonferroni","pval_w.test","padj.w_bonferroni")
#feeding_duration
sample_IDs=c("mutant1","mutant2","mutant3","mutant4")
metric=rep("feeding_duration",4)
var="feeding duration"
m1=mut1$feeding_duration
m2=mut2$feeding_duration
m3=mut3$feeding_duration
m4=mut4$feeding_duration
mutant1=data.frame(group=rep("mutant1",30),value=m1)
mutant2=data.frame(group=rep("mutant2",30),value=m2)
mutant3=data.frame(group=rep("mutant3",30),value=m3)
mutant4=data.frame(group=rep("mutant4",30),value=m4)
df=rbind(mutant1,mutant2,mutant3,mutant4)
hist(df[,2], main=var, xlab="Ratio")
dev.copy(png,"Q2_feeding_duration_hist.png")
## quartz_off_screen
## 3
dev.off()
## quartz_off_screen
## 2
ggplot(df,aes(x=group,y=value,color=group)) + scale_x_discrete(labels=c("1","2","3","4"))+labs(title=var,x=NULL,y="Ratio",color="BDNF mutant") + scale_color_manual(labels=c("1","2","3","4"),values=c("red","blue","orange","green")) +geom_point()+ geom_point(stat="summary",color="black",size=4, shape=18)+geom_errorbar(stat="summary",width=.5)+theme(text=element_text(size=24,face="bold"))
## No summary function supplied, defaulting to `mean_se()`
## No summary function supplied, defaulting to `mean_se()`
dev.copy(png,"Q2_feeding_duration.png")
## quartz_off_screen
## 3
dev.off()
## quartz_off_screen
## 2
shapiro.test(df$value)
##
## Shapiro-Wilk normality test
##
## data: df$value
## W = 0.89064, p-value = 6.706e-08
m1=tidy(shapiro.test(mutant1$value))
m2=tidy(shapiro.test(mutant2$value))
m3=tidy(shapiro.test(mutant3$value))
m4=tidy(shapiro.test(mutant4$value))
pval_shapiro.test=rbind(m1$p.value,m2$p.value,m3$p.value,m4$p.value)
m1=tidy(t.test(mutant1$value,mu=wt))
m2=tidy(t.test(mutant2$value,mu=wt))
m3=tidy(t.test(mutant3$value,mu=wt))
m4=tidy(t.test(mutant4$value,mu=wt))
pval_t.test=rbind(m1$p.value,m2$p.value,m3$p.value,m4$p.value)
m1=tidy(wilcox.test(mutant1$value,mu=wt,exact=FALSE))
m2=tidy(wilcox.test(mutant2$value,mu=wt, exact=FALSE))
m3=tidy(wilcox.test(mutant3$value,mu=wt, exact=FALSE))
m4=tidy(wilcox.test(mutant4$value,mu=wt, exact=FALSE))
pval_w.test=rbind(m1$p.value,m2$p.value,m3$p.value,m4$p.value)
number.tests=rep(4,4)
padj.t_bonferroni=pval_t.test*number.tests
padj.w_bonferroni=pval_w.test*number.tests
feeding_duration_pvalues=cbind(sample_IDs,metric,pval_shapiro.test,pval_t.test,padj.t_bonferroni,pval_w.test,padj.w_bonferroni)
colnames(feeding_duration_pvalues)=c("sample_ID","metric","pval_shapiro.test","pval_t.test","padj.t_bonferroni","pval_w.test","padj.w_bonferroni")
#concatenate results
Q2.result=rbind(weight_pvalues,length_pvalues,feeding_amount_pvalues,feeding_duration_pvalues)
Q2.result.csv=write.csv(Q2.result, "~/Desktop/Q2.result.csv")