Current best practices in single-cell RNA-seq analysis: a tutorial

摘要

這篇文章主要闡述了一個典型的單細胞 RNA-seq 分析的詳細流程，包括了以下流程：

1. pre-processing -數據預處理

   • quality control -質量控制
   • normalization -數據標准化
   • data correction -數據修正（batch correction, noise correction, etc.）
   • feature selection -特征基因選擇
   • dimensionality reduction -降維

2. cell-level downstream analysis -細胞水平的下游分析

3. gene-level downstream analysis -基因水平的下游分析

1. Pre-processing and visualization

experimental workflow

single-cell dissociation, single-cell isolation, library construction, and sequencing 分別為以上四步走

• As a first step, a single‐cell suspension is generated in a process called single‐cell dissociation in which the tissue is digested.
• plate-based & droplet-based 但是都存在問題，In both cases, errors can occur that lead to multiple cells being captured together (doublets or multiplets), non‐viable cells being captured, or no cell being captured at all (empty droplets/wells)
• Each well or droplet contains the necessary chemicals to break down the cell membranes and perform library construction. Furthermore, many experimental protocols also label captured molecules with a unique molecular identifier (UMI).

What is UMI?

UMIs allow us to distinguish between amplified copies of the same mRNA molecule and reads from separeate mRNA molecules transcribed from the same gene

• Libraries are labelled with cellular barcodes and pooled together(multiplexed) for sequencing

  Illumina 單細胞測序工作流程： https://links.jianshu.com/go?to=http%3A%2F%2Fweb.illumina.com.cn%2Flanding%2Fproducts_view.asp%3Fnewsid%3D324

1. Raw data

原始數據根據是否在單細胞文庫構建方法中使用 UMIs一般分為count matrice 或read matrices。

Raw data processing pipelines (alignment and quantification):

• Cell Ranger (Zheng et al, 2017)
• indrops (Klein et al, 2015)
• SEQC (Azizi et al, 2018)
• zUMIs (Parekh et al, 2018)

2. Quality control

在對單細胞基因表達數據分析之前，我們必須確保所有的細胞 barcode 數據對應那些活着的細胞。即一個 barcode 是細胞的唯一身份證。

QC操作一般建立在三個 QC covariates 上：

• 每個 barcode 對應的基因數量
• 每個 barcode 對應的count數量（count depth）
• 每個 barcode 線粒體基因占 count 數的比例

通過一定的閾值，來篩選掉那些過低count（細胞已經死亡發生裂解了）或過高count（細胞的雙峰或多峰）

但是需要注意的是僅僅單一考慮以上一種指標可能會對細胞信號出現一些誤解，因此需要把數據和生物學知識結合起來

例如，相對高的線粒體基因占比可能是細胞正處於呼吸過程中；或者是低 count 或 低 gene 數的細胞對應着靜態細胞群落，並且細胞有着更大 count 數的可能具有更大的體積。

因此需要協同考慮多種指標來進行全局的閾值選擇，以此避免那些 viable 細胞被不小心篩選掉。

目前有一些 doublet 檢測工具提供更好的解決方案如:

DoubletDecon: preprint: DePasquale et al, 2018

Scrublet: Wolock et al, 2019

Doublet Finder: McGinnis et al, 2018

除了以上對細胞完整性的檢查之外，還必須對各個細胞在轉錄組水平上執行 QC 步驟。原始計數矩陣通常包含了 20,000 個基因，通過濾除大部分細胞均表達的基因，這些基因也正是那些對細胞異質性貢獻度不那么大的基因，可以大大減少原始矩陣的基因數目。

設置此閾值的准則是目前預期想要得到最小細胞群落的大小，並為dropout效果留出余地。例如，去掉那些少於 20 個細胞表達的基因可能會使得更難檢測出那些少於 20 個細胞的細胞群落。

進一步的質量控制能夠直接進行再計數矩陣中，Ambient gene expression 指代那些 count 數並不源於 barcode 細胞的，卻來自於在文庫構建之前由於細胞裂解導致的細胞懸液增加了污染的 mRNA 數量，從而使得 count增大。但是這種問題能夠利用 droplet-based scRNA-seq 通過那些大量的空液滴作為空白對照進行矯正。

