8.1 背景介紹
limma軟件包使用稱(chēng)為線性模型的方法來(lái)分析設(shè)計(jì)的微陣列實(shí)驗(yàn)稳捆。這種方法允許分析非常一般的實(shí)驗(yàn)档礁，就像分析簡(jiǎn)單的重復(fù)實(shí)驗(yàn)一樣容易卓箫。這種方法在[42, 55]中進(jìn)行了概述闯睹。該方法需要指定一個(gè)或兩個(gè)矩陣罩扇。第一個(gè)是設(shè)計(jì)矩陣，表明有每個(gè)陣列使用了哪些RNA樣本摔踱。第二個(gè)是對(duì)比矩陣虐先，它規(guī)定了想要在RNA樣品之間進(jìn)行哪些比較怨愤。對(duì)于非常簡(jiǎn)單的實(shí)驗(yàn)派敷，可能不需要指定對(duì)比矩陣。
方法的理念如下，你必須從將數(shù)據(jù)擬合到一個(gè)線性模型開(kāi)始篮愉，完全建模數(shù)據(jù)的系統(tǒng)部分腐芍。該模型由設(shè)計(jì)矩陣指定。設(shè)計(jì)矩陣的每一行對(duì)應(yīng)于實(shí)驗(yàn)的陣列试躏，每一列對(duì)應(yīng)于描述實(shí)驗(yàn)中RNA源的系數(shù)猪勇。Affymetrix或單通道數(shù)據(jù)，或者是有共同參照的雙色數(shù)據(jù)颠蕴，將需要許多系數(shù)泣刹，和RNA源的數(shù)量一致。直接設(shè)計(jì)的雙色數(shù)據(jù)犀被，需要的系數(shù)比RNA源數(shù)目少一個(gè)椅您，除非你想要估計(jì)每個(gè)基因的染料效應(yīng)，在這種情況下RNA源的數(shù)量和系數(shù)的數(shù)量相同寡键。任何一套獨(dú)立的系數(shù)都可以掀泳，只要它們足以描述實(shí)驗(yàn)。該步驟的主要目的是估計(jì)數(shù)據(jù)的變異性西轩，因此系統(tǒng)部分需要進(jìn)行建模员舵，從而可以和隨機(jī)變異進(jìn)行區(qū)分。
在實(shí)踐中藕畔，從可能想要回答的問(wèn)題角度來(lái)說(shuō)马僻，要求系數(shù)與RNA源數(shù)目完全相同太過(guò)苛刻。你可能比較感興趣的是更多或更少的RNA源之間的比較注服。因此巫玻，提供對(duì)比步驟，以便可以利用初始系數(shù)祠汇，并以想要回答任何問(wèn)題的方式進(jìn)行比較仍秤，不管這些可能有多少。
如果有來(lái)自Affymetrix實(shí)驗(yàn)的數(shù)據(jù)可很，比如來(lái)自單通道點(diǎn)微陣列或來(lái)自使用共同參照的點(diǎn)微陣列诗力，那么線性建模與普通方差分析或多重回歸分析相同，只是需要為每個(gè)基因擬合一個(gè)模型我抠。具有這種類(lèi)型的數(shù)據(jù)苇本，可以像單變量數(shù)據(jù)普通建模一樣創(chuàng)建設(shè)計(jì)矩陣。如果你使用直接設(shè)計(jì)的點(diǎn)狀微陣列的數(shù)據(jù)菜拓，即沒(méi)有共同參照的連接設(shè)計(jì)瓣窄，那么線性建模方法非常強(qiáng)大，但是設(shè)計(jì)矩陣的創(chuàng)建可能需要更多的統(tǒng)計(jì)知識(shí)纳鼎。
對(duì)于統(tǒng)計(jì)分析和評(píng)估差異表達(dá)俺夕，limma使用經(jīng)驗(yàn)貝葉斯方法以減輕估計(jì)的對(duì)數(shù)倍變化的標(biāo)準(zhǔn)誤差裳凸。這導(dǎo)致更穩(wěn)定的推導(dǎo)和改進(jìn)的效力，特別是對(duì)于具有少量陣列的實(shí)驗(yàn)[42, 26]劝贸。對(duì)于陣列內(nèi)部有重復(fù)點(diǎn)的陣列姨谷，limma使用合并的相關(guān)方法充分利用重復(fù)點(diǎn)[39]。
8.2單通道設(shè)計(jì)
Affymetrix數(shù)據(jù)通常使用affy包進(jìn)行標(biāo)準(zhǔn)化映九。我們?cè)谶@里假設(shè)數(shù)據(jù)可用作名為eset的ExpressionSet對(duì)象梦湘。這樣的對(duì)象將包含一個(gè)插槽，表示每個(gè)陣列上每個(gè)基因的對(duì)數(shù)表達(dá)值件甥，可以使用exprs(eset)提取捌议。Affymetrix和其他單通道微陣列數(shù)據(jù)非常可能用普通線性模型或anova模型分析引有。與微陣列數(shù)據(jù)的差異在于它幾乎總是需要提取特定的感興趣的對(duì)比禁灼，因此在R中為因子提供的標(biāo)準(zhǔn)參數(shù)通常不足。
limma中有很多方法可以分析一個(gè)復(fù)雜實(shí)驗(yàn)轿曙。一個(gè)直接的策略是設(shè)置最簡(jiǎn)單的設(shè)計(jì)矩陣弄捕，然后從擬合中提取感興趣的對(duì)比。
假設(shè)有三個(gè)RNA源需要進(jìn)行比較导帝。假設(shè)前三個(gè)陣列與RNA1雜交守谓，接下來(lái)兩個(gè)與RNA2雜交，其次三個(gè)與RNA3雜交您单。假設(shè)RNA源之間的所有兩兩比較都是值得關(guān)注的斋荞。我們假設(shè)數(shù)據(jù)已經(jīng)被標(biāo)準(zhǔn)化并存儲(chǔ)在ExpressionSet對(duì)象中，例如通過(guò)
> data <- ReadAffy()
> eset <- rma(data)
下面的命令可以創(chuàng)建適當(dāng)?shù)脑O(shè)計(jì)矩陣虐秦，并使用線性模型進(jìn)行擬合
> design <- model.matrix(~ 0+factor(c(1,1,1,2,2,3,3,3)))
> colnames(design) <- c("group1", "group2", "group3")
> fit <- lmFit(eset, design)
為了三組之間的兩兩比較平酿，可以創(chuàng)建合適的對(duì)比矩陣
> contrast.matrix <- makeContrasts(group2-group1, group3-group2, group3-group1, levels=design)
> fit2 <- contrasts.fit(fit, contrast.matrix)
> fit2 <- eBayes(fit2)
第2組與第1組中差異表達(dá)的基因列表排名可以通過(guò)下面的命令獲得
> topTable(fit2, coef=1, adjust="BH")
每個(gè)假設(shè)檢驗(yàn)的結(jié)果可以使用
> results <- decideTests(fit2)
代表每個(gè)比較中顯著基因數(shù)目的維恩圖可以通過(guò)下面的命令獲得
> vennDiagram(results)
8.3 共同參照設(shè)計(jì)
現(xiàn)在考慮雙色微陣列實(shí)驗(yàn)，其中陣列已經(jīng)使用了共同參照悦陋。這樣的實(shí)驗(yàn)可以非常類(lèi)似于Affymetrix實(shí)驗(yàn)進(jìn)行分析蜈彼，只是必須進(jìn)行染料交換。最簡(jiǎn)單的方法是使用modelMatrix()函數(shù)和目標(biāo)文件設(shè)計(jì)矩陣俺驶。例如幸逆，我們考慮Jo?lle Michaud，Catherine Carmichael和Hamish Scott博士在沃爾特伊麗莎豪爾研究所比較轉(zhuǎn)錄因子在人類(lèi)細(xì)胞系中的作用實(shí)驗(yàn)的一部分暮现。目標(biāo)文件如下：
> targets <- readTargets("runxtargets.txt")
> targets
SlideNumber Cy3 Cy5
1 2144 EGFP AML1
2 2145 EGFP AML1
3 2146 AML1 EGFP
4 2147 EGFP AML1.CBFb
5 2148 EGFP AML1.CBFb
6 2149 AML1.CBFb EGFP
7 2158 EGFP CBFb
8 2159 CBFb EGFP
9 2160 EGFP AML1.CBFb
10 2161 AML1.CBFb EGFP
11 2162 EGFP AML1.CBFb
12 2163 AML1.CBFb EGFP
13 2166 EGFP CBFb
14 2167 CBFb EGFP
在實(shí)驗(yàn)中还绘，綠色熒光蛋白（EGFP）被用作共同參照。利用一個(gè)腺病毒系統(tǒng)將各種轉(zhuǎn)錄因子轉(zhuǎn)運(yùn)到HeLa細(xì)胞的細(xì)胞核中栖袋。在這里我們考慮轉(zhuǎn)錄因子AML1拍顷，CBFbeta或兩者都考慮。生成一個(gè)簡(jiǎn)單的設(shè)計(jì)矩陣和線性模型：
> design <- modelMatrix(targets,ref="EGFP")
> design
AML1 AML1.CBFb CBFb
1 1 0 0
2 1 0 0
3 -1 0 0
4 0 1 0
5 0 1 0
6 0 -1 0
7 0 0 1
8 0 0 -1
9 0 1 0
10 0 -1 0
11 0 1 0
12 0 -1 0
13 0 0 1
14 0 0 -1
> fit <- lmFit(MA, design)
將每個(gè)轉(zhuǎn)錄因子與EGFP進(jìn)行比較塘幅，并且分別比較組合轉(zhuǎn)錄因子與AML1和CBFb是有意義的昔案。形成的適當(dāng)?shù)膶?duì)比矩陣如下：
> contrast.matrix <- makeContrasts(AML1,CBFb,AML1.CBFb,AML1.CBFb-AML1,AML1.CBFb-CBFb,
+ levels=design)
> contrast.matrix
AML1 CBFb AML1.CBFb AML1.CBFb - AML1 AML1.CBFb - CBFb
AML1 1 0 0 -1 0
AML1.CBFb 0 0 1 1 1
CBFb 0 1 0 0 -1
現(xiàn)在可以擴(kuò)展線性模型擬合尿贫，并計(jì)算經(jīng)驗(yàn)貝葉斯統(tǒng)計(jì)量：
> fit2 <- contrasts.fit(fit, contrasts.matrix)
> fit2 <- eBayes(fit2)
8.4 雙色直接設(shè)計(jì)
沒(méi)有共同參照的雙色設(shè)計(jì)需要最多的統(tǒng)計(jì)知識(shí)來(lái)選擇適當(dāng)?shù)脑O(shè)計(jì)矩陣。直接設(shè)計(jì)是沒(méi)有單一RNA源與每個(gè)陣列雜交的設(shè)計(jì)爱沟。例如帅霜，我們考慮由沃爾特伊麗莎豪爾研究所的Mireille Lahoud博士進(jìn)行的比較三種不同樹(shù)突細(xì)胞群體（DC）的基因表達(dá)實(shí)驗(yàn)匆背。
該實(shí)驗(yàn)涉及三個(gè)染料交換實(shí)驗(yàn)中的6個(gè)cDNA微陣列呼伸，每對(duì)比較兩種DC類(lèi)型。設(shè)計(jì)如上圖所示钝尸。目標(biāo)文件如下：
> targets
SlideNumber FileName Cy3 Cy5
ml12med 12 ml12med.spot CD4 CD8
ml13med 13 ml13med.spot CD8 CD4
ml14med 14 ml14med.spot DN CD8
ml15med 15 ml15med.spot CD8 DN
ml16med 16 ml16med.spot CD4 DN
ml17med 17 ml17med.spot DN CD4
這種實(shí)驗(yàn)的設(shè)計(jì)矩陣有很多有效的選擇括享，沒(méi)有一個(gè)絕對(duì)正確的選擇。我們選擇的設(shè)計(jì)矩陣如下：
> design <- modelMatrix(targets, ref="CD4")
Found unique target names:
CD4 CD8 DN
> design
CD8 DN
ml12med 1 0
ml13med -1 0
ml14med 1 -1
ml15med -1 1
ml16med 0 1
ml17med 0 -1
在這個(gè)設(shè)計(jì)矩陣中珍促，CD8和DN群體已經(jīng)與CD4群體進(jìn)行了比較铃辖。通過(guò)線性模型估計(jì)的系數(shù)將對(duì)應(yīng)于CD8與CD4的對(duì)數(shù)比（第一列）和DN與CD4的對(duì)數(shù)比（第二列）。
在表達(dá)數(shù)據(jù)適當(dāng)標(biāo)準(zhǔn)化后猪叙，使用下面命令擬合線性模型
> fit <- lmFit(MA, design)
線性模型現(xiàn)在可以被用來(lái)回答任何感興趣的問(wèn)題娇斩。對(duì)于這個(gè)實(shí)驗(yàn)，目的是在三個(gè)DC群體之間進(jìn)行兩兩比較穴翩。使用對(duì)比矩陣完成
> contrast.matrix <- cbind("CD8-CD4"=c(1,0),"DN-CD4"=c(0,1),"CD8-DN"=c(1,-1))
> rownames(contrast.matrix) <- colnames(design)
> contrast.matrix
CD8-CD4 DN-CD4 CD8-DN
CD8 1 0 1
DN 0 1 -1
對(duì)比矩陣可用于擴(kuò)展線性模型擬合犬第，然后計(jì)算經(jīng)驗(yàn)貝葉斯統(tǒng)計(jì)：
> fit2 <- contrasts.fit(fit, contrast.matrix)
> fit2 <- eBayes(fit2)
