[BioC] limma Normalization question

[Gordon K Smyth]

> Date: Sun, 20 Nov 2005 17:09:09 -0600
> From: "Cecilia McGregor" <cmcgre1 at lsu.edu>
> Subject: [BioC] limma Normalization question
> To: bioconductor at stat.math.ethz.ch
>
> 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.

This is really pushing the envelope.  I guess that loess can probably handle up to around 20%
differentially expressed, while robustspline can handle around 30%.  With 40% you probably need
some symmetry as well.

> 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)

This is probably not a problem, so doesn't need a fix.

Another alternative would be to increase the robustifying iterations used with "loess", say

> MA <- normalizeWithinArrays(RG, method="loess", iterations=10)

That is more robust than the default "loess" method.

Do you have any useful control probes on your arrays?

Best wishes
Gordon

> 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