重生之我在剑桥大学学习单细胞RNA-seq分析——4. 使用Bioconductor进行scRNA-seq分析（2）

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(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）

重生之我在剑桥大学学习单细胞RNA-seq分析——4. 使用Bioconductor进行scRNA-seq分析（2）

推荐阅读更多精彩内容