limma Normalization question
2
0
Entering edit mode
@gordon-smyth
Last seen 2 hours ago
WEHI, Melbourne, Australia
> 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 ADD COMMENT 0 Entering edit mode @cecilia-mcgregor-1508 Last seen 8.1 years ago Dr Smyth I have some control spots that are not differentially expressed on the array. I tried to use the modifyWeights as on page 21 of the Users Guide (6 Oct 2005). However I realized that it doesn't work, cause I don't use weights in my read.maimages. I used a image analysis program, not recognized by limma, so used the following to read the data. RG <- read.maimages(targets$FileName, columns=list(R="Rf",G="Gf",Rb="Rb",Gb="Gb")) I tried to add "phony" column for weights (all values=1) in the intensity files, so I can modify it, but have not been able to get it to work (yet). Cecilia ----- Original Message ----- From: "Gordon K Smyth" <smyth@wehi.edu.au> To: "Cecilia McGregor" <cmcgre1 at="" lsu.edu=""> Subject: [BioC] limma Normalization question Date: Tue, 22 Nov 2005 07:54:59 +1100 (EST) > > > > 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 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 ADD COMMENT 0 Entering edit mode @gordon-smyth Last seen 2 hours ago WEHI, Melbourne, Australia >[BioC] limma Normalization question >Cecilia McGregor cmcgre1 at lsu.edu >Tue Nov 22 00:54:53 CET 2005 > >Dr Smyth > >I have some control spots that are not differentially expressed on the array. >I tried to use the modifyWeights as on page 21 of the Users Guide (6 Oct >2005). >However I realized that it doesn't work, cause I don't use weights in my >read.maimages. I used a image analysis program, not recognized by limma, >so used >the following to read the data. > >RG <- read.maimages(targets$FileName, >columns=list(R="Rf",G="Gf",Rb="Rb",Gb="Gb")) > >I tried to add "phony" column for weights (all values=1) in the intensity >files, >so I can modify it, but have not been able to get it to work (yet). > >Cecilia Dear Cecilia, There's no need to change the data files (a last resort!), you could just use for example RG\$weights <- array(1, dim(RG)) Best wishes Gordon