### R code from vignette source 'DiffExpressVignette.Snw'

###################################################
### code chunk number 1: DiffExpressVignette.Snw:84-85
###################################################
options(width=95)


###################################################
### code chunk number 2: getData
###################################################
library('ChimpHumanBrainData')
library(affy)
celfileDir = system.file('extdata',package='ChimpHumanBrainData')
celfileNames = list.celfiles(celfileDir)
brainBatch=ReadAffy(filenames=celfileNames,celfile.path=celfileDir,compress=TRUE)



###################################################
### code chunk number 3: createNames
###################################################
sampleNames(brainBatch)
sampleNames(brainBatch)=paste(rep(c("CH","HU"),each=12),rep(c(1:3,1:3),each=4),
  rep(c("Prefrontal","Caudate","Cerebellum","Broca"),6),sep="")


###################################################
### code chunk number 4: DiffExpressVignette.Snw:107-111
###################################################
brain.expr=exprs(brainBatch)
library(hexbin)
plot(hexplom(log2(brain.expr[,paste(rep(c("CH","HU"),each=3),
  c(1:3,1:3),rep("Prefrontal",6),sep="")])))


###################################################
### code chunk number 5: DiffExpressVignette.Snw:121-124
###################################################
trts=factor(paste(rep(c("CH","HU"),each=12),
  rep(c("Prefrontal","Caudate","Cerebellum","Broca"),6),sep=""))
blocks=factor(rep(1:6,each=4))


###################################################
### code chunk number 6: DiffExpressVignette.Snw:127-128
###################################################
brain.rma=rma(brainBatch)


###################################################
### code chunk number 7: DiffExpressVignette.Snw:156-158
###################################################
library(limma)
design.trt=model.matrix(~0+trts)


###################################################
### code chunk number 8: DiffExpressVignette.Snw:170-172
###################################################
library(statmod)
corfit <- duplicateCorrelation(brain.rma, design.trt, block = blocks)


###################################################
### code chunk number 9: DiffExpressVignette.Snw:178-179
###################################################
hist(tanh(corfit$atanh.correlations))


###################################################
### code chunk number 10: DiffExpressVignette.Snw:190-191
###################################################
fitTrtMean <- lmFit(brain.rma, design.trt, block = blocks, cor = corfit$consensus.correlation)


###################################################
### code chunk number 11: DiffExpressVignette.Snw:212-224
###################################################
colnames(design.trt)
contrast.matrix=makeContrasts(
  (trtsCHBroca+trtsCHCaudate+trtsCHCerebellum+trtsCHPrefrontal)/4
     -(trtsHUBroca+trtsHUCaudate+trtsHUCerebellum+trtsHUPrefrontal)/4,
  trtsCHBroca-trtsHUBroca,
  trtsCHCaudate-trtsHUCaudate,
  trtsCHCerebellum-trtsHUCerebellum,
  trtsCHPrefrontal-trtsHUPrefrontal,
  (trtsCHCerebellum-trtsHUCerebellum)-(trtsCHBroca-trtsHUBroca),
 levels=design.trt)
colnames(contrast.matrix)=
  c("ChVsHu","Broca","Caudate","Cerebellum","Prefrontal","Interact")


###################################################
### code chunk number 12: DiffExpressVignette.Snw:233-234
###################################################
fit.contrast=contrasts.fit(fitTrtMean,contrast.matrix)


###################################################
### code chunk number 13: DiffExpressVignette.Snw:242-243
###################################################
efit.contrast=eBayes(fit.contrast)


###################################################
### code chunk number 14: DiffExpressVignette.Snw:256-260
###################################################
par(mfrow=c(2,3))
for (i in 1:ncol(efit.contrast$p.value)) {
hist(efit.contrast$p.value[,i],main=colnames(efit.contrast$p.value)[i])
}


###################################################
### code chunk number 15: DiffExpressVignette.Snw:285-288
###################################################
genes=geneNames(brainBatch)
topTable(efit.contrast,coef=1,adjust.method="BY",n=10,p.value=1e-5,genelist=genes)
topTable(efit.contrast,coef=6,adjust.method="BY",n=10,p.value=1e-5,genelist=genes)


###################################################
### code chunk number 16: DiffExpressVignette.Snw:297-302
###################################################
write.table(file="fits.txt",
  cbind(genes,fitTrtMean$coefficients,efit.contrast$coefficients,efit.contrast$p.value),
  row.names=F,
  col.names=c("GeneID",colnames(fitTrtMean$coefficients),colnames(efit.contrast$p.value), 
  paste("p",colnames(efit.contrast$coefficients))),sep=",")


###################################################
### code chunk number 17: DiffExpressVignette.Snw:308-310
###################################################
library(qvalue)
q.values=apply(efit.contrast$p.value,2, function(x) qvalue(x)$qvalues)


###################################################
### code chunk number 18: sessionInfo
###################################################
toLatex(sessionInfo())


###################################################
### code chunk number 19: DiffExpressVignette.Snw:337-338
###################################################
print(gc())