library(limma)
library(vsn)

# DATA PREPROCESSING
# Download RAW data from GEO experiment
setInternet2(use=FALSE)
download.file('http://www.ncbi.nlm.nih.gov/geo/download/?acc=GSE34571&format=file', 'GSE34571_RAW.tar', mode='wb')
untar('GSE34571_RAW.tar', exdir='gpr')
setwd('gpr')

# Set function to assign weights
f <- function(x) {ok.flags <- x$Flags > (-99)
                  ok.snr <- x[,'SNR 532'] > 3
                  as.numeric(ok.snr & ok.flags)
}

# Read single-channel .gpr files as double-channel data and perform bg correction  
RG <- read.maimages(dir(pattern='*.gpr.gz'), source="genepix",
                    columns=list(R="F532 Mean",G="F532 Mean",Rb="B532",Gb="B532"), 
                    wt.fun=f)
RG.b = backgroundCorrect(RG=RG, method='normexp')

# Delete fake second channel
rownames(RG.b$R) <- RG.b$genes$Name
RG.b$G <- NULL
RG.b$Gb <- NULL

# Specify contrast and design matrices
design <- model.matrix(~ 0+factor(c(1,1,1,1,2,2,2,2)))
colnames(design) <- c("group1", "group2")
contrast.matrix <- makeContrasts(group1-group2, levels=design)

# Peform normalisation and stats
norm.vsn <- normalizeVSN(RG.b$R)
fit.vsn <- lmFit(norm.vsn, design)
fit2.vsn <- contrasts.fit(fit.vsn, contrast.matrix)
fit3.vsn <- eBayes(fit2.vsn)
limma::volcanoplot(fit3.vsn)


RG$printer$ngrid.r = 1
imageplot(log2(RG.b$R[,2]), RG$printer)

# DIAGNOSTIC PLOTS

# maplot.1
par(mfrow=c(3,2))
limma::plotMA(log2(RG$R[,c(1:4)]))
limma::plotMA(log2(RG$R[,c(5:8)]))

limma::plotMA(log2(RG.b$R[,c(1:4)]))
limma::plotMA(log2(RG.b$R[,c(5:8)]))

limma::plotMA(norm.vsn[,c(1:4)])
limma::plotMA(norm.vsn[,c(5:8)])

# maplot.2.final
par(mfrow=c(1,1))
limma::plotMA(fit3.vsn)

# maplot.3.filtered
par(mfrow=c(2,1))
limma::plotMA(log2(RG.b$R[((RG.b$weights[,1] == 1) & (RG.b$weights[,2] == 1)),c(1:2)]))
limma::plotMA(log2(RG.b$R[,c(1:2)]))

boxplot(RG$R)
boxplot(norm.vsn)

# plotdensities
plotDensities(RG, singlechannels=c(1:8))
plotDensities(RG.b, singlechannels=c(1:8))