################################################### ### chunk number 1: setup ################################################### options(width=60)

################################################### ### chunk number 2: design-experiment ################################################### dataDir <- system.file("extdata", package="HGFSFN") library("limma") targets <- readTargets("ArrayDesc", path=dataDir) factors <- targets[["Cy3"]] exptlDesign <- data.frame(cellLine=factor(sub("PEC([0-9]+)_.*", "\\1", factors)), SFN=factor("Low", levels=c("Low", "High")), HGF=factor("Low", levels=c("Low", "High"))) exptlDesign[grep("SFN", targets[["Cy3"]]), "SFN"] <- "High" exptlDesign[grep("HGF", targets[["Cy3"]]), "HGF"] <- "High"

################################################### ### chunk number 3: load-data ################################################### library("HGFSFN") data("M")

################################################### ### chunk number 4: read-pedb ################################################### RG <- read.maimages(targets[["FileName"]], path=dataDir, source="genepix", wt.fun=wtflags())

################################################### ### chunk number 5: normalize ################################################### MA <- normalizeBetweenArrays(RG, "vsn", strata=rep(1:2, each=nrow(RG)/2))

################################################### ### chunk number 6: expreseion-set ################################################### library(convert) M <- as(MA, "ExpressionSet") sampleNames(M) <- targets[["Cy3"]] featureNames(M) <- MA[["genes"]][["ID"]]

################################################### ### chunk number 7: drop-rows ################################################### metaD <- read.csv(file.path(dataDir, "PEDB_ARRAY_annotations.csv"), as.is=TRUE) dropRows <- metaD[["Description"]] %in% c("EMPTY", "Failed Sequencing") M <- M[!dropRows,]

################################################### ### chunk number 8: model-matrix ################################################### X <- model.matrix(~SFN*HGF, exptlDesign)

################################################### ### chunk number 9: model-matrix-cell-line-only ################################################### XcellLine <- model.matrix(~cellLine, exptlDesign)

################################################### ### chunk number 10: first-half-select ################################################### Mh <- M[1:(nrow(M)/2),]

################################################### ### chunk number 11: lmFit-vers1 ################################################### fit <- lmFit(Mh, X) fitE <- eBayes(fit)

################################################### ### chunk number 12: decideTests-eg ################################################### summary(decideTests(fitE))

################################################### ### chunk number 13: data-frame-from-eg ################################################### Mh <- (exprs(M)[1:(nrow(M)/2),] + exprs(M)[1:(nrow(M)/2),]) / 2 observedVals <- as.vector(Mh) fittedVals <- as.vector(fitted(fitE)) residualVals <- as.vector(residuals(fitE, Mh)) df <- data.frame(sampleId=rep(sampleNames(M), each=nrow(Mh)), uniqueId=featureNames(M)[1:(nrow(M)/2)], SFN=rep(exptlDesign[["SFN"]], each=nrow(Mh)), HGF=rep(exptlDesign[["HGF"]], each=nrow(Mh)), observed=observedVals, fitteds=fittedVals, residuals=residualVals)

################################################### ### chunk number 14: residuals-by-sample-eg ################################################### library(lattice) print(plot(bwplot(sampleId~residuals, df)))

################################################### ### chunk number 15: full-design ################################################### Xfull <- model.matrix(~cellLine+SFN*HGF, exptlDesign)

################################################### ### chunk number 16: duplicate-correlations ################################################### dupCor <- duplicateCorrelation(exprs(M), design=Xfull, ndups=2, spacing=nrow(M)/2) dupFit <- lmFit(M, design=Xfull, ndups=2, spacing=nrow(M)/2, correlation=dupCor[["consensus.correlation"]]) dupFitE <- eBayes(dupFit)

################################################### ### chunk number 17: dup-decided-full ################################################### summary(dupDecided <- decideTests(dupFitE))