原文鏈接：

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

主編評(píng)語

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

Abstract

單細(xì)胞RNA-seq使研究者能夠以前所未有的分辨率研究基因表達(dá)圖譜。這一潛力吸引著更多科研工作者應(yīng)用單細(xì)胞分析技術(shù)解決研究問題。隨著可用的分析工具越來越多，如何組合成一個(gè)最新最好的數(shù)據(jù)分析流程也越來越難。我們?cè)敿?xì)闡述了一個(gè)典型的單細(xì)胞轉(zhuǎn)錄組分析各個(gè)步驟的細(xì)節(jié)和注意事項(xiàng)，包括預(yù)處理（質(zhì)控、標(biāo)準(zhǔn)化、數(shù)據(jù)校正、特征選擇、降維）和細(xì)胞/基因水平的下游分析等。基于獨(dú)立的比較研究，我們?yōu)槊恳徊蕉纪扑]了當(dāng)前最好的方法和操作建議。隨后把這些最好的工具整合成一個(gè)分析流程并應(yīng)用于一套公共數(shù)據(jù)集的分析以演示其效果。案例研究具體可見https://www.github.com/theislab/single-cell-tutorial。這篇綜述為這個(gè)領(lǐng)域的新人提供了一份學(xué)習(xí)單細(xì)胞分析的指南，并且也能幫助老用戶更新他們的分析流程。

背景

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

分析標(biāo)準(zhǔn)化面臨的挑戰(zhàn)來源于越來越多的可用分析方法（截至2019年3月7日有385種工具）和數(shù)據(jù)集規(guī)模爆炸性的提高。使得我們一直在尋找新的方法來分析處理我們的數(shù)據(jù)。例如，最近已經(jīng)有方法可以預(yù)測(cè)細(xì)胞分化過程中的命運(yùn)選擇。盡管分析工具的不斷改進(jìn)有助于產(chǎn)生新的科學(xué)推論，但它也使分析流程的標(biāo)準(zhǔn)化變得更為復(fù)雜。

標(biāo)準(zhǔn)化的另一個(gè)挑戰(zhàn)在于軟件技術(shù)方面。用于scRNA-seq數(shù)據(jù)的分析工具是用不同的編程語言編寫的 - 最主要的是R和Python。盡管跨編程語言的支持越來越多，但使用的編程語言確實(shí)影響了對(duì)分析工具的選擇。諸如Seurat，Scater或Scanpy等常用工具提供了集成環(huán)境來開發(fā)流程并包含大量分析工具。然而，出于維護(hù)的需要， 這些平臺(tái)限制了它們只能使用各自的編程語言開發(fā)的工具。引申開來，當(dāng)前可用的scRNA-seq分析教程也存在語言限制，而且許多是圍繞這3大平臺(tái)進(jìn)行講述的 （易生信2020年最新的單細(xì)胞課程會(huì)同時(shí)提供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進(jìn)行單細(xì)胞分析？

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

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

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

圖1. 典型單細(xì)胞轉(zhuǎn)錄組分析流程。

原始測(cè)序數(shù)據(jù)預(yù)處理和比對(duì)后獲得Gene count矩陣，作為本分析流程的開始。這些結(jié)果圖是Count矩陣采用最佳實(shí)戰(zhàn)分析流程進(jìn)行預(yù)處理和下游分析獲得的。數(shù)據(jù)來源于Haber et al 2017腸上皮細(xì)胞文章。

預(yù)處理和可視化

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

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

Box1：scRNA-seq實(shí)驗(yàn)流程的關(guān)鍵要素

從生物樣品到單細(xì)胞數(shù)據(jù)需要多步實(shí)驗(yàn)操作。典型的工作流程包含單細(xì)胞解離、單細(xì)胞分離、文庫構(gòu)建和測(cè)序。我們?cè)谶@里簡要介紹這些步驟。不同protocols的更詳細(xì)的解釋和比較可參考下面三篇文章：

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

單細(xì)胞實(shí)驗(yàn)的輸入材料通常是生物組織樣品。第一步，單細(xì)胞解離：消化組織產(chǎn)生單細(xì)胞懸液。為了分別分析每個(gè)細(xì)胞中的mRNA，必須分離單細(xì)胞。根據(jù)實(shí)驗(yàn)方案不同，單細(xì)胞分離的方式也有所不同。

• 基于平板的技術(shù)將細(xì)胞分到到板上的孔中。

• 基于液滴的方法則依賴于微流體液滴捕獲單個(gè)細(xì)胞。

在這兩種情況下，都可能出現(xiàn)一些問題，如多個(gè)細(xì)胞一起被捕獲（doublets or multiplets）、非活細(xì)胞被捕獲或根本沒有細(xì)胞被捕獲（空液滴/孔）。基于液滴的方法需要通過低的輸入細(xì)胞濃度來保持低的doublets率，因此空液滴是特別常見的 （生信寶典注*：一般beads和細(xì)胞的輸入比例是20:1*）。

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

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

