I wondered if anyone could comment on what is generally considered the most appropriate method of calculating fold change values from Affy data. I have a data set from a test vs. control experiment (n=3 in each group) performed on the MOE430_2 array. I used the Affy package to read in the CEL file data and GCRMA to normalize. Finally I used limma to analyze (lmFit and eBayes). Now the FC coefficients I get back from limma appear to be log2(2^(mean(test replicates in log2))/(2^mean(control replicates in log2)). For example:

> exprs(eset)[18731,1:6]
 TA2.CEL  TA5.CEL  TA7.CEL  TA1.CEL  TA8.CEL  TA9.CEL
3.294215 4.851795 4.403851 4.782934 7.716893 9.284909

> fit$coefficients[18731,1:2] #returns the mean of the above log2 values
  Fed_MT   Fed_WT
4.183287 7.261578

> fit2$coefficients[18731,1]
[1] -3.078291

> -1/2^fit2$coefficients[18731,1]
[1] -8.446136

So we have a probe that is reporting a -8.4 fold downregulation.

The problem is that SAM does not calculate the fold change values in the same manner. It appears to take the average of the unlogged data and then use those values to calculate fold changes. Thus:

> 2^exprs(eset)[18731,1:6]
   TA2.CEL    TA5.CEL    TA7.CEL    TA1.CEL    TA8.CEL    TA9.CEL
  9.809738  28.875927  21.168555  27.530016 210.385700 623.786694

> tmp <- c(mean(2^exprs(eset)[18731,1:3]),mean(2^exprs(eset)[18731,4:6]))
> tmp
[1]  19.95141 287.23414

> tmp[1]/tmp[2]
[1] 0.06946043

> -1/(tmp[1]/tmp[2])
[1] -14.39669

So now we have -14.4 fold downregulation.

The example given is of course an extreme, but using the limma method there are 282 probes with a fold change >2, while using the sam method there are only 159, with 141 probes in common. The probes that were called uniquely by limma were almost all downregulated (-2 to -5 fold).

MY QUESTION IS WHICH METHOD IS THE CORRECT WAY?

Thanks,

Peter

Peter White, Ph.D.
Technical Director
Functional Genomics Core
Department of Genetics
University of Pennsylvania
570 Clinical Research Building
415 Curie Boulevard
Philadelphia, PA 19104-6145
Tel: +1 (215) 898-0773
Fax: +1 (215) 573-2326
E-mail: pwhite at mail.med.upenn.edu
http://www.med.upenn.edu/pdc/cores_fgc.html
http://www.med.upenn.edu/kaestnerlab/
http://www.betacell.org/microarrays
http://www.cbil.upenn.edu/EPConDB/

Dear Peter,

Since I'm the limma author, it will come as no surprise to you that I view the limma fold change (diff of log-intensities) as the more natural, consistent and useful.

I don't know why SAM computes fold change as it does. It seems to me somewhat inconsistent, because SAM uses the limma log-fold-change as numerator in its own test statistic for computing statistical significant. So SAM reports a different fold change measure to that which it uses internally for determining statistical significance.

Best wishes
Gordon

BTW, here's an amusing example of the potential contradictions introduced by the SAM fold change definition. In this example, SAM finds the gene to be significantly up-regulated but the reported fold change is in the opposite direction!

 > y <- matrix(rnorm(100*20),100,20)
 > y[1,] <- rep(c(-2,2),each=10); y[1,1] <- y[1,1]-0.1
 > rownames(y) <- paste("Gene",1:100)
 > y[1,1] <- 6
 > out <- sam(y,cl)
 > summary(out,0.4)
SAM Analysis for the Two-Class Unpaired Case Assuming Unequal
Variances

s0 = 0.8  (The 100 % quantile of the s values.)

  Number of permutations: 100

  MEAN number of falsely called variables is computed.

Delta: 0.4
cutlow: -Inf
cutup: 2
p0: 0.866
Identified Genes: 1
Falsely Called Genes: 0
FDR: 0

Identified Genes (using Delta = 0.4):

   Row d.value stdev  rawp q.value R.fold   Name
1   1       2   0.8 1e-04 0.00866  0.604 Gene 1

Note that Gene1 is found to be significantly up-regulated (as indicated by a positive d.value=regularized t) but the reported fold change is less than one, i.e., indicating down-regulation.