######
#
#Load Limma package.  (Must first be installed.)
#
######

library(limma)

######
#
#Examine the limma users guide.
#
######

limmaUsersGuide()

######
#
#Set the working directory to the directory containing the raw
#data files.  The files are available at 
#http://www.public.iastate.edu/~dnett/microarray/sam3.0/
#You will need to download them to your hard drive
#and change the path accordingly.
#
######

setwd("C:\\z\\Courses\\S416\\sam3.0")

######
#
#Read the Targets.txt file that contains information about
#the slides and what was hybridized to each channel. 
#
######

targets=readTargets("Targets.txt")
targets

######
#
#Read the data files into R and examine content.
#
######

RG=read.maimages(targets,source="genepix")

attributes(RG)

head(RG$R)

######
#
#Examine boxplots of red and green backgrounds.
#
######

boxplot(data.frame(cbind(log2(RG$Gb),log2(RG$Rb))),
        main="Background",col=rep(c("green","red"),each=6))

######
#
#For slide 6, examine an image of the green background and green
#signal corresponding to the slide layout.
#
######

imageplot(log2(RG$Gb[,6]),RG$printer)
x11()
imageplot(log2(RG$G[,6]),RG$printer)

######
#
#Plot signal vs. background for the green channel of slide 6.
#
######

plot(log2(RG$Gb[,6]),log2(RG$G[,6]))
lines(c(-9,99),c(-9,99),col=2)

######
#
#Perform background correction for all slides.
#
######

RG=backgroundCorrect(RG,method="subtract")

attributes(RG)

head(RG$R)

######
#
#Examine side-by-side boxplots of background corrected data.
#
######

boxplot(log2(data.frame(RG$R,RG$G))[,as.vector(rbind(1:6,7:12))],
         col=rep(c("red","green"),6))

######
#
#Examine background corrected signal densities for each channel.
#
######

plotDensities(RG)

######
#
#Examine plots of the difference vs. the average
#log background corrected signal.
#
######

plotMA(RG)
plotMA(RG[,2])
plotMA(RG[,3])
plotMA(RG[,4])
plotMA(RG[,5])
plotMA(RG[,6])

######
#
#Loess normalize and median center the data.
#Examine the resulting Red-Green differences.
#
######

MA=normalizeWithinArrays(RG,method="loess")
MA=normalizeWithinArrays(MA,method="median")

attributes(MA)

plotMA(MA[,6])
plotDensities(MA)
boxplot(data.frame(MA$M))

######
#
#Scale normalize the differences.
#
######

MA=normalizeBetweenArrays(MA,method="scale")
boxplot(data.frame(MA$M))

######
#
#Analyze data.
#
######

#Create a design matrix for a gene-specific vector of 
#differences as the response.

targets

design=cbind(rep(1,6),rep(c(1,-1),each=3))
colnames(design)=c("Cy5minusCy3","PminusS")
design

#Use the limma function lmFit to fit the model.

fit=lmFit(MA, design)

#Set up contrasts and compute p-values using the limma approach.  

contrast.matrix=makeContrasts(PminusS, levels=design)
fit2=contrasts.fit(fit, contrast.matrix)
efit=eBayes(fit2)
attributes(efit)

hist(efit$p.value)