質(zhì)控

在分析單細(xì)胞基因表達(dá)數(shù)據(jù)之前，我們必須確保所有barcode都對(duì)應(yīng)于有效細(xì)胞 (viable cells，有活力的細(xì)胞)。質(zhì)控有3個(gè)指標(biāo)：

• 測(cè)到的轉(zhuǎn)錄本分子總數(shù)

• 測(cè)到的基因總數(shù)

• 來源于線粒體基因的轉(zhuǎn)錄本所占比例

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

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

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

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

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

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

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

陷阱和建議：

• 通過可視化檢測(cè)到的基因數(shù)量、總分子數(shù)和線粒體基因的表達(dá)比例的分布中的異常峰來執(zhí)行QC。聯(lián)合考慮這3個(gè)變量，而不是單獨(dú)考慮一個(gè)變量。

• 盡可能設(shè)置寬松的QC閾值；如果下游聚類無法解釋時(shí)再重新設(shè)定嚴(yán)格的QC閾值。

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

標(biāo)準(zhǔn)化

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

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

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

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

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

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

標(biāo)準(zhǔn)化是對(duì)細(xì)胞計(jì)數(shù)數(shù)據(jù)進(jìn)行縮放處理以使其在細(xì)胞之間可比，也可以在基因?qū)用鎸?duì)基因計(jì)數(shù)進(jìn)行歸一化 (scale)以便于基因內(nèi)部進(jìn)行直接比較。基因歸一化是指一個(gè)基因減去其在所有樣品表達(dá)的均值然后除以其在所有樣品表達(dá)值的標(biāo)準(zhǔn)差。歸一化后，這個(gè)基因在所有樣品表達(dá)值均值為0，用單位方差形式表示其表達(dá)值。歸一化后，所有基因在下游分析時(shí)權(quán)重是一樣的。（生信寶典注*：歸一化后，原始基因的表達(dá)豐度信息就沒了，換成了無標(biāo)度的標(biāo)準(zhǔn)差的表示。）是否對(duì)基因進(jìn)行歸一化目前尚無達(dá)成共識(shí)。盡管流行Seurat教程通常應(yīng)用gene scaling，但Slingshot方法的作者在其教程中選擇了不對(duì)基因進(jìn)行scaling。兩種選擇的爭議點(diǎn)在：所有基因不論表達(dá)高低在進(jìn)行下游分析時(shí)權(quán)重一致 ?(生信寶典注*：**scale后所有基因權(quán)重一致)，還是基因表達(dá)量的絕對(duì)值對(duì)下游分析也有貢獻(xiàn) (生信寶典注*：突出高表達(dá)基因?qū)ο掠蔚母笥绊?) 。為了從數(shù)據(jù)中保留盡可能多的生物相關(guān)信息，在本教程中，我們選擇不進(jìn)行歸一化 （scaling genes）這一步操作。

標(biāo)準(zhǔn)化后，數(shù)據(jù)矩陣通常進(jìn)行l(wèi)og(x+1)轉(zhuǎn)換。此轉(zhuǎn)換具有三個(gè)重要作用。

• 首先，對(duì)數(shù)轉(zhuǎn)換后的表達(dá)式值之間的差值可對(duì)應(yīng)于對(duì)數(shù)轉(zhuǎn)換后的倍數(shù)變化，這是衡量基因表達(dá)變化的常用方法。

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

• 最后，對(duì)數(shù)轉(zhuǎn)換可以減少數(shù)據(jù)的偏態(tài)分布，從而使數(shù)據(jù)近似于正態(tài)分布，更符合許多下游分析工具對(duì)數(shù)據(jù)分布的假設(shè)要求。

盡管scRNA-seq數(shù)據(jù)實(shí)際上不是對(duì)數(shù)正態(tài)分布的，但這三個(gè)效果使對(duì)數(shù)轉(zhuǎn)換成為一種粗略但有用的工具。這一有用性在下游差異表達(dá)分析或批次校正時(shí)有更好的體現(xiàn)。但是應(yīng)該注意的是，數(shù)據(jù)的對(duì)數(shù)轉(zhuǎn)換會(huì)在數(shù)據(jù)中引入虛假的差異表達(dá)結(jié)果。而且如果size factor在組間差異更大時(shí)影響尤其明顯。

陷阱和建議：

• 我們建議使用scran來標(biāo)準(zhǔn)化非全長數(shù)據(jù)集。另一種選擇是通過scone評(píng)估標(biāo)準(zhǔn)化方法，尤其是對(duì)基于板的數(shù)據(jù)集。全長scRNA-seq數(shù)據(jù)可以使用bulk轉(zhuǎn)錄組的方法校正基因長度的影響。

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

• 標(biāo)準(zhǔn)化后的數(shù)據(jù)應(yīng)進(jìn)行l(wèi)og(x+1)轉(zhuǎn)換，以使數(shù)據(jù)更符合正態(tài)分布，滿足下游分析方法的初始假設(shè)。

數(shù)據(jù)校正和整合

