原文鏈接：

https://www.embopress.org/doi/10.15252/msb.20188746

主編評語

這篇文章最好的地方不只在于推薦了工具，提供了一套分析流程，更在于詳細介紹了有哪些方法可用、什么時候用哪些方法和有什么注意事項。這也是我們易生信培訓過程中廣泛討論的問題。

Abstract

單細胞RNA-seq使研究者能夠以前所未有的分辨率研究基因表達圖譜。這一潛力吸引著更多科研工作者應用單細胞分析技術解決研究問題。隨著可用的分析工具越來越多，如何組合成一個最新最好的數據分析流程也越來越難。我們詳細闡述了一個典型的單細胞轉錄組分析各個步驟的細節和注意事項，包括預處理（質控、標準化、數據校正、特征選擇、降維）和細胞/基因水平的下游分析等。基于獨立的比較研究，我們為每一步都推薦了當前最好的方法和操作建議。隨后把這些最好的工具整合成一個分析流程并應用于一套公共數據集的分析以演示其效果。案例研究具體可見https://www.github.com/theislab/single-cell-tutorial。這篇綜述為這個領域的新人提供了一份學習單細胞分析的指南，并且也能幫助老用戶更新他們的分析流程。

背景

近年來，單細胞RNA測序（scRNA-seq）大大提高了我們對生物系統的了解。我們已經能夠在研究斑馬魚、青蛙和渦蟲（zebrafish, frogs and planaria）等生物細胞異質性的同時發現先前未知的細胞群體。這項技術的巨大潛力激勵了計算生物學家開發了一系列分析工具。盡管開發者為了確保單個工具的可用性付出了巨大的努力，但是由于該領域的相對不成熟，對于單細胞數據分析的新手來說，入門的障礙是缺少一份標準指南。在本文中，我們提供了一份scRNA-seq分析的參考教程，并概述了當前的最佳實踐方案，為將來的分析標準化奠定了基礎。

分析標準化面臨的挑戰來源于越來越多的可用分析方法（截至2019年3月7日有385種工具）和數據集規模爆炸性的提高。使得我們一直在尋找新的方法來分析處理我們的數據。例如，最近已經有方法可以預測細胞分化過程中的命運選擇。盡管分析工具的不斷改進有助于產生新的科學推論，但它也使分析流程的標準化變得更為復雜。

標準化的另一個挑戰在于軟件技術方面。用于scRNA-seq數據的分析工具是用不同的編程語言編寫的 - 最主要的是R和Python。盡管跨編程語言的支持越來越多，但使用的編程語言確實影響了對分析工具的選擇。諸如Seurat，Scater或Scanpy等常用工具提供了集成環境來開發流程并包含大量分析工具。然而，出于維護的需要， 這些平臺限制了它們只能使用各自的編程語言開發的工具。引申開來，當前可用的scRNA-seq分析教程也存在語言限制，而且許多是圍繞這3大平臺進行講述的 （易生信2020年最新的單細胞課程會同時提供R和Python版本的最新分析流程）：

• R and bioconductor tools:
  https://github.com/drisso/bioc2016singlecell；
  https://hemberg-lab.github.io/scRNA.seq.course/; https://www.ncbi.nlm.nih.gov/pubmed/27909575?dopt=Abstract；

  如何使用Bioconductor進行單細胞分析？

• Seurat: https://satijalab.org/seurat/get_started.html;

• Scanpy: https://scanpy.readthedocs.io/en/stable/tutorials.html；

考慮到上述挑戰，我們選擇概述當前分析過程中每一步的最佳實踐和獨立于編程語言的常用工具，而不是評估一個完整的分析流程。我們將帶領讀者完成scRNA-seq分析流程的各個步驟，介紹每個步驟的最優方法，并討論分析過程中會遇到的陷阱和當前未解決的問題。當然由于一些工具比較新，并且缺乏系統比較，最好的工具很難決定，我們列出了流行的可用工具。列出的工具從基因count矩陣開始到潛在的分析終點。在我們的Github上有結合了已建立的當前最佳實踐流程的案例研究：https：//github.com/theislab/single-cell-tutorial/。我們在實際的示例工作流程中應用了當前的最佳實踐來分析公共數據集。該分析流程使用Jupyter notebook和rpy2整合了R和Python的工具。結合現有的文檔，這已經可以被視作一個可接受的分析流程模板了。

圖1. 典型單細胞轉錄組分析流程。

原始測序數據預處理和比對后獲得Gene count矩陣，作為本分析流程的開始。這些結果圖是Count矩陣采用最佳實戰分析流程進行預處理和下游分析獲得的。數據來源于Haber et al 2017腸上皮細胞文章。

預處理和可視化

原始測序數據經過處理得到分子計數矩陣（count matrix），或者reads count?（讀數矩陣）。這取決于單細胞文庫構建方案中是否包含唯一分子標識符（UMI，?unique molecular identifiers）（參考下面的Box1. Hemberg-lab單細胞轉錄組數據分析（六）- 構建表達矩陣，UMI介紹）。原始數據處理流程如Cell Ranger?(10X單細胞測序分析軟件:Cell ranger、從拆庫到定量),?indrops,?SEQC,?zUMIs負責測序數據質控，確定reads來源的細胞 (barcodes) （這一步也稱為demultiplexing）和mRNA分子 （生信寶典注*：單細胞測序不只獲得**mRNA，更準確說是帶poly-A尾巴的RNA，包括mRNA和lncRNA等*）、基因組比對和定量。獲得的reads或count矩陣的行數等于barcodes的數目，列數等于基因數目（生信寶典注*：也可能反過來，行為基因，列為**barcodes*）。這里使用術語barcodes而不是cell，因為分配給相同barcode的所有reads可能并不只是來源于同一細胞。barcode可能會錯誤地標記多個細胞（doublet），也可能不會標記任何細胞（空液滴/孔） (Hemberg-lab單細胞轉錄組數據分析（四）- 文庫拆分和細胞鑒定)。

盡管reads data和?count data?的測量噪聲級別不同 （生信寶典注*：基于UMI的數據，獲得的是分子計數，**count data，噪聲更低*），但在分析流程中的處理步驟是相同的。為簡單起見，在本教程中，我們將數據稱為計數矩陣（count data）。在reads和計數矩陣的結果不同的地方，將會特別提到reads矩陣。

Box1：scRNA-seq實驗流程的關鍵要素

從生物樣品到單細胞數據需要多步實驗操作。典型的工作流程包含單細胞解離、單細胞分離、文庫構建和測序。我們在這里簡要介紹這些步驟。不同protocols的更詳細的解釋和比較可參考下面三篇文章：

Ziegenhain?et al?(2017)：https://www.embopress.org/doi/10.15252/msb.20188746#msb188746-bib-0159

Macosko?et al?(2015)：https://www.embopress.org/doi/10.15252/msb.20188746#msb188746-bib-0083

Svensson?et al?(2017)：https://www.embopress.org/doi/10.15252/msb.20188746#msb188746-bib-0126

單細胞實驗的輸入材料通常是生物組織樣品。第一步，單細胞解離：消化組織產生單細胞懸液。為了分別分析每個細胞中的mRNA，必須分離單細胞。根據實驗方案不同，單細胞分離的方式也有所不同。

• 基于平板的技術將細胞分到到板上的孔中。

• 基于液滴的方法則依賴于微流體液滴捕獲單個細胞。

在這兩種情況下，都可能出現一些問題，如多個細胞一起被捕獲（doublets or multiplets）、非活細胞被捕獲或根本沒有細胞被捕獲（空液滴/孔）?；谝旱蔚姆椒ㄐ枰ㄟ^低的輸入細胞濃度來保持低的doublets率，因此空液滴是特別常見的 （生信寶典注*：一般beads和細胞的輸入比例是20:1*）。

每個孔或液滴均包含必要的試劑以裂解細胞膜并進行文庫構建 ?(生信寶典注*：植物單細胞就要注意了，需要提前去除細胞壁*)。文庫構建包括捕獲細胞內mRNA、反轉錄為cDNA分子并進行擴增等過程。因為文庫構建時每個細胞是獨立的，所以每個細胞的mRNA也就特異的標記了孔特異性或液滴特異性細胞barcode。此外，許多實驗方案還使用唯一分子標識符（UMI）標記捕獲的RNA分子。一般在測序之前需要先擴增細胞cDNA以增加其被檢測的可能性。但微量擴增更容易引入PCR偏好性。UMI使我們能夠區分測到的reads是來源于mRNA分子的不同擴增拷貝還是來源于獨立的mRNA分子，從而可以進行更準確的定量。

每個細胞單獨構建的cDNA文庫都帶有cell barcode和/或UMI（取決于protocol），后續將這些文庫混合在一起測序。測序產生的reads數據進行質量控制 (Hemberg-lab單細胞轉錄組數據分析（三）- 原始數據質控)，根據其barcodes序列分組（demultiplexing），并且進行后續比對定量。對基于UMI的protocols，reads的數據可以進一步demultiplexed以得到捕獲的mRNA分子的計數（count data）。也就是本套流程的起始輸入數據。

質控

在分析單細胞基因表達數據之前，我們必須確保所有barcode都對應于有效細胞 (viable cells，有活力的細胞)。質控有3個指標：

• 測到的轉錄本分子總數

• 測到的基因總數

• 來源于線粒體基因的轉錄本所占比例

質控就是檢查這3個指標的分布中是否存在異常峰并設置閾值去除。這些異常的barcodes可能對應于死細胞、細胞膜破損的細胞或doublets。比如，如果某個barcode對應的樣品測到的分子總數低、檢測到的基因數少、線粒體基因所占比例高，則表明該樣品可能存在細胞膜破損導致細胞質RNA漏出，只有線粒體中的RNA保留了下來。相反，如果某個barcode對應的樣品有異常高的總分子數和檢測到的基因數，則有可能這個樣品包含2個或以上細胞（doublets）。因此，高的總分子數通常用來過濾潛在的doublets。3個最近發表的doublets檢測工具提供了更好的解決方案：DoubletDecon,?Scrublet,?Doublet Finder。

