隨著單細胞測序技術的流行，我們對複雜疾病和性狀的理解從patient，tissue的表達譜(bulk RNA-seq)到單個細胞的表達譜(single cell RNA-seq)。究其原因，在於bulk RNA-seq產生的是一個細胞群體的平均讀數，而細胞，特別是癌細胞存在極大的異質性，這些平均信號可能不足以反應這堆組織的真實信息。Prof Aviv Regev (MIT教授，HHMI研究院，人類細胞圖譜計劃項目co-chair)曾經形容這種方法就像水果沙拉，顏色和味道都能提示冰沙的成分，但倘若只有幾個是屬於藍莓的味道，那麼很容易就被一堆草莓的味道所覆蓋（如下圖）。因而在細胞尺度上進行大規模的測序分析以對細胞進行重新分型是很有必要的。

載入超時，點擊重試

圖源：http://fightdipg.org/wp-content/uploads/2015/12/single-cell-sequencing.jpg

然而，類似於傳統的在bulk RNA-seq中的分析，我們依然要對兩個不同的細胞類群cell type (group)進行差異表達分析（Differential expression analysis），尋找具有顯著性表達差異的基因，縮小找marker genes的範圍。但是基因在scRNA-seq data中的表達分布並非完全的負二項分布。

一個比較重要的區別是在scRNA-seq Data 中有更多的0值，即「稀疏性（sparse）」。這些0所包含的信息也不一樣，受到技術上和細胞本身的兩個因素的影響。從技術上講，大多數mRNA在單個細胞的量是非常少的，而絕大多數基因在一個細胞中只有少量的mRNA的拷貝，由於在mRNA capture protocol的一些問題，可能會導致一些基因的mRNA在逆轉錄和cDNA擴增階段存在完全的缺失，進而導致漏測，這種現象被稱為dropout event。而也存在一些基因的表達為0的情況是由於其本身在細胞的表達量就非常低，或者由於mRNA bursting的現象（文獻可參考：Golding, I; Paulsson, J; Zawilski, SM; Cox, EC (2005). "Real-time kinetics of gene activity in individual bacteria". Cell. 123 (6): 1025–36. doi:10.1016/j.cell.2005.09.031. PMID 16360033）。那麼對應的，受這些0的影響，原本的負二項分布（negative binomal distribution）的基因可能更多服從於0膨脹負二項分布（zero-inflated NB）。

基於scRNA-seq數據的特殊性(稀疏，高雜訊，高維度)，為bulk RNA-seq data設計的差異表達方法是否適合於scRNA-seq data也需要做進一步討論。但同時即便目前湧現了很多scRNA-seq data，這些data的分析目的不同，使用的pipelines也相差很大，因而重新使用這些已經經過預處理的公開數據，同時進行不同方法的比較也變得很困難。

我們今天主要介紹Charlotte Soneson等於2018年在《Nature Methods》發表的「Bias, robustness and scalability in single-cell differential expression analysis」這篇文章。這項研究利用相同的pipelines處理，整合了多個公開的scRNA-seq數據，同時對所有的注釋基因和轉錄本進行了丰度估計，以及質量評估和研究分析報告提供給使用者以確定他們所需要的合適的資料庫。

Conquer為文章作者開發出的一個整合性資料庫。目前總共有36個數據集：31個來自於full-length protocols同時5個是3』-end sequencing(UMI) protocol。Conquer對這些數據進行一致的預處理和分析，有利於不同計算方法的比較。同時Conquer還包含了不同的物種和細胞類型可以檢驗某些生物學假說的泛化性(generality)，也就是假說的適用廣度和範圍。

分析方法/Methods

在這篇文章中，共選擇了七個數據集(six full-length and one UMI data)以及兩個額外的UMI count數據對如下的差異分析方法進行評估。

