Kharchenko, Peter V. The triumphs, and limitations of computational methods for scRNA-seq. Nature Methods (2021): 1-10. PMID: 34155396

各位親愛的小伙伴們，本期又是小師弟來和你們嘮一嘮生信的文章了，有沒有誰想念我了呢？哈哈，不知道大家在科研中是不是經(jīng)常接觸單細(xì)胞轉(zhuǎn)錄組數(shù)據(jù)呢？相信大伙兒都能找到各式各樣的教程文章帶你一步一步地做單細(xì)胞分析，然后跑出各式各樣漂亮的圖。但是你知道這背后每一步可能碰到的陷阱和障礙嗎？做科研要求我們不僅知其然還要知其所以然，知己知彼方能百戰(zhàn)百勝。近期這篇發(fā)表于nature methods上的文章就以review的形式探討了單細(xì)胞轉(zhuǎn)錄組數(shù)據(jù)分析中每一步所蘊(yùn)含的假設(shè)以及可能會遇到的陷阱。小師弟覺得十分有用，因此特地寫了這篇推送文章來和大家一起解讀一下這篇文章的要點（強(qiáng)烈建議小伙伴們訪問原文地址進(jìn)行詳細(xì)品讀）。為了方便大家快速了解本文內(nèi)容以及選讀感興趣的部分，我在摘要之后以圖文形式列出了本文的要點。

單細(xì)胞轉(zhuǎn)錄組測序技術(shù)在過去數(shù)十年間得到了飛速發(fā)展，與此相伴而生的是各類計算生物學(xué)方法的大量涌現(xiàn)。隨著檢測技術(shù)的精度和通量的不斷提升，生物學(xué)各方面的復(fù)雜性在新算法的幫助下被逐步揭示，這些復(fù)雜性具體體現(xiàn)在細(xì)胞群體組成、基因調(diào)控以及細(xì)胞動態(tài)演變等各方面上。與此同時，計算生物學(xué)的快速發(fā)展使得人們需重新評估各種算法所依賴的統(tǒng)計模型，所要實現(xiàn)的實驗?zāi)康?，以及所要處理的巨大?shù)據(jù)量。在此篇文章中，作者回顧了單細(xì)胞轉(zhuǎn)錄組數(shù)據(jù)分析的關(guān)鍵步驟，總結(jié)了不同方法所依賴的假設(shè)，并指出這些方法的成功之處和仍然存在的不足之處。由于單細(xì)胞轉(zhuǎn)錄組測序技術(shù)已成為生物學(xué)研究的主流技術(shù)之一，因此本文的討論對于廣大科研人員具有極大的參考意義。

1. 如何以統(tǒng)計的視角看待單細(xì)胞數(shù)據(jù)

單細(xì)胞轉(zhuǎn)錄組數(shù)據(jù)由于其技術(shù)的限制，其所存在的不確定性要大于常規(guī)RNA-seq數(shù)據(jù)，其中的原因主要包括兩點：1）每個細(xì)胞的總RNA量僅有一小部分能被儀器所檢測到；2）單細(xì)胞之間被捕獲的RNA量存在巨大差異。因而測序技術(shù)通常通過提高單細(xì)胞的捕獲數(shù)量來抵消隨機(jī)性。

另一方面，研究者們也嘗試通過概率模型來解釋觀測到的隨機(jī)性。很多研究采用了負(fù)二項分布模型來對基因表達(dá)的分布進(jìn)行建模。需要注意的是，許多研究認(rèn)為隨機(jī)性導(dǎo)致單細(xì)胞數(shù)據(jù)包含了過多的零值，因而這些研究的模型通常包含了一個零膨脹(zero-inflation)的部分來解釋過多的零值。但一項既往研究的實驗結(jié)果顯示：使用UMI(uniquemolecular identifier)技術(shù)的單細(xì)胞數(shù)據(jù)已能夠極大地緩解這個現(xiàn)象[1]，已沒有必要使用包含零膨脹的模型。

2. 比較轉(zhuǎn)錄狀態(tài)

