# Libraries library(ABC2017) library(ggplot2) theme_set(theme_minimal()) library(edgeR) # data data("rat")
rat$Design
# Trim EM <- subset(rat$Expression, rowSums(rat$Expression >= 3) >= 2) # Calculate normalization factors dge <- calcNormFactors(DGEList(EM), method="TMM") # Normalized values logEM <- cpm(dge, prior.count=1, log=TRUE)
plotDensities(logEM, group=rat$Design$Time)
# Perform MDS, but only save output mds <-plotMDS(dge, top=1000, gene.selection="pairwise", plot=FALSE) # Save as a data.frame P <- data.frame(rat$Design, mds$cmdscale.out)
qplot(data=P, x=X1, y=X2, color=protocol, shape=Time, label=rownames(P))
# Relevel to make "untreated" intercept rat$Design$Time <- relevel(factor(rat$Design$Time), "2 weeks") mod <- model.matrix(~protocol*Time, data=rat$Design) mod
# Perform MDS, but only save output ldf <-plotRLDF(logEM, design=mod, nprobes=1000, plot=FALSE, trend=TRUE, robust=TRUE) # Save as a data.frame P <- data.frame(rat$Design, ldf$training)
qplot(data=P, x=X1, y=X2, color=protocol, shape=Time, label=rownames(P))
# Calculate dispersion disp <- estimateDisp(dge, design=mod, robust=TRUE) # Fit models fit <- glmFit(disp, design=mod) # Perform the tests for the coefficients time_effect <- glmLRT(fit, coef="Time2 months") SNL_effect <- glmLRT(fit, coef="protocolL5 SNL") interaction_effect <- glmLRT(fit, coef="protocolL5 SNL:Time2 months")
plotBCV(disp)
nDE <- lapply(list(time=time_effect, SNL=SNL_effect, interaction=interaction_effect), decideTestsDGE) nDE <- sapply(nDE, summary) rownames(nDE) <- c("Down", "Stable", "Up") nDE
plotSmear(SNL_effect, de.tags=rownames(SNL_effect)[decideTestsDGE(SNL_effect) != 0])
topTags(SNL_effect)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.