(悄悄說一句，本次文末有資源，自助獲取)
這是思考問題的熊在果子學生信上發的第一篇帖子，屬於邀稿。
我比較喜歡熊的文字，讀起來循序漸進，漸入佳境，起承轉合，水到渠成。同時，跟他的相處中，我能感受到他高大上的情商，注意不僅僅是高情商，而是高大上的情商。我比較喜歡的一種對情商的定義是這樣的：
情商是指，能夠感受到自己以及他人的情緒改變，並知道這種情緒變化的來源，而且有能力在源頭上處理和疏通。
熊就屬於這個級別的，這種屬性也反應在他的文字當中。這個帖子很長，但是結論很簡單，就是用新版本的Deseq2就可以了，熊考慮到很多朋友讀不了長文，就直接在文章前面加了個太長不看版本，這種為他們考慮的態度，就是情商，好的情商，應該讓人感到舒服，而不是那些容易被識破的小聰明。
既然一句話就把主要思想講完了，那為什麼寫這麼長。這些人都屬於專業的生信工作人員，不僅僅要測試臺式機，伺服器的性能，還要體現多線程和單線程的速度，最終還要講清楚速度快的原因，關鍵是寫得很流暢。這就是專業的體現。
以下是正文太長不看版本
DESeq 最新版本（v1.25.9）針對大樣本分析的速度和之前相比有了質的飛越，果子老師的 1215 個 TCGA 樣本差異分析時間從 1200 分鐘縮短到單線程 20 分鐘，使用多線程的情況下最快 6 分鐘搞定。
前奏前一段時間洲更用統計學中的非參數檢驗方法解決了果子老師關於 1200 個 TCGA 樣本使用 DESeq2 分析耗時 20 小時的問題，然後在那篇文章的最後果子老師如此總結道：
這樣看來，洲更說的是對的。不過，dds 這裡還有一個比較好的點就是可以做 logFC 矯正，就像這個帖子裡面提到的，有一些 count 比較小，但是變化值很大的基因，會對 GSEA 分析產生影響。
從這段話裡可以看出他還是心心念念 DEseq2。念念不忘必有迴響，DEseq2 開發者聽到了果子老師心裡的吶喊，最近一次 GitHub 上 DEseq2 的更新就著手解決了這個問題。這段代碼的貢獻者稱他們解決了 DEseq2 處理大量樣本異常耗時問題，時間隨樣本量的變化趨勢目前已經做到了線性關係。
對比圖如下：
竟效果如何，接下來就用果子老師 1215 個樣本數據去測試。
第一次測試環境我選擇了自己的臺式機，16G 內存。不知道是出於對開發者的信任還是對電腦的信任或者是出於對果子老師代碼的信任。我選擇了直接複製了果子老師的代碼開啟並行模式進行分析。
在果子老師的當年的描述中，如此說道：
最困難的一步來了。燃燒你的小電腦。這裡還使用了並行化處理來解決速度的問題。但是，耗時依然很長，4 個小時以上，有可能是 10 個小時
（後來，據我私下和果子老師交流，他說實際運行時間遠遠高於 4 小時，為了顯示出自己機器性能的優越稍微在寫作手法上做了一些加工，人為縮短了這個時間）
這裡我並不知道果子老師小電腦的性能如何以及用了幾個線程，姑且按照 4 個小時（240 分鐘）來計算。首先要從 GitHub 安裝最新版本的包，安裝好之後查看包的版本是 1.25.9，沒有問題。
下是測試步驟，原始的 Rdata 來自果子老師，用構建好的 dds 直接做DESeq。
devtools::install_github("mikelove/DESeq2@ae7c6bd")
library(DESeq2)
library(BiocParallel)
dds <- DESeq(dds,parallel = T)
DESeq 雖然只是一個函數其實包含了若干個步驟，其中耗時很嚴重的是gene-wise dispersion estimates，這裡也是多線程開啟的地方。R 自動在我的臺式機上開啟了 6 個線程。起初我非常開心，新版本gene-wise dispersion estimates 這一步很快完成，目測只有不到 5 分鐘。看情況不出問題的話，全部分析應該會在 10 分鐘左右完成。
然後我就看了眼電腦的任務管理器，如下圖所示。
下來，就出問題了。因為 6 個線程工作再加上電腦只有 16G 內存，fitting model and testing到一半時我的心態伴隨著內存一起都崩了。到崩潰為止，一共運行了大約 11 分鐘。
> dds <- DESeq(dds,parallel = T)
estimating size factors
estimating dispersions
gene-wise dispersion estimates: 6 workers
mean-dispersion relationship
final dispersion estimates, fitting model and testing: 6 workers
error: arma::memory::acquire(): out of memory
Error: BiocParallel errors
element index: 4
first error: std::bad_alloc
In addition: There were 13 warnings (use warnings() to see them)
> proc.time() - t1
user system elapsed
42.21 253.32 702.48
事已至此，勿忘初心，遷移到伺服器上去試試。在伺服器安裝最新版本的 DEseq2 包，繼續運行同樣內容，並行開啟後，這次自動給我分配了 36 個線程，真棒！
接下來一切正常，運行時間如下所示：
> dds <- DESeq(dds,parallel = T)
estimating size factors
estimating dispersions
gene-wise dispersion estimates: 38 workers
mean-dispersion relationship
final dispersion estimates, fitting model and testing: 38 workers
-- replacing outliers and refitting for 9759 genes
-- DESeq argument 'minReplicatesForReplace' = 7
-- original counts are preserved in counts(dds)
estimating dispersions
fitting model and testing
> proc.time() - t1
user system elapsed
1720.880 207.212 379.153
全程跑完一共用時 6 分鐘多一點，和之前果子老師的多線程 4 個多小時的時間比，速度提升了應該有 40 多倍的樣子。這裡再次強調一下，據我私下和果子老師交流實際運行時間遠遠要高於 4 小時哦。如此看來，文章開頭的運行時間對比圖還比較靠譜，線性時間的複雜度大大縮短了處理大量樣本的時間。
單線程測試隨後，測試如果只開啟一個線程的情況（現在很少有單線程的機器，但是不少人都不知道 DEseq 可以多線程運行），如果不用最新版的話，據一些學員反映這個時間是在 20 到 30 個小時之內，如果用了最新版呢？單線程整個時間縮短到 22 分鐘左右。（從一個側面來說，這裡的多線程加速效果也並不是非常明顯。）
> dds <- DESeq(dds,parallel = F)
estimating size factors
estimating dispersions
gene-wise dispersion estimates
mean-dispersion relationship
final dispersion estimates
fitting model and testing
-- replacing outliers and refitting for 9759 genes
-- DESeq argument 'minReplicatesForReplace' = 7
-- original counts are preserved in counts(dds)
estimating dispersions
fitting model and testing
> proc.time() - t1
user system elapsed
1369.724 18.084 1387.730
從下文會提到的開發者更改內容來看，這次更新提速並沒有更改涉及到 DEseq2 模型和參數的內容，所以新老版本的分析結果是完全一致的，包括差異基因的個數和 p 值。
contrast <- c("sample","cancer","normal")
dd1 <- results(dds, contrast=contrast, alpha = 0.05)
library(dplyr)
library(tibble)
library(tidyr)
res <- dd1 %>%
data.frame() %>%
rownames_to_column("gene_id")
> table(res$padj < 0.01)
FALSE TRUE
20502 25189
回到文章開頭，這次的升級當然不是聽到果子老師遠在他鄉用微信公眾號發出的呼喚，而是開發者為了讓 DEseq2 可以更好的處理單細胞數據而準備的。在傳統的差異分析問題上，我們面對的是高深度和少量幾個樣本，但是在單細胞上，面對的是上千個細胞和相對很淺的測序深度。
代碼貢獻者在自己的博客中提到，在幾年前大家都還在爭論轉錄組測序數據中大量 count 為 0 的數據應該如何處理，但是最近已經有陰性對照的實驗數據表明這些 0 並不會有什麼問題。另外，已經有文章測試表明，傳統差異分析的很多工具在單細胞數據中依舊可以使用，甚至要比一些專門為單細胞開發的工具表現要好。如下圖所示：
是 DESeq2 開發者的實驗室發現，DESeq2 和其他一票軟體相比，在處理大量細胞時速度上完敗，細胞越多敗的越徹底。作為差異分析工具的大哥，當然也是要面子的。於是決定找出問題進行優化。
於代碼改動，可以從 GitHub 的 PR 中略知一二。
以下是 R 腳本的一個優化示例：
在core.R中，舊版本的代碼如下圖所示
版本中對應的代碼如下圖所示：
獨看幾行代碼可能有些難以理解，在這裡盡力解釋一下。
上面代碼提到的 modelMatrix 參數來自於getModelMatrix()這個函數，內容如下：
getModelMatrix <- function(object) {
if (is(design(object), "matrix")) {
design(object)
} else if (is(design(object), "formula")) {
stats::model.matrix.default(design(object), data=as.data.frame(colData(object)))
}
}
其中的參數object就是我們實際數據中的 dds，這裡首先判斷 design(object) 的數據類型，以果子老師的數據為例。
> is(design(dds))
[1] "formula" "oldClass"
果子老師的數據中design(object) 是 formula，那麼就會運行stats::model.matrix.default(design(object), data=as.data.frame(colData(object)))來得到 ModelMatrix
實際運行getModelMatrix()我們得到的 modelMatrix其實是一個 1215*2 的矩陣。
下圖所示，第二列samplecancer就是分組信息，其中 0 代表正常，1 代表癌。
原始的metadata 中的樣本數量如下
> table(metadata$sample)
normal cancer
113 1102
ModelMatrix 中第二列的樣本情況
> table(modelMatrix[,2])
0 1
113 1102
而 nOrMoreInCell函數的作用是返回一個每個樣本一一對應的邏輯向量，來表明是否需要被篩選掉。這裡又需要提到一個參數minReplicatesForReplace。
這個參數，可能大多數人都不會被注意到。他的含義是 the minimum number of replicates required in order to use replaceOutliers on a sample.
那麼replaceOutliers()又是做什麼的呢？我們說 DEseq 在進行標準化的時候會丟掉一些 gene counts 數非常異常的值，轉而使用所有樣本中的 trimmed mean 來進行替代。replaceOutliers()就是做這件事的。
書接上文，minReplicatesForReplace 就指定了在對一個樣本使用 replaceOutliers 的時候需要最小的重複數是多少。假如這個值是 5，如果你的一個樣本只有 3 個重複，那麼就不會執行上面的操作。這個參數默認是 7，假設你不想做這一步（例如針對單細胞數據），你就需要把這個參數設置為 Inf。在實際運行中，如何判斷哪個樣品該進行哪些樣品不該進行這些操作呢？就是進行如下一個比較，下文代碼的第一行也展示了nOrMoreInCell()的作用。
sufficientReps <- any(nOrMoreInCell(attr(object,"modelMatrix"),minReplicatesForReplace))
if (sufficientReps) {
object <- refitWithoutOutliers(object, test=test, betaPrior=betaPrior,
full=full, reduced=reduced, quiet=quiet,
minReplicatesForReplace=minReplicatesForReplace,
modelMatrix=modelMatrix,
modelMatrixType=modelMatrixType)
}
扯遠了，再說回來，nOrMoreInCell() 是如何在新版本中提速的。
nOrMoreInCell_old <- function(modelMatrix, n) {
numEqual <- sapply(seq_len(nrow(modelMatrix)), function(i) {
modelMatrixDiff <- t(t(modelMatrix) - modelMatrix[i,])
sum(apply(modelMatrixDiff, 1, function(row) all(row == 0)))
})
numEqual >= n
}
到這裡很可能你已經看不下去了，但是其實不難理解。
原始函數使用sapply在每個樣品上都要迭代一次，因此如果的設計矩陣很簡單就會出現大量的重複工作。對於大量樣本這個函數佔用 DESeq（）相當一大部分運行時間。對於 1200 個樣本的數據來說，這個過程要重複 1200 次才能得到結果。
如果變成新的函數會有什麼效果呢？
nOrMoreInCell_new <- function(modelMatrix, n){
numEqual <- rep(NA, nrow(modelMatrix))
for(idx in seq_len(nrow(modelMatrix))){
if(is.na(numEqual[idx])){
modelMatrixDiff <- t(t(modelMatrix) - modelMatrix[idx,])
equal_to_idx <- apply(modelMatrixDiff, 1, function(row) all(row == 0))
numEqual[equal_to_idx] <- sum(equal_to_idx)
}
}
numEqual >= n
}
首先給 numEqual 全部複製為NA；當進行第一次循環idx=1時，因為numEqual[idx]是NA，所以一定會執行一次運算。
但是只要進行一次運算，equal_to_idx 就會返回一個新的邏輯值，所有equal_to_idx為真的位置對應的numEqual的位置都會被賦值這個樣本重複的數量，也就是sum(equal_to_idx)。然後第二次循環開始，對於和第一次循環同樣組的樣本來說，只需要進行if(is.na(numEqual[idx])) 的判斷，因為必定會返回False 就不需要再進行後續運算了。
以果子老師的數據為例， 他只有兩組數據，1215 個樣本。只需要進行一次運算，完成相當於原始函數的 1102 次計算（癌有 1102 個），然後循環到下一個是 NA 的樣本時（也就是第一個正常樣本時），再進行一次運算，就完成了原始函數的 113 次計算。
也就是當面對 2000 個樣本的 2 組數據時，原始函數需要運行 2000 次，而新的函數只需要運行 2 次。
在 1215 個樣本的 2 組數據中，前者需要 3.17s，而後者只需要 0.03s。
> t1 <- proc.time()
> samplesForCooks <- nOrMoreInCell_new(modelMatrix,3)
> proc.time() - t1
user system elapsed
0.00 0.02 0.03
>
> t1 <- proc.time()
> samplesForCooks <- nOrMoreInCell_old(modelMatrix,3)
> proc.time() - t1
user system elapsed
3.06 0.08 3.17
當然，直接導致時間複雜度改變的升級是在核心的 C 腳本中，我也只能勉強讀懂一些。
原始的 DESeq2.cpp中有一步需求是提取對角線矩陣。原代碼首先構造了一個 n 行 n 列的矩陣，然後再對整個矩陣進行運算，最後再用diagvec()來提取對角線矩陣，這個思路沒什麼問題，但是隨著樣本數量的增加，對於矩陣的運行時間是成樣本數的平方增加的，遇到大量樣本直接歇菜。
但在新的代碼中，避免了對於這個矩陣的計算。取代構造矩陣而是首先構造了一個 n 維向量，接下來不需要對全部矩陣進行計算就也可以得到對角線矩陣。這一步上的時間複雜度，就從之前的 O(ncol(Y)^2) 變成了 O(ncol(Y))的線性複雜度。這也是我們看到文章開始運行時間變化的主要原因。
進一步了解Bias, robustness and scalability in single-cell differential expression analysis
Need for Speed — DESeq2
果子補充這次更新是意義非凡的，對於那些想要從事數據挖掘的朋友來說，硬體是個難題，線下課上的時候，我說，這裡的Deseq2是檢驗大家電腦能力的時刻，差的電腦根本跑不過去，跑的過去可能需要20個小時，但是現在這些都不是問題了。我們現在提供文中提到的測試數據dds，你可以下載後在自己的電腦上測試一下速度和電腦的極限，如果完成不了就家內存，只要能運行這個數據，你的電腦就可以運行其他數據，暫時就不要糾結硬體的問題了。
在果子學生信公眾號回復"果子Deseq2"自助獲取，下載後放到當前的工作目錄下面，運行以下代碼就可以跟帖子中的教材無縫對接了。
load(file="dds_very_long.Rdata")
關於Deseq2的上下遊分析，可以看看這一個
很有誠意！人人可做的轉錄組數據下遊分析
如果什麼都搞不懂，也沒關係
回復"果子學生信"，獲取線下課程的R語言環境配置視頻教程，裡面還有GEO的視頻導學教程，以及大家喜歡的平臺和R包對應文件platformmap
如果對於GEO非編碼RNA的注釋搞不清楚，
回復"果子非編碼"獲取教程和整理好的文件。
如果對於TCGA的數據挖掘感興趣，
回復"果子TCGA"獲取持續更新的系列教程。