單細胞RNA測序數據分析方法研究進展

張淼 孫祥瑞 徐春明

（北京工商大學輕工科學技術學院，北京 100048）

隨著分子生物學技術的不斷發展，高通量測序已廣泛應用于臨床醫學、基礎醫學、生物醫學等眾多領域，因其具有通量高、速度快、靈敏度高等優勢，可在短時間內檢測大量樣本的基因變異及轉錄水平，具有十分廣闊的應用前景。以二代測序（Next generation sequencing，NGS）為典型代表的高通量測序已經廣泛應用于多種疾病的診斷、治療及預后評估［1］。

NGS可分為DNA測序和RNA測序，DNA測序以檢測基因變異為主［2］，如堿基替換、小片段插入及缺失等，RNA測序主要以檢測基因的mRNA豐度為主［3］。根據所測細胞的群體區分，RNA測序可分為全轉錄組測序（Bulk RNA sequencing，bulk RNA-Seq）和單細胞轉錄組測序（Single cell RNA sequencing，scRNA-Seq）。Bulk RNA-Seq是目前最常用的RNA測序方法，但由于測序樣本普遍具有異質性，bulk RNA-Seq結果僅能代表大細胞群體中每一個基因的平均表達水平，對比較轉錄組學研究有幫助，但不利于異質性研究。2009年scRNA-Seq技術［4］被首次引用，測序樣本異質性問題得到一定程度解決，并逐漸成為近年來的研究熱點。scRNA-Seq原理是分離單個細胞并提取RNA，經PCR擴增后進行高通量測序，主要工作流程包括細胞解離、單細胞分離、文庫構建、上機測序及數據分析［3］，其中數據分析是整個scRNA-Seq過程中最重要的一個環節。

由于scRNA-Seq所使用的數據分析方法有別于bulk RNA-Seq，越來越多的針對scRNA-Seq技術的數據分析方法不斷涌現，但每種分析方法都有各自的優勢及局限性。因此，本文對比了scRNA-Seq與bulk RNA-Seq技術在數據分析上的差異，對scRNASeq數據分析研究進展進行總結，探討每種方法的優勢與局限性，以期能夠對scRNA-Seq數據分析方法進行系統了解。

1 數據預處理

1.1 數據的比對

在scRNA-Seq過程中，測序數據下機后需要進行預處理，將基因轉錄序列轉換為fastq格式并與參考序列比對，鑒定差異表達基因或尋找可變剪切位點。scRNA-Seq數據整體質量評估的重要指標是比對率，即reads在參考基因組中所占的比例。比對率越高說明數據利用率越高［5］。目前應用于數據比對的工具主要有TopHat［6］，STAR［7］或HISAT［8］。TopHat工具以Bowtie作為核心算法針對75 bp以上長度的RNA短序列與參考基因組進行比對，找到匹配的序列，對外顯子進行選擇性拼接，具有內存小、準確性高、容錯率低、可跨內含子比對等優勢。Donbin等［7］研究發現STAR工具的核心算法是Maximal mappable prefix（MMP），其直接選用非連續序列進行比對，運行速度較TopHat工具快，但需要更大的內存。HISAT工具基于Burrows-Wheeler變換（BWT）和Ferragina-Manzini索引（FM）結合的算法進行對齊。Kim等［8］研究發現HISAT是第一個采用分層索引以及自適應策略進行對齊的工具，減少了內存需求，也是目前運行速度最快的工具，具有同其他比對工具相同甚至更高的精度。從總體來看，以上3種軟件在運行速度和結果準確性方面均表現良好，但Engstr?m等［9］進一步研究發現運行速度較快的軟件通常檢測的準確性較低。

1.2 質量控制

數據識別和去除低質量細胞是scRNA-Seq質量控制（Quality control，QC）的關鍵步驟。首先，在細胞捕獲過程中應避免參雜混合或死亡的細胞。其次，使用FastQC［10］等工具檢查測序數據的質量，根據QC值決定其在后續分析中是否被舍棄。異常的cut-off值可以人為定義，也可由程序自動定義，但需考慮被分析組織的多樣性。此外，當細胞受損時，細胞質RNA會丟失，但線粒體RNA會保留在受損細胞中，線粒體RNA含量是QC的另一個指標［11］。除上述方法外，Jiang等［12］最新提出的一種針對單細胞RNA序列質量控制（SinQC）工具，通過整合基因表達模式和數據質量信息來檢測并去除低質量細胞。

