综述：高维单细胞RNA测序数据分析工具（上）

生信技能树jimmy

发布于 2020-05-19 23:52:33

2K0

发布于 2020-05-19 23:52:33

文章被收录于专栏：单细胞天地单细胞天地

呐，等你关注都等出蜘蛛网了~

当你的才华还撑不起你的野心时，请潜下心来，脚踏实地，跟着我们慢慢进步。不知不觉在单细胞转录组领域做知识分析也快两年了，通过文献速递这个栏目很幸运聚集了一些小伙伴携手共进，一起成长。

文献速递栏目通过简短介绍，扩充知识面，每天关注，希望你也能有所收获！

单细胞转录组分析综述

单细胞入门-读一篇scRNA-seq综述

文章信息

Author：Yan Wu and Kun Zhang Date：2020 Mar 27 Journal：Nat Rev Nephrol（IF：19.684）

摘要

在单细胞水平上对转录组进行分析的高通量技术发展取得的突破，有助于生物学家对细胞群体的异质性，疾病状态和发育谱系的理解。然而，对这些单细胞RNA测序(scRNA-seq)技术产生的大量数据进行分析和解释有一些挑战。此外，scRNA-seq数据集通常包含了如由不完全的RNA捕获、PCR扩增偏倚和/或来源于特定病人或样本的批次效应引入的技术噪声。如果不加以处理，这种技术噪音会使数据的分析和解释产生偏差。为了应对这些挑战，一套计算工具应运而生被开发来处理、分析和可视化scRNA-seq数据集。尽管任何给定的scRNA-seq分析的具体步骤可能会根据存在的生物学问题而有所不同，但在大多数分析中一般使用的是一个核心工作流。通常，原始测序reads会处理为基因表达矩阵，然后对其进行标准化和缩放(scale)以去除技术噪音。接下来，根据基因表达模式的相似性对细胞进行分组，概括为两维或三维进行可视化。这些数据可以进一步分析，在感兴趣样本中提供一个细胞类型或发展轨迹的深入观点。

总流程：

1 Data pre-processing

从scRNA-seq平台获得的原始数据必须先经过几个预处理步骤，然后才能用于评估基因表达的生物学相关变化。这些预处理步骤将原始数据转换成更有用的格式，并解决与样本质量、广泛的基因表达水平和差异相关的问题。此外，如果要同时分析多个数据集，这些步骤可以减少技术批处理效果的影响。

1.1 Generating a gene expression matrix

在cDNA组成的scRNA-seq实验中首先生成初始输出FASTQ文件，每个read包含一个RNA序列，一个细胞条形码，用来识别产生read的细胞，一个独特的分子index(UMI)用来识别确切的mRNA分子。scRNA-seq分析的第一步是将这些read数据处理为一个count矩阵，该矩阵汇总了在数据集中每个细胞中检测到的每个基因的count数。计数矩阵作为剩余分析步骤的输入，也是一种有效的存储和共享基因表达信息的方法。创建count矩阵通常需要将每次读取的cDNA序列对参考基因组进行比对，以识别读取来自特定基因的read，然后通过细胞条形码将每次读取分配到原来的细胞。

如下图，展示了从测序read比对参考基因组根据barcode和UMI产生Count matrix的过程：

scRNA-seq利用PCR技术对cDNA分子进行指数级扩增，而UMIs技术可以帮助用户识别并消除在扩增过程中可能产生的重复序列，从而降低技术噪音。值得注意的是，UMI中的测序错误可以人为地提高基因表达，因为应该被消除的重复序列被视为了不同的分子。相反，不同的分子可能会被错误地标记为具有相同的UMI序列，从而被视为同一个分子。

对于大多数测序技术，背景RNA污染和测序错误会导致大量的cell barcode，它们的reads数很低，但与真正的细胞并不对应。可以通过设置最小读取read数或细胞barcode的UMI阈值来检测和删除这些空条形码。更复杂的方法，如dropEst也可用。

