9.1 背景介紹
與早期的微陣列不同郭宝，大多數(shù)數(shù)據(jù)現(xiàn)在都是單通道類型。單通道數(shù)據(jù)來(lái)自流行的微陣列技術(shù)，如Affymetrix巧勤，Illumina或Agilent。RNA-Seq的新技術(shù)也會(huì)生成單通道數(shù)據(jù)弄匕，因此本章中的所有內(nèi)容都可以在使用voom函數(shù)預(yù)處理數(shù)據(jù)后應(yīng)用于RNA-Seq分析[15]颅悉。單通道數(shù)據(jù)可能與普通單變量線性模型或方差分析非常相似。微陣列數(shù)據(jù)的差異在于幾乎總是需要提取特定的感興趣的對(duì)比迁匠，因此為R中的因子提供的標(biāo)準(zhǔn)參數(shù)通常不足剩瓶。
我們將在這里假設(shè)示例數(shù)據(jù)已被適當(dāng)?shù)仡A(yù)處理標(biāo)準(zhǔn)化并且可用作名為eset的ExpressionSet或EList對(duì)象。這樣的對(duì)象會(huì)有一個(gè)位置包含可以使用exprs(eset)提取的每個(gè)陣列上每個(gè)基因的對(duì)數(shù)表達(dá)值城丧。
9.2 兩組
最簡(jiǎn)單的單通道實(shí)驗(yàn)是比較兩組延曙。假設(shè)我們希望比較兩個(gè)野生型（Wt）小鼠與三個(gè)突變型（Mu）小鼠：
| FileName | Target |
|---|---|
| File1 | WT |
| File2 | WT |
| File3 | Mu |
| File4 | Mu |
| File5 | Mu |
形成設(shè)計(jì)矩陣有兩種不同的方式。我們可以二選一：
1.創(chuàng)建一個(gè)設(shè)計(jì)矩陣亡哄，其中包括突變型與野生型差異的系數(shù)枝缔，
2.創(chuàng)建一個(gè)設(shè)計(jì)矩陣，其中包括野生型和突變型小鼠分別的系數(shù)蚊惯，然后將其作為對(duì)比魂仍。
對(duì)于第一種方法拐辽，實(shí)驗(yàn)－對(duì)照參數(shù)化，設(shè)計(jì)矩陣應(yīng)該如下：
> design
WT MUvsWT
Array1 1 0
Array2 1 0
Array3 1 1
Array4 1 1
Array5 1 1
這里第一個(gè)系數(shù)估計(jì)野生型小鼠的平均對(duì)數(shù)表達(dá)值擦酌，并充當(dāng)了截距俱诸。第二個(gè)系數(shù)估計(jì)突變型和野生型之間的差異∩薏埃可以通過(guò)下面的命令發(fā)現(xiàn)差異表達(dá)基因
> fit <- lmFit(eset, design)
> fit <- eBayes(fit)
> topTable(fit, coef="MUvsWT", adjust="BH")
其中eset是包含表達(dá)值對(duì)數(shù)的ExpressionSet或matrix對(duì)象睁搭。第二種方法，設(shè)計(jì)矩陣應(yīng)該是
> design
WT MU
Array1 1 0
Array2 1 0
Array3 0 1
Array4 0 1
Array5 0 1
可以通過(guò)下面的命令發(fā)現(xiàn)差異表達(dá)基因
> fit <- lmFit(eset, design)
> cont.matrix <- makeContrasts(MUvsWT=MU-WT, levels=design)
> fit2 <- contrasts.fit(fit, cont.matrix)
> fit2 <- eBayes(fit2)
> topTable(fit2, adjust="BH")
對(duì)于第一種方法笼平，實(shí)驗(yàn)－對(duì)照參數(shù)化园骆，設(shè)計(jì)矩陣可以通過(guò)下面的命令計(jì)算
> design <- cbind(WT=1,MUvsWT=c(0,0,1,1,1))
或者通過(guò)
> Group <- factor(targets$Target, levels=c("WT","Mu"))
> design <- model.matrix(~Group)
> colnames(design) <- c("WT","MUvsWT")
對(duì)于第二種方法，分組參數(shù)化寓调，設(shè)計(jì)矩陣可以通過(guò)下面的命令計(jì)算
> design <- cbind(WT=c(1,1,0,0,0),MU=c(0,0,1,1,1))
或者通過(guò)
> design <- model.matrix(~0+Group)
> colnames(design) <- c("WT","MU")
9.3 若干組
上述分為兩組的方法很容易擴(kuò)展到任何數(shù)量的組锌唾。假設(shè)有三個(gè)RNA靶標(biāo)進(jìn)行比較，這三個(gè)目標(biāo)稱為 “RNA1”夺英，“RNA2” 和 “RNA3”晌涕，并且targets$Target列表示哪一個(gè)RNA與每個(gè)陣列進(jìn)行雜交⊥疵酰可以使用下面的命令創(chuàng)建適當(dāng)?shù)脑O(shè)計(jì)矩陣
> f <- factor(targets$Target, levels=c("RNA1","RNA2","RNA3"))
> design <- model.matrix(~0+f)
> colnames(design) <- c("RNA1","RNA2","RNA3")
為了使三組之間的所有成對(duì)比較都可以進(jìn)行余黎，需要運(yùn)行
> fit <- lmFit(eset, design)
> contrast.matrix <- makeContrasts(RNA2-RNA1, RNA3-RNA2, RNA3-RNA1,
+ levels=design)
> fit2 <- contrasts.fit(fit, contrast.matrix)
> fit2 <- eBayes(fit2)
RNA2與RNA1的差異基因列表可以通過(guò)下面的命令獲得
> topTable(fit2, coef=1, adjust="BH")
每個(gè)假設(shè)檢驗(yàn)的結(jié)果可以這樣獲得
> results <- decideTests(fit2)
每個(gè)比較中顯著基因數(shù)目的維恩圖可以通過(guò)下面的命令獲得
> vennDiagram(results)
統(tǒng)計(jì)量fit2$F和相應(yīng)的fit2$F.p.value將三對(duì)成對(duì)的比較整合到一次F檢驗(yàn)中。這相當(dāng)于每個(gè)基因的單向ANOVA分析载萌，除了基因間的平均殘差平方和已被調(diào)節(jié)惧财。以任何方式尋找在三種RNA之間變化的基因，需要找到p值小的基因扭仁。找出排名前30的基因：
> topTableF(fit2, number=30)
9.4 其他模型和區(qū)塊化
9.4.1 配對(duì)樣本
配對(duì)樣本通常出現(xiàn)在當(dāng)我們比較兩種實(shí)驗(yàn)條件時(shí)發(fā)生垮衷，每個(gè)樣品給予一種處理，這樣就與給定另一種處理的特定樣品自然配對(duì)乖坠。這是具有兩個(gè)區(qū)塊的區(qū)塊化的一個(gè)特殊情況搀突。與這種情況相關(guān)的經(jīng)典測(cè)試是配對(duì)t檢驗(yàn)。
假設(shè)進(jìn)行一個(gè)實(shí)驗(yàn)以比較新的實(shí)驗(yàn)處理（T）與對(duì)照（C）瓤帚。來(lái)自三個(gè)血親關(guān)系的六條狗被作為實(shí)驗(yàn)材料使用描姚。對(duì)于每個(gè)同胞對(duì)涩赢，一條狗被給予實(shí)驗(yàn)處理戈次，而另一條狗作為對(duì)照。這生成了目標(biāo)框架：
| FileName | SibShip | Treatment |
|---|---|---|
| File1 | 1 | C |
| File2 | 1 | T |
| File3 | 2 | C |
| File4 | 2 | T |
| File5 | 3 | C |
| File6 | 3 | T |
可以通過(guò)在線性模型中允許同胞對(duì)效應(yīng)來(lái)計(jì)算適中配對(duì)t檢驗(yàn)：
> SibShip <- factor(targets$SibShip)
> Treat <- factor(targets$Treatment, levels=c("C","T"))
> design <- model.matrix(~SibShip+Treat)
> fit <- lmFit(eset, design)
> fit <- eBayes(fit)
> topTable(fit, coef="TreatT")
9.4.2 區(qū)塊化
上述用于配對(duì)樣品的方法可以應(yīng)用于有批次效應(yīng)的實(shí)驗(yàn)或進(jìn)行區(qū)塊化實(shí)驗(yàn)的任何情況筒扒。實(shí)驗(yàn)處理可以通過(guò)使用以下形式的模型公式來(lái)區(qū)分塊之間的差異進(jìn)行調(diào)整：
> design <- model.matrix(~Block+Treatment)
在這種類型的分析中怯邪，僅在每個(gè)區(qū)塊內(nèi)進(jìn)行比較實(shí)驗(yàn)處理。
9.5 互動(dòng)模型：2×2因子設(shè)計(jì)
9.5.1 感興趣的問(wèn)題
因子設(shè)計(jì)是那些不止一個(gè)不同的實(shí)驗(yàn)維度和各種實(shí)驗(yàn)條件的組合花墩。假設(shè)從野生型和突變型小鼠提取細(xì)胞悬秉，這些細(xì)胞或者被刺激（S）或者未刺激（U）澄步。來(lái)自被處理細(xì)胞的RNA被提取后與微陣列雜交。我們將簡(jiǎn)單地假設(shè)陣列是單色陣列和泌，如Affymetrix村缸。考慮以下目標(biāo)框架：
| FileName | Strain | Treatment |
|---|---|---|
| File1 | WT | U |
| File2 | WT | S |
| File3 | Mu | U |
| File4 | Mu | S |
| File5 | Mu | S |
這里的兩個(gè)實(shí)驗(yàn)尺度或因素是種系和實(shí)驗(yàn)處理武氓。種系指定提取細(xì)胞的小鼠的基因型梯皿，實(shí)驗(yàn)處理指定細(xì)胞是否被刺激。觀察到所有四種種系和實(shí)驗(yàn)處理的組合县恕，所以這是一個(gè)因子設(shè)計(jì)东羹。使用下面的命令將種系/實(shí)驗(yàn)處理組合收集到一個(gè)矢量中很方便：
> TS <- paste(targets$Strain, targets$Treatment, sep=".")
> TS
[1] "WT.U" "WT.S" "Mu.U" "Mu.S" "Mu.S"
使用因子設(shè)計(jì)來(lái)決定感興趣的對(duì)比是特別重要的。我們?cè)谶@里假設(shè)實(shí)驗(yàn)者對(duì)下面幾點(diǎn)感興趣：
1.哪些基因在野生型細(xì)胞對(duì)刺激作出反應(yīng)忠烛，
2.哪些基因在突變型細(xì)胞對(duì)刺激作出反應(yīng)属提，
3.哪些基因在突變型細(xì)胞中與野生型細(xì)胞相比反應(yīng)程度不同。
因?yàn)檫@些是在分子生物學(xué)背景下最通常相關(guān)的問(wèn)題美尸。第一個(gè)問(wèn)題與WT.SvsWT.U的比較有關(guān)冤议，第二個(gè)與Mu.SvsMu.U相關(guān)，第三個(gè)涉及到差異之間的不同火惊，即(Mu.S-Mu.U) - (WT.S-WT.U)求类，術(shù)語(yǔ)被稱為相互作用。
9.5.2 單因素分析
我們首先描述了使用limma命令分析這個(gè)實(shí)驗(yàn)的簡(jiǎn)單方法屹耐，與兩個(gè)樣本設(shè)計(jì)分析類似尸疆。然后我們將繼續(xù)描述更經(jīng)典的使用因子模型公式的統(tǒng)計(jì)方法。所有考慮的方法是等同的惶岭，并產(chǎn)生相同的基本結(jié)果寿弱。最基本的方法是用一個(gè)有效的模型系數(shù)擬合四因子組合中的每一個(gè)，然后提取感興趣的比較作為對(duì)比：
> TS <- factor(TS, levels=c("WT.U","WT.S","Mu.U","Mu.S"))
> design <- model.matrix(~0+TS)
> colnames(design) <- levels(TS)
> fit <- lmFit(eset, design)
這是一個(gè)具有對(duì)應(yīng)于WT.U按灶，WT.S症革，MuU和MuS的四個(gè)系數(shù)的模型。我們可以提取三種感興趣的對(duì)比
> cont.matrix <- makeContrasts(
+ SvsUinWT=WT.S-WT.U,
+ SvsUinMu=Mu.S-Mu.U,
+ Diff=(Mu.S-Mu.U)-(WT.S-WT.U),
+ levels=design)
> fit2 <- contrasts.fit(fit, cont.matrix)
> fit2 <- eBayes(fit2)
我們可以使用topTable()函數(shù)來(lái)查看三個(gè)對(duì)比中的每一個(gè)的差異表達(dá)基因的列表鸯旁，或者
> results <- decideTests(fit2)
> vennDiagram(results)
同時(shí)查看所有三個(gè)對(duì)比噪矛。
對(duì)于大多數(shù)用戶而言，建議使用此方法铺罢，因?yàn)檎跍y(cè)試的比較是明確形成的艇挨。
9.5.3 嵌套互動(dòng)公式
R中的模型公式是非常靈活的，并提供許多可能的快捷鍵用于設(shè)置設(shè)計(jì)矩陣韭赘。不過(guò)缩滨，他們也需要高水平的統(tǒng)計(jì)知識(shí)來(lái)可靠地使用，并且它們并沒(méi)有在主要的R文檔中完整地描述。如果我們只是想測(cè)試上述前兩個(gè)問(wèn)題脉漏，一個(gè)設(shè)置設(shè)計(jì)矩陣的簡(jiǎn)單方法是使用一個(gè)嵌套相互作用：
> Strain <- factor(targets$Strain, levels=c("WT","Mu"))
> Treatment <- factor(targets$Treatment, levels=c("U","S"))
> design <- model.matrix(~Strain+Strain:Treatment)
> colnames(design)
[1] "(Intercept)" "StrainMu" "StrainWT:TreatmentS" "StrainMu:TreatmentS"
模型公式中的第一項(xiàng)是種系效應(yīng)苞冯。這引入了截距列到設(shè)計(jì)矩陣，用來(lái)估計(jì)野生型未刺激細(xì)胞的平均對(duì)數(shù)表達(dá)水平侧巨，種系列估計(jì)未刺激狀態(tài)下突變型與野生型的差異舅锄。模型公式中的第二項(xiàng)表示刺激和種系之間的相互作用。因?yàn)槟Ｐ椭械膶?shí)驗(yàn)處理沒(méi)有主效應(yīng)司忱，相互作用是嵌套到巧娱。它引入了第三和第四列到設(shè)計(jì)矩陣中，分別代表野生型和突變型小鼠的刺激效應(yīng)烘贴，與之前定義的SvsUinWT和SvsUinMu的對(duì)比完全一樣禁添。
> fit <- lmFit(eset, design)
> fit <- eBayes(fit)
> topTable(fit, coef=3)
會(huì)發(fā)現(xiàn)野生型小鼠響應(yīng)刺激的相關(guān)基因，
> topTable(fit, coef=4)
會(huì)發(fā)現(xiàn)突變型小鼠響應(yīng)刺激的相關(guān)基因桨踪。最后老翘，我們可以提取上面考慮的相互作用對(duì)比Diff
> fit2 <- contrasts.fit(fit, c(0,0,-1,1))
> fit2 <- eBayes(fit2)
> topTable(fit2)
這會(huì)發(fā)現(xiàn)突變型和野生型小鼠刺激后差異反應(yīng)的基因。
9.5.4 經(jīng)典交互模型
對(duì)因子設(shè)計(jì)的分析在統(tǒng)計(jì)學(xué)中有悠久的歷史锻离，一個(gè)因子模型公式體系被開(kāi)發(fā)用來(lái)促進(jìn)復(fù)雜設(shè)計(jì)的分析铺峭。重要的是要了解上述三個(gè)分子生物學(xué)問(wèn)題與任何一個(gè)用于因子設(shè)計(jì)的統(tǒng)計(jì)經(jīng)典參數(shù)化都沒(méi)有關(guān)系。所以我們一般推薦上面已經(jīng)考慮過(guò)的用于微陣列分析的方法汽纠。
假設(shè)我們以通常的統(tǒng)計(jì)方式進(jìn)行卫键，
> design <- model.matrix(~Strain*Treatment)
這創(chuàng)建了一個(gè)符合以下解釋的四個(gè)系數(shù)的設(shè)計(jì)矩陣：
| Coeffcient | Comparison | Interpretation |
|---|---|---|
| Intercept | WT.U | Baseline level of unstimulated WT |
| StrainMu | Mu.U-WT.U | Difference between unstimulated strains |
| TreatmentS | WT.S-WT.U | Stimulation effect for WT |
| StrainMu:TreatmentS | (Mu.S-Mu.U)-(WT.S-WT.U) | Interaction |
這被稱為實(shí)驗(yàn)－對(duì)照參數(shù)化。請(qǐng)注意虱朵，我們感興趣的一個(gè)比較莉炉，Mu.S-Mu.U，沒(méi)有體現(xiàn)碴犬，而Mu.U-WT.U的比較卻包括在內(nèi)絮宁，這可能不是我們直接感興趣的。我們需要使用對(duì)比來(lái)提取所有感興趣的比較：
> fit <- lmFit(eset, design)
> cont.matrix <- cbind(SvsUinWT=c(0,0,1,0),SvsUinMu=c(0,0,1,1),Diff=c(0,0,0,1))
> fit2 <- contrasts.fit(fit, cont.matrix)
> fit2 <- eBayes(fit2)
提取了WT刺激效應(yīng)作為第三系數(shù)服协，相互作用作為第四系數(shù)绍昂。突變刺激效應(yīng)被提取為原始模型第三和第四系數(shù)的總和。該分析產(chǎn)生的結(jié)果與之前的分析完全相同偿荷。
對(duì)于階乘實(shí)驗(yàn)窘游，更經(jīng)典的統(tǒng)計(jì)學(xué)方法是使用至零求和參數(shù)化。在R中是這樣實(shí)現(xiàn)的
> contrasts(Strain) <- contr.sum(2)
> contrasts(Treatment) <- contr.sum(2)
> design <- model.matrix(~Strain*Treatment)
以下解釋定義了四個(gè)系數(shù)：
| Coeffcient | Comparison | Interpretation |
|---|---|---|
| Intercept | (WT.U+WT.S+Mu.U+Mu.S)/4 | Grand mean |
| Strain1 | (WT.U+WT.S-Mu.U-Mu.S)/4 | Strain main effect |
| Treatment1 | (WT.U-WT.S+Mu.U-Mu.S)/4 | Treatment main effect |
| Strain1:Treatment1 | (WT.U-WT.S-Mu.U+Mu.S)/4 | Interaction |
該參數(shù)化具有許多吸引人的數(shù)學(xué)特性跳纳，是經(jīng)典的參數(shù)化忍饰，用于很多實(shí)驗(yàn)設(shè)計(jì)理論中的因子設(shè)計(jì)。但是它只定義了我們直接感興趣的一個(gè)系數(shù)棒旗，即相互作用喘批。三個(gè)我們感興趣的對(duì)比可以使用下面命令提取
> fit <- lmFit(eset, design)
> cont.matrix <- cbind(SvsUinWT=c(0,0,-2,-2),SvsUinMu=c(0,0,-2,2),Diff=c(0,0,0,4))
> fit2 <- contrasts.fit(fit, cont.matrix)
> fit2 <- eBayes(fit2)
結(jié)果與前三種方法相同。
這里描述的解決2×2因子問(wèn)題的各種方法是等價(jià)的铣揉，不同僅存在于線性模型的參數(shù)化選擇中饶深。擬合的模型對(duì)象fit僅在coefficients和相關(guān)組件中有所不同。剩余標(biāo)準(zhǔn)差fit$sigma逛拱，剩余自由度fit$df.residual和fit2的所有組件都是相同的敌厘，即使使用了不同的參數(shù)化。由于方法是等效的朽合，用戶可以自由選擇哪一個(gè)最直觀或最方便俱两。
9.6 時(shí)間過(guò)程實(shí)驗(yàn)
9.6.1 重復(fù)時(shí)間點(diǎn)
時(shí)間過(guò)程實(shí)驗(yàn)是在一些實(shí)驗(yàn)處理或刺激開(kāi)始以后的幾個(gè)時(shí)間點(diǎn)提取RNA的實(shí)驗(yàn)。分析時(shí)間過(guò)程實(shí)驗(yàn)的最佳方法取決于實(shí)驗(yàn)的性質(zhì)曹步，特別是不同時(shí)間點(diǎn)的數(shù)量宪彩。我們首先考慮具有相對(duì)較少重復(fù)時(shí)間點(diǎn)的實(shí)驗(yàn)。簡(jiǎn)單的這種類型的時(shí)間過(guò)程實(shí)驗(yàn)與9.3節(jié)中涉及的多組實(shí)驗(yàn)相似讲婚。
例如尿孔，我們這里考慮一個(gè)雙向?qū)嶒?yàn)，其中將要進(jìn)行兩種基因型的時(shí)間流程分析比較筹麸』詈希考慮下列目標(biāo)框架
| FileName | Target |
|---|---|
| File1 | wt.0hr |
| File2 | wt.0hr |
| File3 | wt.6hr |
| File4 | wt.24hr |
| File5 | mu.0hr |
| File6 | mu.0hr |
| File7 | mu.6hr |
| File8 | mu.24hr |
目標(biāo)是在0，6和24小時(shí)時(shí)間點(diǎn)從野生型和突變型動(dòng)物收集的RNA樣品物赶。這可以被看作是一個(gè)因子實(shí)驗(yàn)白指，但是更簡(jiǎn)單的方法是使用群體平均參數(shù)化。
> lev <- c("wt.0hr","wt.6hr","wt.24hr","mu.0hr","mu.6hr","mu.24hr")
> f <- factor(targets$Target, levels=lev)
> design <- model.matrix(~0+f)
> colnames(design) <- lev
> fit <- lmFit(eset, design)
野生型的哪些基因在6小時(shí)或24小時(shí)內(nèi)響應(yīng)酵紫？我們可以通過(guò)提取野生型時(shí)間之間的對(duì)比找到這些基因告嘲。
> cont.wt <- makeContrasts(
+ "wt.6hr-wt.0hr",
+ "wt.24hr-wt.6hr",
+ levels=design)
> fit2 <- contrasts.fit(fit, cont.wt)
> fit2 <- eBayes(fit2)
> topTableF(fit2, adjust="BH")
三個(gè)時(shí)間點(diǎn)之間的任何兩個(gè)對(duì)比將產(chǎn)生相同的結(jié)果。例如奖地，相同的基因列表將通過(guò)“wt.24hr-wt.0hr”得到而不是“wt.24hr-wt.6hr”状蜗。
哪些基因在突變型中響應(yīng)（即隨時(shí)間發(fā)生變化）？
> cont.mu <- makeContrasts(
+ "mu.6hr-mu.0hr",
+ "mu.24hr-mu.6hr",
+ levels=design)
> fit2 <- contrasts.fit(fit, cont.mu)
> fit2 <- eBayes(fit2)
> topTableF(fit2, adjust="BH")
相對(duì)于野生型鹉动，突變型中哪些基因隨時(shí)間推移響應(yīng)不同轧坎？
> cont.dif <- makeContrasts(
+ Dif6hr =(mu.6hr-mu.0hr)-(wt.6hr-wt.0hr),
+ Dif24hr=(mu.24hr-mu.6hr)-(wt.24hr-wt.6hr),
+ levels=design)
> fit2 <- contrasts.fit(fit, cont.dif)
> fit2 <- eBayes(fit2)
> topTableF(fit2, adjust="BH")
本節(jié)中描述的分析方法用于組蛋白脫乙酰酶抑制劑六點(diǎn)時(shí)間過(guò)程實(shí)驗(yàn)[25]。
9.6.2 多時(shí)間點(diǎn)
現(xiàn)在我們考慮一個(gè)每組有很多時(shí)間點(diǎn)的例子泽示。當(dāng)有很多時(shí)間點(diǎn)時(shí)缸血，假設(shè)表達(dá)隨著時(shí)間平滑地變化而不是從一個(gè)時(shí)間點(diǎn)到另一個(gè)時(shí)間點(diǎn)的離散式變化是合理的。這種類型的時(shí)間過(guò)程可以通過(guò)使用回歸樣條或多項(xiàng)式擬合時(shí)間趨勢(shì)來(lái)分析械筛。
考慮以下32行目標(biāo)框架：
| FileName | Group | Time |
|---|---|---|
| File1 | Control | 1 |
| File2 | Control | 2 |
| . | . | . |
| . | . | . |
| . | . | . |
| File16 | Control | 16 |
| File17 | Treat | 1 |
| File18 | Treat | 2 |
| . | . | . |
| . | . | . |
| . | . | . |
| File32 | Treat | 16 |
利用適量節(jié)點(diǎn)的三次樣條曲線表示特定條件下特定基因的時(shí)間過(guò)程可能是合理的捎泻。選擇有效的自由度在范圍3-5內(nèi)是合理的。為自然回歸樣條設(shè)置基礎(chǔ)：
> library(splines)
> X <- ns(targets$Time, df=5)
然后為對(duì)照組和實(shí)驗(yàn)組單獨(dú)擬合曲線：
> Group <- factor(targets$Group)
> design <- model.matrix(~Group*X)
> fit <- lmFit(y, design)
> fit <- eBayes(fit)
創(chuàng)建了一個(gè)具有12個(gè)參數(shù)的模型埋哟，最后5個(gè)對(duì)應(yīng)于相互作用笆豁，也就是組間曲線的差異郎汪。檢測(cè)對(duì)照組和實(shí)驗(yàn)組不同時(shí)間趨勢(shì)的基因：
> topTable(fit, coef=8:12)
對(duì)5 df上的每個(gè)基因進(jìn)行了一個(gè)調(diào)節(jié)F檢驗(yàn)，可以檢測(cè)出非常普遍的實(shí)驗(yàn)和對(duì)照曲線之間的差異闯狱。
注意煞赢，對(duì)于這種分析沒(méi)有必要進(jìn)行重復(fù)，也沒(méi)有必要在相同的時(shí)間點(diǎn)觀察兩個(gè)實(shí)驗(yàn)組哄孤。
9.7 多層次實(shí)驗(yàn)
我們考慮了配對(duì)比較照筑，考慮了兩個(gè)獨(dú)立組的比較。然而瘦陈，有實(shí)驗(yàn)結(jié)合了這兩種類型的比較凝危。
考慮使用以下目標(biāo)框架的單通道實(shí)驗(yàn)：
| FileName | Subject | Condition | Tissue |
|---|---|---|---|
| File01 | 1 | Diseased | A |
| File02 | 1 | Diseased | B |
| File03 | 2 | Diseased | A |
| File04 | 2 | Diseased | B |
| File05 | 3 | Diseased | A |
| File06 | 3 | Diseased | B |
| File07 | 4 | Normal | A |
| File08 | 4 | Normal | B |
| File09 | 5 | Normal | A |
| File10 | 5 | Normal | B |
| File11 | 6 | Normal | A |
| File12 | 6 | Normal | B |
該實(shí)驗(yàn)涉及6名受試者，其中3名患有該疾病的患者和3名正常受試者晨逝。每個(gè)受試者蛾默，我們有兩種組織類型的表達(dá)譜，A和B捉貌。
在分析本實(shí)驗(yàn)時(shí)趴生，我們要比較兩種組織類型。這個(gè)比較可以在受試者內(nèi)進(jìn)行昏翰，因?yàn)槊總€(gè)受試者都會(huì)產(chǎn)生兩種組織的值苍匆。我們也想比較患病對(duì)象和正常受試者，但這種比較是受試者之間的棚菊。
如果我們只想比較兩種組織類型浸踩，我們可以做一個(gè)配對(duì)的樣本比較。如果我們只想比較生病的和正常的受試者统求，我們可以做一個(gè)普通的兩組比較检碗。由于我們需要在受試者本身和之間進(jìn)行比較，因此需要把患者作為隨機(jī)效應(yīng)码邻≌厶辏可以使用duplicateCorrelation函數(shù)在limma中完成。
兩個(gè)實(shí)驗(yàn)因素條件和組織可以以多種方式處理像屋。在這里怕犁，我們假設(shè)將兩者合并成一個(gè)組合因素是方便的：
> Treat <- factor(paste(targets$Condition,targets$Tissue,sep="."))
> design <- model.matrix(~0+Treat)
> colnames(design) <- levels(Treat)
然后我們估計(jì)同一受試者的測(cè)量數(shù)據(jù)之間的相關(guān)性：
> corfit <- duplicateCorrelation(eset,design,block=targets$Subject)
> corfit$consensus
然后將該對(duì)象內(nèi)相關(guān)性輸入到線性模型擬合中：
> fit <- lmFit(eset,design,block=targets$Subject,correlation=corfit$consensus)
現(xiàn)在我們可以通過(guò)常規(guī)方式對(duì)實(shí)驗(yàn)條件進(jìn)行比較，例如：
> cm <- makeContrasts(
+ DiseasedvsNormalForTissueA = Diseased.A-Normal.A,
+ DiseasedvsNormalForTissueB = Diseased.B-Normal.B,
+ TissueAvsTissueBForNormal = Normal.B-Normal.A,
+ TissueAvsTissueBForDiseased = Diseased.B-Diseased.A,
+ levels=design)
然后計(jì)算這些對(duì)比和緩和t檢驗(yàn)：
> fit2 <- contrasts.fit(fit, cm)
> fit2 <- eBayes(fit2)
然后
> topTable(fit2, coef="DiseasedvsNormalForTissueA")
將找到在A組織類型中正常和患病受試者之間差異表達(dá)的那些基因己莺，等等奏甫。
這個(gè)實(shí)驗(yàn)有兩個(gè)級(jí)別的差異性。首先凌受，人與人之間存在差異阵子，我們稱之為受試者之間的層次。那么同一受試者重復(fù)測(cè)量的可變性稱為受試者內(nèi)層次胜蛉。受試者間層次通常認(rèn)為大于受試者內(nèi)層次挠进，因?yàn)楹笳弑徽{(diào)整為受試者之間的基線差異色乾。在這里，組織類型之間的比較可以在受試者內(nèi)進(jìn)行领突，因此應(yīng)該比病人和正常人之間的比較要精確得多暖璧。
