## R Code for assignment 3 library(limma) targets <- readTargets("SwirlSample.txt") swirlRG <- read.maimages(targets$Names, source = "spot") swirlRG$genes <- readGAL("fish.gal") swirlRG$printer <- getLayout(swirlRG$genes) names(swirlRG) print(swirlRG$targets) # Examine the information stored here head(swirlRG$genes) # Check the first few rows only # Proceed similarly, if required, for other named elements # Exercise 1 names(swirlRG) # Exercise 2 -- Spatial Plots slide.image(swirlRG$R[, 1], layout = swirlRG$printer) # Example 3 - MA Plots & Normalization # No normalization swirlMAn <- normalizeWithinArrays(swirlRG, method = "none") plotPrintTipLoess(swirlMAn) rm(swirlMAn) # Apply loess normalization, and check the MA plots: \begin{verbcode} swirlMAw <- normalizeWithinArrays(swirlRG) plotPrintTipLoess(swirlMAw) # Check whether normalization seems required between arrays: boxplot(swirlMAw$M ~ col(swirlMAw$M), names = colnames(swirlMAw$M)) # normalize between arrays swirlMA <- normalizeBetweenArrays(swirlMAw) rm(swirlMAw) boxplot(swirlMA$M ~ col(swirlMA$M), names = colnames(swirlMA$M)) # Exercise 4 -- Checks for Differential Expression # First, fit a simple statistical model that allows for the possible # effect of the dye swap, and that can be used as the basis for # checks for differential expression. The calculations require a # design vector that has -1 wherever there is a dye swap: design <- c(-1, 1, -1, 1) swirlfit <- lmFit(swirlMA, design) # From the fit object, calculate empirical Bayes moderated $t$-statistics. efit <- eBayes(swirlfit) ## You might like to examine a qq-plot qqt(efit$t, df = efit$df.prior + efit$df.residual, pch = 16, cex = 0.2) # Print out a table that shows the top 25 differentially expressed # genes, in order of the values of the moderated$t$-statistic: options(digits = 3) topTable(efit, number = 25, adjust="fdr") # Repeat, but now omitting: adjust="fdr" # Repeat, but now set the prior estimate of the proportion to be 0.1. # (For this, specify: eBayes(swirlfit, proportion=0.1, adjust="fdr") ## Optional extra -- not part of the assignment plot(efit$coef, efit$lods, pch = 16, cex = 0.2, xlab = "log(fold change)", ylab = "log(odds)") ord <- order(efit$lods, decreasing = TRUE) top8 <- ord[1:8] text(efit$coef[top8], efit$lods[top8], labels = swirlRG\$genes[top8, "Name"], cex = 0.8, col = "blue")