利用Bioconductor实现scRNA-seq数据的分析

Posted on | In
Bioconductor提供了完成单细胞转录组分析所需的大量软件包，本综述向我们展示了如何利用Bioconductor实现整个分析流程

Introduction

Bioconductor 是一个基于 R 语言的生物学开源软件包仓库，吸引了来自不同领域，包括基因组学，蛋白质组学，代谢组学，转录组学等的大量开发者和用户。随着 scRNA-seq 技术的广泛应用，其产生的数据集的两大特性需要我们有效地解决：一是样本量规模的大幅增加，另一个则是数据稀疏度的增加。数据复杂度以及稀疏度的增加需要我们在数据获取，管理，分析等方面进行根本性地改变，并开发出新的方法来分析这些数据，而 Bioconductor 恰好提供了最新的分析软件和流程，从而实现：

  1. 高效的数据导入和表示
  2. 统一的数据容器，用于存储和软件包之间的交互
  3. 快速且强健的方法，将原始的单细胞数据转换成可用于后续分析的加工后数据
  4. 交互式数据可视化
  5. 下游分析，注释以及生物学解释

本文后续将对整个分析流程以及用到的软件进行说明。此外，作者还提供了在线电子书，囊括了更详细的说明和更多的相关资源

测序数据预处理

得到测序数据后，通常我们会将 reads 比对到参考基因组，量化，构建表达矩阵。R/Bioconductor 上提供了 scPipe 包，能够利用 Rsubread 包在 R 中实现比对。对于基于液滴的测序方法如 10X Genomics 等的数据，DropletUtils 包能够读入通过 CellRanger 流程产生的 UMI 计数矩阵
此外，一些新的 pseudoaligners 使得比对能够在 PC 上完成，如 Salmon 和 Kallisto。tximport 包可以将这些比对软件的结果导入到 R 会话中，并在 tximeta 包的配合下，构建后续分析所需要的实例

数据结构

Bioconductor 分析流程的一大优势就是使用了统一的数据结构，从而实现模块化以及软件包之间的交互。Bioconductor 利用 S4 面向对象编程风格，对数据的存储和获取进行标准化，并在单个数据容器中纳入丰富的元信息注释，便于分析。下面我们将着重介绍 SingleCellExperiment 类
SingleCellExperiment 类衍生自 SummarizedExperiment 类，因此继承了其适于大规模数据分析的优势，并增加了单细胞数据分析特需的方法和结构。不同信息部分称为 slot，以下是常见的一些 slot：

  • assay slot，对应最核心的数据，即表达矩阵。可以是 read 或 UMI 计数，且可以有一至多个矩阵，比如分别表示原始以及标准化后的计数。通常矩阵的行对应基因（或转录本等 feature），列对应细胞（sample）
  • rowData slot，行和列都会有丰富的注释信息，该 slot 中放置的即为行相关的元信息，比如基因的 ID，GC 含量等，特别地，如果行对应的是基因组特征，那么会有 rowRanges slot 存储基因组坐标信息
  • colData slot，列元信息，即样本层面的信息，可以自行增加列元信息，比如每个细胞中测到的 reads 数量等，适用于质控时的筛选
  • sizeFactors slot，指向列（细胞），包含标准化所需的信息
  • reducedDims slot，包含数据在低维度空间中的表示，可以有多种降维方法，如 PCA，t-SNE，UMAP

另外，最近涌现了许多同时测量细胞中基因，表观遗传以及转录信息的方法，MultiAssayExperiment 包能够整合异质的数据类型，比如 SingleCellExperiment，DelayedArray等，以便进行后续分析

质控和标准化

构建好 SingleCellExperiment 对象后，第一步需要做的就是鉴定、移除、校正低质量的 feature 和 sample，整个步骤包括：过滤低质量的细胞，选择能提供有用信息的特征，细胞和基因层面的标准化以移除细胞和基因特异性偏好性，以及针对已知协变量和隐变量进行调整

细胞和基因质控

