概述
本章概述了典型的scRNA-seq分析工作流程的框架（圖1）。 隨后將更詳細(xì)地描述每個(gè)分析步驟。
實(shí)驗(yàn)設(shè)計(jì)
在開始分析本身之前蛉拙，對實(shí)驗(yàn)設(shè)計(jì)進(jìn)行一些討論可能會有所幫助。最明顯的問題是技術(shù)的選擇彻亲，大致可以分為：
1.基于液滴的：10XGenomics孕锄，inDrop, Dropseq
2.基于板的獨(dú)特分子識別符（UMI）：CEL-seq吮廉，MARS-seq
3.基于板的讀取：Smart-seq2
4.其他：sci-RNA-seq畸肆，Seq-Well
這些方法中的每一種都有其優(yōu)點(diǎn)和缺點(diǎn)宦芦，其他地方對此進(jìn)行了廣泛討論（Mereu等人(http://bioconductor.org/books/release/OSCA/overview.html#ref-mereu2019benchmarking); Ziegenhain等人2017(http://bioconductor.org/books/release/OSCA/overview.html#ref-ziegenhain2017comparative)）。實(shí)際上轴脐，基于液滴的技術(shù)由于其吞吐量和每個(gè)細(xì)胞的低成本而成為當(dāng)前的事實(shí)上的標(biāo)準(zhǔn)调卑。基于板的方法可以捕獲其他表型信息（例如形態(tài)）大咱，并且更適合定制恬涧。基于reads的方法提供了完整的轉(zhuǎn)錄本覆蓋范圍碴巾，這在某些應(yīng)用（例如剪接溯捆，外顯子組突變）中很有用; 基于UMI的方法將減輕PCR擴(kuò)增噪聲的影響，因此更為流行厦瓢。方法的選擇取決于具體情況——但是以下我們分析流程中的大多數(shù)方面都與技術(shù)無關(guān)提揍。
下一個(gè)問題是應(yīng)該捕獲多少個(gè)細(xì)胞，以及應(yīng)該對它們進(jìn)行測序的深度煮仇。簡短的答案是“盡可能多地花錢”劳跃。長答案是，這取決于分析的目的浙垫。如果我們旨在發(fā)現(xiàn)稀有的細(xì)胞亞群刨仑，那么我們需要更多的細(xì)胞。如果我們旨在刻畫細(xì)微的差異夹姥，那么我們需要更多的測序深度贸人。目前，對文獻(xiàn)的非正式調(diào)查表明佃声，典型的基于液滴的實(shí)驗(yàn)將捕獲10,000至100,000個(gè)細(xì)胞艺智，每個(gè)細(xì)胞以1,000至10,000個(gè)UMI進(jìn)行測序（通常與細(xì)胞數(shù)量成反比）』鳎基于液滴的方法還需要在通量和雙峰速率之間進(jìn)行權(quán)衡十拣，從而影響測序的真正效率。
對于涉及多個(gè)樣品或條件的研究志鹃，設(shè)計(jì)注意事項(xiàng)與批量RNA-seq實(shí)驗(yàn)的考慮相同夭问。每個(gè)條件應(yīng)有多個(gè)生物學(xué)重復(fù)，并且條件不應(yīng)與批次混淆曹铃。請注意缰趋，單個(gè)細(xì)胞不是重復(fù)單元。相反，我們指的是來自重復(fù)供體或培養(yǎng)物的樣品秘血。
獲取表達(dá)矩陣
來自scRNA-seq實(shí)驗(yàn)的測序數(shù)據(jù)必須轉(zhuǎn)換成可用于統(tǒng)計(jì)分析的表達(dá)矩陣味抖。考慮到測序數(shù)據(jù)的離散性灰粮，通常是一個(gè)計(jì)數(shù)矩陣仔涩，其中包含映射到每個(gè)細(xì)胞中每個(gè)基因的UMI或讀數(shù)的數(shù)量。量化表達(dá)的過程往往取決于技術(shù)：
1.對于10X Genomics數(shù)據(jù)粘舟，CellRanger軟件包提供了一個(gè)自定義管道來獲取計(jì)數(shù)矩陣熔脂。這使用STAR將reads與參考基因組比對，然后計(jì)算映射到每個(gè)基因的獨(dú)特UMI的數(shù)量柑肴。
2.偽比對方法（例如alevin）可更高效地獲取計(jì)數(shù)矩陣霞揉。這避免了精準(zhǔn)對齊的需要，從而減少了計(jì)算時(shí)間和內(nèi)存使用量晰骑。
3.對于其他高度復(fù)用的協(xié)議适秩，scPipe軟件包提供了更通用的管道來處理scRNA-seq數(shù)據(jù)。這使用Rsubread比對reads些侍，然后對每個(gè)基因的UMI進(jìn)行計(jì)數(shù)。
4.對于CEL-seq或CEL-seq2數(shù)據(jù)政模，scruff軟件包提供了專用的量化管道岗宣。
5.對于基于reads的方法，我們通沉苎可以重復(fù)使用處理大量RNA-seq數(shù)據(jù)相同的管道耗式。
6.對于涉及spike-in轉(zhuǎn)錄本的任何數(shù)據(jù)，在比對和定量過程中應(yīng)將spike-in序列作為參考基因組的一部分趁猴。
量化后刊咳，我們將計(jì)數(shù)矩陣導(dǎo)入R并創(chuàng)建SingleCellExperiment對象。這可以通過基本方法（例如read.table（））來完成儡司，然后再應(yīng)用SingleCellExperiment（）函數(shù)構(gòu)造娱挨。另外，對于特定的文件格式捕犬，我們可以使用DropletUtils（用于10X數(shù)據(jù)）或tximport / tximeta包（用于偽對齊方法）等專用方法跷坝。根據(jù)數(shù)據(jù)的來源，需要注意以下幾點(diǎn)：
1.某些feature計(jì)數(shù)工具會在計(jì)數(shù)矩陣中報(bào)告映射統(tǒng)計(jì)信息（例如碉碉，未對齊或未分配的reads數(shù)）柴钻。盡管這些值可用于質(zhì)量控制，但如果將其視為基因表達(dá)值垢粮，則會產(chǎn)生誤導(dǎo)贴届。因此，在進(jìn)行進(jìn)一步分析之前，應(yīng)將其刪除（或至少移至colData）毫蚓。
2.小心使用^ ERCC正則表達(dá)式來檢測人類數(shù)據(jù)中的spike-in行占键，其中計(jì)數(shù)矩陣的行名稱是基因符號。 ERCC基因家族實(shí)際上存在于人類注釋中绍些，因此這將導(dǎo)致錯(cuò)誤地將基因識別為spike-in轉(zhuǎn)錄本捞慌。通過使用帶有標(biāo)準(zhǔn)標(biāo)識符（例如Ensembl，Entrez）的計(jì)數(shù)矩陣柬批，可以避免此問題啸澡。
數(shù)據(jù)處理與下游分析
在最簡單的情況下，工作流程具有以下形式：
1.我們計(jì)算質(zhì)量控制指標(biāo)氮帐，以去除會干擾下游分析的低質(zhì)量細(xì)胞嗅虏。這些細(xì)胞在處理過程中可能已經(jīng)損壞，或者可能沒有被測序方案完全捕獲上沐。常用指標(biāo)包括每個(gè)細(xì)胞的總計(jì)數(shù)皮服，spike-in或線粒體reads的比例以及檢測到的feature的數(shù)量。
2.我們將計(jì)數(shù)轉(zhuǎn)換為標(biāo)準(zhǔn)化的表達(dá)值参咙，以消除特定于細(xì)胞的偏倚（例如捕獲效率）龄广。這使我們能夠在下游諸如聚類等步驟中在細(xì)胞間執(zhí)行準(zhǔn)確的比較。我們還應(yīng)用了一個(gè)轉(zhuǎn)換（通常是對數(shù)）來調(diào)整均值-方差關(guān)系蕴侧。
3.我們執(zhí)行feature選擇以選擇有興趣的特征子集進(jìn)行下游分析择同。這是通過對每個(gè)基因的細(xì)胞間差異建模并保留高度可變的基因來完成的。目的是減少不必要的基因的計(jì)算和噪聲净宵。
4.我們應(yīng)用降維來壓縮數(shù)據(jù)并進(jìn)一步降低噪聲敲才。通常使用主成分分析來獲得初始的低階表示，以進(jìn)行更多的計(jì)算工作择葡，然后再采用更具激進(jìn)的方法紧武，例如t-隨機(jī)鄰居嵌入可視化。
5.我們根據(jù)其（標(biāo)準(zhǔn)化）表達(dá)譜的相似性將細(xì)胞分組敏储。這旨在獲得用作不同生物學(xué)狀態(tài)的經(jīng)驗(yàn)代表分組阻星。我們通常通過識別細(xì)胞群之間差異表達(dá)的標(biāo)記基因來解釋這些分組。
諸如數(shù)據(jù)整合和細(xì)胞注釋之類的其他步驟以后再進(jìn)行討論已添。
快速上手的示例
在這里迫横，我們使用Macosko等人的基于液滴的視網(wǎng)膜數(shù)據(jù)集。（http://bioconductor.org/books/release/OSCA/overview.html#ref-macosko2015highly）酝碳，在scRNAseq包中有提供矾踱。 這個(gè)例子從計(jì)數(shù)矩陣開始，并以聚類結(jié)束疏哗，為生物學(xué)解釋做準(zhǔn)備呛讲。
library(scRNAseq)
sce <- MacoskoRetinaData()
# Quality control.
library(scater)
is.mito <- grepl("^MT-", rownames(sce))
qcstats <- perCellQCMetrics(sce, subsets=list(Mito=is.mito))
filtered <- quickPerCellQC(qcstats, percent_subsets="subsets_Mito_percent")
sce <- sce[, !filtered$discard]
# Normalization.
sce <- logNormCounts(sce)
# Feature selection.
library(scran)
dec <- modelGeneVar(sce)
hvg <- getTopHVGs(dec, prop=0.1)
# Dimensionality reduction.
set.seed(1234)
sce <- runPCA(sce, ncomponents=25, subset_row=hvg)
sce <- runUMAP(sce, dimred = 'PCA', external_neighbors=TRUE)
# Clustering.
g <- buildSNNGraph(sce, use.dimred = 'PCA')
colLabels(sce) <- factor(igraph::cluster_louvain(g)$membership)
# Visualization.
plotUMAP(sce, colour_by="label")
