#########################################################
## THIS EXAMPLE R SCRIPT COMPUTES THE RESULTS OF
## THE PROPOSED DIFFERENTIAL EXPRESSION LIKELIHOOD 
## APPROACH FOR THE DATA OF HUNG ET AL. (2002)
## JOURNAL OF BIOLOGICAL CHEMISTRY
## 277:40309-40323.  IT COMPARES EXPRESSION IN TWO DIFFERENT
## E. COLI STRAINS LRP+ AND LRP- (FOUR AFFYMETRIX ARRAYS EACH).
## THE EXPERIMENT IS DESCRIBED CONCISELY
## IN THE LIMMA USER'S MANUAL, PP 54-57,
## DATED JUNE 6 2005,
## http://bioinf.wehi.edu.au/limma/usersguide.pdf
##
## WE ASSUME A BASIC FAMILIARITY WITH R/BIOCONDUCTOR,
## AND THAT ALL FILES AT
## http://www.bios.unc.edu/~fwright/Affyshrink
## ARE IN THE R WORKING DIRECTORY.
## 
## THE FIRST SET OF COMMANDS BELOW COMPUTES EXPRESSION
## ESTIMATES FOR THE ARRAYS USING RMA AND EXPRESSION
## DIFFERENCES USING LIMMA.  YOU CAN SKIP THIS PART BY INSTEAD
## TYPING THE COMMAND
## "x<-read.table("expression.log2scale",header=T)"
## AND GOING STRAIGHT TO THE SECTION "THE PROPOSED APPROACH"
##################################################

library(limma)
library(affy)
Data<-ReadAffy() # ASSUMES THE 8 CEL FILES ARE IN THE WORKING DIRECTORY
eset<-rma(Data) # ASSUMES ecolicdf PACKAGE HAS BEEN INSTALLED
pData(eset)
strain<-c("lrp-","lrp-","lrp-","lrp-","lrp+","lrp+","lrp+","lrp+")
design<-model.matrix(~factor(strain))
colnames(design)<-c("lrp-","lrp+vs-")
design
fit<-lmFit(eset,design)
fit<-eBayes(fit)
options(digits=5) # WE WILL SHOW MORE DIGITS THAN IN THE LIMMA MANUAL
topTable(fit, coef=2, n=40, adjust="fdr")

limma.t<-fit$t[,2] # THIS IS THE LIMMA MODERATED T-STATISTIC FOR DIFFERENTIAL
                   # EXPRESSION
x<-exprs(eset)   #####RMA expression intensity estimates####

###################################################
## THE PROPOSED APPROACH
###################################################
source("allfunctions.txt") # HERE YOU READ IN ALL THE FUNCTIONS FOR OUR METHOD
original.x<-2^x # THIS IS THE EXPRESSION MATRIX ON ORIGINAL
                # (NOT LOG2) SCALE
x1<-original.x[,5:8] # EXPRESSION MATRIX FOR GROUP 1
x2<-original.x[,1:4] # EXPRESSION MATRIX FOR GROUP 2

SIMPLE.REGRESSION(x1) # RESULT SHOWS XI^2 ESTIMATE POSITIVE
SIMPLE.REGRESSION(x2) # RESULT SHOWS XI^2 ESTIMATE POSITIVE

# RESULTS ABOVE IMPLY THAT WE SHOULD FIT THE LIKELIHOOD
# MODEL WITH THE ETA TERMS INCLUDED

output1eta<-DO.ALL.TRUEINTEGRAL(x1,iter=2)
output2eta<-DO.ALL.TRUEINTEGRAL(x2,iter=2)

meanvar.t<-(output1eta$muhat-output2eta$muhat)/sqrt(output1eta$sigma2.predicted/4+output2eta$sigma2.predicted/4)

showresults.table<-cbind(meanvar.t,limma.t)
showresults.table[order(-abs(meanvar.t)),][1:10,]  # NOTE THAT lrp IS RANKED
                                           # FIRST BY THE MEAN-VARIANCE METHOD