对于基于液滴的测序方法，DropletUtils 包可用于鉴定空液滴以及减弱 barcode swapping 带来的影响。scater 包能够自动计算一系列质控指标，其中文库大小（library size，单个细胞中 read 或 UMI 总量）、spike-ins 或线粒体基因上计数所占比例、检测到的基因数等是常用的细胞筛选指标；而基因的平均表达量或检测到的频率可用于移除那些不提供信息的基因。simpleSingleCell 流程中展示了如何利用这些包进行质控

基因表达标准化

细胞和基因水平的偏好可以通过比例因子（scaling factor），也称量化因子（size factor）进行校正，从而使得不同的细胞能够进行比较。SCnorm，scran 和 scater 包均提供了标准化的方法，其中 SCnorm 包能够估算并消除由于测序深度，GC含量以及基因长度等引起的基因水平差异，scran 首先将表达谱相似的细胞聚集在一起，然后计算细胞特异性量化因子用于标准化，scater 包则是基因文库大小计算量化因子。需要注意的是，以上方法只能计算能直接从数据中推导出来的内部因子，无法针对实验水平比如批次效应进行校正
部分包除了能够进行标准化外，还能够进行降维或差异表达分析，如 BASiCS，zinbwave 和 MAST 等，均提出了特有的统计学框架用于 scRNA-seq 数据的分析。这些框架除了能够针对内部技术性因素进行校正之外，还能够处理像时间、处理以及批次效应等
scone 包可用于比较不同标准化策略的效果，该报可探索不同策略和参数的结果。如果想对 scRNA-seq 数据标准化方法进一步了解，建议阅读 McDavid et al. 和 Vallejos et al. 的综述

针对细胞周期进行调整

除非进行了细胞周期同步化处理，不然群体中的细胞一般处于不同的细胞周期阶段，可能会掩盖主要的生物学效应或者与之高度关联，从而影响分析结果。因此 scRNA-seq 数据分析特有的标准化步骤就是针对细胞周期差异进行校正，比如通过哑变量模型来消除这些异质性。scran 包基于人类和小鼠已知的细胞周期基因训练了一个分类器，能够给细胞赋上 G1，G2/M 和 S 期的分支，在后续分析中将这些分支作为协变量以消除细胞周期的效应。此外，Oscope 包能够在单个细胞的 mRNA 表达量随着细胞周期上下波动时鉴定出振荡器，或者震荡基因

合并数据集

随着 scRNA-seq 技术的广泛应用，整合不同来源的数据集将成为常态。虽然线性或广义线性模型等方法能够整合不同的数据集，但对于 scRNA-seq 数据而言，这些方法的表现可能不佳，因为它们通常默认不同细胞群的组成是已知或者相同的，但这不一定是事实，因此我们需要开发新的方法
一种方法是利用相互最近邻居（mutual nearest neighbors, MNNs）方法找出批次间生物学层面最相似的细胞对，该方法在 scran 和 batchelor 包中均有实现。MNN 校正后的数据可以直接用于下游的分析，类似于其他的降维方法。scMerge 包采用了类似的方法，不过换成找最相似的分支，但是 scMerge 包需要用户提供分支的数量或者已经注释好细胞类型作为先验知识。还有就是 scmap 包，通过近似最近邻居搜索的方法，将查询数据集中的细胞或分支映射至另一个数据集的细胞或分支上。其他的方法包括 scAlign 包中的无监督式深度学习方法以及 mixOmics 包中多种针对混池数据的整合方法
完成整合后，需要注意后续进一步验证以确保整合结果的可靠性。除了消除批次效应外，整合方法的另一应用是将新产生的数据与参考数据集如 Human Cell Atlas，Tabula Muris 等资源进行比较，从而帮助注释这些新数据

特征选择