如前所述，數(shù)據(jù)標(biāo)準(zhǔn)化是去除實(shí)驗(yàn)過程中隨機(jī)性的影響 (count sampling)。但是，標(biāo)準(zhǔn)化后的數(shù)據(jù)可能仍然包含有與研究不相干的因素帶來的影響。數(shù)據(jù)校正的目的就是進(jìn)一步去除技術(shù)因素和非關(guān)注的生物學(xué)混雜因素，例如批次、dropout或細(xì)胞周期不同帶來的影響 (如何火眼金睛鑒定那些單細(xì)胞轉(zhuǎn)錄組中的混雜因素)。需要注意的是這些混淆因素并不總是需要校正。相反，考慮在分析中移除哪些因素取決于下游分析目的。我們建議分別考慮對(duì)生物學(xué)和技術(shù)混雜因素 (covariates)的校正，因?yàn)樗鼈冇糜诓煌哪康牟⒋聿煌姆治鎏魬?zhàn)。

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

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

在校正特定的生物因素影響之前，應(yīng)考慮幾個(gè)方面。首先，校正生物學(xué)因素影響并不總是有助于解釋scRNA-seq數(shù)據(jù)。雖然消除細(xì)胞周期影響可以改善對(duì)發(fā)育軌跡的推斷，但細(xì)胞周期信號(hào)也可以提供有意義的生物學(xué)信息。例如，可以根據(jù)細(xì)胞周期評(píng)分確定增殖細(xì)胞群 (在Github的案例研究中有講)。同樣，必須根據(jù)具體研究的對(duì)象理解生物學(xué)信號(hào)。假設(shè)生物過程發(fā)生在同一有機(jī)體內(nèi)，則這些過程之間可能存在依賴性。因此，校正一個(gè)過程的影響可能會(huì)無意間掩蓋另一個(gè)過程的信號(hào)。另外，也有研究者提出細(xì)胞大小對(duì)轉(zhuǎn)錄組的影響通常也歸因于細(xì)胞周期不同。因此，通過標(biāo)準(zhǔn)化或?qū)Ｓ霉ぞ?#xff08;例如cgCorrect）在校正細(xì)胞大小的同時(shí)也部分校正了scRNA-seq數(shù)據(jù)中細(xì)胞周期的影響。

去除技術(shù)影響

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

消除測(cè)序深度影響的另一個(gè)策略是使用更嚴(yán)格的標(biāo)準(zhǔn)化程序，例如downsampling或非線性歸一化方法 (見前面標(biāo)準(zhǔn)化部分)。這些方法可能更適用于基于板的scRNA-seq數(shù)據(jù)集，因?yàn)榧?xì)胞之間較大的測(cè)序深度差異可能掩蓋細(xì)胞之間的異質(zhì)性。

批次效應(yīng)和數(shù)據(jù)整合

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

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

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

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

缺失值填充

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

陷阱和建議：

• 僅在進(jìn)行軌跡推斷和校正的生物學(xué)混雜因素不影響感興趣的生物學(xué)過程時(shí)才校正這些因素的影響。

• 如果校正的話，所有因素同時(shí)校正而不是分別校正技術(shù)和非關(guān)注的生物因素變量。

• 基于板的數(shù)據(jù)集預(yù)處理時(shí)需要校正count數(shù)的影響，建議采用非線性標(biāo)準(zhǔn)化方法或downsampling方法進(jìn)行標(biāo)準(zhǔn)化。

• 當(dāng)批次之間的細(xì)胞類型和狀態(tài)組成一致時(shí)，建議通過ComBat執(zhí)行批次校正。

• 數(shù)據(jù)整合和批次校正應(yīng)通過不同的方法進(jìn)行。數(shù)據(jù)整合工具可能會(huì)過度校正簡單的批次效應(yīng)。

• 用戶需要對(duì)只在缺失值填充后才能發(fā)現(xiàn)的信號(hào)格外注意。探索性分析時(shí)最好不進(jìn)行缺失值填充操作。

特征選擇、降維和可視化

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

特征選擇

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

圖EV1. 不同大小數(shù)據(jù)集選擇HVG基因數(shù)量的分布展現(xiàn)。**（數(shù)據(jù)來源于手動(dòng)整理的最近發(fā)表的多篇scRNA-seq文章，顏色代表不同的技術(shù)方式）

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

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

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

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

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

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

可視化

