################################################### ### chunk number 1: setup ################################################### library("RbcBook1") library("Biobase") library("multtest") library("limma") library("statmod") library("genefilter") library("annotate") library("GO") library("hgu95av2") set.seed(4711) quantilePlot <- function (px, py, quant=0.99, pch = ".", ...) { n <- length(px) stopifnot(length(px)==n) rpx <- rank(px, na.last = FALSE) dm <- 0.05 mids <- seq(dm, 1 - dm, by = dm) within <- function(x, x1, x2) { x >= x1 & x <= x2 } quantilewind <- function(mp) { quantile(py[within(rpx/n, mp - dm, mp + dm)], quant, na.rm = TRUE) } plot(rpx, py, pch = pch, col="#999999", ...) q <- sapply(mids, quantilewind) lines(mids * n, q, col = "black", type = "b", pch = 19) } ################################################### ### chunk number 2: dataALL ################################################### library("ALL") data(ALL) pdat <- pData(ALL) subset <- intersect(grep("^B", as.character(pdat$BT)), which(pdat$mol %in% c("BCR/ABL","NEG"))) eset <- ALL[,subset] ################################################### ### chunk number 3: prefiltering ################################################### library("genefilter") f1 <- pOverA(.25, log2(100)) f2 <- function(x)(IQR(x)>0.5) ff <- filterfun(f1, f2) selected <- genefilter(eset, ff) sum(selected) esetSub <- eset[selected,] ################################################### ### chunk number 4: multtestDo ################################################### cl <- as.numeric(esetSub$mol== "BCR/ABL") ## permutation test, using the Welch statistic resT <- cache("multtest-resT", mt.maxT(exprs(esetSub), classlabel=cl, B=10000)) ord <- order(resT$index) ## the original gene order rawp <- resT$rawp[ord] ## raw permutation $p$-values #names(resT) names(rawp) <- geneNames(esetSub) ################################################### ### chunk number 5: multtest eval=FALSE ################################################### ## cl <- as.numeric(esetSub$mol== "BCR/ABL") ## ## permutation test, using the Welch statistic ## resT <- mt.maxT(exprs(esetSub), classlabel=cl, B=10000)#1e+04) ## ord <- order(resT$index) ## the original gene order ## rawp <- resT$rawp[ord] ## raw permutation $p$-values #names(resT) ## names(rawp) <- geneNames(esetSub) ################################################### ### chunk number 6: histpraw ################################################### hist(rawp, breaks=30, main="", xlab="", col="#B2DF8A") ################################################### ### chunk number 7: WY ################################################### sum(resT$adjp<0.05) ################################################### ### chunk number 8: BH ################################################### res <- mt.rawp2adjp(rawp, proc = "BH") ################################################### ### chunk number 9: ################################################### res <- mt.rawp2adjp(rawp, proc = "BH") sum(res$adjp[,"BH"]<0.05) ################################################### ### chunk number 10: ################################################### rownames(res$adjp) <- names(rawp) ################################################### ### chunk number 11: calc.presel ################################################### IQRs <- esApply(eset, 1, IQR) intensityscore <- esApply(eset, 1, function(x) quantile(x, 0.75)) abs.t <- abs(mt.teststat(exprs(eset), classlabel=cl)) ## Welch statistic ################################################### ### chunk number 12: plotpresel ################################################### par(mfrow=c(1,2)) quantilePlot(IQRs, abs.t, quant=0.95, xlab="rank(IQR)", ylab="abs.t") quantilePlot(intensityscore, abs.t, quant=0.95, xlab="rank(intensity score)", ylab="abs.t") par(mfrow=c(1,1)) ################################################### ### chunk number 13: GOp1 ################################################### tykin <- unique(lookUp("GO:0004713", "hgu95av2", "GO2ALLPROBES")) length(tykin) ################################################### ### chunk number 14: GO2 ################################################### gN <- geneNames(esetSub) whSel <- match(tykin, gN) whSel <- whSel[!is.na(whSel)] ################################################### ### chunk number 15: GO3 ################################################### resTGO <- cache("AvHDS-resTGO", mt.maxT(exprs(esetSub[whSel,]), classlabel=cl, B=10000)) ordGO <- order(resTGO$index) #$ the original gene order adjpGO <- resTGO$adjp[ordGO] #$ raw permutation $p$-values #names(resT) names(adjpGO) <- geneNames(esetSub[whSel,]) otherp <- resT$adjp[ord] #$ names(otherp) <- geneNames(esetSub) op2 <- otherp[names(adjpGO)] ndiff <- sum(adjpGO < 0.05) ################################################### ### chunk number 16: GO4 ################################################### print("GO analysis") sort(adjpGO)[1:ndiff] print("All Genes") sort(op2)[1:ndiff] ################################################### ### chunk number 17: preprocessKidney ################################################### library("vsn") library("kidpack") data(qua) data(spotanno) data(hybanno) data(cloneanno) data(esetSpot) ################################################### ### chunk number 18: kidney.data ################################################### pdat <- pData(esetSpot) table(pdat$type) ################################################### ### chunk number 19: model ################################################### design <- model.matrix(~ -1 + factor(pdat$type)) colnames(design) <- c("ccRCC", "chRCC", "pRCC") ################################################### ### chunk number 20: dupcorDo ################################################### dupcor <- cache("AvHDS-dupcor", duplicateCorrelation(exprs(esetSpot),design=design,ndups=2,spacing=4224)) fit <- cache("AvHDS-fit", lmFit(esetSpot, design=design, ndups=2, spacing=4224, correlation=dupcor$cor)) ################################################### ### chunk number 21: dupcorShow1 eval=FALSE ################################################### ## dupcor <- duplicateCorrelation(exprs(esetSpot),design=design,ndups=2,spacing=4224) ## fit <- lmFit(esetSpot, design=design, ndups=2, spacing=4224, correlation=dupcor$cor) ################################################### ### chunk number 22: dupcorShow2 ################################################### dupcor$cor ################################################### ### chunk number 23: fit.contrasts ################################################### contrast.matrix <- makeContrasts(ccRCC-chRCC, ccRCC-pRCC, chRCC-pRCC, levels=design) contrast.matrix fit2 <- contrasts.fit(fit, contrast.matrix) ################################################### ### chunk number 24: tmod ################################################### fit3 <- eBayes(fit2) ################################################### ### chunk number 25: topTable ################################################### topTable(fit3, coef=3, n=8, adjust.method="fdr") ################################################### ### chunk number 26: decideTests ################################################### clas <- decideTests(fit3, method="nestedF", adjust.method="fdr", p=0.05) ##get the numbers of significant genes per contrast: colSums(abs(clas)) ################################################### ### chunk number 27: compareDupCor ################################################### nclones <- 4224 datAvDup <- (exprs(esetSpot)[1:nclones,] + exprs(esetSpot)[nclones+1:nclones,])/2 fitAvDup <- lmFit(datAvDup, design=design) fit2AvDup <- contrasts.fit(fitAvDup, contrast.matrix) fit3AvDup <- eBayes(fit2AvDup) ################################################### ### chunk number 28: AvDup ################################################### plot(log10(fit3AvDup$p.value[,3]), log10(fit3$p.value[,3]), xlab="log10(p_avdup)", ylab="log10(p_dupcor)", xlim=c(-42,0), ylim=c(-42,0), pch=19, cex=0.3) abline(c(0,1),col="#E31A1C") ################################################### ### chunk number 29: xxx ################################################### gc() ################################################### ### chunk number 30: loadandnormalizeestrogendata ################################################### library("estrogen") library("limma") library("hgu95av2cdf") datadir <- system.file("extdata", package="estrogen") targets <- readTargets("phenoData.txt",path=datadir,sep="") covdesc <- list("present or absent","10 or 48 hours") names(covdesc) <- names(targets)[-1] pdata <- new("phenoData",pData=targets[,-1],varLabels=covdesc) rownames(pData(pdata)) <- targets[,1] gc() esAB <- ReadAffy(filenames=file.path(datadir,targets$filename),phenoData=pdata) esEset <- rma(esAB) ################################################### ### chunk number 31: ################################################### esEset pData(esEset) ################################################### ### chunk number 32: outlierplots ################################################### par(mfrow=c(1,2)) par(las=2) for (i in c("728_at","33379_at")){ expvals <- exprs(esEset)[i,] plot(expvals,axes=F,cex=1.5, xlab="Conditions",ylab="log 2 Expression Estimate") if(i=="728_at") points(7:8,expvals[7:8],pch=16,cex=1.5) if(i=="33379_at") points(1:2,expvals[1:2],pch=16,cex=1.5) axis(1,at=1:8,labels=c("et1","et2","Et1","Et2","eT1","eT2","ET1","ET2")) axis(2) title(i) } par(mfrow=c(1,1)) ################################################### ### chunk number 33: outliers ################################################### library("factDesign") op1 <- outlierPair(exprs(esEset)["728_at",],INDEX=pData(esEset)) op1 madOutPair(exprs(esEset)["728_at",],op1[[3]]) op2 <- outlierPair(exprs(esEset)["33379_at",],INDEX=pData(esEset)) madOutPair(exprs(esEset)["33379_at",],op2[[3]]) ################################################### ### chunk number 34: model ################################################### pdat <- pData(esEset) design <- model.matrix(~factor(estrogen)*factor(time.h),pdat) colnames(design) <- c("Intercept","ES","T48","ES:T48") fit <- lmFit(esEset,design) fit$coefficients[1:3,] ################################################### ### chunk number 35: ################################################### contM <- cbind(es10=c(0,1,0,0),es48=c(0,1,0,1)) fitC <- contrasts.fit(fit,contM) fitC <- eBayes(fitC) esClas <- classifyTestsF(fitC, p=0.00001) print(colSums(abs(esClas))) ################################################### ### chunk number 36: heatmap1 ################################################### both1048 <- which(esClas[,"es10"] != 0 & esClas[,"es48"] != 0) heatmap(exprs(esEset)[both1048,],Colv=NA,col=cm.colors(256),labCol=targets[,"filename"]) ################################################### ### chunk number 37: heatmap2 ################################################### library("stats") only10 <- which(esClas[,"es10"] != 0 & esClas[,"es48"] == 0) heatmap(exprs(esEset)[only10,],Colv=NA,col=cm.colors(256),labCol=targets[,"filename"],cexRow=1,cexCol=1) title("Estrogen target at 10 hours only") ################################################### ### chunk number 38: heatmap3 ################################################### only48 <- which(esClas[,"es10"] == 0 & esClas[,"es48"] != 0) heatmap(exprs(esEset)[only48,],Colv=NA,col=cm.colors(256),labCol=targets[,"filename"]) title("Estrogen target at 48 hours only")