Current best practices in single-cell RNA-seq analysis: a tutorial

Current best practices in single-cell RNA-seq analysis: a tutorial

摘要

这篇文章主要阐述了一个典型的单细胞 RNA-seq 分析的详细流程，包括了以下流程：

1. pre-processing -数据预处理

   • quality control -质量控制
   • normalization -数据标准化
   • data correction -数据修正（batch correction, noise correction, etc.）
   • feature selection -特征基因选择
   • dimensionality reduction -降维

2. cell-level downstream analysis -细胞水平的下游分析

3. gene-level downstream analysis -基因水平的下游分析

1. Pre-processing and visualization

experimental workflow

single-cell dissociation, single-cell isolation, library construction, and sequencing 分别为以上四步走

• As a first step, a single‐cell suspension is generated in a process called single‐cell dissociation in which the tissue is digested.
• plate-based & droplet-based 但是都存在问题，In both cases, errors can occur that lead to multiple cells being captured together (doublets or multiplets), non‐viable cells being captured, or no cell being captured at all (empty droplets/wells)
• Each well or droplet contains the necessary chemicals to break down the cell membranes and perform library construction. Furthermore, many experimental protocols also label captured molecules with a unique molecular identifier (UMI).

What is UMI?

UMIs allow us to distinguish between amplified copies of the same mRNA molecule and reads from separeate mRNA molecules transcribed from the same gene

• Libraries are labelled with cellular barcodes and pooled together(multiplexed) for sequencing

  Illumina 单细胞测序工作流程： https://links.jianshu.com/go?to=http%3A%2F%2Fweb.illumina.com.cn%2Flanding%2Fproducts_view.asp%3Fnewsid%3D324

1. Raw data

原始数据根据是否在单细胞文库构建方法中使用 UMIs一般分为count matrice 或read matrices。

Raw data processing pipelines (alignment and quantification):

• Cell Ranger (Zheng et al, 2017)
• indrops (Klein et al, 2015)
• SEQC (Azizi et al, 2018)
• zUMIs (Parekh et al, 2018)

2. Quality control

在对单细胞基因表达数据分析之前，我们必须确保所有的细胞 barcode 数据对应那些活着的细胞。即一个 barcode 是细胞的唯一身份证。

QC操作一般建立在三个 QC covariates 上：

• 每个 barcode 对应的基因数量
• 每个 barcode 对应的count数量（count depth）
• 每个 barcode 线粒体基因占 count 数的比例

通过一定的阈值，来筛选掉那些过低count（细胞已经死亡发生裂解了）或过高count（细胞的双峰或多峰）

但是需要注意的是仅仅单一考虑以上一种指标可能会对细胞信号出现一些误解，因此需要把数据和生物学知识结合起来

例如，相对高的线粒体基因占比可能是细胞正处于呼吸过程中；或者是低 count 或 低 gene 数的细胞对应着静态细胞群落，并且细胞有着更大 count 数的可能具有更大的体积。

因此需要协同考虑多种指标来进行全局的阈值选择，以此避免那些 viable 细胞被不小心筛选掉。

目前有一些 doublet 检测工具提供更好的解决方案如:

DoubletDecon: preprint: DePasquale et al, 2018

Scrublet: Wolock et al, 2019

Doublet Finder: McGinnis et al, 2018

除了以上对细胞完整性的检查之外，还必须对各个细胞在转录组水平上执行 QC 步骤。原始计数矩阵通常包含了 20,000 个基因，通过滤除大部分细胞均表达的基因，这些基因也正是那些对细胞异质性贡献度不那么大的基因，可以大大减少原始矩阵的基因数目。

设置此阈值的准则是目前预期想要得到最小细胞群落的大小，并为dropout效果留出余地。例如，去掉那些少于 20 个细胞表达的基因可能会使得更难检测出那些少于 20 个细胞的细胞群落。

