Hi Edwin, sorry for the delay -- had to leave early on friday to enjoy the Swiss national holiday... :-) W.r.t. normalization, I am NOT a statistician and there are people on this list who know much better. I am rather following the general recommendation to normalize only as much as necessary and with the least aggressive method that will remove the issues from your diagnostic plots. We are using vsn for that and are so happy with the results that I have not looked into normalization for quite some time and felt no need to use any of loess, quantile, or control spot-based methods but I'm pretty sure that at least for vsn and loess you could fit a model on a subset of spots and predict for all spots, and that subset could be spike-ins, invariants, house-keeping genes or whatever you believe in... I also seem to recall, that for control spots the recommendation was to have a dilution series placed (i.e., spotted) on the array to avoid bias from manually pipetting the spike-ins... Sorry but I don't know the Yang and Thorne paper and cannot comment on that. If you tell your colleagues there's too much variation, how do they like that answer? :-) Seriously, I've gone through that, too, and one result is a function for fancy lattice dotplots showing individual "expression levels" according to another general recommendation: not to look at summary statistics if you can look at the real data... and with only three replicates I think we can do that... If you are interested in having a look at this function let me know and I'll try to convince the relevant people to share it. Cheers, - axel Axel Klenk Research Informatician Actelion Pharmaceuticals Ltd / Gewerbestrasse 16 / CH-4123 Allschwil / Switzerland "Edwin Groot" <edwin.groot at="" biol="" ogie.uni-freiburg="" to="" .de=""> "Bioconductor List" Sent by: <bioconductor at="" stat.math.ethz.ch=""> bioconductor-boun cc ces at stat.math.eth Axel.Klenk at actelion.com z.ch Subject Re: [BioC] LIMMA: Why does eBayes expression differ from observed? 31.07.2009 17:26 Hello Axel, Thanks for looking over my code. Yes, the arrays are Agilent, and the Feature Extraction program uses either Linear or LOWESS methods for normalization. Are you suggesting that I use control transcripts for normalization, when there are large changes in the transcriptome? The hybridizations include Agilent SpikeIn controls, but I do not know how to implement this in LIMMA. What do you think of only normalizing between arrays, using a Quantile method, as illustrated in Yang and Thorne (2003)? Do you have any idea of how to tell topTable() to output AveExpr and LogFC for only the contrast specified? This point is bothering my collegue, who wanted to know why certain transcripts were not at the top of the list. He wanted to see the replicate intensities. My short answer to this quibble, is that the expression was too variable, of course! Regards, Edwin --- On Fri, 31 Jul 2009 15:50:21 +0200 Axel.Klenk at Actelion.Com wrote: > > Dear Edwin, > > as far as I understand what you're trying to do, there are two > problems: > > First, according to ?topTable, AveExpr is the average expression over > all > (18) > channels whereas logFC is the fold change for the three mut1 vs, WT > only > and > you would need to also use the average fold change over all > conditions to > backcalculate the mean of your original R and G. > > Second, I think the formula you should use for that is the one used > in > RG.MA(): > R <- 2^(A + M/2) > G <- 2^(A - M/2) > and not the one you're using. > > Last, I'm guessing from your original column names that these are > Agilent > arrays > and AFAIK [gr]ProcessedSignal are loess-normalized by default (among > others) > which may be inappropriate in a common reference design if your > common > reference > is very different from your second channel as you have indicated -- > if you > want to avoid > that you should start from raw data using source = "agilent" in > limma's > read.maimages() > and then apply the preprocessing of your choice. > > Hope this helps, > > - axel > > > Axel Klenk > Research Informatician > Actelion Pharmaceuticals Ltd / Gewerbestrasse 16 / CH-4123 Allschwil > / > Switzerland > > > > > > "Edwin Groot" > > <edwin.groot at="" biol=""> > ogie.uni-freiburg > To > .de> "Bioconductor List" > > Sent by: > <bioconductor at="" stat.math.ethz.ch=""> > bioconductor-boun > cc > ces at stat.math.eth > > z.ch > Subject > Re: [BioC] LIMMA: Why does > eBayes > expression differ from > observed? > 31.07.2009 11:51 > > > > > > > > > > > > > > > > On Thu, 30 Jul 2009 14:57:43 +0200 > Axel.Klenk at Actelion.Com wrote: > > > > Dear Edwin, > > > > first of all, the code you have posted cannot work because it never > > assigns > > the results > > of your computations to any variables. In order to help you, we > need > > to > > know at least > > 1) what code you have really run, > > 2) your refdesign matrix, > > 3) the values from topTable() and RG.MA() for at least one example > > probe, > > 4) HOW you backcalculated R, G, and FC from topTable(), and > > 5) the obligatory output of sessionInfo() although this doesn't > > really > > sound like a version issue. :-) > > > > --- > Here is the code, with minimal commenting: > > refdesign <- modelMatrix(targets, ref="WT") > > refdesign > cf1 mut1 mut1_cf1 > [1,] 0 1 0 > [2,] 1 0 0 > [3,] 0 0 1 > [4,] 0 1 0 > [5,] 1 0 0 > [6,] 0 0 1 > [7,] 0 1 0 > [8,] 1 0 0 > [9,] 0 0 1 > > RG <- read.maimages(targets$FileName, columns = list(G = > "gProcessedSignal", Gb = "gBGMeanSignal", R = "rProcessedSignal", Rb > = > "rBGMeanSignal"), annotation= c("FeatureNum", "Row", "Col", > "ProbeUID", > "ControlType", "ProbeName", "Description", "GeneName", > "SystematicName")) > > i <- RG$genes$ControlType==0 > > RGnc <- RG[i,] > #RGnc already contains BG-subtracted and normalized intensities > > MA <- MA.RG(RGnc, bc.method="none") > #Intensity spread among arrays dissimilar > > MAAq <- normalizeBetweenArrays(MA, method="Aquantile") > > fit <- lmFit(MAAq, refdesign) > > > fitWT <- eBayes(fit) > > RGAq <- RG.MA(MAAq) > > options(digits=2) > > out.mut1 <- topTable(fitWT, coef="mut1", > genelist=cbind(fit$genes$GeneName,RGnc$G[,1],RGnc$R[,1],RGAq$G[,1:9],R Gnc$G[,1],RGAq$R[,1],RGAq$R[,4],RGAq$R[,7],RGnc$R[,1])) > > #Average WT (green) intensity > > out.mut1[,2] <- (2*2^out.mut1$AveExpr)/(2^out.mut1$logFC+1) > #Average mut1 (red) intensity > > out.mut1[,3] <- (2*2^out.mut1$AveExpr)/(2^-out.mut1$logFC+1) > > out.mut1[,2:17] <- as.numeric(as.matrix(out.mut1[,2:17])) > #Calculated mean WT intensity > > for (g in 1:10) out.mut1[g,13] <- > mean(as.numeric(as.matrix(out.mut1[g,4:12]))) > #Calculated mean mut1 intensity > > for (g in 1:10) out.mut1[g,17] <- > mean(as.numeric(as.matrix(out.mut1[g,14:16]))) > #Fold Change > > out.mut1$logFC <- 2^out.mut1$logFC > > colnames(out.mut1)[1:18] <- > c("GeneName","AvWT","Avmut1","WT1","WT2","WT3","WT4","WT5","WT6","WT7" , > "WT8","WT9","CalcWT","mut11","mut12","mut13","Calcmut1","FoldChange") > > out.mut1[9,] > GeneName AvWT Avmut1 WT1 WT2 WT3 WT4 WT5 WT6 WT7 WT8 WT9 > CalcWT > 10943 AT1G65370 31 298 64 57 37 56 47 59 69 64 69 > 58 > mut11 > 10943 423 > mut12 mut13 Calcmut1 FoldChange AveExpr t P.Value adj.P.Val > B > 10943 683 727 611 9.5 7.4 12 1.2e-06 0.0055 > 4.7 > --- > > Pardon the hack with mean(). It did not like the factor levels in the > data frame. > The WT and mut1 replicates, "observed RG", from the MAAq object are > in > WT1 to Wt9 and mut11 to mut13. > In the output table out.mut1, the expression calculated from topTable > (AvWT and Avmut1) is about half that of the expression calculated > from > the MAAq object (CalcWT and Calcmut1). The question is why are they > so > different? > The FC calculated from the MAAq object (10.5) also differs from the > FC > from topTable (9.5). > > Thanks in advance for your help, > Edwin > --- > > > AFAIK, eBayes() will affect the t, p, and B statistics computed > from > > your M > > values but not > > the AveExpr and logFC. > > > > RG.MA() backcalculates background-adjusted and normalized R and G > > values, > > so I am > > not sure what you mean by "observed" RG -- you don't get back the > > original > > R, Rb, G, Gb > > from an MAList via RG.MA(). > > > > Cheers, > > > > - axel > > > > > > Axel Klenk > > Research Informatician > > Actelion Pharmaceuticals Ltd / Gewerbestrasse 16 / CH-4123 > Allschwil > > / > > Switzerland > > > > > > > > > > > > > "Edwin Groot" > > > > <edwin.groot at="" biol=""> > > > ogie.uni-freiburg > > To > > .de> > bioconductor at stat.math.ethz.ch > > > > Sent by: > > cc > > bioconductor-boun > > > > ces at stat.math.eth > > Subject > > z.ch [BioC] LIMMA: Why does > eBayes > > > > expression differ from > > observed? > > > > > > > 30.07.2009 12:38 > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > I am getting odd results from a common reference design analysis of > > two-colour data in LIMMA. Previously I had analyzed only > > simple-comparison designs. Hopefully you can help restore > credibility > > of Bioconductor to my supervisor. > > Why does the AveExpr and logFC reported in topTable() differ from > the > > replicates in my MA object? > > > > The analysis is a textbook example of comparing a series of mutants > > to > > the wild type. Wild type is always green. The summary is as > follows: > > > > lmfit(MA, refdesign) > > eBayes(fit) > > topTable(fit, coef="mut1") > > #Regenerate the RG from MA > > RG.MA(MA) > > > > Backcalculating the topTable AveExpr and logFC to green, red and FC > > gives the following expressions as an example: > > WT: 31 mut1: 298 FC: 9.5 > > Compare that to the average of 9 WT and 3 mut1 in the RG.MA(MA) > list: > > WT: 58 mut1: 611 FC: 10.5 > > > > A survey of other probes gives an over and underestimate of the > > observed RG from 1.5 to 5 times. > > What is the explanation for that? > > What should I troubleshoot, besides looking at the usual BG and FG > > distributions and MA plots? > > > > Regards, > > Edwin > > --- > > Dr. Edwin Groot, postdoctoral associate > > AG Laux > > Institut fuer Biologie III > > Schaenzlestr. 1 > > 79104 Freiburg, Deutschland > > +49 761-2032945 > > > > _______________________________________________ > > Bioconductor mailing list > > Bioconductor at stat.math.ethz.ch > > https://stat.ethz.ch/mailman/listinfo/bioconductor > > Search the archives: > > http://news.gmane.org/gmane.science.biology.informatics.conductor > > _______________________________________________ > Bioconductor mailing list > Bioconductor at stat.math.ethz.ch > https://stat.ethz.ch/mailman/listinfo/bioconductor > Search the archives: > http://news.gmane.org/gmane.science.biology.informatics.conductor > > > > > The information of this email and in any file transmitted with it is > strictly confidential and may be legally privileged. > It is intended solely for the addressee. If you are not the intended > recipient, any copying, distribution or any other use of this email > is prohibited and may be unlawful. In such case, you should please > notify the sender immediately and destroy this email. > The content of this email is not legally binding unless confirmed by > letter. > Any views expressed in this message are those of the individual > sender, except where the message states otherwise and the sender is > authorised to state them to be the views of the sender's company. For > further information about Actelion please see our website at > http://www.actelion.com > Dr. Edwin Groot, postdoctoral associate AG Laux Institut fuer Biologie III Schaenzlestr. 1 79104 Freiburg, Deutschland +49 761-2032945 _______________________________________________ Bioconductor mailing list Bioconductor at stat.math.ethz.ch https://stat.ethz.ch/mailman/listinfo/bioconductor Search the archives: http://news.gmane.org/gmane.science.biology.informatics.conductor The information of this email and in any file transmitted with it is strictly confidential and may be legally privileged. It is intended solely for the addressee. If you are not the intended recipient, any copying, distribution or any other use of this email is prohibited and may be unlawful. In such case, you should please notify the sender immediately and destroy this email. The content of this email is not legally binding unless confirmed by letter. Any views expressed in this message are those of the individual sender, except where the message states otherwise and the sender is authorised to state them to be the views of the sender's company. For further information about Actelion please see our website at http://www.actelion.com