在建立了模型之后，下一步就是利用模型對基因的轉(zhuǎn)錄狀態(tài)進(jìn)行比較，這就是常見的差異表達(dá)分析，其目的是比較一個基因在兩個細(xì)胞群體間表達(dá)量的差異。通過統(tǒng)計模型可知基因在兩個細(xì)胞群體中的分布形狀，那么兩個分布之間的重疊部分就是基因表達(dá)在兩個群體中相同的概率（見下圖），可理解為統(tǒng)計檢驗中的P值（紅字為個人觀點，有待商榷）。在一些研究中，人們使用參數(shù)模型（比如負(fù)二項分布，高斯分布等）來完成分布的估計和差異表達(dá)分析，然而當(dāng)單細(xì)胞數(shù)量增多后，參數(shù)模型不僅計算緩慢并且效果不佳。因而在許多研究中，研究者傾向于采用計算速度較快的非參數(shù)分布來完成差異表達(dá)分析(例如 Wilcoxon秩和檢驗)。不過需要注意的是，如果細(xì)胞群體足夠大，幾乎所有的基因都會呈現(xiàn)統(tǒng)計學(xué)差異(見下圖)，因此引入其他標(biāo)準(zhǔn)來輔助差異基因的篩選就變得十分重要(例如設(shè)定一個fold-change閾值，使用更小的細(xì)胞亞群進(jìn)行差異表達(dá)分析等)。

除了進(jìn)行差異表達(dá)分析外，比較轉(zhuǎn)錄狀態(tài)還包括了計算細(xì)胞與細(xì)胞間的距離或相似性（此步驟是下游其他分析諸如細(xì)胞聚類，可視化以及軌跡分析的基礎(chǔ)）。若將單細(xì)胞的每個基因視作一個維度，那么毫無疑問單細(xì)胞數(shù)據(jù)是一種高維數(shù)據(jù)，在高維空間中計算細(xì)胞間的距離將不可避免的面臨“維數(shù)災(zāi)難”(curse of dimensionality)的問題。簡而言之，在高維空間中，常用的距離度量方式(歐式距離)將無法區(qū)分距離近的細(xì)胞和距離遠(yuǎn)的細(xì)胞。所幸的是，可通過降維的方式來緩解這一問題，下一小節(jié)將會集中探討單細(xì)胞數(shù)據(jù)降維的問題。

3. 高維度數(shù)據(jù)降維

單細(xì)胞數(shù)據(jù)的有效維度數(shù)量要比其原始維度低得多。何為有效維度？有效維度是表征單細(xì)胞數(shù)據(jù)變化程度所需的最低維度。從生物學(xué)意義上也可以很好地理解這個問題。基因組中的許多基因通常存在功能的高度相關(guān)性，因此原始的高維數(shù)據(jù)冗余程度極高，而降維則等同于尋找最簡約的表征數(shù)據(jù)的方式。

第二大特點是單細(xì)胞數(shù)據(jù)的稀疏性。由于PCA的原理為求協(xié)方差矩陣的矩陣分解，因此只有當(dāng)數(shù)據(jù)分布呈對稱形狀時PCA才能達(dá)到最佳效果。然而，如前所述，單細(xì)胞數(shù)據(jù)含有很多零值，導(dǎo)致數(shù)據(jù)分布形狀向零值方向偏斜（見下圖），整個分布呈極為不對稱的形狀。針對此，既往文獻(xiàn)中采用了一些校正方式，可供參考[2-4]。

PCA之所以無法克服以上兩大困難，原因在于其計算矩陣分解的過程中沒有考慮到單細(xì)胞數(shù)據(jù)自身的統(tǒng)計分布。針對此，既往研究采用的思路為把數(shù)據(jù)分布建模和矩陣分解結(jié)合到一起，例如ZINB-WaVE就是其中的代表之一，詳情可見[5]。

PCA歸根結(jié)底仍只是一種高維到低維空間的線性投影(linear projection)。例如照片就是一種高維到低維的線性投影方式，把三維空間的物體投影到維度更低的二維空間（在此小師弟又要安利一下三體小說，關(guān)于維度的描述實在太棒）。然而投影方式不一定要局限于線性投影，非線性投影很多時候能夠更好地反應(yīng)數(shù)據(jù)本身的結(jié)構(gòu)，現(xiàn)今有關(guān)單細(xì)胞分析的文章中常用的非線性投影方式有UMAP、t-SNE以及autoender等。(小師弟個人觀點，UMAP速度快，同時能夠保留遠(yuǎn)距離細(xì)胞群體間的信息，對下游分析諸如軌跡分析和細(xì)胞聚類等很有幫助)。

