Skip to content.

bioconductor.org

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

Sections

synthesis.R

################################################### ### chunk number 1: init ################################################### options(continue = " ")

################################################### ### chunk number 2: simulated ###################################################

k <- 7 # number of experiments n <- 10 # number of samples in each group

mu1 <- 0; mu2 <- 1; sigma <- 2.5

x <- matrix(rnorm(k n, mean = mu1, sd = sigma), n, k) y <- matrix(rnorm(k n, mean = mu2, sd = sigma), n, k)

################################################### ### chunk number 3: ttest ################################################### t.test(x[, 1], y[, 1])

################################################### ### chunk number 4: ttestpvals ################################################### pvals <- sapply(1:k, function(i) t.test(x[, i], y[, i])$p.value) pvals.less <- sapply(1:k, function(i) t.test(x[, i], y[, i], alternative = "less")$p.value) direction <- sapply(1:k, function(i) sign(mean(x[, i]) - mean(y[, i]))) pvals.less direction

################################################### ### chunk number 5: votes ################################################### sum(pvals.less < 0.05)

################################################### ### chunk number 6: pgamma ################################################### pgamma(sum(-log(pvals.less)), shape = k, rate = 1, lower.tail = FALSE)

################################################### ### chunk number 7: loaddata ################################################### library("nlme") library("GeneMeta") library("GeneMetaEx")

data("NevinsER") data("VantER") data("NevinsLN") data("HolstegeLN")

## test to make sure that we have matching data stopifnot(all(featureNames(VantER) == featureNames(NevinsER))) stopifnot(all(featureNames(NevinsLN) == featureNames(HolstegeLN)))

################################################### ### chunk number 8: zscores ###################################################

eSER <- list(NevinsER, VantER) eCER <- list(NevinsER$ERstatus, VantER$ERstatus) eCER <- lapply(eCER, function(x) ifelse(x == "pos", 1, 0))

wSFEM <- zScores(eSER, eCER, useREM=FALSE) wSREM <- zScores(eSER, eCER, useREM=TRUE)

################################################### ### chunk number 9: zscoresOutput ################################################### t(head(wSFEM, 4))

################################################### ### chunk number 10: plotT ###################################################

library("geneplotter") smoothScatter(wSFEM[,"Effect_Ex_1"], wSFEM[,"Effect_Ex_2"], xlab="Nevins", ylab="van't Veer", main = "per gene effect sizes in the two ER studies") abline(0, 1)

################################################### ### chunk number 11: qqplots ###################################################

library("lattice") plot(qqmath(wSFEM[, "Qvals"], type = c("p", "g"), distribution = function(p) qchisq(p, df = 1), panel = function(...) { panel.abline(0, 1) panel.qqmath(...) }))

################################################### ### chunk number 12: muestimates ################################################### plot(wSFEM[, "MUvals"], wSREM[, "MUvals"], pch = ".") abline(0, 1, col = "red")

################################################### ### chunk number 13: combining ###################################################