大部分实验中，只有一部分基因的表达存在异质性。特征选择（feature selection），即指找出这部分导致生物学差异的，能提供有用信息的特征或基因。通过特征选择，能够有效减少计算资源负荷，降低噪音（指那些变化没有规律或者基本无变化的基因）干扰
如果已知细胞分组，那么可采用基于相关性的方法或鉴定组间差异表达基因来挑选特征。但实际中我们一般无法获得这些先验知识，因此可采用的策略包括以下几种：

  • 基于表达量，挑选整体平均表达量最高的那些基因
  • 基于方差，挑选高可变基因（highly variable gene）
  • 基于 droupout，挑选 droupout 较多的基因
  • 基于异常行为（deviance），即根据基因拟合常量表达的空模型（null model）的情况来挑选基因
  • 上述策略的混合

不同方法间的比较可以看一下 Yip et al. 和 Andrews et al. 的综述

降维

测序的细胞和基因数量大幅度上升为后续分析带来了巨大的挑战，容易导致“维度灾难”。特征选择可以稍微减轻这一问题，但还远远不够，因此我们需要通过降维在低维度空间中表示这些数据，并保留有意义的结构。降维结果可用于可视化以及聚类和轨迹推断等下游分析
如前所述，SingleCellExperiment 类提供了专用的 reducedDims slot 用于存放降维结果。降维方法包括 PCA、t-SNE（Rtsne 包中实现）、UMAP（umap 包）和传播图（destiny 包）等。此外，zinbwave 包基于针对 zero-inflated 计数数据的 ZINB-WaVE 模型实现了另一种降维策略。对于开发者而言，BiocSingular 包提供了准确和近似奇异值分解（SVD）的实现，有助于他们在自己的包中实现各种形式的 SVD

下游分析

聚类

单细胞数据能够帮助研究者揭示组织异质性，鉴定新的细胞类型，尤其是那些混池测序无法检测到的罕见细胞群体。通过无监督聚类，将相似的细胞聚集在一起，是解析单细胞数据异质性的基础。scRNA-seq 数据的特性给聚类带来了新的挑战，因此一些新的软件包整合了最近在最近邻居算法方面的进展，例如通过近似方法牺牲一些准确性来提升计算效率。此外，近似方法能够减轻噪音和稀疏性，可能能更好地拟合数据
SC3 和 clusterExperiment 包是两大无监督式聚类框架；SIMLR 包利用 kernels 获得细胞间的距离度量，适用于充满噪音的数据集；对于大型数据集，mbkmeans 包实现了小批量的 k 均值聚类，能够高效地复现经典 k 均值聚类的结果；对于使用了 spike-in 作为对照基因的实验，BEARscc 包能够估算由于技术噪音导致的分支变异。SC3，clusterExperiment 和 clustree 包能够在视觉上以及定量地评估聚类结果，从而比较不同方法和参数的表现
每种方法有它适用的场景，因此在选择时需要注意。此外，结果也需要通过各种实例进行验证。

差异表达

差异表达分析可以确定能区分每个分支的特征，从而用于确定细胞的身份；或是比较不同条件下（处理和对照）的差异，确定参与的基因。差异表达分析的方法主要有两大类。一类是对单细胞数据集进行改造，使之适用于最初设计用于混池数据的分析框架如 edgeR，DESeq2 和 limma 等。以 zinbwave 包为例，它能对基于 zero-inflated 负二项分布（zero-inflated negative binomial distribution, ZiNB）产生的单细胞数据建模，从而解释其固有的稀疏性，换句话说，它能在散度分布和模型拟合步骤中对过量的零值进行降权，从而提升差异表达分析的效果
另一类就是专门设计用于单细胞数据分析的方法，这些方法显式地将基因表达分为两部分，一个是离散部分，描述一个二进制成分（零 vs 非零）的频率，另一个是连续部分，即测量的基因表达水平。MAST 包使用了栅栏模型框架，而 scDD，BASiCS 和 SCDE 包则使用了贝叶斯混合和分层模型

轨迹分析

不同于之前的离散分组，有些情况下异质性用连续谱解释可能更好，因此有了降维的一个特殊应用——轨迹分析，也称拟时间推测，从而鉴定新的细胞亚型，分化过程或者导致分叉的事件。这类分析是单细胞特有的，自2014年起涌现出大量的软件，Saelens et al. 的综述对这些方法进行了比较，并给出了实用软件选择建议。需要记住的是每种方法都有其适用性，可以尝试使用多种方法以获得更加可靠的结果