分開考慮這三個QC變量中的任何一個都可能導致對細胞狀態的錯誤解讀。例如，線粒體計數相對較高的細胞可能是呼吸過程比較活躍?（生信寶典注*：如心臟細胞總**mRNA的30%是線粒體基因，具體見對一篇單細胞RNA綜述的評述：細胞和基因質控參數的選擇*）。同樣，其他QC變量也具有相應的生物學意義。總分子數和/或基因數低的細胞可能是處于靜息狀態的細胞群體；總分子數和/或基因數高的細胞也可能是細胞自身體積較大。實際上，細胞之間的總分子數自身也可能有很大差異（具體見Github上的案例研究）。因此，在做出是否過濾的單閾值決策時，應聯合考慮3個QC變量，并且應將這些閾值設置的盡可能寬松，以避免無意間濾除有效的細胞群。將來，考慮多元QC依賴的過濾模型可能會提供更敏感的QC選擇方式。

包含不同類型細胞混合體的數據集的每個QC指標可能會呈現多個分布聚集峰。例如，圖2D顯示了具有不同QC分布的兩個細胞群。如果未執行先前的過濾步驟（請注意，Cell Ranger也會執行細胞QC），則應僅將最低總分子數和基因數的barcode視為無效細胞。另一個閾值設定準則是考慮所選閾值過濾掉的細胞比例。比如過濾高分子總數細胞時，過濾掉的細胞比例不應超過預期的doublet率。

除了檢查細胞的完整性外，還必須在轉錄本水平上執行QC步驟。原始計數矩陣通常包含超過20,000個基因?？梢酝ㄟ^濾除在多數細胞中不表達的基因 (通常認為這些基因對細胞異質性分析貢獻不大)，大大減少該數目。設置此閾值的準則是使用感興趣的最小細胞亞群大小，并為dropout事件留出一些余地 (生信寶典注*：比最小細胞亞群大小數字再小一點*)。例如，濾除掉在少于20個細胞中表達的基因可能會使鑒定少于20個細胞的細胞亞群變得困難。對于dropout比率高的數據集，此閾值也可能會使不太大的細胞亞群難以被檢測到。閾值的選擇應隨待分析數據集中的細胞數量和既定的下游分析而定。

另外可以直接對計數矩陣 (count matrix)執行進一步的質量控制。輸入型基因表達（Ambient gene expression）是指某個barcode檢測到的轉錄本是源自其他細胞在建庫前發生破裂而釋放到細胞懸液中的mRNA。這些增加的外來計數會影響下游分析，如Marker基因鑒定或其他差異表達分析過程。由于基于液滴的scRNA-seq數據集中存在大量空液滴，因此可以通過空液滴建模分析細胞懸液中的RNA構成和豐度來校正這一影響。最近開發的SoupX使用這種方法直接校正count matrix。另外，在下游分析中直接忽略這些有強影響的輸入型基因也是處理這個問題的一個實用方法。

