Analyzing RNA-seq data with DESeq2翻譯（3）

多因素分析設計

可以使用包含附加變數的設計公式分析影響count值的多個因素。實際上，DESeq2包可以分析任何能用固定效應項表示的實驗設計，如多因素、具有交互作用的設計、具有連續變數的設計等分析都是可以的。

通過向設計公式中添加變數，可以分析影響count值的其他因素。例如，如果對不同處理的樣本進行實驗批次的平衡，即在設計公式中包含batch factor，就可以提高發現由處理條件引起的差異的靈敏度。分析附加變數有很多種方法。

pasilla包中的數據（如下所示）包含處理條件（condition），以及測序類型的信息（type）。

我們用多因素設計公式對這些數據進行重新分析

ddsMF <-ddslevels(ddsMF$type)##[1]"paried-end" "single-red"levels(ddsMF$type) <- sub("-.* " ,"", levels(ddsMF$type))levels(ddaMF$type)##[1]"paired" "single"

我們可以把不同測序類型的影響計算在內，得出其對處理的貢獻。

**因為condition是我們感興趣的變數，所以要放在設計公式的最後**這樣result函數進行比較的時候會默認根據condition比較，除非使用contrast或者name參數進行了其他指定。

我們重新進行DESeq分析：

design(ddsMF) <- formula(~ type + condition) #這個公式代表平衡type的影響後分析condtion的差異ddsMF <- DESeq(ddsMF)resMF <- results(ddsMF)head(resMF)

可以同樣提取type變數的log2 fold change，p-value以及adjusted p-value。results函數的contrast參數可以設定三個字元串向量：變數名稱，log2比率中作為分子的factor levle名，log2比率中作為分母的factor level名。 contrast也可以採用其他形式，具體請參考results函數的幫助信息。

resMFType <- results(ddsMF,contrast=c("type", "single", "paired"))head(resMFType)

如果變數是連續的或交互項（請參閱交互部分），則可以使用results函數中的name參數提取結果，其中變數名稱是resultsNames（dds）返回的元素之一。

數據轉換及可視化

count數據轉換

為了下游數據分析及可視化的方便，count數據通常要進行轉換。最常用的轉換是進行對數轉換。由於基因的count值在某些條件下為0，其他條件下又不為0，因此一些人提倡使用偽計數（pseudocounts），轉換形式為：

[ y = log_2(n + n_0) ]

其中n是count 值，n_0是一個正的常數。

這裡我們介紹另外兩種轉換方法，它們與上述（n_0 ）參數等效。其中一種是利用方差穩定變換（variance stabilizing transformations，VST）的概念（Tibshirani 1988; Huber et al. 2003; Anders and Huber 2010)另一種是正則化的對數或rlog（Love， Huber和Anders 2014）。這兩種轉換都在log2尺度上產生轉換數據，其已經根據庫的大小或其他標準化因子進行了歸一化。

VST和rlog轉換的重點，是消除方差對平均值的依賴性，尤其是平均值低時計數數據的對數的高方差。VST和rlog都使用實驗範圍內的均值方差趨勢，以便轉換數據已消除實驗範圍內的趨勢。請注意，我們不需要或者要求有所的基因在轉換後具有完全相同的方差。事實上，如下圖所示，在轉換後，具有相同平均值的基因沒有完全相同的標準偏差，但是實驗範圍內的趨勢已經趨於平緩。正是那些行方差高於趨勢的基因是我們能夠將樣本聚類到感興趣的分組中。

盲分散估計

vst和rlog這兩個函數有一個參數叫blind，用來確定變換是否對設計公式指定的信息可知。當blind=TRUE（默認情況），函數將僅使用截距重新估計色散。這樣設置的好處是在比較樣本時以完全不偏於實驗組信息的方式進行比較。例如對樣本進行質量控制時（下面會詳細介紹）。