可視化時(shí)一般使用非線性降維方法（圖4）。scRNA-seq數(shù)據(jù)可視化的最常見的降維方法是t‐SNE?（?t‐distributed stochastic neighbour embedding） （還在用PCA降維？快學(xué)學(xué)大牛最愛的t-SNE算法吧, 附Python/R代碼）。t‐SNE的維度著重于以犧牲全局結(jié)構(gòu)為代價(jià)來保留局部相似性 (生信寶典注*：PCA則是盡可能多的保留全局差異*)。因此，這些可視化可能會(huì)夸大細(xì)胞群體之間的差異，而忽略群體之間的潛在聯(lián)系 (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)。另一個(gè)困難是對(duì)參數(shù)perplexity parameter的選擇，因?yàn)閠-SNE圖會(huì)因?yàn)檫@個(gè)參數(shù)值不同而顯示出明顯不同的分簇?cái)?shù)。t‐SNE的常見替代方法是UMAP?(Uniform Approximation and Projection method)或基于圖的工具，例如SPRING。UMAP和SPRING的力導(dǎo)向布局算法ForceAtlas2(*force‐directed layout algorithm*)可以說是數(shù)據(jù)潛在拓?fù)浣Y(jié)構(gòu)展示最好的方法（Wolf et al，2019）。在此比較中，使UMAP與眾不同的是它的速度快和能應(yīng)用于更大規(guī)模數(shù)據(jù)的能力（Becht et al，2018）。因此，在沒有特定生物學(xué)問題限制的情況下，我們將UMAP視為探索性數(shù)據(jù)可視化分析的最佳實(shí)踐。此外，UMAP還可以把數(shù)據(jù)降維到二維以上的新數(shù)據(jù)。盡管我們不知道UMAP`在數(shù)據(jù)匯總中的任何應(yīng)用，但它可能是PCA的一個(gè)合適的替代方案。

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

陷阱和建議：

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

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

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

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

• UMAP+PAGA是可視化特別復(fù)雜的數(shù)據(jù)集的合適替代方法。

數(shù)據(jù)預(yù)處理步驟

雖然我們?cè)谏厦姘错樞蚋攀隽藄cRNA‐seq中常見的預(yù)處理步驟，但下游分析通常需要不同級(jí)別的預(yù)處理數(shù)據(jù)，建議根據(jù)下游應(yīng)用調(diào)整預(yù)處理的各個(gè)步驟。為了向新用戶說明這種情況，我們將預(yù)處理分為五個(gè)數(shù)據(jù)處理階段：（i）原始數(shù)據(jù)，（ii）標(biāo)準(zhǔn)化后的數(shù)據(jù)，（iii）校正后的數(shù)據(jù)，（iv）特征選擇后的數(shù)據(jù)，以及（v）降維后的數(shù)據(jù)。數(shù)據(jù)處理的這些階段可以歸類為三個(gè)預(yù)處理層：測(cè)量的數(shù)據(jù)，校正的數(shù)據(jù)和降維的數(shù)據(jù)。細(xì)胞和基因的質(zhì)量控制篩選是必須執(zhí)行的步驟，因此在此處略過。盡管層的順序代表了scRNA-seq分析的典型工作流程，但也可以跳過某一步或在處理階段的順序中稍作更改。例如，單批次數(shù)據(jù)集可能不需要批次校正。在表1中，我們總結(jié)了預(yù)處理數(shù)據(jù)的每一層所適用的下游處理程序。

表1

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

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

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

基因表達(dá)的統(tǒng)計(jì)比較需要基于測(cè)量數(shù)據(jù)層。去噪、校正批次或其他噪聲源的方法總不是太完美。因此，數(shù)據(jù)校正方法不可避免地會(huì)對(duì)數(shù)據(jù)進(jìn)行過度或不足校正，因此會(huì)以意想不到的方式改變至少某些基因表達(dá)譜的方差。基因表達(dá)的統(tǒng)計(jì)檢驗(yàn)依賴于將背景方差評(píng)估為數(shù)據(jù)噪聲的空模型。隨著數(shù)據(jù)校正趨于減少背景變化（圖EV2），背景變化被數(shù)據(jù)校正方法過度校正的基因?qū)⒏锌赡鼙辉u(píng)估為顯著差異表達(dá)基因。此外，某些數(shù)據(jù)校正方法（例如ComBat）將不吻合于實(shí)驗(yàn)設(shè)計(jì)的表達(dá)信號(hào)定義為噪聲，隨后將其從數(shù)據(jù)中刪除。這一基于實(shí)驗(yàn)設(shè)計(jì)的優(yōu)化方法除了可能造成數(shù)據(jù)中的噪聲被低估還會(huì)導(dǎo)致對(duì)效應(yīng)大小 （effect size）的高估。鑒于這些考慮，與校正后的數(shù)據(jù)相比，基于測(cè)量數(shù)據(jù)進(jìn)行差異檢驗(yàn)是一種更為保守的方法。使用測(cè)量的數(shù)據(jù)時(shí)，可以并且應(yīng)該在差異統(tǒng)計(jì)分析模型中考慮技術(shù)協(xié)變量 （參考DESeq2差異基因分析和批次效應(yīng)移除）。

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

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

陷阱和建議：

• 原始測(cè)量的數(shù)據(jù)用于差異基因統(tǒng)計(jì)檢驗(yàn)分析；

• 校正后的數(shù)據(jù)用于數(shù)據(jù)的可視化比較；

• 降維后的數(shù)據(jù)用于其他基于數(shù)據(jù)形態(tài)（biological data manifold）的下游分析。

下游分析

經(jīng)過預(yù)處理后，我們稱之為下游分析的方法指應(yīng)用于生物學(xué)發(fā)現(xiàn)并描述潛在的生物學(xué)系統(tǒng)的方法。通過將可以解釋的模型擬合到數(shù)據(jù)中獲得相應(yīng)的結(jié)論，比如有相似基因表達(dá)譜的細(xì)胞群代表一個(gè)細(xì)胞簇、相似細(xì)胞之間基因表達(dá)的微小變化指示連續(xù)（分化）軌跡；或具有相關(guān)表達(dá)趨勢(shì)的基因指示共調(diào)控等。

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

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

在此，我們?cè)谠敿?xì)描述獨(dú)立于這些細(xì)胞結(jié)構(gòu)進(jìn)行的基因水平分析之前，先描述細(xì)胞水平和基因水平的聚類和軌跡分析工具。

聚類分析

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

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

社群檢測(cè)方法是一種圖分類算法，依賴于單細(xì)胞數(shù)據(jù)的網(wǎng)絡(luò)圖表示。圖是使用K最近鄰方法（KNN圖）獲得的。細(xì)胞在圖中表示為節(jié)點(diǎn)。每個(gè)細(xì)胞與其K個(gè)最相似的細(xì)胞相連，通常對(duì)主成分降維后的數(shù)據(jù)計(jì)算歐氏距離作為相似性度量。根據(jù)數(shù)據(jù)集的大小，K通常設(shè)置為5到100個(gè)最近的鄰居。獲得的KNN結(jié)果圖捕獲了表達(dá)圖譜的基礎(chǔ)拓?fù)浣Y(jié)構(gòu)（Wolf et al., 2019）。表達(dá)譜相似的細(xì)胞集合對(duì)應(yīng)于圖的緊密連接區(qū)域，然后使用社群檢測(cè)方法檢測(cè)圖中這些緊密連接區(qū)域。社群檢測(cè)方法通常比聚類方法要快，因?yàn)橹挥邢噜彽募?xì)胞對(duì)會(huì)被視為屬于同一群集。因此，這種方法大大減少了可能的細(xì)胞簇的搜索空間。

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

陷阱和建議：

• 我們建議在單細(xì)胞KNN圖上執(zhí)行Louvain社群檢測(cè)算法獲得聚類簇。

• 細(xì)胞聚簇不一定使用單一分辨率進(jìn)行。對(duì)特定細(xì)胞簇進(jìn)行細(xì)化聚類是鑒定數(shù)據(jù)集中亞結(jié)構(gòu)的有效方法，可以發(fā)現(xiàn)細(xì)胞亞群。

細(xì)胞簇注釋

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

盡管把單細(xì)胞數(shù)據(jù)中鑒定到的簇歸類為一個(gè)有生物意義的細(xì)胞類型是一個(gè)誘人的結(jié)果，但是細(xì)胞身份的確定不是那么簡單（Wagner et al., 2016；Clevers et al., 2017）。首先，并不總是清楚細(xì)胞類型定義的精確程度。例如，“T細(xì)胞”可能是某些人滿意的細(xì)胞類型標(biāo)記，但其他人可能會(huì)在數(shù)據(jù)集中尋找T細(xì)胞亞型并區(qū)分CD4+T細(xì)胞和CD8+T細(xì)胞?。另外，同一細(xì)胞類型的不同狀態(tài)可能會(huì)被分到不同的細(xì)胞簇。基于上述原因，最好使用術(shù)語cell identities而不是cell types。在對(duì)細(xì)胞進(jìn)行聚類和注釋之前，用戶必須確定要注釋到的細(xì)胞信息精確水平，從而確定合適的聚簇分辨率。

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

圖6. 聚類分析結(jié)果

(A) UMAP可視化Louvain算法鑒定的細(xì)胞簇，并標(biāo)記注釋的細(xì)胞類型名稱。(B) 通過細(xì)胞類型Marker基因的表達(dá)標(biāo)記細(xì)胞簇所代表的細(xì)胞類型，如干細(xì)胞 (Slc12a2), 腸上皮細(xì)胞 (Arg2), 杯狀細(xì)胞 (Tff3) 和潘氏細(xì)胞 (Defa24)。灰色代表低表達(dá)，紅色代表高表達(dá)。Marker基因也可能在其它細(xì)胞簇有表達(dá)，比如Tff3在杯狀細(xì)胞和潘氏細(xì)胞都有表達(dá)，只是在杯狀細(xì)胞更高。? 腸上皮細(xì)胞區(qū)域細(xì)胞類型組成熱圖（上面是近端區(qū)域，下面是遠(yuǎn)端區(qū)域）。相對(duì)細(xì)胞密度高的區(qū)域用深紅色表示。

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

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

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

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

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

陷阱和建議：

• 不要使用標(biāo)記基因的P-value值來驗(yàn)證細(xì)胞簇身份，尤其是當(dāng)檢測(cè)到的標(biāo)記基因無益于注釋細(xì)胞簇時(shí)。P值可能會(huì)被夸大。

• 請(qǐng)注意，由于數(shù)據(jù)集的細(xì)胞類型和狀態(tài)組成，同一細(xì)胞簇的標(biāo)記基因可能在數(shù)據(jù)集之間有所不同。

• 如果存在相關(guān)的參考數(shù)據(jù)集，我們建議綜合使用自動(dòng)聚類注釋和基于數(shù)據(jù)的標(biāo)記基因的手動(dòng)注釋來注釋細(xì)胞類型。

細(xì)胞組成分析?（Compositional analysis）。在細(xì)胞層面，我們可以根據(jù)其組成結(jié)構(gòu)分析聚類數(shù)據(jù)。細(xì)胞組成數(shù)據(jù)分析圍繞不同樣品落入每個(gè)細(xì)胞簇的細(xì)胞比例進(jìn)行分析。這一比例可能在疾病狀態(tài)下發(fā)生變化。例如，研究表明沙門氏菌感染會(huì)增加小鼠腸上皮區(qū)域腸上皮細(xì)胞（enterocytes）的比例（Haber et al.，2017）。

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

陷阱和建議：

• 統(tǒng)計(jì)檢驗(yàn)分析時(shí)需要考慮到樣本之間的細(xì)胞簇比例的變化是互相依賴的（非獨(dú)立的）

Trajectory analysis

**軌跡推斷。**細(xì)胞多樣性不能通過離散的分類系統(tǒng)（例如細(xì)胞聚類）充分描述。驅(qū)動(dòng)觀察到的細(xì)胞異質(zhì)性發(fā)展的生物進(jìn)程是一個(gè)連續(xù)過程（Tanay＆Regev，2017）。因此，為了捕獲細(xì)胞身份之間的過渡狀態(tài)、不同的分化分支或生物學(xué)功能的漸進(jìn)式非同步變化，我們需要?jiǎng)討B(tài)的基因表達(dá)模型。這類方法稱為軌跡推斷（trajectory inference，TI）。

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

圖7. 小鼠腸上皮細(xì)胞軌跡分析結(jié)果展示

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

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

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

推斷的軌跡不一定要完全對(duì)應(yīng)生物發(fā)育過程。首先，推斷的軌跡僅表示轉(zhuǎn)錄相似性。很少有TI方法在其模型中包括不確定性評(píng)估（Griffiths et al., 2018）。因此，需要更多的信息來驗(yàn)證是否確實(shí)捕獲了生物過程。這些信息可以來源于干擾實(shí)驗(yàn)、推斷的調(diào)控基因動(dòng)力學(xué)以及RNA velocity數(shù)據(jù)的支持等。

陷阱和建議：

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

• 推斷的軌跡不需要完全對(duì)應(yīng)生物過程， 應(yīng)該收集更多的證據(jù)來解釋軌跡。

**基因表達(dá)動(dòng)力學(xué)。**推斷的軌跡也可能是轉(zhuǎn)錄噪聲擬合的結(jié)果，一個(gè)排除方式是在基因水平上對(duì)鑒定出的軌跡進(jìn)行進(jìn)一步分析。在整個(gè)偽時(shí)間范圍內(nèi)連續(xù)平穩(wěn)變化的基因是決定軌跡的關(guān)鍵基因，可用于識(shí)別潛在的生物學(xué)過程。此外，這一組與軌跡相關(guān)的基因應(yīng)該包含調(diào)控相應(yīng)生物進(jìn)程的基因。調(diào)節(jié)基因可以幫助我們理解該生物過程是如何和為什么被觸發(fā)的，并且可能是潛在的藥物靶標(biāo)（Gashaw et al., 2012）。

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

當(dāng)前幾乎沒有專用于分析基因表達(dá)時(shí)間動(dòng)力學(xué)的分析工具。BEAM是整合于Monocle?軌跡推斷方法中的一個(gè)工具，可以分析分支特異的基因表達(dá)動(dòng)力學(xué)。除了這個(gè)，用戶可以選擇使用LineagePulse?(https://github.com/YosefLab/LineagePulse)，這個(gè)工具考慮了dropout噪音，但仍然在開發(fā)中。或者基于標(biāo)準(zhǔn)R包或Limma包寫作自己的分析測(cè)試工具。在Slingshot`的使用幫助中可以找到一個(gè)這樣應(yīng)用的例子供參考。

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