（例如 SoupX preprint: Young & Behjati, 2018 能夠直接應用於 count data）

我們知道質量控制是為下游分析服務的，一個數據質量好壞是由下游分析結果的性能緊密相關的（例如聚類中細胞群落的標注）。因此，我們在數據分析時往往需要多次進行質量控制的再回顧。從一個較為寬松的 QC 閾值開始，再逐漸嚴格通常來講是比較有益的。

該方法對於包含異類細胞群體的數據集特別有用，在這些數據集中，細胞類型或狀態可能會被誤解為低質量的異常細胞。在低質量數據集中，可能需要嚴格的QC閾值。 QC閾值不宜用於改善統計測試的結果。相反，可以根據數據集可視化和聚類中QC協變量的分布來評估QC。注意 data peeking 的發生（Data peeking goes something like this: you run some subjects, then do an analysis. If it comes out significant, you stop. If not, you run some more subjects and try again）

3. Normalization

歸一化和標准化的區別

https://cloud.tencent.com/developer/article/1486102

• 常用的歸一化Normalize處理目的是將離散程度很大的數據集中化, 對數據轉換能夠讓同一基因在不同樣本具有可比性（RPKM/TPM） ，在Seurat中 LogNormalize 便是利用log1p(value/colSums[cell-idx] *scale_factor)
• 常用的標准化Scale是基於之前歸一化結果，再添加z-score計算，考慮到了不同樣本對表達量的影響，消除了表達的平均水平和偏離度的影響，即統一數據的均值和方差

歸一化使用場景

• 對表達量的范圍有要求
• 表達量穩定，無極端值
• 數據不太符合正態分布時，可以使用歸一化
• 機器學習算法（SVM,KNN,神經網絡等）要求input數據為歸一化

標准化使用場景

• 離散程度大，存在異常值和較多噪音，用標准化可以避免異常值和極端值的影響

• 在分類、聚類、PCA算法中，使用z-score值結果更好

• pheatmap(dat) # scale之前
  n = t(scale(t(dat)))
  n[n>2]=2 # up-limitation
  n[n<-2]=-2 #down-limitation
  pheatmap(n,show_colnames = F,show_rownames = F)
  

3.1 Normalize cellular count data (to make it comparable between cells)

Count 在計數矩陣中有兩個維度，一個是細胞數目，一個是基因數目。表達矩陣中，每個count值都需要細胞捕獲、反轉錄、測序這三個步驟成功完成而記錄下來。相同細胞的測序深度（檢測出的基因數目）可能由於每個步驟中操作差異而不同。因此，在比較兩個細胞時，其中的任何差異都有可能是由於實驗測序誤差產生的，而不是真實的生物學差異。

在這里首先討論如何實現細胞之間的可比較性。通過歸一化能夠消除這種因為采樣等操作造成的偏差，例如scaling 計數數據，以利用縮放后的數據反映細胞內部基因相對的表達豐度，利用細胞內部的相對值放在同一水平比較不同細胞的基因表達量能夠減少實驗測序誤差，突出生物學差異。

最常見的標准化方法是count depth scaling，也被稱作counts per million或者CPM normalization。這個方法常用於bulk轉錄組，它會根據每個細胞的總表達量計算一個 size factor ，然后對其中各個基因表達量進行normalize。在CPM方法中，有一個前提假設，即數據集中的所有細胞最初都含有相同數量的mRNA分子並且計數深度的差異僅由於取樣操作。

由於單細胞數據是具有不同大小及分子數量的異質細胞群，更復雜的標准化方法應該適用。例如，Weinreb et al (2018) 利用了一種CPM拓展方法，在計算他們的size factor時，排除了那些占細胞分子總數至少5%的基因。主要是預防少量極高表達量基因的存在，從而顯現其他基因的差異。Scran包有個pooling‐based size factor estimation方法，允許更高的細胞異質性存在；另外Scran包在批次矯正和差異分析環節也比其他歸一化方法表現更好(Buttner et al, 2019)