有几种工具可用于read处理(表1)，包括CellRanger，它与10X genomics Chromium scRNA seq平台配套使用。CellRanger处理cDNA read，运行序列比对，通过UMIs处理重复read和进行QC然后输出count matrix。CellRanger还可以执行二级分析，如聚类(即根据细胞基因表达模式的相似性对细胞进行分组)和可视化(稍后将详细讨论)，尽管使用的是相当基本的工作流。但是，CellRanger可能相当慢且占用大量内存，使用最大30 GB的RAM和花22小时处理7.84亿个read(相当于~50,000个细胞,每个细胞15,000 reads的深度)。然而，CellRanger与另一款10X基因组学软件Loupe Cell Browser的集成，为非专业用户提供了一个交互式浏览器，可用于可视化聚类结果和标记基因的表达。

表1

在过去的几年里，研究人员开发出了scRNA seq方法，通过组合检索，可以在一次实验中分析数十万到数百万个细胞。这些方法包括基于分裂池连接的转录组测序和单细胞组合索引RNA seq。考虑到这些技术上的进步，并考虑CellRanger所需的内存和处理时间，还开发了将cDNA读入单细胞基因表达count的替代计算流程。例如The dropEst流程，比CellRanger更短的运行时长和更低的内存使用情况，并通过纠正细胞条形码和UMIs中的测序错误来提供更准确的基因表达估计，DropEst还通过使用机器学习模型识别空条形码来改进数据恢复，从而能够恢复大小小于平均大小的细胞类型，以及可能被排除在分析之外的RNA含量较低的细胞类型。UMI-Tools 是另一种纠正细胞barcode和UMIs中测序错误的流程，以提供更准确的基因表达定量。

CellRanger流程中最慢的步骤之一是将cDNA read比对到参考基因组。The Kallisto pseudoaligner，使用BUStools作为存储和操作scRNA-seq数据的方法，是一种高效的替代方法。因为它为每个read创建了一个兼容的转录本列表(伪比对)，而不是将单个read调整到转录组中的一个精确位置(比对)。当处理50,000个cells时，合并后的Kallisto-BUStools方法比CellRanger 快51倍，使用了最大~ 12gb的RAM。然而，Kallisto-BUStools不会删除空的cell barcode。STARSolo和Alevin分别是两种比对和伪比对方法的扩展，也可以用来处理scRNA-seq数据。STARSolo和Alevin的运行时间都明显快于CellRanger，但STARSolo的最大内存使用量更高。

综上所述，scRNA-seq分析的第一步是将raw reads处理成count matrix。对于10X genomics scRNA-seq平台的用户，CellRanger为这一处理提供了一种方便的方法，尽管速度慢且需要大量内存。CellRanger还运行基本的聚类和标记基因分析，可以用Loupe细胞浏览器可视化。DropEst、Kallisto-BUStools、UMI-Tools、STARSolo和Alevin都是可选的read处理方法，它们对运行时长和内存进行了实质的改进，允许用户处理他们的scRNA-seq数据而不必在计算基础设施上投入太多。此外，与CellRanger相比，DropEst、UMITools和Kallisto-BUS Tools提供的对UMI和cell barcode错误的增强校正可以改善基因表达估计。

1.2 Quality control and doublet detection

所有的scRNA-seq方法都会产生技术偏差和噪声—一些基本的QC在下游分析之前就解决了这些问题。例如，用于单细胞分离和测序的协议可能会引起细胞应激并导致细胞死亡，这将影响基因表达，并可能导致下游分析中出现假的死亡细胞簇。过滤掉cDNA read数较低或UMI数较低的细胞，以及UMIs总数量中线粒体read数较多的细胞(也称为线粒体分数)，可以帮助去除死亡细胞。与胞质RNA不同，线粒体RNA的存在提示细胞死亡。read数或线粒体read数比例的适当阈值取决于数据集中存在的细胞类型和使用的scRNA-seq方法。设置检测到基因的最小细胞数的阈值，也可以帮助排除那些只在少数细胞中表达，而且不太可能提供信息的基因。但是，用户应该确保这个阈值不是太高，否则在下游分析中可能会遗漏罕见的细胞类型。

另外一篇文章的过滤标准（Current best practices in single-cell RNA-seq analysis: a tutorial）： Cell QC is commonly performed based on three QC covariates:

