R help - Mar 2010 - How to calculate fold change and its confidence interval

Dear all,

I am reading the book "statistic analysis of gene expression microarray 
data" by Terry Speed .
On page16-17, there is an example that I can not repeat.

The content can be found using below linkage:
http://books.google.com/books?id=MBzcRjez9ccC&pg=PA16&lpg=PA16&dq=confidence+interval+for+fold+change++statistic+analysis+of+gene+expression+microarray+data&source=bl&ots=pSfB1_-sb4&sig=zaUycGxgD__GMHZtaSnOwqLUnsk&hl=en&ei=k6iyS_DGGZLcM9CjpaEE&sa=X&oi=book_result&ct=result&resnum=4&ved=0CBIQ6AEwAw#v=onepage&q=&f=false


I want to know how to calculate the confidence interval for fold change 
of microarray data in R.
My questions in details are:
1) why the rowMeans() gives different results as simple function mean()
2) How to calculate the CI for the fold changes.


For example, we can get a matrix of 4 genes of 4 samples in in two groups:

 > tb <- 
matrix(c(860,406,284,46,42,31,29,64,348,164,115,19,36,44,18,85),nrow=4)
 > colnames(tb) <-
c("A","A","B","B")
 > tb
       A  A   B  B
[1,] 860 42 348 36
[2,] 406 31 164 44
[3,] 284 29 115 18
[4,]  46 64  19 85

The fold change can be calculated as
 > fc <- apply(tb[,1:2],1,mean) / apply(tb[,3:4],1,mean)
 > fc
[1] 2.348958 2.100962 2.353383 1.057692

And the paired fold change should be:
 > fc.perSample <- cbind(tb[,1]/tb[,3], tb[,2]/tb[,4])
 > fc.perSample
         [,1]      [,2]
[1,] 2.471264 1.1666667
[2,] 2.475610 0.7045455
[3,] 2.469565 1.6111111
[4,] 2.421053 0.7529412

 > rowMeans(fc.perSample)
[1] 1.818966 1.590078 2.040338 1.586997


But if you manually check every gene one by one, totally different 
results were got. for example:
 > mean(2.471264, 1.1666667)
[1] 2.471264


Many thanks for any help.

Xinan

Hi,

Firstly, you should bring your bioinformatics related questions over
to the bioconductor mailing list. Please sign up and post that portion
there, see:

http://www.bioconductor.org/docs/mailList.html

But, getting to your last question:
>> fc.perSample <- cbind(tb[,1]/tb[,3], tb[,2]/tb[,4])
>> fc.perSample
> ? ? ? ?[,1] ? ? ?[,2]
> [1,] 2.471264 1.1666667
> [2,] 2.475610 0.7045455
> [3,] 2.469565 1.6111111
> [4,] 2.421053 0.7529412
>
>> rowMeans(fc.perSample)
> [1] 1.818966 1.590078 2.040338 1.586997
>
> But if you manually check every gene one by one, totally different results
> were got. for example:
>> mean(2.471264, 1.1666667)
> [1] 2.471264
In short: you're doing it wrong.

You are passing 1.166 as the second parameter to the mean function.
See its help page to figure out what the params are: ?mean

You have to pass it in a *vector* of numbers to take the mean from:

R> mean(c(2.471264, 1.1666667))
[1] 1.818965

-steve

-- 
Steve Lianoglou
Graduate Student: Computational Systems Biology
 | Memorial Sloan-Kettering Cancer Center
 | Weill Medical College of Cornell University
Contact Info: http://cbio.mskcc.org/~lianos/contact