# Analysis of the LMS controls from the QC data set

library(limma)

# Load data
load("qc.Rdata")

# Normalise controls using an intensity based loess fit, with an increased span
MAlms <- normalizeWithinArrays(RGlms, method="loess", span=0.75)

# Remove arrays for which spike-in controls were not added
spikedIn <- c(1:7,10,15,21,50) 

# Estimate array weights
lms.array.wts <- arrayWeights(MAlms[,-spikedIn])

# Set up array numbers for plotting
array.nums <- rep("", 100)
array.nums[1] <- 1
array.nums[100] <- 100
for(i in 2:10) {
	array.nums[(i-1)*10] <- (i-1)*10
	}

# Plot array weights
barplot(lms.array.wts, xlab="Array", names.arg=array.nums, ylab="Weight", col="white", las=2)
abline(h=1, lwd=1, lty=2)

# Lowest and highest weights from the data set
round(lms.array.wts[c(19,91)],2)
# [1] 0.11 3.68

# MA plots for arrays with highest and lowest weights
par(mfrow=c(2,1))
i <- 91
plot(MAlms$A[status[indexControls]=="U10",-spikedIn][,i], MAlms$M[status[indexControls]=="U10",-spikedIn][,i], pch=1, xlab="A", ylab="M", main=paste("Array", i, ", w = ", round(lms.array.wts[i], 2), sep=" "), col=gray(0.8), xlim=c(6,16), ylim=c(-4,4))
points(MAlms$A[status[indexControls]=="U03",-spikedIn][,i], MAlms$M[status[indexControls]=="U03",-spikedIn][,i], pch=2, col=gray(0.6))
points(MAlms$A[status[indexControls]=="DR",-spikedIn][,i], MAlms$M[status[indexControls]=="DR",-spikedIn][,i], pch=3, col=gray(0.4))
points(MAlms$A[status[indexControls]=="D03",-spikedIn][,i], MAlms$M[status[indexControls]=="D03",-spikedIn][,i], pch=4, col=gray(0.2))
points(MAlms$A[status[indexControls]=="D10",-spikedIn][,i], MAlms$M[status[indexControls]=="D10",-spikedIn][,i], pch=5, col=gray(0))
abline(h=log(10,2), col=gray(0.8), lty=2, lwd=1)
abline(h=log(3,2), col=gray(0.6), lty=2, lwd=1)
abline(h=0, col=gray(0.4), lty=2, lwd=1)
abline(h=-log(3,2), col=gray(0.2), lty=2, lwd=1)
abline(h=-log(10,2), col=gray(0), lty=2, lwd=1)
legend(6, 4, legend=c("U10", "U03", "DR", "D03", "D10"), col=gray(c(0.8, 0.6, 0.4, 0.2, 0)), pch=seq(1:5), bg="white")
i <- 19
plot(MAlms$A[status[indexControls]=="U10",-spikedIn][,i], MAlms$M[status[indexControls]=="U10",-spikedIn][,i], pch=1, xlab="A", ylab="M", main=paste("Array", i, ", w = ", round(lms.array.wts[i], 2), sep=" "), col=gray(0.8), xlim=c(6,16), ylim=c(-6,4))
points(MAlms$A[status[indexControls]=="U03",-spikedIn][,i], MAlms$M[status[indexControls]=="U03",-spikedIn][,i], pch=2, col=gray(0.6))
points(MAlms$A[status[indexControls]=="DR",-spikedIn][,i], MAlms$M[status[indexControls]=="DR",-spikedIn][,i], pch=3, col=gray(0.4))
points(MAlms$A[status[indexControls]=="D03",-spikedIn][,i], MAlms$M[status[indexControls]=="D03",-spikedIn][,i], pch=4, col=gray(0.2))
points(MAlms$A[status[indexControls]=="D10",-spikedIn][,i], MAlms$M[status[indexControls]=="D10",-spikedIn][,i], pch=5, col=gray(0))
abline(h=log(10,2), col=gray(0.8), lty=2, lwd=1)
abline(h=log(3,2), col=gray(0.6), lty=2, lwd=1)
abline(h=0, col=gray(0.4), lty=2, lwd=1)
abline(h=-log(3,2), col=gray(0.2), lty=2, lwd=1)
abline(h=-log(10,2), col=gray(0), lty=2, lwd=1)
legend(6, 4, legend=c("U10", "U03", "DR", "D03", "D10"), col=gray(c(0.8, 0.6, 0.4, 0.2, 0)), pch=seq(1:5), bg="white")

# Linear model analysis
# Fit linear model using equal arrays weights
out.no.wts <- lmFit(MAlms[,-spikedIn])
out.no.wts <- eBayes(out.no.wts)

# Fit linear model with empirical array weights 
# Make a matrix of weights for lmFit
lms.array.wts <- matrix(rep(lms.array.wts, each=dim(MAlms[,-spikedIn])[1]), 
                         dim(MAlms[,-spikedIn])[1], dim(MAlms[,-spikedIn])[2])

out.array.wts <- lmFit(MAlms[,-spikedIn], weights=lms.array.wts)
out.array.wts <- eBayes(out.array.wts)

# Calculate ordinary t statistics
# standard errors
sigma.lms <- cbind(out.no.wts$sigma*out.no.wts$stdev.unscaled,
                   out.array.wts$sigma*out.array.wts$stdev.unscaled)

# t statistics
ord.tstat.lms <- cbind(out.no.wts$coef/sigma.lms[,1], 
                   out.array.wts$coef/sigma.lms[,2])
	
# Plot t statistics
names.spot.class <- rep(c("D03 ", "D10 ", "DR ", "U03 ", "U10 "), each=2)
weight.names <- c("Equal", "Array")

par(mfrow=c(1,1))
boxplot(ord.tstat.lms[,1]~as.factor(status[indexControls]), at=1:5-0.2,col="gray", boxwex=0.4,
	las=2, ylab="Ordinary t statistics", pch=".", ylim=c(-150, 125), medlwd=1) 
boxplot(ord.tstat.lms[,2]~as.factor(status[indexControls]),boxwex=0.4, 
        at=1:5+0.2,col="white",add=TRUE,yaxt="n",xaxt="n",medlwd=1, pch=".")
abline(h=0, col="black", lty=2, lwd=1)
legend(0.5, 100, legend=weight.names, fill=c("gray", "white"), cex=0.8)