makeDf <- function(i, expr1, expr2, cov1, cov2, experiment.names = c("1", "2")) { y <- c(expr1[i, ], expr2[i, ]) # i-th gene treatment <- c(as.character(cov1), as.character(cov2)) experiment <- rep(experiment.names, c(length(cov1), length(cov2))) ans <- data.frame(y = y, treatment = factor(treatment), experiment = factor(experiment)) rownames(ans) <- NULL ans }

makeDf.ER <- function(i) { makeDf(i, exprs(NevinsER), exprs(VantER), NevinsER$ERstatus, VantER$ERstatus, experiment.names = c("Nevins", "Vant")) }

df3 <- makeDf.ER(3) summary(df3)

################################################### ### chunk number 14: combining2 ###################################################

makeDf2 <- function(i, ...) { args <- list(...) exprlist <- lapply(args, function(x) x$exprs[i, ]) covlist <- lapply(args, function(x) as.character(x$cov)) nsamples <- sapply(covlist, length) # number of samples experiment.names <- names(args) ans <- data.frame(y = unlist(exprlist), treatment = factor(unlist(covlist)), experiment = factor(rep(experiment.names, nsamples))) rownames(ans) <- NULL ans }

################################################### ### chunk number 15: scaling ################################################### scaled.exprs <- function(eset) { x <- exprs(eset) t(scale(t(x))) }

scaled.exprs.2 <- function(eset) { x <- exprs(eset) vx <- as.vector(x) x[] <- scale(vx, center = median(vx), scale = mad(vx)) x }

################################################### ### chunk number 16: scaleERdata ################################################### NevinsER.info <- list(exprs = scaled.exprs(NevinsER), cov = NevinsER$ERstatus) VantER.info <- list(exprs = scaled.exprs(VantER), cov = VantER$ERstatus) df3 <- makeDf2(3, Nevins = NevinsER.info, Vant = VantER.info)

################################################### ### chunk number 17: scaleLNdata eval=FALSE ################################################### ## NevinsLN.info <- list(exprs = scaled.exprs(NevinsLN), ## cov = NevinsLN$LNstatus) ## HolstegeLN.info <- list(exprs = scaled.exprs(HolstegeLN), ## cov = HolstegeLN$LNstatus)

################################################### ### chunk number 18: exampledf ################################################### str(df3)

################################################### ### chunk number 19: nlme eval=FALSE ################################################### ## ## ## ngenes <- nrow(exprs(NevinsER)) # takes very long ## ngenes <- 5 ## ## ERpvals <- numeric(ngenes) ## ## for(i in 1:ngenes) ## { ## cat(sprintf("\r%5g / %5g", i, ngenes)) ## dfi <- makeDf2(i, Nevins = NevinsER.info, Vant = VantER.info) ## fm.null <- ## lme(y ~ 1 + treatment, data = dfi, ## random = ~ 1 | experiment, ## method = "ML") ## fm.full <- ## lme(y ~ 1 + treatment, data = dfi, ## random = ~ 1 | experiment/treatment, ## method = "ML") ## ERpvals[i] <- anova(fm.null, fm.full)[["p-value"]][2] ## }

################################################### ### chunk number 20: ERpvs ###################################################

data("ERpvs") ERpvals <- unlist(eapply(ERpvs, function(x) x[2]))

plot(qqmath(~ ERpvals, outer = TRUE, main = "p-values for interaction effect", xlab = "Quantiles of Uniform", ylab = "Observed Quantiles", distribution = qunif, pch = "."))

################################################### ### chunk number 21: diffexp eval=FALSE ################################################### ## ## ## ngenes <- nrow(exprs(NevinsER)) # slow ## ngenes <- 20 ## ## ER.meandiff.combined <- data.frame(matrix(NA, nrow = ngenes, ncol = 4)) ## ER.meandiff.Nevins <- data.frame(matrix(NA, nrow = ngenes, ncol = 4)) ## ER.meandiff.Vant <- data.frame(matrix(NA, nrow = ngenes, ncol = 4)) ## ## colnames(ER.meandiff.combined) <- ## colnames(ER.meandiff.Nevins) <- ## colnames(ER.meandiff.Vant) <- ## c("Value", "Std.Error", "t.value", "p.value") ## ## for(i in 1:ngenes) ## { ## cat(sprintf("\r%5g / %5g", i, ngenes)) ## dfi <- makeDf2(i, Nevins = NevinsER.info, Vant = VantER.info) ## fm.lme <- ## lme(y ~ 1 + treatment, data = dfi, ## random = ~ 1 | experiment, ## method = "ML") ## fm.Nevins <- ## lm(y ~ 1 + treatment, data = dfi, ## subset = (experiment == "Nevins")) ## fm.Vant <- ## lm(y ~ 1 + treatment, data = dfi, ## subset = (experiment == "Vant")) ## ## ER.meandiff.combined[i, ] <- summary(fm.lme)$tTable[2, -3] ## ER.meandiff.Nevins[i, ] <- summary(fm.Nevins)$coefficients[2, ] ## ER.meandiff.Vant[i, ] <- summary(fm.Vant)$coefficients[2, ] ## } ##

################################################### ### chunk number 22: pvalqq ###################################################

load("ER.estimates.RData")

comb.df <- make.groups(Combined = ER.meandiff.combined$p.value, Nevins = ER.meandiff.Nevins$p.value, Vant = ER.meandiff.Vant$p.value)

plot(qqmath(~ data | which, comb.df, pch = ".", distribution = qunif, aspect = "iso"))

################################################### ### chunk number 23: coefpairs ################################################### cvERB1 <- ER.meandiff.combined$t.value cvERB2 <- ER.meandiff.Nevins$t.value cvERB3 <- ER.meandiff.Vant$t.value

pairs(data.frame(combined = cvERB1, Nevins = cvERB2, Vant = cvERB3), pch = ".")

################################################### ### chunk number 24: vennDiag ###################################################

pvERB1 <- ER.meandiff.combined$p.value pvERB2 <- ER.meandiff.Nevins$p.value pvERB3 <- ER.meandiff.Vant$p.value

vennDiag <- function(p1, p2, p3, c1, c2, c3, pthresh = 0.01, labels = TRUE) { p <- cbind(p1,p2,p3)0 e <- sign(cbind(c1,c2,c3)[sel, ]) p[sel, ] ordfun <- function(x) { order(rowSums(x matrix(2^(ncol(x):1), ncol = ncol(x), nrow = nrow(x), byrow = TRUE))) } image(y = 1:nrow(e), x = 1:ncol(e), z = t(e[ordfun(e), ]), ylab = paste(nrow(e), "probes"), xlab="", xaxt="n", col=c("skyblue","white","firebrick1"), main = "comparison of\nsignificant gene sets") axis(1, at = 1:ncol(e), labels = labels) }

vennDiag(pvERB1, pvERB3, pvERB2, cvERB1, cvERB3, cvERB2, labels=c("C", "V", "N"))

################################################### ### chunk number 25: eval=FALSE ################################################### ## openVignette("GeneMetaEx")

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.