进一步的质量控制能够直接进行再计数矩阵中，Ambient gene expression 指代那些 count 数并不源于 barcode 细胞的，却来自于在文库构建之前由于细胞裂解导致的细胞悬液增加了污染的 mRNA 数量，从而使得 count增大。但是这种问题能够利用 droplet-based scRNA-seq 通过那些大量的空液滴作为空白对照进行矫正。

（例如 SoupX preprint: Young & Behjati, 2018 能够直接应用于 count data）

我们知道质量控制是为下游分析服务的，一个数据质量好坏是由下游分析结果的性能紧密相关的（例如聚类中细胞群落的标注）。因此，我们在数据分析时往往需要多次进行质量控制的再回顾。从一个较为宽松的 QC 阈值开始，再逐渐严格通常来讲是比较有益的。

该方法对于包含异类细胞群体的数据集特别有用，在这些数据集中，细胞类型或状态可能会被误解为低质量的异常细胞。在低质量数据集中，可能需要严格的QC阈值。 QC阈值不宜用于改善统计测试的结果。相反，可以根据数据集可视化和聚类中QC协变量的分布来评估QC。注意 data peeking 的发生（Data peeking goes something like this: you run some subjects, then do an analysis. If it comes out significant, you stop. If not, you run some more subjects and try again）

3. Normalization

归一化和标准化的区别

https://cloud.tencent.com/developer/article/1486102

• 常用的归一化Normalize处理目的是将离散程度很大的数据集中化, 对数据转换能够让同一基因在不同样本具有可比性（RPKM/TPM） ，在Seurat中 LogNormalize 便是利用log1p(value/colSums[cell-idx] *scale_factor)
• 常用的标准化Scale是基于之前归一化结果，再添加z-score计算，考虑到了不同样本对表达量的影响，消除了表达的平均水平和偏离度的影响，即统一数据的均值和方差

归一化使用场景

• 对表达量的范围有要求
• 表达量稳定，无极端值
• 数据不太符合正态分布时，可以使用归一化
• 机器学习算法（SVM,KNN,神经网络等）要求input数据为归一化

标准化使用场景

• 离散程度大，存在异常值和较多噪音，用标准化可以避免异常值和极端值的影响

• 在分类、聚类、PCA算法中，使用z-score值结果更好

• pheatmap(dat) # scale之前
  n = t(scale(t(dat)))
  n[n>2]=2 # up-limitation
  n[n<-2]=-2 #down-limitation
  pheatmap(n,show_colnames = F,show_rownames = F)
  

3.1 Normalize cellular count data (to make it comparable between cells)

Count 在计数矩阵中有两个维度，一个是细胞数目，一个是基因数目。表达矩阵中，每个count值都需要细胞捕获、反转录、测序这三个步骤成功完成而记录下来。相同细胞的测序深度（检测出的基因数目）可能由于每个步骤中操作差异而不同。因此，在比较两个细胞时，其中的任何差异都有可能是由于实验测序误差产生的，而不是真实的生物学差异。

在这里首先讨论如何实现细胞之间的可比较性。通过归一化能够消除这种因为采样等操作造成的偏差，例如scaling 计数数据，以利用缩放后的数据反映细胞内部基因相对的表达丰度，利用细胞内部的相对值放在同一水平比较不同细胞的基因表达量能够减少实验测序误差，突出生物学差异。

最常见的标准化方法是count depth scaling，也被称作counts per million或者CPM normalization。这个方法常用于bulk转录组，它会根据每个细胞的总表达量计算一个 size factor ，然后对其中各个基因表达量进行normalize。在CPM方法中，有一个前提假设，即数据集中的所有细胞最初都含有相同数量的mRNA分子并且计数深度的差异仅由于取样操作。

由于单细胞数据是具有不同大小及分子数量的异质细胞群，更复杂的标准化方法应该适用。例如，Weinreb et al (2018) 利用了一种CPM拓展方法，在计算他们的size factor时，排除了那些占细胞分子总数至少5%的基因。主要是预防少量极高表达量基因的存在，从而显现其他基因的差异。Scran包有个pooling‐based size factor estimation方法，允许更高的细胞异质性存在；另外Scran包在批次矫正和差异分析环节也比其他归一化方法表现更好(Buttner et al, 2019)

