1 文章信息题目：Current best practices in single-cell RNA-seq analysis: a tutorial 发表日期：2019年6月19日 杂志：Mol Syst Biol 文章在：https://www./doi/10.15252/msb.20188746 DOI：https:///10.15252/msb.20188746  图1 2 摘要单细胞领域日新月异，大量的工具被开发出来，但很难去判断是否好用，而且如何组建一个分析流程是一个难点。本文将详细介绍单细胞转录组数据分析的步骤，包括预处理（质控、归一化标准化、数据矫正、挑选基因、降维）以及细胞和基因层面的下游分析。并且作者将整个流程应用在了一个公共数据集作为展示（详细说明在：https://www.github.com/theislab/single-cell-tutorial），目的是帮助新入坑用户建立一个知识体系，已入坑用户更新知识体系。 3 前言
现在已经可以利用scRNA研究斑马鱼、青蛙、涡虫的细胞异质性(Briggs et al, 2018; Plass et al, 2018; Wagner et al, 2018) ，重新理解以前的细胞群体，但这个领域面临的一个问题就是没有成熟的标准化流程。标准化之路的困难有：大量分析方法和工具的诞生（截止2019.3.7 已经有385种工具）、爆炸式增长的数据量(Angerer et al, 2017; Zappia et al, 2018)。另外根据不同研究目的，各种分支也突显，例如在细胞分化过程中预测细胞命运(La Manno et al, 2018)。在我们眼界大开的同时，分析流程标准化就变得更加困难。 在未来分析流程标准化之路上，困难还会存在于技术整合层面。比如现在大量的scRNA工具都是用R和Python写的，跨平台分析需求在增长，而对编程语言的喜好也决定了工具的选择。很多好用的分析工具将自己限制在用各自的编程语言开发的环境中，例如Seurat、Scater、Scanpy。 接下来，就一起看看作者列出了哪些他认为比较好的软件和流程吧 先上一个scRNA分析总体流程图：  图2 4 预处理和可视化4.1 首先看一下实验过程比较详细的介绍可以看：Ziegenhain et al (2017); Macosko et al (2015); Svensson et al (2017). 原文描述的关键点是：
感觉原文描述的还没有illumina给出的详细，那么就看看illumina的图文并茂版：
 图3  图4 原始测序数据要经过处理得到表达矩阵，注意这里有两种表述方式：molecular counts (count matrices) 【也即是使用UMI的】和 read counts (read matrices)，取决于是否使用UMI。而作者介绍的流程中，默认使用 count matrices，除非readmatrices和 count matrices得到的结果存在差异，才会特别介绍read matrices