scRNA-seq測序技術分為全長測序及3'端富集測序；全長測序考慮到了基因的長度，而后者則沒有考慮。TPM歸一化常用於全長測序數據。

3.2 Normalize (scale) gene counts (to make it comparable between the same genes expressed in different cells)

基因的歸一化是將gene counts的數據通過放縮(scaling)創造出其均值為0，方差為1(z scores)的一組數據。這樣處理的好處是所有的基因權重都是一樣的，在下游分析中能夠與其他基因平等的進行比較。目前對基因是否進行縮放尚未達成一致。（Seurat 選擇進行縮放，而Slingshot則不選擇進行縮放）對基因進行縮放意味着在下游分析中所有的基因都被平均權重，而不進行縮放則將基因表達的量級作為基因重要性的一種信息保留下來。為了盡可能保留多的生物學信息，本文避免對基因進行放縮。

對數據進行歸一化之后，數據矩陣通常進行log(x+1)轉化。這樣的轉換有三種重要的影響:

• 第一，對數轉換后的基因表達量之間的距離反應出了對數倍數改變量，即

  \[logFOLD=logx-logy \]

  它能夠很好的反應出表達量發生的變化

• 第二，對數轉換能夠減緩單細胞數據中均值-偏差的關系

• 第三，對數轉換可以減少數據的偏斜度（skewness），從而近似於許多下游分析工具對數據呈正態分布的假設

盡管scRNA-seq數據實際上不是呈對數正態分布，這三個影響使得對數轉換成為稍許粗糙但非常有用的工具。

對數轉換對於下游分析中檢測差異表達基因及批次效應矯正(Finak et al, 2015; Ritchie et al, 2015)。但是也應該注意到對數轉換的一個弊端，即給數據引入了虛假的差異表達效應(preprint: Lun, 2018)當標准化大小因子分布在測試組之間差異很大時，此效果尤其明顯。

4. Data correction and integration

上述的歸一化處理（normalization）試圖消除采樣導致的計數偏差，但是歸一化后的數據可能仍然包含一些偏差。數據矯正的目的是進一步消除技術性或生物性的偏差，如批次效應、dropout、細胞周期效應。這些因素並不總是能得到矯正，因此需要根據預期的下游分析來考慮矯正哪個因素。本文建議分別考慮生物學和技術性變量的矯正，因為它們用於不同的目的。

4.1 Regressing out biological effects

雖然矯正技術性偏差對於發現潛在生物學信號可能至關重要，矯正生物學變量卻能夠揭示出特定的生物學信號。最常見的生物學數據矯正是在轉錄組中除去細胞周期的影響。可以通過Scanpy和Seurat平台中實現的細胞周期評分進行簡單的線性回歸或再具有更復雜混合模型的專用程序包中（如scLVM 或 fscLVM）。可以從文獻中獲得用於計算細胞周期評分的標記基因列表(Macosko et al, 2015)。這些方法能夠用於消退其他已知的生物學效應，例如線粒體基因表達，這被解釋為細胞壓力的反應。

還有一些方面需要在矯正生物影響偏差之前考慮到。

1. 矯正生物影響並不能總是對解釋單細胞測序數據有幫助。盡管除去了細胞周期的影響，提高了對發育周期推測的合理性，但是細胞周期信號對於生物學同樣也能提供有益的信息。例如，可以根據細胞周期評分以確定增殖細胞群（proliferating cell populations）。
2. 除此之外，生物學信號必須在相應的環境背景下才能夠理解。假定生物學過程發生在同一生物體內，則這些過程之間存在依賴性。因此，矯正了一個過程則可能誤將另一個信號掩蓋住。
3. 最后，細胞大小的變化也與細胞周期造成的轉錄組效應有關。因此，通過歸一化或專用工具（如cgCorrect Blasi et al 2017 ）也能部分矯正細胞周期對單細胞測序數據的影響。

4.2 Regressing out technical effects