這些差異表達分析方法包括了對bulk RNA-seq data的經典分析方法 (SAMseq, edgeR, limma, DESeq2等), 專門開發出來用於scRNA-seq data的分析方法 (SCDE, DEsingle，Seurat，MAST等)以及一些常見的統計方法 (Wilcoxon, t-test等)

數據準備/data preparation：

Signal data:

我們將之前選擇的九個資料庫只保留已經注釋好的兩個cell group，從t-SNE的結果看，在低維流形上絕大多數數據對應的兩個cell group在表達上具有明顯差異。對於每個數據集，我們生成一個『maximal』大小的實例 (舉例：即每個group的細胞數量等於原始的較小那個cell group的細胞數量)，以及從maximal大小的數據集中隨機採樣(random sampling) 生成若干個小的數據集，這些數據集都作為signal data。

null data:

null data用於後續的type 1 error分析，是在signal data中單個細胞群體中隨機採樣得到新的數據作為null data。在這些null data中理論上不存在差異表達基因，因為都是來源於一個細胞群體。

Simulation data：

基於powsim R package（鏈接： https://github.com/bvieth/powsimR）和三個原始數據集生成simulation data進行分析(表示不太懂powsim這個包，不過貌似是生成和原始數據在分布和dispersion等等參數上類似的模擬數據，但選擇了特定的10%的基因為差異表達基因)，這些simulation data是用於後續的演算法性能在TPR和AUROC等等方面進行評估的。

Number of differentially expressed and non-tested genes

文章對九個scRNA-seq data對應的signal data計算和統計了差異表達基因數量。對於full-length數據集而言，在沒有做prefiltering （prefiltering為過度掉低表達基因）的情況下，SeuratBimod檢測到了最多存在顯著表達差異的基因，而經過prefiltering以後，大多數方法的檢測到的差異表達基因數量有所下降，edgeR/QLF檢測到的差異表達基因數量下降是最多的。且隨著兩個group細胞數量的上升，發現的差異表達基因數量越多。

大多數差異表達方法並不是在全部基因上都做了test，除了方法內置有filtering步驟，一些基因也會運行報錯，可能是在擬合的步驟無法收斂而去除的一些基因，這在bulk RNA-seq。

（觀察右邊的圖可以發現相比於左邊的沒有filtering，返回NA的基因數量明顯降低了)

Type 1 error control

Type 1 error也是很重要的，其實這一部分也就是看false positive rate(FPR)在各個數據集各個方法上的情況，對於8個原始的數據中通過對單個細胞組(cell group)進行隨機採樣得到一個null data，那麼在null data中理論上不存在差異表達基因(因為都是取自同一種細胞)。因而我們希望在這些data上運行差異表達分析方法以後返回後的差異表達基因佔全部檢測的基因數量的佔比要小於0.05(對應的就是type 1 error也就是FPR)。

對於unfiltered的數據集而言，ROTS以及SeuratTobit的性能是最好的。但對於edgeR/QLF，SeuratBimod兩種方法我們可以看到結果有很多假陽性，經過filtering以後會好很多，但是會少很多檢測顯著的基因。特別是SeuratBimodIsExpr2返回的DE的基因數量就不是很多

從UMI數據集看，monocle使用TPM數據反而效果不佳，當在full-length數據中，TPM數據的效果比read count要更好，特別是篩選一些低表達基因以後，相比於其他方法， voom-limma，ROTSvoom以及edgeR/QLF這些基於線性模型的bulk 分析方法的性能有了很大提升同時結果更加穩定。

除了FPR，從p value的分布上也有一些很有意思的現象，在沒有gene prefiltering的時候，每個方法在p value的分布上並不是均勻的，特別是在SeuratBimod中p value較小的基因占絕大多數，相反的monocle中p value較大的基因占絕大多數。其他基因也存在類似的情況，但經過prefiltering以後都有所改善(趨向於均勻分布)，特別的，即便有prefiltering，SeuratBimod依然傾向於給出較小的p value，考慮到之前FPR的結果，SeuratBimod的假陽性就顯得有點過高了。