4. 利用最近鄰圖構(gòu)建流形

在先前的小節(jié)里已經(jīng)提到了高維數(shù)據(jù)存在極大冗余性，因此需要采用降維的手段來分析數(shù)據(jù)。除此之外，我們還可以通過流形在低維空間中表征高維空間的數(shù)據(jù)。何為流形？流形簡而言之就是一種低維空間的幾何表面，其所代表的區(qū)域能夠逼近單細(xì)胞數(shù)據(jù)的分布。現(xiàn)今學(xué)術(shù)界普遍采用的方式是用最近鄰圖(nearest-neighbor graph)來估計流形的形狀。最近鄰圖的構(gòu)建簡單明了，只需通過一定手段計算細(xì)胞與細(xì)胞間的距離，再將距離作為邊，細(xì)胞作為節(jié)點即可構(gòu)建最近鄰圖（也稱為最近鄰“網(wǎng)絡(luò)”，見下圖）。最近鄰圖在單細(xì)胞數(shù)據(jù)分析中的應(yīng)用十分廣泛，前一小節(jié)提到的UMAP和t-SNE降維方法實際上就是最近鄰圖在二維空間上的可視化，與PCA相比，它們不受表達(dá)值量級和分布形狀的影響，可以很好地表示高維空間中細(xì)胞與細(xì)胞間的局部距離和數(shù)據(jù)的總體結(jié)構(gòu)。

5. 細(xì)胞聚類

細(xì)胞聚類也是單細(xì)胞數(shù)據(jù)分析中一個重要步驟。聚類的結(jié)果是將細(xì)胞分為了不同的群體，以方便下游的進(jìn)一步分析，比如在群體之間進(jìn)行差異表達(dá)分析。通常來說，其依據(jù)的基礎(chǔ)是利用最近鄰圖來鑒定不同的細(xì)胞群體，使得群體中細(xì)胞與細(xì)胞的聯(lián)系比群體間細(xì)胞與細(xì)胞的聯(lián)系更為緊密。常用的聚類算法有Louvain聚類和Leiden聚類等。不過，這些算法所依據(jù)的只是細(xì)胞轉(zhuǎn)錄組數(shù)據(jù)的相似性，實則并沒有考慮到真正的生物學(xué)意義。因此，對聚類結(jié)果的解釋一定要結(jié)合實際生物學(xué)問題。這使得在實踐中，對聚類的分辨率大小的選擇變得十分重要。例如在有些研究中需要將T淋巴細(xì)胞視作一個大類群而另一些研究中則需要將T細(xì)胞細(xì)分為不同亞群，這時候就需要給予聚類算法不同的分辨率等級。一些層次聚類算法，例如walkTrap[6]就能很好捕捉到不同分辨率等級下的細(xì)胞類群，為分析提供了靈活性。

6. 細(xì)胞的動態(tài)演變過程

單細(xì)胞轉(zhuǎn)錄組數(shù)據(jù)是探究器官組織發(fā)育以及細(xì)胞類型分化演變的利器之一。現(xiàn)有的轉(zhuǎn)錄組測序技術(shù)可對細(xì)胞群體在某個時間點的狀態(tài)拍一個“快照”。在對原始數(shù)據(jù)降維之后，在低維空間中可通過這個靜態(tài)的“快照”推斷細(xì)胞群體的演變軌跡，這就是常見的單細(xì)胞軌跡分析。Monocle是目前基于這一原理最常被使用的軌跡分析工具(最新版本為Monocle3)[7]。需要注意的是，并非所有的細(xì)胞都參與了演變或分化的過程，但是此類軌跡分析工具通常將所有細(xì)胞都作為輸入用于構(gòu)建軌跡。

此外，此類軌跡分析還隱含一大假設(shè)，即只要測序得到的細(xì)胞樣本數(shù)量足夠大，那么單一時間點所捕獲的細(xì)胞“快照”能夠涵蓋細(xì)胞所有可能的狀態(tài)。但是此假設(shè)對于某些生物學(xué)過程來說顯然有誤。例如，細(xì)胞群體對外部刺激的響應(yīng)過程就無法使用單一時間點的“快照”衡量。雖然有一些軌跡分析工具可以結(jié)合多個時間點的單細(xì)胞數(shù)據(jù)[8,9]，但大部分主流工具仍然只能針對單一時間點進(jìn)行軌跡分析。