the number of counts per barcode (count depth)：每个细胞的测序深度，gene count数总和
the number of genes per barcode：每个细胞内表达基因的数目
the fraction of counts from mitochondrial genes per barcode：每个细胞中线粒体基因的比例

对于大多数scRNAseq方法，doublets的产生，即当两个或两个以上的细胞被分配到相同的细胞条形码，会在下游分析中产生假的cluster，因为合并两种不同细胞类型的基因表达模式可能会产生一种独特的表达特征，这在任何真正的细胞类型中都找不到。然而，手动的从真实的cluster中区分出doublet cluster会有很大的挑战，特别是对有多种细胞类型的大数据集。识别doublets的一种常见策略是通过组合数据集中不同cluster中的细胞来生成模拟doublets，并评估哪些细胞具有与模拟doublets细胞相似的表达谱。然而，这种策略只有在数据集包含离散的细胞类型而不是连续的细胞轨迹时才可行。

数据集之间的QC阈值可能不同，一些探索性的数据分析，如每个细胞或基因UMIs分布的直方图，可以帮助设置每个数据集的阈值。在某些情况下，例如当下游分析中出现一组假的死亡或正在死亡的细胞时，在运行整个分析流程后修改这些阈值并重复分析后也会有所帮助（下图）。Seurat和SCANPY是scRNA-seq分析流程包，包括计算QC指标的功能，如每个细胞表达的基因的比例、线粒体比例和总count数；用户确定过滤数据集中的基因和细胞的阈值。Scater也提供了一套用于计算关键QC指标的工具。

1.3 Data normalization

在每个细胞中捕获的RNA的比例可能因逆转录效率、引物捕获效率以及与collapsing的UMIs相关的错误等因素而有所不同。因此，每个细胞的UMIs或reads的总数差异可能是由技术因素而不是生物学变异造成的（下图）。

如果没有标准化，总的UMIs或reads的技术差异可能会主导下游分析。例如，总UMIs或reads数量相似的细胞聚集在一起，而不是基因表达模式相似的细胞。因此，normalization对于揭示数据集的真正生物异质性至关重要。大多数归一化方法试图估计每个cell的偏差(也称为细胞size因子)。所有细胞的UMI或read计数可以通过除以大小因子来标准化，从而可以比较不同细胞间的基因表达水平。总count数归一化是一种简单的归一化策略，其中大小因子很由每个cell中的UMIs或read总数组成。然而，总count的归一化可能由高表达基因主导，当存在很强的细胞类型特异性基因表达时，即细胞或者组织类型差异非常大时，会导致size因子的偏倚估计。此外，一些细胞类型更大，比其他细胞有更多的RNA分子，这样的一种生物学因素，当简单地用UMIs或reads的数量除以总counts数时，它就被掩盖了。

scran包在估计大小因子之前将具有相似表达模式的细胞汇集在一起，因此解决了由于细胞类型特异性基因表达或UMI计数而导致的标准化问题。但是，使用相同大小因子对高表达和低表达的基因进行标准化会导致低表达的基因如转录因子的过校正，高表达的基因如管家基因的低校正。SCnorm通过将与总UMI或read count数有相似依赖关系的基因汇集起来，并在每个池中计算大小因子来解决这个问题。sctransform(在Seurat包中实现)使用概率模型来计算总UMI或每个基因read计数的效应，这也使它能够稳定基因变异(稍后将详细讨论)，并识别过分散的基因（over-dispersed genes ）。

总的来说，某种类型的归一化对scRNA-seq分析是至关重要的，尽管总count归一化成功地缓解了技术上的偏差，但它在一定程度上掩盖了真正的生物异质性。使用专门的标准化方法，如SCnorm和sctransform，可以解锁数据集中额外的异质性。

1.4 Variance stabilization

基因表达水平可以有很大的差异，基因的平均表达(或量)与基因的变异密切相关，这种效应被称为均值-方差关系（如下图）。

