limma_normalisation.R

library(limma)
Cy5 <- "F635 Median"
Cy5b <- "B635 Mean"
targets <- readTargets("Grade2_targets.csv")
f <- factor(targets$Condition, levels = unique(targets$Condition))
design <- model.matrix(~0 + f)
cont.matrix <- makeContrasts(Grade2vsControl=Grade2-Control, levels=design)
colnames(design) <- levels(f)

RG <- read.maimages( targets$FileName,
                     source="genepix.custom", green.only=TRUE, 
                     columns=list(G=Cy5,Gb=Cy5b))
RG.bc=backgroundCorrect(RG, method="normexp", offset=15)
controls <- grep("CONTROL", RG.bc$genes[,"ID"])
other.genes <- grep("empty|Empty|CONTROL", invert=TRUE, RG.bc$genes[,"ID"])
RG.nba <- normalizeBetweenArrays(RG.bc, method="quantile")
RG.final <- RG.nba[other.genes, ]
rownames(RG.final) <- RG.final$genes$ID
RG.final <- RG.final[order(rownames(RG.final)), ]

corfit <- duplicateCorrelation(RG.final, design, ndups=2, spacing=1)
fit <- lmFit(RG.final, design, ndups=2, correlation=corfit$consensus)
fit2 <- contrasts.fit(fit, cont.matrix)


fit2 <- eBayes(fit2)
#topTable(fit2, adjust="BH")