質控是用于保證下游分析時數據質量足夠好。由于不能先驗地確定什么是足夠好的數據質量，所以只能基于下游分析結果（例如，細胞簇注釋）來對其進行判斷。因此，在分析數據時可能需要反復調整參數進行質量控制 ?(生信寶典注*：反復分析，多次分析，是做生信的基本要求。這也是為啥需要上易生信培訓班而不是單純依賴公司的分析。*)。從寬松的QC閾值開始并研究這些閾值的影響，然后再執行更嚴格的QC，通常是有益的。這種方法對于細胞種類混雜的數據集特別有用，在這些數據集中，特定細胞類型或特定細胞狀態可能會被誤解為低質量的異常細胞。在低質量數據集中，可能需要嚴格的質量控制閾值。數據集的質量可以通過實驗質量控制指標來確定 (see Appendix Supplementary Text S2，簡單來說就是原始數據測序質量(FastQC的結果、序列比對率和比對到已經注釋的基因區的reads比率、實際檢測到的細胞數和預估細胞數是否吻合)。在QC閾值迭代優化過程中，要避免數據挑選 (data peeking)。QC閾值不應用于改善統計檢驗的結果。相反，可以根據數據集可視化和聚類中QC變量的分布來評估QC選取的閾值是否合理。

圖2. 質控變量的分布圖 (小鼠腸道上皮細胞數據集)
(A) 每個細胞檢測到的總分子數的分布用直方圖展示。右上角的子圖展示了總分子數小于4,000的部分的分布。因為在count數為1200處有個峰，所以閾值設置為了1,500。(B) 每個細胞檢測到的基因數的直方圖展示。在400個基因處有個小峰（噪音峰），閾值設置為700。（生信寶典注*：這閾值設置的隨意性也沒誰了~~~*）? 每個細胞檢測到的總分子數從高到底繪制rank plot，類似于Cell Ranger檢測并過濾空液滴的log-log plot。在總分子數為1500時，存在一個快速降低的拐點，則1500為篩選閾值。(D) 同時展示檢測到的基因數（縱軸）、總分子數（橫軸）和線粒體基因的比例 (顏色)。線粒體基因只在低基因數和低總分子數的細胞中比例高。這些細胞已經被上面設置的總分子數和基因數閾值過濾掉了。質控參數聯合可視化顯示更低的基因數篩選閾值可能也是足夠的。

陷阱和建議：

• 通過可視化檢測到的基因數量、總分子數和線粒體基因的表達比例的分布中的異常峰來執行QC。聯合考慮這3個變量，而不是單獨考慮一個變量。

• 盡可能設置寬松的QC閾值；如果下游聚類無法解釋時再重新設定嚴格的QC閾值。

• 如果樣品之間的QC變量分布不同（存在多個強峰），則需要考慮樣品質量差異，應按照Plasschaert et al. (2018)的方法為每個樣品分別確定QC閾值。

標準化

計數矩陣中的每個數值代表細胞中一個mRNA分子被成功捕獲、逆轉錄和測序（Box 1）。由于每個操作步驟固有的可變性，即便同一個細胞測序兩次獲得的計數深度也可能會有所不同。因此，當基于原始計數數據比較細胞之間的基因表達時，得到的差異可能來自于技術原因。Normalization可以通過調整計數數據 (scaling count data)等解決這一問題，以獲得細胞之間可比的相對基因表達豐度。

普通轉錄組分析有很多可用的標準化方法。盡管其中一些方法已應用于scRNA-seq分析，但單細胞數據特有的變異特征（例如technical dropouts:?zero counts due to sampling?）催生了scRNA-seq特異性標準化方法的發展。(什么？你做的差異基因方法不合適？)

最常用的標準化方法是測序深度標準化，也稱為“每百萬計數”或CPM normalization。該方法來自普通轉錄組表達分析，使用每個細胞的測序深度作為size factor對計數數據進行標準化。CPM標準化假設數據集中的所有細胞最初都包含相等數量的mRNA分子，并且計數深度差異來源于技術問題。這一方法的變體有：把size factor放大10^6倍如CPM，或size factor乘以數據集中所有細胞的測序深度的中位數。downsampling方法也基于同樣的假設，從原始數據中隨機抽取預設的等量reads以保證所有細胞測序深度相同。downsampling方法在扔掉一些數據的同時增加了technical dropout rates，而CPM或其它全局的標準化方式則不會影響。因此downsampling方法可能更真實地反應了相同測序深度下不同細胞的表達譜特征。

由于單細胞數據集通常由大小和分子數不同的異質細胞群體組成，因此通常需要更復雜的標準化方法。例如，Weinreb et al對CPM算法進行了擴展，在計算size factors時排除在任何細胞中總計數大于5％的基因。這一方法屏蔽掉少數高表達基因對總體表達變化的影響。軟件包Scran的pooling‐based size factor方法對細胞異質性的影響處理更好。首先把細胞合并到一起避免technical dropout效應，然后基于基因表達的線性回歸模型估算size factor。這一方法允許細胞有少于50%的差異表達基因，并且在不同的測試評估研究中這一標準化方法都表現最好。評估發現，scran比其他標準化方法對后續批次校正和差異表達分析效果更好。并且在scran作者小范圍的比較中也展現出這個方法穩定性比較好。

CPM,?high‐count filtering CPM和scran?使用線性全局縮放對計數數據進行標準化。還有非線性標準化方式可以處理更多的混雜因素的影響。許多此類方法涉及對計數數據進行參數化建模。例如，?Mayer et al?使用技術相關協變量（例如測序深度和每個基因的計數）作為模型參數對count data擬合負二項模型。模型的殘差作為標準化后的基因表達值。這種方法可以在標準化的同時校正技術差異和生物來源的差異（例如批次校正或對細胞周期影響的校正）。已表明，非線性標準化方法優于全局標準化方法，尤其是在有較強批次效應影響的情況下。因此，非線性標準化方法特別適用于基于板的scRNA-seq數據，因為不同板之間往往會產生批次效應。此外，基于板的數據與基于液滴的數據相比，每個細胞的測序深度變化更大。理論上非線性標準化方法或諸如downsampling的放法更適合于板的數據 (plate)，但仍需要進行比較研究以確認這一觀點。在本教程中，我們傾向于將標準化和數據校正（批處理校正、噪聲校正等）步驟分開，以強調數據的不同處理階段。因此，我們專注于全局縮放標準化方法。

我們不能期望一個標準化方法適用于所有類型的scRNA-seq數據。例如，Vieth et al?表明，reads data和count data適用不同的模型。確實，Cole et al?發現不同的標準化方法只對適合的數據集表現最佳，并認為應使用他們的scone工具對數據集進行評估以選擇適用的標準化方法。此外，scRNA-seq技術可分為全長和3'富集方法。全長方案的數據可能會受益于考慮了基因長度的標準化方法，而3'富集測序數據則不然。TPM歸一化是全長scRNA-seq數據常用的歸一化方法，它來自bulk RNA-seq分析。

標準化是對細胞計數數據進行縮放處理以使其在細胞之間可比，也可以在基因層面對基因計數進行歸一化 (scale)以便于基因內部進行直接比較。基因歸一化是指一個基因減去其在所有樣品表達的均值然后除以其在所有樣品表達值的標準差。歸一化后，這個基因在所有樣品表達值均值為0，用單位方差形式表示其表達值。歸一化后，所有基因在下游分析時權重是一樣的。（生信寶典注*：歸一化后，原始基因的表達豐度信息就沒了，換成了無標度的標準差的表示。）是否對基因進行歸一化目前尚無達成共識。盡管流行Seurat教程通常應用gene scaling，但Slingshot方法的作者在其教程中選擇了不對基因進行scaling。兩種選擇的爭議點在：所有基因不論表達高低在進行下游分析時權重一致 ?(生信寶典注*：**scale后所有基因權重一致)，還是基因表達量的絕對值對下游分析也有貢獻 (生信寶典注*：突出高表達基因對下游的更大影響*) 。為了從數據中保留盡可能多的生物相關信息，在本教程中，我們選擇不進行歸一化 （scaling genes）這一步操作。

標準化后，數據矩陣通常進行log(x+1)轉換。此轉換具有三個重要作用。

• 首先，對數轉換后的表達式值之間的差值可對應于對數轉換后的倍數變化，這是衡量基因表達變化的常用方法。

• 其次，對數轉換可減輕（但不能消除）單細胞數據的均值-方差關系 (mean-variance relationship) （生信寶典注*：均值-方差關系展示數據在特征空間的分布關系。方差越大數據分布越廣，后續采用線性回歸算法時效果越差。*）。

• 最后，對數轉換可以減少數據的偏態分布，從而使數據近似于正態分布，更符合許多下游分析工具對數據分布的假設要求。

盡管scRNA-seq數據實際上不是對數正態分布的，但這三個效果使對數轉換成為一種粗略但有用的工具。這一有用性在下游差異表達分析或批次校正時有更好的體現。但是應該注意的是，數據的對數轉換會在數據中引入虛假的差異表達結果。而且如果size factor在組間差異更大時影響尤其明顯。

陷阱和建議：

• 我們建議使用scran來標準化非全長數據集。另一種選擇是通過scone評估標準化方法，尤其是對基于板的數據集。全長scRNA-seq數據可以使用bulk轉錄組的方法校正基因長度的影響。

• 基因表達歸一化（scale）到均值為0和單位方差尚無共識。我們傾向于不進行scale操作。

• 標準化后的數據應進行log(x+1)轉換，以使數據更符合正態分布，滿足下游分析方法的初始假設。

數據校正和整合

如前所述，數據標準化是去除實驗過程中隨機性的影響 (count sampling)。但是，標準化后的數據可能仍然包含有與研究不相干的因素帶來的影響。數據校正的目的就是進一步去除技術因素和非關注的生物學混雜因素，例如批次、dropout或細胞周期不同帶來的影響 (如何火眼金睛鑒定那些單細胞轉錄組中的混雜因素)。需要注意的是這些混淆因素并不總是需要校正。相反，考慮在分析中移除哪些因素取決于下游分析目的。我們建議分別考慮對生物學和技術混雜因素 (covariates)的校正，因為它們用于不同的目的并代表不同的分析挑戰。

去除與研究不相干的生物因素的影響

盡管校正實驗中技術因素的影響可能對揭示潛在的生物學信號至關重要，但對不關注的生物學因素的校正可用來挑選出特定的感興趣的目標生物學信息。最常見的生物因素校正是消除細胞周期對轉錄組中基因表達的影響。可以使用Scanpy和Seurat對每個細胞的細胞周期評分進行簡單的線性回歸校正或通過應用了更復雜的混合模型的專用程序包如scLVM或f-scLVM進行校正。用于計算細胞周期評分的標記基因列表可在文獻中獲取 (Seurat亮點之細胞周期評分和回歸)。這些方法還可用于校正其他已知的生物因素的影響，例如線粒體基因表達，通常是細胞應激的標識。

在校正特定的生物因素影響之前，應考慮幾個方面。首先，校正生物學因素影響并不總是有助于解釋scRNA-seq數據。雖然消除細胞周期影響可以改善對發育軌跡的推斷，但細胞周期信號也可以提供有意義的生物學信息。例如，可以根據細胞周期評分確定增殖細胞群 (在Github的案例研究中有講)。同樣，必須根據具體研究的對象理解生物學信號。假設生物過程發生在同一有機體內，則這些過程之間可能存在依賴性。因此，校正一個過程的影響可能會無意間掩蓋另一個過程的信號。另外，也有研究者提出細胞大小對轉錄組的影響通常也歸因于細胞周期不同。因此，通過標準化或專用工具（例如cgCorrect）在校正細胞大小的同時也部分校正了scRNA-seq數據中細胞周期的影響。

去除技術影響

用于移除不相干生物因素影響的回歸模型也可應用于校正實驗技術因素的影響。單細胞數據中最突出的技術影響因素是測序深度和批次。盡管標準化可以使得細胞之間的基因定量數據可比，但測序深度的影響仍保留在數據中。這種測序深度效應既可以是生物學自身的影響，也可能是技術帶來的影響。例如，細胞大小可能不同，因此mRNA分子數量也可能不同。而且，標準化后技術帶來的計數影響可能依然存在，因為沒有辦法推斷由于實驗過程中的隨機性帶來的未檢測到的基因的表達。對測序深度進行回歸校正可以提高軌跡推斷算法的性能，該算法依賴于查找細胞之間的過渡（具體見Github上的案例研究）。在校正多個協變量（干擾因素）（例如細胞周期和測序深度）時，應在單個步驟中對所有協變量（干擾因素）執行回歸，以解決協變量（干擾因素）之間的依賴性。

消除測序深度影響的另一個策略是使用更嚴格的標準化程序，例如downsampling或非線性歸一化方法 (見前面標準化部分)。這些方法可能更適用于基于板的scRNA-seq數據集，因為細胞之間較大的測序深度差異可能掩蓋細胞之間的異質性。

批次效應和數據整合

當將細胞分組操作時可能會帶來批次效應（DESeq2差異基因分析和批次效應移除），比如不同芯片上的細胞、不同測序通道中的細胞或在不同時間點收集的細胞都歸類于不同的組。實驗操作過程中細胞所經歷的不同環境可能會影響轉錄組的測量結果或甚至影響細胞自身的轉錄變化。所產生的影響存在多個層面：同一實驗不同的細胞組、同一實驗室的不同實驗或不同實驗室的數據集之間。在這里，我們把第一種情況與后面兩種情況區分開。校正同一實驗中樣品或細胞之間的批次效應是bulk RNA測序批次效應的一種經典方案。我們將其與整合來自多個實驗的數據（稱為數據整合）區分開。通常批次效應校正使用線性方法，而非線性方法則用于數據整合。

最近對經典批次校正方法的比較顯示，ComBat在中低到中復雜度的單細胞實驗中也表現良好 (收藏 Combat開發者、北大生信平臺” 單細胞分析、染色質分析” 視頻和PPT分享)。ComBat構建了基因表達的線性模型，其中批次貢獻在數據的均值和方差中均得到校正（圖3）。與采用什么計算方法無關，批次校正的最佳方法是通過巧妙的實驗設計來避免存在不同批次。通過合并不同實驗條件和樣品的細胞一起進行后續操作可以避免批次效應。使用細胞標記等策略或根據基因組變異信息，可以對實驗中合并的細胞進行拆分。（生信寶典注*：最不建議的實驗設計方式是對照組樣品是一個批次，處理組樣品是一個批次；如果這么做，將不能確定差異是來源于批次，還是真的存在生物差異。*）

圖3. 批次校正前后UMAP可視化結果展示
細胞按其來源的樣本上色。在批次校正前可以看到顏色區分很明顯（細胞所來源的樣品對細胞分群影響大），表明批次影響比較大。使用ComBat校正批次后，同一個顏色分布比較散了，細胞所來源的樣品對細胞分群影響變小了。

與批次校正相比，數據整合面臨的另一個挑戰是數據集之間的組成差異 (compositional differences)。估計批次效應時，ComBat使用同一批次的所有細胞來擬合批次校正參數。這種方法將批次效應與數據集之間非共有的細胞類型或狀態之間的生物學差異混淆。數據整合方法（Cell 深度| 一套普遍適用于各類單細胞測序數據集的錨定整合方案），例如典型相關分析（Canonical Correlation Analysis,?CCA），相互最近鄰（Mutual Nearest Neighbours,?MNN),?Scanorama,?RISC,?scGen,?LIGER,?BBKNN和Harmony已經開發用于解決這個問題。盡管數據整合方法也可應用于簡單的批次校正處理，但是鑒于非線性數據整合方法會增大自由度，使用時需要注意過度校正問題。例如，在較簡單的批次校正設置中，ComBat的表現優于MNN（Buttner et al.，2019）。需要進一步進行數據整合和批次校正方法之間的比較研究，才可以評估這些方法的適用范圍。

缺失值填充

技術數據校正的另一種類型是缺失值填充（也稱為降噪或插補,?denoising or imputation）。單細胞轉錄組的數據包含各種噪聲。這種噪音的一個特別突出的來源是dropout。推斷dropouts事件，用推斷出的合適的表達值替換這些零以減少數據集中的噪聲成為幾種最新工具的目標 ?(MAGIC,?DCA,?scVI,?SAVE,?scImpute)。已證明進行缺失值填充可改善基因與基因相關性的估計。此外，這一步也可以與標準化、批次校正和其他下游分析整合，就像在scVI工具中實現的那樣。盡管大多數數據校正方法都將標準化后的數據作為輸入，但是某些缺失值填充方法是基于預期的負二項噪聲分布，需要基于原始計數矩陣進行操作。在應用缺失值填充時，應考慮到沒有一種方法是完美的。因此，任何方法都可能會對數據中的噪聲進行過高校正或校正不足。確實，已有報道表明缺失值填充可能引入錯誤的相關信號。鑒于在實際應用中難以評估缺失值填充是否得當，用戶選擇是否應用這一方法也是很大的挑戰。當前缺失值填充方法是否能拓展應用到大數據集還是一個問題。鑒于這些考慮，目前尚無關于應如何使用缺失值填充的共識。謹慎的方法是僅在視覺展示數據時使用缺失值填充，而非在探索性數據分析過程中基于填充的數據作出推論或假設。全面的實驗驗證在這里尤為重要。

陷阱和建議：

• 僅在進行軌跡推斷和校正的生物學混雜因素不影響感興趣的生物學過程時才校正這些因素的影響。

