################################################################
#
# This is the script to demonstrate the data analysis on a
# multiple factor experiment.
#
################################################################

# clear the work space
rm(list=ls())
# close all figures
graphics.off()

# load library and read in data
library(maanova)
nperm= 5; source('name.txt')
# read data, note that data has technical replicate 
paigen.raw <- read.madata("paigen.txt", designfile="paigendesign.txt", probeid=1, metarow=2, metacol=3, row=4, col=5, intensity=6, arrayType='twoColor', log.trans=T, spotflag=F,n.rep=2)

############################################
# STEP I: Data Quality Check on Raw Data
############################################
# grid check

# RI plot on raw data
riplot(paigen.raw)
# close figures
graphics.off()

# Array View on raw data
arrayview(paigen.raw)
# close figures
graphics.off()

paigen <- transform.madata(paigen.raw, method="rlowess")
graphics.off()

##############################
# STEP I: fixed model analysis
##############################
# Fit ANOVA model - full model; 'spot' can test technical replicate
anova.full.fix <- fitmaanova(madata = paigen, formula=~Array+Dye+Spot+Strain+Diet+Strain:Diet+Sample)
anova.noSample.fix <- fitmaanova(madata=paigen, formula=~Array+Dye+Spot+Strain+Diet+Strain:Diet)
# model without interaction
anova.noint.fix <- fitmaanova(madata=paigen, formula=~Array+Dye+Spot+Strain+Diet)
# model without strain effect
anova.nostrain.fix <- fitmaanova(madata=paigen, formula=~Array+Dye+Spot+Diet)
# model without diet effect
anova.nodiet.fix <- fitmaanova(madata=paigen, formul=~Array+Dye+Spot+Strain)

# permutation tests - all test use non-restricted residual shuffling
# test for Sample
test.Sample.fix <- matest(paigen, anova.full.fix, term="Sample", 
	n.perm=nperm, shuffle.method="resid",  test.method=rep(1,2))
#save(test.Sample.fix, file="testsamplefix.RData")
#load("testsamplefix.RData")

# test for interaction effect
test.int.fix <- matest(paigen, anova.noSample.fix, term="Strain:Diet", 
	n.perm=nperm, shuffle.method="resid",  test.method=rep(1,2))
#save(test.int.fix, file="testintfix.RData")
#load("testintfix.RData")
 
# test for strain effect
test.strain.fix <- matest(paigen, anova.noint.fix, term="Strain", n.perm=nperm,
        shuffle.method="resid", test.method=rep(1,2))
#save(test.strain.fix, file="teststrainfix.RData")
#load("teststrainfix.RData")

# test for diet effect
test.diet.fix <- matest(paigen, anova.noint.fix, term="Diet", n.perm=nperm,
        shuffle.method="resid", test.method=rep(1,2))
#save(test.diet.fix, file="testdietfix.RData")
#load("testdietfix.RData")

# volcano plot
idx.sample.fix <- volcano(test.Sample.fix, title="Test for Sample - fixed model")
x11();idx.int.fix <- volcano(test.int.fix, title="Interaction test - fixed model")
x11(); idx.strain.fix <- volcano(test.strain.fix, title="Strain test - fixed model")
x11(); idx.diet.fix <- volcano(test.diet.fix, title="Diet test - fixed model")

# T-test pairwise comparison of Strain
C <- matrix(c(1,-1,0,1,0,-1, 0,1,-1), nrow=3, byrow=T)
ttest.strain.fix <- matest(paigen, anova.noint.fix, term="Strain", Contrast=C, 
	n.perm=nperm, test.method=rep(1,2))
#save(ttest.strain.fix, file="tteststrainfix.RData")
#load("tteststrainfix.RData")
# volcano plot
volcano(ttest.strain.fix)

#################################
# STEP II: Mixed model analysis
#################################
# make severl model objects and fit ANOVA
# full model
anova.full.mix <- fitmaanova(madata=paigen, 
	formula=~Array+Spot+Dye+Strain+Diet+Strain:Diet+Sample,
	random=~Array+Spot+Sample)
varplot(anova.full.mix)

# test interaction effect
ftest.int.mix <- matest(data=paigen, anovaobj=anova.full.mix, term="Strain:Diet",n.perm=nperm,test.method=rep(1,2))

# volcano plot
x11(); idx.int.mix <- volcano(ftest.int.mix, title="Interaction test - mixed model")

# model without interaction
anova.noint.mix <- fitmaanova(madata=paigen, formula=~Dye+Array+Spot+Strain+Diet+Sample, 
        random=~Array+Spot+Sample)
varplot(anova.noint.mix)
ftest.strain.mix <- matest(data=paigen, anovaobj=anova.noint.mix, term="Strain", n.perm=nperm, test.method=rep(1,2))

# volcano plot
x11();idx.strain.mix <- volcano(ftest.strain.mix, title="Strain test - mixed model")