Background

From: Luecken, Malte D, and Fabian J Theis. 2019. “Current Best Practices in Single‐Cell RNA‐Seq Analysis: a Tutorial.” Molecular Systems Biology 15(6): 360–23.

Experimental workflow

Typical workflows incorporate single-cell dissociation, single-cell isolation, library construction, and sequencing.

– Dissociation: As a first step, a single-cell suspension is generated in a process called single-cell dissociation in which the tissue is digested.

– Isolation: To profile the mRNA in each cell separately, cells must be isolated. Single-cell isolation is performed differently depending on the experimental protocol.While plate-based techniques isolate cells into wells on a plate, droplet-based methods rely on capturing each cell in its own microfluidic droplet.

– Library construction: Each well or droplet contains the necessary chemicals to break down the cell membranes and perform library construction. Library construction is the process in which the intracellular mRNA is captured, reverse-transcribed to cDNA molecules and amplified.As cells undergo this process in isolation, the mRNA from each cell can be labelled with a well- or droplet-specific cellular barcode. Further- more, many experimental protocols also label captured molecules with a unique molecular identifier (UMI).

– Sequencing: These libraries are pooled together (multiplexed) for sequencing. Sequencing produces read data, which undergo quality control, grouping based on their assigned barcodes (demultiplexing) and alignment in read processing pipelines. For UMI-based protocols, read data can be further demultiplexed to produce counts of captured mRNA molecules (count data).

Data analysis workflow

Overview of the data analysis workflow:

Pre-processing layer content

Downstream analysis

Pre-processing content and method

QC

细胞层面QC最主要需要考虑的三个方面：

1）通过基因数量、计数深度和线粒体reads分数的异常峰来执行QC。需要同时考虑这些因素，而不是单独考虑。

2）尽可能考虑宽松的QC阈，如果下游聚类无法解释，则重新QC。

3）如果QC协变量在样品之间的分布不同，则应针对每个样品分别确定QC阈值，以解释样品质量差异，如Plasschaert等(2018)文章提到的。

检查这几个QC的指标分布，有助于得到outlier，并通过阈值进行过滤。这些outlier的barcode主要是源于死细胞，膜破损细胞，或者一个系统中多细胞。比如，一个barcode的count很少，基因检测数量也很少，但是线粒体基因比例很高，就很有可能是因为其膜破裂了，导致细胞质中的mRNA漏出了，而存在于线粒体中的遗传物质仍然得到了保留。相对的，如果一个barcode的count和检出的基因都异常很多，则可能是由于一个体系中存在多个细胞。

单独考虑这三种QC条件中的某一种，可能会导致对细胞信号的误判。细胞间的分子count可能会差异比较大，因此QC条件需要联合应用，并且这些条件应该设置的尽可能宽松，以免无意中过滤掉了活细胞群体。比如：细胞具有相对高的线粒体count，可能是由于参与呼吸过程；细胞具有低的count与gene数，可能是由于其为休眠的细胞群落；细胞具有高的count可能是由于细胞体积更大。

转录本层面也需要QC：

原始数据初始检出基因一般大于2w，但是通过过滤掉不在几个细胞中表达的基因（因此种基因也对于细胞异质性分析提供不了有效信息），最终的基因数就会大幅减少。设定最少细胞数量的阈值，一般是考虑兴趣点，同时保留一些dropout effect。比如，如果过滤掉了少于20个细胞表达的基因，就会使检出少于20cell的细胞群落变的很困难，结果中可能就会遗漏这种小群落的信息。但是对于high dropout rate的数据，20细胞的阈值也可能使鉴定大细胞群变的困难。因此阈值的选择，要基于数据集的实际情况，以及考虑后续分析的重点。

也可以基于基因表达计数文件做一些其他QC：Ambient gene expression（环境基因表达）是指不是来自barcode细胞，而是来自其他溶解细胞的count，这些细胞的 mRNA 在文库构建之前污染了细胞悬液，这些污染会歪曲下游分析结果，如标记基因鉴定或其他差异表达检测（特别是在样本水平变化时）等。在基于液滴的scRNA-seq中，这样的情况是可以被矫正的，由于存在大量的空载液滴，因此这些可以被用来建模出环境污染的基因表达谱。soupx，Pragmatic 就是利用这种方法开发出来的软件。