1.3 數據標準化

數據標準化是數據預處理的關鍵步驟。在bulk RNA-Seq中，DESeq2［13］和TMM［14］是常用的標準化處理算法。然而，DEseq2算法并不適用于scRNA-Seq數據分析，因為DEseq2算法假設所有樣本RNA總量相等，reads數僅與測序深度有關，根據不同樣本的reads數來計算比例，但對于單個細胞可能會受到零值和高變異性的影響，結果不穩定。因此，Bacher等［15］提出SCnorm工具，其使用分位數回歸的方法對測序數據進行標準化，可以避免針對bulk RNA-Seq的傳統標準化方法對scRNA-Seq數據進行標準化時所引入的錯誤，改善主成分分析和差異表達基因的識別，可以用于scRNA-Seq數據標準化。最常用的數據歸一化方法為 Count depth scaling，又稱為Counts per million（CPM），它會根據每個細胞的總表達量計算一個size factor，然后對其中各個基因表達量進行標準化。除CPM外，非線性標準化方法使用測序深度擬合的負二項模型，可解釋更復雜的異質性。因此，數據標準化方法的使用需要根據不同細胞特性進行選擇，單一的標準化方法不能適用于所有類型的scRNA-seq數據。

2 插補

在測序過程中，很多低表達或中度表達的基因無法有效檢測到，導致表達值為零或減少，影響下游后續分析，增加細胞間變異，甚至不能獲得完整的單個細胞轉錄組信息，這種情況稱之為Dropout。在實際操作過程中Dropouts［16］對數據分析影響較大，合適的插補方法可以彌補Dropout產生的影響。目前針對插補開發的算法有MAGIC［17］、ScImpute［18］、SAVER［19］、DrImpute［20］和AutoImpute［21］，各種算法的計算原理不同。MAGIC算法使用基于Markov親和力的矩陣確定細胞間的相似性，對高度相似細胞中基因表達進行聚集，以估算基因表達量。ScImpute算法通過擬合Gamma-Normal混合模型估算基因缺失概率，根據相似細胞信息估算可能的Dropout。SAVER和MAGIC算法分析可能會造成未受Dropout影響的基因表達發生變化，但ScImpute算法可以利用其他類似細胞中不太可能受Dropout影響的相同基因信息，在不引入新偏差情況下計算缺失值。有研究表明，MAGIC和scImpute算法都依賴于相似細胞基因數據，這會消除細胞間的隨機性，而Huang等［19］提出的SAVER算法來接收具有唯一分子索引的矩陣后，假定每個基因都遵循Poisson-Gamma模型，使用多元廣義泊松回歸模型的貝葉斯分析還原基因表達水平，消除技術差異的同時還可保留不同細胞間的生物學差異。DrImpute是一種集群分析算法，通過使用Spearman和Pearson相關系數計算距離矩陣，可將Dropout從真正的零值中有效地分離出來。Gong等［20］將DrImpute算法與多種插補算法進行性能比較發現，與MAGIC和scImpute算法相比，DrImpute可以恢復更多的缺失值，提高后續細胞類型識別和擬時間推斷的準確性。受到上述軟件的啟發，Talwar等［21］提出了AutoImpute自編碼分析算法，通過學習scRNA-Seq數據的固有分布和模式來尋找缺失值。通過與現有的9種獨立數據集的插補算法進行比較，AutoImpute被證實是唯一能對最大數據集進行插補而不會消耗內存的算法。

3 批次效應校正

