## 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")