scRNA-seq测序技术分为全长测序及3'端富集测序；全长测序考虑到了基因的长度，而后者则没有考虑。TPM归一化常用于全长测序数据。

3.2 Normalize (scale) gene counts (to make it comparable between the same genes expressed in different cells)

基因的归一化是将gene counts的数据通过放缩(scaling)创造出其均值为0，方差为1(z scores)的一组数据。这样处理的好处是所有的基因权重都是一样的，在下游分析中能够与其他基因平等的进行比较。目前对基因是否进行缩放尚未达成一致。（Seurat 选择进行缩放，而Slingshot则不选择进行缩放）对基因进行缩放意味着在下游分析中所有的基因都被平均权重，而不进行缩放则将基因表达的量级作为基因重要性的一种信息保留下来。为了尽可能保留多的生物学信息，本文避免对基因进行放缩。

对数据进行归一化之后，数据矩阵通常进行log(x+1)转化。这样的转换有三种重要的影响:

• 第一，对数转换后的基因表达量之间的距离反应出了对数倍数改变量，即

  [logFOLD=logx-logy ]

  它能够很好的反应出表达量发生的变化

• 第二，对数转换能够减缓单细胞数据中均值-偏差的关系

• 第三，对数转换可以减少数据的偏斜度（skewness），从而近似于许多下游分析工具对数据呈正态分布的假设

尽管scRNA-seq数据实际上不是呈对数正态分布，这三个影响使得对数转换成为稍许粗糙但非常有用的工具。

对数转换对于下游分析中检测差异表达基因及批次效应矫正(Finak et al, 2015; Ritchie et al, 2015)。但是也应该注意到对数转换的一个弊端，即给数据引入了虚假的差异表达效应(preprint: Lun, 2018)当标准化大小因子分布在测试组之间差异很大时，此效果尤其明显。

4. Data correction and integration

上述的归一化处理（normalization）试图消除采样导致的计数偏差，但是归一化后的数据可能仍然包含一些偏差。数据矫正的目的是进一步消除技术性或生物性的偏差，如批次效应、dropout、细胞周期效应。这些因素并不总是能得到矫正，因此需要根据预期的下游分析来考虑矫正哪个因素。本文建议分别考虑生物学和技术性变量的矫正，因为它们用于不同的目的。

4.1 Regressing out biological effects

虽然矫正技术性偏差对于发现潜在生物学信号可能至关重要，矫正生物学变量却能够揭示出特定的生物学信号。最常见的生物学数据矫正是在转录组中除去细胞周期的影响。可以通过Scanpy和Seurat平台中实现的细胞周期评分进行简单的线性回归或再具有更复杂混合模型的专用程序包中（如scLVM 或 fscLVM）。可以从文献中获得用于计算细胞周期评分的标记基因列表(Macosko et al, 2015)。这些方法能够用于消退其他已知的生物学效应，例如线粒体基因表达，这被解释为细胞压力的反应。

还有一些方面需要在矫正生物影响偏差之前考虑到。

1. 矫正生物影响并不能总是对解释单细胞测序数据有帮助。尽管除去了细胞周期的影响，提高了对发育周期推测的合理性，但是细胞周期信号对于生物学同样也能提供有益的信息。例如，可以根据细胞周期评分以确定增殖细胞群（proliferating cell populations）。
2. 除此之外，生物学信号必须在相应的环境背景下才能够理解。假定生物学过程发生在同一生物体内，则这些过程之间存在依赖性。因此，矫正了一个过程则可能误将另一个信号掩盖住。
3. 最后，细胞大小的变化也与细胞周期造成的转录组效应有关。因此，通过归一化或专用工具（如cgCorrect Blasi et al 2017 ）也能部分矫正细胞周期对单细胞测序数据的影响。

4.2 Regressing out technical effects

