[BioC] LIMMA: Why does eBayes expression differ from observed?
Axel.Klenk at Actelion.Com
Axel.Klenk at Actelion.Com
Fri Jul 31 15:50:21 CEST 2009
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],RGnc$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
More information about the Bioconductor
mailing list