寫在前面的廢話
不研究不知道，一研究嚇一跳，原來DESeq2這么復雜，這10000多的引用量真不是吹的……
廢話超多系列
DESeq2的差異表達分析涉及多個步驟，具體步驟參見下面流程圖中的藍色部分
簡單地說，DESeq2將對原始reads進行建模，使用標準化因子(scale factor)來解釋庫深度的差異。然后，DESeq2估計基因的離散度，并縮小這些估計值以生成更準確的離散度估計，從而對reads count進行建模。最后，DESeq2擬合負二項分布的模型，并使用Wald檢驗或似然比檢驗進行假設檢驗。
為什么說DESeq2復雜呢？因為上一篇文章講了七個步驟，也僅僅只是完成了這個流程圖中Estimate size facors這一步。
在使用DESeq2進行基因表達差異分析之前，最重要的是明確我們的研究目的，了解數(shù)據(jù)中的變異來源。一旦我們了解了數(shù)據(jù)的主要變異來源，就可以在分析之前提前移除它們，或者通過將這些變量包含在統(tǒng)計模型的公式中對它們進行分析。
DESeq2在進行分析之前需要我們自己書寫公式，以便讓軟件知道變異來源以及在差異表達分析中，我們感興趣的地方。以下面這個數(shù)據(jù)為例：
如果我們知道性別是數(shù)據(jù)中一個比較顯著的變異來源，那么我們就需要將sex寫入到統(tǒng)計模型的公式中。公式應該包含數(shù)據(jù)中的所有因素，這些因素解釋了數(shù)據(jù)中主要的變化來源，其中公式中的最后一個因素，應為我們最為關注的因素
比如，我們想要知道treatment的影響，其中sex和age是主要的變異來源，那么我們的公式則應該為design <- ~sex + age + treatment
公式中的波浪線應該在所有的代數(shù)式之前，從而告訴DESeq2在進行差異表達分析時，使用后面的公式。而公式中代數(shù)的名稱應該與數(shù)據(jù)框中的列名相匹配。
此外，DESeq2還允許我們研究變異之間的交互作用。比如，我們想知道sex對于treatment的影響，那么我們的公式就應該是design <- ~ sex + age + treatment + sex:treatment
此處需要注意，因為我們關注的是sex對于treatment的交互作用，因此
sex:treatment應該放在公式的最末尾
接下來就是無腦運行軟件，進行差異表達分析。
首先創(chuàng)建一個DESeq2Dataset對象
dds <- DESeqDataSetFromTximport(txi, colData = meta, design = ~ sampletype)
- txi是reads count的矩陣，每一行是一個基因，每一列是一個樣本
- colData則是一個因子數(shù)據(jù)框，每一個因子表示一個樣本，相同處理方式的樣本采用同樣的factor
- design就是剛剛上面所介紹的統(tǒng)計模型的公式
txi和colData的描述可能有點抽象，這里舉一個例子進一步說明：
下面就是txi應該有的格式：
倘若txi如上圖所示，則colData則應通過下述代碼得到
meta <- factor(rep(c('WT','KO'), each=3))
meta <- data.frame(row.names=colnames(txi), meta)
接下來進行差異表達分析，我們調(diào)用DESeq()函數(shù)即可
dds <- DESeq(dds)
這一步通過調(diào)用DESeq()，將軟件的運行結果重新賦給了dds對象。雖然我們僅僅用了一個命令，但是其中涉及到了多個步驟。軟件運行的輸出信息打印出了它所執(zhí)行的各個步驟：
estimating size factors
estimating dispersions
gene-wise dispersion estimates
mean-dispersion relationship
final dispersion estimates
fitting model and testing
除了這種無腦式的一鍵調(diào)用，DESeq2還提供了一些單獨的功能，可以讓我們一步一步地執(zhí)行工作流中的每一步。接下來我們詳細看看這幾個步驟的原理
1. estimating size factors
這一步也就是上一篇文章所說的文庫矯正，通過取log，找中位數(shù)，減少異常值對scale factor的影響，從而找到一個合適的scale factor。
我們要想單獨運行這一步，可以使用函數(shù)estimateSizeFactors()，示例如下：
dds <- DESeqDataSetFromTximport(txi, colData = meta, design = ~ sampletype)
## 計算sizefactors
dds.sizefactor <- estimateSizeFactors(dds)
## 如果想要知道具體的sizefactor是多少，可以使用sizeFactors()函數(shù)
sizeFactors(dds.sizefactor )
sizeFactors()函數(shù)除了可以查看表達矩陣評估得到的具體的sizefactor，還可以給一個DESeqDataSet對象的sizefactor賦值，這樣DESeq2在對DESeqDataSet對象進行差異表達分析時，就可以使用這個賦值的sizefactor進行后續(xù)分析。
dds <- DESeqDataSetFromTximport(txi, colData = meta, design = ~ sampletype)
# WT和KO兩種處理，六個樣本的sizefactor分別是2，2，2，8，9，10
sizefactor <- c(2,2,2,8,9,10)
names(sizefactor) <- c('WT_1','WT_2','WT_3','KO_1','KO_2','KO_3')
sizeFactors(dds) <- sizefactor
所以，如果你的樣本加有spike-in，你通過各種方式最終得到了一個
scale factor，也可以通過這種方式賦給你的基因表達數(shù)據(jù)集
2. Estimate gene-wise dispersion
差異表達分析的第二步是對數(shù)據(jù)離散程度的評估，在RNA-seq的reads count數(shù)據(jù)中，我們需要知道兩點：
- 為了確定差異表達的基因，我們需要根據(jù)組內(nèi)(重復之間)的方差來確定基因的表達值在組間是否有顯著差異
- 組內(nèi)(重復之間)的變異需要考慮到方差隨表達量的平均值增加的情況，如下圖所示(每個黑點是一個基因)。
為了更加準確的確定差異表達基因，DESeq2需要解釋方差和均值的關系。從上圖可知，在低表達的基因中，它們的方差也更低，因此DESeq2不希望差異表達基因都是低表達基因。
DESeq2使用離散度(dispersion)作為方差的度量方式，離散度既可以解釋基因表達值的方差也可以解釋基因的平均表達值。其具體公式為：Var = μ + α*μ^2。其中Var表示方差，μ表示均值，α表示離散度。因此我們可以得到這么一個關系
| 離散度 | |
|---|---|
| 方差增加 | 離散度增加 |
| 平均值增加 | 離散度降低 |
那么在表達水平較高的基因中，離散度的平方根 就等于方差系數(shù)
。其中σ是標準差，μ是平均值。那么α=0.01就意味著，在樣本生物重復之間存在著10%的標準差。
表達水平較高時，μ對于公式的影響顯著小于
因此，具有相同平均值的基因的離散度僅根據(jù)其方差而存在差異，離散度反映了一個給定平均值的基因表達的差異程度。
那么接下來一個比較重要的問題，就是如何將離散度與我們的模型建立聯(lián)系呢？為了更精確的為我們的數(shù)據(jù)建模，我們需要知道每個基因組內(nèi)方差的精確評估值。
但是，生物這個行當，一般3個重復就了不起了，6個重復就頂破天了?？雌饋硭坪跬Χ嗟模沁h遠不夠，這就導致我們得到的組內(nèi)方差是相當?shù)牟豢煽俊?/p>
6個重復都不一定夠，所以那些做RNA-seq一個處理只有一個重復的同學，你們是想搞事情么？
針對這個問題，DESeq2使用一種叫作shrinkage的方法，共享基因之間的表達信息，根據(jù)基因的表達水平生成更為準確的方差估計。DESeq2假設具有相似表達水平的基因也具有相似的方差。這樣DESeq2就可以基于基因的平均表達水平和離散度來評估每個基因的離散度。
在下文里我就暫且把
shrinkage翻譯為收縮
3. Fit curve to gene-wise dispersion estimates
DESeq2的第三步就是根據(jù)基因的離散度擬合一個曲線。那么為什么要做擬合呢？不同的基因生物學重復中存在不同的方差，但是，在所有的基因中，將會有一個合理的離散分布。
這個曲線如下圖紅線所示，其中紅線的橫坐標是基因的表達強度，縱坐標是理論離散值。而每一個黑點的橫坐標是基因的平均表達水平，縱坐標是經(jīng)過最大似然估計的離散值。
4. Shrink gene-wise dispersion estimates toward the values predicted by the curve
有了擬合曲線，接下來就是對基因表達水平的離散度進行矯正，即：將基因的實際離散度向紅線收縮(shrink)。當樣本量較小時，該曲線可以讓我們更為準確的識別差異表達基因。既然知道要將基因的離散度向紅色曲線收縮，那么收縮多少比較合適呢？有兩點需要考慮：
- 基因的離散度距離紅色曲線的距離
- 樣本量（樣本量越大，則收縮的越少）
這種方法在差異表達分析時，可以極大的減少數(shù)據(jù)的假陽性。離散度較低的基因朝著理論值收縮，從而得到一個更為準確的離散值。而那些離散度較高的基因，則不能無腦朝著理論值收縮。
這是因為這類基因可能不遵循建模假設，并且由于生物學或技術原因使得這類基因與其他基因具有更高的方差。DESeq2識別這類基因后，將不采用shrink方法對它們進行處理。下圖中藍色圓圈圈住的便是這類基因。
5. GLM fit for each gene
DESeq2的第五步，對每個基因使用廣義線性模型進行擬合，我自己也沒搞明白，就不在這里獻丑了
寫在后面的話
也許你會問既然DESeq()函數(shù)可以一鍵實現(xiàn)這么多操作，為什么還要一步一步單拎出來做呢？主要是我自己有點強迫癥，知其然不知其所以然的感覺太痛苦了……更何況，有了spike-in的教訓，我還敢無腦運行軟件么？
參考資料
- https://hbctraining.github.io/ ：Introduction to DGE
- Beginner's guide to using the DESeq2 package
