概述

本章概述了典型的scRNA-seq分析工作流程的框架（图1）。 随后将更详细地描述每个分析步骤。

图1：典型的scRNA-seq分析工作流程的示意图。 每个阶段（由虚线分隔）由许多特定步骤组成，其中许多步骤都在SingleCellExperiment实例上进行操作和修改。

实验设计

在开始分析本身之前，对实验设计进行一些讨论可能会有所帮助。最明显的问题是技术的选择，大致可以分为：

1.基于液滴的：10XGenomics，inDrop, Dropseq
2.基于板的独特分子识别符（UMI）：CEL-seq，MARS-seq
3.基于板的读取：Smart-seq2
4.其他：sci-RNA-seq，Seq-Well
这些方法中的每一种都有其优点和缺点，其他地方对此进行了广泛讨论（Mereu等人(http://bioconductor.org/books/release/OSCA/overview.html#ref-mereu2019benchmarking); Ziegenhain等人2017(http://bioconductor.org/books/release/OSCA/overview.html#ref-ziegenhain2017comparative)）。实际上，基于液滴的技术由于其吞吐量和每个细胞的低成本而成为当前的事实上的标准。基于板的方法可以捕获其他表型信息（例如形态），并且更适合定制。基于reads的方法提供了完整的转录本覆盖范围，这在某些应用（例如剪接，外显子组突变）中很有用; 基于UMI的方法将减轻PCR扩增噪声的影响，因此更为流行。方法的选择取决于具体情况——但是以下我们分析流程中的大多数方面都与技术无关。

下一个问题是应该捕获多少个细胞，以及应该对它们进行测序的深度。简短的答案是“尽可能多地花钱”。长答案是，这取决于分析的目的。如果我们旨在发现稀有的细胞亚群，那么我们需要更多的细胞。如果我们旨在刻画细微的差异，那么我们需要更多的测序深度。目前，对文献的非正式调查表明，典型的基于液滴的实验将捕获10,000至100,000个细胞，每个细胞以1,000至10,000个UMI进行测序（通常与细胞数量成反比）。基于液滴的方法还需要在通量和双峰速率之间进行权衡，从而影响测序的真正效率。

对于涉及多个样品或条件的研究，设计注意事项与批量RNA-seq实验的考虑相同。每个条件应有多个生物学重复，并且条件不应与批次混淆。请注意，单个细胞不是重复单元。相反，我们指的是来自重复供体或培养物的样品。

获取表达矩阵

来自scRNA-seq实验的测序数据必须转换成可用于统计分析的表达矩阵。考虑到测序数据的离散性，通常是一个计数矩阵，其中包含映射到每个细胞中每个基因的UMI或读数的数量。量化表达的过程往往取决于技术：

1.对于10X Genomics数据，CellRanger软件包提供了一个自定义管道来获取计数矩阵。这使用STAR将reads与参考基因组比对，然后计算映射到每个基因的独特UMI的数量。
2.伪比对方法（例如alevin）可更高效地获取计数矩阵。这避免了精准对齐的需要，从而减少了计算时间和内存使用量。
3.对于其他高度复用的协议，scPipe软件包提供了更通用的管道来处理scRNA-seq数据。这使用Rsubread比对reads，然后对每个基因的UMI进行计数。
4.对于CEL-seq或CEL-seq2数据，scruff软件包提供了专用的量化管道。
5.对于基于reads的方法，我们通常可以重复使用处理大量RNA-seq数据相同的管道。
6.对于涉及spike-in转录本的任何数据，在比对和定量过程中应将spike-in序列作为参考基因组的一部分。
量化后，我们将计数矩阵导入R并创建SingleCellExperiment对象。这可以通过基本方法（例如read.table（））来完成，然后再应用SingleCellExperiment（）函数构造。另外，对于特定的文件格式，我们可以使用DropletUtils（用于10X数据）或tximport / tximeta包（用于伪对齐方法）等专用方法。根据数据的来源，需要注意以下几点：

1.某些feature计数工具会在计数矩阵中报告映射统计信息（例如，未对齐或未分配的reads数）。尽管这些值可用于质量控制，但如果将其视为基因表达值，则会产生误导。因此，在进行进一步分析之前，应将其删除（或至少移至colData）。
2.小心使用^ ERCC正则表达式来检测人类数据中的spike-in行，其中计数矩阵的行名称是基因符号。 ERCC基因家族实际上存在于人类注释中，因此这将导致错误地将基因识别为spike-in转录本。通过使用带有标准标识符（例如Ensembl，Entrez）的计数矩阵，可以避免此问题。

数据处理与下游分析

在最简单的情况下，工作流程具有以下形式：
1.我们计算质量控制指标，以去除会干扰下游分析的低质量细胞。这些细胞在处理过程中可能已经损坏，或者可能没有被测序方案完全捕获。常用指标包括每个细胞的总计数，spike-in或线粒体reads的比例以及检测到的feature的数量。
2.我们将计数转换为标准化的表达值，以消除特定于细胞的偏倚（例如捕获效率）。这使我们能够在下游诸如聚类等步骤中在细胞间执行准确的比较。我们还应用了一个转换（通常是对数）来调整均值-方差关系。
3.我们执行feature选择以选择有兴趣的特征子集进行下游分析。这是通过对每个基因的细胞间差异建模并保留高度可变的基因来完成的。目的是减少不必要的基因的计算和噪声。
4.我们应用降维来压缩数据并进一步降低噪声。通常使用主成分分析来获得初始的低阶表示，以进行更多的计算工作，然后再采用更具激进的方法，例如t-随机邻居嵌入可视化。
5.我们根据其（标准化）表达谱的相似性将细胞分组。这旨在获得用作不同生物学状态的经验代表分组。我们通常通过识别细胞群之间差异表达的标记基因来解释这些分组。
诸如数据整合和细胞注释之类的其他步骤以后再进行讨论。

快速上手的示例

在这里，我们使用Macosko等人的基于液滴的视网膜数据集。（http://bioconductor.org/books/release/OSCA/overview.html#ref-macosko2015highly），在scRNAseq包中有提供。 这个例子从计数矩阵开始，并以聚类结束，为生物学解释做准备。

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")

视网膜数据集的UMAP图，其中每个点都是一个细胞，并通过群集标识进行着色。