希望大家能有所收获!
序言
通过Bioconductor进行单细胞分析(一)
本教程是翻译自 Orchestrating Single-Cell Analysis with Bioconductor(https://osca.bioconductor.org/), 由于前三章是序言和软件介绍等无关紧要的内容,所以我们从第四章数据基础框架开始
正
文
第四章 数据基础架构
使用Bioconductor软件包的优势之一是它们利用了通用的数据基础架构,从而使分析可以在各种软件包之间实现互操作。此外,要使此基础架构强大且可扩展,需要付出大量工程努力。在这里,我们将详细描述
使用Bioconductor软件包的优势之一是它们利用了通用的数据基础架构,从而使分析可以在各种软件包之间实现互操作。此外,要使此基础架构强大且可扩展,需要付出大量工程努力。在这里,我们将详细描述SingleCellExperiment对象(或简写成sce`形式),以描述如何构造,在下游分析中使用该对象以及如何存储各种类型的主数据和元数据。
4.1
先决条件
按以下方式安装SingleCellExperiment包:
代码语言:javascript复制BiocManager::install('SingleCellExperiment')
此外,还会用到scater
,scran
和CRAN包uwot
。安装命令如下:
BiocManager::install(c('scater', 'scran', 'uwot'))
对于本节,需要加载的是SingleCellExperiment包:
代码语言:javascript复制library(SingleCellExperiment)
4.2
SingleCellExperiment Class
图1 Overview of the SingleCellExperiment class object
4.2.1 Primary Data:The `assays` Slot
SingleCellExperiment(sce
)对象是基于在Bioconductor的单细胞分析应用的基础。该sce
对象是一个S4对象(https://adv-r.hadley.nz/s4.html),与R中其他可用的方法相比,它本质上为数据的构造和访问提供了一种更为形式化的方法。S4对象最大最重要的用处就是可以安全有效检查数据并且能够扩展数据内容。
假设的sce
是一艘船,sce
的slot可被认为是货箱——可以作为sce
中独立的实体。并且每个slot都以自己的格式保存数据。可以理解为水果和砖头需要不同的货箱。在sce
中,某些slot需要数字矩阵,而其他slot可能需要数据框。
要构建基本sce
对象,我们需要的是一个slot:
assays
slot:包含基本数据,例如列表中的counts,列表中的每个条目均采用矩阵格式,其中行对应于特征(基因),列对应于样本(cell)(图1A,蓝色框)
让我们从生成3个细胞,10个基因的简单数据开始:
代码语言:javascript复制counts_matrix <- data.frame(cell_1 = rpois(10, 10),
cell_2 = rpois(10, 10),
cell_3 = rpois(10, 30))
rownames(counts_matrix) <- paste0("gene_", 1:10)
counts_matrix <- as.matrix(counts_matrix) # must be a matrix object!
现在可以构建我们的第一个SingleCellExperiment对象,使用已定义的构造函数,SingleCellExperiment()
。请注意,我们以命名列表的形式提供数据,并且列表的每个条目都是一个矩阵。在这里,我们将counts_matrix
命名为counts
。
sce <- SingleCellExperiment(assays = list(counts = counts_matrix))
要查看对象,我们可以简单地在控制台输入sce
以查看一些相关信息,这些信息将显示可供我们使用的各种slot的概述(可能有也可能没有数据)。
sce
## class: SingleCellExperiment
## dim: 10 3
## metadata(0):
## assays(1): counts
## rownames(10): gene_1 gene_2 ... gene_9 gene_10
## rowData names(0):
## colnames(3): cell_1 cell_2 cell_3
## colData names(0):
## reducedDimNames(0):
## spikeNames(0):
## altExpNames(0):
要访问我们刚刚提供的计数数据,我们可以执行以下任一操作:
assay(sce, "counts")
这是最通用的方法,我们可以提供测定的名称作为第二个参数。counts(sce)
与上述相同,但仅适用于"counts"
。
counts(sce)
## cell_1 cell_2 cell_3
## gene_1 10 5 37
## gene_2 7 11 30
## gene_3 5 12 34
## gene_4 12 10 27
## gene_5 12 11 32
## gene_6 12 8 28
## gene_7 7 13 22
## gene_8 16 11 32
## gene_9 17 10 29
## gene_10 10 8 29
代码语言:javascript复制## assay(sce, "counts") ## same as above in this special case
4.2.2 Extending the `assays ` Slot
使assay
slot特别强大的原因是它可以容纳主要数据的多种表示形式。这对于存储数据的原始版本和规范化版本特别有用。我们可以使用scran
和scater
包来计算初始主数据的规范化和用对数转换表示初始数据。
请注意,您需要将结果附加到sce
对象上。
sce <- scran::computeSumFactors(sce)
sce <- scater::normalize(sce)
再次查看该对象,我们看到这些函数添加了一些新条目:
代码语言:javascript复制sce
## class: SingleCellExperiment
## dim: 10 3
## metadata(1): log.exprs.offset
## assays(2): counts logcounts
## rownames(10): gene_1 gene_2 ... gene_9 gene_10
## rowData names(0):
## colnames(3): cell_1 cell_2 cell_3
## colData names(0):
## reducedDimNames(0):
## spikeNames(0):
## altExpNames(0):
具体来说,我们看到该assays
slot 已由两个条目组成:counts
(初始数据)和logcounts
(规范化数据)。与counts
相似,该logcounts
名称也是一个可以直接访问的特殊名称logcounts(sce)
。
logcounts(sce)
## cell_1 cell_2 cell_3
## gene_1 4.06 3.25 4.45
## gene_2 3.58 4.31 4.16
## gene_3 3.14 4.43 4.33
## gene_4 4.31 4.18 4.02
## gene_5 4.31 4.31 4.25
## gene_6 4.31 3.87 4.07
## gene_7 3.58 4.54 3.74
## gene_8 4.70 4.31 4.25
## gene_9 4.79 4.18 4.12
## gene_10 4.06 3.87 4.12
代码语言:javascript复制## assay(sce, "logcounts") ## same as above
要查看中所有可用的检测方法sce
,我们可以输入:
assays(sce)
## List of length 2
## names(2): counts logcounts
虽然上面的功能演示了sce
的assays
可以自动添加,但是在某些情况下,我们可能希望执行自己的计算并将结果保存到assays
slot中。让我们附加偏移的数据,在原来的基础上 100
。
counts_100 <- assay(sce, "counts") 100
assay(sce, "counts_100") <- counts_100 # assign a new entry to assays slot
然后,我们可以使用访问器assays()
(请注意,这是复数!)来查看到目前为止assay
在slot中的所有条目。请注意,要查看所有测定,我们使用复数assays()
访问器,并使用单数访问器assay()
检索单个测定条目(作为矩阵),并提供想要检索的测定名称。
assays(sce)
## List of length 3
## names(3): counts logcounts counts_100
这些条目也可以在默认视图中看到sce
:
sce
## class: SingleCellExperiment
## dim: 10 3
## metadata(1): log.exprs.offset
## assays(3): counts logcounts counts_100
## rownames(10): gene_1 gene_2 ... gene_9 gene_10
## rowData names(0):
## colnames(3): cell_1 cell_2 cell_3
## colData names(0):
## reducedDimNames(0):
## spikeNames(0):
## altExpNames(0):
图1B(深蓝色框)assays
以图形方式表示了slot的这种扩展,显示了将logcounts
矩阵添加到assays
slot中。
以类似的方式,sce
许多slot可通过分配而扩展,从而允许与面向单细胞的生物导体封装之外的功能互操作所需的多种定制功能。
4.2.3 Column (Meta)Data: `colData` Slot
为了进一步注释我们的sce
对象,第一个也是最有用的信息是添加描述主要数据列的元数据,例如实验的样本或cell 。此数据输入到colData
slot:
colData
slot:提供描述样本(cell)的元数据,以data.frame
(或DataFrame
)形式保存,其中行对应于cell,列对应于样本(cell)元数据功能(例如id,批处理,作者等)(图1A ,橙色框)。
因此,让我们为cell提供一些元数据,从批处理变量开始,其中cell1和2在批处理1中,而cell3在批处理2中。
代码语言:javascript复制cell_metadata <- data.frame(batch = c(1, 1, 2))
rownames(cell_metadata) <- paste0("cell_", 1:3)
现在,我们可以采用两种方法将cell_metadata
追加到我们现有的sce
,或者通过SingleCellExperiment()
构造函数从头开始并从一开始就提供。现在,我们将从头开始,但还将展示如何添加数据:
## From scratch:
sce <- SingleCellExperiment(assays = list(counts = counts_matrix),
colData = cell_metadata)
## Appending to existing object (requires DataFrame() coercion)
## colData(sce) <- DataFrame(cell_metadata)
与相似assays
,我们可以看到colData默认填充到sce
:
sce
## class: SingleCellExperiment
## dim: 10 3
## metadata(0):
## assays(1): counts
## rownames(10): gene_1 gene_2 ... gene_9 gene_10
## rowData names(0):
## colnames(3): cell_1 cell_2 cell_3
## colData names(1): batch
## reducedDimNames(0):
## spikeNames(0):
## altExpNames(0):
通过colData()
访问数据:
colData(sce)
## DataFrame with 3 rows and 1 column
## batch
## <numeric>
## cell_1 1
## cell_2 1
## cell_3 2
最后,一些软件包会自动添加到colData
slot中,例如,该scater
中的calculateQCMetrics()
函数,该函数会附加许多质量控制数据。在这里,我们显示附加了质量控制指标的前五列。
sce <- scater::calculateQCMetrics(sce)
colData(sce)[, 1:5]
## DataFrame with 3 rows and 5 columns
## batch is_cell_control total_features_by_counts log10_total_features_by_counts
## <numeric> <logical> <integer> <numeric>
## cell_1 1 FALSE 10 1.04139268515822
## cell_2 1 FALSE 10 1.04139268515822
## cell_3 2 FALSE 10 1.04139268515822
## total_counts
## <integer>
## cell_1 108
## cell_2 99
## cell_3 300
4.2.3.1 Using `colData` for Subsetting
colData
最常用的功能是提取子集。一种简单的访问colData
的方式是通过使用$
运算符,这是访问colData
slot内变量的捷径:
sce$batch
## [1] 1 1 2
代码语言:javascript复制## colData(sce)$batch # same as above
如果只需要批次1中的cell,则可以按如下方式对sce对象进行取子集(我们可以在列上取子集,是因为在这里是通过cell/样本进行过滤的)。
代码语言:javascript复制sce[, sce$batch == 1]
## class: SingleCellExperiment
## dim: 10 2
## metadata(0):
## assays(1): counts
## rownames(10): gene_1 gene_2 ... gene_9 gene_10
## rowData names(7): is_feature_control mean_counts ... total_counts log10_total_counts
## colnames(2): cell_1 cell_2
## colData names(10): batch is_cell_control ... pct_counts_in_top_200_features
## pct_counts_in_top_500_features
## reducedDimNames(0):
## spikeNames(0):
## altExpNames(0):
4.2.4 Feature Metadata: `rowData`/`rowRanges`
最后,这些行还具有自己的元数据slot,用于存储与sce
对象功能有关的信息:
rowData
slot:包含data.frame
(DataFrame
)格式的数据,该数据描述与主要数据行相对应的方面(图1A,绿色框)。
此外,还有一个特殊的slot,适用于具有基因组坐标的特征:
rowRanges
slot:以GRangesList
(其中每个条目均为GenomicRanges
格式)的数据形式描述染色体,开始和结束坐标,特征(基因,基因组区域)等。
这两种方法都能通过各自的存取访问,rowRanges()
和rowData()
。在我们的例子中,rowRanges(sce)
产生一个空列表:
rowRanges(sce) # empty
## GRangesList object of length 10:
## $gene_1
## GRanges object with 0 ranges and 0 metadata columns:
## seqnames ranges strand
## <Rle> <IRanges> <Rle>
## -------
## seqinfo: no sequences
##
## $gene_2
## GRanges object with 0 ranges and 0 metadata columns:
## seqnames ranges strand
## <Rle> <IRanges> <Rle>
## -------
## seqinfo: no sequences
##
## $gene_3
## GRanges object with 0 ranges and 0 metadata columns:
## seqnames ranges strand
## <Rle> <IRanges> <Rle>
## -------
## seqinfo: no sequences
##
## ...
## <7 more elements>
上一节中调用calculateQCMetrics(sce)
的结果保存在sce
对象的rowData
slot中,如下所示(为简便起见,仅显示了前三列):
rowData(sce)[, 1:3]
## DataFrame with 10 rows and 3 columns
## is_feature_control mean_counts log10_mean_counts
## <logical> <numeric> <numeric>
## gene_1 FALSE 17.3333333333333 1.26324143477458
## gene_2 FALSE 16 1.23044892137827
## gene_3 FALSE 17 1.25527250510331
## gene_4 FALSE 16.3333333333333 1.23888208891514
## gene_5 FALSE 18.3333333333333 1.28630673884327
## gene_6 FALSE 16 1.23044892137827
## gene_7 FALSE 14 1.17609125905568
## gene_8 FALSE 19.6666666666667 1.31527043477859
## gene_9 FALSE 18.6666666666667 1.29373075692248
## gene_10 FALSE 15.6666666666667 1.22184874961636
以类似于colData
slot的方式,可以在创建SingleCellExperiment对象时在开始时提供此类功能元数据,我们将其留给读者作为练习。
4.2.4.1 Subsetting by Rows
要sce
在特征/基因级别上将对象细分,我们可以通过提供数字索引或名称向量来进行类似于其他R对象的行设置操作:
sce[c("gene_1", "gene_4"), ]
## class: SingleCellExperiment
## dim: 2 3
## metadata(0):
## assays(1): counts
## rownames(2): gene_1 gene_4
## rowData names(7): is_feature_control mean_counts ... total_counts log10_total_counts
## colnames(3): cell_1 cell_2 cell_3
## colData names(10): batch is_cell_control ... pct_counts_in_top_200_features
## pct_counts_in_top_500_features
## reducedDimNames(0):
## spikeNames(0):
## altExpNames(0):
代码语言:javascript复制## sce[c(1, 4), ] # same as above in this case
4.2.5 Size Factors Slot: `sizeFactors`
简而言之,我们已经通过调用scran::computeSumFactors(sce)
(添加了一个sizeFactors slot)
:
sizeFactors
slot:在数字向量中包含有关样本/细胞归一化因子的信息,该因子用于产生归一化的数据(图1B,棕色框)
sce <- scran::computeSumFactors(sce)
sce <- scater::normalize(sce)
sizeFactors(sce)
## [1] 0.639 0.586 1.775
4.3
A Brief Recap: From `se` to `sce`
到目前为止,我们已经涵盖了assays
(原始数据),colData
(样本元数据),rowData
/ rowRanges
(特征元数据)和SingleCellExperiment的sizeFactors
slot。
需要注意的是,SingleCellExperiment类是从SummarizedExperiment
(se
)类导出,并继承assays
,colData
和rowData
/ rowRanges
slot。大多数SummarizedExperiment
功能保留在SingleCellExperiment中。这使与之配合使用的现有方法在SingleCellExperiment对象上类似地工作。
那么SingleCellExperiment类又有什么新消息呢?对于我们的讨论,最重要的更改是添加了一个名为reducedDims
的新slot。
4.4
The `reducedDims` Slot
reducedDims
slot是一个新增功能,专门用于存储通过PCA,tSNE,UMAP等方法获得的原始数据的降维表示。
reducedDims
slot:包含数字matrix
条目的列表,这些条目描述了降维的原始数据表示,因此行代表原始数据的列(样本/cell),而列则代表维度
最重要的是,就像assays
slot一样,该reducedDims
slot可以容纳许多条目的列表。因此,它可以在reducedDims
slot内保存给定数据集的PCA,TSNE和UMAP表示形式。
在我们的示例中,我们可以使用scater
中的函数runPCA()
计算数据的PCA。我们看到,sce
现在显示了一个新的reducedDim
,并且通过访问器reducedDim()
对来自的规范化数据logcounts(sce)
产生了运行PCA的结果。
sce <- scater::runPCA(sce)
reducedDim(sce, "PCA")
## PC1 PC2
## cell_1 -1.126 -0.309
## cell_2 0.991 -0.455
## cell_3 0.135 0.764
## attr(,"percentVar")
## [1] 71.9 28.1
由此,我们还可以使用scater
函数runTSNE()
来计算tSNE表示形式,并且可以在默认视图中sce
看到它:
sce <- scater::runTSNE(sce, perplexity = 0.1)
## Perplexity should be lower than K!
代码语言:javascript复制reducedDim(sce, "TSNE")
## [,1] [,2]
## cell_1 1231 -5557
## cell_2 4197 3859
## cell_3 -5428 1699
我们可以通过访问器reducedDims()
查看reducedDims
slot中所有条目的名称(请注意,该名称是复数形式,与reducedDim()
不同:
reducedDims(sce)
## List of length 2
## names(2): PCA TSNE
现在,假设我们要尝试使用一种不同的降维算法,该算法尚未实现与SingleCellExperiment的直接兼容性,例如,使用scater
。为了适应这种情况(或者,当我们想直接对数据本身而不是通过包装器运行降维方法时),我们可以直接添加到reducedDims
slot中。这类似于我们之前assays
使用自定义条目counts_100
扩展slot的方式。
下面,我们展示如何直接通过uwot
包中的umap()
,而不是通过scater
关联包装函数runUMAP()
来运行,并保存中间结果,然后将它们添加到sce
我们先前拥有的对象中。
u <- uwot::umap(t(logcounts(sce)), n_neighbors = 2)
reducedDim(sce, "UMAP_uwot") <- u
reducedDim(sce, "UMAP_uwot")
## [,1] [,2]
## cell_1 0.526 -0.402
## cell_2 -0.266 0.632
## cell_3 -0.260 -0.230
## attr(,"scaled:center")
## [1] 9.12 15.20
现在,当我们查看访问reducedDims()
器输出时,我们还可以看到其条目:
reducedDims(sce)
## List of length 3
## names(3): PCA TSNE UMAP_uwot
4.5
The `metadata` Slot
一些分析产生的结果不适合上述slot。值得庆幸的是,有一个slot仅用于这种类型的混乱数据,实际上它可以容纳任何类型的数据,只要它在命名列表中即可:
metadata
slot:一个命名的条目列表,列表中的每个条目都可以是您想要的任何内容
例如,假设我们有一些喜欢的基因,例如高度可变的基因,我们希望将其保存在内部sce
以便以后使用。我们可以简单地通过如下所示将其附加到元数据slot来做到这一点:
my_genes <- c("gene_1", "gene_5")
metadata(sce) <- list(favorite_genes = my_genes)
metadata(sce)
## $favorite_genes
## [1] "gene_1" "gene_5"
同样,我们可以通过$
运算符附加更多信息:
your_genes <- c("gene_4", "gene_8")
metadata(sce)$your_genes <- your_genes
metadata(sce)
## $favorite_genes
## [1] "gene_1" "gene_5"
##
## $your_genes
## [1] "gene_4" "gene_8"
4.6
About Spike-Ins
您可能已经注意到sce
默认视图中带有spikeNames
。SingleCellExperiment对象对通过掺入(ERCC)对照 的实验数据预留了通道。我们将此留给有兴趣的读者,以在SingleCellExperiment入门小插图中(https://bioconductor.org/packages/release/bioc/vignettes/SingleCellExperiment/inst/doc/intro.html)了解更多信息。
4.7
Recording Analyses in SingleCellExperiment
在随后的部分中,我们将显示一个示例工作流,该工作流使用SingleCellExperiment对象作为其基础,并且与上面的SingleCellExperiment类的演练类似,将连续追加新条目以保存分析结果。因此,SingleCellExperiment可以以此方式用作分析记录。这使得它对于协作特别有用,因为可以通过图形用户界面(例如iSEE)来传输,分析甚至可视化对象。
4.8
结论
Bioc-verse单细胞相关的软件包,与SingleCellExperiment class可以连接。正是这种与SingleCellExperiment的连接,使得许多这些程序包在scRNA-seq分析的整个过程中都易于互操作和模块化。此外,它允许任何人以SingleCellExperiment为基础,建立自己的scRNA-SEQ分析工具,促进发展。