完整的RNA測序流程包括細胞分離、RNA提取、文庫構建、上機測序及測序后數據分析等多個環節，但不同實驗室、不同時間以及不同人員操作會造成批次效應，影響結果可靠性。批次效應也成為scRNA-Seq技術中常見的變異來源［22］。由于scRNASeq與bulk RNA-Seq在數據特征上具有差異，常規用于 bulk RNA-Seq的批次校正算法，如RUVseq［23］和svaseq［24］等算法可能并不適用于scRNA-Seq。但是，在scRNA-seq研究中，批次之間的種群組成通常并不相同，即使假設每個批次中存在相同的細胞類型，數據集中每種細胞類型的豐度也會根據細胞培養或組織提取、解離等過程中的細微差異而變化，因此造成變異的因子并非僅考慮技術性因素。為了校正單細胞測序中的批次效應，多種scRNASeq數據校正算法被開發，包括ComBat［25］、相互最 近 鄰（Mutual nearest neighbours，MNN）［26］和Scanorama［27］等。當批次信息可用時，ComBat算法使用參數和非參數經驗貝葉斯框架通過批處理效應變量的加法組合來描述基因表達。Haghverdi等［26］提出的MNN算法是計算成對細胞的余弦歸一化表達譜之間的歐氏距離，再根據每個批次中共享種群的偏差來調整批次效果。盡管MNN和ComBat是常用的分析算法，但研究發現MNN算法性能要優于ComBat。Hie等［27］最近提出的Scanorama是采用一種可識別并準確整合數據集合的算法，利用匹配的信息進行批次效應校正，相比于MNN算法，該技術不需要依賴于數據集的順序，將鄰近搜索優化為低維嵌入的基因表達譜，極大減少了搜索時間。

4 降維分析

高維性是scRNA-Seq數據的顯著特點，數據分析時常常要用到降維分析法。主成分分析（Principal component analysis，PCA）作為一種經典的無監督降維算法，借助正交變換使線性維數減少，產生一組不相關的分量，通過最大化投影數據的方差，將高維數據投影到低維線性空間上。其具有兩大主要優勢：第一，PCA通過正交線性投影可以消除基因間的冗余，被用作多種降維方法的預處理步驟。第二，PCA可將高維數據投影到低維線性空間上，可以預測多維數據的相關性。研究表明［28］通過分散表達水平來過濾基因，然后選擇數百個最具可變性的基因來捕獲整個種群的重要特征。PCA已成功應用于scRNA-Seq數據分析中［29-31］，以捕獲細胞異質性的整體結構，其局限性在于無法可視化細胞聚類和細胞類型識別所必須的局部結構。

為了彌補PCA無法可視化的局限性，t分布隨機領域嵌入（t-distributed stochastic neighbor embedding，t-SNE）算法被引入單細胞測序分析。Alexander等［32］提出的t-SNE是一種用于高維數據可視化的非線性分析算法，通過捕獲局部結構，將原始高維空間中不相似單元以大距離建模，而相似單元則以小距離建模，在不丟失數據點間相對距離的基礎上，將高維數據嵌入到二維或三維空間中進行可視化。通過降維與最近鄰網絡相結合來考慮數據點之間的局部距離，目的是分離不同的群集。t-SNE可以通過構造概率分布來描述數據集，相似的單元格分配概率高，相異的單元格分配概率低，在高維空間中相似的單元將在低維空間中聚集在一起。t-SNE在維持相似細胞群集能力方面優于PCA。目前，t-SNE還不能很好地捕獲全局結構，如群集之間的距離。盡管t-SNE在scRNA-Seq數據可視化方面取得了成功，但仍存在兩種算法的缺陷［33］。首先，由于t-SNE的隨機性，同一數據集在不同的運行中可能產生不同的可視化效果。為了獲得對種群結構的認識，可能需要對同一數據集進行多次t-SNE運行。其次，雖然t-SNE將原始空間中相似單元格放置在低維空間中來維持簇，但原始空間中不相似單元格不一定會在低維空間中按比例放置。最近，一種基于黎曼幾何和代數拓撲理論的UMAP工具［34］被開發出來，其性能和效率均優于t-SNE。UMAP工具能夠沿著分化軌跡排列簇并保留瞬時細胞的分化連續體，通過在二維或三維圖上覆蓋標記基因的表達或與生物過程有關的一組基因的活性，捕獲scRNASeq數據中局部和全局結構。

5 細胞亞型鑒定

