################################################### ### chunk number 1: ################################################### #line 59 "vignettes/ecolitk/inst/doc/ecolitk.Rnw" library(ecolitk) ################################################### ### chunk number 2: ################################################### #line 71 "vignettes/ecolitk/inst/doc/ecolitk.Rnw" data(ecoligenomeSYMBOL2AFFY) data(ecoligenomeCHRLOC) ################################################### ### chunk number 3: ################################################### #line 82 "vignettes/ecolitk/inst/doc/ecolitk.Rnw" lac.i <- grep("^lac", ls(ecoligenomeSYMBOL2AFFY)) lac.symbol <- ls(ecoligenomeSYMBOL2AFFY)[lac.i] lac.affy <- unlist(lapply(lac.symbol, get, envir=ecoligenomeSYMBOL2AFFY)) beg.end <- lapply(lac.affy, get, envir=ecoligenomeCHRLOC) beg.end <- matrix(unlist(beg.end), nc=2, byrow=TRUE) ################################################### ### chunk number 4: ################################################### #line 93 "vignettes/ecolitk/inst/doc/ecolitk.Rnw" par(mfrow=c(2,2)) n <- 10 thetas <- rev(seq(0, 2 * pi, length=n)) rhos <- rev(seq(1, n) / n) xy <- polar2xy(rhos, thetas) colo <- heat.colors(n) plot(0, 0, xlim=c(-2, 2), ylim=c(-2, 2), type="n") for (i in 1:n) linesCircle(rhos[i]/2, xy$x[i], xy$y[i]) plot(0, 0, xlim=c(-2, 2), ylim=c(-2, 2), type="n") for (i in 1:n) polygonDisk(rhos[i]/2, xy$x[i], xy$y[i], col=colo[i]) plot(0, 0, xlim=c(-2, 2), ylim=c(-2, 2), type="n", xlab="", ylab="") for (i in 1:n) polygonArc(0, thetas[i], rhos[i]/2, rhos[i], center.x = xy$x[i], center.y = xy$y[i], col=colo[i]) plot(0, 0, xlim=c(-2, 2), ylim=c(-2, 2), type="n", xlab="", ylab="") for (i in (1:n)[-1]) { linesCircle(rhos[i-1], col="gray", lty=2) polygonArc(thetas[i-1], thetas[i], rhos[i-1], rhos[i], col=colo[i], edges=20) arrowsArc(thetas[i-1], thetas[i], rhos[i] + 1, col=colo[i], edges=20) } ################################################### ### chunk number 5: ################################################### #line 135 "vignettes/ecolitk/inst/doc/ecolitk.Rnw" cPlotCircle(main.inside = "E. coli - K12") ################################################### ### chunk number 6: ################################################### #line 147 "vignettes/ecolitk/inst/doc/ecolitk.Rnw" lac.o <- order(beg.end[, 1]) lac.i <- lac.i[lac.o] lac.symbol <- lac.symbol[lac.o] lac.affy <- lac.affy[lac.o] beg.end <- beg.end[lac.o, ] ################################################### ### chunk number 7: ################################################### #line 157 "vignettes/ecolitk/inst/doc/ecolitk.Rnw" cPlotCircle(main.inside = "E. coli - K12", main = "lac genes") polygonChrom(beg.end[, 1], beg.end[, 2], ecoli.len, 1, 1.4) ################################################### ### chunk number 8: ################################################### #line 163 "vignettes/ecolitk/inst/doc/ecolitk.Rnw" l <- data.frame(x=c(0.470, 0.48), y=c(0.87, 0.90)) cPlotCircle(xlim=range(l$x), ylim=range(l$y), main = "lac genes") polygonChrom(beg.end[, 1], beg.end[, 2], ecoli.len, 1, 1.007, col=rainbow(4)) legend(0.47, 0.9, legend=lac.symbol, fill=rainbow(4)) ################################################### ### chunk number 9: ################################################### #line 178 "vignettes/ecolitk/inst/doc/ecolitk.Rnw" library(Biobase) data(ecoligenomeBNUM2STRAND) data(ecoligenomeBNUM) data(ecoligenomeBNUM2MULTIFUN) data(ecoligenomeCHRLOC) affyids <- ls(ecoligenomeCHRLOC) affypos <- mget(affyids, ecoligenomeCHRLOC, ifnotfound=NA) ## 'unlist' as the mapping is one-to-one bnums <- unlist(mget(affyids, ecoligenomeBNUM, ifnotfound=NA)) strands <- unlist(mget(bnums, ecoligenomeBNUM2STRAND, ifnotfound=NA)) ## multifun <- mget(bnums, ecoligenomeBNUM2MULTIFUN, ifnotfound=NA) ## select the entries in the categorie "1.6" ## ("Macromolecules (cellular constituent) biosynthesis") f <- function(x) { if (all(is.na(x))) return(FALSE) length(grep("^1\\.6", x) > 0) } is.selected <- unlist(lapply(multifun, f)) cPlotCircle(main.inside="E.coli K12") beg.end <- matrix(unlist(affypos), nc=2, byrow=TRUE) ## plot 'bnums'... strand + good <- strands == ">" & is.selected linesChrom(beg.end[good, 1], beg.end[good, 2], ecoli.len, 1.4, col="red", lwd=3) ## plot 'bnums'... strand - good <- strands == "<" & is.selected linesChrom(beg.end[good, 1], beg.end[good, 2], ecoli.len, 1.5, col="blue", lwd=3) ################################################### ### chunk number 10: ################################################### #line 281 "vignettes/ecolitk/inst/doc/ecolitk.Rnw" library(Biobase) data(ecoligenome.operon) data(ecoligenomeSYMBOL2AFFY) tmp <- mget(ecoligenome.operon$gene.name, ecoligenomeSYMBOL2AFFY, ifnotfound=NA) ecoligenome.operon$affyid <- unname(unlist(tmp)) # clean up NAs ecoligenome.operon <- subset(ecoligenome.operon, !is.na(affyid)) ################################################### ### chunk number 11: ################################################### #line 299 "vignettes/ecolitk/inst/doc/ecolitk.Rnw" attach(ecoligenome.operon) affyoperons <- split(affyid, operon.name) detach(ecoligenome.operon) ## a sample of 5 operons affyoperons[18:22] ################################################### ### chunk number 12: ################################################### #line 309 "vignettes/ecolitk/inst/doc/ecolitk.Rnw" library(affy) ################################################### ### chunk number 13: ################################################### #line 314 "vignettes/ecolitk/inst/doc/ecolitk.Rnw" library(ecoliLeucine) data(ecoliLeucine) ################################################### ### chunk number 14: ################################################### #line 321 "vignettes/ecolitk/inst/doc/ecolitk.Rnw" abatch.nqt <- normalize(ecoliLeucine, method="quantiles") ################################################### ### chunk number 15: ################################################### #line 326 "vignettes/ecolitk/inst/doc/ecolitk.Rnw" ## the operon taken as an example: names(affyoperons)[18] #colnames(abatch.nqt@exprs) <- NULL eset <- computeExprSet(abatch.nqt, pmcorrect.method="pmonly", summary.method="medianpolish", ids = affyoperons[[18]]) ################################################### ### chunk number 16: computeExprSetOperon ################################################### #line 338 "vignettes/ecolitk/inst/doc/ecolitk.Rnw" operons.eset <- computeExprSet(abatch.nqt, pmcorrect.method="pmonly", summary.method="medianpolish", ids = unlist(affyoperons)) ################################################### ### chunk number 17: ################################################### #line 349 "vignettes/ecolitk/inst/doc/ecolitk.Rnw" library(multtest) my.ttest <- function(x, i, j) { pval <- t.test(x[i], x[j])$p.value return(pval) } is.lrpplus <- pData(operons.eset)$strain == "lrp+" is.lrpmoins <- pData(operons.eset)$strain == "lrp-" operons.ttest <- esApply(operons.eset, 1, my.ttest, is.lrpplus, is.lrpmoins) ## adjustment for multiple testing. operons.ttest.adj <- mt.rawp2adjp(operons.ttest, "BY")$adjp ## flag whether or not the probeset can be considered differentially expressed operons.diff.expr <- operons.ttest.adj < 0.01 ################################################### ### chunk number 18: ################################################### #line 371 "vignettes/ecolitk/inst/doc/ecolitk.Rnw" operons.i <- split(seq(along=operons.ttest), ecoligenome.operon$operon.name)