方差稳定可以调整数据，消除基因表达量对基因方差的影响。这一步骤确保下游分析集中于生物学上最相关的基因(即数据集中在特定细胞类型中表达的基因)，而不是简单地集中于表达量最高的基因。例如，方差稳定可能通过使具有低平均表达水平的基因(如转录因子)，有助于分析，从而促进一个发育前体细胞的两个亚群的分离，否则这两个亚群可能合并。尽管这些基因在细胞内的总体表达水平较低，但它们在揭示细胞命运方面可能很重要。

一种简单的稳定变异的方法是归一化的count数进行对数变换，它可以减少高表达和低表达基因之间的差异。可以用来消除基因平均表达对基因方差影响的流程包括Seurat、Pagoda2和SCANPY。它明确地符合均值-方差关系并应用了一个比例因子（scaling factor）。ZINB-Wave、singlecell variational inference (scVI)和deep count autoencoder (DCA)使用不同方法的替代方法(负二项分布)来模拟单细胞count数据。

总的来说，虽然方差稳定对于scRNA-seq分析并不是严格必需的，但是调整数据集以适应平均基因表达的大范围变化可以增强生物学相关基因对下游分析的贡献。这种方法消除了在所有细胞中都有丰富的表达但表达水平相似的基因的影响，例如管家基因，这种基因同时对细胞异质性的研究没有帮助。在方差稳定之后，识别和选择高度可变的基因可以提高下游分析中细胞类型的分辨率，特别是在被分析的细胞类型相当相似的情况下。这一可选的处理步骤包括在调整平均基因表达差异后选择残差最大的基因。

1.5 Batch effects and data integration

对来自不同的患者或样本，或从不同的实验中获得的样本，使用不同的技术产生的多个scRNA seq数据集的联合分析，增加了分析的细胞总数。该方法可以提高细胞亚型的分辨率和罕见细胞表型的检测，还可以对患者、样本或技术进行直接比较。然而，由于基因表达的技术差异可以掩盖相关的生物学现象，即批次效应，因此这种类型的分析常常具有挑战性。批次内的变异通常是由于细胞类型和生物学相关因素的差异造成的，而批次间的变异也可能是技术因素造成的。这些批次效应可能来自于患者、样本或影响RNA捕获效率或细胞活力的协议(包括操作员驱动的变化)的可变性。批次效应的强度取决于数据集的类型，在进行分析之前很难预测。

消除技术批次效应的一个简单方法是，基本假设每个批次必须在所有细胞中具有相同的平均基因表达，并使用回归模型消除批次之间的差异。虽然这种批次校正方法适用于bulk RNAseq数据，但当不同批次的数据不一致时，它可能会过度校正。具体来说，如果细胞类型的比例在批次之间不同，这种类型的粗批校正可能会对细胞类型的解析能力产生影响。例如，如果一个肾脏样本比另一个样本含有更多的近端肾小管细胞的收集管，那么由于细胞类型组成的不同，每个样本中细胞间的平均基因表达是不同的。应用批次校正，强制两个样本细胞间的平均基因表达量相同，降低了收集导管和小管间基因表达量的差异，降低了分辨这些细胞类型的能力。

专门为scRNA-seq数据整合定制的方法能够在消除批次效应的同时保留批次之间的细胞类型比例的差异。这些方法中的大多数依赖于在不同批次中找到对应于相同细胞类型或状态的细胞对的概念（如下图）。一旦这些配对(也称为互斥的最近邻居（mutual nearest neighbours，MNNs))，被确定，MNNs之间任何剩余的基因表达差异都被认为是由于批次效应造成的，然后可以被纠正。MNN方法的一个优点是，如果一个细胞类型或状态对于特定批次是唯一的，则不会将其标识为MNN，从而保留每个批次的独特生物学特性。例如，当使用MNN方法整合来自两只小鼠的肾脏scRNA-seq数据时，一只是野生型对照组，另一只是缺乏足细胞的基因敲除，即来自对照小鼠的足细胞在集成后将保持在一个单独的cluster中。值得注意的是，如果敲除导致所有足细胞的基因表达发生了一致的变化，而不是发生了足细胞的丢失，那么这种变化在整合后可能会消失，因为它与批次效应难以区分。