Characteristics of false-positive genes

那麼這些false-positive genes和其他True-negative genes在表達的特徵上有什麼區別？這個研究主要估計了各個基因的四個表徵(characteristic)：平均表達量(average of CPM)，方差(variance of CPM)，變異係數(variance and coefficient of variantion of CPM)以及表達量為0的細胞在所有細胞中的佔比(Fraction Zeros)。同時計算了信噪比統計量進行比較(signal-to-noise，其實可以理解為scaling，方便不同方法之間進行比較)。結果發現對於NODEs，ROTS，SAMseq，SeuratBimod這些方法而言，假陽性的基因表達量為0的情況比較少，同時表達量較高，變異係數較小。相反的，對於edgeR/QLF，SeuratTobit，MAST以及metagenomeSeq有較多的表達量為0的值，平均表達量較小，變異係數較大。

Between-method similarity

對於一個差異表達分析方法而言，可重複性是很重要的，因而每次分析出的差異基因的一致性(concordance)不應受到樣本量的影響。或者說，如果對於兩個population中有真實的潛在的差異信號，不管取這兩個population中多大的樣本量我們都可以發現這個signal。因而假如我們隨機抽取樣本做差異表達分析，每次做出的結果應該有很高的一致性。相反的，對於兩個population中不存在差異表達信號的(null data)，我們每次做出的結果理應是隨機的，不一致的，倘若存在一致性。那麼這個數據集或者方法或者說兩者的組合就會存在一定的bias。

因而基於此，為了測量（measure）這種一致性，我們可以依據顯著性(pvalue)進行rank，我們選擇top100個具有顯著性差異的基因，在每一對具有相同樣本量的數據實例(data instance)的比較中，我們看top1-100有多少個基因是一致的，比如Top20個基因中有五個基因是一樣的，那麼就記為5。從橫坐標1-100對應的作出縱坐標相應的一致基因的數量，就可以得到一條一致性曲線(concordance curve)，類似於AUC之於ROC，我們計算AUCC，也就是一致曲線的下面積。同時除以100^2/2進行標準化(100^2/2相當於從top1-100都是一致的)。

從結果上看，我們可以看到顯著的signal data中存在較高的similarity(concordance)，相反的null data的一致性就不是很高了。一個良好的數據集在signal data與null data中一致性差異應該是比較大的，而EMTAB2805這個數據集signal data與null data差異卻不是很高，而10XMonoCytoT的差異就很明顯了。有意思的是，一些無參數方法(Wilcoxon test,SAMseq)以及一些log-like transformation處理的方法(limma-trend, ROTSvoom, voom-limma)相比於基於count數據分析的方法在signal data和null data之間的差異更大，說明在魯棒性上，對於scRNA-seq數據的一些小的變化，基於count數據可能更加敏感。

利用AUCC，我們也可以比較不同方法之間結果的相似性，從聚類的結果可以看出來，依據數據的輸入格式，transformation的方式以及modeling的方式，都會對方法的一致性產生影響。但同時需要指出的是，如果在進行一致性比較時兩個實例中的細胞數量與AUCC值存在一定的相關性。

FDR control and power

最後，這項研究也基於simulated data比較了FDR control，TPR(true positive rate)以及AUC的情況。一些方法如基於Census data的voom/limma，ROTSpm，MAST，以及有prefiltering的SeuratTobit，SeuratBimod和SAMseq都能控制FDR在0.05這個水平左右，而SCDE，scDD和t test，limma-trend以及wilcoxon test等等方法可以控制FDR低於0.05這個水平。edgeR/QLF的FDR水平較高，經過filtering以後，edgeR/QLF的FDR顯著下降。同時大多數方法提高樣本量都能降低FDR的水平。

