尽管readsdata和countdata的测量噪声级别不同（生信宝典注：基于UMI的数据，获得的是分子计数，countdata，噪声更低），但在分析流程中的处理步骤是相同的。为简单起见，在本教程中，我们将数据称为计数矩阵（countdata）。在reads和计数矩阵的结果不同的地方，将会特别提到reads矩阵。
在分析单细胞基因表达数据之前，我们必须确保所有barcode都对应于有效细胞(viablecells，有活力的细胞)。质控有3个指标：
质控就是检查这3个指标的分布中是否存在异常峰并设置阈值去除。这些异常的barcodes可能对应于死细胞、细胞膜破损的细胞或doublets。比如，如果某个barcode对应的样品测到的分子总数低、检测到的基因数少、线粒体基因所占比例高，则表明该样品可能存在细胞膜破损导致细胞质RNA漏出，只有线粒体中的RNA保留了下来。相反，如果某个barcode对应的样品有异常高的总分子数和检测到的基因数，则有可能这个样品包含2个或以上细胞（doublets）。因此，高的总分子数通常用来过滤潜在的doublets。3个最近发表的doublets检测工具提供了更好的解决方案：DoubletDecon,Scrublet,DoubletFinder。
包含不同类型细胞混合体的数据集的每个QC指标可能会呈现多个分布聚集峰。例如，图2D显示了具有不同QC分布的两个细胞群。如果未执行先前的过滤步骤（请注意，CellRanger也会执行细胞QC），则应仅将最低总分子数和基因数的barcode视为无效细胞。另一个阈值设定准则是考虑所选阈值过滤掉的细胞比例。比如过滤高分子总数细胞时，过滤掉的细胞比例不应超过预期的doublet率。
除了检查细胞的完整性外，还必须在转录本水平上执行QC步骤。原始计数矩阵通常包含超过20,000个基因。可以通过滤除在多数细胞中不表达的基因(通常认为这些基因对细胞异质性分析贡献不大)，大大减少该数目。设置此阈值的准则是使用感兴趣的最小细胞亚群大小，并为dropout事件留出一些余地(生信宝典注：比最小细胞亚群大小数字再小一点)。例如，滤除掉在少于20个细胞中表达的基因可能会使鉴定少于20个细胞的细胞亚群变得困难。对于dropout比率高的数据集，此阈值也可能会使不太大的细胞亚群难以被检测到。阈值的选择应随待分析数据集中的细胞数量和既定的下游分析而定。
另外可以直接对计数矩阵(countmatrix)执行进一步的质量控制。输入型基因表达（Ambientgeneexpression）是指某个barcode检测到的转录本是源自其他细胞在建库前发生破裂而释放到细胞悬液中的mRNA。这些增加的外来计数会影响下游分析，如Marker基因鉴定或其他差异表达分析过程。由于基于液滴的scRNA-seq数据集中存在大量空液滴，因此可以通过空液滴建模分析细胞悬液中的RNA构成和丰度来校正这一影响。最近开发的SoupX使用这种方法直接校正countmatrix。另外，在下游分析中直接忽略这些有强影响的输入型基因也是处理这个问题的一个实用方法。
图2.质控变量的分布图(小鼠肠道上皮细胞数据集)(A)每个细胞检测到的总分子数的分布用直方图展示。右上角的子图展示了总分子数小于4,000的部分的分布。因为在count数为1200处有个峰，所以阈值设置为了1,500。(B)每个细胞检测到的基因数的直方图展示。在400个基因处有个小峰（噪音峰），阈值设置为700。（生信宝典注：这阈值设置的随意性也没谁了~~~）(C)每个细胞检测到的总分子数从高到底绘制rankplot，类似于CellRanger检测并过滤空液滴的log-logplot。在总分子数为1500时，存在一个快速降低的拐点，则1500为筛选阈值。(D)同时展示检测到的基因数（纵轴）、总分子数（横轴）和线粒体基因的比例(颜色)。线粒体基因只在低基因数和低总分子数的细胞中比例高。这些细胞已经被上面设置的总分子数和基因数阈值过滤掉了。质控参数联合可视化显示更低的基因数筛选阈值可能也是足够的。
计数矩阵中的每个数值代表细胞中一个mRNA分子被成功捕获、逆转录和测序（Box1）。由于每个操作步骤固有的可变性，即便同一个细胞测序两次获得的计数深度也可能会有所不同。因此，当基于原始计数数据比较细胞之间的基因表达时，得到的差异可能来自于技术原因。Normalization可以通过调整计数数据(scalingcountdata)等解决这一问题，以获得细胞之间可比的相对基因表达丰度。
由于单细胞数据集通常由大小和分子数不同的异质细胞群体组成，因此通常需要更复杂的标准化方法。例如，Weinrebetal对CPM算法进行了扩展，在计算sizefactors时排除在任何细胞中总计数大于5％的基因。这一方法屏蔽掉少数高表达基因对总体表达变化的影响。软件包Scran的pooling‐basedsizefactor方法对细胞异质性的影响处理更好。首先把细胞合并到一起避免technicaldropout效应，然后基于基因表达的线性回归模型估算sizefactor。这一方法允许细胞有少于50%的差异表达基因，并且在不同的测试评估研究中这一标准化方法都表现最好。评估发现，scran比其他标准化方法对后续批次校正和差异表达分析效果更好。并且在scran作者小范围的比较中也展现出这个方法稳定性比较好。
标准化后，数据矩阵通常进行log(x+1)转换。此转换具有三个重要作用。
尽管scRNA-seq数据实际上不是对数正态分布的，但这三个效果使对数转换成为一种粗略但有用的工具。这一有用性在下游差异表达分析或批次校正时有更好的体现。但是应该注意的是，数据的对数转换会在数据中引入虚假的差异表达结果。而且如果sizefactor在组间差异更大时影响尤其明显。
在校正特定的生物因素影响之前，应考虑几个方面。首先，校正生物学因素影响并不总是有助于解释scRNA-seq数据。虽然消除细胞周期影响可以改善对发育轨迹的推断，但细胞周期信号也可以提供有意义的生物学信息。例如，可以根据细胞周期评分确定增殖细胞群(在Github的案例研究中有讲)。同样，必须根据具体研究的对象理解生物学信号。假设生物过程发生在同一有机体内，则这些过程之间可能存在依赖性。因此，校正一个过程的影响可能会无意间掩盖另一个过程的信号。另外，也有研究者提出细胞大小对转录组的影响通常也归因于细胞周期不同。因此，通过标准化或专用工具（例如cgCorrect）在校正细胞大小的同时也部分校正了scRNA-seq数据中细胞周期的影响。
用于移除不相干生物因素影响的回归模型也可应用于校正实验技术因素的影响。单细胞数据中最突出的技术影响因素是测序深度和批次。尽管标准化可以使得细胞之间的基因定量数据可比，但测序深度的影响仍保留在数据中。这种测序深度效应既可以是生物学自身的影响，也可能是技术带来的影响。例如，细胞大小可能不同，因此mRNA分子数量也可能不同。而且，标准化后技术带来的计数影响可能依然存在，因为没有办法推断由于实验过程中的随机性带来的未检测到的基因的表达。对测序深度进行回归校正可以提高轨迹推断算法的性能，该算法依赖于查找细胞之间的过渡（具体见Github上的案例研究）。在校正多个协变量（干扰因素）（例如细胞周期和测序深度）时，应在单个步骤中对所有协变量（干扰因素）执行回归，以解决协变量（干扰因素）之间的依赖性。
消除测序深度影响的另一个策略是使用更严格的标准化程序，例如downsampling或非线性归一化方法(见前面标准化部分)。这些方法可能更适用于基于板的scRNA-seq数据集，因为细胞之间较大的测序深度差异可能掩盖细胞之间的异质性。
人单细胞RNA-seq数据集可包含多达25,000个基因的表达值。对于一个给定的scRNA-seq数据集，其中有许多基因都不能提供有用信息，并且大多只包含零计数。即使在QC步骤中滤除了这些零计数基因后，单细胞数据集的特征空间也可能超过15,000个维度（即还会剩余15,000多基因）。为了减轻下游分析工具的计算负担、减少数据中的噪声并方便数据可视化，可以使用多种方法来对数据集进行降维。
scRNA-seq数据集降维的第一步通常是特征选择。在此步骤中，对数据集基因进行过滤仅保留对数据的变异性具有信息贡献的基因（在数据中变异大的基因）。这些基因通常被定义为高变化基因（HVG，highlyvariablegenes）。根据任务和数据集的复杂性，通常选择1,000到5,000个HVG用于下游分析。Kleinetal.的初步结果表明，下游分析对HVG的数量不太敏感。在HVG数量从200到2,400之间选择不同的数目时，评估显示PCA结果相差不大。基于此结果，我们宁愿选择更多的HVG用于下游分析（erronthesideof：可以借鉴的英语短语）。
在Scanpy和Seurat中都实现了一种简单而流行的选择HVG的方法。在这里，基因按其均值表达进行分组，将每个组内方差/均值比最高的基因选为每个分组的HVG。该算法在不同软件中输入不同，Seurat需要原始countdata；CellRanger需要对数转换的数据。最佳地，应在技术等干扰因素校正之后选择HVG，以避免选择仅由于例如批次效应等引入的高可变基因。Yipetal.综述了其他HVG选择方法。
特征选择后，可以通过专用的降维算法进一步对单细胞表达矩阵进行降维。这些算法将表达式矩阵映射到低维空间中，同时以尽可能少的维数捕获数据中所有的信息。鉴于单细胞RNA测序数据固有的低维性特征，这一方法是合适的。也就是说，细胞表达图谱构成的生物形态(biologicalmanifold)可以使用远少于基因数目的维度信息来展示。降维旨在找到这些维度。
在细胞水平上经典可视化方法的替代方法是PAGA(partition‐basedgraphabstraction）。事实证明，此工具可以充分展示数据的拓扑结构，同时使用聚簇生成粗粒度的可视化图像。结合以上任何一种可视化方法，PAGA可以生成粗粒度的可视化图像（coarse‐grainedvisualizations），从而可以简化单细胞数据的解释，尤其是在测序细胞量大或整合了大量细胞的情况下。
虽然我们在上面按顺序概述了scRNA‐seq中常见的预处理步骤，但下游分析通常需要不同级别的预处理数据，建议根据下游应用调整预处理的各个步骤。为了向新用户说明这种情况，我们将预处理分为五个数据处理阶段：（i）原始数据，（ii）标准化后的数据，（iii）校正后的数据，（iv）特征选择后的数据，以及（v）降维后的数据。数据处理的这些阶段可以归类为三个预处理层：测量的数据，校正的数据和降维的数据。细胞和基因的质量控制筛选是必须执行的步骤，因此在此处略过。尽管层的顺序代表了scRNA-seq分析的典型工作流程，但也可以跳过某一步或在处理阶段的顺序中稍作更改。例如，单批次数据集可能不需要批次校正。在表1中，我们总结了预处理数据的每一层所适用的下游处理程序。
表1
表1中的预处理阶段分为三组：测量的数据，校正后的数据和降维后的数据。我们将测量的数据定义为原始数据和处理后但保留了表达值为零的基因的数据。应用细胞特异的sizefactor的全局标准化方法即使在log(x+1)转换后仍保留零表达值。相比之下，校正后的数据会去除不需要的变化信息进而替换零表达值。校正后的数据层表示数据的“最干净”版本，它是数据代表的生物信号的最近似值。我们称最终预处理层为降维后的数据（reduceddata）。此数据层强调数据的主要差异，可以使用降维后的特征集来描述。
前述特性决定了预处理数据对特定下游应用的适用性。作为最后的预处理阶段，降维后的数据是广泛应用的数据层。但是，差异表达分析只能在基因空间内进行生物学解释，而无法（或无法完全）体现在降维后的数据中。降维后数据的优势在于生物数据信息的汇总(summarization)和影响生物数据的噪声的降低。因此，降维后的数据应用于后续的探索性方法（如可视化、邻域图推断、聚类）以及计算复杂的下游分析工具（如轨迹推断）。实际上，许多轨迹推断方法在工具本身中都包含了降维功能。
图EV2.批次校正和降噪后变异系数的变化负值代表数据校正后变异系数降低。第一行展示的是ComBat校正(A)mouseintestinalepithelium(mIE)和(B)mouseembryonicstemcell(mESC)数据变异系数的变化。第二行展示的是DCA降噪后(C)mIE和(D)mESC数据变异系数的变化。
下游分析可分为细胞水平和基因水平的方法，如图5所示。细胞水平分析通常着重于两种结构的描述：簇和轨迹。这些结构又可以在细胞和基因水平上进行分析，即簇分析和轨迹分析方法。广义地讲，簇分析方法试图将细胞分类为离散的多个组来解释数据的异质性。相反，在轨迹分析中，数据被视为细胞连续变化动态过程的一个个快照，轨迹分析方法研究了这一连续变化过程。
图5.下游分析方法总览（蓝色背景的方法是基因水平的分析方法）。
在此，我们在详细描述独立于这些细胞结构进行的基因水平分析之前，先描述细胞水平和基因水平的聚类和轨迹分析工具。
社群检测方法是一种图分类算法，依赖于单细胞数据的网络图表示。图是使用K最近邻方法（KNN图）获得的。细胞在图中表示为节点。每个细胞与其K个最相似的细胞相连，通常对主成分降维后的数据计算欧氏距离作为相似性度量。根据数据集的大小，K通常设置为5到100个最近的邻居。获得的KNN结果图捕获了表达图谱的基础拓扑结构（Wolfetal.,2019）。表达谱相似的细胞集合对应于图的紧密连接区域，然后使用社群检测方法检测图中这些紧密连接区域。社群检测方法通常比聚类方法要快，因为只有相邻的细胞对会被视为属于同一群集。因此，这种方法大大减少了可能的细胞簇的搜索空间。
在开创性的PhenoGraph方法（Levineetal.，2015）出现之后，聚类单细胞数据集的标准方法已演变成多分辨率模块优化算法(multi‐resolutionmodularityoptimization），即在单细胞KNN图中应用Louvain算法。此方法是在Scanpy和Seurat单细胞分析平台中应用的默认聚类方法。评估发现，这一方法应用于单细胞RNA测序数据以及流式细胞和质谱数据分析时优于其他聚类方法。从概念上讲，Louvin算法将组内细胞连接数高于预期的一组细胞视作一个簇。(生信宝典注：Louvainalgorithm算法采用迭代方式，先计算每个点加入到其邻居节点所在社区带来的模块性评估收益，然后将其加入收益最大的邻居节点所在社区，对所有未归类点重复进行这一过程，直至结果稳定。然后再把每个社区作为一个节点，重复上一步。)优化后的模块鉴定函数包括一个分辨率参数，该参数允许用户确定簇的近似大小。通过对KNN图取子集，还可以只聚类一部分特殊的细胞簇。这种子集聚簇模式可以允许用户识别某个细胞类型簇内的细胞状态，但也可能会发现仅由数据中的噪声带来的聚类模式。
尽管把单细胞数据中鉴定到的簇归类为一个有生物意义的细胞类型是一个诱人的结果，但是细胞身份的确定不是那么简单（Wagneretal.,2016；Cleversetal.,2017）。首先，并不总是清楚细胞类型定义的精确程度。例如，“T细胞”可能是某些人满意的细胞类型标记，但其他人可能会在数据集中寻找T细胞亚型并区分CD4+T细胞和CD8+T细胞。另外，同一细胞类型的不同状态可能会被分到不同的细胞簇。基于上述原因，最好使用术语cellidentities而不是celltypes。在对细胞进行聚类和注释之前，用户必须确定要注释到的细胞信息精确水平，从而确定合适的聚簇分辨率。
图6.聚类分析结果
(A)UMAP可视化Louvain算法鉴定的细胞簇，并标记注释的细胞类型名称。(B)通过细胞类型Marker基因的表达标记细胞簇所代表的细胞类型，如干细胞(Slc12a2),肠上皮细胞(Arg2),杯状细胞(Tff3)和潘氏细胞(Defa24)。灰色代表低表达，红色代表高表达。Marker基因也可能在其它细胞簇有表达，比如Tff3在杯状细胞和潘氏细胞都有表达，只是在杯状细胞更高。(C)肠上皮细胞区域细胞类型组成热图（上面是近端区域，下面是远端区域）。相对细胞密度高的区域用深红色表示。
鉴定标记基因时应注意两个方面。首先，用于鉴定标记基因的P值基于以下假设：即鉴定出的细胞簇是有生物学意义的。如果细胞簇鉴定过程中存在不确定性，则必须在统计检验中考虑簇分配与标记基因鉴定之间的关系。这个关系起因于鉴定细胞簇和标记基因都是基于同一套基因表达数据。差异基因检测的零假设（nullhypothesis）是两组细胞整体基因的表达值具有相同的分布。然而，由于这两个聚类组是基于基因表达变化的聚类结果得到的，其基因表达谱从本质上肯定存在差异。因此即使在应用splatter随机生成的数据的聚类结果中也可以鉴定出显著的Marker基因(Zappiaetal,2017)。为了在聚类数据中获得合理的显著性度量，可以使用置换检验减少聚类步骤带来的影响。补充介绍S3中对这一检验有详细描述。最近的一个差异表达工具也专门解决了这个问题（Preprint：Zhangetal，2018）。
在当前的研究中，P值通常会被夸大，这可能导致对标记基因数量的高估。但是，基于P值对基因进行排序则不受影响。假设聚类结果是有生物学意义的，那么排名最高的标记基因仍将是最佳候选标记基因。首先，我们可以通过可视化展示粗略地查验获得的标记基因。我们强调，特别是通过无监督聚类方法定义细胞聚类簇时，会导致夸大的P值。当改为通过单个基因的表达确定细胞簇的身份时，P值可以解释为相对于其他基因的期望值。尽管很常见，但是这种单变量的聚类簇注释方法除了在特定情况下不建议使用（例如，β细胞中的胰岛素或红细胞中的血红蛋白）。其次，标记基因在数据集中将一个细胞簇与其他细胞簇区分开，不仅依赖于细胞簇鉴定结果，还依赖于数据集的组成。如果数据集组成不能准确表示背景基因表达，则检测到的标记基因将偏向在其它细胞中未检测到的基因。特别是在细胞多样性低的数据集中计算标记基因时，必须特别考虑这一方面。
最近，自动细胞簇注释已可用。通过直接将注释好的参考簇的基因表达谱与单个细胞进行比较，使用诸如scmap（Kiselevetal.,2018b）或Garnett（Preprint：Plineretal.,2019）之类的工具可以在参考集和自己的数据集之间比较注释。因此，这些方法可以同时执行细胞注释和细胞簇鉴定，而无需数据驱动的细胞聚类步骤。但由于不同实验条件下细胞类型和状态组成不同（Segerstolpeetal.，2016;Tanay＆Regev，2017），基于参考数据的聚类不应取代数据驱动的聚类过程。
聚类、聚类注释、重新聚类或子聚类以及重新注释的迭代是很耗时的过程。自动化的细胞簇注释方法大大加快了此过程。但是，自动和手动方法各自存在优点和局限性，速度的提高与灵活性的降低并存，因此很难有一种全局适应的方法。如前所述，参考数据集很难与正在研究的数据集包含完全一样的细胞类型。因此，不应放弃用于手动注释标记基因的计算。特别是对于包含许多细胞簇的大型数据集，当前的最佳实践是两种方法（自动+手动）的组合。为了提高速度，可以使用自动化细胞类型注释工具粗略地标记细胞并确定可能需要鉴定细胞子群的簇。随后，应针对细胞簇计算标记基因，并将其与参考数据集或文献中的已知标记基因集进行比较。对于较小的数据集和缺少参考数据库的数据集，手动注释就足够了。
细胞组成分析（Compositionalanalysis）。在细胞层面，我们可以根据其组成结构分析聚类数据。细胞组成数据分析围绕不同样品落入每个细胞簇的细胞比例进行分析。这一比例可能在疾病状态下发生变化。例如，研究表明沙门氏菌感染会增加小鼠肠上皮区域肠上皮细胞（enterocytes）的比例（Haberetal.，2017）。
轨迹推断。细胞多样性不能通过离散的分类系统（例如细胞聚类）充分描述。驱动观察到的细胞异质性发展的生物进程是一个连续过程（Tanay＆Regev，2017）。因此，为了捕获细胞身份之间的过渡状态、不同的分化分支或生物学功能的渐进式非同步变化，我们需要动态的基因表达模型。这类方法称为轨迹推断（trajectoryinference，TI）。
图7.小鼠肠上皮细胞轨迹分析结果展示
(A)Slingshot鉴定出的远端和近端肠上皮细胞分化轨迹。远端谱系细胞按pseudotime从红色到蓝色上色。数据集中其它细胞用灰色点展示。(B)PCA空间中Slingshot轨迹和细胞簇比较展示。簇名字简写如下：cuEP—enterocyteprogenitors;Imm.Ent.—immatureenterocytes;Mat.Ent.—matureenterocytes;Prox.—proximal;Dist.—distal.(C)远端肠上皮细胞轨迹中每个pseudotime对应的细胞密度分布。柱子按每个pseudotime区间主要的细胞簇上色。(D)数据集投射到UMAP空间的抽象图展示。每个簇用不同颜色点展示。出现在其它轨迹中的簇特别标记处理用于比较分析。TA代表短暂增殖下拨(transitamplifyingcells)。(E)R包GAM展示肠上皮细胞中基因表达随pseudotime的动态变化。
细胞水平联合分析。聚类和轨迹推断代表单细胞数据的两个不同视角，并可以在粗粒度(coarse‐grained)的图形展示中进行统一。用点代表单细胞簇，轨迹作为簇之间相连的边，可以同时展示数据的静态和动态属性。这一联合分析在PAGA工具中首先提出。PAGA提出一个统计模型进行细胞簇互作分析，将细胞相似度高于期望值的两个簇之间用线连接。在最近的综述中(Saelensetal,2018)，PAGA相比于其它轨迹分析方法整体表现更优。PAGA是唯一一个讨论到的可以处理不相连的拓扑结构和包含环形结构的复杂轨迹图的一个工具。这一特征使得PAGA在可视化整个数据集并进行探索分析时很有用。
基因水平分析。到目前为止，虽然我们专注于表征细胞结构的基因水平分析方法，但是单细胞数据的基因水平分析有更广泛的内容。差异表达分析、基因集分析和基因调控网络分析都可以直接研究数据中的分子信号。这些方法不是描述细胞异质性，而是把这种异质性作为理解基因表达差异的背景。
我们这儿描述的场景中，实验条件协变量是在实验设计中决定的。因此在同一簇内基于这一协变量的差异基因分析是独立于聚类过程的。这将基于实验条件的差异基因分析和基于细胞簇的差异基因分析区别开来。基于实验条件的差异基因检测获得的P-values是显著性检验统计指标，需要进一步做多重检验校正。为了降低多重检验校正的压力，不感兴趣的基因通常需要预先排除掉。虽然假基因和非编码基因也可以提供有用信息，但在分析中也通常会被忽略掉。
基因集分析。基因水平的分析通常会获得一长串难以解释其生物含义的基因列表。比如，处理组和对照样品通常会有数以千计差异表达基因。我们可以通过基因共有的特征对基因进行分组并检验这些特征在候补基因列表中的出现是否显著富集(overrepresented)以便更好的解释结果。
分析领域富集分析的一个新操作是基于配对基因标记的配体-受体分析。通过受体和共轭配体之间的表达关系推断细胞簇之间的互作。CellPhoneDB提供了受体-配体对信息，并用统计模型解释簇之间的高表达基因是否可以配对(Zeppetal,2017;Zhouetal,2017;Cohenetal,2018;Vento‐Tormoetal,2018)。
虽然存在专门针对scRNA-seq数据开发的GRN推断方法（SCONE：Matsumotoetal.,2017;PIDC：Chanetal.,2017;SCENIC：Aibaretal.,2017），但最近的比较分析表明普通和单细胞转录组的方法在这些数据表现都不好（Chen＆Mar，2018）。GRN推断方法仍可提供识别生物过程的因果调节子的有用信息，但我们建议谨慎使用这些方法。
分析平台。单细胞分析流程是多个独立开发的工具的集合。为了促进这些工具之间的数据传递，单细胞工具开发时都考虑使用一致的数据格式。这些工具为构建分析流程提供了基础。目前可用的分析平台有基于R（McCarthyetal.,2017;Butleretal.,2018）或Python（Wolfetal.,2018）的命令行，或具有图形用户界面（GUI）的本地应用程序（Patel，2018;Scholzetal.,2018）或Web服务（Gardeuxetal.,2017;Zhuetal.,2017）。Zhuetal.（2017）和Zappiaetal.（2018）对这些分析平台进行了综述。
推荐阅读
LueckenMD,TheisFJ.Currentbestpracticesinsingle-cellRNA-seqanalysis:atutorial.MolSystBiol.2019Jun19;15(6):e8746.doi:10.15252