在特定條件下，對組織中的細胞亞群進行鑒定是scRNA-Seq數據分析的關鍵目標之一，其結果可以揭示細胞異質性［35］。結合降維分析方法，通過聚類分析實現細胞亞群的鑒定。在無監督聚類分析中，主要以分層聚類和K-means聚類為主。分層聚類無需預先定義聚類數量，以聚集或分裂的方式進行連續合并或拆分，目前常用的工具包括SINCERA［36］和bigSCale［37］。其中，Iacono等［37］提出的bigSCale工具框架構建了一個概率模型來定義所有可變性成對細胞之間的表型距離。與在簡單或混合概率模型中假設負二項式，伽馬或泊松分布的其他方法相比，bigSCale工具構建了一個高精度、全面的噪聲數值模型，通過將P值分配給每個基因來量化細胞間距離。而K-means聚類［11］則是先確定簇中心，再將細胞分配到最近的簇中心，迭代優化質心位置，將細胞分為 k個簇，根據質心聚類中細胞的平均值重新計算質心，工作速度快于分層聚類。以上兩種傳統聚類方法都會受到數據規模和噪聲的影響。為此，Lin等［38］在聚類前通過插補和降維進行聚類（CIDR）使用非線性最小二乘回歸擬合數據，并對零值進行插補來減弱Dropout影響。該算法可識別并評估Dropout與基因表達水平之間的關系，計算基因表達譜之間的差異。通過實驗證實CIDR運算速度遠快于傳統算法。近年來，新的聚類算法不斷被開發，如graph-based 聚類包括SNN［39］和RaceID2［40］，這些算法將單元格嵌入圖形中，每個邊代表兩個單元格之間的相似度，將圖形劃分為高度互連的模塊，具有高效性和穩定性。Seurat基于共享近鄰（Shared nearest neighbor，SNN）聚類算法來識別細胞簇，通過差異表達或方差分析來識別不同亞群標記物，基于表達水平的相似度的不同構建共享近鄰網絡。為了證明SNN算法的有效性，Xu等［39］通過在不同結構的數據集上進行測試發現，與原始數據的研究結論相同。為開發出一種可靠的方式推斷分化軌跡，Grün等［40］提出RaceID2工具，該軟件適合測試微分動力學，在增加集群數后可通過識別集群內的飽和點確定亞群數量，使數據比K-means聚類更可靠，且已在動物實驗中證實。除上述方法外，單細胞一致性聚類（Single cell consensus clustering，SC3）是特別為scRNA-Seq數據開發的聚類算法，通過共識方法將多個聚類算法組合在一起，具有高度的準確性和魯棒性［41］，相比于K-means，SNN和SINCERA算法分析，SC3缺點是運行時間長，但準確性最高［42］。

6 差異表達分析

差異表達基因分析可以檢測不同細胞類型、不同細胞亞群間的mRNA豐度，通過組間比對，獲得不同樣本或不同處理方法對基因表達水平的影響，或上調或下調［43］，進而可對差異表達基因進行功能分析，如通過基因本體分析（Gene ontology，GO），確定基因所參與的生物學過程、分子功能及細胞組分，通過KEGG分析差異基因參與的信號通路。盡管scRNA-Seq結果可以鑒定差異表達基因，但其也存在一定的局限性。首先，由于單細胞測序數據通常具較高的背景噪音，很多低表達或者中等表達水平的基因不能被有效檢測到。所以，針對bulk RNASeq數據開發的差異表達檢測算法，并不完全適用于scRNA-Seq。針對scRNA-Seq的差異表達算法被陸 續 開 發，如SCDE［44］，MAST［45］，Census［46］和BCseq［47］等。其中，SCDE是一種運用貝葉斯算法，從單個測量中獲得的不確定信息，使用泊松過程來解釋Dropouts，通過對比分析證明SCDE算法具有比傳統方法更高的靈敏度［44］。MAST算法采用線性模型對轉錄陽性表達的平均值進行建模，同時控制模型的離散性和技術因素。其采用廣義相加模型（Generalized additive models，GAMS）與Tobit模型進行正態分布。Finak等［45］通過比較發現SCDE算法檢測的差異表達基因數量高于MAST，但MAST算法的特異性更高。Qiu等［46］研究發現Census算法可將常規的相對表達量轉換為相對轉錄本計數，與標準的讀取計數相比，使用回歸技術更容易建模，而且顯著提高準確度。BCseq算法無需指定偏差的來源或格式即可校正表達量化中的偏差，即以自適應的方式糾正固有偏差，有效降低技術噪音。Chen等［47］通過對比發現BCseq算法在細胞類型分類性能上優于MAST和SCDE。盡管多種算法被用于差異表達基因分析，但不同算法處理后的數據結果仍存在一定偏差。早期一項針對36種基因差異表達分析算法有效性的研究發現，不同算法得到的差異表達基因在特征及數量上均存在顯著差異［48］。因此，測序后數據分析算法的優化仍然是今后的一項重要工作。