用於以上生物學偏差的回歸模型的變體也能應用於技術性偏差。單細胞數據中最主要的技術性偏差是count depth和 batch。盡管歸一化對數據進行放縮，使得細胞之間的基因count數能夠在同一水平進行比較，count depth效應依舊存在於數據中。例如，細胞在大小上具有差異，因此它們的mRNA分子數也不盡相同。但是，技術效應在歸一化后仍可能存在，因為沒有一種放縮方法能夠推測那些由於不良取樣過程而未檢測到的基因表達量。對計數深度影響進行回歸能夠提高軌跡推測算法（找到細胞之間的轉換）。

一種基於回歸策略消除count effect的替代方法是使用更嚴格的歸一化過程，例如下采樣或非線性歸一化方法（請參見“歸一化”部分）。這些方法可能與基於板的scRNA-seq數據集特別相關，其中每個細胞計數深度的較大變化可以掩蓋細胞之間的異質性。

4.3 Batch effects and data integration

Batch effect

當將細胞分為不同的組別時，批次效應將會發生。這些組中的細胞可能來自不同芯片、不同測序通道或在不同時間節點收集得到的。細胞所經歷的不同環境可能會影響轉錄組的測序或轉錄組本身。所產生的影響存在多個層面：實驗中的細胞組之間，同一實驗室中進行的實驗之間或不同實驗室的數據集之間。本文將對前兩種情況作出區分。

校正在同一實驗中的樣品或細胞之間的批次效應是一種經典情況，稱為bulk RNA-seq 中的批次效應。我們將其與來自多個實驗的數據整合（data integration）區分開。通常使用線性方法校正批處理效果，而將非線性方法用於數據集成。

一個經典的批次效應矯正方法ComBat(Johnson et al, 2006)在中低復雜度的單細胞實驗中性能也良好。**ComBat包含了基因表達的線性模型，其中在數據的均值和方差中均考慮到了批次效應的影響。無論采用何種計算方法，最好的批效應校正方法都是通過巧妙地設計實驗來規避。可以通過將不同實驗條件和樣品中的細胞進行合並來避免批次效應。利用細胞標記(preprint: Gehring et al, 2018)，或者通過遺傳變異(Kang et al, 2018)，也能夠解構實驗中合並的細胞。

Data integration

與批次效應相比，數據整合方法的挑戰圍繞數據集之間的差異。估計批處理效果時，ComBat使用在同一批次中的所有單細胞來對應批次參數。這種方法將批次效應與細胞培訓或狀態之間的生物學差異混淆。數據整合方法，例如:

Canonical Correlation Analysis CCA	Butler et al, 2018
Mutual Nearest Neighbours MNN	Haghverdi et al, 2018
Scanorama	preprint: Hie et al, 2018
RISC	preprint: Liu et al, 2018
scGen	preprint: Lotfollahi et al, 2018
LIGER	preprint: Welch et al, 2018
BBKNN	preprint: Park et al, 2018
Harmony	preprint: Korsunsky et al, 2018

盡管數據整合方法也能夠被應用於簡單的批次效應校正問題，考慮到非線性數據整合方法的自由度增加，我們建議要注意到過度校正的出現。例如，在更簡單的批次校正設置中，MNN的表現

第一張圖片每個細胞都按照樣品來源進行着色，可以很明顯看出樣品與樣品之間有分層，說明了該數據集有很明顯的批次效應。而利用ComBat在小鼠腸上皮細胞數據集上進行批效應校正后，可以看見各類樣品在UMAP投影中分布更均勻，批效應明顯減弱。

Expression recovery

技術校正的進一步形式是表達恢復（expression recovery, also denoising or imputation）。單細胞轉錄組sexual包含了各種噪聲源(Gru¨n et al, 2014; Kharchenko et al, 2014; Hicks et al, 2017).噪聲的最主要方面主要是dropout造成的。現有的一些工具通過推測dropout事件，用合適的表達量去替代這些0值，進而降低數據集的噪音。

例如：MAGIC: van Dijk et al, 2018; DCA: Eraslan et al, 2019; scVI: Lopez et al, 2018; SAVER: Huang et al, 2018; scImpute: Li & Li, 2018 。