• 如果校正的話，所有因素同時校正而不是分別校正技術和非關注的生物因素變量。

• 基于板的數據集預處理時需要校正count數的影響，建議采用非線性標準化方法或downsampling方法進行標準化。

• 當批次之間的細胞類型和狀態組成一致時，建議通過ComBat執行批次校正。

• 數據整合和批次校正應通過不同的方法進行。數據整合工具可能會過度校正簡單的批次效應。

• 用戶需要對只在缺失值填充后才能發現的信號格外注意。探索性分析時最好不進行缺失值填充操作。

特征選擇、降維和可視化

人單細胞RNA-seq數據集可包含多達25,000個基因的表達值。對于一個給定的scRNA-seq數據集，其中有許多基因都不能提供有用信息，并且大多只包含零計數。即使在QC步驟中濾除了這些零計數基因后，單細胞數據集的特征空間也可能超過15,000個維度（即還會剩余15,000多基因）。為了減輕下游分析工具的計算負擔、減少數據中的噪聲并方便數據可視化，可以使用多種方法來對數據集進行降維。

特征選擇

scRNA-seq數據集降維的第一步通常是特征選擇。在此步驟中，對數據集基因進行過濾僅保留對數據的變異性具有信息貢獻的基因（在數據中變異大的基因）。這些基因通常被定義為高變化基因（HVG，highly variable genes）。根據任務和數據集的復雜性，通常選擇1,000到5,000個HVG用于下游分析。Klein et al.的初步結果表明，下游分析對HVG的數量不太敏感。在HVG數量從200到2,400之間選擇不同的數目時，評估顯示PCA結果相差不大?；诖私Y果，我們寧愿選擇更多的HVG用于下游分析 （err on the side of：可以借鑒的英語短語）。

圖EV1. 不同大小數據集選擇HVG基因數量的分布展現。**（數據來源于手動整理的最近發表的多篇scRNA-seq文章，顏色代表不同的技術方式）

在Scanpy和Seurat中都實現了一種簡單而流行的選擇HVG的方法。在這里，基因按其均值表達進行分組，將每個組內方差/均值比最高的基因選為每個分組的HVG。該算法在不同軟件中輸入不同，Seurat需要原始count data；Cell Ranger需要對數轉換的數據。最佳地，應在技術等干擾因素校正之后選擇HVG，以避免選擇僅由于例如批次效應等引入的高可變基因。Yip et al.綜述了其他HVG選擇方法。

特征選擇后，可以通過專用的降維算法進一步對單細胞表達矩陣進行降維。這些算法將表達式矩陣映射到低維空間中，同時以盡可能少的維數捕獲數據中所有的信息。鑒于單細胞RNA測序數據固有的低維性特征，這一方法是合適的。也就是說，細胞表達圖譜構成的生物形態 (biological manifold)可以使用遠少于基因數目的維度信息來展示。降維旨在找到這些維度。

降維有兩個主要目標：可視化和信息匯總（summarization）?？梢暬菄L試在二維或三維空間最優地展示數據集。降維后的維度值就是數據在新的空間進行可視化如繪制散點圖時的坐標值。信息匯總沒有規定輸出的維數；但更高的維數對表示原有數據的差異越來越不重要（生信寶典注*：可以理解為PCA中各個主成分對于原始數據差異的解釋依次降低*）。匯總技術可通過計算數據的固有維數來將數據降維到基本組成（主）成分，從而有助于下游分析。雖然不應使用二維可視化數據來匯總數據集，但匯總方法獲得的降維數據可用來可視化數據。另外，專用的可視化技術通常會更好地展示數據的原始結構和變異性 （還在用PCA降維？快學學大牛最愛的t-SNE算法吧, 附Python/R代碼）。

降維后的維度是通過對基因表達向量進行線性或非線性組合生成的 (PCA主成分分析實戰和可視化 附R代碼和測試數據)。特別是在非線性情況下，降維后的數據難以解釋其生物含義。圖4中顯示了一些常用降維方法的示例應用。隨著可供選擇的方法的增加，詳細討論這些方法已超出了本教程的范圍。相反，我們簡要概述了可能有助于用戶在常見的降維方法之間進行選擇時需要考慮的因素。Moon et al.?提供了有關單細胞分析降維的更詳細的綜述 (Manifold learning‐based methods for analyzing single‐cell RNA‐sequencing data. Curr Opin Syst Biol 7: 36–46)。

圖4. scRNA‐seq數據常用的可視化方法。
(A) PCA, (B) t‐SNE, ? diffusion maps, (D) UMAP and (E) A force‐directed graph layout via ForceAtlas2. 細胞根據測序深度上色 (F) 前31個主成分解釋的原始數據的差異。圖中的拐點( “elbow”)用于選擇下游分析相關的主成分，位于PCs 5 和PCs 7之間。

主成分分析（PCA）和diffusion maps是兩種常用的降維方法，在單細胞分析中也很流行。主成分分析是一種線性方法，通過最大化每個可能維度中捕獲的殘差 (residual variance)來進行降維。盡管PCA不能像非線性方法那樣在更少的維度捕獲原始數據更多的信息，但是它是許多當前可用的聚類或軌跡推斷工具的基礎。實際上，PCA通常用作非線性降維方法的預處理步驟。通常，PCA通過其前N個主成分來代表原始數據集，其中N可以通過elbow算法（參見圖4F）或基于置換檢驗的jackstraw方法確定。PCA簡單線性化的優勢是：降維空間中的距離在該空間的所有區域具有一致的解釋。因此，我們可以將感興趣的統計量與主成分進行關聯分析，以評估其重要性。例如，可以將主成分投影到技術干擾協變量上，以研究QC、數據校正和標準化過程的效果（Buttner et al.，2019），或顯示基因在數據集中的重要性 (PCA bi-plot)。diffusion map是一種非線性數據降維技術。由于diffusion component強調數據的轉換，因此主要用于諸如細胞分化之類的連續過程。通常，每個diffusion component（即diffusion map 維度）突出顯示不同細胞群體的異質性。

可視化