7 擬時序分析

擬時序分析法是指根據單個細胞的基因表達模式推斷出細胞發育或分化的動態路徑［49］。與bulk RNA-Seq不同，單細胞測序可以沿著一個連續發展的過程對細胞進行排序，在軌跡的開始、中間和結束狀態對細胞類型進行識別，進展越少越接近原始細胞狀態，進展越多越接近終點細胞狀態。針對scRNA-Seq擬時序分析開發的算法有Monocle［50］，Waterfall［51］，TSCAN［52］，Sincell［53］，SLICER［54］和Wishbone［55］。Monocle算法將無監督的數據與反向圖形嵌入結合在一起，通過分化進程對細胞進行排序，揭示關鍵調控因子表達中的變化及細胞分化的新型調控因子。Pere?íni 等［50］將Monocle算法可應用于骨骼肌分化過程中，明確了一系列形態學和轉錄組動力學。該算法將每個細胞的表達譜表示為高維空間的一個點，高度相似的單元格之間添加連接邊，構建最小生成樹圖，找到最長路徑的對應轉錄序列，即可找到分化過程中單個細胞的動態路徑。當對體內連續生物學過程進行scRNA-Seq分析時，由于缺少足夠的信息不適合使用Monocle等傳統方法，因此Jaehoon等［51］開發了一種更為通用的算法Waterfall，它可以對連續生物學過程的多位單細胞數據集進行無偏差統計分析，該算法使用Hidden markov模型以無偏差方式確定擬時序上每個基因的表達狀態，并量化為隨時間變化的分子級聯圖，最后將基因表達水平與擬時序相關聯。盡管Waterfall算法中曾考慮過細胞聚類的影響，但并未對細胞聚類對細胞有序性影響進行系統評估，在此基礎上，Ji等［52］提出TSCAN工具，基于構建最小生成樹之前對細胞進行聚類可以降低復雜性，解決Monocle由于高復雜性造成的穩定性差。Juliá 等［53］開發的Sincell實現了兩種算法，以區分穩定細胞和噪聲影響的層次結構，第一種依賴于基因重采樣程序；第二種系統隨機生成的復制細胞，這些復制細胞遵循原始細胞的隨機模式。Sincell可以提供細胞狀態層次結構，同時考慮scRNA-Seq序列中的隨機因素。上述方法無法推斷對于非線性基因表達的變化和與過程無關的基因的分析，Welch等［54］開發SLICER使用局部線性嵌入來重建細胞軌跡，該算法可以推斷非線性的軌跡，無需了解過程即可選擇基因，并自動確定分支的位置和數量，通過對小鼠非細胞和神經干細胞的驗證，證實了此算法的有效性。如果要為多細胞構建分支軌跡，上述方法在分辨率和準確性上較為不足，Manu等［55］提出的Wishbone算法解決了此類問題，關鍵技術在于通過重復采樣邊緣子集來確定軌跡，并用動物實驗證明了Wishbone的準確性。因此，選擇合適的方法應主要依賴于數據集特點。

8 結語

隨著多種單細胞測序數據分析方法的開發與應用，在一定程度上促進了單細胞轉錄組學的發展，改善了高背景噪音和高變異性對數據產生的影響，為細胞異質性研究奠定了分子基礎。但單細胞測序數據分析仍面臨著新的問題和挑戰。首先，隨著單細胞測序數據集激增，如何提高軟件運行速度和儲存效率是目前需要解決的一項重要問題。其次，由于不同實驗室在實驗方案和數據處理流程方面存在差異，結果的室間比對較為困難。因此，我們仍然有必要對現有的分析方法進行優化，不斷開發新的高效的數據分析方法，進一步提升單細胞測序結果的準確性和可靠性。