[BioC] limma Normalization question

[Cecilia McGregor]

Hi Everyone

I've described my experiment in an earlier message, on Nov 16. After running the commands
(see end of message for commands) I looked at the plotMA(fit2)plot. I attached the plot.
It seems to me that at high A-values, the plot is going more in the direction of positive
M-values. I tested a few genes in the circled area (for high A-values)with Q-RT-PCR and it
confirms that these genes should have M-values around 0. The fact that more genes
are down-regulated, than up-regulated is expected from our knowledge of the experiment.
In the limma Users Guide (6 Oct 2005) p21, it says that loess does not assume equal number
of genes up- and down-regulated, so this should not be a problem. I have about 2500 unique
genes on the array.However it does say that most genes need to be not differentially
expressed. I estimate about 60% of genes are not differentially expressed. In a message
on this board (Sept 3, 2003), Gordon Smyth suggests the use of the following command if
a large number of genes are differentially expressed.

normalizeWithinArrays(RG, method="robustspline", robust="MM")

However I get 18 of the following warning messages:
Warning messages:
1: rlm failed to converge in 20 steps in: rlm.default(x, y, weights = w, method = method)

I believe Paul Boutros suggested a fix for this in a message on this board(Aug 11 2004),
but since I don't know much about the subject, I'm not very eager to change the limma
code (I think that is what he suggest).

How can I fix this problem I'm having at high A-Values? I get the correct results at low
A-values. I have analised this data also with limmaGUI (I just treated all replicates
as biological)and I get the similar results. I used several different within and between
normalization methods, but cannot get better results.

I even tried normalizing with Genespring, altough I know it is not apropriate for loop
designs. With Genesspring I can get an acceptable MA plot if I use "Per array normalizations
to the 20th percentile". However I'm not very comfortable with this, since inspecting the
MA plots before normalization, I can see that I have an intensity dependent bias on all my
arrays.

Any suggestions would be much appreciated.

Code I Used: (I get the same results with the code suggested by Steen Krogsgaard on Nov 17)

library(limma)
setwd("C:\\Program Files\\R_JSM\\rw2011\\RFlimma")
targets <- readTargets()
show(targets)
RG <- read.maimages(targets$FileName,
columns=list(R="Rf",G="Gf",Rb="Rb",Gb="Gb"))
RG$genes <- readGAL()
spottypes <-readSpotTypes()
RG$genes$Status <-controlStatus(spottypes, RG)
MA <- normalizeWithinArrays(RG, method="loess")
MA <- normalizeBetweenArrays(MA, method="Aquantile")
design <-modelMatrix(targets,ref="S1")
cor <- duplicateCorrelation(MA,design,ndups=3,spacing=3120)
cor$consensus.correlation
fit <-lmFit(MA,design,ndups=3,spacing=3120,correlation=0.3128828)
cont.matrix <- makeContrasts(fvss=(F1+F2+F3-S2-S3)/3, levels=design)
fit2 <- contrasts.fit(fit, cont.matrix)
fit2 <- eBayes(fit2)
plotMA(fit2)

Thanks
Cecilia McGregor

PhD Student
Sweetpotato Breeding and Genetics Lab
JC Miller Hall room 236
Louisiana State University
Baton Rouge
LA, 70803
USA

Phone: (225) 578 2173