由于“足够好”的数据的定义难以一开始就给出，所以最终需要看下游分析的表现，比如cluster注释，来看数据质量如何，所以有时候可能是要做多次数据质控的。从允许的QC阈值开始并研究这些阈值的影响，然后再执行更严格的QC往往是有益的。这样的方法尤其与那种含有异质性的cell群体相关，其中的cell类型和状态，可能会被错误的认为是低质量的outlier。而对于低质量的数据集，严格的质控阈值是必须的。数据集的质量，可以被experimental QC matrix确定。在这种迭代 QC 优化中，应该注意data peeking。不应调整 QC 阈值以改善统计检验的结果。相反，可根据数据集可视化和聚类中的 QC 协变量分布评价 QC。

Normalization

1）推荐使用scan用于非全长转录组的数据归一化。也可以选择使用scone，特别是基于芯片/小孔原理的平台数据产生的数据。全长转录组可以矫正基因长度，工具选择的范围很广。2）没有共识认为一定要把数据缩放到均值为0，具有单位方差的情况（z-score）。此文章更倾向不去缩放基因表达。3）标准化的数据需要形成log(x+1)的形式，便于后续假设数据为正太分布的方法的使用。

最常用的标准化方法是利用read count 做归一化，也就是CPM【count per million normalization：原始的表达量除以该样本表达量的总和(此方法的假设前提是：最开始在所有细胞的mRNA量是相同的，count depth不同是因为取样产生的)，在乘以一百万就得到了CPM值 。从公式可以看出， CPM其实就是相对丰度，只不过考虑到测序的reads总量很多，所以总的reads数目以百万为单位。】这种方法使用与每个细胞的count depth成比例的叫做size factor的值来做归一化。对CPM做一个简单的扩展应用：在计算size factor的时候，去除那些在任何细胞中都占总数5%以上的基因。这样就能考虑到少数高表达基因的分子计数的差异性。

Scran的方法，使能分析更多的细胞异质性，此方法是通过在细胞被pooled的情况下（排除技术dropouts），利用在基因间的线性回归模型，去预估size factor。该方法将在不同细胞内差异表达的基因的variability限制在50%以下。

某一个归一化的方法是不能适应所有类型的单细胞数据的。非线性的模型，在plate based（芯片，孔原理）的平台上更具有优势，其能更好的消除批次效应。全长转录组的数据可能受益于考虑到基因长度的归一化方法，但是3’测序的数据却不能。全长转录组测序的归一化方法一般是用TPM进行归一化。

细胞内数量归一化（按理解这里其实是标准化，z score标准化）使细胞间可比，基因数量归一化使基因数量标度化以及在基因间可比。基因数量归一化使其形成均值为0，方差为单位方差的数据集（z score）。

基因间的归一化是否应该做目前还没有共识，比如Seurat是做了基因的scaling的，Slingshot是没有做的。为了保留更多的数据中的生物学信息，此文章的流程是选择没有做基因间scaling的指导流程。

归一化之后，一般要进行log转化。After normalization, data matrices are typically log(x+1)-trans- formed. This transformation has three important effects：1）数据能体现出来foldchange;2)缓和（但是却没有消除）方差均值比；3）更适用于下游分析的工具（假设数据正太分布）

Data correction and integration

归一化试图消除计数采样的影响，但是数据矫正是为了消除一些技术和生物学上的差异，比如批次效应、脱落或细胞周期效应。但是这些变量并不总是需要被矫正。是否需要矫正，取决于预期的下游分析。并且建议分别考虑技术与生物学因素，因为这两种矫正是用于不同的目的，以及存在独特的考虑点/挑战。