用于以上生物学偏差的回归模型的变体也能应用于技术性偏差。单细胞数据中最主要的技术性偏差是count depth和 batch。尽管归一化对数据进行放缩，使得细胞之间的基因count数能够在同一水平进行比较，count depth效应依旧存在于数据中。例如，细胞在大小上具有差异，因此它们的mRNA分子数也不尽相同。但是，技术效应在归一化后仍可能存在，因为没有一种放缩方法能够推测那些由于不良取样过程而未检测到的基因表达量。对计数深度影响进行回归能够提高轨迹推测算法（找到细胞之间的转换）。

一种基于回归策略消除count effect的替代方法是使用更严格的归一化过程，例如下采样或非线性归一化方法（请参见“归一化”部分）。这些方法可能与基于板的scRNA-seq数据集特别相关，其中每个细胞计数深度的较大变化可以掩盖细胞之间的异质性。

4.3 Batch effects and data integration

Batch effect

当将细胞分为不同的组别时，批次效应将会发生。这些组中的细胞可能来自不同芯片、不同测序通道或在不同时间节点收集得到的。细胞所经历的不同环境可能会影响转录组的测序或转录组本身。所产生的影响存在多个层面：实验中的细胞组之间，同一实验室中进行的实验之间或不同实验室的数据集之间。本文将对前两种情况作出区分。

校正在同一实验中的样品或细胞之间的批次效应是一种经典情况，称为bulk RNA-seq 中的批次效应。我们将其与来自多个实验的数据整合（data integration）区分开。通常使用线性方法校正批处理效果，而将非线性方法用于数据集成。

一个经典的批次效应矫正方法ComBat(Johnson et al, 2006)在中低复杂度的单细胞实验中性能也良好。**ComBat包含了基因表达的线性模型，其中在数据的均值和方差中均考虑到了批次效应的影响。无论采用何种计算方法，最好的批效应校正方法都是通过巧妙地设计实验来规避。可以通过将不同实验条件和样品中的细胞进行合并来避免批次效应。利用细胞标记(preprint: Gehring et al, 2018)，或者通过遗传变异(Kang et al, 2018)，也能够解构实验中合并的细胞。

Data integration

与批次效应相比，数据整合方法的挑战围绕数据集之间的差异。估计批处理效果时，ComBat使用在同一批次中的所有单细胞来对应批次参数。这种方法将批次效应与细胞培训或状态之间的生物学差异混淆。数据整合方法，例如:

Canonical Correlation Analysis CCA	Butler et al, 2018
Mutual Nearest Neighbours MNN	Haghverdi et al, 2018
Scanorama	preprint: Hie et al, 2018
RISC	preprint: Liu et al, 2018
scGen	preprint: Lotfollahi et al, 2018
LIGER	preprint: Welch et al, 2018
BBKNN	preprint: Park et al, 2018
Harmony	preprint: Korsunsky et al, 2018

尽管数据整合方法也能够被应用于简单的批次效应校正问题，考虑到非线性数据整合方法的自由度增加，我们建议要注意到过度校正的出现。例如，在更简单的批次校正设置中，MNN的表现

第一张图片每个细胞都按照样品来源进行着色，可以很明显看出样品与样品之间有分层，说明了该数据集有很明显的批次效应。而利用ComBat在小鼠肠上皮细胞数据集上进行批效应校正后，可以看见各类样品在UMAP投影中分布更均匀，批效应明显减弱。

Expression recovery

技术校正的进一步形式是表达恢复（expression recovery, also denoising or imputation）。单细胞转录组sexual包含了各种噪声源(Gru¨n et al, 2014; Kharchenko et al, 2014; Hicks et al, 2017).噪声的最主要方面主要是dropout造成的。现有的一些工具通过推测dropout事件，用合适的表达量去替代这些0值，进而降低数据集的噪音。

例如：MAGIC: van Dijk et al, 2018; DCA: Eraslan et al, 2019; scVI: Lopez et al, 2018; SAVER: Huang et al, 2018; scImpute: Li & Li, 2018 。

