4.2.3 矩阵统计
因为SingleCellExperiment对象中存储的主要数据是一个矩阵，所以覆盖矩阵行或列间汇总指标的函数很有用。有几个函数可以做到这一点，详细信息请参见下面的信息框。
例如，要计算数据集中每个细胞（列）的平均计数：
> colMeans(counts(tung))
我们可以将此信息作为新列添加到我们的colData中：
> colData(tung)$mean_counts <- colMeans(counts(tung))
如果我们查看colData，我们可以看到已添加新列：
> colData(tung)
DataFrame with 864 rows and 6 columns
individual replicate well batch sample_id
<character> <character> <character> <character> <character>
NA19098.r1.A01 NA19098 r1 A01 NA19098.r1 NA19098.r1.A01
NA19098.r1.A02 NA19098 r1 A02 NA19098.r1 NA19098.r1.A02
NA19098.r1.A03 NA19098 r1 A03 NA19098.r1 NA19098.r1.A03
NA19098.r1.A04 NA19098 r1 A04 NA19098.r1 NA19098.r1.A04
NA19098.r1.A05 NA19098 r1 A05 NA19098.r1 NA19098.r1.A05
... ... ... ... ... ...
NA19239.r3.H08 NA19239 r3 H08 NA19239.r3 NA19239.r3.H08
NA19239.r3.H09 NA19239 r3 H09 NA19239.r3 NA19239.r3.H09
NA19239.r3.H10 NA19239 r3 H10 NA19239.r3 NA19239.r3.H10
NA19239.r3.H11 NA19239 r3 H11 NA19239.r3 NA19239.r3.H11
NA19239.r3.H12 NA19239 r3 H12 NA19239.r3 NA19239.r3.H12
mean_counts
<numeric>
NA19098.r1.A01 3.32801
NA19098.r1.A02 3.36238
NA19098.r1.A03 2.29306
NA19098.r1.A04 2.83397
NA19098.r1.A05 3.72560
... ...
NA19239.r3.H08 3.097861
NA19239.r3.H09 2.933568
NA19239.r3.H10 0.316813
NA19239.r3.H11 2.678615
NA19239.r3.H12 4.068482
有一些函数可用于计算矩阵（或稀疏矩阵）的行或列的汇总指标——例如平均值、中位数、方差等。
例如：
#row (feature) summaries
rowSums(counts(tung)) # sum
rowMeans(counts(tung)) # mean
rowSds(counts(tung)) # standard deviation
rowVars(counts(tung)) # variance
rowIQRs(counts(tung)) # inter-quartile range
rowMads(counts(tung)) # mean absolute deviation
#column (sample) summaries
colSums(counts(tung)) # sum
colMeans(counts(tung)) # mean
colSds(counts(tung)) # standard deviation
colVars(counts(tung)) # variance
colIQRs(counts(tung)) # inter-quartile range
colMads(counts(tung)) # mean absolute deviation
4.2.4 SCE对象取子集
与R中的data.frame和矩阵对象类似，我们可以使用[运算符通过行（基因）或列（细胞）对我们的SingleCellExperiment进行取子集。一般语法是：sce[rows_of_interest, columns_of_interest]。
例如：
> tung[1:3, ] # the first 3 genes, keep all cells
> tung[, 1:3] # the first 3 cells, keep all genes
> tung[1:3, 1:2] # the first 3 genes and first 2 cells
# subset by name
> tung[c("ENSG00000069712", "ENSG00000237763"), ]
> tung[, c("NA19098.r1.A01", "NA19098.r1.A03")]
> tung[c("ENSG00000069712", "ENSG00000237763"), c("NA19098.r1.A01", "NA19098.r1.A03")]
虽然手动对对象进行取子集有时很有用，但更多时候我们希望根据TRUE/FALSE逻辑语句进行条件取子集。这对于过滤我们的数据非常有用。让我们看一些实际的例子。
假设我们想要保留平均计数大于0.01的基因。要计算每个基因的平均计数（SCE对象中的行），可以使用rowMeans()函数：
> gene_means <- rowMeans(counts(tung))
> gene_means[1:10]
ENSG00000237683 ENSG00000187634 ENSG00000188976 ENSG00000187961 ENSG00000187583
0.28240741 0.03009259 2.63888889 0.23842593 0.01157407
ENSG00000187642 ENSG00000188290 ENSG00000187608 ENSG00000188157 ENSG00000237330
0.01273148 0.02430556 1.89583333 3.32523148 0.02199074
使用逻辑运算符将其转换为TRUE/FALSE向量：
> gene_means[1:10] > 0.01
ENSG00000237683 ENSG00000187634 ENSG00000188976 ENSG00000187961 ENSG00000187583
TRUE TRUE TRUE TRUE TRUE
ENSG00000187642 ENSG00000188290 ENSG00000187608 ENSG00000188157 ENSG00000237330
TRUE TRUE TRUE TRUE TRUE
在[中使用这样的逻辑向量来过滤我们的数据，它将只返回值为TRUE的情况：
> tung[gene_means > 0.01, ]
class: SingleCellExperiment
dim: 16052 864
metadata(0):
assays(2): counts logcounts
rownames(16052): ENSG00000237683 ENSG00000187634 ... ERCC-00170
ERCC-00171
rowData names(0):
colnames(864): NA19098.r1.A01 NA19098.r1.A02 ... NA19239.r3.H11
NA19239.r3.H12
colData names(6): individual replicate ... sample_id mean_counts
reducedDimNames(0):
altExpNames(0):
这样生成的SCE对象比原始对象具有更少的基因。
另一个常见例子是保留基因表达量超过一定阈值的细胞。对于这个问题，我们需要将问题分解成几个部分。
首先让我们检查计数矩阵中哪些基因的表达超过了某个阈值：
> counts(tung) > 0
NA19098.r1.A01 NA19098.r1.A02 NA19098.r1.A03 NA19098.r1.A04
ENSG00000237683 FALSE FALSE FALSE TRUE
ENSG00000187634 FALSE FALSE FALSE FALSE
ENSG00000188976 TRUE TRUE TRUE TRUE
ENSG00000187961 FALSE FALSE FALSE FALSE
ENSG00000187583 FALSE FALSE FALSE FALSE
ENSG00000187642 FALSE FALSE FALSE FALSE
ENSG00000188290 FALSE FALSE FALSE FALSE
ENSG00000187608 TRUE TRUE FALSE TRUE
ENSG00000188157 TRUE TRUE TRUE TRUE
ENSG00000237330 FALSE FALSE FALSE FALSE
可以看到矩阵现在仅由TRUE/FALSE值组成。由于TRUE/FALSE被编码为1/0，我们可以使用colSums()来计算每个细胞中高于此阈值的基因总数：
> total_detected_per_cell <- colSums(counts(tung) > 0)
> total_detected_per_cell[1:10]
NA19098.r1.A01 NA19098.r1.A02 NA19098.r1.A03 NA19098.r1.A04 NA19098.r1.A05
8368 8234 7289 7985 8619
NA19098.r1.A06 NA19098.r1.A07 NA19098.r1.A08 NA19098.r1.A09 NA19098.r1.A10
8659 8054 9429 8243 9407
最后，将这个向量应用到最终条件中，例如我们想要至少有5000个检测到的基因的细胞：
> tung[, total_detected_per_cell > 5000]
class: SingleCellExperiment
dim: 19027 849
metadata(0):
assays(2): counts logcounts
rownames(19027): ENSG00000237683 ENSG00000187634 ... ERCC-00170
ERCC-00171
rowData names(0):
colnames(849): NA19098.r1.A01 NA19098.r1.A02 ... NA19239.r3.H11
NA19239.r3.H12
colData names(6): individual replicate ... sample_id mean_counts
reducedDimNames(0):
altExpNames(0):
这里新的SCE对象比原来的SCE对象拥有更少的细胞。
以下是一些常见过滤方式的语法摘要：
4.3 数据可视化
有多种方法可以将数据进行可视化。我们将使用ggplot2包和Bioconductor包scater来绘制图表，该包有一些辅助函数，可用于从SingleCellExperiment对象中检索可用于可视化的数据。
ggplot的基本组成部分是：
- 包含要绘制的数据的data.frame
- 将映射到图形的不同样式（例如轴、颜色、形状等）的变量（数据框的列）
- 将在图表上绘制的几何图形（例如点、线、箱线图、小提琴图等）
这转化为以下基本语法：
ggplot(data = <data.frame>,
mapping = aes(x = <column of data.frame>, y = <column of data.frame>)) +
geom_<type of geometry>()
例如，让我们直观地看一下每个批次中每个细胞的总数分布情况。该信息存储在colData中，因此我们需要从对象中提取它并将其转换为标准data.frame：
> colData(tung)$total_counts <- colSums(counts(tung))
> cell_info <- as.data.frame(colData(tung))
> head(cell_info)
individual replicate well batch sample_id mean_counts
NA19098.r1.A01 NA19098 r1 A01 NA19098.r1 NA19098.r1.A01 3.328008
NA19098.r1.A02 NA19098 r1 A02 NA19098.r1 NA19098.r1.A02 3.362380
NA19098.r1.A03 NA19098 r1 A03 NA19098.r1 NA19098.r1.A03 2.293057
NA19098.r1.A04 NA19098 r1 A04 NA19098.r1 NA19098.r1.A04 2.833973
NA19098.r1.A05 NA19098 r1 A05 NA19098.r1 NA19098.r1.A05 3.725600
NA19098.r1.A06 NA19098 r1 A06 NA19098.r1 NA19098.r1.A06 3.576549
total_counts
NA19098.r1.A01 63322
NA19098.r1.A02 63976
NA19098.r1.A03 43630
NA19098.r1.A04 53922
NA19098.r1.A05 70887
NA19098.r1.A06 68051
现在我们准备开始作图：
> library(ggplot2)
> ggplot(data = cell_info, aes(x = batch, y = total_counts)) +
geom_violin(fill = 'brown') + theme_bw() +
theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1))
如果我们想要直观地看到每个批次中特定基因的表达分布该怎么办？现在情况变得有点复杂，因为基因表达信息存储在我们的SCE的counts中，而批次信息存储在colData中。为了将这两条信息整合在一起，我们需要进行大量的数据操作，将其全部放入一个数据框中。这是scater包非常有用的地方，因为它提供了ggcells()函数，让我们可以为图表指定所有这些信息。
例如，可以直接从tung SCE对象绘制与上述相同的图：<
> library(scater)
> ggcells(tung, aes(x = batch, y = total_counts)) +
geom_violin(fill = 'orange') + theme_bw() +
theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1))
如果我们想绘制某个基因的表达情况，可以这样做：
> ggcells(tung, aes(x = batch, y = ENSG00000198938), exprs_values = "logcounts") +
geom_violin(fill = 'coral2') + theme_bw() +
theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1))
这里指定了我们想要用于表达值的assay来源（exprs_values）。默认值为“logcounts”，因此在这种情况下我们不必指定它，但如果您想从不同的assay中可视化表达，则需要了解这一点。scater包提供的功能远不止绘图，它还包括用于质量控制的函数，我们将在下一章中介绍。
4.4 总结
- SingleCellExperiment（SCE）对象用于存储表达数据以及与细胞（列）和基因（行）有关的信息。
- 要创建新的SCE对象，我们可以使用
SingleCellExperiment()函数。要读取cellranger的输出，我们可以使用专用函数DropletUtils::read10xCounts()。 - 该对象主要包括：
- assay —— 一个或多个表达定量矩阵
- “counts”包含所有其他assay所基于的原始计数。
- rowData —— 与基因有关的信息。
- colData —— 与细胞有关的信息。
- reducedDim —— 一个或多个降维信息。
- assay —— 一个或多个表达定量矩阵
- 可以使用同名函数访问此对象的所有部分。例如，
assay(sce, "counts")从对象中检索计数矩阵。 - 可以使用赋值运算符
<-添加/修改此对象的各个部分。例如，assay(sce, "logcounts") <- log2(counts(sce) + 1)将向对象添加一个名为“logcounts”的新assay。 - 矩阵汇总指标对于探索表达数据的属性非常有用。一些更有用的函数包括
rowSums()/colSums()和rowMeans()/colMeans()。这些函数可用于汇总assay的信息，例如colSums(counts(sce))会提供每个细胞（列）中的总计数。 - 将矩阵汇总与条件运算符（
>、<、==、!=）相结合可用于使用[进行条件取子集。 - 可以使用
ggcells()函数（来自scater包）直接从SCE对象生成ggplot样式的图。
往期内容：
重生之我在剑桥大学学习单细胞RNA-seq分析——3. 单细胞分析中的R/Bioconductor简介（3）
重生之我在剑桥大学学习单细胞RNA-seq分析——3. 单细胞分析中的R/Bioconductor简介（4）
重生之我在剑桥大学学习单细胞RNA-seq分析——4. 使用Bioconductor进行scRNA-seq分析（1）