**亞穩(wěn)態(tài)。**軌跡的細(xì)胞水平分析是指研究偽時(shí)間軸上的細(xì)胞密度。假設(shè)以無偏方式對(duì)細(xì)胞進(jìn)行隨機(jī)采樣，軌跡中細(xì)胞密集的區(qū)域表示優(yōu)選的轉(zhuǎn)錄狀態(tài)。當(dāng)將軌跡解釋為時(shí)間過程時(shí)，這些密集區(qū)域可能代表例如發(fā)育中的亞穩(wěn)態(tài)（Haghverdi et al., 2016）。我們可以通過繪制偽時(shí)間坐標(biāo)的直方圖來找到這些亞穩(wěn)態(tài)（圖7C）。

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

**基因水平分析。**到目前為止，雖然我們專注于表征細(xì)胞結(jié)構(gòu)的基因水平分析方法，但是單細(xì)胞數(shù)據(jù)的基因水平分析有更廣泛的內(nèi)容。差異表達(dá)分析、基因集分析和基因調(diào)控網(wǎng)絡(luò)分析都可以直接研究數(shù)據(jù)中的分子信號(hào)。這些方法不是描述細(xì)胞異質(zhì)性，而是把這種異質(zhì)性作為理解基因表達(dá)差異的背景。

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

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