已被證實進行表達恢復能夠提高基因與基因相關性的估計(van Dijk et al, 2018; Eraslan et al, 2019)。不僅如此，這個步驟也能夠與歸一化、批次效應校正、其他在scVI工具中實現的下游分析進行整合。盡管大多數據校正方法將歸一化后的數據作為輸入，一些表達恢復方法卻基於預期的負二項式噪音分布，故輸入為原始計數數據。

當應用表達恢復時，應該注意到沒有一個方法是完美的。因此，任何方法都可能過校正或校正不足。實際上，已報道由於表達恢復造成了錯誤的相關信號。鑒於在實際應用中難以評估表達恢復的准確度，這種情況對用戶思考是否對數據進行去噪處理構成了挑戰。由於這些考量，目前尚無關於應如何對數據去噪的共識。謹慎而言，我們將僅在數據可視化時才使用表達式恢復，而在探索性數據分析過程中不使用假設的數據。因此全面的實驗驗證相當重要。

5. Feature selection, dimensionality reduction and visualization

一個人類的單細胞rna測序數據集能夠包含高達25,000個基因的表達量。對於一個給定的scRNA-seq數據集，許多基因並不能提供有用的信息，並且大多基因會包含0值。盡管在QC步驟中除去了這些zero count gene，對於一個細胞向量而言，也有高達15,000個維度（基因）。為了減少下游分析工具的計算負擔，減少數據中的噪聲並使數據可視化，可以使用多種方法來減少數據集的維數。

5.1 Feature selection

降維的第一步通常是基因挑選。在這個步驟中，數據集將進行篩選以保留那些可以顯示出數據集中差異的那一部分“informative”基因。因此經常使用highly variable genes(HVGs)。根據任務和數據集的復雜性，通常選擇1000-5000個HVG用於下游分析。數量更多的HVG會呈現出更好的PCA結果。

在Scanpy和Seurat中都實現了一種簡單而流行的選擇HVG的方法。在此，通過基因的均值表達對基因進行分類，並在每個分類中選擇方差均值最高的基因作為HVG。該算法存在不同的風格，對於Seurat則使用count data，而對於Cell Ranger 則使用對數轉換過的數據。不過最好還是在技術偏差校正之后選擇HVG，以避免選擇僅由於例如批次效應導致的高度可變基因。Other methods for HVG selection are reviewed in Yip et al (2018).

5.2 Dimensionality reduction

在基因選擇過后，單細胞表達矩陣的維度能夠通過專門的降維算法進一步降低。這些算法能夠使得表達矩陣嵌入在低維空間中，該空間旨在以盡可能少的維數捕獲數據中的深層結構。其原理是單細胞RNA-seq數據本身具有低維性(Heimberg et al, 2016)；或者是可以用比基因數目少得多的維數來充分描述存在在細胞表達譜圖上的生物學特征。降維就是找到這些特征維度。

降維的兩個目標是為了可視化（visualization）和匯總（summarization）。可視化是試圖將數據集以二維或三維投影的最佳方式展現出來。這些降維后的特征值作為散點圖上的坐標，以獲得數據的直觀表示。匯總沒有規定輸出成分的數量，反而對於描述數據中存在的可變性，較高的分量變得不再重要。匯總技術可用於通過查找數據的固有維數將數據縮減為其基本組成部分，從而有助於進行下游分析。盡管對輸出結果不進行二維可視化，但可以使用匯總方法利用top的簡化后的成分進行可視化。但是專用的可視化技術通常能夠更好地表示差異性。

降低的維度是通過特種空間維度（基因表達向量）的線性或非線性組合生成的。特別是在非線性情況下，在此過程可能會降低維度的可解釋性。圖4中顯示了一些常用的降維方法的示例應用。本文主要闡述如何在常見的降維方法之間進行選擇。A more detailed review of dimensionality reduction for single-cell analysis can be found in Moon et al (2018).

PCA

主成分分析時一種線性方法，通過最大化與其他每個維度中的殘差來生成縮減的維度，即主成分。盡管PCA並沒有在幾個維度上或非線性方法中獲得數據的結構，但它是許多聚類或軌跡推斷分析工具的基礎。事實上，PCA通常用作非線性降維方法的預處理步驟。PCA通過其前N個主成分來總結數據集，其中N的選定可以用"elbow"heuristics 或基於置換測試的jackstraw方法確定。PCA的簡單線性性質有以下優勢：

