###1###
#
#Read a subset of the soybean data and examine first few rows.
#
#######

d=read.delim("http://www.public.iastate.edu/~dnett/microarray/soybeanloopdata.txt")
d[1:4,]

###2###
#
#Load some functions for normalization and analysis.
#
#######

source("http://www.public.iastate.edu/~dnett/microarray/rfunctions.txt")

###3###
#
#Lowess normalize the data within each slide and median center all channels.
#
#######

norm.d=msn(d[,-1])

###4###
#
#Examine normalized data.
#
#######

boxplot(norm.d)
slide.comp(norm.d,1)
slide.comp(norm.d,2)
slide.comp(norm.d,3)
slide.comp(norm.d,4)
slide.comp(norm.d,5)
slide.comp(norm.d,7)
slide.comp(norm.d,8)
slide.comp(norm.d,9)
slide.comp(norm.d,10)
slide.comp(norm.d,1)

#There is some evidence of problems, especially on slide 1.
#Could adjust frac to "straighten" data.

norm.d=msn(d[,-1],frac=.2)
slide.comp(norm.d,1)

#However, we will just go with the standard frac=.4

norm.d=msn(d[,-1])

###5###
#
#Prepare data for mixed model analysis in SAS.
#
#######

#Examine and read the design file.

design=read.delim("http://www.public.iastate.edu/~dnett/microarray/design.txt")
design

#Examine the function make.SAS.data

makeSASdata

#Set the directory in which you would like to save the SAS data.

setwd("C:/z/Courses/Smicroarray07/SoybeanLoop")

#Write the SAS data to a file in the directory specified above.

makeSASdata(cbind(geneID=d[,1],norm.d),design)

###6###
#
#Now use SAS to analyze the data (see analysis.sas).
#Read the results file created by SAS.
#
#######

results=read.delim("results.txt")
results[1:4,]

###7###
#
#Separate the dye and daf p-values and analyze.
#
#######

p=matrix(results[,6],byrow=T,ncol=2)
p.dye=p[,1]
p.daf=p[,2]

#Examine the histograms of p-values.

hist(p.daf)
hist(p.dye)

#Estimate the number of genes whose expression is
#constant accross the time points profiled.

estimate.m0(p.daf)

#Compute q-values for daf

q.daf=q.value(p.daf)

#Find number of genes with q-values <= 0.05.

sum(q.daf<=0.05)

#Obtain IDs for the most significant genes.

(d[order(p.daf),1])[1:10]

#Plot data for most significant genes.

plot.soy.gene(61)
plot.soy.gene(1559)
plot.soy.gene(1140)
plot.soy.gene(914)
plot.soy.gene(1882)
plot.soy.gene(1595)
plot.soy.gene(153)
plot.soy.gene(777)
plot.soy.gene(776)
plot.soy.gene(1523)

#Obtain IDs for genes with most significant dye effect.

(d[order(p.dye),1])[1:10]

#Plot data for genes with most significant dye effect.

plot.soy.gene(1567)
plot.soy.gene(859)
plot.soy.gene(1606)
plot.soy.gene(693)