## ----include=FALSE------------------------------------------------------- knitr::opts_chunk$set(tidy=FALSE) ## ----style, eval=TRUE, echo=FALSE, results="asis"------------------------ BiocStyle::latex() options(bitmapType = "cairo") ## ----eval=FALSE---------------------------------------------------------- ## ## library(knitr) ## purl(system.file("doc/RefactoredAffycoretools.Rnw", package="affycoretools")) ## ## ----start--------------------------------------------------------------- suppressMessages(library(affycoretools)) data(sample.ExpressionSet) eset <- sample.ExpressionSet eset ## ----include=FALSE------------------------------------------------------- featureNames(eset) <- gsub("/", "_", featureNames(eset)) ## ----model--------------------------------------------------------------- suppressMessages(library(limma)) pd <- pData(phenoData(eset)) design <- model.matrix(~0+type+sex, pd) colnames(design) <- gsub("type|sex", "", colnames(design)) contrast <- matrix(c(1,-1,0)) colnames(contrast) <- "Case vs control" fit <- lmFit(eset, design) fit2 <- contrasts.fit(fit, contrast) fit2 <- eBayes(fit) topTable(fit2, 1)[,1:4] ## ----annotate------------------------------------------------------------ suppressMessages(library(hgu95av2.db)) gns <- select(hgu95av2.db, featureNames(eset), c("ENTREZID","SYMBOL","GENENAME")) ## There are one-to many mappings here, so we just ## removed duplicates in a very naive way. gns <- gns[!duplicated(gns[,1]),] fit2$genes <- gns topTable(fit2, 1)[,1:3] ## ----output-------------------------------------------------------------- suppressMessages(library(ReportingTools)) htab <- HTMLReport("afile", "My cool results") publish(topTable(fit2, 1), htab) finish(htab) ## ----rtools-------------------------------------------------------------- htab <- HTMLReport("afile2", "My cool results, ReportingTools style") publish(fit2, htab, eset, factor = pd$type, coef = 1, n = 10) finish(htab) ## ----output2------------------------------------------------------------- d.f <- topTable(fit2, 2) out <- makeImages(df = d.f, eset = eset, grp.factor = pd$type, design = design, contrast = contrast, colind = 1, repdir = ".") htab <- HTMLReport("afile3", "My cool results, affycoretools style") publish(out$df, htab, .mofifyDF = list(entrezLinks, affyLinks)) finish(htab) ## ----param2-------------------------------------------------------------- grps <- factor(apply(pd[,1:2], 1, paste, collapse = "_")) design <- model.matrix(~0+grps) colnames(design) <- gsub("grps", "", colnames(design)) contrast <- matrix(c(1,-1,0,0, 0,0,1,-1, 1,-1,-1,1), ncol = 3) colnames(contrast) <- c("Female_Case vs Female_Control", "Male_Case vs Male_Control", "Interaction") fit <- lmFit(eset, design) fit2 <- contrasts.fit(fit, contrast) fit2 <- eBayes(fit2) fit2$genes <- gns ## ----output3------------------------------------------------------------- ## get a list containing the output for each comparison out <- lapply(1:3, function(x) topTable(fit2, x)) ## process the output to add images htab <- lapply(1:3, function(x){ tmp <- HTMLReport(gsub("_", " ", colnames(contrast)[x]), colnames(contrast)[x], "./reports") tmp2 <- makeImages(out[[x]], eset, grps, design, contrast, x) publish(tmp2$df, tmp, .modifyDF = list(affyLinks, entrezLinks)) finish(tmp) return(tmp) }) ## Now make an index.html file to contain links to the various comps idx <- HTMLReport("index", "Links to our stuff") publish(hwriter::hwrite("Univariate comparisons", br = TRUE, header = 2), idx) publish(Link(htab), idx) finish(idx) ## ----makevenn------------------------------------------------------------ collist <- list(1:2) venn <- makeVenn(fit2, contrast, design, collist = collist, adj.meth = "none") vennlnk <- vennPage(venn, "venn_diagram", "Venn diagram") ## ----venn, fig.cap="Venn diagram"---------------------------------------- vennDiagram(venn[[1]]$venncounts, cex = 0.9) ## ----addvenn------------------------------------------------------------- idx <- HTMLReport("index","Links to our stuff") publish(hwriter::hwrite("Univariate comparisons", br = TRUE, header = 2), idx) publish(Link(htab), idx) publish(hwriter::hwrite("Venn diagrams", br = TRUE, header = 2), idx) publish(Link("Venn1", vennlnk), idx) finish(idx) ## ----sessioninfo, results="asis"----------------------------------------- toLatex(sessionInfo())