## ----style, echo=FALSE, results='asis'----------------------------------- BiocStyle::markdown() suppressPackageStartupMessages({ library(edgeR) library(goseq) library(org.Hs.eg.db) library(GO.db) }) ## ----prostate-edgeR-input------------------------------------------------ library(edgeR) path <- system.file(package="goseq", "extdata", "Li_sum.txt") table.summary <- read.table(path, sep='\t', header=TRUE, stringsAsFactors=FALSE) counts <- table.summary[,-1] rownames(counts) <- table.summary[,1] grp <- factor(rep(c("Control","Treated"), times=c(4,3))) summarized <- DGEList(counts, lib.size=colSums(counts), group=grp) ## ----prostate-edgeR-de--------------------------------------------------- disp <- estimateCommonDisp(summarized) tested <- exactTest(disp) topTags(tested) ## ----prostate-edgeR-padj------------------------------------------------- padj <- with(tested$table, { keep <- logFC != 0 value <- p.adjust(PValue[keep], method="BH") setNames(value, rownames(tested)[keep]) }) genes <- padj < 0.05 table(genes) ## ----prostate-edgeR-pwf-------------------------------------------------- pwf <- nullp(genes,"hg19","ensGene") head(pwf) ## ----prostate-goseq-wall------------------------------------------------- GO.wall <- goseq(pwf, "hg19", "ensGene") head(GO.wall) ## ----prostate-goseq-nobias----------------------------------------------- GO.nobias <- goseq(pwf,"hg19","ensGene",method="Hypergeometric") ## ----prostate-goseq-compare, fig.width=5, fig.height=5------------------- idx <- match(GO.nobias$category, GO.wall$category) plot(log10(GO.nobias[, "over_represented_pvalue"]) ~ log10(GO.wall[idx, "over_represented_pvalue"]), xlab="Wallenius", ylab="Hypergeometric", xlim=c(-5, 0), ylim=c(-5, 0)) abline(0, 1, col="red", lwd=2)