可視化時一般使用非線性降維方法（圖4）。scRNA-seq數據可視化的最常見的降維方法是t‐SNE?（?t‐distributed stochastic neighbour embedding） （還在用PCA降維？快學學大牛最愛的t-SNE算法吧, 附Python/R代碼）。t‐SNE的維度著重于以犧牲全局結構為代價來保留局部相似性 (生信寶典注*：PCA則是盡可能多的保留全局差異*)。因此，這些可視化可能會夸大細胞群體之間的差異，而忽略群體之間的潛在聯系 (t‐SNE dimensions focus on capturing local similarity at the expense of global structure. Thus, these visualizations may exaggerate differences between cell populations and overlook potential connections between these populations)。另一個困難是對參數perplexity parameter的選擇，因為t-SNE圖會因為這個參數值不同而顯示出明顯不同的分簇數。t‐SNE的常見替代方法是UMAP?(Uniform Approximation and Projection method)或基于圖的工具，例如SPRING。UMAP和SPRING的力導向布局算法ForceAtlas2(*force‐directed layout algorithm*)可以說是數據潛在拓撲結構展示最好的方法（Wolf et al，2019）。在此比較中，使UMAP與眾不同的是它的速度快和能應用于更大規模數據的能力（Becht et al，2018）。因此，在沒有特定生物學問題限制的情況下，我們將UMAP視為探索性數據可視化分析的最佳實踐。此外，UMAP還可以把數據降維到二維以上的新數據。盡管我們不知道UMAP`在數據匯總中的任何應用，但它可能是PCA的一個合適的替代方案。

在細胞水平上經典可視化方法的替代方法是PAGA?(partition‐based graph abstraction）。事實證明，此工具可以充分展示數據的拓撲結構，同時使用聚簇生成粗粒度的可視化圖像。結合以上任何一種可視化方法，PAGA可以生成粗粒度的可視化圖像 （coarse‐grained visualizations），從而可以簡化單細胞數據的解釋，尤其是在測序細胞量大或整合了大量細胞的情況下。

陷阱和建議：

• 我們建議根據數據集的復雜程度選擇1,000至5,000個高可變基因 (HVG)。

• 當將基因表達值歸一化為均值為0和單位方差時，或者將模型擬合的殘差用作標準化表達值時，不能使用基于基因表達均-方差的特征選擇方法。因此，在選擇HVG之前，必須考慮要執行哪些預處理。

• 信息匯總 (summarization，類比于挑選重要的主成分)和可視化應使用不同的降維方法。

• 我們建議使用UMAP進行探索性分析可視化；使用PCA做為通用數據降維方法；diffusion maps可以在軌跡推斷時替代PCA。

• UMAP+PAGA是可視化特別復雜的數據集的合適替代方法。

數據預處理步驟

雖然我們在上面按順序概述了scRNA‐seq中常見的預處理步驟，但下游分析通常需要不同級別的預處理數據，建議根據下游應用調整預處理的各個步驟。為了向新用戶說明這種情況，我們將預處理分為五個數據處理階段：（i）原始數據，（ii）標準化后的數據，（iii）校正后的數據，（iv）特征選擇后的數據，以及（v）降維后的數據。數據處理的這些階段可以歸類為三個預處理層：測量的數據，校正的數據和降維的數據。細胞和基因的質量控制篩選是必須執行的步驟，因此在此處略過。盡管層的順序代表了scRNA-seq分析的典型工作流程，但也可以跳過某一步或在處理階段的順序中稍作更改。例如，單批次數據集可能不需要批次校正。在表1中，我們總結了預處理數據的每一層所適用的下游處理程序。

表1

表1中的預處理階段分為三組：測量的數據，校正后的數據和降維后的數據。我們將測量的數據定義為原始數據和處理后但保留了表達值為零的基因的數據。應用細胞特異的size factor的全局標準化方法即使在log(x+1)轉換后仍保留零表達值。相比之下，校正后的數據會去除不需要的變化信息進而替換零表達值。校正后的數據層表示數據的“最干凈”版本，它是數據代表的生物信號的最近似值。我們稱最終預處理層為降維后的數據（reduced data）。此數據層強調數據的主要差異，可以使用降維后的特征集來描述。

前述特性決定了預處理數據對特定下游應用的適用性。作為最后的預處理階段，降維后的數據是廣泛應用的數據層。但是，差異表達分析只能在基因空間內進行生物學解釋，而無法（或無法完全）體現在降維后的數據中。降維后數據的優勢在于生物數據信息的匯總(summarization)和影響生物數據的噪聲的降低。因此，降維后的數據應用于后續的探索性方法（如可視化、鄰域圖推斷、聚類）以及計算復雜的下游分析工具（如軌跡推斷）。實際上，許多軌跡推斷方法在工具本身中都包含了降維功能。

單個基因的表達譜只能在基因空間中進行比較，這一信息存在于測量的數據和校正后的數據中。表達譜可以進行視覺和統計學比較。我們認為視覺比較和統計比較應該在不同的數據層上進行。對于目測基因表達，校正的數據是最合適的，有助于結果解釋。如果對原始數據進行視覺比較，則要求用戶牢記數據中的偏差以準確解釋結果。但是，此處應分別考慮技術和生物協變量的校正數據。雖然對生物協變量的校正可能會增加特定生物信號的強度，但其對潛在的生物學意義表示可能不甚準確，并且會掩蓋可能相關的其他生物信號。因此，經過生物學不相關因素校正的數據主要適用于專注于特定生物學過程（例如軌跡推斷方法）的分析工具。

基因表達的統計比較需要基于測量數據層。去噪、校正批次或其他噪聲源的方法總不是太完美。因此，數據校正方法不可避免地會對數據進行過度或不足校正，因此會以意想不到的方式改變至少某些基因表達譜的方差?；虮磉_的統計檢驗依賴于將背景方差評估為數據噪聲的空模型。隨著數據校正趨于減少背景變化（圖EV2），背景變化被數據校正方法過度校正的基因將更有可能被評估為顯著差異表達基因。此外，某些數據校正方法（例如ComBat）將不吻合于實驗設計的表達信號定義為噪聲，隨后將其從數據中刪除。這一基于實驗設計的優化方法除了可能造成數據中的噪聲被低估還會導致對效應大小 （effect size）的高估。鑒于這些考慮，與校正后的數據相比，基于測量數據進行差異檢驗是一種更為保守的方法。使用測量的數據時，可以并且應該在差異統計分析模型中考慮技術協變量 （參考DESeq2差異基因分析和批次效應移除）。

圖EV2. 批次校正和降噪后變異系數的變化
負值代表數據校正后變異系數降低。第一行展示的是ComBat校正(A) mouse intestinal epithelium (mIE) 和 (B) mouse embryonic stem cell (mESC) 數據變異系數的變化。第二行展示的是DCA降噪后? mIE 和 (D) mESC 數據變異系數的變化。

最近一篇單細胞差異表達分析方法的文章也支持以上觀點, 它只使用原始或標準化后的矩陣作為輸入 (Soneson & Robinson, 2018)。這篇研究中標準化后的數據僅限于全局標準化方法。然而當前可用的非線性標準化方法模糊了標準化和數據校正的界限（具體見標準化部分）。這樣標準化之后的數據已不適合用于差異基因分析。

陷阱和建議：

• 原始測量的數據用于差異基因統計檢驗分析；

• 校正后的數據用于數據的可視化比較；

• 降維后的數據用于其他基于數據形態（biological data manifold）的下游分析。

下游分析

經過預處理后，我們稱之為下游分析的方法指應用于生物學發現并描述潛在的生物學系統的方法。通過將可以解釋的模型擬合到數據中獲得相應的結論，比如有相似基因表達譜的細胞群代表一個細胞簇、相似細胞之間基因表達的微小變化指示連續（分化）軌跡；或具有相關表達趨勢的基因指示共調控等。

下游分析可分為細胞水平和基因水平的方法，如圖5所示。細胞水平分析通常著重于兩種結構的描述：簇和軌跡。這些結構又可以在細胞和基因水平上進行分析，即簇分析和軌跡分析方法。廣義地講，簇分析方法試圖將細胞分類為離散的多個組來解釋數據的異質性。相反，在軌跡分析中，數據被視為細胞連續變化動態過程的一個個快照，軌跡分析方法研究了這一連續變化過程。

圖5. 下游分析方法總覽（藍色背景的方法是基因水平的分析方法）。

在此，我們在詳細描述獨立于這些細胞結構進行的基因水平分析之前，先描述細胞水平和基因水平的聚類和軌跡分析工具。

聚類分析

聚類。將細胞聚類成簇通常是任何單細胞分析的第一個中間結果。聚類成簇使我們可以推斷成員細胞的身份。簇是通過基于細胞基因表達譜的相似性將細胞分組得到的。表達譜相似性是通過對將降維的數據進行距離度量確定的。相似性評分的一個常見示例是在主成分降維后的表達空間上計算的歐氏距離。存在兩種根據這些相似性分數生成細胞簇的方法：聚類算法和基于圖的社群檢測方法?（community detection methods）。

聚類是直接基于距離矩陣的經典無監督機器學習問題。通過最小化簇內距離或在降維后的表達空間中鑒定密集區域，將細胞分配給簇。流行的k均值聚類算法通過確定聚類質心并將細胞分配給最近的聚類質心，將細胞分為k個類。質心位置經過迭代優化。這種方法需要輸入預期的簇數，但這通常是未知的，必須進行啟發式校準。k均值在單細胞數據中的應用因所使用的距離度量而異。標準距離度量是使用歐氏距離，替代距離包括余弦相似度、基于相關的距離度量或SIMLR方法，該方法使用Gaussian kernels為每個數據集計算距離度量。最近的比較表明，與k均值一起使用或作為Gaussian kernels分析的基礎時，基于相關的距離可能會勝過其他距離度量方法。（基因共表達聚類分析及可視化）

社群檢測方法是一種圖分類算法，依賴于單細胞數據的網絡圖表示。圖是使用K最近鄰方法（KNN圖）獲得的。細胞在圖中表示為節點。每個細胞與其K個最相似的細胞相連，通常對主成分降維后的數據計算歐氏距離作為相似性度量。根據數據集的大小，K通常設置為5到100個最近的鄰居。獲得的KNN結果圖捕獲了表達圖譜的基礎拓撲結構（Wolf et al., 2019）。表達譜相似的細胞集合對應于圖的緊密連接區域，然后使用社群檢測方法檢測圖中這些緊密連接區域。社群檢測方法通常比聚類方法要快，因為只有相鄰的細胞對會被視為屬于同一群集。因此，這種方法大大減少了可能的細胞簇的搜索空間。

在開創性的PhenoGraph方法（Levine et al.，2015）出現之后，聚類單細胞數據集的標準方法已演變成多分辨率模塊優化算法 (multi‐resolution modularity optimization），即在單細胞KNN圖中應用Louvain算法。此方法是在Scanpy和Seurat單細胞分析平臺中應用的默認聚類方法。評估發現，這一方法應用于單細胞RNA測序數據以及流式細胞和質譜數據分析時優于其他聚類方法。從概念上講，Louvin算法將組內細胞連接數高于預期的一組細胞視作一個簇。(生信寶典注*：**Louvain algorithm算法采用迭代方式，先計算每個點加入到其鄰居節點所在社區帶來的模塊性評估收益，然后將其加入收益最大的鄰居節點所在社區，對所有未歸類點重復進行這一過程，直至結果穩定。然后再把每個社區作為一個節點，重復上一步。*) 優化后的模塊鑒定函數包括一個分辨率參數，該參數允許用戶確定簇的近似大小。通過對KNN圖取子集，還可以只聚類一部分特殊的細胞簇。這種子集聚簇模式可以允許用戶識別某個細胞類型簇內的細胞狀態，但也可能會發現僅由數據中的噪聲帶來的聚類模式。

陷阱和建議：

• 我們建議在單細胞KNN圖上執行Louvain社群檢測算法獲得聚類簇。

• 細胞聚簇不一定使用單一分辨率進行。對特定細胞簇進行細化聚類是鑒定數據集中亞結構的有效方法，可以發現細胞亞群。

細胞簇注釋

在基因水平上，對每個簇的分析是鑒定其標記基因 (Marker genes)。這些所謂的標記基因代表了細胞簇的特征，用于給細胞簇一個有生物學意義的標簽（Celaref | 單細胞測序細胞類型注釋工具）。該標簽表示集群中細胞的身份。由于任何聚類算法都會聚類出細胞簇，因此聚類獲得的生物簇的準確性只能通過其生物學注釋進行衡量?(生信寶典注*：**這也是前面和易生信課程中反復強調的，細胞過濾時標準盡量松一些，根據聚類結果回看之前的參數設置是否合理。*)。

盡管把單細胞數據中鑒定到的簇歸類為一個有生物意義的細胞類型是一個誘人的結果，但是細胞身份的確定不是那么簡單（Wagner et al., 2016；Clevers et al., 2017）。首先，并不總是清楚細胞類型定義的精確程度。例如，“T細胞”可能是某些人滿意的細胞類型標記，但其他人可能會在數據集中尋找T細胞亞型并區分CD4+T細胞和CD8+T細胞?。另外，同一細胞類型的不同狀態可能會被分到不同的細胞簇?；谏鲜鲈?#xff0c;最好使用術語cell identities而不是cell types。在對細胞進行聚類和注釋之前，用戶必須確定要注釋到的細胞信息精確水平，從而確定合適的聚簇分辨率。

識別和注釋細胞簇依賴于外部信息，這些信息描述了每個細胞類群預期的表達譜。隨著小鼠腦圖集（Zeisel et al., 2018）、人類細胞圖集（Regev et al., 2017）和其他工作的繼續，可用的參考數據庫越來越多。這些數據庫極大地方便了細胞簇注釋。在沒有相關參考數據庫的情況下，可以通過將細胞聚類簇的Marker基因與文獻中的Marker基因進行比較來注釋細胞身份（請參閱Github上的案例研究）或通過直接對文獻報道的標記基因在自己研究的數據中進行表達值可視化來確定細胞身份（圖6B）。應當注意的是，后一種方法將用戶限制于普通轉錄組數據鑒定出的經典細胞類型，而不是對細胞身份 (cell identities)的理解。而且，研究表明常用的細胞表面標志基因在定義細胞身份上是效果有限的（Tabula Muris Consortium et al., 2018）。

圖6. 聚類分析結果

(A) UMAP可視化Louvain算法鑒定的細胞簇，并標記注釋的細胞類型名稱。(B) 通過細胞類型Marker基因的表達標記細胞簇所代表的細胞類型，如干細胞 (Slc12a2), 腸上皮細胞 (Arg2), 杯狀細胞 (Tff3) 和潘氏細胞 (Defa24)?；疑淼捅磉_，紅色代表高表達。Marker基因也可能在其它細胞簇有表達，比如Tff3在杯狀細胞和潘氏細胞都有表達，只是在杯狀細胞更高。? 腸上皮細胞區域細胞類型組成熱圖（上面是近端區域，下面是遠端區域）。相對細胞密度高的區域用深紅色表示。

有兩種方法可以使用參考數據庫信息注釋細胞簇：使用計算出的Marker基因或使用完整的基因表達譜。可以通過對目標簇的細胞和數據集中的所有其他細胞進行差異表達（DE）分析來鑒定標記基因集。通常，我們關注在目標簇中上調的基因。由于預期標記基因表達變化幅度較大，因此通常使用簡單的統計檢驗（例如Wilcoxon秩和檢驗或t檢驗）就可以對兩組細胞的基因進行差異檢驗分析。根據統計檢驗結果，選擇排名最高的基因視為標記基因?？梢酝ㄟ^比較自己的數據中鑒定出的標記基因和來自參考數據集的標記基因對細胞簇進行注釋。一些在線工具如www.mousebrain.org (Zeisel et al, 2018) 或http://dropviz.org/ (Saunders et al, 2018)允許用戶在參考數據集中可視化用戶鑒定出的Marker基因從而方便細胞簇的注釋。

鑒定標記基因時應注意兩個方面。首先，用于鑒定標記基因的P值基于以下假設：即鑒定出的細胞簇是有生物學意義的。如果細胞簇鑒定過程中存在不確定性，則必須在統計檢驗中考慮簇分配與標記基因鑒定之間的關系。這個關系起因于鑒定細胞簇和標記基因都是基于同一套基因表達數據。差異基因檢測的零假設（null hypothesis）是兩組細胞整體基因的表達值具有相同的分布。然而，由于這兩個聚類組是基于基因表達變化的聚類結果得到的，其基因表達譜從本質上肯定存在差異。因此即使在應用splatter隨機生成的數據的聚類結果中也可以鑒定出顯著的Marker基因 (Zappia et al, 2017)。為了在聚類數據中獲得合理的顯著性度量，可以使用置換檢驗減少聚類步驟帶來的影響。補充介紹S3中對這一檢驗有詳細描述。最近的一個差異表達工具也專門解決了這個問題（Preprint：Zhang et al，2018）。

在當前的研究中，P值通常會被夸大，這可能導致對標記基因數量的高估。但是，基于P值對基因進行排序則不受影響。假設聚類結果是有生物學意義的，那么排名最高的標記基因仍將是最佳候選標記基因。首先，我們可以通過可視化展示粗略地查驗獲得的標記基因。我們強調，特別是通過無監督聚類方法定義細胞聚類簇時，會導致夸大的P值。當改為通過單個基因的表達確定細胞簇的身份時，P值可以解釋為相對于其他基因的期望值。盡管很常見，但是這種單變量的聚類簇注釋方法除了在特定情況下不建議使用（例如，β細胞中的胰島素或紅細胞中的血紅蛋白）。其次，標記基因在數據集中將一個細胞簇與其他細胞簇區分開，不僅依賴于細胞簇鑒定結果，還依賴于數據集的組成。如果數據集組成不能準確表示背景基因表達，則檢測到的標記基因將偏向在其它細胞中未檢測到的基因。特別是在細胞多樣性低的數據集中計算標記基因時，必須特別考慮這一方面。

最近，自動細胞簇注釋已可用。通過直接將注釋好的參考簇的基因表達譜與單個細胞進行比較，使用諸如scmap（Kiselev et al., 2018b）或Garnett（Preprint：Pliner et al., 2019）之類的工具可以在參考集和自己的數據集之間比較注釋。因此，這些方法可以同時執行細胞注釋和細胞簇鑒定，而無需數據驅動的細胞聚類步驟。但由于不同實驗條件下細胞類型和狀態組成不同（Segerstolpe et al.，2016;Tanay＆Regev，2017），基于參考數據的聚類不應取代數據驅動的聚類過程。

聚類、聚類注釋、重新聚類或子聚類以及重新注釋的迭代是很耗時的過程。自動化的細胞簇注釋方法大大加快了此過程。但是，自動和手動方法各自存在優點和局限性，速度的提高與靈活性的降低并存，因此很難有一種全局適應的方法。如前所述，參考數據集很難與正在研究的數據集包含完全一樣的細胞類型。因此，不應放棄用于手動注釋標記基因的計算。特別是對于包含許多細胞簇的大型數據集，當前的最佳實踐是兩種方法（自動+手動）的組合。為了提高速度，可以使用自動化細胞類型注釋工具粗略地標記細胞并確定可能需要鑒定細胞子群的簇。隨后，應針對細胞簇計算標記基因，并將其與參考數據集或文獻中的已知標記基因集進行比較。對于較小的數據集和缺少參考數據庫的數據集，手動注釋就足夠了。

陷阱和建議：

• 不要使用標記基因的P-value值來驗證細胞簇身份，尤其是當檢測到的標記基因無益于注釋細胞簇時。P值可能會被夸大。

• 請注意，由于數據集的細胞類型和狀態組成，同一細胞簇的標記基因可能在數據集之間有所不同。

• 如果存在相關的參考數據集，我們建議綜合使用自動聚類注釋和基于數據的標記基因的手動注釋來注釋細胞類型。

細胞組成分析?（Compositional analysis）。在細胞層面，我們可以根據其組成結構分析聚類數據。細胞組成數據分析圍繞不同樣品落入每個細胞簇的細胞比例進行分析。這一比例可能在疾病狀態下發生變化。例如，研究表明沙門氏菌感染會增加小鼠腸上皮區域腸上皮細胞（enterocytes）的比例（Haber et al.，2017）。

要研究單細胞數據的細胞組成變化，需要足夠的細胞數量來穩健地評估細胞簇的比例，并需要足夠的樣本數量來評估細胞簇組成中的背景變化。由于合適的數據集直到最近才可用，因此專用工具尚待開發。在上述小鼠研究中，細胞身份計數使用泊松分布建模，使用實驗條件作為協變量和檢測到的細胞總數作為偏移量。在此，可以對回歸系數進行統計檢驗，以評估特定細胞群體的出現頻率是否發生了顯著變化。但是，對同一數據集中其他細胞簇的統計檢驗分析并非彼此獨立。如果一個細胞聚類簇的比例發生變化，那么所有其他細胞聚類簇的比例也會發生變化。因此，使用該模型無法評估整體構成是否發生了顯著變化。在沒有專用工具的情況下，細胞構成數據的可視化展示會給不同樣品之間細胞組成變化提供有用的信息（圖6C）。該領域的未來發展可能會借鑒質譜流式細胞技術（Mass Cytometry）（例如Tibshirani et al., 2002；Arvaniti＆Claassen., 2017；Lun et al., 2017；Weber et al., 2018）或微生物組的方法（Gloor et al., 2017），組成數據分析在這些領域已經受到了更多關注。

陷阱和建議：

• 統計檢驗分析時需要考慮到樣本之間的細胞簇比例的變化是互相依賴的（非獨立的）

Trajectory analysis

**軌跡推斷。**細胞多樣性不能通過離散的分類系統（例如細胞聚類）充分描述。驅動觀察到的細胞異質性發展的生物進程是一個連續過程（Tanay＆Regev，2017）。因此，為了捕獲細胞身份之間的過渡狀態、不同的分化分支或生物學功能的漸進式非同步變化，我們需要動態的基因表達模型。這類方法稱為軌跡推斷（trajectory inference，TI）。

軌跡推斷方法將單細胞數據視為連續過程的一個個快照（教你如何定義新亞群?| 在單細胞水平上解析人肝硬化的纖維化微環境）。這一過程通過最小化相鄰細胞之間的轉錄改變構建細胞空間的轉換路徑（圖7A和B）。這些路徑上的細胞排序由偽時間變量 （pseudotime variable）描述。雖然此變量是基于距離根細胞的轉錄距離計算的，但通常被解釋為發育時間的代名詞（Moignard et al.，2015; Haghverdi et al.，2016; Fischer et al.，2018; Griffiths et al.，2018）。

圖7. 小鼠腸上皮細胞軌跡分析結果展示

(A) Slingshot鑒定出的遠端和近端腸上皮細胞分化軌跡。遠端譜系細胞按pseudotime從紅色到藍色上色。數據集中其它細胞用灰色點展示。(B) PCA空間中Slingshot 軌跡和細胞簇比較展示。簇名字簡寫如下：cuEP—enterocyte progenitors;?Imm. Ent.—immature enterocytes;?Mat. Ent.—mature enterocytes;?Prox.—proximal;?Dist.—distal. ? 遠端腸上皮細胞軌跡中每個pseudotime對應的細胞密度分布。柱子按每個pseudotime 區間主要的細胞簇上色。(D) 數據集投射到UMAP空間的抽象圖展示。每個簇用不同顏色點展示。出現在其它軌跡中的簇特別標記處理用于比較分析。TA代表短暫增殖下撥 (transit amplifying cells)。(E)R包GAM展示腸上皮細胞中基因表達隨pseudotime的動態變化。

自Monocle（Trapnell et al. 2014）和Wanderlust（Bendall et al. 2014）建立了TI?（trajectory inference）領域以來，可用的TI方法數量激增。當前可用的TI方法的差別在于構建的發育軌跡模型拓撲結構復雜性不同，從簡單的線性軌跡或二分支軌跡到復雜樹形軌跡、多分支軌跡或組合多種拓撲結構軌跡。在最近對TI方法進行的全面比較中（Saelens et al.，2018）發現沒有一種單獨的方法可以在所有類型的軌跡分析中都表現最優 （NBT|45種單細胞軌跡推斷方法比較，110個實際數據集和229個合成數據集）。相反，應根據預期軌跡的復雜性選擇TI（軌跡分析）方法，研究比較表明Slingshot（Street et al.，2018）在簡單軌跡分析如從線性軌跡到二分支和多分軌跡表現最佳。如果預期數據對應更復雜的軌跡，作者建議使用PAGA（Wolf et al.，2019）。如果知道精確的軌跡模型，則可以選擇使用更特定的方法來提高性能（Saelenset al.，2018）。通常，應使用多種方法來確定評估推斷出的軌跡，以避免方法偏差。

在典型的分析流程中，軌跡推斷（TI）方法應用于降維后的數據。如果使用的TI工具自帶了降維功能，則基于校正后的數據進行分析。由于通常在細胞內同時發生多種生物學過程，因此消除其他生物過程的影響對鑒定預期軌跡可能很有用。例如，T細胞在成熟過程中可能會經歷細胞周期轉換（Buettner et al.，2015）。此外，由于幾種性能最好的TI方法依賴于聚類后的數據，因此TI通常在聚類之后執行。軌跡中的細胞簇可能表示穩態或亞穩態細胞（請參見“亞穩態”；圖7B和C）。隨后，可以將RNA velocities?（RNA速度，或RNA表達動力學）疊加到軌跡上確認發育方向（La Manno et al.，2018）。（生信寶典注*：*新生轉錄本成熟過程中需要進行剪接操作。對于一個穩定表達的基因，總會在細胞中找到存在一定比例的未剪接的非成熟RNA形式，用于補充老的轉錄本的降解。如果一個基因剛被激活，短時間內將會有高比例的未成熟轉錄本。相反，當一個基因被抑制時，轉錄過程會早于轉錄本降解過程而被抑制，未成熟轉錄本的比例會降低。因此對于細胞中每個基因，未剪接的mRNA相對于剪接的mRNA的比例（RNA velocity）可以推斷瞬時表達動力學，進一步推演組織內發生的細胞轉變。https://www.nature.com/articles/d41586-018-05882-8）

推斷的軌跡不一定要完全對應生物發育過程。首先，推斷的軌跡僅表示轉錄相似性。很少有TI方法在其模型中包括不確定性評估（Griffiths et al., 2018）。因此，需要更多的信息來驗證是否確實捕獲了生物過程。這些信息可以來源于干擾實驗、推斷的調控基因動力學以及RNA velocity數據的支持等。

陷阱和建議：

• 我們建議使用Saelens et al.（2018）的綜述（NBT|45種單細胞軌跡推斷方法比較，110個實際數據集和229個合成數據集）作為指南 。

• 推斷的軌跡不需要完全對應生物過程， 應該收集更多的證據來解釋軌跡。

**基因表達動力學。**推斷的軌跡也可能是轉錄噪聲擬合的結果，一個排除方式是在基因水平上對鑒定出的軌跡進行進一步分析。在整個偽時間范圍內連續平穩變化的基因是決定軌跡的關鍵基因，可用于識別潛在的生物學過程。此外，這一組與軌跡相關的基因應該包含調控相應生物進程的基因。調節基因可以幫助我們理解該生物過程是如何和為什么被觸發的，并且可能是潛在的藥物靶標（Gashaw et al., 2012）。

早期尋找軌跡相關基因的方法是沿軌跡在細胞簇之間進行差異基因分析（DE）。現在通過將基因表達與偽時間 (pusedotime)信息進行回歸來鑒定隨軌跡顯著變化的基因。為了保證表達值隨偽時間變量平滑變化，首先通過擬合樣條曲線 (fit a spline)或其他局部回歸方法（例如loess）對偽時間進行平滑處理。回歸框架的區別有兩點：一是其假設的噪聲模型不同，二是構建的表達-偽時間函數的類別不同。通過基因對偽時間的依賴性進行模型篩選獲得潛在的調控基因。這種基于偽時間回歸進行的DE基因測試也會受到軌跡推斷方法的影響，就像不同細胞簇之間的差異基因分析會受到聚類方法的影響一樣（參見“細胞簇注釋”部分）。因此，在此步驟中獲得的P值不應視為顯著性的評估。

當前幾乎沒有專用于分析基因表達時間動力學的分析工具。BEAM是整合于Monocle?軌跡推斷方法中的一個工具，可以分析分支特異的基因表達動力學。除了這個，用戶可以選擇使用LineagePulse?(https://github.com/YosefLab/LineagePulse)，這個工具考慮了dropout噪音，但仍然在開發中?；蛘呋跇藴蔙包或Limma包寫作自己的分析測試工具。在Slingshot`的使用幫助中可以找到一個這樣應用的例子供參考。

