Microarray: 如何讓 limma 使用 Tab delimited text table 計算 p 值

在統計上用來區分兩群族有無差異,常用的方法是 t-test 之類的方法, t-test 的原理是用觀察兩群樣本的分布來作判定是否有差異,而為了顯現族群的分布,適當的取樣是必需的(至少上百個sample), 但是microarray 實驗,通常沒有多少人可以拿出三五百萬來作一個實驗, 所以 limma 就被發展出來, 用以補足Microarray 傳統上 t-test 統計方法上的不足, limma 是一個線性的模型的方法概念上有點像ANOVA,不過原理我就不詳述了,而在應用上,目前我也只知道 R 上有 package, limma, 而且還出了 GUI 的界面,讓不會寫R 的使用者可以用選單的方法, raw data (.gpr or .cel) 餵進去, p-value 跑出來, 算是很成熟的演算法了.  不過也許太過 user freindly 了,反而讓我很傷腦筋. 最近我接了一個任務,要作一組分析,客戶指定要用 limma 來挑基因,不過,不是從 raw data 開始,我拿到的是一組已經處理過的數據(tab delimited text table), 可是 limma 可沒支援這種半生不熟的data 格式, 所以我花了一點時間,想了一個方法騙過 limma, 從中間插入.如果你也有相同的問題,希望我的方法,對你有幫助:

load limma 的 package
>library(limma)

接下來,讀取 data file,格式像這樣

            S1  S2  S3   ...    S7
gene1    1    2    2    ...    1
gene3    2    3    1    ...    0

以 Tab 分開. s1-s7是sample name

 >raw <- read.table (file="c:/yourdata.txt",sep="\t",row.names=1,header=T)

然後把他假裝成 marray 產生的物件 ExpresssionSet, limma就可以把資料吃進去了!
> eset<-new("ExpressionSet",exprs=as.matrix(raw))

接下來就一樣了,先設定實驗設計模型, 然後套入linear model 及 contrast model
>design <- model.matrix (~-1+factor(c(1,1,1,1,2,2,2,2))
>colnames(design) <- c("test","control")
> fit <- lmFit(eset,design=design)
>contrast.matrix <-makeContrasts(test-control, levels=design)
>fit2<-contrasts.fit(fit,contrast.matrix)
然後用 empirical Bayes smoothing 的方法,求 p-values
>P <- eBayes(fit2)

這個就是你要的啦!
>P$p.value
##ReadMore##