因?yàn)橛糜诓町惢蚍治龅氖菢?biāo)準(zhǔn)化后但未進(jìn)行技術(shù)變量和非相關(guān)生物變量校正的數(shù)據(jù)，因此在分析過程中考慮這些混雜變量對(duì)獲得穩(wěn)定的結(jié)果是關(guān)鍵的。通常差異基因分析工具都允許用戶考慮混雜變量（confounders），但是把哪些變量納入模型需要謹(jǐn)慎評(píng)估。比如，在多數(shù)單細(xì)胞數(shù)據(jù)集中，樣本和實(shí)驗(yàn)處理?xiàng)l件是互斥的，因?yàn)閹缀醪豢赡軐?duì)一個(gè)樣品進(jìn)行多個(gè)實(shí)驗(yàn)處理。如果我們?cè)谀Ｐ椭型瑫r(shí)考慮樣品和實(shí)驗(yàn)條件兩個(gè)協(xié)變量，跟這些協(xié)變量相關(guān)的變化則不會(huì)被明確確定。因此在考慮不同實(shí)驗(yàn)條件下的差異分析時(shí)，通常不在模型中引入樣品信息作為協(xié)變量。當(dāng)校正多個(gè)確定的批次分類變量時(shí)，可視化展示協(xié)變量的混雜影響將會(huì)變得更困難。在這個(gè)情況下，需要先確認(rèn)模型設(shè)計(jì)矩陣是否滿秩 （the model design matrix is full rank，full rank是指模型設(shè)計(jì)矩陣中的每一列都不能通過對(duì)其它列的線性組合得來，即它們之間不存在線性相關(guān)）。即便一些設(shè)計(jì)矩陣不是滿秩，差異基因分析工具通常也不會(huì)給出警告而是繼續(xù)運(yùn)行。這時(shí)獲得的結(jié)果將可能不是預(yù)期的分析方向。