由于可用的工具很少，因此尚不能確定研究基因時間表達動力學的最佳方法。使用上述所有方法肯定可以對基因表達動力學進行探索性研究。將來，高斯過程 (Gaussian processes)可能會提供一個natural model來研究基因的時間動力學。此外，對調節模塊整體而不是單個基因進行統計檢驗分析可能會提高信噪比并更方便生物學解釋。

**亞穩態。**軌跡的細胞水平分析是指研究偽時間軸上的細胞密度。假設以無偏方式對細胞進行隨機采樣，軌跡中細胞密集的區域表示優選的轉錄狀態。當將軌跡解釋為時間過程時，這些密集區域可能代表例如發育中的亞穩態（Haghverdi et al., 2016）。我們可以通過繪制偽時間坐標的直方圖來找到這些亞穩態（圖7C）。

**細胞水平聯合分析。**聚類和軌跡推斷代表單細胞數據的兩個不同視角，并可以在粗粒度(coarse‐grained)的圖形展示中進行統一。用點代表單細胞簇，軌跡作為簇之間相連的邊，可以同時展示數據的靜態和動態屬性。這一聯合分析在PAGA工具中首先提出。PAGA提出一個統計模型進行細胞簇互作分析，將細胞相似度高于期望值的兩個簇之間用線連接。在最近的綜述中 (Saelens et al, 2018)，PAGA相比于其它軌跡分析方法整體表現更優。PAGA是唯一一個討論到的可以處理不相連的拓撲結構和包含環形結構的復雜軌跡圖的一個工具。這一特征使得PAGA在可視化整個數據集并進行探索分析時很有用。