原始数据处理工具主要有：CellRanger、indrops、SEQC、zUMIs 它们主要做了这么几件事：
得到的矩阵行是转录本，列是barcodes【这里用barcodes而不是直接叫细胞，是因为不同细胞的reads也可能属于同一个barcode =》如果出现一孔/液滴多细胞（doublet情况），那么barcode在多个细胞都是一样的】当然也会出现有barcode但实际没有细胞的情况（一个孔/液滴没有细胞即droplet，但这个孔/液滴也会赋予barcode）
关于10X实验环节，可以看我之前写的：https://mp.weixin.qq.com/s/0DEybX7GnuDFhfY1uj9t9A  图5 4.2 质控在正式分析之前，先要确定barcode是不是对应真正的细胞（上面已经了解了barcode和细胞的关系），也就是进行Cell QC，主要考虑三个因素（这几个因素也就是现在流程中常用的过滤指标）：
先看图A：其中这个小的直方图就是把count depth小于4000的放大，这里设定了一个阈值1500，也就是一个barcode中至少有1500的表达量 图B：每个细胞中包含的基因数直方图。可以看到横坐标有一个小的峰在400附近，这里设定的阈值是700 图C：依旧是看count depth。从高到低排列count depth值，可以过滤一些空的液滴（empty droplets），看到从”肘部“也就是纵坐标1500左右开始迅速下降 图D：看线粒体比例。如果占比很高并且细胞类型不是线粒体特别丰富的那种（如心肌细胞），可能说明这个细胞本身的基因数不多并且总体表达量也不高  图6 以上三个指标固然重要，但如果只关注其中某一个，也会产生误导作用，所以作者建议看问题一定要全面，并且要把数据和生物学知识结合起来。作者举了个例子：比如线粒体表达量相对较高的细胞也可能参与了呼吸过程。细胞总体表达量低或者基因数量少，也可能是因为当时取的细胞处于静止；细胞表达量很高，也可能因为本身细胞体积就比较大。的确，细胞与细胞之间的总表达量还是存在较大差异的。未来也许QC会提供更多的选择。 除了检查细胞完整度，QC还要进行转录本层面上的检查。原始的count矩阵一般包含超过20000个基因。这里一般要根据在细胞中有表达的数量进行过滤，但这个阈值要根据总体细胞数和预计的分群情况来灵活调整。比如有的细胞类型本身就数量比较少（也许就50个），那么如果我们要设定”在少于50个细胞中有表达的基因“这种条件，那么可能会丢失那些总共就50个细胞中的marker基因，最终导致鉴定的细胞亚群会缺失。 质控的目的就是给下游提供更高质量的数据，但一开始谁也不知道这个质量高不高，只能先进行下游分析，看看结果（比如细胞分群结果）再判断。尤其是针对异质性高的细胞群体
小结
4.3 归一化/标准化
表达矩阵中的每个count值都表示成功的细胞捕获、成功的反转录、成功的测序。但即使是相同类型的细胞，它们的count depth（也就是每个细胞的全部表达量）也会有变化，变化的来源就在于上面说的那三步。因此在比较两个细胞时，任何差异都可能由于实验测序误差产生，而不是真的生物学差异。归一化就是解决这个问题，它把要比较的两个count值根据各自身处的环境求出一个相对丰度，也就是放在了一个水平上考虑，减少实验测序误差，突出更多的生物学差异。 最常用的归一化方法就是：count depth scaling，也称为counts per million（CPM），这个方法常用于bulk转录组，它会根据每个细胞的总表达量计算一个 size factor ，然后对其中各个基因表达量进行normalize。
单细胞测序中使用的归一化方法由于细胞种类和基因错综复杂，有人就在bulk的基础上进行了改动。例如：Weinreb et al (2018) 先排除了表达量超过总体5%的基因，然后再计算size factor，主要是预防少量极高表达量基因的存在；Scran包有个pooling‐based size factor estimation方法，允许更高的细胞异质性存在；另外Scran包在批次矫正和差异分析环节也比其他归一化方法表现更好(Buttner et al, 2019)。 在单细胞RNA测序领域，目前有三种常用方法：其一是以10x Genomics为代表的微滴(droplet-based)测序；其二是以Namocell为代表的PCR板(plate-based)测序；其三是以BD Rhapsody为代表的微孔(micro-well-based)测序。就测序长度来说，Smart-seq/C1和Smart-seq2基于full length的测序方案，CEL-seq2, Drop-seq, MARS-seq, SCRBseq是基于UMI的测序方案。 不能指望某一种方法适用于所有类型的scRNA数据，(Cole et al, 2019)就发现不同的归一化方法对于不同类型数据集表现不同，使用scone工具可以帮助选择合适的方法。 一般在归一化后，数据都会变成 使用log转换的一个好处就是：让数据更加集中，减少数据的偏斜度，从而近似于许多下游分析工具对数据为正态分布的假设（尽管scRNA数据并不是真正的符合正态分布），比如在差异表达分析和批次矫正环节 小结
4.4 数据矫正与整合数据矫正的对象种技术和生物因素都有，例如：不同批次、捕获失败（dropout）、不同细胞周期。这些在之前的归一化中没有被矫正，但这些差异因素都可能会后面的分析产生影响，它们现在都是导致差异的”嫌疑人“之一。这里要做的就是把这些差异来源去掉（Regressing out 《=》【专门查的词典】 同义词partialling out ：剔除） 4.4.1 首先是生物因素最常见的生物矫正因素就是：转录组中的细胞周期信息。简单一点的方式就像Scanpy和Seurat对细胞周期评分进行简单线性回归；复杂点的方式就像scLVM和f‐scLVM。用来计算细胞周期分数的marker基因可以从文献中获得 (Macosko et al, 2015)。另外，这些方法还能用来去除其他已知的生物因素，例如线粒体基因表达量（可以作为细胞应激的标记）。 需要注意的是：
4.4.2 然后是技术因素最常见的技术矫正因素就是：样本测序深度、批次、噪音。 去除测序深度的影响，可以促进轨迹推断算法的表现，因为它需要在细胞之间找变化的路径，只要放在同一水平才能看到更准确的总体表达高低。 批次的来源可能是：细胞捕获的时期不同、文库制备使用的芯片不同、测序使用的lane不同。由此产生的效应存在于多个层面:一次实验中各个细胞群之间、同一实验室中进行的不同实验之间、或来自不同实验室的数据集之间。这里主要介绍第一种和最后一种情况：
看一下Combat矫正前后的差别：其中颜色表示不同样本  图7 去噪也是矫正的一种类型。单细胞数据的一个特点就是含有许多噪音来源，其中一个就是dropout。一些工具就用来推断dropout，用适当的表达量来替代0，例如：MAGIC、DCA、scVI、SAVER、scImpute。去噪可以提高基因间相关性的估计。这一步可以和归一化、批次矫正及其他下游分析整合起来，例如基于Python的scVI工具。但任何方法都可能导致矫正过度或不足。 4.4.3 小结
4.5 挑选基因、降维、可视化人类的scRNA数据中可能会包含25000个基因，但其中许多基因并非能提供有用信息，还有很多基因表达量直接为0。即使在QC阶段去掉这些表达量为0的基因，一个单细胞数据的基因空间依然会有超过15000个维度（一个基因表示一个维度），因此需要降低维度 4.5.1 首先挑选基因就是挑那些真正”具有情报价值“的基因，也就是会数据变化起作用的基因。因此我们这里会挑选名为HVG的基因，也就是highly variable genes。根据数据集的复杂程度不同，HVGs一般会有1000-5000个（如下图就对不同数据集的HVGs做了个统计）  图8 之前有研究表明，HVGs数量从200到2400，它们降维后的表现差不多（Klein et al (2015)，作者建议先尽量多选一些HVGs。 比较流行的挑选HVGs的方法有Scanpy和Seurat，而且最好是在去除技术因素后挑选，避免因为批次、测序等因素导致错误挑选HVG。当然还有其他挑选的方法，看Yip et al (2018). 4.5.2 接着降维挑出来HVGs后，就是降维了，力求在最少的维度中捕捉到最多的数据特征。 常用的降维方法：A-F分别是：PCA、t-SNE、diffusion maps、UMAP、ForceAtlas2（force‐directed graph）、Variance explained by the first 31 principal components (PCs)。关于单细胞数据的降维方法，详细可以看：Moon et al (2018)  图9 其中两个应用比较广的方法是：PCA（Pearson, 1901）和diffusion maps (Coifman et al, 2005) 【diffusion maps 于2015年在单细胞领域走红 Haghverdi et al (2015) 】
4.5.3 最后可视化可视化一般使用非线性降维的方法。最常用的就是2008年提出的t-SNE（ t‐distributed stochastic neighbour embedding）。t-SNE的一个特性就是关注局部而忽视整体，因此带来的一个影响就是：可视化结果可能夸大了细胞群之间的差异，忽略了这些细胞群之间的潜在联系 另外，使用t-SNE的一大难点就是 除了t-SNE，还有2018年推出的UMAP和SPRING可以用，在缺乏明确的生物学问题时，可以用UMAP作为不错的数据探索。 小结
4.6 「总结」 预处理的各个阶段作者贴心将预处理比作5种类型数据的处理： 原始数据（raw data）、归一化数据（normalized data）、矫正后的数据（corrected data）、挑选后的数据（feature‐selected data）、降维后的数据（dimensionality‐reduced data） 这5个阶段又分成3个层次：
其中每个步骤适时调整，例如单一批次的数据集，就可以跳过矫正批次这一步  图10 5 下游分析之细胞层面下游分析的目的是解释生物问题，例如根据表达量将细胞划分成不同的类型；相似细胞间表达量的微小变化也会体现连续的分化路径；基因表达量之间的相关性可能与基因共表达有关... 下游分析也是有细胞层面和基因层面：
 图11