已被证实进行表达恢复能够提高基因与基因相关性的估计(van Dijk et al, 2018; Eraslan et al, 2019)。不仅如此，这个步骤也能够与归一化、批次效应校正、其他在scVI工具中实现的下游分析进行整合。尽管大多数据校正方法将归一化后的数据作为输入，一些表达恢复方法却基于预期的负二项式噪音分布，故输入为原始计数数据。

当应用表达恢复时，应该注意到没有一个方法是完美的。因此，任何方法都可能过校正或校正不足。实际上，已报道由于表达恢复造成了错误的相关信号。鉴于在实际应用中难以评估表达恢复的准确度，这种情况对用户思考是否对数据进行去噪处理构成了挑战。由于这些考量，目前尚无关于应如何对数据去噪的共识。谨慎而言，我们将仅在数据可视化时才使用表达式恢复，而在探索性数据分析过程中不使用假设的数据。因此全面的实验验证相当重要。

5. Feature selection, dimensionality reduction and visualization

一个人类的单细胞rna测序数据集能够包含高达25,000个基因的表达量。对于一个给定的scRNA-seq数据集，许多基因并不能提供有用的信息，并且大多基因会包含0值。尽管在QC步骤中除去了这些zero count gene，对于一个细胞向量而言，也有高达15,000个维度（基因）。为了减少下游分析工具的计算负担，减少数据中的噪声并使数据可视化，可以使用多种方法来减少数据集的维数。

5.1 Feature selection

降维的第一步通常是基因挑选。在这个步骤中，数据集将进行筛选以保留那些可以显示出数据集中差异的那一部分“informative”基因。因此经常使用highly variable genes(HVGs)。根据任务和数据集的复杂性，通常选择1000-5000个HVG用于下游分析。数量更多的HVG会呈现出更好的PCA结果。

在Scanpy和Seurat中都实现了一种简单而流行的选择HVG的方法。在此，通过基因的均值表达对基因进行分类，并在每个分类中选择方差均值最高的基因作为HVG。该算法存在不同的风格，对于Seurat则使用count data，而对于Cell Ranger 则使用对数转换过的数据。不过最好还是在技术偏差校正之后选择HVG，以避免选择仅由于例如批次效应导致的高度可变基因。Other methods for HVG selection are reviewed in Yip et al (2018).

5.2 Dimensionality reduction

在基因选择过后，单细胞表达矩阵的维度能够通过专门的降维算法进一步降低。这些算法能够使得表达矩阵嵌入在低维空间中，该空间旨在以尽可能少的维数捕获数据中的深层结构。其原理是单细胞RNA-seq数据本身具有低维性(Heimberg et al, 2016)；或者是可以用比基因数目少得多的维数来充分描述存在在细胞表达谱图上的生物学特征。降维就是找到这些特征维度。

降维的两个目标是为了可视化（visualization）和汇总（summarization）。可视化是试图将数据集以二维或三维投影的最佳方式展现出来。这些降维后的特征值作为散点图上的坐标，以获得数据的直观表示。汇总没有规定输出成分的数量，反而对于描述数据中存在的可变性，较高的分量变得不再重要。汇总技术可用于通过查找数据的固有维数将数据缩减为其基本组成部分，从而有助于进行下游分析。尽管对输出结果不进行二维可视化，但可以使用汇总方法利用top的简化后的成分进行可视化。但是专用的可视化技术通常能够更好地表示差异性。

降低的维度是通过特种空间维度（基因表达向量）的线性或非线性组合生成的。特别是在非线性情况下，在此过程可能会降低维度的可解释性。图4中显示了一些常用的降维方法的示例应用。本文主要阐述如何在常见的降维方法之间进行选择。A more detailed review of dimensionality reduction for single-cell analysis can be found in Moon et al (2018).

PCA

主成分分析时一种线性方法，通过最大化与其他每个维度中的残差来生成缩减的维度，即主成分。尽管PCA并没有在几个维度上或非线性方法中获得数据的结构，但它是许多聚类或轨迹推断分析工具的基础。事实上，PCA通常用作非线性降维方法的预处理步骤。PCA通过其前N个主成分来总结数据集，其中N的选定可以用"elbow"heuristics 或基于置换测试的jackstraw方法确定。PCA的简单线性性质有以下优势：

