## ----setup, include=FALSE----------------------------------------------------- knitr::opts_chunk$set(echo = TRUE) library(limma) library(TeachingDemos) options(digits=2) set.seed(20190429) col.fit <- "#56B4E9" ## ----echo=FALSE, out.width='80%'---------------------------------------------- knitr::include_graphics(here::here("vignettes", "Table1.png")) ## ----echo=FALSE, out.width='80%'---------------------------------------------- knitr::include_graphics(here::here("vignettes", "Table2.png")) ## ----echo=FALSE, out.width='80%'---------------------------------------------- knitr::include_graphics(here::here("vignettes", "Table3.png")) ## ----echo=FALSE, out.width='80%'---------------------------------------------- knitr::include_graphics(here::here("vignettes", "Table4.png")) ## ----echo=FALSE, out.width='80%'---------------------------------------------- knitr::include_graphics(here::here("vignettes", "Table5.png")) ## ----echo=FALSE, out.width='80%'---------------------------------------------- knitr::include_graphics(here::here("vignettes", "Table6.png")) ## ---- echo=FALSE-------------------------------------------------------------- id = paste0("MOUSE", 1:6) mouse <- id age = c(1,2,3,4,5,6) expression = age/2 + 2 + rnorm(n=6, sd=0.1) data <- data.frame(expression, mouse, age) options(digits=3) data options(digits=2) ## ---- eval=TRUE, echo=FALSE--------------------------------------------------- ploti <- 3 ## ---- design3, eval=FALSE, echo=FALSE----------------------------------------- # cat("Output for Figure", ploti) # model.matrix(~age) # fit <- lm(expression~age) # fit ## ---- include=FALSE, echo=FALSE----------------------------------------------- cat("Output for Figure", ploti) model.matrix(~age) fit <- lm(expression~age) fit plot.name <- paste0("fig", ploti, ".png") png(plot.name) model.name <- paste0("E(y)=", round(fit$coef[1],2), "+", round(fit$coef[2],2), "x") plot(age, expression, pch=16, ylab="", xlab="", xlim=c(0,7), ylim=c(0,7)) abline(fit$coef[1], fit$coef[2], col=col.fit, lty=1) abline(v=0, col="lightgrey") title(ylab="expression (y)", xlab="age (x)", line=2) title(main=model.name) dev.off() ## ---- design3, eval=TRUE, echo=TRUE, include=TRUE----------------------------- cat("Output for Figure", ploti) model.matrix(~age) fit <- lm(expression~age) fit ## ---- eval=TRUE, echo=FALSE--------------------------------------------------- ploti <- ploti+1 ## ---- design4, eval=FALSE, echo=FALSE----------------------------------------- # cat("Output for Figure", ploti) # model.matrix(~0+age) # fit <- lm(expression~0+age) # fit ## ---- include=FALSE, echo=FALSE----------------------------------------------- cat("Output for Figure", ploti) model.matrix(~0+age) fit <- lm(expression~0+age) fit plot.name <- paste0("fig", ploti, ".png") png(plot.name) model.name <- paste0("E(y)=", round(fit$coef[1],2), "x") plot(age, expression, pch=16, ylab="", xlab="", xlim=c(0,7), ylim=c(0,7)) abline(0, fit$coef[1], col=col.fit, lty=1) abline(v=0, col="lightgrey") title(ylab="expression (y)", xlab="age (x)", line=2) title(main=model.name) dev.off() ## ---- design4, eval=TRUE, echo=TRUE, include=TRUE----------------------------- cat("Output for Figure", ploti) model.matrix(~0+age) fit <- lm(expression~0+age) fit ## ---- echo=FALSE-------------------------------------------------------------- group <- as.factor(rep(c("HEALTHY", "SICK"), each=3)) data <- data.frame(expression=round(expression,2), mouse=id, group=as.character(group)) options(digits=3) data options(digits=2) ## ---- eval=TRUE, echo=FALSE--------------------------------------------------- ploti <- ploti+1 ## ---- design5, eval=FALSE, echo=FALSE----------------------------------------- # cat("Output for Figure", ploti) # model.matrix(~group) # fit <- lm(expression~group) # fit ## ---- include=FALSE, echo=FALSE----------------------------------------------- cat("Output for Figure", ploti) model.matrix(~group) fit <- lm(expression~group) fit plot.name <- paste0("fig", ploti, ".png") png(plot.name) model.name <- paste0("E(y)=", round(fit$coef[1],2), "+", round(fit$coef[2],2), "x") plot(as.numeric(group), expression, pch=16, ylab="", xlab="", xlim=c(0,3), ylim=c(0,7), xaxt="n") abline(h=fit$coef[1], col=col.fit, lty=1) abline(h=fit$coef[1]+fit$coef[2], col=col.fit, lty=2) parameters <- c("", "(x=1)") axis(side=1, at=1:nlevels(group), labels=paste(levels(group), parameters)) title(xlab="group", line=3) title(ylab="expression (y)", line=2) title(main=model.name) dev.off() ## ---- design5, eval=TRUE, echo=TRUE, include=TRUE----------------------------- cat("Output for Figure", ploti) model.matrix(~group) fit <- lm(expression~group) fit ## ---- eval=TRUE, echo=FALSE--------------------------------------------------- ploti <- ploti+1 ## ---- design6, eval=FALSE, echo=FALSE----------------------------------------- # cat("Output for Figure", ploti) # model.matrix(~0+group) # fit <- lm(expression~0+group) # fit ## ---- include=FALSE, echo=FALSE----------------------------------------------- cat("Output for Figure", ploti) model.matrix(~0+group) fit <- lm(expression~0+group) fit plot.name <- paste0("fig", ploti, ".png") png(plot.name) model.name <- paste0("E(y)=", round(fit$coef[1],2), "x1+", round(fit$coef[2],2), "x2") plot(as.numeric(group), expression, pch=16, ylab="", xlab="", xlim=c(0,3), ylim=c(0,7), xaxt="n") abline(h=fit$coef[1], col=col.fit, lty=1) abline(h=fit$coef[2], col=col.fit, lty=1) parameters <- c("(x1=1)", "(x2=1)") axis(side=1, at=1:nlevels(group), labels=paste(levels(group), parameters), cex.axis=0.7) title(xlab="group", line=3) title(ylab="expression (y)", line=2) title(main=model.name) dev.off() ## ---- design6, eval=TRUE, echo=TRUE, include=TRUE----------------------------- cat("Output for Figure", ploti) model.matrix(~0+group) fit <- lm(expression~0+group) fit ## ----------------------------------------------------------------------------- design <- model.matrix(~0+group) makeContrasts(groupSICK-groupHEALTHY, levels=colnames(design)) ## ---- echo=FALSE-------------------------------------------------------------- treatment = rep(c("CTL", "I", "II", "III"), each=3) treatment <- factor(treatment, levels=unique(treatment)) expression <- rep(c(0,1,2,4), each=3) + 1 expression = expression + rnorm(n=length(treatment), sd=0.1) id = paste0("MOUSE", 1:length(treatment)) data <- data.frame(expression=round(expression,2), mouse=id, treatment=treatment) options(digits=3) data options(digits=2) ## ---- eval=TRUE, echo=FALSE--------------------------------------------------- ploti <- ploti+1 ## ---- design7, eval=FALSE, echo=FALSE----------------------------------------- # cat("Output for Figure", ploti) # model.matrix(~treatment) # fit <- lm(expression~treatment) # fit ## ---- include=FALSE, echo=FALSE----------------------------------------------- cat("Output for Figure", ploti) model.matrix(~treatment) fit <- lm(expression~treatment) fit plot.name <- paste0("fig", ploti, ".png") png(plot.name) model.name <- paste0("E(y)=", round(fit$coef[1],2), "+", round(fit$coef[2],2), "x1+", round(fit$coef[3],2), "x2+", round(fit$coef[4],2), "x3") plot(as.numeric(treatment), expression, pch=16, ylab="", xlab="", xlim=c(0.5,4.5), ylim=c(0,6), xaxt="n") abline(h=fit$coef[1], col=col.fit, lty=1) abline(h=fit$coef[1]+fit$coef[2], col=col.fit, lty=2) abline(h=fit$coef[1]+fit$coef[3], col=col.fit, lty=2) abline(h=fit$coef[1]+fit$coef[4], col=col.fit, lty=2) parameters <- c("", "(x1=1)", "(x2=1)", "(x3=1)") axis(side=1, at=1:nlevels(treatment), labels=paste(levels(treatment), parameters), cex.axis=0.7) title(xlab="treatment", line=3) title(ylab="expression (y)", line=2) title(main=model.name) dev.off() # save this fit for later use fit.mr <- fit ## ---- design7, eval=TRUE, echo=TRUE, include=TRUE----------------------------- cat("Output for Figure", ploti) model.matrix(~treatment) fit <- lm(expression~treatment) fit ## ---- eval=TRUE, echo=FALSE--------------------------------------------------- ploti <- ploti+1 ## ---- design8, eval=FALSE, echo=FALSE----------------------------------------- # cat("Output for Figure", ploti) # model.matrix(~0+treatment) # fit <- lm(expression~0+treatment) # fit ## ---- include=FALSE, echo=FALSE----------------------------------------------- cat("Output for Figure", ploti) model.matrix(~0+treatment) fit <- lm(expression~0+treatment) fit plot.name <- paste0("fig", ploti, ".png") png(plot.name) model.name <- paste0("E(y)=", round(fit$coef[1],2), "x0+", round(fit$coef[2],2), "x1+", round(fit$coef[3],2), "x2+", round(fit$coef[4],2), "x3") plot(as.numeric(treatment), expression, pch=16, ylab="", xlab="", xlim=c(0.5,4.5), ylim=c(0,6), xaxt="n") abline(h=fit$coef[1], col=col.fit, lty=1) abline(h=fit$coef[2], col=col.fit, lty=1) abline(h=fit$coef[3], col=col.fit, lty=1) abline(h=fit$coef[4], col=col.fit, lty=1) parameters <- c("(x0=1)", "(x1=1)", "(x2=1)", "(x3=1)") axis(side=1, at=1:nlevels(treatment), labels=paste(levels(treatment), parameters), cex.axis=0.7) title(xlab="treatment", line=3) title(ylab="expression (y)", line=2) title(main=model.name) dev.off() ## ---- echo=FALSE-------------------------------------------------------------- design <- model.matrix(~0+treatment) ## ----------------------------------------------------------------------------- contrasts <- makeContrasts( treatmentI-treatmentCTL, treatmentII-treatmentCTL, treatmentIII-treatmentCTL, treatmentII-treatmentI, treatmentIII-treatmentI, treatmentIII-treatmentII, levels=colnames(design)) colnames(contrasts) <- abbreviate(colnames(contrasts)) contrasts ## ---- design8, eval=TRUE, echo=TRUE, include=TRUE----------------------------- cat("Output for Figure", ploti) model.matrix(~0+treatment) fit <- lm(expression~0+treatment) fit ## ----------------------------------------------------------------------------- makeContrasts((treatmentI+treatmentII+treatmentIII)/3-treatmentCTL, levels=colnames(design)) ## ----------------------------------------------------------------------------- makeContrasts((treatmentCTL+treatmentIII)/2-(treatmentI+treatmentII)/2, levels=colnames(design)) ## ---- warning=FALSE----------------------------------------------------------- design <- model.matrix(~treatment) makeContrasts(treatmentIII-treatmentI-treatmentII, levels=colnames(design)) ## ---- echo=FALSE-------------------------------------------------------------- treat1 <- as.factor(c(0,0,0,1,1,1,0,0,0,1,1,1)) treat2 <- as.factor(c(0,0,0,0,0,0,1,1,1,1,1,1)) levels(treat1) <-levels(treat2) <- c("NO", "YES") data <- data.frame(expression=expression, mouse=id, treat1=treat1, treat2=treat2) options(digits=3) data options(digits=2) ## ---- eval=TRUE, echo=FALSE--------------------------------------------------- ploti <- ploti+1 ## ---- design9, eval=FALSE, echo=FALSE----------------------------------------- # cat("Output for Figure", ploti) # model.matrix(~treat1*treat2) # fit <- lm(expression~treat1*treat2) # fit ## ---- include=FALSE, echo=FALSE----------------------------------------------- cat("Output for Figure", ploti) model.matrix(~treat1*treat2) fit <- lm(expression~treat1*treat2) fit plot.name <- paste0("fig", ploti, ".png") png(plot.name) model.name <- paste0("E(y)=", round(fit$coef[1],2), "+", round(fit$coef[2],2), "x1+", round(fit$coef[3],2), "x2+", round(fit$coef[4],2), "x1x2") plot(as.numeric(treatment), expression, pch=16, ylab="", xlab="", xlim=c(0.5,4.5), ylim=c(0,6), xaxt="n") abline(h=fit$coef[1], col=col.fit, lty=1) abline(h=fit$coef[1]+fit$coef[2], col=col.fit, lty=2) abline(h=fit$coef[1]+fit$coef[3], col=col.fit, lty=2) abline(h=fit$coef[1]++fit$coef[2]+fit$coef[3]+fit$coef[4], col="darkblue", lty=3) parameters <- c("", "(x1=1)", "(x2=1)", "(x1=1,x2=1)") axis(side=1, at=1:nlevels(treatment), labels=paste(c("CTL", "I", "II", "I+II"), parameters), cex.axis=0.7) title(xlab="treatment", line=3) title(ylab="expression (y)", line=2) title(main=model.name) dev.off() ## ---- design9, eval=TRUE, echo=TRUE, include=TRUE----------------------------- cat("Output for Figure", ploti) model.matrix(~treat1*treat2) fit <- lm(expression~treat1*treat2) fit ## ---- eval=TRUE, echo=FALSE--------------------------------------------------- ploti <- ploti+1 ## ---- design10, eval=FALSE, echo=FALSE---------------------------------------- # cat("Output for Figure", ploti) # model.matrix(~treat1+treat2) # fit <- lm(expression~treat1+treat2) # fit ## ---- include=FALSE, echo=FALSE----------------------------------------------- cat("Output for Figure", ploti) model.matrix(~treat1+treat2) fit <- lm(expression~treat1+treat2) fit plot.name <- paste0("fig", ploti, ".png") png(plot.name) model.name <- paste0("E(y)=", round(fit$coef[1],2), "+", round(fit$coef[2],2), "x1+", round(fit$coef[3],2), "x2") plot(as.numeric(treatment), expression, pch=16, ylab="", xlab="", xlim=c(0.5,4.5), ylim=c(0,6), xaxt="n") abline(h=fit$coef[1], col=col.fit, lty=1) abline(h=fit$coef[1]+fit$coef[2], col=col.fit, lty=2) abline(h=fit$coef[1]+fit$coef[3], col=col.fit, lty=2) abline(h=fit$coef[1]+fit$coef[2]+fit$coef[3], col="darkblue", lty=3) parameters <- c("", "(x1=1)", "(x2=1)", "(x1=1,x2=1)") axis(side=1, at=1:nlevels(treatment), labels=paste(c("CTL", "I", "II", "I+II"), parameters), cex.axis=0.7) title(xlab="treatment", line=3) title(ylab="expression (y)", line=2) title(main=model.name) dev.off() ## ---- design10, eval=TRUE, echo=TRUE, include=TRUE---------------------------- cat("Output for Figure", ploti) model.matrix(~treat1+treat2) fit <- lm(expression~treat1+treat2) fit ## ---- echo=FALSE-------------------------------------------------------------- tissue <- treat1 levels(tissue) <- c("LUNG", "BRAIN") cells <- treat2 levels(cells) <- c("B", "T") options(digits=3) data.frame(expression, id, tissue, cells) options(digits=2) ## ---- echo=FALSE-------------------------------------------------------------- group <- paste(tissue, cells, sep="_") group <- factor(group, levels=unique(group)) options(digits=3) data.frame(expression, id, tissue, cells, group) options(digits=2) design <- model.matrix(~0+group) ## ----------------------------------------------------------------------------- contrasts <- makeContrasts( BVsT=(groupLUNG_B+groupBRAIN_B)/2-(groupLUNG_T+groupBRAIN_T)/2, LungVsBrain=(groupLUNG_B+groupLUNG_T)/2-(groupBRAIN_B+groupBRAIN_T)/2, BVsT_Lung=groupLUNG_B-groupLUNG_T, BVsT_Brain=groupBRAIN_B-groupBRAIN_T, LungVsBrain_B=groupLUNG_B-groupBRAIN_B, LungVsBrain_T=groupLUNG_T-groupBRAIN_T, levels=colnames(design)) rownames(contrasts) <- gsub("group", "", rownames(contrasts)) contrasts ## ---- eval=TRUE, echo=FALSE--------------------------------------------------- ploti <- ploti+1 ## ---- design11, eval=FALSE, echo=FALSE---------------------------------------- # cat("Output for Figure", ploti) # model.matrix(~0+group) # fit <- lm(expression~0+group) # fit ## ---- include=FALSE, echo=FALSE----------------------------------------------- cat("Output for Figure", ploti) model.matrix(~0+group) fit <- lm(expression~0+group) fit plot.name <- paste0("fig", ploti, ".png") png(plot.name) model.name <- paste0("E(y)=", round(fit$coef[1],2), "x1+", round(fit$coef[2],2), "x2+", round(fit$coef[3],2), "x3+", round(fit$coef[4],2), "x4") plot(as.numeric(group), expression, pch=16, ylab="", xlab="", xlim=c(0.5,4.5), ylim=c(0,6), xaxt="n") abline(h=fit$coef[1], col=col.fit, lty=1) abline(h=fit$coef[2], col=col.fit, lty=1) abline(h=fit$coef[3], col=col.fit, lty=1) abline(h=fit$coef[4], col=col.fit, lty=1) parameters <- c("(x1=1)", "(x2=1)", "(x3=1)", "(x4=1)") axis(side=1, at=1:nlevels(group), labels=paste(levels(group), parameters), cex.axis=0.6) title(xlab="group", line=3) title(ylab="expression (y)", line=2) title(main=model.name) dev.off() C <- fit$coef%*%contrasts ## ---- design11, eval=TRUE, echo=TRUE, include=TRUE---------------------------- cat("Output for Figure", ploti) model.matrix(~0+group) fit <- lm(expression~0+group) fit ## ---- echo=FALSE-------------------------------------------------------------- group <- as.factor(rep(LETTERS[1:4], each=3)) lane <- rep(c("L1", "L2"), c(6,6)) lane <- as.factor(sample(lane, length(group), replace=FALSE)) technician <- c("I", "II") technician <- as.factor(sample(technician, length(group), replace=TRUE)) options(digits=3) data.frame(expression, id, group, lane, technician) options(digits=2) ## ---- eval=TRUE, echo=FALSE--------------------------------------------------- ploti <- ploti+1 ## ---- design12, eval=FALSE, echo=FALSE---------------------------------------- # cat("Output for Figure", ploti) # model.matrix(~0+group+lane+technician) # fit <- lm(expression~0+group+lane+technician) # fit ## ---- include=FALSE, echo=FALSE----------------------------------------------- cat("Output for Figure", ploti) model.matrix(~0+group+lane+technician) fit <- lm(expression~0+group+lane+technician) fit plot.name <- paste0("fig", ploti, ".png") png(plot.name) model.name <- paste0("E(y)=", round(fit$coef[1],2), "x1+", round(fit$coef[2],2), "x2+", round(fit$coef[3],2), "x3+", round(fit$coef[4],2), "x4+", round(fit$coef[5],2), "x5+", round(fit$coef[6],2), "x6") plot(jitter(as.numeric(group)), expression, pch=16, ylab="", xlab="", xlim=c(0.5,4.5), ylim=c(0,6), xaxt="n", type="n") text(jitter(as.numeric(group)), expression, label=lane, col=as.numeric(technician)) abline(h=fit$coef[1], col=col.fit, lty=1) abline(h=fit$coef[2], col=col.fit, lty=1) abline(h=fit$coef[3], col=col.fit, lty=1) abline(h=fit$coef[4], col=col.fit, lty=1) parameters <- c("(x1=1)", "(x2=1)", "(x3=1)", "(x4=1)") axis(side=1, at=1:nlevels(group), labels=paste(levels(group), parameters), cex.axis=0.6) title(xlab="group", line=3) title(ylab="expression (y)", line=2) title(main=model.name) dev.off() ## ---- design12, eval=TRUE, echo=TRUE, include=TRUE---------------------------- cat("Output for Figure", ploti) model.matrix(~0+group+lane+technician) fit <- lm(expression~0+group+lane+technician) fit ## ---- echo=FALSE-------------------------------------------------------------- batch <- rep(c("B1", "B2"), c(6,6)) options(digits=3) data.frame(expression, id, group, batch) options(digits=2) ## ---- echo=FALSE-------------------------------------------------------------- treatment3 <- treat2 levels(treatment3) <- c("X", "Y") df <- data.frame(expression, id, treatment=treatment3) df <- cbind(df, timepoint=paste0("T", rep(c(1,2), each=3))) df$id <- paste0("MOUSE", c(rep(1:3, 2), rep(4:6, 2))) id <- df$id treatment <- df$treatment timepoint <- df$timepoint options(digits=3) df options(digits=2) ## ----------------------------------------------------------------------------- design <- model.matrix(~0+id) design <- cbind(design, X= treatment=="X" & timepoint=="T2") design <- cbind(design, Y= treatment=="Y" & timepoint=="T2") ## ---- eval=TRUE, echo=FALSE--------------------------------------------------- ploti <- ploti+1 ## ---- design13, eval=FALSE, echo=FALSE---------------------------------------- # cat("Output for Figure", ploti) # design # fit <- lm(expression~0+design) # fit ## ---- include=FALSE, echo=FALSE----------------------------------------------- cat("Output for Figure", ploti) design fit <- lm(expression~0+design) fit plot.name <- paste0("fig", ploti, ".png") png(plot.name) model.name <- paste0("E(y)=", round(fit$coef[1],2), "x1+", round(fit$coef[2],2), "x2+", round(fit$coef[3],2), "x3+", round(fit$coef[4],2), "x4+", round(fit$coef[5],2), "x5+", round(fit$coef[6],2), "x6+", round(fit$coef[7],2), "x7+", round(fit$coef[8],2), "x8") plot(as.numeric(group), expression, ylab="", xlab="", xlim=c(0.5,4.5), ylim=c(0,6), xaxt="n", type="n") abline(h=mean(fit$coef[1:3]), col="salmon", lwd=3) abline(h=mean(fit$coef[4:6]), col="turquoise", lwd=3) abline(h=fit$coef[1], col=1, lty=1) abline(h=fit$coef[2], col=1, lty=1) abline(h=fit$coef[3], col=1, lty=1) abline(h=fit$coef[4], col=1, lty=1) abline(h=fit$coef[5], col=1, lty=1) abline(h=fit$coef[6], col=1, lty=1) abline(h=mean(fit$coef[1:3])+fit$coef[7], col="salmon", lwd=3, lty=2) abline(h=mean(fit$coef[4:6])+fit$coef[8], col="turquoise", lwd=3, lty=2) text(jitter(as.numeric(group)), expression, label=substr(id, 6,6)) group <- as.factor(paste(treatment, timepoint)) axis(side=1, at=1:4, labels=levels(group)) title(ylab="expression (y)", line=2) title(main=model.name, cex.main=0.8) dev.off() ## ---- design13, eval=TRUE, echo=TRUE, include=TRUE---------------------------- cat("Output for Figure", ploti) design fit <- lm(expression~0+design) fit ## ---- echo=FALSE-------------------------------------------------------------- group <- as.factor(paste(treatment, timepoint, sep="_")) df <- data.frame(df, group) options(digits=3) df options(digits=2) design <- model.matrix(~0+group) ## ---- echo=FALSE-------------------------------------------------------------- options(digits=1) ## ----------------------------------------------------------------------------- cor <- duplicateCorrelation(expression, design, block=id) cor$consensus.correlation ## ---- echo=FALSE-------------------------------------------------------------- options(digits=2) ## ----------------------------------------------------------------------------- fit <- lmFit(object=expression, design=design, block=id, correlation=cor$consensus.correlation) ## ---- eval=TRUE, echo=FALSE--------------------------------------------------- ploti <- ploti+1 ## ---- design14, eval=FALSE, echo=FALSE---------------------------------------- # cat("Output for Figure", ploti) # model.matrix(~0+group) # fit$coef ## ---- include=FALSE, echo=FALSE----------------------------------------------- cat("Output for Figure", ploti) model.matrix(~0+group) fit$coef plot.name <- paste0("fig", ploti, ".png") png(plot.name) model.name <- paste0("E(y)=", round(fit$coef[1],2), "x1+", round(fit$coef[2],2), "x2+", round(fit$coef[3],2), "x3+", round(fit$coef[4],2), "x4") plot(as.numeric(group), expression, ylab="", xlab="", xlim=c(0.5,4.5), ylim=c(0,6), xaxt="n", type="n") abline(h=fit$coef[1], col=col.fit, lty=1) abline(h=fit$coef[2], col=col.fit, lty=1) abline(h=fit$coef[3], col=col.fit, lty=1) abline(h=fit$coef[4], col=col.fit, lty=1) text(jitter(as.numeric(group)), expression, label=substr(id, 6,6)) group <- as.factor(paste(treatment, timepoint)) axis(side=1, at=1:4, labels=levels(group)) title(ylab="expression (y)", line=2) title(main=model.name, cex.main=0.8) dev.off() ## ---- include=FALSE, echo=FALSE----------------------------------------------- options(digits=3) ## ---- design14, eval=TRUE, echo=TRUE, include=TRUE---------------------------- cat("Output for Figure", ploti) model.matrix(~0+group) fit$coef ## ---- include=FALSE, echo=FALSE----------------------------------------------- options(digits=2) ## ----------------------------------------------------------------------------- contrasts <- makeContrasts( X_T2vsT1=groupX_T2-groupX_T1, Y_T2vsT1=groupY_T2-groupY_T1, XvsY=(groupX_T2-groupX_T1)-(groupY_T2-groupY_T1), levels=colnames(design)) contrasts ## ---- eval=FALSE-------------------------------------------------------------- # fit <- contrasts.fit(fit, contrasts) ## ---- echo=FALSE-------------------------------------------------------------- df <- data.frame(expression, id, treatment, timepoint) df$id <- paste0("MOUSE", 1:12) df$timepoint <- rep(1:2, each=3) colnames(df)[4] <- "time" options(digits=3) df options(digits=2) expression <- df$expression treatment <- df$treatment time <- df$time ## ---- eval=TRUE, echo=FALSE--------------------------------------------------- ploti <- ploti+1 ## ---- design15, eval=FALSE, echo=FALSE---------------------------------------- # cat("Output for Figure", ploti) # model.matrix(~0+treatment+treatment:time) # fit <- lm(expression~0+treatment+treatment:time) # fit ## ---- include=FALSE, echo=FALSE----------------------------------------------- cat("Output for Figure", ploti) model.matrix(~0+treatment+treatment:time) fit <- lm(expression~0+treatment+treatment:time) fit plot.name <- paste0("fig", ploti, ".png") png(plot.name) model.name <- paste0("E(y)=", round(fit$coef[1],2), "x1+", round(fit$coef[2],2), "x2+", round(fit$coef[3],2), "x1.t+", round(fit$coef[4],2), "x2.t") plot(time, expression, ylab="", xlab="", xlim=c(0.5,2.5), ylim=c(0,6), xaxt="n", type="n") text(time, expression, label=treatment) abline(fit$coef[1], fit$coef[3], col=col.fit) abline(fit$coef[2], fit$coef[4], col=col.fit) axis(side=1, at=1:2, labels=unique(time)) title(xlab="time (t)", line=2) title(ylab="expression (y)", line=2) title(main=model.name, cex.main=0.8) dev.off() ## ---- design15, eval=TRUE, echo=TRUE, include=TRUE---------------------------- cat("Output for Figure", ploti) model.matrix(~0+treatment+treatment:time) fit <- lm(expression~0+treatment+treatment:time) fit ## ---- echo=FALSE-------------------------------------------------------------- TIME <- rep(1:10, each=2) EXPRS <- 2*(sin(TIME-1)+1) EXPRS <- EXPRS + rnorm(n=length(EXPRS), sd=0.2) ## ---- echo=FALSE-------------------------------------------------------------- i <- TIME<7 time <- TIME[i] expression <- EXPRS[i] data <- data.frame(expression, mouse=paste0("MOUSE", 1:sum(i)), time) data ## ---- eval=TRUE, echo=FALSE--------------------------------------------------- ploti <- ploti+1 ## ---- design16, eval=FALSE, echo=FALSE---------------------------------------- # cat("Output for Figure", ploti) # model.matrix(~time) # fit <- lm(expression~time) # fit ## ---- include=FALSE, echo=FALSE----------------------------------------------- cat("Output for Figure", ploti) model.matrix(~time) fit <- lm(expression~time) fit plot.name <- paste0("fig", ploti, ".png") png(plot.name) model.name <- paste0("E(y)=", round(fit$coef[1],2), "+", round(fit$coef[2],2), "t") plot(time, expression, pch=16, ylab="", xlab="", xlim=c(0.5,6.5), ylim=c(-0.5,4.5)) time_grid <- seq(from=0, to=11, by=0.1) preds <- predict(fit, newdata=list(time=time_grid)) lines(time_grid, preds, col=col.fit) title(xlab="time (t)", line=2) title(ylab="expression (y)", line=2) title(main=model.name, cex.main=0.8) dev.off() ## ---- design16, eval=TRUE, echo=TRUE, include=TRUE---------------------------- cat("Output for Figure", ploti) model.matrix(~time) fit <- lm(expression~time) fit ## ---- eval=TRUE, echo=FALSE--------------------------------------------------- ploti <- ploti+1 ## ---- design17, eval=FALSE, echo=FALSE---------------------------------------- # cat("Output for Figure", ploti) # design <- model.matrix(~poly(time, degree=2, raw=TRUE)) # colnames(design)[2:3] <- c("time", "time2") # design # fit <- lm(expression~poly(time, degree=2, raw=TRUE)) # fit ## ---- include=FALSE, echo=FALSE----------------------------------------------- cat("Output for Figure", ploti) design <- model.matrix(~poly(time, degree=2, raw=TRUE)) colnames(design)[2:3] <- c("time", "time2") design fit <- lm(expression~poly(time, degree=2, raw=TRUE)) fit plot.name <- paste0("fig", ploti, ".png") png(plot.name) model.name <- paste0("E(y)=", round(fit$coef[1],2), "+", round(fit$coef[2],2), "t+", round(fit$coef[3],2), "t^2") plot(time, expression, pch=16, ylab="", xlab="", xlim=c(0.5,6.5), ylim=c(-0.5,4.5)) preds <- predict(fit, newdata=list(time=time_grid)) lines(time_grid, preds, col=col.fit) title(xlab="time (t)", line=2) title(ylab="expression (y)", line=2) title(main=model.name, cex.main=0.8) dev.off() peak.time <- -fit$coef[2]/(2*fit$coef[3]) peak.exprs <- fit$coef[1]+fit$coef[2]*peak.time+fit$coef[3]*peak.time^2 ## ---- design17, eval=TRUE, echo=TRUE, include=TRUE---------------------------- cat("Output for Figure", ploti) design <- model.matrix(~poly(time, degree=2, raw=TRUE)) colnames(design)[2:3] <- c("time", "time2") design fit <- lm(expression~poly(time, degree=2, raw=TRUE)) fit ## ---- echo=FALSE-------------------------------------------------------------- i <- TIME<9 time <- TIME[i] expression <- EXPRS[i] data <- data.frame(expression, mouse=paste0("MOUSE", 1:sum(i)), time) data ## ---- eval=TRUE, echo=FALSE--------------------------------------------------- ploti <- ploti+1 ## ---- design18, eval=FALSE, echo=FALSE---------------------------------------- # cat("Output for Figure", ploti) # design <- model.matrix(~poly(time, degree=3, raw=TRUE)) # colnames(design)[2:4] <- c("time", "time2", "time3") # design # fit <- lm(expression~poly(time, degree=3, raw=TRUE)) # fit ## ---- include=FALSE, echo=FALSE----------------------------------------------- cat("Output for Figure", ploti) design <- model.matrix(~poly(time, degree=3, raw=TRUE)) colnames(design)[2:4] <- c("time", "time2", "time3") design fit <- lm(expression~poly(time, degree=3, raw=TRUE)) fit plot.name <- paste0("fig", ploti, ".png") png(plot.name) model.name <- paste0("E(y)=", round(fit$coef[1],2), "+", round(fit$coef[2],2), "t+", round(fit$coef[3],2), "t^2+", round(fit$coef[4],2), "t^3") plot(time, expression, pch=16, ylab="", xlab="", xlim=c(0.5,8.5), ylim=c(-0.5,4.5)) preds <- predict(fit, newdata=list(time=time_grid)) lines(time_grid, preds, col=col.fit) title(xlab="time (t)", line=2) title(ylab="expression (y)", line=2) title(main=model.name, cex.main=0.8) dev.off() ## ---- design18, eval=TRUE, echo=TRUE, include=TRUE---------------------------- cat("Output for Figure", ploti) design <- model.matrix(~poly(time, degree=3, raw=TRUE)) colnames(design)[2:4] <- c("time", "time2", "time3") design fit <- lm(expression~poly(time, degree=3, raw=TRUE)) fit ## ---- echo=FALSE-------------------------------------------------------------- i <- TIME<11 time <- TIME[i] expression <- EXPRS[i] data <- data.frame(expression, mouse=paste0("MOUSE", 1:sum(i)), time) data ## ----------------------------------------------------------------------------- cycle <- 6 sinphase <- sin(2*pi*time/cycle) cosphase <- cos(2*pi*time/cycle) ## ----------------------------------------------------------------------------- design <- model.matrix(~sinphase+cosphase) ## ---- eval=TRUE, echo=FALSE--------------------------------------------------- ploti <- ploti+1 ## ---- design19, eval=FALSE, echo=FALSE---------------------------------------- # cat("Output for Figure", ploti) # design # fit <- lm(expression~sin(2*pi*time/cycle)+cos(2*pi*time/cycle)) # fit ## ---- include=FALSE, echo=FALSE----------------------------------------------- cat("Output for Figure", ploti) design fit <- lm(expression~sin(2*pi*time/cycle)+cos(2*pi*time/cycle)) fit plot.name <- paste0("fig", ploti, ".png") png(plot.name) model.name <- paste0("E(y)=", round(fit$coef[1],2), "+", round(fit$coef[2],2), "sin(pi/3 t)+", round(fit$coef[3],2), "cos(pi/3 t)") plot(time, expression, pch=16, ylab="", xlab="", xlim=c(0.5,10.5), ylim=c(-0.5,4.7)) preds <- predict(fit, newdata=list(time=time_grid)) lines(time_grid, preds, col=col.fit) abline(h=fit$coef[1], col="lightgrey") title(xlab="time (t)", line=2) title(ylab="expression (y)", line=2) title(main=model.name, cex.main=0.8) dev.off() ## ---- design19, eval=TRUE, echo=TRUE, include=TRUE---------------------------- cat("Output for Figure", ploti) design fit <- lm(expression~sin(2*pi*time/cycle)+cos(2*pi*time/cycle)) fit ## ---- eval=TRUE, echo=FALSE--------------------------------------------------- expression <- expression + time/4 + rnorm(n=length(expression), sd=0.2) ploti <- ploti+1 ## ---- design20, eval=FALSE, echo=FALSE---------------------------------------- # cat("Output for Figure", ploti) # model.matrix(~time+sinphase+cosphase) # fit <- lm(expression~time+sin(2*pi*time/cycle)+cos(2*pi*time/cycle)) # fit ## ---- include=FALSE, echo=FALSE----------------------------------------------- cat("Output for Figure", ploti) model.matrix(~time+sinphase+cosphase) fit <- lm(expression~time+sin(2*pi*time/cycle)+cos(2*pi*time/cycle)) fit plot.name <- paste0("fig", ploti, ".png") png(plot.name) model.name <- paste0("E(y)=", round(fit$coef[1],2), "+", round(fit$coef[2],2), "t+", round(fit$coef[3],2), "sin(pi/3 t)+", round(fit$coef[4],2), "cos(pi/3 t)") plot(time, expression, pch=16, ylab="", xlab="", xlim=c(0.5,10.5), ylim=c(0.5,6.8)) preds <- predict(fit, newdata=list(time=time_grid)) lines(time_grid, preds, col=col.fit) abline(fit$coef[1], fit$coef[2], col="lightgrey") title(xlab="time (t)", line=2) title(ylab="expression (y)", line=2) title(main=model.name, cex.main=0.8) dev.off() ## ---- design20, eval=TRUE, echo=TRUE, include=TRUE---------------------------- cat("Output for Figure", ploti) model.matrix(~time+sinphase+cosphase) fit <- lm(expression~time+sin(2*pi*time/cycle)+cos(2*pi*time/cycle)) fit ## ----------------------------------------------------------------------------- sessionInfo() ## ----------------------------------------------------------------------------- contrasts <- cbind( AvsC=c(1,0,-1,0), BvsC=c(0,1,-1,0), ABvsCD=c(0.5,0.5,-0.5,-0.5)) rownames(contrasts) <- LETTERS[1:4] contrasts ## ---- eval=FALSE-------------------------------------------------------------- # group <- as.factor(c(1,1,1,2,2,2)) # design <- model.matrix(~0+group) # contrasts <- makeContrasts(group1-group2, levels=colnames(design)) # v <- voom(counts, design) # fit <- lmFit(v, design) # fit <- contrasts.fit(fit, contrasts) # fit <- eBayes(fit) # topTable(fit)