注释

高通量数据分析中常遇到的一个挑战就是如何将数据表征为我们熟悉的术语，因此基因特征（gene signature）的概念应运而生，即在一组基因的环境中描述单个基因的表达变化。所有基因特征的基本结构是它们提供一组具有共同生物背景的基因，这可以是通过代谢或分子通路数据库。更复杂的版本可以包含定量或定性指标，或基因特征间关系等额外的信息。最近一些新开发的方法绕过了基因特征的使用，转而采取以数据为中心的策略，基于参考数据集注释新数据

基因特征富集

目前已经开发了大量公共数据库用于下游富集分析，常规的如 KEGG，Reactome 和 Gene Ontology（GO），不同免疫学模块条件下差异表达的 MsigDB 等
类似于单细胞数据中的差异表达分析，富集分析方法也可分为两大类。第一类是对用于芯片和 bulk 数据的方法基于权重进行调整，以解释单细胞数据中大量的零值

  • GSEA（或预计算过的 fqsea 包）
  • goseq
  • PADOG
  • EnrichmenBrowser 包提供了10种不同方法，而且提供选项用于方法间结果的排序

第二类方法则是专门用于 scRNA-seq 数据：

  • MAST：利用栅栏模型实现了一个竞争性基因集测试，用于解释基因间关联
  • AUCell：利用基于排名的打分方法对基因集的活跃水平打分，并计算每个细胞的基因集活性得分
  • slalom：利用因子花大细胞因变量模型解释单细胞数据集中的变异

以数据为中心的富集方法

以数据为中心的方法通过使用更加全面的量化特征定义来实现更佳的注释，目前其中一个特殊的应用是基于已知的细胞类型将未知的细胞进行分类。基于是给每个细胞还是（预先定义的）细胞群贴标签，采用的方法也会不同。在分支水平，可以克服稀疏性，降低特征量化的不确定性，但是对分支的定义本身内在就存在经验性的过程，因而可能导致后续结果也会存在偏好性。目前可用的包包括

  • celaref
  • scmap
  • cellassign

可重复的分析

交互式数据可视化

iSEE 包使用户能够通过浏览器对 scRNA-seq 数据集进行交互式可视化，它能够直接与 SingleCellExperiment 对象交互，如可视化降维结果等，还允许用户将生成交互图表的代码输出，有利于探索性分析的可重复性

生成报告

为确保结果具有良好的文档记录，可分享，可重复，我们需要用到 rmarkdown 和 bookdown 包。此外，Bioconductor 还发表了 BiocStyle 包，提供包的说明文档的标准化格式。在质控报告方面，countsimQC 和 batchQC 包可对多个样本的特征，如文库大小，基因数目等进行可视化

发表以及模拟数据

包含有单细胞数据的包有：

  • ExperimentHub：提供了对发表数据进行索引的标准接口
  • TENxPBMCData
  • TENxBrainData
  • HCAData
  • HCABrowser
  • conquer：相较其他的数据是由各自作者提供，处理上可能存在差异，该报提供的数据以一致的方式处理，且有质控步骤

基准比较数据集：

  • CellBench：包含来自不同平台和细胞系的数据集
  • Tung et al 基于平板方法得到的数据
  • DuoClustering2018：包含各种来源的数据，本来是用于比较不同聚类方法的

模拟数据集：

  • splatter：可模拟不同细胞类型，批次效应，droupout 水平，差异表达等等

总结

Bioconductor 为单细胞数据的分析提供了全面的生态环境，而且作为 R 生态的一部分，可以很方便地调用各类 R 包用于统计和数据科学探索。此外，对其他的 R 包如 Seurat或者其他语言如 Python 中的 scanpy 包的结果，也提供了相互转换的方法

REF

  1. Robert A. Amezquita, et al. Orchestrating Single-Cell Analysis with Bioconductor. Preprint at bioRxiv https://www.biorxiv.org/content/10.1101/590562v1 (2019)