在降维空间中的距离与该空间所有区域中具有一致的解释性。因此我们可以将感兴趣的量与主成分相关联，以评估其重要性。例如可以将主成分投影到technical nuisance covariates （如 percent_mt, feature count）以研究QC的性能，数据矫正和归一化步骤。

Diffusion maps

扩散图是非线性数据汇总技术。由于扩散成分强调数据中的转移，通常，每个扩散成分（扩散图的维度）突出显示不同细胞群体的异质性。

5.3 Visualization

出于可视化的目的，标准实践手段是使用非线性降维方法。最常见的则为t分布随机邻近嵌入（t-SNE）。t-SNE维度着重于全局结构为代价获取局部相似性。因此，这些可视化会夸大细胞群体之间的差异，而忽略掉这些群体之间的潜在联系。另一个困难是其perplexity参数的选择，因为t-SNE图可能会根据其值显示出明显不同的簇数。t-SNE常见的替代方法是，Uniform Approximation adn Projection method(UMAP), 或者是基于图的工具（graph-based tool），例如SPRING(Weinreb et al, 2018)。UMAP和SPRING的force-directed 布局算法ForceAtlas2可以说是底层拓扑的最佳近似。

但是UMAP的优势在于它的速度和对大数据集的兼容。因此，在没有特定生物学问题的情况下，UMAP是最佳的探索性数据可视化工具。此外，UMAP还可以汇总二维以上的数据。

6. Stages of pre-processed data

尽管上述内容按顺序概述了scRNA-seq中常见的预处理步骤，但下游分析通常更喜欢采用不同层次的预处理数据，并且建议根据下游应用来适应预处理。

预处理分为5个数据处理阶段

1. raw data
2. normalized data
3. corrected data
4. feature-selected data
5. dimensionality-reduced data

这些数据处理的不同阶段被分为三个层次：measured data 用于统计检验, corrected data 用于数据比较可视化, reduced data 用于下游分析

细胞和基因的质量控制应始终执行。但后续操作可能不需要操作，比如单批次的数据不需要数据校正。

6.1 Downstream analysis

经过预处理后，使用称之为下有分析的方法去提取生物学信息并描述其潜在的生物学系统。这些描述通过将可解释模型拟合到数据而获得的。例如具有相似基因表达谱的细胞群，它们代表同一种细胞类型的群落。在相似细胞中基因表达的微小改变能够表示细胞间的连续（分化）轨迹；或者基因具有相关联的表达模式表明它们具有共调控关系。

下游分析可以分为细胞水平和基因水平两种方法

细胞层次的分析一般着重两个结构的描绘：聚类和轨迹分析。广义地讲，聚类分析方法试图根据细胞分类为组别来解释数据的异质性。而在轨迹分析中则恰恰相反，数据被认为是动态过程中的某一时刻的快照，可以帮助理解某个动态发育过程。基因层面是：差异分析，富集分析，互作网络

7. Cluster analysis

7.1 clustering

聚类是任何单细胞分析中第一个中间结果。聚类能够使得我们推测细胞的身份。通过基于细胞基因表达谱图的相似性给细胞分组能够得到细胞的类。表达谱图相似性通过距离矩阵来决定，距离度量通常将降维的表示形式作为输入。相似性评分一般是在PC降维后的表达空间上计算欧氏距离。存在两种根据这些相似性分数生成细胞簇的方法：clustering algorithms & community detection methods

Clustering

