单细胞转录组数据分析流程

回顾一下：什么是单细胞转录组测序技术？

细胞异质性是生物组织的普遍特征。由于传统的转录组测序(RNA-Seq)技术的测序水平是在个体或群体水平上对数万个细胞进行转录组测序，因此传统转录组测序技术的测序结果就只能检测到个体间或者群体间的转录组差异，而细胞间的转录差异则无法精确地检测到。而单细胞转录组技术则提供了一种在单个细胞水平进行高通量转录组测序的一项新技术，能够有效解决细胞间转录组异质性以细胞群间转录组异质性的难题。

单细胞转录组数据分析的难点主要在于细胞的质量不确定，细胞的数量大，从单细胞测序技术诞生至今，测到的细胞通量越来越高，现在一次单细胞转录组测到的细胞数可达100K~200K[1]。因而，对分析人员的要求也越来越高。

虽然单细胞转录组的分析不容易，但依然是有清晰的流程哒（见下图）:

接下来我们一起看看，每一步都需要做些啥。

01测序原始数据的处理

测序原始数据通常指测序下机得到的fastq文件，需要经过一定的处理，将其中我们需要的信息，如barcode，UMI以及基因的序列等，给提取出来，方便下一步分析。

最初处理原始数据常用的是perl脚本，后来有了更方便的软件或工具。目前我们常用的是fastp、 cutadapt、 trimmomatic等分析工具。这步处理主要是为了去除测序时引入的连续的N、低质量reads、以及建库时引入的接头序列等。

通过这步分析，我们可以得到关注的barcode、UMI以及基因的序列。

02 获得表达矩阵

处理完fastq之后，我们需要从中分析出每个细胞中基因表达的信息，即获得表达矩阵。对于这一步处理，我们常采用的是STAR或者salmon，kallisto等比对工具，将测得的序列片段比对到参考基因组或者转录组。同时根据建库时的barcode白名单对每个真实捕获到的细胞barcode进行比对，分出每个细胞的基因表达矩阵。

表达矩阵示意图

表达矩阵中包含了每个细胞转录组中各个基因表达水平的信息，是我们后续各类分析的基础。

这样的分析之后，我们可以统计得到细胞的个数，各个细胞表达的基因数等信息。同时，通过对这些信息的统计分析，我们还可以判断单细胞测序数据整体的质量，为后面的分析步骤提供依据和参考。

单细胞测序数据质控的指标有很多，这里我们来重点看看3个最为常见的指标。

细胞数 Number of Cells
即捕获到的细胞数，是通过分析与细胞关联的条形码的数目计算出来的。根据这个值，我们可以知道这次单细胞测序捕获了多少细胞。

中值UMI数 Median UMI Counts per Cell
这个指标代表的是每个细胞中被检测到UMI数据的中位数。UMI是目前许多高通量单细胞测序平台用到的一种分子标签，会给细胞中每个被捕获的mRNA分子打上一个独特的标签，用来在分析中校准基因的表达量。通过这个指标，我们可以了解到每个高质量细胞中大概有多少个mRNA分子被捕获到。

中值基因数 Median Genes per Cell
这个指标代表的是每个细胞中被检测到基因数目的中位数。虽然人体一共有约2万个基因，但由于转录水平的不同和测序量的限制，每个细胞中能测到的基因只是这2万个中的一部分——当然，我们希望能测到的基因越多越好。这个指标可以让我们了解到，在这次单细胞测序实验中，每个细胞中大概有多少个基因被测到。

03 细胞过滤

虽然上一步中我们得到了所有细胞中基因表达的信息，但并不是每个细胞中信息的质量都符合我们后续分析的标准，因此，我们需要对细胞进行过滤，以便获得相对完好的细胞。那么，怎样进行过滤呢？

在单细胞测序分析中，过滤的标准往往是某些特定基因的表达量，用来鉴别出质量欠佳的细胞，将其过滤掉。其中最重要的参考标准是基因数以及线粒体基因表达情况。

以下3幅小提琴图，分别展示了基因数，mRNA分子总数、线粒体基因占比这三个常用的过滤指标。

细胞过滤参考的指标

首先可以通过基因数、mRNA分子数、线粒体基因占比三个参数进行质控去除质量差的细胞。

• nFeature_RNA 是每个细胞中检测到的基因数量。
• nCount_RNA 是细胞内检测到的mRNA分子总数。
• percent.mt 是细胞内线粒体基因表达量占所有基因表达量的比例。