在降維空間中的距離與該空間所有區域中具有一致的解釋性。因此我們可以將感興趣的量與主成分相關聯，以評估其重要性。例如可以將主成分投影到technical nuisance covariates （如 percent_mt, feature count）以研究QC的性能，數據矯正和歸一化步驟。

Diffusion maps

擴散圖是非線性數據匯總技術。由於擴散成分強調數據中的轉移，通常，每個擴散成分（擴散圖的維度）突出顯示不同細胞群體的異質性。

5.3 Visualization

出於可視化的目的，標准實踐手段是使用非線性降維方法。最常見的則為t分布隨機鄰近嵌入（t-SNE）。t-SNE維度着重於全局結構為代價獲取局部相似性。因此，這些可視化會誇大細胞群體之間的差異，而忽略掉這些群體之間的潛在聯系。另一個困難是其perplexity參數的選擇，因為t-SNE圖可能會根據其值顯示出明顯不同的簇數。t-SNE常見的替代方法是，Uniform Approximation adn Projection method(UMAP), 或者是基於圖的工具（graph-based tool），例如SPRING(Weinreb et al, 2018)。UMAP和SPRING的force-directed 布局算法ForceAtlas2可以說是底層拓撲的最佳近似。

但是UMAP的優勢在於它的速度和對大數據集的兼容。因此，在沒有特定生物學問題的情況下，UMAP是最佳的探索性數據可視化工具。此外，UMAP還可以匯總二維以上的數據。

6. Stages of pre-processed data

盡管上述內容按順序概述了scRNA-seq中常見的預處理步驟，但下游分析通常更喜歡采用不同層次的預處理數據，並且建議根據下游應用來適應預處理。

預處理分為5個數據處理階段

1. raw data
2. normalized data
3. corrected data
4. feature-selected data
5. dimensionality-reduced data

這些數據處理的不同階段被分為三個層次：measured data 用於統計檢驗, corrected data 用於數據比較可視化, reduced data 用於下游分析

細胞和基因的質量控制應始終執行。但后續操作可能不需要操作，比如單批次的數據不需要數據校正。

6.1 Downstream analysis

經過預處理后，使用稱之為下有分析的方法去提取生物學信息並描述其潛在的生物學系統。這些描述通過將可解釋模型擬合到數據而獲得的。例如具有相似基因表達譜的細胞群，它們代表同一種細胞類型的群落。在相似細胞中基因表達的微小改變能夠表示細胞間的連續（分化）軌跡；或者基因具有相關聯的表達模式表明它們具有共調控關系。

下游分析可以分為細胞水平和基因水平兩種方法

細胞層次的分析一般着重兩個結構的描繪：聚類和軌跡分析。廣義地講，聚類分析方法試圖根據細胞分類為組別來解釋數據的異質性。而在軌跡分析中則恰恰相反，數據被認為是動態過程中的某一時刻的快照，可以幫助理解某個動態發育過程。基因層面是：差異分析，富集分析，互作網絡

7. Cluster analysis

7.1 clustering

聚類是任何單細胞分析中第一個中間結果。聚類能夠使得我們推測細胞的身份。通過基於細胞基因表達譜圖的相似性給細胞分組能夠得到細胞的類。表達譜圖相似性通過距離矩陣來決定，距離度量通常將降維的表示形式作為輸入。相似性評分一般是在PC降維后的表達空間上計算歐氏距離。存在兩種根據這些相似性分數生成細胞簇的方法：clustering algorithms & community detection methods

Clustering

聚類是一個經典的非監督機器學習問題，直接借助距離矩陣進行分析。細胞通過最小化類別間距離或在減少的表達空間中找到密集區域，被歸到相應的簇中。最常用的k均值聚類算法（k-means clustering）通過確定聚類中心並將細胞分配給最近的聚類中心，將細胞分為k個類別。質心位置經過迭代優化。這個方法需要輸入預期的簇數，這通常是未知的往往需要啟發式校准。k均值算法在單細胞數據中的應用在所使用的距離矩陣中有所不同。標准的歐氏距離的替代包括cosine similarity(余弦相似性)， correlation-based distance metrics(Kim et al, 2018), 或者SIMLR方法，這利用了Gaussian kernels(Wang et al, 2017)為每個數據集學習距離度量。最近的比較表明，correlation-based distance 表現更佳。