聚类是一个经典的非监督机器学习问题，直接借助距离矩阵进行分析。细胞通过最小化类别间距离或在减少的表达空间中找到密集区域，被归到相应的簇中。最常用的k均值聚类算法（k-means clustering）通过确定聚类中心并将细胞分配给最近的聚类中心，将细胞分为k个类别。质心位置经过迭代优化。这个方法需要输入预期的簇数，这通常是未知的往往需要启发式校准。k均值算法在单细胞数据中的应用在所使用的距离矩阵中有所不同。标准的欧氏距离的替代包括cosine similarity(余弦相似性)， correlation-based distance metrics(Kim et al, 2018), 或者SIMLR方法，这利用了Gaussian kernels(Wang et al, 2017)为每个数据集学习距离度量。最近的比较表明，correlation-based distance 表现更佳。

Community detection methods

社区检测方法是图形划分算法，因此依赖于单细胞数据的图形表示。使用K最邻近方法（KNN图）获得此图表示。细胞在图中表示为节点（node），每个细胞都与其K个最相似的细胞相连，这些细胞通常使用PC减少的表达空间上的欧式距离获得。根据数据集的大小，K通常设置为5到100个最近的邻居之间。结果图捕获了表达数据的拓扑结构。表达空间的密集采样区域表示为图的密集连接区域。使用社区检测方法检测这些密集区域，社区检测通常比聚类方法要快，因为只有相邻细胞对才被视为同种细胞类型。因此，这种方法大大减少了可能的群集的搜索空间，这是在Louvain算法（Blondel）中实现的等人，2008年）。

7.2 Cluster annotation

在基因水平上，通过找到每个群落的标记基因聚类数据能得以被分析。这些称之为marker genes 表征了每个群落，并被用来用有意义的生物学标记对其进行标记。该标签表示了群落中细胞的类型（身份）。由于任何聚类算法都会产生数据分区，因此只能通过其代表的生物学特性的成功注释来确定已识别聚类的有效性。

虽然很可能会尝试假定单细胞数据中的簇代表真实的细胞类型，但是有多个变量决定了细胞的类别。

1. 首先，并不总是清楚什么构成了一个细胞类别。例如，尽管“T细胞”可能是某些人满意的细胞类别标记，但其他人可能会在数据集中寻找T细胞亚型并区分CD4+和CD8+ T细胞
2. 此外，相同类型的细胞处于不同转台可能会在不同的群落中检测出

由于上述原因，细胞类型的定义存在争议。最好使用术语"cell identities"而不是“cell types”。因为就像上述提到的T细胞，T细胞本身是T细胞，但是在发育阶段和状态不一样，身份会发生改变 ，但不能简答理解为它的类型发生改变了。

识别和注释cluster依赖于外部信息，即那些描绘了单个细胞identities的预期表达模式。即reference database。这些数据库大大的促进了细胞identities的标注。在没有相关参考数据库的情况下，可以通过将数据衍生的标记基因与文献中的标记基因比较来注释细胞身份，或者直接将文献来源的标记基因的表达值可视化。应当注意的是，后一种方法将用户限制于对从bulk expression研究中得出的细胞类型的经典理解，而不是cell identities。此外，已经表明，常用的细胞标志物在限定细胞identities方面能力有限。

利用reference database 标注细胞群落的两种策略：

• 利用data-derived 标记基因
• 利用全基因表达谱

标记基因集能通过应用差异表达检测找出（differential expression DE testing）。一般而言，我们更关注于那些群落中上调基因。由于标记基因被认为是强有力的差异表达因素，简单统计检验（如Wilcoxon 秩和检验或t检验等）通过两个群落中基因表达差异的程度去排序基因。各自检测中排序靠前的基因被视作marker gene。可以通过比较来自数据集的标记基因和来自参考数据集的标记基因（通过富集测试，Jaccard索引或其他重叠统计数据）来注释聚类。

检测标记基因需要注意两个方面

1. 由于DE测试中，0假设使基因在两组之间具有相同的表达值分布，但是这两组是通过聚类得到的，本身即是存在差异的。即是是随机数据的聚类时，也能发现标记基因。因此制定合适的P值很关键，
2. 并且单一变量的簇注释方法在特定情况下不建议使用。

7.3 Compositional analysis