從真陽性率TPR上去看，經過filtering後，如edgeR/QLF，SAMseq，DEsingle還有voom-limma均有良好的表現。DESeq2在不管有無intern filtering後均有良好的表現，同時從AUROC上去看的話，最好的是edgeR，MAST，limma以及SCDE等等，經過prefiltering後，大多數方法的AUROC值都有了相對的提高。文章最後還比較了各個演算法的運行速度和這些方法的開源性以及documentation，不作詳細的介紹。

結論

從這篇研究我們可以看到，過去一些研究認為的在bulk RNA-seq中分析的方法不太適用於scRNA-seq數據，究其原因在於scRNA-seq由於Drop-out event以及單個細胞的transcript太少，估計出的單個基因在多個細胞中的表達水平會存在大量的0值，因而這時候的基因應該服從0膨脹的負二項分布分布(zero-inflated negative binomial distribution，可以理解為在0這一點的狄拉克分布和負二項分布的混合分布)，不同於在bulk RNA-seq中的負二項分布。

而這一研究通過對不同方法在不同數據以及不同輸入格式(CPM, TPM, Census Count等)進行比較，發現prefiltering確實會顯著影響到DE演算法的性能，同時傳統的一些limma，edgeR在scRNA-seq data中一樣能有很好的表現。但值得注意的是，如edgeR和limma對一些0值較多的表達量比較低的基因檢驗顯著，因而在prefiltering這一步相比於bulk RNA-seq要更加嚴格。

其實懂得No Free Lunch theory（1、一種演算法（演算法A）在特定數據集上的表現優於另一種演算法（演算法B）的同時，一定伴隨著演算法A在另外某一個特定的數據集上有著不如演算法B的表現；2、具體問題具體分析。）便會知道，即便我們看到了edgeR和limma這些在bulk RNA-seq中較好的軟體可能在scRNA-seq一樣也適用，但是依然得從數據的分布出發，假如我們手頭的scRNA-seq數據絕大多數基因表達都比較稀疏，其實使用有prefiltering的一些方法可能會錯過很多signal。因而我們可以畫出 scRNA-seq data整體表達情況的分布圖，同時從分布上去看，是否服從0膨脹分布的基因的佔比也會影響到演算法的robust也是值得討論的。

除了在方法的選擇上看，對於今後凡是做演算法性能比較的real data（在做出新的演算法時為了和過去演算法進行比較，需要真實的數據進行性能比較），這項研究也給出了很有意義的參考。我們可以看到一些scRNA-seq data分析出的差異表達基因的可重複性容易受到比較的時候細胞數量的影響，而選擇更穩定的，受隨機因素影響更小的scRNA-seq data也許在與其他演算法進行比較時會更好，比如10XMonoCytoT。為了給用戶提供適合他們使用目的的scRNA-seq數據，conquer還提供了探索性分析報告。該報告計算並可視化細胞的各種質量測量，並提供細胞的低維表示，由不同的表型注釋著色等等信息。

(http://imlspenticton.uzh.ch:3838/conquer/)

同時，在單細胞測序數據分析中，不僅僅在差異表達分析上存在一些需要討論的問題，在做依據表達譜進行細胞分型也是需要做進一步標準化資料庫，方法等等，放在同一個框架下討論各個方法的適用範圍。之前聽的一個seminar上就說目前幾十種聚類演算法，能找到幾十種不同的分型結果，富集到各種不同的pathways，究竟哪一種更make sense；以及，在bulk RNA seq data中用於聚類的經典演算法(PCA, K-means等)能否一樣適用於scRNA-seq data，可能也要像這項研究一樣進行更加細緻的討論。

參考文獻：

[1] :Soneson C, Robinson M D. Bias, robustness and scalability in single-cell differential expression analysis[J]. Nature Methods, 2018.

[2] Miao Z, Deng K, Wang X, et al. DEsingle for detecting three types of differential expression in single-cell RNA-seq data.[J]. Bioinformatics, 2018.

推薦閱讀：

相关文章