我們這兒描述的場(chǎng)景中，實(shí)驗(yàn)條件協(xié)變量是在實(shí)驗(yàn)設(shè)計(jì)中決定的。因此在同一簇內(nèi)基于這一協(xié)變量的差異基因分析是獨(dú)立于聚類過程的。這將基于實(shí)驗(yàn)條件的差異基因分析和基于細(xì)胞簇的差異基因分析區(qū)別開來。基于實(shí)驗(yàn)條件的差異基因檢測(cè)獲得的P-values是顯著性檢驗(yàn)統(tǒng)計(jì)指標(biāo)，需要進(jìn)一步做多重檢驗(yàn)校正。為了降低多重檢驗(yàn)校正的壓力，不感興趣的基因通常需要預(yù)先排除掉。雖然假基因和非編碼基因也可以提供有用信息，但在分析中也通常會(huì)被忽略掉。

陷阱和建議：

• DE測(cè)試不應(yīng)基于校正后的數(shù)據(jù)（去噪，批次校正等），而應(yīng)基于原始測(cè)量數(shù)據(jù)并在模型中引入干擾協(xié)變量。

• 用戶不應(yīng)依賴差異基因檢測(cè)工具校正帶有混雜協(xié)變量的模型。設(shè)計(jì)矩陣模型應(yīng)當(dāng)保證滿秩。

• 我們建議使用MAST或limma進(jìn)行差異基因分析。

• 39個(gè)轉(zhuǎn)錄組分析工具，120種組合評(píng)估(轉(zhuǎn)錄組分析工具哪家強(qiáng)-導(dǎo)讀版)

• DESeq2差異基因分析和批次效應(yīng)移除

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

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

不同用途的基因集信息可以在校正過的基因注釋數(shù)據(jù)庫中獲取。為了解釋差異基因結(jié)果，通常按共有的生物進(jìn)程對(duì)基因進(jìn)行分組。生物進(jìn)程信息存儲(chǔ)于MSigDB、Gene Ontology，通路數(shù)據(jù)存儲(chǔ)于KEGG和Reactome數(shù)據(jù)庫。基因列表的功能注釋富集分析有多種工具可用，具體見Huang et al (2009)和?Tarca et al (2013)的綜述。

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

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

• 沒錢買KEGG怎么辦？REACTOME開源通路更強(qiáng)大

• 一文掌握GSEA，超詳細(xì)教程