在细胞层次上，我们可以根据数据组成结构来分析聚类数据。成分数据分析主要围绕归入每个cell-identity 类别的细胞比例。这些比例可能会由于疾病而发生改变。例如，沙门氏菌感染已显示会增加小鼠肠上皮中肠上皮细胞的比例。

在单细胞数据集分析成分变化需要足够多的细胞数量以稳健地评估细胞类别的比例，并且足够的样本数量以评估细胞身份簇组成中预期的背景变化。

8. Trajectory analysis

8.1 Trajectory inference

细胞多样性不能充分被离散分类系统（如聚类）描述。驱动观测到的异质性发展的生物学过程是连续过程。因此，为了获得细胞身份之间的转换，分支分化过程或生物学功能的缓慢非同步变化，我们需要动态的基因表达模型，这类方法称之为轨迹推断（TI）

轨迹推断方法将单细胞数据解释为连续过程的快照（snapshot），通过寻找穿过细胞空间的路径来重建该过程，该路径将相邻细胞之间的转录变化降低至最低。细胞沿着这个路径的顺序被称之为pseudotime。尽管该变量与距离root细胞的转录距离有关，但通常被解释为发育时间的代名词。

从Monocle和Wanderlust建立TI领域以来，可用的方法数量激增。当前可用的TI方法在建模路径的复杂性方面有所不同。模型的范围从简单的线性或分叉轨迹到复杂的图形，树或多分叉轨迹。没有一种单独的方法对所有类型的轨迹都能表现出最佳性能。相反，应根据预期轨迹的复杂性选择TI方法。对比显示，Slingshot（Street等人，2018）优于其他方法（从线性模型到双叉模型和多叉模型）的简单轨迹。如果期望更复杂的轨迹，作者建议使用PAGA（Wolf等人，2019）。如果确切的轨迹模型是已知的，则可以选择使用更专业的方法来提高性能（Saelens等，2018）。通常，应使用替代方法来确定任何推断的轨迹，以避免方法偏差。

在典型的工作流程中，TI方法应该用于降维或校正后的数据。由于在细胞内同时发生多种生物学过程，因此消退其他过程的生物学效应以区别我们预期的轨迹是有效果的。

例如T细胞在成熟过程中可能正在经历细胞周期转换，此外由于几种性能高的TI方法依赖于聚类数据，因此TI通常在聚类之后进行。推断轨迹中的类可能表示稳定或亚稳定。随后，可以将RNA速度叠加到轨迹上以增加方向性。

需要注意的是推断的轨迹并不代表真实的生物过程，这些仅表示转录相似性。因此需要进一步信息对其验证。

8.2 Gene expression dynamic

一个推断的轨迹不是拟合转录噪音的结果而是在基因水平上分析该轨迹。在伪时间内平滑变化的基因表征了轨迹，可用于识别潜在的生物学过程。此外，这组与轨迹相关的基因被认为是包含调节建模过程的基因。Regulator 基因帮助我们了解生物过程的触发方式和原因，并代表潜在的药物靶标。

8.3 Metastable states

亚稳态是一个相对稳定但又会变化的状态

亚稳态状态轨迹的细胞水平分析研究了伪时间的细胞密度。假定细胞采样的过程是无偏的，沿轨迹的密集区域表示了优选的转录组状态。当将轨迹解释为时间过程时，这些密集区域可能代表例如发展中的亚稳态。我们可以通过绘制伪时间坐标的直方图来找到这些亚稳态。下图中B中，有很多中状态，但C中直方图认为这几个密集的区域才属于亚稳态。

整合分群与轨迹分析

clustering 和 trajectory 可以整合在一起展示结果。将分群的结果当成节点（node），将轨迹当成节点之间的桥梁（edge），所以动态与静态的数据结合在了一起，利用partition-based graph abstraction PAGA 工具可以得到图D效果。

结语

文章阐述了单细胞数据从获得到预处理到分析等等一整套流程，虽然精读了一遍，也学了很多背景知识，还是发现许多遗漏的点，等以后实战中遇上了再查漏补缺建立自己的体系认知。