scRNA单细胞试验记录

SingleCellExperiment

2019-11-06  本文已影响0人  麒麟991

Bioconductor软件包SingleCellExperiment提供了SingleCellExperiment类以供使用。当使用依赖于SingleCellExperiment类的任何程序包和加载程序包时,可以按以下方式显式安装(并加载)程序包:

BiocManager::install('SingleCellExperiment')
BiocManager::install(c('scater', 'scran', 'uwot'))
library(SingleCellExperiment)

所述SingleCellExperimentsce)对象是基于在Bioconductor的单细胞分析应用的基础。该sce对象是一个S4对象,与R中其他可用的方法相比,它本质上为数据的构造和访问提供了一种更为形式化的方法。

如果我们想象sce对象是一艘船,slots可被认为是个别货箱——为sce对象内单独的实体。此外,每个slot都包含自己的格式数据。打个比喻,我们可以想象水果和砖头需要不同的货箱变化形式。在某些情况下sce,一些slot包含矩阵,而其他slot则是数据框。

image
要构建基本sce对象,我们只需要一个slot:
assays slot:包含基本数据的一个列表,列表中的每个条目均采用矩阵格式,其中行对应于特征(基因),列对应于样本(单元格)(图1A,蓝色框)

让我们从生成10个基因的3个细胞计数数据开始

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

要访问我们刚刚提供的计数数据,我们可以执行以下任一操作:

assay(sce, "counts") ——这是最通用的方法,我们可以在其中提供测定的名称作为第二个参数。
counts(sce)——与上述相同,但仅适用于具有特殊名称的测定"counts"。

counts(sce)

扩展assays slot

使assay slot特别强大的原因是它可以容纳多种数据。这对于存储数据的原始版本和normalized版本特别有用。我们可以如下所示,使用scranscater包来计算初始主数据的对数转换形式。

请注意,这里,sce在将结果重新分配给我们之前,我们将覆盖之前的内容sce——这是因为这些函数返回了SingleCellExperiment对象。某些功能(尤其是那些面向单细胞的Bioconductor封装之外的功能)没有,在这种情况下,您需要将结果附加到sce对象上。

sce <- scran::computeSumFactors(sce)
sce <- scater::normalize(sce)

再次查看该对象,我们看到这些函数添加了一些新条目:

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.43   4.56   4.14
## gene_2    4.83   4.10   4.54
## gene_3    4.04   4.23   4.19
## gene_4    4.31   3.62   4.43
## gene_5    3.70   3.62   4.23
## gene_6    4.43   4.35   3.53
## gene_7    3.26   3.96   3.90
## gene_8    4.54   4.46   4.31
## gene_9    3.70   4.56   4.00
## gene_10   3.50   3.62   4.05
## assay(sce, "logcounts") ## same as above

请注意,之前的数据在细胞1、2与3之间的计数存在严重差异,并且通过归一化可以缓解这种差异。

要查看sce中所有可用的assay,我们可以输入:

assays(sce)
## List of length 2
## names(2): counts logcounts

虽然上面的功能演示了向我们的sce对象自动添加,但是在某些情况下,我们可能希望执行自己的计算并将结果保存到assays中。特别对于使用不返回SingleCellExperiment对象的函数很重要。

让我们附加一个修改后的数据新版本+100

counts_100 <- assay(sce, "counts") + 100
assay(sce, "counts_100") <- counts_100 # assign a new entry to assays slot

然后,我们可以使用访问器assays()(请注意,这是复数!)来查看assay到目前为止的所有条目。请注意,要查看所有条目,我们使用复数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以图形方式表示了slots的这种扩展,显示了将矩阵添加到slot中。

以类似的方式,sce如上所述,许多slot可通过分配而扩展,从而允许与面向单细胞的封装功能包之外的功能互相操作。

列数据:colData

为了进一步注释我们的sce对象,第一个也是最有用的信息之一是添加描述我们原始数据的列,例如实验的样本或细胞名。此数据输入到colData

因此,让我们为细胞提供一些注释数据,从批次变量开始,其中细胞1和2在批次1中,而细胞3在批次2中。

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目录中,例如,scater软件包具有一个功能calculateQCMetrics(),该功能会附加许多质量控制数据。在这里,我们显示colData(sce)前五列,并附加了质量控制指标。

sce <- scater::calculateQCMetrics(sce)
colData(sce)[, 1:5]
## DataFrame with 3 rows and 5 columns
##            batch is_cell_control total_features_by_counts
##        <numeric>       <logical>                <integer>
## cell_1         1           FALSE                       10
## cell_2         1           FALSE                       10
## cell_3         2           FALSE                       10
##        log10_total_features_by_counts total_counts
##                             <numeric>    <integer>
## cell_1               1.04139268515822           98
## cell_2               1.04139268515822          104
## cell_3               1.04139268515822          303

使用colData的子集

colData的一个常见的操作是它在子集选取中的使用。一种简单的访问方式是通过使用$运算符,这是访问colData内变量的捷径:

sce$batch
## [1] 1 1 2
## colData(sce)$batch # same as above

如果只需要批次1中的单元格,则可以按如下方式对sce对象进行子集化(请记住,在这种情况下,我们在上具有子集,因为这里是按细胞/样本进行过滤)。

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):

行数据:rowData/rowRanges

最后,这些行还具有自己的注释信息单元,用于存储与sce对象功能有关的信息:

此外,还有一个特殊的slot,储存基因组坐标特征:

这两种方法都能通过各自的存取访问,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)填充在对象的rowData中,如下所示(为简便起见,仅显示了前三列):

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 18.6666666666667  1.29373075692248
## gene_2               FALSE               22  1.36172783601759
## gene_3               FALSE               17  1.25527250510331
## gene_4               FALSE 18.3333333333333  1.28630673884327
## gene_5               FALSE 15.3333333333333  1.21307482530885
## gene_6               FALSE 14.3333333333333  1.18563657696191
## gene_7               FALSE               13  1.14612803567824
## gene_8               FALSE               20  1.32221929473392
## gene_9               FALSE               16  1.23044892137827
## gene_10              FALSE 13.6666666666667  1.16633142176653

以类似于colData的方式,可以在创建SingleCellExperiment对象时在开始时就提供此类功能数据,我们将其留给读者作为练习。

按行取子集

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):
## sce[c(1, 4), ] # same as above in this case

sizeFactors

我们已经通过scran::computeSumFactors(sce)调用添加了一个sizeFactorsslot:

sce <- scran::computeSumFactors(sce)
sce <- scater::normalize(sce)
sizeFactors(sce)
## [1] 0.582 0.618 1.800

简要回顾:从sesce

到目前为止,我们已经涵盖了assays(原始数据),colData(样本注释数据),rowData/ rowRanges(特征注释数据)和SingleCellExperimentsizeFactors

需要注意的是,SingleCellExperiment是从它的前身SummarizedExperimentse)导出的,特别是继承了assayscolDatarowData/ rowRanges。这样,大多数SummarizedExperiment功能保留在SingleCellExperiment中。这使与之配合使用的现有方法可以`在SingleCellExperiment对象上类似地工作。

那么SingleCellExperiment又有什么创新呢?对于我们的讨论,最重要的更改是添加了一个reducedDims

reducedDims slot

reducedDims是一个新增功能,专门用于存储通过PCA,tSNE,UMAP等方法获得的原始数据的降维信息。

最重要的是,就像assays一样,该reducedDims是可以容纳许多条目的列表。因此,它可以在reducedDims插槽内保存给定数据集的PCA,TSNE和UMAP表示形式。

在我们的示例中,我们可以使用scaterpackage函数如下计算数据的PCA表示形式runPCA()。我们看到,sce现在显示了一个新的reducedDim,并且reducedDim()对来自logcounts(sce)的规范化数据产生了运行PCA的结果。

sce <- scater::runPCA(sce)
reducedDim(sce, "PCA")
##           PC1    PC2
## cell_1  0.194  0.818
## cell_2 -0.887 -0.258
## cell_3  0.693 -0.560
## attr(,"percentVar")
## [1] 55.4 44.6

由此,我们还可以使用scaterpackagerunTSNE()函数来计算tSNE表示形式,并且可以在默认视图中sce或通过access 看到它:

sce <- scater::runTSNE(sce, perplexity = 0.1)
## Perplexity should be lower than K!
reducedDim(sce, "TSNE")
##         [,1]  [,2]
## cell_1  1341  5533
## cell_2 -5458 -1597
## cell_3  4117 -3936

我们可以通过访问器查看reducedDims中所有条目的名称reducedDims()(请注意,该名称是复数形式)

reducedDims(sce)
## List of length 2
## names(2): PCA TSNE

现在,假设我们要尝试使用一种不同的降维算法,该算法尚未实现与SingleCellExperiment的直接兼容。为了适应这种情况(或者,当我们想直接对数据本身而不是通过包运行降维方法时),我们可以直接添加到reducedDims中。这类似于我们之前assays使用自定义条目扩展counts_100的方式。

下面,我们展示如何直接umap()uwot包中实现,而不是通过在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.413  0.0368
## cell_2  0.243  0.5494
## cell_3  0.170 -0.5862
## attr(,"scaled:center")
## [1] -2.32 -1.64

现在,当我们查看reducedDims()输出时,我们还可以看到其条目:

reducedDims(sce)
## List of length 3
## names(3): PCA TSNE UMAP_uwot

metadata

一些分析产生的结果不适合上述slot。值得庆幸的是,有一个slot仅用于这种类型的混乱数据,并且实际上,可以容纳任何类型的数据,只要它在命名列表中即可:

例如,假设我们有一些喜欢的基因,例如高度可变的基因,我们希望将其保存在内部sce以便以后使用。我们可以简单地通过如下所示将其附加到metadata 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"

SingleCellExperiment记录分析

在随后的部分中,我们将显示一个示例工作流,该工作流使用SingleCellExperiment对象作为其基础,并且与上面的SingleCellExperiment类的演练类似,将连续追加新条目以保存分析结果。因此,SingleCellExperiment可以以此方式用作分析记录。这使得它对于合作特别有用,因为可以通过iSEE等图形用户界面来传输,分析甚至可视化对象。

结论

这种与SingleCellExperiment连接,使得许多这些程序包在scRNA-seq分析的整个过程中都易于互操作和模块化。此外,它允许任何人在SingleCellExperiment里面实施并建立自己的scRNA-SEQ分析工具。

原始学习网站:https://osca.bioconductor.org/data-infrastructure.html

上一篇下一篇

猜你喜欢

热点阅读