Community detection methods

社區檢測方法是圖形划分算法，因此依賴於單細胞數據的圖形表示。使用K最鄰近方法（KNN圖）獲得此圖表示。細胞在圖中表示為節點（node），每個細胞都與其K個最相似的細胞相連，這些細胞通常使用PC減少的表達空間上的歐式距離獲得。根據數據集的大小，K通常設置為5到100個最近的鄰居之間。結果圖捕獲了表達數據的拓撲結構。表達空間的密集采樣區域表示為圖的密集連接區域。使用社區檢測方法檢測這些密集區域，社區檢測通常比聚類方法要快，因為只有相鄰細胞對才被視為同種細胞類型。因此，這種方法大大減少了可能的群集的搜索空間，這是在Louvain算法（Blondel）中實現的等人，2008年）。

7.2 Cluster annotation

在基因水平上，通過找到每個群落的標記基因聚類數據能得以被分析。這些稱之為marker genes 表征了每個群落，並被用來用有意義的生物學標記對其進行標記。該標簽表示了群落中細胞的類型（身份）。由於任何聚類算法都會產生數據分區，因此只能通過其代表的生物學特性的成功注釋來確定已識別聚類的有效性。

雖然很可能會嘗試假定單細胞數據中的簇代表真實的細胞類型，但是有多個變量決定了細胞的類別。

1. 首先，並不總是清楚什么構成了一個細胞類別。例如，盡管“T細胞”可能是某些人滿意的細胞類別標記，但其他人可能會在數據集中尋找T細胞亞型並區分CD4+和CD8+ T細胞
2. 此外，相同類型的細胞處於不同轉台可能會在不同的群落中檢測出

由於上述原因，細胞類型的定義存在爭議。最好使用術語"cell identities"而不是“cell types”。因為就像上述提到的T細胞，T細胞本身是T細胞，但是在發育階段和狀態不一樣，身份會發生改變 ，但不能簡答理解為它的類型發生改變了。

識別和注釋cluster依賴於外部信息，即那些描繪了單個細胞identities的預期表達模式。即reference database。這些數據庫大大的促進了細胞identities的標注。在沒有相關參考數據庫的情況下，可以通過將數據衍生的標記基因與文獻中的標記基因比較來注釋細胞身份，或者直接將文獻來源的標記基因的表達值可視化。應當注意的是，后一種方法將用戶限制於對從bulk expression研究中得出的細胞類型的經典理解，而不是cell identities。此外，已經表明，常用的細胞標志物在限定細胞identities方面能力有限。

利用reference database 標注細胞群落的兩種策略：

• 利用data-derived 標記基因
• 利用全基因表達譜

標記基因集能通過應用差異表達檢測找出（differential expression DE testing）。一般而言，我們更關注於那些群落中上調基因。由於標記基因被認為是強有力的差異表達因素，簡單統計檢驗（如Wilcoxon 秩和檢驗或t檢驗等）通過兩個群落中基因表達差異的程度去排序基因。各自檢測中排序靠前的基因被視作marker gene。可以通過比較來自數據集的標記基因和來自參考數據集的標記基因（通過富集測試，Jaccard索引或其他重疊統計數據）來注釋聚類。

檢測標記基因需要注意兩個方面

1. 由於DE測試中，0假設使基因在兩組之間具有相同的表達值分布，但是這兩組是通過聚類得到的，本身即是存在差異的。即是是隨機數據的聚類時，也能發現標記基因。因此制定合適的P值很關鍵，
2. 並且單一變量的簇注釋方法在特定情況下不建議使用。

7.3 Compositional analysis