• 如果nFeature_RNA 过低，表示该细胞可能已经死亡或将要死亡或者可能是空液滴。

• 如果nFeature_RNA 与 nCount_RNA 数值过高，表示细胞在形成油包水的结构制备过程中，两个或者多个细胞被包裹在一个液滴中。

• 如果线粒体基因占比较高，则说明细胞的质量较差。这是因为线粒体基因会在受损或凋亡细胞表达升高，因而线粒体基因占比较高，表明细胞可能已经受损或者正处于凋亡过程中。

不过，每种细胞或组织类型如何设定线粒体阈值，要依实际情况而定。比如某些细胞的呼吸作用很旺盛，其线粒体基因的比例就会可能很高，而不是因为细胞破裂或者细胞状态不好引起的。而有些细胞本来基因的表达数就很少，比如中性粒细胞。所以这三个参数的设置要根据细胞类型而设置。

04 降维和聚类

拿到过滤后的细胞后，我们就可以进行进一步的分析，了解样本中有哪些类型的细胞，每个细胞分别属于哪种细胞类型，甚至细胞亚型。

要做到这一点，我们首先要知道哪些细胞是属于同一类的，这就需要进行降维和聚类。

所谓降维，就是把多维度的复杂数据用更少的维度展示出来，同时尽量保留原始数据中的主要信息。比如照片和地图，就是对三维物体和真实世界的一种降维展示。

从三维的地球到二维的世界地图，就是一种“降维”

而聚类的概念就比较简单了，顾名思义，就是把相似的类别聚在一起。

单细胞测序分析的降维聚类图，就是将各个细胞的基因表达情况在二维平面上展示出来，并且将基因表达特征近似的细胞聚在一起。

在降维聚类图中，细胞间的距离是由它们表达谱的相似程度决定的。表达谱相似的细胞会聚在一起，被标记为同一种颜色，提示它们可能属于同一种细胞类型，为后续判断细胞类型提供分析基础。

聚类后UMAP可视化结果

05 找到细胞簇的 Maker 基因

对于第四步中发现的每一个细胞簇（cluster，即降维聚类图中聚在一起的一群细胞），我们可以通过分析找到在其中特异表达的cluster marker基因，用于后续的细胞类型注释分析。

在通常情况下，我们会将某一个cluster与其他所有cluster相比的差异基因作为这个cluster的marker基因。当然，如果需要的话，也可以计算两实验组间或者两cluster间的差异基因来作为marker。这些都可以用Seurat软件包内的FindMarkers函数来实现。

06 细胞类型注释

在得到细胞簇以及它们的marker基因后，我们就要对这些细胞簇的细胞类型进行判定，这一步就是细胞类型注释。

细胞类型注释是基于不同细胞类型中特异表达的marker基因来进行的。在第五步中，我们找到了每个细胞簇的marker基因，如果某个细胞簇的marker和某个细胞类型的marker基因相符，就可以被判定为对应的细胞类型。

这一步是单细胞分析中非常重要的环节，有一些细胞自动注释软件可以帮助我们定义细胞类型，比如singleR或者scCATCH。

当然受限于前期实验设计或数据分析的差异，自动注释的结果有时并不能与预期相符，我们还可以通过单细胞公共数据库（比如CellMarker、PangLaoDB、CancerSCEM、SingleCellPortal等）或者已发表文章，来寻找自己感兴趣的单细胞注释参考数据集或已知的细胞类型marker，以提高注释准确度。

比如，对于外周血单个核细胞（PBMC）数据集，我们可以用第五步中的方法计算出每个细胞簇的marker（下表中第二列），然后基于这些marker基因，就可以找到对应的细胞类型（下表中第三列），于是就能轻松地进行细胞类型注释啦！

进行了注释后，我们在降维聚类图上看到的，就不再是以数字编号的细胞簇，而是有名有姓的具体细胞类型：

细胞类型注释结果样例

当我们获得了完整的细胞类型注释后，就可以开始进行下游的深入分析啦，比如不同细胞类型的差异基因、通路富集，也可以进行拟时序分析、细胞通讯分析等等，对样本中各类细胞的功能、状态和相互作用进行更加深入详细的分析。

其他

继续介绍一下转录本定量分析、实验设计、批次效应和混杂因素。🤒

