差異分析流程示例與資料
基因芯片的差異表達(dá)分析主要有 構(gòu)建基因表達(dá)矩陣、構(gòu)建實(shí)驗(yàn)設(shè)計(jì)矩陣重归、構(gòu)建對(duì)比模型（對(duì)比矩陣）米愿、線性模型擬合、貝葉斯檢驗(yàn)和生成結(jié)果報(bào)表 六個(gè)關(guān)鍵步驟鼻吮。
下面是模擬的一個(gè)示例：
# Simulate gene expression data for 100 probes and 6 microarraexprSets
# MicroarraexprSet are in two groups
# First two probes are differentiallexprSet expressed in second group
# Std deviations varexprSet between genes with prior df=4
# 構(gòu)建模擬的表達(dá)矩陣育苟，實(shí)際處理時(shí)換成自己的表達(dá)矩陣即可
sd <- 0.3*sqrt(4/rchisq(100,df=4))
exprSet <- matrix(rnorm(100*6,sd=sd),100,6)
rownames(exprSet) <- paste("Gene",1:100)
colnames(exprSet) <- c(paste0("con-",1:3), paste0("G3-",1:3))
exprSet[1:2,4:6] <- exprSet[1:2,4:6] + 2
library(limma)
# 構(gòu)建實(shí)驗(yàn)設(shè)計(jì)矩陣
group_list = c(rep("con",3), rep("G3",3))
# 這里根據(jù)實(shí)際的情況設(shè)置（表型）分組，對(duì)應(yīng)表達(dá)矩陣的列：樣本
design <- model.matrix(~0+factor(group_list))
design
colnames(design)=levels(factor(group_list))
rownames(design)=colnames(exprSet)
design
# 實(shí)驗(yàn)設(shè)計(jì)矩陣的每一行對(duì)應(yīng)一個(gè)樣品的編碼椎木，
# 每一列對(duì)應(yīng)樣品的一個(gè)特征违柏。這里只考慮了一個(gè)因素兩個(gè)水平，
# 如果是多因素和多水平的實(shí)驗(yàn)設(shè)計(jì)香椎，會(huì)產(chǎn)生更多的特征漱竖，需要參考文檔設(shè)計(jì)。
# 構(gòu)建對(duì)比模型畜伐，比較兩個(gè)實(shí)驗(yàn)條件下表達(dá)數(shù)據(jù)
contrast.matrix<-makeContrasts(G3-con,levels = design)
#contrast.matrix<-makeContrasts(paste0(unique(group_list),collapse = "-"),levels = design)
contrast.matrix ##這個(gè)矩陣聲明馍惹，我們要把G3組跟con組進(jìn)行差異分析比較
##### 差異分析
##step1 線性模型擬合
fit <- lmFit(exprSet,design)
##step2 根據(jù)對(duì)比模型進(jìn)行差值計(jì)算
fit2 <- contrasts.fit(fit, contrast.matrix)
##step3 貝葉斯檢驗(yàn)
fit2 <- eBayes(fit2)
##step4 生成所有基因的檢驗(yàn)結(jié)果報(bào)告
tempOutput = topTable(fit2, coef=1, n=Inf)
##step5 用P.Value進(jìn)行篩選，得到全部差異表達(dá)基因
dif <- tempOutput[tempOutput[, "P.Value"]<0.01,]
# 顯示一部分報(bào)告結(jié)果
head(dif)
參考：
- 用limma對(duì)芯片數(shù)據(jù)做差異分析
- Bioconductor分析基因芯片數(shù)據(jù)
-
limFit函數(shù)文檔?limFit()
更新資料：
相關(guān)問(wèn)題請(qǐng)下面留言討論玛界。
差異分析比較的是什么
在使用limma時(shí)万矾，我一直對(duì)對(duì)比的事物存有疑惑，特別是當(dāng)你可能看到下面這種分析結(jié)果相同時(shí)：
#1:
library(CLL)
data(sCLLex)
library(limma)
design=model.matrix(~factor(sCLLex$Disease))
fit=lmFit(sCLLex,design)
fit=eBayes(fit)
options(digits = 4)
#topTable(fit,coef=2,adjust='BH')
> topTable(fit,coef=2,adjust='BH')
logFC AveExpr t P.Value adj.P.Val B
39400_at 1.0285 5.621 5.836 8.341e-06 0.03344 3.234
36131_at -0.9888 9.954 -5.772 9.668e-06 0.03344 3.117
33791_at -1.8302 6.951 -5.736 1.049e-05 0.03344 3.052
1303_at 1.3836 4.463 5.732 1.060e-05 0.03344 3.044
36122_at -0.7801 7.260 -5.141 4.206e-05 0.10619 1.935
36939_at -2.5472 6.915 -5.038 5.362e-05 0.11283 1.737
41398_at 0.5187 7.602 4.879 7.824e-05 0.11520 1.428
32599_at 0.8544 5.746 4.859 8.207e-05 0.11520 1.389
36129_at 0.9161 8.209 4.859 8.212e-05 0.11520 1.389
37636_at -1.6868 5.697 -4.804 9.355e-05 0.11811 1.282
#2:
library(CLL)
data(sCLLex)
library(limma)
design=model.matrix(~0+factor(sCLLex$Disease))
colnames(design)=c('progres','stable')
fit=lmFit(sCLLex,design)
cont.matrix=makeContrasts('progres-stable',levels = design)
fit2=contrasts.fit(fit,cont.matrix)
fit2=eBayes(fit2)
options(digits = 4)
topTable(fit2,adjust='BH')
logFC AveExpr t P.Value adj.P.Val B
39400_at -1.0285 5.621 -5.836 8.341e-06 0.03344 3.234
36131_at 0.9888 9.954 5.772 9.668e-06 0.03344 3.117
33791_at 1.8302 6.951 5.736 1.049e-05 0.03344 3.052
1303_at -1.3836 4.463 -5.732 1.060e-05 0.03344 3.044
36122_at 0.7801 7.260 5.141 4.206e-05 0.10619 1.935
36939_at 2.5472 6.915 5.038 5.362e-05 0.11283 1.737
41398_at -0.5187 7.602 -4.879 7.824e-05 0.11520 1.428
32599_at -0.8544 5.746 -4.859 8.207e-05 0.11520 1.389
36129_at -0.9161 8.209 -4.859 8.212e-05 0.11520 1.389
37636_at 1.6868 5.697 4.804 9.355e-05 0.11811 1.282
上述代碼資料來(lái)自jimmy
為什么第一種方式?jīng)]有做對(duì)比矩陣慎框，結(jié)果一致良狈！
大家運(yùn)行一下這些代碼就知道，兩者結(jié)果是一模一樣的笨枯。
而差異比較矩陣的需要與否薪丁，主要看分組矩陣如何制作的！design=model.matrix(~factor(sCLLex$Disease)) design=model.matrix(~0+factor(sCLLex$Disease))有本質(zhì)的區(qū)別Ｏ诰Ｑ鲜取！
前面那種方法已經(jīng)把需要比較的組做出到了一列洲敢，需要比較多次阻问，就有多少列，第一列是截距不需要考慮沦疾，第二列開(kāi)始往后用coef這個(gè)參數(shù)可以把差異分析結(jié)果一個(gè)個(gè)提取出來(lái)称近。
而后面那種方法，僅僅是分組而已哮塞，組之間需要如何比較刨秆，需要自己再制作差異比較矩陣，通過(guò)makeContrasts函數(shù)來(lái)控制如何比較忆畅！
--- jimmy
另一個(gè)問(wèn)題：這兩種方法哪一種更可取呢衡未？
在我沒(méi)有實(shí)際操作之前，我覺(jué)得簡(jiǎn)單的更清爽家凯，適用缓醋，但適用后我的結(jié)論是第二種各種可取。在前幾天的一次分析中绊诲，我將差異比較的兩個(gè)水平分為：High和Low送粱，結(jié)果分析的差異基因fold change反了！在沒(méi)有顯式指定的情況下掂之，我們難以真正確定我們比對(duì)的結(jié)果是High-Low還是Low-High抗俄。另外，后一種方法更利于將差異的比較過(guò)程程序化世舰。
最后动雹，再?gòu)?qiáng)調(diào)一下流程：
基因芯片的差異表達(dá)分析主要有 構(gòu)建基因表達(dá)矩陣、構(gòu)建實(shí)驗(yàn)設(shè)計(jì)矩陣跟压、構(gòu)建對(duì)比模型（對(duì)比矩陣）胰蝠、線性模型擬合、貝葉斯檢驗(yàn)和生成結(jié)果報(bào)表 六個(gè)關(guān)鍵步驟震蒋。
