概述了典型的scRNA-seq分析工作流程的框架（图1）。

image.png

图1：典型的scRNA-seq分析流程的示意图。每个阶段（用虚线分隔）由许多特定步骤组成，其中许多步骤可操作并修改。

实验设计

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

• 基于液滴的：10X基因组学，inDrop，Drop-seq
• 基于板的独特分子识别符（UMI）：CEL-seq，MARS-seq
• 基于板的读取：Smart-seq2
• 其他：sci-RNA-seq，Seq-Well

这些方法中的每一种都有其优点和缺点，在其他地方都进行了广泛讨论（Mereu等人，2019； Ziegenhain等人，2017）。在实际应用中，基于微滴的技术是当前由于其吞吐量和每个单元的低成本而成为标准配置。基于板的方法可以捕获其他表型信息（例如形态），并且更适合定制。基于读取的方法可提供完整的转录本覆盖，这在某些应用程序（例如剪接，外显子组突变）中很有用；否则，基于UMI的方法将减轻PCR扩增噪声的影响，因此更为流行。

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

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

获取计数矩阵

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

• 对于10X Genomics数据，该CellRanger软件套件提供了自定义管道以获取计数矩阵。它使用STAR将读段与参考基因组比对，然后计算映射到每个基因的独特UMI的数量。
• 伪对齐方法（例如）alevin可用于从相同数据中以更高的效率获得计数矩阵。这避免了显式对齐的需要，从而减少了计算时间和内存使用量。
• 对于其他高度复用的协议，scPipe软件包提供了更通用的管道来处理scRNA-seq数据。这使用Rsubread对齐器对齐读取，然后对每个基因的UMI进行计数。
• 对于CEL-seq或CEL-seq2数据，scruff软件包提供了专用的量化管道。
• 对于基于读取的协议，我们通常可以重复使用相同的管道来处理大量RNA-seq数据。
• 对于涉及刺入转录本的任何数据，在比对和定量过程中应将刺入序列作为参考基因组的一部分。

量化后，我们将计数矩阵导入R并创建一个SingleCellExperiment对象。这可以通过基本方法（例如read.table()）完成，然后应用SingleCellExperiment()构造函数。另外，对于特定的文件格式，我们可以使用DropletUtils（用于10X数据）或tximport / tximeta包（用于伪对齐方法）中的专用方法。根据数据的来源，这需要保持警惕：

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

数据处理和下游分析

在最简单的情况下，工作流程具有以下形式：

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

快速入门

在这里，我们使用Macosko等人的droplet-based的视网膜数据集。（2015），在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')
sce$clusters <- factor(igraph::cluster_louvain(g)$membership)

# Visualization.
plotUMAP(sce, colour_by="clusters")

image.png

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