setwd("~/Desktop")
library(ggplot2)
library(broom)
#load data experiment 1
exp1_untreated=read.table(text="clone1 17.737 17.395 18.035
clone2 19.637 19.294 NA
clone3 13.469 13.035 13.401
clone4 10.979 10.893 11.023
clone5 27.588 27.192 28.025
clone6 10.388 9.293 10.284
clone7 7.280 7.193 7.011
clone8 9.338 9.039 8.923
clone9 14.932 14.385 14.923
clone10 16.956 16.242 NA",header=F,sep="\t")
exp1_treated=read.table(text="clone11 NA 7.349 8.125
clone12 19.932 19.984 19.395
clone13 17.865 17.395 18.022
clone14 21.075 21.233 21.087
clone15 24.668 24.563 24.035
clone16 20.349 NA 20.010
clone17 30.855 30.932 30.496
clone18 15.070 15.149 15.439
clone19 14.185 14.122 14.550
clone20 21.281 20.843 21.434",header=F,sep="\t")
exp1=rbind(exp1_untreated,exp1_treated)
colnames(exp1)=c("cloneID","r1","r2","r3")
exp1$cloneID=c("clone01","clone02","clone03","clone04","clone05","clone06","clone07","clone08","clone09","clone10","clone11","clone12","clone13","clone14","clone15","clone16","clone17","clone18","clone19","clone20")
exp1$mean=rowMeans(exp1[,2:4], na.rm=TRUE)
exp1$treatment=c(rep("untreated",10),rep("H2O2_treated",10))
#load data experiment 2
exp2_untreated=read.table(text="clone1 18.954 18.435 19.323
clone2 33.488 33.934 33.023
clone3 29.700 29.303 29.454
clone4 7.881 6.955 NA
clone5 18.149 19.143 18.943
clone6 11.844 12.320 12.433
clone7 8.476 8.343 8.749",header=F,sep="\t")
exp2_treated=read.table(text="clone1 22.267 22.758 22.039
clone2 35.028 NA 34.888
clone3 31.509 32.039 31.349
clone4 8.989 9.343 8.984
clone5 22.267 22.434 22.344
clone6 11.944 11.439 12.329
clone7 NA 10.485 10.221",header=F,sep="\t")
exp2=rbind(exp2_untreated,exp2_treated)
colnames(exp2)=c("cloneID","r1","r2","r3")
exp2$cloneID=c("clone01","clone02","clone03","clone04","clone05","clone06","clone07","clone01","clone02","clone03","clone04","clone05","clone06","clone07")
exp2$mean=rowMeans(exp2[,2:4], na.rm=TRUE)
exp2$treatment=c(rep("untreated",7),rep("H2O2_treated",7))
#experiment 1
treatment=exp1$treatment
mean1=exp1$mean
clone=exp1$cloneID
experiment1=data.frame(group=clone,treatment=treatment,value=mean1)
ggplot(experiment1, aes(x=treatment,y=value,fill=treatment))+geom_boxplot(outlier.shape = NULL,outlier.color = NA)+geom_jitter(color="black",shape=16,position=position_jitter(0.1),size=2)+labs(title="Experiment 1",x=NULL,y="DNA fragmentation",color="H2O2 treatment")+theme(text=element_text(size=18,face="bold"))
dev.copy(png,"Q3_experiment1.png")
## quartz_off_screen
## 3
dev.off()
## quartz_off_screen
## 2
#experiment 1 effect size
mean1_untreated=mean(exp1$mean[1:10])
mean1_treated=mean(exp1$mean[11:20])
sd1_untreated=sd(exp1$mean[1:10])
sd1_treated=sd(exp1$mean[11:20])
sd1_pooled=sqrt((sd1_untreated^2+sd1_treated^2)/2)
cohen1=(mean1_treated-mean1_untreated)/sd1_pooled #cohen1=0.74987
#experiment 1 power
E1_power=tidy(power.t.test(10,type="two.sample",sig.level = 0.05,power=NULL,sd=1,delta=cohen1)) #power=0.35480
#experiment 1 distribution
E1_untreated.shapiro=tidy(shapiro.test(experiment1$value[1:10]))
print(E1_untreated.shapiro)
## # A tibble: 1 x 3
## statistic p.value method
## <dbl> <dbl> <chr>
## 1 0.938 0.530 Shapiro-Wilk normality test
E1_treated.shapiro=tidy(shapiro.test(experiment1$value[11:20]))
print(E1_treated.shapiro)
## # A tibble: 1 x 3
## statistic p.value method
## <dbl> <dbl> <chr>
## 1 0.970 0.895 Shapiro-Wilk normality test
shapiro.test(experiment1$value)
##
## Shapiro-Wilk normality test
##
## data: experiment1$value
## W = 0.97133, p-value = 0.7828
hist(experiment1$value[1:10], main="experiment 1 untreated", xlab="DNA fragmentation")
dev.copy(png,"Q3_experiment1_untreated_hist.png")
## quartz_off_screen
## 3
dev.off()
## quartz_off_screen
## 2
hist(experiment1$value[11:20], main="experiment 1 treated", xlab="DNA fragmentation")
dev.copy(png,"Q3_experiment1_treated_hist.png")
## quartz_off_screen
## 3
dev.off()
## quartz_off_screen
## 2
#experiment 1 unpaired t test
E1_unpaired_t=tidy(t.test(experiment1$value[1:10],experiment1$value[11:20], paired=F))
#experiment 2
treatment=exp2$treatment
mean2=exp2$mean
clone=exp2$cloneID
experiment2=data.frame(group=clone,treatment=treatment,value=mean2)
ggplot(experiment2, aes(x=treatment,y=value,fill=treatment))+geom_boxplot(outlier.shape = NULL,outlier.color = NA)+geom_jitter(color="black",shape=16,position=position_jitter(0.1),size=2)+labs(title="Experiment 2",x=NULL,y="DNA fragmentation",color="H2O2 treatment")+theme(text=element_text(size=18,face="bold"))
dev.copy(png,"Q3_experiment2.png")
## quartz_off_screen
## 3
dev.off()
## quartz_off_screen
## 2
#experiment 2 effect size
mean2_untreated=mean(exp2$mean[1:7])
mean2_treated=mean(exp2$mean[8:14])
sd2_untreated=sd(exp2$mean[1:7])
sd2_treated=sd(exp2$mean[8:14])
sd2_pooled=sqrt((sd2_untreated^2+sd2_treated^2)/2)
cohen2=(mean2_treated-mean2_untreated)/sd2_pooled #cohen2=0.19426
#experiment 2 power
E2_power=tidy(power.t.test(7,type="paired",sig.level = 0.05,power=NULL,sd=1,delta=cohen2)) #power=0.06401229
#experiment 2 distribution
E2_untreated.shapiro=tidy(shapiro.test(experiment2$value[1:7]))
print(E2_untreated.shapiro)
## # A tibble: 1 x 3
## statistic p.value method
## <dbl> <dbl> <chr>
## 1 0.912 0.412 Shapiro-Wilk normality test
E2_treated.shapiro=tidy(shapiro.test(experiment2$value[8:14]))
print(E2_treated.shapiro)
## # A tibble: 1 x 3
## statistic p.value method
## <dbl> <dbl> <chr>
## 1 0.893 0.293 Shapiro-Wilk normality test
shapiro.test(experiment2$value)
##
## Shapiro-Wilk normality test
##
## data: experiment2$value
## W = 0.89963, p-value = 0.1114
hist(experiment2$value[1:7], main="experiment 2 untreated", xlab="DNA fragmentation")
dev.copy(png,"Q3_experiment2_untreated_hist.png")
## quartz_off_screen
## 3
dev.off()
## quartz_off_screen
## 2
hist(experiment2$value[8:14], main="experiment 2 treated", xlab="DNA fragmentation")
dev.copy(png,"Q3_experiment2_treated_hist.png")
## quartz_off_screen
## 3
dev.off()
## quartz_off_screen
## 2
#experiment 2 paired t test
E2_paired_t=tidy(t.test(experiment2$value[1:7],experiment2$value[8:14],paired=T))
#concatenate results
means=c(mean1_treated,mean1_untreated,mean2_treated,mean2_untreated)
cohen=c(cohen1,cohen1,cohen2,cohen2)
power=c(E1_power$power,E1_power$power,E2_power$power,E2_power$power)
pval_shapiro.test=rbind(E1_treated.shapiro$p.value,E1_untreated.shapiro$p.value,E2_untreated.shapiro$p.value,E2_untreated.shapiro$p.value)
pval_t.test=rbind(E1_unpaired_t$p.value,E1_unpaired_t$p.value,E2_paired_t$p.value,E2_paired_t$p.value)
experiment=c(rep("experiment1",2),rep("experiment2",2))
condition=c("H202_treated","untreated","H202_treated","untreated")
Q3.result=cbind(experiment,condition,means,cohen,power,pval_shapiro.test,pval_t.test)
colnames(Q3.result)=c("experiment","treatment","mean","cohen_d","power","pval_shapiro.test","pval_t.test")
Q3.result.csv=write.csv(Q3.result, "~/Desktop/Q3.result.csv")
print(Q3.result)
## experiment treatment mean cohen_d
## [1,] "experiment1" "H202_treated" "19.2453166666667" "0.749870957754811"
## [2,] "experiment1" "untreated" "14.66515" "0.749870957754811"
## [3,] "experiment2" "H202_treated" "20.379380952381" "0.194263781256146"
## [4,] "experiment2" "untreated" "18.3937142857143" "0.194263781256146"
## power pval_shapiro.test pval_t.test
## [1,] "0.354795853146098" "0.894579913119539" "0.110878321823104"
## [2,] "0.354795853146098" "0.529718595252342" "0.110878321823104"
## [3,] "0.0640122891892976" "0.412088064777498" "0.00713927463179746"
## [4,] "0.0640122891892976" "0.412088064777498" "0.00713927463179746"