### R code from vignette source 'EnrichmentBrowser.Rnw' ### Encoding: UTF-8 ################################################### ### code chunk number 1: style ################################################### BiocStyle::latex() ################################################### ### code chunk number 2: read-eset ################################################### library(EnrichmentBrowser) data.dir <- system.file("extdata", package="EnrichmentBrowser") exprs.file <- file.path(data.dir, "exprs.tab") pdat.file <- file.path(data.dir, "pData.tab") fdat.file <- file.path(data.dir, "fData.tab") eset <- read.eset(exprs.file, pdat.file, fdat.file) ################################################### ### code chunk number 3: help (eval = FALSE) ################################################### ## ?read.eset ## ?ExpressionSet ################################################### ### code chunk number 4: load-ALL ################################################### library(ALL) data(ALL) ################################################### ### code chunk number 5: subset-ALL ################################################### ind.bs <- grep("^B", ALL$BT) ind.mut <- which(ALL$mol.biol %in% c("BCR/ABL", "NEG")) sset <- intersect(ind.bs, ind.mut) all.eset <- ALL[, sset] ################################################### ### code chunk number 6: show-ALL ################################################### dim(all.eset) exprs(all.eset)[1:4,1:4] ################################################### ### code chunk number 7: probe2gene ################################################### all.eset <- probe.2.gene.eset(all.eset) head(featureNames(all.eset)) ################################################### ### code chunk number 8: show-probe2gene ################################################### head(fData(eset)) ################################################### ### code chunk number 9: load-airway ################################################### library(airway) data(airway) ################################################### ### code chunk number 10: preproc-airway ################################################### expr <- assays(airway)[[1]] expr <- expr[grep("^ENSG", rownames(expr)),] expr <- expr[rowMeans(expr) > 10,] air.eset <- new("ExpressionSet", exprs=expr, annotation="hsa") dim(air.eset) exprs(air.eset)[1:4,1:4] ################################################### ### code chunk number 11: norm-ma ################################################### before.norm <- exprs(all.eset) all.eset <- normalize(all.eset, norm.method="quantile") after.norm <- exprs(all.eset) ################################################### ### code chunk number 12: plot-norm ################################################### par(mfrow=c(1,2)) boxplot(before.norm) boxplot(after.norm) ################################################### ### code chunk number 13: norm-rseq ################################################### norm.air <- normalize(air.eset, norm.method="quantile") ################################################### ### code chunk number 14: lgc (eval = FALSE) ################################################### ## ids <- head(featureNames(air.eset)) ## lgc <- EDASeq::getGeneLengthAndGCContent(ids, org="hsa", mode="biomart") ################################################### ### code chunk number 15: norm2-rseq ################################################### lgc.file <- file.path(data.dir, "air_lgc.tab") fData(air.eset) <- read.delim(lgc.file) norm.air <- normalize(air.eset, within=TRUE) ################################################### ### code chunk number 16: sample-groups-ALL ################################################### pData(all.eset)$GROUP <- ifelse(all.eset$mol.biol == "BCR/ABL", 1, 0) table(pData(all.eset)$GROUP) ################################################### ### code chunk number 17: sample-groups-airway ################################################### pData(air.eset)$GROUP <- ifelse(colData(airway)$dex == "trt", 1, 0) table(pData(air.eset)$GROUP) ################################################### ### code chunk number 18: sample-blocks ################################################### pData(air.eset)$BLOCK <- colData(airway)$cell table(pData(air.eset)$BLOCK) ################################################### ### code chunk number 19: DE-ana-ALL ################################################### all.eset <- de.ana(all.eset) head(fData(all.eset), n=4) ################################################### ### code chunk number 20: plot-DE ################################################### par(mfrow=c(1,2)) pdistr(fData(all.eset)$ADJ.PVAL) volcano(fData(all.eset)$FC, fData(all.eset)$ADJ.PVAL) ################################################### ### code chunk number 21: DE-exmpl ################################################### fData(all.eset)[ which.min(fData(all.eset)$ADJ.PVAL), ] ################################################### ### code chunk number 22: DE-ana-airway ################################################### air.eset <- de.ana(air.eset, de.method="edgeR") head(fData(air.eset), n=4) ################################################### ### code chunk number 23: idmap-idtypes ################################################### id.types("hsa") ################################################### ### code chunk number 24: idmap-airway ################################################### head(featureNames(air.eset)) air.eset <- map.ids(air.eset, from="ENSEMBL", to="ENTREZID") head(featureNames(air.eset)) ################################################### ### code chunk number 25: get-kegg-gs (eval = FALSE) ################################################### ## kegg.gs <- get.kegg.genesets("hsa") ################################################### ### code chunk number 26: get-go-gs (eval = FALSE) ################################################### ## go.gs <- get.go.genesets(org="hsa", onto="BP", mode="GO.db") ################################################### ### code chunk number 27: parseGMT ################################################### gmt.file <- file.path(data.dir, "hsa_kegg_gs.gmt") hsa.gs <- parse.genesets.from.GMT(gmt.file) length(hsa.gs) hsa.gs[1:2] ################################################### ### code chunk number 28: sbea-methods ################################################### sbea.methods() ################################################### ### code chunk number 29: sbea ################################################### sbea.res <- sbea(method="ora", eset=all.eset, gs=hsa.gs, perm=0, alpha=0.05) gs.ranking(sbea.res) ################################################### ### code chunk number 30: ea-browse (eval = FALSE) ################################################### ## ea.browse(sbea.res) ################################################### ### code chunk number 31: dummy-sbea ################################################### dummy.sbea <- function(eset, gs, alpha, perm) { sig.ps <- sample(seq(0,0.05, length=1000),5) insig.ps <- sample(seq(0.1,1, length=1000), length(gs)-5) ps <- sample(c(sig.ps, insig.ps), length(gs)) names(ps) <- names(gs) return(ps) } ################################################### ### code chunk number 32: sbea2 ################################################### sbea.res2 <- sbea(method=dummy.sbea, eset=all.eset, gs=hsa.gs) gs.ranking(sbea.res2) ################################################### ### code chunk number 33: dwnld-pwys (eval = FALSE) ################################################### ## pwys <- download.kegg.pathways("hsa") ################################################### ### code chunk number 34: compile-grn ################################################### pwys <- file.path(data.dir, "hsa_kegg_pwys.zip") hsa.grn <- compile.grn.from.kegg(pwys) head(hsa.grn) ################################################### ### code chunk number 35: nbea-methods ################################################### nbea.methods() ################################################### ### code chunk number 36: nbea ################################################### nbea.res <- nbea(method="ggea", eset=all.eset, gs=hsa.gs, grn=hsa.grn) gs.ranking(nbea.res) ################################################### ### code chunk number 37: ggea-graph ################################################### par(mfrow=c(1,2)) ggea.graph( gs=hsa.gs[["hsa05217_Basal_cell_carcinoma"]], grn=hsa.grn, eset=all.eset) ggea.graph.legend() ################################################### ### code chunk number 38: combine ################################################### res.list <- list(sbea.res, nbea.res) comb.res <- comb.ea.results(res.list) ################################################### ### code chunk number 39: browse-comb (eval = FALSE) ################################################### ## ea.browse(comb.res, graph.view=hsa.grn, nr.show=5) ################################################### ### code chunk number 40: all-in-one (eval = FALSE) ################################################### ## ebrowser( meth=c("ora", "ggea"), ## exprs=exprs.file, pdat=pdat.file, fdat=fdat.file, ## org="hsa", gs=hsa.gs, grn=hsa.grn, comb=TRUE, nr.show=5)