背景介绍
Bioconductor项目的主要优势之一在于使用通用数据基础结构，该数据基础结构可实现跨越R包的相互操作性。 用户应该能够使用来自不同Bioconductor软件包的功能来分析其数据，而无需在格式之间进行转换。 为此，SingleCellExperiment类（来自SingleCellExperiment包）数据结构可以当做在70多个单细胞相关的Bioconductor包进行数据分析的“通用货币”。 此类实现了一个数据结构，该数据结构存储了我们单细胞数据的各个方面-逐个基因表达数据，每个细胞meta数据和每个基因注释（图1），并以同步方式对其进行操作。
图1：SingleCellExperiment类的结构概述。 分析的每一行对应于rowData的一行（粉红色阴影），而分析的每一列对应于colData和reduceDims的列（黄色阴影）
使用任何依赖于SingleCellExperiment类的软件包时，都会调用预先安装的SingleCellExperiment软件包，但也可以如下方式安装（和加载）该软件包：
BiocManager :: install（'SingleCellExperiment'）
另外，我们还会使用scater和scran软件包以及CRAN软件包uwot的某些功能（也可以通过BiocManager :: install方便地安装它们）。这些功能将根据需要通过<package> :: <function>进行调用。
BiocManager :: install（c（'scater'，'scran'，'uwot'））
通常，我们将SingleCellExperiment包加载到我们的R会话中。这样就避免了在函数调用前加上::的前缀，特别是对于在整个分析工作流程中大量使用的包。
libarary(SingleCellExperiment)
SingleCellExperiment对象中的每条（元）数据都由一个单独的“插槽”表示。 （此术语来自S4类系统，但现在并不重要。）如果我们将SingleCellExperiment对象想象为一艘货船，则可以将这些插槽视为具有不同内容的单个货箱，例如，某些插槽装载数字矩阵，而其他矩阵可能装载数据框。在本章的其余部分，我们将讨论可用的插槽，它们装载的数据格式以及如何与之交互。
存储初始的实验数据
1填充assays插槽
要构建基本的SingleCellExperiment对象，我们只需要填充assays槽即可。其中包含主要数据，例如测序计数矩阵，其中行对应于特征（基因），列对应于样本（细胞）（图1，蓝色框）。让我们从十个基因中生成三个单元格的计数数据开始，从简单开始：
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）＃必须是矩阵对象！
由此，我们现在可以使用SingleCellExperiment（）函数构造第一个SingleCellExperiment对象。请注意，我们将数据作为命名列表提供，其中列表的每个条目都是一个矩阵。在这里，我们将counts_matrix条目简称为“ counts”。
sce <-SingleCellExperiment（assays = list（counts = counts_matrix））
要查看对象，我们只需在控制台中键入sce即可查看一些相关信息，这些信息将显示可供我们使用的各种插槽的概述（可能有或没有任何数据）。
sce
##class：SingleCellExperiment
##dim：10 3
##metadata（0）：
## assays（1）：counts
## rownames（10）：gene_1gene_2 ... gene_9gene_10
## rowData names（0）：
## colnames（3）：cell_1 cell_2 cell_3
## colData names（0）：
## reduceDimNames（0）：
## altExpNames（0）：
要访问我们刚刚提供的计数数据，我们可以执行以下任一操作：
-
assay（sce，“ counts”)-这是最通用的方法，我们可以在其中提供assay名称作为第二个参数。 -
counts（sce）-这是上面的快捷方式，但仅适用于特殊名称“ counts”。
counts（sce）
## cell_1 cell_2 cell_3
## gene_1 8 8 33
## gene_2 13 9 21
## gene_3 15 8 33
## gene_4 7 9 40
## gene_5 8 8 27
## gene_6 10 10 39
## gene_7 11 8 21
## gene_8 14 13 41
## gene_9 8 10 45
## gene_10 12 10 31
2添加更多的“assays”
使“assays”槽特别强大的原因在于，它可以容纳原始数据的多种表示形式。 这对于存储原始计数矩阵以及数据的normalized版本特别有用。 我们可以如下所示使用scater包来执行原始计数的normalize，以算出初始数据的规范化和对数转换表示形式。
sce <-scater :: logNormCounts（sce）
请注意，在每个步骤中，我们都会通过将结果重新赋给sce来覆盖之前的sce。 因为这些特定函数返回SingleCellExperiment对象，该对象除包含原始数据外还包含结果。 （某些函数-尤其是那些面向单细胞的Bioconductor软件包的函数是不会有这种操作的-在这种情况下，需要将结果附加到sce对象中（请参阅下面的示例。）再次查看该对象，我们看到这些函数添加了一些新条目：
sce
## class: SingleCellExperiment
## dim: 10 3
## metadata(0):
## 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(1): sizeFactor
## reducedDimNames(0):
## altExpNames(0):
具体来说，我们看到“assays”槽已增长为包含两个条目：“counts”（我们的初始数据）和“logcounts”（对数转换的归一化数据）。 与“ counts”类似，“ logcounts”可以使用logcounts（sce）方便地访问，尽管普通版本也可以使用。
logcounts(sce)
## cell_1 cell_2 cell_3
## gene_1 3.841302 4.017667 4.218261
## gene_2 4.502500 4.177661 3.609809
## gene_3 4.700440 4.017667 4.218261
## gene_4 3.662965 4.177661 4.482167
## gene_5 3.841302 4.017667 3.945877
## gene_6 4.142958 4.321670 4.447296
## gene_7 4.273018 4.017667 3.609809
## gene_8 4.604862 4.683435 4.516216
## gene_9 3.841302 4.321670 4.644902
## gene_10 4.392317 4.321670 4.133056
要查看sce中所有可用的assays，我们可以使用assays（）访问。 相比之下，assay（）仅返回单个感兴趣的测定。
assays(sce)
## List of length 2
## names(2): counts logcounts
尽管以上函数会自动添加assays到我们的sce对象中，但在某些情况下，我们可能希望执行自己的计算并将结果保存到assays槽中。 当使用不返回SingleCellExperiment对象的函数时，这通常是必需的。 为了说明这一点，让我们附加一个新版本的数据，该数据已通过向所有值加100来生成。
counts_100 <- counts(sce) + 100
assay(sce, "counts_100") <- counts_100 # assign a new entry to assays slot
assays(sce) # new assay has now been added.
## List of length 3
## names(3): counts logcounts counts_100
处理meta信息数据
1在列上处理
为了进一步注释我们的SingleCellExperiment对象，我们可以添加meta信息来描述主要数据的列，例如实验的样本或细胞。 此数据输入到colData插槽，形成一个data.frame或DataFrame对象，其中每一行对应于单个细胞，列对应于meta信息字段，例如原产批次，处理条件（图1，橙色框）。 现在，我们为细胞提供一些meta信息，从批次变量开始，其中细胞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（）```构造函数从头开始。 现在，我们将从头开始：
sce <- SingleCellExperiment(assays = list(counts = counts_matrix),
colData = 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):
## altExpNames(0):
我们可以使用colData（）函数访问列数据：
colData(sce)
## DataFrame with 3 rows and 1 column
## batch
## <numeric>
## cell_1 1
## cell_2 1
## cell_3 2
更简单地说，我们可以使用$快捷方式提取单个字段：
sce$batch
## [1] 1 1 2
某些函数通过返回额外字段到colData插槽中而自动添加列的meta信息数据。 例如，scater程序包包含addPerCellQC（）函数，该函数附加了许多质量控制数据。 在这里，我们显示colData（sce）的前五列，并附加了质量控制指标。
sce <- scater::addPerCellQC(sce)
colData(sce)
## DataFrame with 3 rows and 4 columns
## batch sum detected total
## <numeric> <numeric> <numeric> <numeric>
## cell_1 1 106 10 106
## cell_2 1 93 10 93
## cell_3 2 331 10 331
另外，我们可能想手动向列的meta信息数据添加更多字段：
sce$more_stuff <- runif(ncol(sce))
colnames(colData(sce))
## [1] "batch" "sum" "detected" "total" "more_stuff"
对colData的常见操作是使用其值进行挑选子集。 例如，如果我们只希望批次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(0):
## colnames(2): cell_1 cell_2
## colData names(5): batch sum detected total more_stuff
## reducedDimNames(0):
## altExpNames(0):
2.在行上处理
为了存储“feature-level”的注释信息，SingleCellExperiment的rowData插槽包含一个DataFrame，其中的每一行都与一个基因相对应，并包含诸如转录本长度或基因符号之类的注释。 此外，还有一个特殊的rowRanges插槽以GRanges或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>
当前，rowData插槽也为空。 但是，类似于上一节中对addPerCellQC（）的调用，addPerFeatureQC（）函数将在sce对象的rowData插槽中插入值：
sce <- scater::addPerFeatureQC(sce)
rowData(sce)
## DataFrame with 10 rows and 2 columns
## mean detected
## <numeric> <numeric>
## gene_1 16.3333 100
## gene_2 14.3333 100
## gene_3 18.6667 100
## gene_4 18.6667 100
## gene_5 14.3333 100
## gene_6 19.6667 100
## gene_7 13.3333 100
## gene_8 22.6667 100
## gene_9 21.0000 100
## gene_10 17.6667 100
类似于colData插槽的方式，可以在开始创建SingleCellExperiment对象时提供此类meta信息数据。 确切的方式取决于对齐和定量过程中可用的生物体和注释。 例如，给定Ensembl标识符，我们可以使用AnnotationHub资源库中下载Ensembl注释对象并提取基因体以存储在SingleCellExperiment的rowRanges中。
library(AnnotationHub)
edb <- AnnotationHub()[["AH73881"]] # Human, Ensembl v97.
genes(edb)[,2]
## GRanges object with 67667 ranges and 1 metadata column:
## seqnames ranges strand | gene_name
## <Rle> <IRanges> <Rle> | <character>
## ENSG00000223972 1 11869-14409 + | DDX11L1
## ENSG00000227232 1 14404-29570 - | WASH7P
## ENSG00000278267 1 17369-17436 - | MIR6859-1
## ENSG00000243485 1 29554-31109 + | MIR1302-2HG
## ENSG00000284332 1 30366-30503 + | MIR1302-2
## ... ... ... ... . ...
## ENSG00000224240 Y 26549425-26549743 + | CYCSP49
## ENSG00000227629 Y 26586642-26591601 - | SLC25A15P1
## ENSG00000237917 Y 26594851-26634652 - | PARP4P1
## ENSG00000231514 Y 26626520-26627159 - | CCNQP2
## ENSG00000235857 Y 56855244-56855488 + | CTBP2P1
## -------
## seqinfo: 424 sequences from GRCh38 genome
要在feature/gene 级别上对SingleCellExperiment对象进行子集化，我们可以通过提供数字索引或名称向量来进行类似于其他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(2): mean detected
## colnames(3): cell_1 cell_2 cell_3
## colData names(5): batch sum detected total more_stuff
## reducedDimNames(0):
## altExpNames(0):
sce[c(1, 4), ] # same as above in this case
## class: SingleCellExperiment
## dim: 2 3
## metadata(0):
## assays(1): counts
## rownames(2): gene_1 gene_4
## rowData names(2): mean detected
## colnames(3): cell_1 cell_2 cell_3
## colData names(5): batch sum detected total more_stuff
## reducedDimNames(0):
## altExpNames(0):
3.其他meta信息数据
一些分析包含的结果和注释不适合存储于上述插槽，例如研究设计本身的meta数据。 值得庆幸的是，有一个插槽专门用于这种混乱数据的存储-metadata插槽，一个命名的列表，列表中的每个条目都可以是您想要的任何内容。 例如，假设我们有一些喜欢的基因（例如，高度可变的基因），我们希望将其存储在sce中，以备后用。 我们可以简单地通过添加到metadata数据插槽来做到这一点，如下所示：
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"
单细胞特有的一些字段
1.背景
到目前为止，我们已经介绍了SingleCellExperiment类的assays（原始数据），colData（细胞注释数据），rowData / rowRanges（feature注释数据）和其他注释数据插槽（metadata slots）。 这些插槽实际上是与SummarizedExperiment父类继承的（有关详细信息，请参见[此处](https://bioconductor.org/packages/3.12/SummarizedExperiment/vignettes/SummarizedExperiment.html），因此，适用于SummarizedExperiment的任何方法也将适用于SingleCellExperiment对象。 但是为什么我们需要一个单独的SingleCellExperiment类？ 这是出于简化某些特定于单细胞数据操作的需求，我们将在本节的其余部分中进行讨论。
2.降维结果
reduceDims插槽专门设计用于存储通过PCA和t-SNE等方法获得的原始数据的降维结果。 该插槽包含低维的原始数据的数字矩阵列表，其中的行表示原始数据的列（即细胞），而列表示维度。 由于此插槽包含一个列表，因此我们可以存储同一数据集的多个PCA /t-SNE 等方式降维的结果。
在下面的示例中，我们可以使用scater中的runPCA（）函数来计算数据的PCA表示形式。 我们看到该sce现在显示了一个新的reduceDim，可以使用reduceDim（）进行访问。
sce <- scater::logNormCounts(sce)
sce <- scater::runPCA(sce)
reducedDim(sce, "PCA")
## PC1 PC2
## cell_1 -0.9220510 0.1440957
## cell_2 0.0821027 -0.3350241
## cell_3 0.8399483 0.1909284
## attr(,"varExplained")
## [1] 0.78121599 0.08472916
## attr(,"percentVar")
## [1] 90.215412 9.784588
## attr(,"rotation")
## PC1 PC2
## gene_2 -0.49697176 -0.36357578
## gene_4 0.46744137 -0.09457483
## gene_9 0.45724971 -0.04428599
## gene_3 -0.29509275 0.80659316
## gene_7 -0.36995311 -0.24239829
## gene_1 0.21191718 0.07604183
## gene_6 0.17299986 -0.01042145
## gene_10 -0.14308909 -0.15243691
## gene_5 0.06548595 -0.23085318
## gene_8 -0.04352653 -0.25521831
我们还可以使用scater包函数runTSNE（）计算原始数据的tSNE表示形式：
sce <- scater::runTSNE(sce, perplexity = 0.1)
## Perplexity should be lower than K!
reducedDim(sce, "TSNE")
## [,1] [,2]
## cell_1 -2518.011 -5105.8218
## cell_2 5683.687 368.9928
## cell_3 -3165.676 4736.8290
我们可以通过reduceDims（）来查看ReducedDims插槽中所有条目的名称。 请注意，这是复数形式，并返回所有结果的列表，而reduceDim（）仅返回单个结果。
reducedDims(sce)
## List of length 2
## names(2): PCA TSNE
我们还可以手动将内容添加到reduceDims（）插槽中，就像我们之前将矩阵添加到assays插槽中的方式一样。 为了说明这一点，我们直接用uwot包运行umap（）函数来生成UMAP坐标矩阵，该矩阵将添加到sce对象的reduceDims中。 （实际上，scater具有一个runUMAP（）包装函数，该函数为我们添加了结果，但是出于演示目的，我们将在此处手动调用umap（）。
u <- uwot::umap(t(logcounts(sce)), n_neighbors = 2)
reducedDim(sce, "UMAP_uwot") <- u
reducedDims(sce) # Now stored in the object.
## List of length 3
## names(3): PCA TSNE UMAP_uwot
reducedDim(sce, "UMAP_uwot")
## [,1] [,2]
## cell_1 -0.4999728 -0.4328734
## cell_2 0.3765395 -0.2061922
## cell_3 0.1234334 0.6390656
## attr(,"scaled:center")
## [1] 3.059862 5.467025
3.alternative实验
SingleCellExperiment类提供“alternative实验”的概念，其中我们具有一组不同的feature但具有相同的样本/细胞的数据。 经典的应用程序是存储每个细胞的spike-in转录本计数。 这使我们可以保留此数据供下游使用，但可以将其与保存内源基因计数的assay分开。 这种分离尤其重要，因为此类feature通常需要单独处理。后续会单独阐述。
如果我们有feature集的数据，则可以将其存储在SingleCellExperiment中作为alternative实验。 例如，如果我们有一些spike-in的数据，我们首先创建一个单独的SummarizedExperiment对象：
spike_counts <- cbind(cell_1 = rpois(5, 10),
cell_2 = rpois(5, 10),
cell_3 = rpois(5, 30))
rownames(spike_counts) <- paste0("spike_", 1:5)
spike_se <- SummarizedExperiment(list(counts=spike_counts))
spike_se
## class: SummarizedExperiment
## dim: 5 3
## metadata(0):
## assays(1): counts
## rownames(5): spike_1 spike_2 spike_3 spike_4 spike_5
## rowData names(0):
## colnames(3): cell_1 cell_2 cell_3
## colData names(0):
然后，我们通过altExp（）将此SummarizedExperiment存储在sce对象中。 像assays（）和reduceDims（）一样，我们还可以使用altExps（）检索所有可用的替代实验。
altExp(sce, "spike") <- spike_se
altExps(sce)
## List of length 1
## names(1): spike
alternative实验概念可确保将单个单细胞数据集的所有相关方面保存在单个对象中。 这也很方便，因为它可以确保我们的spike-in数据与内源基因的数据同步。 例如，如果我们对sce进行了子集化，则spike-in数据将被子集化以匹配：
sub <- sce[,1:2] # retain only two samples.
altExp(sub, "spike")
## class: SummarizedExperiment
## dim: 5 2
## metadata(0):
## assays(1): counts
## rownames(5): spike_1 spike_2 spike_3 spike_4 spike_5
## rowData names(0):
## colnames(2): cell_1 cell_2
## colData names(0):
任何SummarizedExperiment对象都可以存储为替代实验，包括另一个SingleCellExperiment！ 这使高级用户可以进行更为复杂的数据分析。
4.尺寸因素
sizeFactors（）函数使我们能够获取或设置用于归一化的每个细胞缩放因子的数值矢量。 通常由归一化功能自动添加，如下所示，用于基于scran基于反卷积的size因子：
sce <- scran::computeSumFactors(sce)
sizeFactors(sce)
## [1] 0.6000000 0.5264151 1.8735849
另外，我们可以手动添加size factor，如下所示，以获取库大小衍生的因子：
sizeFactors(sce) <- scater::librarySizeFactors(sce)
sizeFactors(sce)
## cell_1 cell_2 cell_3
## 0.6000000 0.5264151 1.8735849
从技术上讲，sizeFactors概念并不是单细胞分析所独有的。 尽管如此，我们还是在这里提到它，因为它是对SummarizedExperiment父类提供的扩展。
5.列标签
colLabels（）函数使我们能够获取或设置每个细胞标签的向量或因子，通常对应于无监督聚类分配的分组或分类算法中预测的细胞类型标识。
colLabels(sce) <- LETTERS[1:3]
colLabels(sce)
## [1] "A" "B" "C"
这是一个方便设置的字段，因为几个函数（例如scran :: findMarkers）将尝试通过colLabels（）自动检索标签。 因此，我们可以避免一些额外的操作去指定集群。
小结：
SingleCellExperiment类的广泛使用为Bioconductor生态系统中与单细胞相关的软件包之间的相互操作性奠定了基础。 一个包生成的SingleCellExperiment对象可以用作另一个包的输入，从而促进了协同工作，使我们的分析大于其各个部分的总和。 分析的每个步骤还将在assays，colData，reduceDims等中添加新条目，这意味着最终的SingleCellExperiment对象有效地充当了分析的独立记录。 这很方便，因为可以将对象保存起来以备将来使用，也可以将其转移给协作者以进行进一步分析。
网友评论