1）去除生物因素协变量只是为了轨迹推断，或者其他感兴趣的生物过程没有被去掉的生物协变量所覆盖。2）最好同时考虑去除生物学和技术协变量，而不是分开考虑。3）plate based数据预处理需要回归技术，一般通过非线性的归一化或者向下采样的方法。4）当批次间的cell类型和状态成分一致时，推荐使用Combat做批次效应矫正。5）不同实验的数据的整合，和同一实验的批次效应矫正，应该使用不同的方法（因其需要解决的问题不同）。数据整合的工具，可能会过度纠正简单的批次效应。6）需要对只出现于表达恢复/填补后的显著信号保持怀疑。最好在探索性分析的时候不要使用表达恢复/填补的步骤。

Regressing out bilogical effects

最常见的生物学因素矫正是去除细胞周期的影响，scanpy或者seurat，或者一些其他的复杂方法，比如scLVM，可以使用。生物学影响的数据矫正之前需要明确：1）并不是总是对数据整合有帮助。比如矫正细胞周期，对于推断发育轨迹有帮助，但是其本身也有生物学意义，比如增生细胞群可以通过细胞周期分数鉴定。2）生物信号必须在背景/环境中理解。3）有人认为，细胞大小的变化解释了通常归因于细胞周期的转录组效应。

Regressing out technical effects

用回归的方法消除生物因素的协变量的回归模型，也同样适用于消除技术协变量。单细胞数据中最突出的技术协变量为count depth和批次差别。削弱count depth的影响有助于轨迹推断类的算法的表现，这种算法一般依赖于细胞之间的过渡进行推算（应该是表达差异从少到多的路径推算吧？）。当涉及多个协变量的时候，应该在一步处理中同时考虑多个协变量，以考虑变量之间的依赖性。其他基于回归的消除count effect的方法是用更严格的方法，比如downsampling或者非线性回归。这些在基于plate平台的单细胞数据中尤为适用，因为细胞间的计数深度差异过大会掩盖细胞间的异质性。

Batch effects and data integration

细胞用不同的分组处理时，可能发生批次效应。这些组可以是不同芯片上的细胞、不同测序泳道中的细胞或不同时间收获的细胞。即不同平台，不同时间，都可能造成同样的细胞产生批次效应。仅对来源于统一实验室的不同实验、或者不同实验室的数据，叫做数据整合（同一实验的就叫做批次效应矫正）。批次效应一般用线性方法矫正，而数据整合一般用非线性的方法进行。ComBat，在批次效应梳理表现就良好，有基因表达的线性模型组成，在数据的平均值和方差中考虑了批次贡献。不考虑计算方法的话，最好的批次效应的矫正方法好还是通过巧妙的实验设计预先消除及避免影响。通过合并实验间的条件和细胞（使所有条件一样呗），可以避免批次效应。使用例如细胞标记，或者遗传变异，可能能够区分被合并到一起的细胞。与批次校正相比，数据集成方法面临的另一个挑战是整合不同的数据集。估计批效应时，ComBat 使用一批中的所有细胞来拟合批次参数。这种方法将混淆批次效应与细胞类型或数据集之间不相同状态之间的生物学差异。为克服该问题，开发了经典相关性分析CCA、相互最近邻 MNN、Scanorama、RISC、scGen、LIGER、BBKNN和 Harmony等数据整合方法。数据整合的方法，也可以应用于一些简单的批次效应的矫正，但是要注意非线性模型自由度增高后的过拟合。

Expression recovery

表达恢复/去噪/填补。单细胞数据里面比较突出的一个噪音就是dropout（没有测到数据信息），推断什么是dropout并用合适的数据去替换它，是最新的几个工具在干的事情。表达恢复有助于改善基因间相关性的估计。并且这一步可以被整合到其他步骤，比如归一化，或者批次效应矫正或者scVI工具的下游分析中。虽然大多数数据校正方法以归一化数据作为输入，但一些表达恢复方法是基于预期负二项分布噪声，因此在原始计数数据上运行。没有方法是完美的，方法都可能导致数据过矫正，或者矫正不足，而且也有报道表达恢复导致了假的相关性信号。

Featuture selection

