download code
## Chunk 1
library("limma")
library("CCl4")
dataPath = system.file("extdata", package="CCl4")
## Chunk 2
## dir(dataPath)
## Chunk 3
adf = read.AnnotatedDataFrame("samplesInfo.txt", 
    path=dataPath)
adf
## Chunk 4
targets = pData(adf)
targets$FileName = row.names(targets)
## Chunk 5
RG = read.maimages(targets, path=dataPath, source="genepix")
## Chunk 6
head(RG$genes)
## Chunk 7
par(mfrow=c(5,1))
imageplot(log2(RG$Rb[,1]), RG$printer, low="white", 
    high="red")
imageplot(log2(RG$Gb[,1]), RG$printer, low="white", 
    high="green")
imageplot(rank(RG$Rb[,1]), RG$printer, low="white", 
    high="red")
imageplot(rank(RG$Gb[,1]), RG$printer, low="white", 
    high="green")
imageplot(rank(log(RG$R[,1])+log(RG$G[,1])), 
    RG$printer, low="white", high="blue")
## Chunk 8
MA = normalizeWithinArrays(RG, method="none", 
    bc.method="none")
## Chunk 9
library("geneplotter")
smoothScatter(MA$A[, 1], MA$M[, 1], xlab="A", ylab="M")
abline(h=0, col="red")
## Chunk 10
plotformula = log2(RG$G)~col(RG$G) 
boxplot(plotformula, ylim=c(5,9), outline=FALSE, 
    col="forestgreen", xlab="arrays", 
    ylab=expression(log[2]~G), main="boxplot")
## Chunk 11
multidensity(plotformula, xlim=c(5,9), 
    main="densities", xlab=expression(log[2]~G))
## Chunk 12
rin = with(MA$targets, ifelse(Cy5=="CCl4", RIN.Cy5, 
    RIN.Cy3))
rin
select = (rin == max(rin))
RGgood = RG[, select]
adfgood = adf[select, ]
## Chunk 13
## library("vsn")
## ccl4 = justvsn(RGgood, backgroundsubtract=TRUE)
## Chunk 14
library("vsn")
## Chunk 15
ccl4 = justvsn(RGgood, backgroundsubtract=TRUE)
## Chunk 16
r = assayData(ccl4)$R
g = assayData(ccl4)$G
meanSdPlot(cbind(r, g))
## Chunk 17
rownames(pData(adfgood)) = sub("\\.gpr$", "", 
    rownames(pData(adfgood)))
pData(adfgood)
## Chunk 18
varMetadata(adfgood)$channel = factor(c("G", "R", "G", "R"),
    levels = c("G", "R", "_ALL_"))
## Chunk 19
phenoData(ccl4) = adfgood
validObject(ccl4)
## Chunk 20
ccl4AM = ccl4
assayData(ccl4AM) = assayDataNew(A=(r+g)/2, M=r-g)
## Chunk 21
varMetadata(phenoData(ccl4AM))$channel[] = "_ALL_"
validObject(ccl4AM)
## Chunk 22
## smoothScatter(assayData(ccl4AM)$A[,2], 
##               assayData(ccl4AM)$M[,2])
## abline(h=0, col="red")
## Chunk 23
design = modelMatrix(pData(ccl4AM), ref="DMSO")
## Chunk 24
fit = lmFit(assayData(ccl4AM)$M, design)
## Chunk 25
fit = eBayes(fit)
## Chunk 26
class(fit)
names(fit)
## Chunk 27
hist(fit$p.value, 1000)
## Chunk 28
fit$genes = pData(featureData(ccl4AM))
topTable(fit, number=10, adjust="BH")
## Chunk 29
plot(fit$coefficients, -log10(fit$p.value), pch=".")
## Chunk 30
## write.fit(fit, file="fit.tab")