在細胞層次上，我們可以根據數據組成結構來分析聚類數據。成分數據分析主要圍繞歸入每個cell-identity 類別的細胞比例。這些比例可能會由於疾病而發生改變。例如，沙門氏菌感染已顯示會增加小鼠腸上皮中腸上皮細胞的比例。

在單細胞數據集分析成分變化需要足夠多的細胞數量以穩健地評估細胞類別的比例，並且足夠的樣本數量以評估細胞身份簇組成中預期的背景變化。

8. Trajectory analysis

8.1 Trajectory inference

細胞多樣性不能充分被離散分類系統（如聚類）描述。驅動觀測到的異質性發展的生物學過程是連續過程。因此，為了獲得細胞身份之間的轉換，分支分化過程或生物學功能的緩慢非同步變化，我們需要動態的基因表達模型，這類方法稱之為軌跡推斷（TI）

軌跡推斷方法將單細胞數據解釋為連續過程的快照（snapshot），通過尋找穿過細胞空間的路徑來重建該過程，該路徑將相鄰細胞之間的轉錄變化降低至最低。細胞沿着這個路徑的順序被稱之為pseudotime。盡管該變量與距離root細胞的轉錄距離有關，但通常被解釋為發育時間的代名詞。

從Monocle和Wanderlust建立TI領域以來，可用的方法數量激增。當前可用的TI方法在建模路徑的復雜性方面有所不同。模型的范圍從簡單的線性或分叉軌跡到復雜的圖形，樹或多分叉軌跡。沒有一種單獨的方法對所有類型的軌跡都能表現出最佳性能。相反，應根據預期軌跡的復雜性選擇TI方法。對比顯示，Slingshot（Street等人，2018）優於其他方法（從線性模型到雙叉模型和多叉模型）的簡單軌跡。如果期望更復雜的軌跡，作者建議使用PAGA（Wolf等人，2019）。如果確切的軌跡模型是已知的，則可以選擇使用更專業的方法來提高性能（Saelens等，2018）。通常，應使用替代方法來確定任何推斷的軌跡，以避免方法偏差。

在典型的工作流程中，TI方法應該用於降維或校正后的數據。由於在細胞內同時發生多種生物學過程，因此消退其他過程的生物學效應以區別我們預期的軌跡是有效果的。

例如T細胞在成熟過程中可能正在經歷細胞周期轉換，此外由於幾種性能高的TI方法依賴於聚類數據，因此TI通常在聚類之后進行。推斷軌跡中的類可能表示穩定或亞穩定。隨后，可以將RNA速度疊加到軌跡上以增加方向性。

需要注意的是推斷的軌跡並不代表真實的生物過程，這些僅表示轉錄相似性。因此需要進一步信息對其驗證。

8.2 Gene expression dynamic

一個推斷的軌跡不是擬合轉錄噪音的結果而是在基因水平上分析該軌跡。在偽時間內平滑變化的基因表征了軌跡，可用於識別潛在的生物學過程。此外，這組與軌跡相關的基因被認為是包含調節建模過程的基因。Regulator 基因幫助我們了解生物過程的觸發方式和原因，並代表潛在的葯物靶標。

8.3 Metastable states

亞穩態是一個相對穩定但又會變化的狀態

亞穩態狀態軌跡的細胞水平分析研究了偽時間的細胞密度。假定細胞采樣的過程是無偏的，沿軌跡的密集區域表示了優選的轉錄組狀態。當將軌跡解釋為時間過程時，這些密集區域可能代表例如發展中的亞穩態。我們可以通過繪制偽時間坐標的直方圖來找到這些亞穩態。下圖中B中，有很多中狀態，但C中直方圖認為這幾個密集的區域才屬於亞穩態。

整合分群與軌跡分析

clustering 和 trajectory 可以整合在一起展示結果。將分群的結果當成節點（node），將軌跡當成節點之間的橋梁（edge），所以動態與靜態的數據結合在了一起，利用partition-based graph abstraction PAGA 工具可以得到圖D效果。

結語

文章闡述了單細胞數據從獲得到預處理到分析等等一整套流程，雖然精讀了一遍，也學了很多背景知識，還是發現許多遺漏的點，等以后實戰中遇上了再查漏補缺建立自己的體系認知。