[BioC] RE: RMA normalization

[w.huber at dkfz-heidelberg.de]

On Wed, Sep 15, 2004 at 10:20:42AM +0200, Matthew  Hannah wrote:

>> I've tried to get a discussion on this several times but have got very
>> few responses.
>>
>> I'm looking at some data where the treatment has a very BIG effect, but
>> I don't think this is unusual it's just that a lot of people don't
>> realise it or ignore it.

Hi Matthew,

Robert Gentleman and myself just talked about the issue that you brought
up and here are some of our points:

If by big effect you mean that a lot of genes change expression under
different conditions, then I think that many people are aware of the
issue, but that it has not been widely discussed. Frank Holstege and his
group in Utrecht have worked quite systematically on this, and I would
recommend looking at one of their recent papers:

Monitoring global messenger RNA changes in externally controlled
microarray experiments. van de Peppel J, Kemmeren P, van Bakel H, Radonjic
M, van Leenen D, Holstege FC. EMBO Rep. 2003 Apr;4(4):387-93.  PMID:
12671682

You seemed to indicate that having a large fraction of genes changing
being only a problem for RMA and GCRMA, but our understanding is that it
is a problem for all "intrinsic" methods of normalizing, i.e. those do not
use external spike-in controls. We believe that the problem is largely to
do with normalization, and not with computing expression estimates.

As long as microarray experiments are not trying to measure absolute
molecule abundances, but rather just "relative expression", then we think
that the problem of interpreting situations in which most genes change
will always remain hard.

On the other hand side, if you use an "extrinsic" method, you need to
decide whether you want to measure number of molecules per cell, or per
total RNA, or per what ... so that's a conceptual issue that needs to be
worked out.

This is also known as a research *opportunity*.

>> If we take the average Pearson correlation of treated versus untreated
>> as a crude indication of the number of changes (is this valid?) then in
>> our experiments this is 0.96. Comparing 25 treated versus 25 untreated
>> replicates (GCRMA, LIMMA gene-wise fdr corrected p<0.001) we get c.30%
>> of transcripts on the chip changing!

You are too imprecise in your description for us to comment on whether the
method is valid, but why don't you use good old-fashioned statistics: for
each gene, calculate a p-value from e.g. t-test or an appropriate linear
model generalization, and look at the histogram of p-values. The empirical
p-values at the right end of the histogram should be approximately
uniformly distributed, and the number of non-differentially expressed
genes can be estimated by 2 times the genes with p>0.5.

>> Looking at a couple of public datasets I don't think our treatment
>> effect (as indicated by the Pearson) is that unusual, it's just that we
>> have the statistical power to detect the changes. Also looking at the
>> changes, and considering the biology it seems reasonable to get these
>> changes.

It depends a lot on what the factors are. Robert has some collaborators
who use treatments that greatly change things, and others that use
treatments that are so specific that the number of changes genes is under
10. No global statement that can be made here. For the former, they need
to be warned that their experiment lies outside of the currently available
technology, and then you do the best you can. With the latter, we should
be able to do a good job, although one may never find the signal if
p-value corrections are applied in a naive fashion (but that is a
different story).

>> In the discussions on RMA/GCRMA there are 2 assumptions discussed
>> 1)few genes changing - obviously not
>> 2)equal # up and down - despite the huge amount of changes there are
>> only 20 more transcripts going up compared to down - so yes.

As I said above, I do not believe this is specific to RMA or GCRMA, but
rather a general problem for all normalization methods.

Also, these are sufficient conditions, but not necessary. I.e. if they are
fulfilled, (GC)RMA can be guaranteed to work, but if they aren't, the
results from (GC)RMA, or other normalization methods, may still be valid
to a sufficient degree!

>> I've also looked at a number of control genes and can't find any real
>> bias, in fact there is quite a bit of (random?) variation, so if you
>> normalised on a few of these then you may get strange results...

 I think the problem with using control genes is that there are
 typically few of them and they do not necessarily span the range of
 intensities so they provide a poor basis for normalization. Although,
 in the present case they may be better than the alternatives.

>> Finally, I will at some point try separate GCRMAs and then scaling. If
>> anyone has any scripts for mean, robust mean or median scaling a series
>> of separate exprs sets then I'd appreciate it.

I hope that this is a rather simple use of lapply or similar
(depending on how you have your exprSets stored).

  Regards,
    Robert, Wolfgang


+------------------------------------------------------------------------+
| Robert Gentleman phone : (617) 632-5250                                |
| Associate Professor fax: (617) 632-2444                                |
| Department of Biostatistics office: M1B20                              |
| Harvard School of Public Health email: rgentlem at jimmy.harvard.edu      |
+------------------------------------------------------------------------+