• 這個(gè)只需一步就可做富集分析的網(wǎng)站還未發(fā)表就被CNS等引用超過350次

• 單基因GSEA怎么做？

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

**基因調(diào)控網(wǎng)絡(luò)。**基因不是獨(dú)立發(fā)揮作用的。相反，基因的表達(dá)水平是由與其他基因和小分子之間的復(fù)雜調(diào)控決定的。揭示這些調(diào)控作用是基因調(diào)控網(wǎng)絡(luò)（GRN）推斷方法的目標(biāo)（SCENIC | 從單細(xì)胞數(shù)據(jù)推斷基因調(diào)控網(wǎng)絡(luò)和細(xì)胞類型）。

基因調(diào)控網(wǎng)絡(luò)推斷是基于對(duì)基因共表達(dá)如相關(guān)性、互信息（mutual information）或回歸模型分析進(jìn)行測(cè)量（Chen＆Mar，2018）。如果在考慮了其他所有基因都作為混雜因素后，兩個(gè)目標(biāo)基因之間依然存在共表達(dá)關(guān)系則可以推斷這兩個(gè)基因之間存在因果調(diào)控關(guān)系。基因調(diào)控關(guān)系鑒定與軌跡相關(guān)的調(diào)控基因檢測(cè)相關(guān)。確實(shí)，幾種單細(xì)胞GRN推論方法使用的機(jī)械微分方程模型中都引入了軌跡信息（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）。

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

陷阱和建議：

• 用戶應(yīng)警惕推斷的調(diào)控關(guān)系中的不確定性。基因調(diào)控關(guān)系豐富的基因模塊比調(diào)控信息稀疏的模塊更可靠。

**分析平臺(tái)。**單細(xì)胞分析流程是多個(gè)獨(dú)立開發(fā)的工具的集合。為了促進(jìn)這些工具之間的數(shù)據(jù)傳遞，單細(xì)胞工具開發(fā)時(shí)都考慮使用一致的數(shù)據(jù)格式。這些工具為構(gòu)建分析流程提供了基礎(chǔ)。目前可用的分析平臺(tái)有基于R（McCarthy et al., 2017; Butler et al., 2018）或Python（Wolf et al., 2018）的命令行，或具有圖形用戶界面（GUI）的本地應(yīng)用程序（Patel，2018; Scholz et al., 2018）或Web服務(wù)（Gardeux et al., 2017; ?Zhu et al., 2017）。Zhu et al.（2017）和Zappia et al.（2018）對(duì)這些分析平臺(tái)進(jìn)行了綜述。

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

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

結(jié)論與展望

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

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

隨著單細(xì)胞組學(xué)技術(shù)的進(jìn)步，對(duì)組學(xué)分析流程開發(fā)的需求將會(huì)增長（Tanay＆Regev，2017）。未來的單細(xì)胞分析平臺(tái)將必須能夠處理不同的數(shù)據(jù)源，例如DNA甲基化（Smallwood et al., 2014），染色質(zhì)可及性（Buenrostro et al., 2015）或蛋白質(zhì)豐度數(shù)據(jù)（Stoeckius et al., 2017），以及能對(duì)這些數(shù)據(jù)進(jìn)行整合分析。如果是這個(gè)目的，將不能只使用我們流程起點(diǎn)的單個(gè)reads或計(jì)數(shù)矩陣作為輸入。而且，整合RNA velocity分析的平臺(tái)已經(jīng)在適應(yīng)多組學(xué)（多模態(tài)）數(shù)據(jù)結(jié)構(gòu)，輸入的是未剪接和剪接的reads數(shù)據(jù)（La Manno et al，2018）。單細(xì)胞多組學(xué)整合可以通過一致性聚類方法、多組學(xué)因子分析（Argelaguet et al., 2018）或多組學(xué)基因調(diào)控網(wǎng)絡(luò)推斷（Colomé-Tatché＆Theis，2018）來實(shí)現(xiàn)。具有這些功能的分析流程將是未來開發(fā)的方向。我們預(yù)計(jì)這種多組學(xué)分析流程將建立在我們?yōu)閟cRNA-seq分析奠定的基礎(chǔ)上。

文章解讀：苑曉梅
文章校對(duì)：生信寶典

推薦閱讀

• cellassign：用于腫瘤微環(huán)境分析的單細(xì)胞注釋工具（9月Nature）

• NC |SCALE準(zhǔn)確鑒定單細(xì)胞ATAC-seq數(shù)據(jù)中染色質(zhì)開放特征

• Cell子刊 | 植物單細(xì)胞轉(zhuǎn)錄組綜述·植物功能基因組學(xué)的高分辨率研究方法

• 七龍珠|召喚一份單細(xì)胞數(shù)據(jù)庫匯總

• RNA-seq最強(qiáng)綜述名詞解釋&思維導(dǎo)圖|關(guān)于RNA-seq，你想知道的都在這（續(xù)）

參考文獻(xiàn)：

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

總結(jié)

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

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔推薦給好友。