特征选择的目的是保留最具信息的变量（最少的变量代表足够的信息那种），通常使用高度可变基因HVGs，一般1000-5000个，用于下游分析（否则原始数据维度太高计算压力很大）。最常用的HVG的选择已经被整合在了Scanpy和Seurat中。HVG的选择方法为：按基因的均值表达对基因进行分类，并选择方差/均值比最高的基因作为每个分类中的HVG。并且HVG需要在技术矫正后选出，以免被选出的基因为批次效应导致。

Dimensionality reduction

特征选择之后就是数据降维，数据降维可以将数据鞠正转化到一个低维空间，目的是在尽可能少的维度捕获数据结构。这种方法的原理是scRNA-seq数据本来就是低维的。

1）特征选择的时候，推荐根据数据的复杂性，选用1000-5000个高可变基因。2）使用基因表达的均值和方差的特征选择的方法不能在数据被标准化为均为为0具有标准方法时，或者模型拟合的残差被归一化表达值时使用，因此必须要考虑在选择HVG走之前做什么预处理是合适的。3）数据降维与可视化的方法应该分别考虑。两种常用的降维方法为PCA和diffusion map。PCA是线性的，diffusion map是非线性的。

Visualization

1）推荐使用UMAP来做探索性的可视化。使用PCA来进行一般的数据概括。在进行轨迹推断时使用diffusion maps作为PCA的替代选项。2）具有UMAP的PAGA是可视化特别复杂的数据是的合适替代方案。

可视化一般使用非线性的降维方法，最常用的是t-SNE降维以全局结构为代价来获取局部相似性。因此，这些可视化可能夸大细胞群体之间的差异，并忽略这些群体之间的潜在联系T-SNE的替代方案是umap和spring，UMAP在扩展至大量细胞时的速度和能力都很好，在没有特殊生物学问题的情况下，最好选择UMAP。

Downstream analysis

下游分析是提取生物学的insight以及描述潜在的生物学系统。比如具有相似基因表达谱的细胞群代表细胞类型簇；相似细胞之间基因表达的微小变化表示连续（分化）轨迹；或具有相关表达谱的基因表明其共同调节作用。

Cluster analysis

Cluster

细胞分类是靠识别那些具有相似的基因表达的细胞，相似性是靠降维后的距离矩阵来计算（最常用的就是欧几里得距离）。主要两类算法：聚类算法和社区检测算法。1）建议在单细胞KNN图上通过Louvain社区检测进行聚类。2）聚类不需要在单个分辨率下执行。特定的cell群是关注数据集中更详细的子结构的有效方法。

Cluster annotation

只有被成功注释了，cluster才算有效，所以需要群落的细胞类型注释。需要注意的是1）不要用marker gene的p值来证明细胞身份，特别是marker gene不能用于注释时。2）注意相同的细胞检测出来的marker gene可能由于不同的数据集，或者细胞集的细胞类型与状态不同导致。3）如果存在参考数据，推荐使用自动注释联合手动注释。

Compositional analysis

组分分析是为了弄清楚细胞组成比例，这种比例可能会由于疾病发生改变。比如沙门菌感染可以增加小鼠肠上皮中肠上皮细胞的比例。计算整个群体中，分类后的每种细胞的占比。但样本间的细胞群体的比例变化是相互依赖的，某一个细胞的比例变化了，其他的都会跟着变化。为了分析组分的变化情况，单细胞数据需要足够的细胞数和样本数，足够的细胞数用于获得细胞比例，足够的样本数用于评估群体组成变化的背景（有背景才能看出来变化）。

Trajectory analysis

发展成为能观测到的异质性的生物学过程是一个连续性的过程。为了观测到细胞间的转变/变化，分支分化的过程，或者渐进的不同步的变化，我们需要能动态分析基因表达的模型。这类方法就叫轨迹推断（TI）。

Trajectory inference

1）推荐使用saelens的文章作为TI分析的参考。2）轨迹推断并不一定代表了生物过程，为了解释轨迹需要进一步的证据（比如扰动实验，推断调控基因动力学，或者从RNA速度获得支持？）。模型的范围从简单的线性或分叉轨迹，到复杂的图形、树或多分支轨迹都有。在最近的 TI 方法综合比较中，得出的结论是，没有一种方法的性能，能在所有类型的轨迹分析中表现都最佳。

