library(lattice)
library(nlme)
library(MASS)

#Comet Assays
mydata1<-read.table("cometgill.csv",header=T, sep=",")

mydata2<-read.table("comethaemo.csv",header=T, sep=",")

lm(tail ~ (treatment:group) + (run:group) + (treatment:group)*(run:group), data = mydata1)
anova(lm(tail ~ (treatment:group) + (run:group) + (treatment:group)*(run:group), data = mydata1))

lm(tail ~ (treatment:group) + (run:group) + (treatment:group)*(run:group), data = mydata2)
anova(lm(tail ~ (treatment:group) + (run:group) + (treatment:group)*(run:group), data = mydata2))

#Oxidative Stress
mydatasod<-read.table("sod.csv",header=T, sep=",")

mydatatbars<-read.table("tbar.csv",header=T, sep=",")

mydataGPx<-read.table("GPx.txt",header=T,sep="\t")

mydataGSH<-read.table("GSH.txt",header=T,sep="\t")

lm(sod ~ treatment + run + treatment*run, data = mydatasod)
anova(lm(sod ~ treatment + run + treatment*run, data = mydatasod))

lm(tbars ~ treatment + run + treatment*run, data = mydatatbars)
anova(lm(tbars ~ treatment + run + treatment*run, data = mydatatbars))

t.test(GSH~Treatment, data = mydataGSH)

t.test(GPx~Treatment, data = mydataGPx)

#Oxygen Consumption
resp<-read.csv("respiration2.csv",header=TRUE)

head(resp)
str(resp)

resp$fTank<-as.factor(resp$Tank)

xyplot(O2~Time|Treatment, groups=fTank, type = c("p","r"), data=resp, auto.key=FALSE, xlab="Time", ylab="oxygen saturation (%)")

model1<-lme(O2~Treatment*Time, random=~1 | fTank, data=resp)

plot(model1, which=1)

summary(model1)
anova(model1)

#Algal Clearance
cell_counts<-read.csv("raw_cell_counts.csv",header=TRUE)

head(cell_counts)
str(cell_counts)

cell_counts$fTank<-as.factor(cell_counts$Tank)

xyplot(mean_cell_count~Time|Treatment, groups=fTank, type = c("p","r"), data=cell_counts, auto.key=FALSE, xlab="Time", ylab="cell count")

model4<-lme(log(mean_cell_count+1)~Treatment*Time, random=~1|fTank, weights=varIdent(form=~1|fTank), data=cell_counts)

plot(model4, which=1)

anova(model4)

#Valve Movement
mydatamovement<-read.table("valve movement.txt",header=T,sep="\t")

wilcox.test(gape~treatment, data = mydatamovement)

wilcox.test(time~treatment, data = mydatamovement)