Skip to content.

bioconductor.org

Bioconductor is an open source and open development software project
for the analysis and comprehension of genomic data.

Sections

ComplexAnalysisLecture.R

################################################### ### chunk number 1: setup ################################################### options(width=60) dataDir <- system.file("extdata", package="arrayQ") workDir <- tempdir()

################################################### ### chunk number 2: expt-design-ours-setup ################################################### library("limma") description <- readTargets("ArrayDesc", path=dataDir) factors <- description[["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", description[["Cy3"]]), "SFN"] <- "High" exptlDesign[grep("HGF", description[["Cy3"]]), "HGF"] <- "High"

################################################### ### chunk number 3: expt-design-ours ################################################### head(exptlDesign) table(exptlDesign[,-1])

################################################### ### chunk number 4: formula-a eval=FALSE ################################################### ## Response ~ Dependent

################################################### ### chunk number 5: formula-our eval=FALSE ################################################### ## ~1 + SFN + HGF + SFN:HGF ## ~SFN*HGF

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

################################################### ### chunk number 7: model-fit-setup ################################################### RG <- read.maimages(description$FileName, path=dataDir, source="genepix", wt.fun=wtflags())

MA <- normalizeBetweenArrays(RG, "vsn", strata=rep(1:2, each=nrow(RG)/2))

library(convert) M <- as(MA, "ExpressionSet") sampleNames(M) <- description[["Cy3"]] featureNames(M) <- MA[["genes"]][["ID"]]

metaD <- read.csv(file.path(dataDir, "PEDB_ARRAY_annotations.csv"), as.is=TRUE) dropRows <- metaD[["Description"]] %in% c("EMPTY", "Failed sequencing") M <- M[!dropRows,]

dupCor <- duplicateCorrelation(exprs(M), ndup=2, spacing=nrow(M)/2)

################################################### ### chunk number 8: M-eg ################################################### M

################################################### ### chunk number 9: M-half-eg ################################################### Mh <- M[1:(nrow(M)/2),]

################################################### ### chunk number 10: limma-lmFit-eg ################################################### fit <- lmFit(Mh, X) efit <- eBayes(fit)

################################################### ### chunk number 11: limma-lmFit-eg-summary ################################################### names(efit) head(efit[["coefficients"]])

################################################### ### chunk number 12: expected-value ################################################### head(cbind(X, X %*% efit[["coefficients"]][1,]))

################################################### ### chunk number 13: decideTests-eg ################################################### decided <- decideTests(efit) summary(decided)

################################################### ### chunk number 14: topTable-eg ################################################### topTable(efit, coef=2, n=sum(decided[,2]!=0))

################################################### ### chunk number 15: duplicate-correlations-eg ################################################### dupCor <- duplicateCorrelation(exprs(M), design=X, ndups=2, spacing=nrow(M)/2)

################################################### ### chunk number 16: ################################################### summary(dupCor[["atanh.correlations"]]) dupCor[["consensus"]]

################################################### ### chunk number 17: duplicate-analysis ################################################### dupFit <- lmFit(M, design=X, ndups=2, spacing=nrow(M)/2, correlation=dupCor[["consensus.correlation"]]) dupFitE <- eBayes(dupFit) dupDecided <- decideTests(dupFitE)

################################################### ### chunk number 18: ################################################### summary(dupDecided)

################################################### ### chunk number 19: cell-line-eg ################################################### levels(exptlDesign[["cellLine"]])

################################################### ### chunk number 20: cell-line-design ################################################### head(model.matrix(~cellLine, exptlDesign))

################################################### ### chunk number 21: full-model-eg ################################################### X <- model.matrix(~cellLine+SFN*HGF, exptlDesign) dupCor <- duplicateCorrelation(exprs(M), design=X, ndups=2, spacing=nrow(M)/2) dupFit <- lmFit(M, design=X, ndups=2, spacing=nrow(M)/2, correlation=dupCor[["consensus.correlation"]]) dupFitE <- eBayes(dupFit) dupDecided <- decideTests(dupFitE)

################################################### ### chunk number 22: full-model-highlights ################################################### dupCor[["consensus.correlation"]] summary(dupDecided)

################################################### ### chunk number 23: full-model-topTable ################################################### topTable(dupFitE, coef=6)

################################################### ### chunk number 24: fitted ################################################### fittedVals <- fitted(dupFitE) dim(fittedVals)

################################################### ### chunk number 25: ave-exprs ################################################### observedVals <- (exprs(M)[1:(nrow(M)/2),]+exprs(M)[nrow(M)/2+1:(nrow(M)/2),])/2 residualVals <- residuals(dupFitE, observedVals)

################################################### ### chunk number 26: as-data-frame ################################################### df <- data.frame(sampleId=rep(1:20, each=nrow(observedVals)), SFN=rep(exptlDesign[["SFN"]], each=nrow(observedVals)), HGF=rep(exptlDesign[["HGF"]], each=nrow(observedVals)), observed=as.vector(observedVals), fitted=as.vector(fittedVals), residuals=as.vector(residualVals))

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

News
2009-10-26

BioC 2.5, consisting of 352 packages and designed to work with R 2.10.z, was released today.

2009-01-07

R, the open source platform used by Bioconductor, featured in a series of articles in the New York Times.