最后，軌跡分析的一大缺陷還在于無法給出細(xì)胞群體演變的方向。于是本文作者在這個方向上做了努力，在一項已發(fā)表的工作中提出了RNA速度(RNA velocity)這一概念。簡單來說，RNA速度衡量的是前體RNA和成熟mRNA之間的比例，如比例達(dá)到平衡態(tài)則基因表達(dá)在下一時刻可能不發(fā)生變化；如果前體RNA數(shù)量遠(yuǎn)大于成熟mRNA，則下一時刻基因有可能上調(diào)表達(dá)，反之則下調(diào)。利用這一特性可以預(yù)測該時間點每個細(xì)胞在下一時刻的狀態(tài)，進(jìn)而推斷細(xì)胞群體演變的方向。（小師弟云：RNA速度也是一個很酷的概念，相應(yīng)的分析工具叫做velocyto，大家可以在這里下載安裝嘗試一下：http:///。另外一個更新的基于RNA速度的工具叫scVelo，可從這里獲得：https://scvelo.）。

7. 單細(xì)胞轉(zhuǎn)錄組分析在未來將走向何方？

目前為止，單細(xì)胞轉(zhuǎn)錄組數(shù)據(jù)分析工具的主要努力方向在于解決數(shù)據(jù)的高噪聲和高稀疏性。今后，隨著測序技術(shù)的不斷完善，研究的重心也會轉(zhuǎn)移到對其他類型分子的測量和分析上：如DNA甲基化、染色質(zhì)可及性以及蛋白質(zhì)豐度等。此外，空間轉(zhuǎn)錄組也是未來的一大熱門方向。如何整合這些不同的數(shù)據(jù)類型將會是計算生物學(xué)領(lǐng)域所要面對的挑戰(zhàn)。

本文是單細(xì)胞轉(zhuǎn)錄組分析領(lǐng)域的一篇極具指導(dǎo)意義的綜述文章，基本涵蓋了單細(xì)胞數(shù)據(jù)分析的所有重要步驟以及每一步蘊(yùn)含的假設(shè)和可能踩到的“坑”，值得一讀再讀。同時，文中列出的許多分析工具對于廣大科研人員的數(shù)據(jù)分析也是相當(dāng)實用，大家可以參考我文末附上的參考文獻(xiàn)。當(dāng)然有時間的同學(xué)可以品讀一下原文，由于推文篇幅所限，無法一一展現(xiàn)其中精華內(nèi)容。

以上就是本期全部的推送內(nèi)容，我是愛你們的小師弟，我們下期再見啦！

1.Svensson,Valentine. Droplet scRNA-seq is not zero-inflated. NatureBiotechnology 38.2 (2020): 147-150.

2.Fan, Jean, et al. Characterizing transcriptional heterogeneity through pathway and gene set overdispersion analysis. Nature methods 13.3 (2016): 241-244.

3.Eling, Nils, et al. Correcting the mean-variance dependency for differential variability testing using single-cell RNA sequencing data. Cell Systems 7.3 (2018): 284-294.

4.Hafemeister, Christoph, and Rahul Satija. Normalization and variance stabilization of single-cell RNA-seq data using regularized negative binomial regression. Genome Biology 20.1 (2019): 1-15.

5.Risso, Davide, et al. A general and flexible method for signal extraction from single-cell RNA-seq data. Nature communications9.1 (2018): 1-17.

6.Pons, Pascal, and Matthieu Latapy. Computing communities in large networks using random walks. International symposium on computer and information sciences. Springer, Berlin, Heidelberg, 2005.

7. Cao, Junyue, et al. The single-cell transcriptional landscape of mammalian organogenesis. Nature 566.7745 (2019):496-502.

8.Schiebinger, Geoffrey, et al. Optimal-transport analysis of single-cell gene expression identifies developmental trajectories in reprogramming. Cell 176.4 (2019): 928-943.

9.Tran, Thinh N., and Gary D. Bader. Tempora: cell trajectory inference using time-series single-cell RNA sequencing data. PLoS computational biology 16.9 (2020): e1008205.

注：詳細(xì)信息請參見原文！

https://www./articles/s41592-021-01171-x