我们先思考几个问题，如下：
Q1: 不同protocol有什么区别，优缺点是什么？
Q2: 在进行scRNA-seq的实验设计时，要考虑哪些问题？
Q3: 与bulk RNA-seq的数据相比，scRNA-seq数据有什么不同？

1. 定量方法

目前我们常见的转录本定量方法有两种，full-length和tag。full-length实现整个转录本的count，而tag的只capture5’或3’端。

1.1 full-length
scRNA-seq的full-length文库构建与bulk RNA-seq相似，如SMART-seq2。从理论上讲，full-length应该可以提供一个均匀的转录本coverage，但有时在coverage上还是有一定的偏差。full-length一大优势就是可以检测到不同剪接体(splice variants)。

1.2 tag
如果使用tag的方法进行scRNA-seq，则只对转录本的一端（3'或5'）进行测序。目前大多数scRNA-seq都是基于tag的，如10x Chromium，

优点：可以与UMI（unique molecular identifiers）结合，提高定量的准确性。
缺点: 由于只限于转录本的一端，无法区分isoforms。

Note! 这个图展示了不同细胞中average coverage的情况，有明显的3' bias。

而且3个细胞群明显离群，可能是RNA降解导致的。

1.3 为什么使用UMI
由于在PCR的过程中，扩增是指数级的，可能会导致扩增不均，从而高估基因的表达量。为了解决这个问题，cell barcodes会标记上一段随机核苷酸序列（UMI），而这个UMI是唯一的。在读取count时，将UMI纳入，从而更准确的计算转录本的丰度。

1.4 选3’ 还是5’ tag
这个可能要根据大家具体的实验目的来进行选择，常用的就是3’的方法。但5'也有其优势，如可以获得有关转录起始位点（TSS）的信息，从而探索不同细胞之间是否存在不同的TSS。

2. 实验设计

**那么多方法怎么选？
首先我们要明确的就是选择不同方法还是要基于你的科学问题，你的研究目的。

低通量的方法与高通量的方法相比具有更高的灵敏度，如10x Chromium。

另一方面，低通量方法很难capture到样本中一些比较稀有的细胞类型，导致细胞群的特征不完整。

scRNA-seq数据的不同之处
测序完成后，每个library代表一个细胞，而不是一群细胞。所以，每个细胞都是独一无二的，在单细胞水平上没有办法进行 “生物学重复”。我们一般需要进行相似性聚类，然后在相似细胞群之间进行比较。

批次效应
批次效应（batch effects）是一定要考虑到的问题，即使用不同的技术对相同的样本进行scRNA-seq，也会有批次效应，可以通过normalise来减少批次效应。

混杂因素
整个scRNA-seq的过程中，应避免实验因素（如治疗、表型或疾病等）、准备样品时间、测序时间等对结果的影响。

举个栗子
假设我们准备对10个病人的control和diseased组织进行scRNA-seq，如果每天只能处理10个样本，最好是每天做5个control和5个diseased的样本，而不是一天准备所有control的样本，另一天准备所有diseased的样本。
另一个需要考虑到的就是样本的可重复性。

当从一个器官收集组织时，最好从器官的不同部位采集多个样本。
由于基因表达可能受昼夜节律（circadian changes）的影响，我们最好也在同一个时间点进行取样。

参考文献

[1] Svensson V, Vento-Tormo R, Teichmann S A. Exponential scaling of single-cell RNA-seq in the past decade[J]. Nature Protocols, 2018, 13(4):599-604.

[2] Malte D L., Fabian J T.. Current best practices in single‐cell RNA‐seq analysis: a tutorial. Molecular Systems Biology. 2019 Jun; 15(6): e8746.

[3] Macosko, E. Z. , Basu, A. , Satija, R. , Nemesh, J. , & Mccarroll, S. A. . Highly parallel genome-wide expression profiling of individual cells using nanoliter droplets. Cell, 2015, 161(5), 1202-1214.

[4] Butler, A. , Hoffman, P. , Smibert, P. , Papalexi, E. , & Satija, R. . Integrating single-cell transcriptomic data across different conditions, technologies, and species. Nature Biotechnology, 2018, 36(5).

[5] Papalexi E, Satija R. Single-cell RNA sequencing to explore immune cell heterogeneity. Nat Rev Immunol. 2018;18(1):35-45.

参考：
https://zhuanlan.zhihu.com/p/532134856
https://blog.csdn.net/m0_72224305/article/details/127148666