limma Normalization question
2
0
Entering edit mode
@gordon-smyth
Last seen 4 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
Genetics Normalization limma GeneSpring limmaGUI Genetics Normalization limma GeneSpring • 1.2k views
ADD COMMENT
0
Entering edit mode
@cecilia-mcgregor-1508
Last seen 9.6 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 4 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
ADD COMMENT

Login before adding your answer.

Traffic: 955 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6