Gene expression dynamics

一种支持推断轨迹不是拟合转录噪声结果的方法是在基因水平分析轨迹。轨迹相关的基因可能包括了调节模型过程的基因，调控基因可以帮助我们理解生物学过程为什么被触发，并且这些调节基因可能作为一些药物靶标。以前轨迹推断的方法主要是顺着轨迹利用cluster细胞的的差异表达来进行测试推断，现在主要是利用回归的方法检测在拟时轨迹过程中基因表达的变化。为了使表达沿该协变量平滑变化，通过拟合曲线或通过额外的局部回归步骤（例如 loess）平滑假时间。拟时分析过程中的差异表达基因会有轨迹推断算法的混淆在里面，同样的cluster间的差异表达基因也有聚类算法的混淆在里面，因此这些过程中获得的P值不应该被认为就是显著的（必须要其他证据？）。由于目前可用的基因时间动力学分析的软件很少，因此目前研究基因时间动态的最佳方法还不能确定。以后高斯模型可能会是研究时间动态的一个自然模型，此外，研究一个调控模块，相对于研究一个基因来说，应该会有跟高的信噪比以及能更加促进生物学的解释。

Metastable states

Cell-level analysis unification

细胞层面的数据整合：聚类和轨迹推断代表了两种不同的view of single cell 数据，这两种数据可以被整合，常用方法为PAGA，可以利用细胞间交互/相互作用的统计模型，在比预期更相似的cluster 节点放置一个edge。以此来能应付被断开的拓扑异构结构以及包括细胞周期的复杂的图表。可以可视化整个数据集的拓扑异构结构，以便进行探索性分析。

Gene-level analysis

Differential expression testing

1）DE分析需要在measure data层面考虑技术协变量进行，不能使用矫正后的数据。2）不能依赖DE工具来矫正具有混杂协变量的数据。模型的规范性需要被保证以获得全帙矩阵。3）推荐使用MAST（更快，综合效果也可）或者limma进行DE分析。

利用异质性表达来理解基因表达的背景（差异原因）。虽然单细胞有一些技术噪音，比如dropout cell与cell间的高可变性，也有一些专门针对单细胞这些问题开发的工具，但是传统的转录组的DE分析的一些工具，也被证明在单细胞上表型良好，特别是在测试中引入基因权重时。根据测试，表现最好的是结合ZINB-wave权重时，DESseq2和EdgeR软件。

Gene set analysis

基因水平差异分析往往产生大量的候选基因，难以解释。因此可以利用共有的特征进行整合/分组这些基因到基因集，并集检验这些特征是不是在候选基因集中过度表现出来。整合DE结果，常利用一些常见的生物学过程，这些生物过程一般存储在一些数据库中，比如MSigDB,GO,KEGG,Reactome等。单细胞领域最近建立了一个利用配对基因来进行配体受体分析的领域。利用细胞簇间的交互/相互作用，是利用受体及配对的同源配体表达推断。配体受体信息可以从CellPhoneDB获得，并可以用于利用统计模型整合不同细胞簇间的高表达基因。

Gene regulatory networks

基因并不是单独组总有的，解释基因间、与小分子间的相互作用，是基因调控网络方法的主要目的，这种方法一般基于基因的共表达，互斥表达，或者回归模型进行的。目前的GRN分析的工具表现都还不太好，因此建议谨慎使用以便获得有价值的结果。需要警惕推断出的生物关系的不确定性，因调控关系富集的基因模块可能比单个基因更可靠。

Analysis platforms

Scater(R), Seurat(R), Scanpy(python)…

Tools

https://hemberg-lab.github.io/scRNA.seq.course

https://github.com/drisso/bioc2016singlecell

Scanpy: https://scanpy.readthedocs.io/en/stable/tutorials.html

Seurat: https://satijalab.org/seurat/get_started.html

https://scrnaseq-course.cog.sanger.ac.uk/website/index.html