但是，**如果預計多數基因的count值差異很大，並且這種差異可以通過實驗設計來解釋，並且希望為下游分析轉換數據，則盲分散不是合適的選擇**。在這種情況下，使用盲色散估計會導致大量的色散估計，因為它將本來歸因於實驗設計的差異歸因為雜訊，並且會導致轉換後的數據過度收縮。將blind設置為FALSE，已經估計的色散將被用於執行變換，如果它不存在，則將使用當前的設計公式來估計它們。請注意，只有平均色散趨勢線擬合的色散估計值用於轉換（整個實驗中色散平均值的整體依賴性）。因此，在轉換中將blind設置為FALSE仍然適用於大部分不使用哪個樣本在哪個實驗組中這種信息的情況。

提取轉換後的值

轉換函數返回一個DESeqTransform類的對象，它是RangedSummarizedExperiment的一個子類。 對於在新創建的DESeqDataSet上運行的約20個樣本，rlog轉換可能需要30秒，而VST轉換隻需要不到1秒。 當使用blind = FALSE並且DESeq函數已經運行時，運行時間會更短，因為不需要重新估計色散值。 assay函數用於提取標準化值的矩陣。

vsd <- vst(dds, blind=FALSE)rld <-rlog(dds, blind=FALSE)head(assay(vsd), 3)

VST轉換（方差穩定轉換variance stabilizing transformation）

上面例子中的代碼使用了一個參數來擬合分散。 在這種情況下，vst函數使用方差穩定變換的閉合表達式。 如果使用局部擬合（選項fitType =「locfit」來估計分散體），則使用數值積分。 轉換後的數據應近似方差穩定，並且還包括對尺寸因子或歸一化因子進行校正。 轉換後的數據在log2範圍內的大計數。

rlog轉換（規範化的log轉換）

函數rlog代表正則化log轉換，通過擬合一個模型，將每個樣本的一個項和一個根據數據估計的係數的先驗分布進行擬合，將原始計數數據轉換為log2尺度。這與DESeq和nbinomWaldTest所使用的對數倍數變化是相同的收縮（有時稱為正則化或調節）。結果數據包含的元素定義如下：

[ log_2（q_ {ij}）= beta_ {i0} + beta_ {ij} ]

其中（q_ {ij} ）是與基因i和樣本j的片段的預期真實濃度成比例的參數（參見下面的公式），（ beta_ {i0} ）是不經歷收縮的截距，和（ beta_ {ij} ）是基於整個數據集上的色散平均趨勢向零收縮的樣本特定效應。這種趨勢典型地表現為低計數的高分散性，因此這些基因表現出更高的收縮率。請注意，由於（q_ {ij} ）表示大小因子（s_j ）被分開之後的平均值（ mu_ {ij} ）的部分，因此很明顯rlog轉換具有固有的考慮到測序深度的差異。如果沒有先驗，這個設計矩陣將導致一個非獨特的解決方案，然而，在非截距測試中增加一個先驗的解決方案可以找到一個獨特的解決方案。

轉換對方差的影響

下圖使用移位對數轉換（shifted logarithm transformation）、正則化對數轉換和方差穩定轉換，繪製了轉換後的數據在樣本間的平均值與均值的標準差。 移位對數在較低的計數範圍內具有較高的標準偏差，正則化的對數在較小程度上較高，而對於方差穩定的數據，標準偏差在整個動態範圍內大致恆定。

請注意，這些圖中的垂直軸是所有樣本方差的平方根，因此包括由於實驗條件導致的方差。 雖然平均值方差平方根曲線可能看起來像是這種變換的目標，但由於實驗條件，在數據集有很多真實差異的情況下，這可能是不合理的。

# this gives log2(n + 1)ntd <- normTransform(dds)library("vsn")meanSdPlot(assay(ntd))

meanSdPlot(assay(vsd))

meanSdPlot(assay(rld))

PS.好多統計理論，背景知識太多。