Limma

library(limma)

setwd("d:/agti/BioCat/swirl")

targets <- readTargets("SwirlSample.txt")

# source <- c( "genepix", "genepix.median", "imagene", "quantarray", "smd.old", "smd", "spot", "spot.close.open")

RG <- read.maimages(targets$FileName, source="spot")

RG$genes <- readGAL("fish.gal")

spottypes <- readSpotTypes("SpotTypes.txt") # information spot types

RG$genes$Status <- controlStatus(spottypes, RG)

RG$printer <- getLayout(RG$genes)


# Normalize Within
MAw <- normalizeWithinArrays(RG) # bacground correction method = "subtract"

# method <- c("none", "subtract", "half", "minimum", "movingmin", "edwards", "normexp")
RGb <- backgroundCorrect(RG, method="subtract")
MAw <- normalizeWithinArrays(RGb) # bacground correction method = "subtract"

# Normalize Between
MAwb <- normalizeBetweenArrays(MAw)

### Unnormalized Data
plotMA(RG[,2]) # Color-Coded M A Plot (for one slide) for the unnormalized data
plotPrintTipLoess(RG[,2]) # Print-Tip Group M A Plot (for one slide)
imageplot(log2(RG$R[,1]), RG$printer, low="white", high="red") # image array plot
boxplot(RG$R~col(RG$R))

### Normalized Data
plotMA(MAw[,2])
plotPrintTipLoess(MAw)
imageplot(MAw$M[,1], RG$printer, zlim=c(-3,3))
boxplot(MAw$M~col(MAw$M), names=colnames(MAw$M))
boxplot(MAwb$M~col(MAwb$M), names=colnames(MAwb$M))

### Linear Model
design <- c(-1, 1, -1, 1)
fit <- lmFit(MAwb, design)

# ordinary.t <- fit$coef / fit$stdev.unscaled / fit$sigma
fiteb <- eBayes(fit)

options(digits=3)

### Number of genes in table
ngene = c(10,30,50,100,nrow(fit2))
mgene = ngene[menu(ngene, graphics=F, title="Number of genes in table?")]

### Sort by
sortby = c("M","A","t","P","B")
msortby = sortby[menu(sortby, graphics=F, title="Sort by?")]

### Adjust method
adjust = c("bonferroni","holm","hochberg","hommel","fdr","none")
madjust = adjust[menu(adjust, graphics=F, title="Adjust method?")]
topTable(fiteb, coef=1, number=mgene, sort.by=msortby, adjust=madjust)

# MA Plot (using fitted M values)
plotMA(fiteb)
ord <- order(fiteb$lods, decreasing=T)
top <- ord[1:mgene]
text(fiteb$Amean[top], fiteb$coef[top], labels=fiteb$genes[top,"Name"], cex=0.8, col="blue")

# Log Odds Plot
volcanoplot(fiteb, highlight=8, names=fiteb$genes$Name)

# Quantile-Quantile t statistic plot
qqt(fiteb$t, df=fiteb$df.prior+fiteb$df.residual, pch=16, cex=0.2)
abline(0,1)