## ----include=FALSE------------------------------------------------------- knitr::opts_chunk$set(tidy=FALSE) ## ----style, eval=TRUE, echo=FALSE, results="asis"--------------------------------------- BiocStyle::latex() options(bitmapType = "cairo") ## ----echo = FALSE, include = FALSE------------------------------------------------------ eset <- matrix(rnorm(200), ncol = 2) Treat <- factor(c("Treated","Control")) design <- model.matrix(~Treat) ## ----eval = FALSE----------------------------------------------------------------------- # # fit <- lmFit(eset, design) # fit2 <- eBayes(fit) # Error in ebayes(fit = fit, proportion = proportion, stdev.coef.lim = stdev.coef.lim, : # No residual degrees of freedom in linear model fits ## --------------------------------------------------------------------------------------- design ## --------------------------------------------------------------------------------------- head(eset) ## ----eval = FALSE----------------------------------------------------------------------- # # fit <- lmFit(eset, design) # fit2 <- eBayes(fit) # Error in ebayes(fit = fit, proportion = proportion, stdev.coef.lim = stdev.coef.lim, : # No residual degrees of freedom in linear model fits # ## --------------------------------------------------------------------------------------- d.f <- data.frame(Samples = paste("Sample", 1:10), Type = factor(rep(c("Tumor","Control"), each = 5))) d.f ## --------------------------------------------------------------------------------------- design1 <- model.matrix(~0 + Type, d.f) design1 ## --------------------------------------------------------------------------------------- design2 <- model.matrix(~ Type, d.f) design2 ## --------------------------------------------------------------------------------------- data.frame(Type = d.f$Type, design1) ## --------------------------------------------------------------------------------------- data.frame(Type = d.f$Type, design2) ## --------------------------------------------------------------------------------------- library(limma) set.seed(0xabeef) y <- matrix(rnorm(10000), 1000) ## note that for design we need to use a contrasts matrix contrast <- makeContrasts(TypeTumor - TypeControl, levels = design1) fit <- lmFit(y, design1) fit2 <- contrasts.fit(fit, contrast) fit2 <- eBayes(fit2) fit3 <- lmFit(y, design2) fit3 <- eBayes(fit3) all.equal(topTable(fit2,1), topTable(fit3,2)) ## --------------------------------------------------------------------------------------- d.f2 <- data.frame(Patient = paste("Patient", rep(1:5, 2)), Type = d.f$Type) design3 <- model.matrix(~Type + Patient, d.f2) colnames(design3) <- gsub("Type|^Patient", "", colnames(design3)) design3 ## --------------------------------------------------------------------------------------- d.f3 <- data.frame(Treatment = rep(c("Treat","Cont","Treat","Cont"), c(4,4,3,3)), Batch = rep(c("A","B"), c(8,6))) design4 <- model.matrix(~Treatment + Batch, d.f3) design4 ## ----interaction, echo = FALSE, fig.cap = "Interaction example", out.width = "0.6\\linewidth", fig.pos = "!htbp"---- library(ggplot2) set.seed(1) tmp <- data.frame(GeneX = c(rnorm(21, 4), rnorm(7, 7)), Type = factor(rep(c("WT Control","WT Treated","KO Control","KO Treated"), each = 7))) ggplot(tmp, aes(Type, GeneX)) + geom_point() ## --------------------------------------------------------------------------------------- ## let's assume that we have seven replicates per group d.f4 <- data.frame(Treatment = rep(c("Control","Treated"), times = 2, each = 7), Genotype = rep(c("WT","KO"), each = 14)) design5 <- model.matrix(~Treatment * Genotype, d.f4) design5 ## --------------------------------------------------------------------------------------- newgrp <- factor(apply(d.f4, 1, paste, collapse = "_")) design5.1 <- model.matrix(~ 0 + newgrp) colnames(design5.1) <- gsub("newgrp", "", colnames(design5.1)) design5.1 ## ----eval = FALSE----------------------------------------------------------------------- # # library(limma) # debug(lmFit) # lmFit(y, design) # ## ----eval = FALSE----------------------------------------------------------------------- # # library("affycoretools") # vennPage # drawVenn # debug(drawVenn) # ## access unexported functions with the ::: operator # affycoretools:::drawVenn # debug(affycoretools:::drawVenn) # ## ----eval = FALSE----------------------------------------------------------------------- # # eBayes(lmFit(eset, design)) # traceback() # ## --------------------------------------------------------------------------------------- print head(methods(print)) ## ----eval = FALSE----------------------------------------------------------------------- # # mod <- aov(yield ~ block + N * P + K, npk) # class(mod) # print(mod) # print.aov # ## no idea where this function lives. Use getAnywhere() # getAnywhere(print.aov)$where # debug(stats:::print.aov) # ## can now step through debugger by doing print(mod) # ## ----eval = FALSE----------------------------------------------------------------------- # # library("Biobase") # show # showMethods(show) # ## we can look at individual functions by saying what class we # ## want and setting includeDefs to TRUE # showMethods(show, class = "AnnotatedDataFrame", includeDefs = TRUE) # ## ----eval = FALSE----------------------------------------------------------------------- # # .showAnnotatedDataFrame # ## Hmmm # getAnywhere(.showAnnotatedDataFrame) # ## can we use debug()? # # df <- data.frame(x=1:6, # y=rep(c("Low", "High"),3), # z=I(LETTERS[1:6]), # row.names=paste("Sample", 1:6, sep="_")) # # ad <- AnnotatedDataFrame(data=df) # show(ad) # debug(Biobase:::.showAnnotatedDataFrame) # ## show(ad) to run through the debugger # ## ----eval = FALSE----------------------------------------------------------------------- # # showMethods(show, class = "eSet", includeDefs = TRUE) # ## just a function there. Can we debug() directly? # debug(show) # data(sample.ExpressionSet) # sample.ExpressionSet # ## ----eval = FALSE----------------------------------------------------------------------- # # trace(show, tracer = browser, signature = c(x = "eSet"), where = getNamespace("Biobase")) # sample.ExpressionSet # # ## ----sessioninfo, results="asis"-------------------------------------------------------- toLatex(sessionInfo())