**基因水平分析。**到目前為止，雖然我們專注于表征細胞結構的基因水平分析方法，但是單細胞數據的基因水平分析有更廣泛的內容。差異表達分析、基因集分析和基因調控網絡分析都可以直接研究數據中的分子信號。這些方法不是描述細胞異質性，而是把這種異質性作為理解基因表達差異的背景。

**差異表達分析。**表達數據分析的一個常見問題是兩個實驗條件之間是否有基因發生了差異表達。差異基因分析源于大量細胞基因表達分析（Scholtens＆von Heydebreck，2005），已經有過很好的描述。相對于bulk差異分析，單細胞的優勢是可以基于細胞簇進行細胞異質性評估。分析某個類型細胞簇在不同實驗條件下的轉錄響應（Kang et al.，2018）。雖然用于解決同樣的問題，普通和單細胞差異基因分析工具在方法上差別很大。普通轉錄組的差異基因分析方法解決的是少量樣品中精確評估基因表達變化的問題，而單細胞數據中不存在這一問題。另一方面，單細胞數據包含特有的技術噪音如dropout和高細胞變異性 (Hicks et al, 2017; Vallejos et al, 2017)。這些影響因素在設計單細胞專用的差異分析方法時都要有考慮。然而，在最近一個大規模的差異基因分析方法評估中發現普通轉錄組差異基因檢測工具與表現最好的單細胞差異分析工具性能相當 (Soneson & Robinson, 2018)。進一步地，當在普通轉錄組的差異分析工具中引入基因權重對單細胞數據建模時，這些工具的表現優于同類的單細胞差異分析工具 (Van den Berge et al, 2018)。基于這個比較分析，DESeq2?(Love et al, 2014) 和EdgeR?(Robinson et al, 2010) 加ZINB‐wave?(Risso et al, 2018)估計出的權重信息表現最好。不過仍需加權的普通轉錄組差異分析的獨立比較研究進一步確認這個結果。

轉錄組差異分析工具的加權建模后效果的改善是以犧牲計算效率為代價的。隨著單細胞數據中細胞量越來越大，算法的運行時間在方法選擇中越來越重要。因此單細胞工具MAST被視作加權普通轉錄組差異分析工具的替代方法。MAST應用hurdle model消除dropout的影響，并構建基因表達相對于實驗條件和技術協變量的變化模型。這是在前述評估研究中表現最好的單細胞差異基因分析方法 (Soneson & Robinson, 2018)，并且在一個單一數據集的小范圍研究中表現優于其它普通或單細胞差異基因分析方法 (Vieth et al, 2017)。MAST相比于加權的普通轉錄組差異分析工具有10-100倍的速度提升，還可以使用limma–voom模型進一步獲得10倍的速度提升。雖然limma是普通轉錄組基因差異分析方法，但limma–voom與MAST性能相當。

因為用于差異基因分析的是標準化后但未進行技術變量和非相關生物變量校正的數據，因此在分析過程中考慮這些混雜變量對獲得穩定的結果是關鍵的。通常差異基因分析工具都允許用戶考慮混雜變量（confounders），但是把哪些變量納入模型需要謹慎評估。比如，在多數單細胞數據集中，樣本和實驗處理條件是互斥的，因為幾乎不可能對一個樣品進行多個實驗處理。如果我們在模型中同時考慮樣品和實驗條件兩個協變量，跟這些協變量相關的變化則不會被明確確定。因此在考慮不同實驗條件下的差異分析時，通常不在模型中引入樣品信息作為協變量。當校正多個確定的批次分類變量時，可視化展示協變量的混雜影響將會變得更困難。在這個情況下，需要先確認模型設計矩陣是否滿秩 （the model design matrix is full rank，full rank是指模型設計矩陣中的每一列都不能通過對其它列的線性組合得來，即它們之間不存在線性相關）。即便一些設計矩陣不是滿秩，差異基因分析工具通常也不會給出警告而是繼續運行。這時獲得的結果將可能不是預期的分析方向。