5.1 细胞分群5.1.1 先是：分群方法
将细胞分群基本就是任何单细胞分析的必经之路。群的划分就是根据细胞中基因表达谱的相似性，表达谱的相似性是由于欧几里得距离量度决定的，而距离量度又是利用的降维的数据。一般有两种方法计算：clustering algorithms、community detection methods
5.1.2 然后是：分群后的注释这个过程主要是基因层面的操作，为每个cluster找marker gene（也就是能代表这个cluster的基因，而这个基因又和已知的细胞类型有关）。任何的分群算法和参数设置都会将一整团细胞分成多个群，但这些群是否真的有意义，就要靠这一步来和生物背景结合起来。 我们希望看到的是存在很多类型的细胞，来说明细胞异质性的问题，但这里关于细胞类型这个定义还是存在争议。首先，细胞类型的划分怎样算是清楚，对于一些人来说，”T cells“这个名称可以叫一个细胞类型，但还有人认为，必须继续深入，像”CD4+ T cells“、”CD8+ T cells“才叫细胞类型；另外，即使是同一种细胞类型的细胞也会有不同的发育状态，因此它们也会显示不同的分群结果。但不管如何，它们都是当时细胞的一种身份（identity）
因此，我们将分群的结果称为不同身份的细胞（cell identities）会比不同类型的细胞（cell types）要好一些【即每个亚群可能并不是真的不同类型细胞，只是显示了此时此刻的细胞身份】 对于不同细胞身份的注释，近年来也随之细胞图谱的研究而加速，例如小鼠脑细胞图谱 (Zeisel et al, 2018) 、人类细胞图谱 (Regev et al, 2017)的发现，产生了许多参考数据库。在缺乏相关背景的情况下，我们可以借用数据库中已发现的细胞marker 基因套入我们的细胞，帮助判断细胞身份。需要注意：通常使用的细胞表面marker基因在细胞身份鉴定方面存在局限性(Tabula Muris Consortium et al, 2018)
图12
利用差异分析，分成两组：某个cluster中的细胞、数据集中其余全部的细胞。然后重点关注这个cluster中上调的基因，因为marker基因一般具有更强的表达作用。差异分析也会使用简单的统计检验，例如Wilcoxon rank‐sum test、t-test，将基因的差异大小排个序，选出排名靠前的基因来作为marker基因
将数据集中选出的marker基因和参考数据集进行比对，统计方法可以是：enrichment tests、the Jaccard index、other overlap statistics 参考数据集可以是网页工具： www.mousebrain.org、 http:///，可以将选出的marker基因在参考数据集中进行可视化，帮助判断这个marker基因是什么细胞身份
细胞分群、分群注释、重分群、重注释...这个循环很耗费时间。自动化注释方法加快了这个过程，例如scmap (Kiselev et al, 2018b) 、Garnett (preprint: Pliner et al, 2019) ，但这样的方法有利有弊。自动化提高了速度，但相比手动注释也降低了灵活性。毕竟自动化工具使用的参考数据集中可能并不包含我们数据中的这样细胞。因此，有自动化工具也不能完全抛弃手动挑选，尤其针对大型数据集中多种多样的细胞。自动化的过程可以先帮我们粗略地给细胞加个标记，如果有需要，我们可以继续手动对这种细胞继续划分子细胞。对于小型数据集或者缺乏参考基因集的，手动注释就足够了。 5.1.3 注意
5.1.4 细胞分群衍生——细胞组成分析（Compositional analysis）就像上面的图12中的C图，显示的是近端(上图)和远端(下图)肠上皮区域的细胞身份组成图（颜色越深细胞密度越大）。研究细胞组成的变化也是一个新方向，例如沙门氏菌感染已被证明会增加小鼠肠上皮细胞的比例 (Haber et al, 2017)。 这个分析既需要足够多的细胞数量来推断各个cluser的占比，又需要足够的样本数量来证明是单纯一个样本得cluster数量这样变还是总体都会这样变。相关的分析工具还没有太多，未来的开发可能会借鉴单细胞质谱流式（mass cytometry）或者是宏基因组分析【单细胞与宏基因组的结合...】 5.2 轨迹分析5.2.1 轨迹推断Trajectory inference轨迹推断就是为了找到不同细胞身份、分化或者生物过程中渐进式非同步的变化，构建出的一个动态模型。它认为单细胞数据实际上就是一个连续过程中的快照（snapshot），这个过程可以通过在细胞空间中寻找最小化相邻细胞间转录变化的路径来重建
图13 2014年Monocle和Wanderlust先推出了轨迹推断，之后诞生的分析方法更加丰富，它们在建模路径的复杂性上有所不同，从简单的linear or bifurcating(分叉) trajectories，到复杂的graphs, trees, or multifurcating(多叉) trajectories。Saelens et al, 2018)进行过轨迹推断方法的比较，结论是没有一种方法对所有类型的轨迹推断有效，应该根据预期轨迹的复杂度来选择。不过，Slingshot在简单轨迹推断中优于其他方法(Street et al, 2018) 。如果期望得到更复杂的轨迹，PAGA值得推荐。轨迹推断是一个不确定的过程，可以用多种方法来进行佐证。
5.2.2 基因表达量的动态变化在拟时序（pseudotime）中变化的基因描述了轨迹，这组与轨迹相关的基因有望包含调控建模过程的基因，可以用来识别潜在的生物过程。 目前很少有专门分析基因表达动态变化的工具。BEAM将Monocle的轨迹推断整合进来，允许检测在轨迹分支过程中相关基因的动态变化。另外还有LineagePulse (https://github.com/YosefLab/LineagePulse)考虑了dropout技术噪音但还在开发中。 下面这样的图在Slingshot的帮助文档就有提及：https:///packages/release/bioc/vignettes/slingshot/inst/doc/vignette.html 【4.1：Identifying temporally expressed genes】 Slingshot基因表达量的动态变化 5.2.3 细胞亚稳态分析 Metastable states
拟时序分析会展示出不同阶段细胞数量的多少。假设细胞以无偏的方式采样，其中轨迹中的稠密区域就表示转录时首选的方案。当把轨迹理解为一条时间线时（例如在发育这个时间线），这些密集的区域可能代表细胞的亚稳态，可以结合拟时间坐标来绘制直方图，找到这些亚稳态【因此看到B图中很多种状态，但C中直方图认为这几个密集的区域才属于亚稳态】 Metastable states 5.2.4 整合分群与轨迹分析
将分群的结果当成节点（node），将轨迹当成节点之间的桥梁（edge），所以将动静数据结合在了一起。利用partition‐based graph abstraction（PAGA）这个工具就能得到类似下面这个图。
整合分群与轨迹分析 6 下游分析之基因层面之前都是对细胞进行分析，但细胞中的基因分析会提供更多的信息。例如差异表达分析、基因集分析和基因调控网络推断，不是表面上研究细胞异质性，而是基于异质性探索基因表达相关的原因 6.1 差异表达分析
这个方法也是常规bulk转录组中经常做的。不过单细胞相比于bulk转录组的一个优势就是：可以深入一个层次，原来bulk只是看一块组织的平均表达量，但现在经过分群后，能得到一块组织中各种各样的亚群，再结合差异分析，对理解异质性问题更有帮助。 虽然都是朝着一个方向前进，但单细胞和bulk转录组的差异分析方法还是不同的。
但最近(Soneson & Robinson, 2018)研究表明，基于大批量的差异分析，bulk分析方法的性能与最好的单细胞分析方法相当。当bulk方法进行改进，加入基因权重分析后，表现要好于单细胞原有工具。例如：bulk差异分析工具DESeq2/EdgeR + ZINB‐wave工具估算的权重。 不过，bulk差异分析工具的性能虽然好，但是计算的效率很难提升。毕竟单细胞数据样本数量越来越多，程序跑的时间长短也成了衡量工具优劣的重要因素。单细胞工具MAST脱颖而出。在单个数据集的小范围比较中，完胜bulk和其他单细胞方法(Vieth et al, 2017)。而且MAST比bulk方法快了10到100倍 (Van den Berge et al, 2018) 。 小结
6.2 基因集分析
例如差异分析我们往往能得到上千基因，为了比较方便解读，一般会把有共同特性的基因归为一组，然后检查我们归类的可靠性 【grouping the genes into sets based on shared characteristics and testing whether these characteristics are overrepresented in the candidate gene list.】 我们一般关注基因在生物过程（biological processes, BP）中的富集，可以使用MSigDB、GO、KEGG pathway、Reactome数据库 另外，单细胞中的一个新进展就是利用成对基因标签进行配体受体分析（ ligand–receptor analysis）
配体-受体成对标签可以从：CellPhoneDB数据库获得，然后用来解释cluster之间高表达基因的联系 例如，利用celltalker 就可以做 Celltalker分析 6.3 基因调控网络 gene regulatory network (GRN)
方法例如：SCONE、PIDC、SCENIC (Single-Cell rEgulatory Network Inference and Clustering)，但发展还不是很完善，推断的调控关系不是很稳定【谨慎使用】 7 分析平台现在开发了很多平台，整合了一套分析流程，有基于R的(McCarthy et al, 2017; Butler et al, 2018) ，python的 (Wolf et al, 2018)，本地的(Patel, 2018; preprint: Scholz et al, 2018) ，网页版带可视化的(Gardeux et al, 2017; Zhu et al, 2017) Seurat是使用最广泛的，Scater在QC和预处理中表现优异；除此以外，基于Python的scanpy也逐渐发展起来，它对于大量细胞的标准化方面表现不错 如果不使用命令行，可视化界面也有，只不过用户只能跑人家已经写好的脚本，操作灵活性不足。这样的平台更多的用处是在可视化探索上，例如Granatum、ASAP。未来 Human Cell Atlas（HCA）会在数据可视化探索上迅速发展： https://www./data-sharing 8 结语8.1 作者的结语作者把流程测试和说明都放在了：https://github.com/theislab/single-cell-tutorial 感兴趣的可以跟着走一遍，比较一下不同的工具。作者希望这一篇能代表单细胞领域目前发展的一个最新动向。他也提到，新方法层出不穷，本文介绍的大量的方法是经过实践比较、验证过的。目前可用的方法不管是运行效率还是易用性可能都不如最新开发的方法，但要注意：新方法在未被大量验证之前都需小心使用。而且新方法一般都是针对单个层面（比如降维、分群、轨迹推断等），大体的分析流程基本固定了。 未来整合深度学习和单细胞多组学是两个重要的发展方向，流程化运行更是趋势。 随着文库制备和测序技术的进步，未来的单细胞平台必将可以处理多种类型数据：DNA甲基化、蛋白丰度等等。 8.2 刘小泽的结语
三天的时间，基本每天都会花半天时间在阅读这篇综述上。从第一眼看到它的文章逻辑，就感觉：嗯是它，没错了！连午觉都不想睡了。 一开始想强迫自己看下去，没想到，越看越精彩。尤其是将整个流程和自己的知识结合起来，就看得比较顺畅。为了更加易读，我在其中加了很多注释，包括之前自己写的一些推文和网上一些好的资源，可以帮助梳理知识点。 最后，希望看完本文对你有帮助🤓！ 欢迎关注我们的公众号～_～ Welcome to our bioinfoplanet! |
|
|