我們這兒描述的場景中，實驗條件協變量是在實驗設計中決定的。因此在同一簇內基于這一協變量的差異基因分析是獨立于聚類過程的。這將基于實驗條件的差異基因分析和基于細胞簇的差異基因分析區別開來。基于實驗條件的差異基因檢測獲得的P-values是顯著性檢驗統計指標，需要進一步做多重檢驗校正。為了降低多重檢驗校正的壓力，不感興趣的基因通常需要預先排除掉。雖然假基因和非編碼基因也可以提供有用信息，但在分析中也通常會被忽略掉。

陷阱和建議：

• DE測試不應基于校正后的數據（去噪，批次校正等），而應基于原始測量數據并在模型中引入干擾協變量。

• 用戶不應依賴差異基因檢測工具校正帶有混雜協變量的模型。設計矩陣模型應當保證滿秩。

• 我們建議使用MAST或limma進行差異基因分析。

• 39個轉錄組分析工具，120種組合評估(轉錄組分析工具哪家強-導讀版)

• DESeq2差異基因分析和批次效應移除

• 什么？你做的差異基因方法不合適？

**基因集分析。**基因水平的分析通常會獲得一長串難以解釋其生物含義的基因列表。比如，處理組和對照樣品通常會有數以千計差異表達基因。我們可以通過基因共有的特征對基因進行分組并檢驗這些特征在候補基因列表中的出現是否顯著富集 (overrepresented)以便更好的解釋結果。

不同用途的基因集信息可以在校正過的基因注釋數據庫中獲取。為了解釋差異基因結果，通常按共有的生物進程對基因進行分組。生物進程信息存儲于MSigDB、Gene Ontology，通路數據存儲于KEGG和Reactome數據庫?；蛄斜淼墓δ茏⑨尭患治鲇卸喾N工具可用，具體見Huang et al (2009)和?Tarca et al (2013)的綜述。

• 無需寫代碼的高顏值富集分析神器

• 去東方，最好用的在線GO富集分析工具

• 沒錢買KEGG怎么辦？REACTOME開源通路更強大

• 一文掌握GSEA，超詳細教程

• 這個只需一步就可做富集分析的網站還未發表就被CNS等引用超過350次

• 單基因GSEA怎么做？

分析領域富集分析的一個新操作是基于配對基因標記的配體-受體分析。通過受體和共軛配體之間的表達關系推斷細胞簇之間的互作。CellPhoneDB提供了受體-配體對信息，并用統計模型解釋簇之間的高表達基因是否可以配對 (Zepp et al, 2017; Zhou et al, 2017;?Cohen et al, 2018;?Vento‐Tormo et al, 2018)。

**基因調控網絡。**基因不是獨立發揮作用的。相反，基因的表達水平是由與其他基因和小分子之間的復雜調控決定的。揭示這些調控作用是基因調控網絡（GRN）推斷方法的目標（SCENIC | 從單細胞數據推斷基因調控網絡和細胞類型）。

基因調控網絡推斷是基于對基因共表達如相關性、互信息（mutual information）或回歸模型分析進行測量（Chen＆Mar，2018）。如果在考慮了其他所有基因都作為混雜因素后，兩個目標基因之間依然存在共表達關系則可以推斷這兩個基因之間存在因果調控關系。基因調控關系鑒定與軌跡相關的調控基因檢測相關。確實，幾種單細胞GRN推論方法使用的機械微分方程模型中都引入了軌跡信息（Inferring gene regulatory relationships is related to the detection of trajectory‐associated regulatory genes. Indeed, several single‐cell GRN inference methods use trajectories with mechanistic differential equation models）（Ocone et al., 2015; Matsumoto et al., 2017）。

雖然存在專門針對scRNA-seq數據開發的GRN推斷方法（SCONE：Matsumoto et al., 2017;?PIDC：Chan et al., 2017;?SCENIC：Aibar et al., 2017），但最近的比較分析表明普通和單細胞轉錄組的方法在這些數據表現都不好（Chen＆Mar，2018）。GRN推斷方法仍可提供識別生物過程的因果調節子的有用信息，但我們建議謹慎使用這些方法。

陷阱和建議：

• 用戶應警惕推斷的調控關系中的不確定性?；蛘{控關系豐富的基因模塊比調控信息稀疏的模塊更可靠。

**分析平臺。**單細胞分析流程是多個獨立開發的工具的集合。為了促進這些工具之間的數據傳遞，單細胞工具開發時都考慮使用一致的數據格式。這些工具為構建分析流程提供了基礎。目前可用的分析平臺有基于R（McCarthy et al., 2017; Butler et al., 2018）或Python（Wolf et al., 2018）的命令行，或具有圖形用戶界面（GUI）的本地應用程序（Patel，2018; Scholz et al., 2018）或Web服務（Gardeux et al., 2017; ?Zhu et al., 2017）。Zhu et al.（2017）和Zappia et al.（2018）對這些分析平臺進行了綜述。

在命令行平臺中，Scater（McCarthy et al.，2017）和Seurat（Butler et al.，2018）可以輕松地與R Bioconductor項目提供的各種分析工具進行交互（Huber et al.，2015）（讓你的單細胞數據動起來！|iCellR(二)）。Scater在質量控制和數據預處理方面具有特殊優勢，而Seurat可以說是最受歡迎和最全面的分析平臺，提供了大量分析模塊和教程。scanpy是一個基于Python的不斷發展的新的單細胞分析平臺，該平臺可以擴展應用于分析更大規模的單細胞數據集。它可以整合Python寫作的很多工具，尤其是在機器學習中流行的工具。

圖形用戶界面工具可以幫助非專業用戶構建單細胞分析工作流。通常會基于固定的分析流程引導用戶，這些分析流程有助于簡化分析，但同時也限制了用戶的靈活性，在探索性分析時特別有用。諸如Granatum（Zhu et al.，2017）和ASAP（Gardeux et al.，2017）之類的平臺在集成的工具方面有所不同，其中Granatum包含了更多的方法。作為Web服務器，這兩個平臺都直接可用，但是其不同的計算框架將限制它們擴展到大規模數據集的能力。ASAP的測試數據集只包含92個細胞?；赪eb的GUI平臺的替代產品是FASTGenomics（Scholz et al., 2018）、iSEE（Rue-Albrecht et al., 2018）、IS-CellR（Patel，2018）和Granatum等軟件包，它們在本地服務器上運行。這些平臺和GUI封裝工具可以隨著本地電腦計算力的增強進行性能擴展。將來，Human Cell Atlas（https://www.humancellatlas.org/data-sharing）的不斷發展將促進功能更強大的可應用于大規模細胞數據的可視化探索工具的開發。

結論與展望

我們綜述了典型的scRNA‐seq分析流程的各個步驟，并在案例研究教程（https://www.github.com/theislab/single-cell-tutorial）中實現了這些步驟。本教程旨在通過對當前可用分析方法的比較研究確定當前最佳分析實踐流程。雖然匯總單個最佳實踐工具不能保證流程最佳，但我們希望我們的分析流程能夠代表單細胞分析領域最新技術現狀的一個掠影。它為新來者提供了進入該領域的合適切入點，并希望能有助于人類細胞圖譜的scRNA-seq分析流程構建（Regev et al.，2018）。應該注意的是，方法比較必然落后于最新的方法開發。因此，我們提到了可能尚未被獨立評估的新工具。隨著新工具和更好工具的進一步發展和比較研究，文中推薦的每個單獨的工具也需要一直更新，但是有關數據處理階段的一般注意事項應保持不變 （生信寶典注*：**工具一直在變，但核心原理不變。來易生信，掌握核心知識。*）。

特別受關注并且可能改變現有分析流程的兩個開發方向是深度學習分析和單細胞多組學整合分析。由于深度學習可以很容易拓展到大規模數據分析，已經廣泛用于計算機視覺到自然語言處理等多個領域，而且在基因組領域應用也越來越多。機器學習首先應用于scRNA-seq的降維和降噪分析領域，例如scVis：Ding et al.，2018;?scGen：Preprint：Lotfollahi et al.，2018;?DCA：Eraslan et al.，2019）。最近，深度學習已形成完善的分析流程用于擬合數據、對數據進行去噪并在模型框架內執行下游分析如聚類和差異表達等（scVI：Lopez et al., 2018）。在這一分析方法中，可以在保留數據變化的準確估計的同時，把噪聲和批次效應的影響納入下游統計分析模型中。諸如此類的整合建模方法有可能取代當前雜糅了多個工具的分析流程。

隨著單細胞組學技術的進步，對組學分析流程開發的需求將會增長（Tanay＆Regev，2017）。未來的單細胞分析平臺將必須能夠處理不同的數據源，例如DNA甲基化（Smallwood et al., 2014），染色質可及性（Buenrostro et al., 2015）或蛋白質豐度數據（Stoeckius et al., 2017），以及能對這些數據進行整合分析。如果是這個目的，將不能只使用我們流程起點的單個reads或計數矩陣作為輸入。而且，整合RNA velocity分析的平臺已經在適應多組學（多模態）數據結構，輸入的是未剪接和剪接的reads數據（La Manno et al，2018）。單細胞多組學整合可以通過一致性聚類方法、多組學因子分析（Argelaguet et al., 2018）或多組學基因調控網絡推斷（Colomé-Tatché＆Theis，2018）來實現。具有這些功能的分析流程將是未來開發的方向。我們預計這種多組學分析流程將建立在我們為scRNA-seq分析奠定的基礎上。

文章解讀：苑曉梅
文章校對：生信寶典

推薦閱讀

• cellassign：用于腫瘤微環境分析的單細胞注釋工具（9月Nature）

• NC |SCALE準確鑒定單細胞ATAC-seq數據中染色質開放特征

• Cell子刊 | 植物單細胞轉錄組綜述·植物功能基因組學的高分辨率研究方法

• 七龍珠|召喚一份單細胞數據庫匯總

• RNA-seq最強綜述名詞解釋&思維導圖|關于RNA-seq，你想知道的都在這（續）

參考文獻：

Luecken MD, Theis FJ. Current best practices in single-cell RNA-seq analysis:a tutorial. Mol Syst Biol. 2019 Jun 19;15(6):e8746\. doi: 10.15252

總結

以上是生活随笔為你收集整理的重磅综述：三万字长文读懂单细胞RNA测序分析的最佳实践教程 （原理、代码和评述）的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。