最近在分析单细胞数据，因此和Seurat包打交道了不少。由于seurat的对象设计太复杂了，遂花了一些时间去学习其数据结构以及属性访问方法，整理为笔记特放置于此，谨供需要的朋友们参考。
一、背景知识：R语言的S3,R4,RC,R6对象
另外参考博客往期文章： 《R语言的对象系统（S3-R6）》
R语言脱胎于S语言，因此在设计上继承了许多S语言的特性。近年来，随着编程技术的发展，R语言又吸收了许多新的设计思想，因此在原先的特性的基础上产生了许多扩展。这就是R语言对象系统的复杂性，其按照出现的先后顺序分为S3,R4,RC,R6。 在本文中我们将要详细介绍的seurat就属于R4对象系统。
- S3 对象系统：R语言最原始的对象系统，基于一种特殊的泛型函数实现。使用
attr(obj,"class") 进行对象创建和数据设置，使用 UseMethod() 创建泛型函数。Seurat当中几乎不涉及S3对象系统。
- R4对象系统：这是标准的 R 语言面向对象实现方式，比 S3 的定义更加严格，S4 对象有专门的函数用于定义类（
setClass）、泛型函数（setGeneric）、方法（setMethod）以及实例化对象（new），提供了参数检查，多重继承功能。另外， S4 有一个重要的组件 slot，它是对象的属性组件，可以使用专门的运算符 @来访问（对于Seurat对象，我们还会经常碰到这个运算符）。
- RC对象系统：全名是Reference Classes，是在 R 2.12 版本开始引入的新一代的面向对象系统，通过类似
C++ 语言的风格进行类定义和方法封装。使用 $ 符号来调用方法，获取和修改对象的属性，调用方法或设置属性的值会修改对象，这种方式不同于常用的函数式编程模型。
- R6对象系统：这是一个第三方的 R 面向对象编程的实现，比内置的 RC 类更简单，更快，更轻量级。
二、Seurat Object的结构
使用 str() 函数可以查看一个对象的结构。下面我们以Seurat官方文档的示例数据为例，展示一下其内部结构：
读取数据：
1
2
3
4
5
6
7
8
| library(Seurat)
data_dir = "D:/linux/R/Seurat-demo/filtered_gene_bc_matrices/hg19"
pbmc.data <- Read10X(data.dir = data_dir)
pbmc <- CreateSeuratObject(counts = pbmc.data, project = "pbmc3k", min.cells = 3, min.features = 200)
str(pbmc)
|
下面是使用 str() 展示的数据结构。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
| > str(pbmc)
Formal class 'Seurat' [package "SeuratObject"] with 13 slots
..@ assays :List of 1
.. ..$ RNA:Formal class 'Assay5' [package "SeuratObject"] with 8 slots
.. .. .. ..@ layers :List of 1
.. .. .. .. ..$ counts:Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
.. .. .. .. .. .. ..@ i : int [1:2282976] 29 73 80 148 163 184 186 227 229 230 ...
.. .. .. .. .. .. ..@ p : int [1:2701] 0 779 2131 3260 4220 4741 5522 6304 7094 7626 ...
.. .. .. .. .. .. ..@ Dim : int [1:2] 13714 2700
.. .. .. .. .. .. ..@ Dimnames:List of 2
.. .. .. .. .. .. .. ..$ : NULL
.. .. .. .. .. .. .. ..$ : NULL
.. .. .. .. .. .. ..@ x : num [1:2282976] 1 1 2 1 1 1 1 41 1 1 ...
.. .. .. .. .. .. ..@ factors : list()
.. .. .. ..@ cells :Formal class 'LogMap' [package "SeuratObject"] with 1 slot
.. .. .. .. .. ..@ .Data: logi [1:2700, 1] TRUE TRUE TRUE TRUE TRUE TRUE ...
.. .. .. .. .. .. ..- attr(*, "dimnames")=List of 2
.. .. .. .. .. .. .. ..$ : chr [1:2700] "AAACATACAACCAC-1" "AAACATTGAGCTAC-1" "AAACATTGATCAGC-1" "AAACCGTGCTTCCG-1" ...
.. .. .. .. .. .. .. ..$ : chr "counts"
.. .. .. .. .. ..$ dim : int [1:2] 2700 1
.. .. .. .. .. ..$ dimnames:List of 2
.. .. .. .. .. .. ..$ : chr [1:2700] "AAACATACAACCAC-1" "AAACATTGAGCTAC-1" "AAACATTGATCAGC-1" "AAACCGTGCTTCCG-1" ...
.. .. .. .. .. .. ..$ : chr "counts"
.. .. .. ..@ features :Formal class 'LogMap' [package "SeuratObject"] with 1 slot
.. .. .. .. .. ..@ .Data: logi [1:13714, 1] TRUE TRUE TRUE TRUE TRUE TRUE ...
.. .. .. .. .. .. ..- attr(*, "dimnames")=List of 2
.. .. .. .. .. .. .. ..$ : chr [1:13714] "AL627309.1" "AP006222.2" "RP11-206L10.2" "RP11-206L10.9" ...
.. .. .. .. .. .. .. ..$ : chr "counts"
.. .. .. .. .. ..$ dim : int [1:2] 13714 1
.. .. .. .. .. ..$ dimnames:List of 2
.. .. .. .. .. .. ..$ : chr [1:13714] "AL627309.1" "AP006222.2" "RP11-206L10.2" "RP11-206L10.9" ...
.. .. .. .. .. .. ..$ : chr "counts"
.. .. .. ..@ default : int 1
.. .. .. ..@ assay.orig: chr(0)
.. .. .. ..@ meta.data :'data.frame': 13714 obs. of 0 variables
.. .. .. ..@ misc :List of 1
.. .. .. .. ..$ calcN: logi TRUE
.. .. .. ..@ key : chr "rna_"
..@ meta.data :'data.frame': 2700 obs. of 3 variables:
.. ..$ orig.ident : Factor w/ 1 level "pbmc3k": 1 1 1 1 1 1 1 1 1 1 ...
.. ..$ nCount_RNA : num [1:2700] 2419 4903 3147 2639 980 ...
.. ..$ nFeature_RNA: int [1:2700] 779 1352 1129 960 521 781 782 790 532 550 ...
..@ active.assay: chr "RNA"
..@ active.ident: Factor w/ 1 level "pbmc3k": 1 1 1 1 1 1 1 1 1 1 ...
.. ..- attr(*, "names")= chr [1:2700] "AAACATACAACCAC-1" "AAACATTGAGCTAC-1" "AAACATTGATCAGC-1" "AAACCGTGCTTCCG-1" ...
..@ graphs : list()
..@ neighbors : list()
..@ reductions : list()
..@ images : list()
..@ project.name: chr "pbmc3k"
..@ misc : list()
..@ version :Classes 'package_version', 'numeric_version' hidden list of 1
.. ..$ : int [1:3] 5 0 2
..@ commands : list()
..@ tools : list()
>
|
我们可以看到，在pbmc对象的第一层，是几个数据槽，包括@assays , @meta.data 等。 其中 @assays 槽是最为关键的数据槽，其中存储了表达矩阵的信息（即 $RNA ）。
另外，毋庸置疑的，当我们继续运行seurat的分析流程，对pbmc对象进行更深层次的分析（如归一化、降维聚类、细胞分群与注释）以后， str(pbmc) 给出的列表还会继续增长：因为Seurat的R4对象会将一切数据都存储在对象结构内部。
那么，如何访问呢？这是下一小节的内容。
三、Seurat对象内部属性的访问方法
3.1 使用R4数据槽操作符 @ 进行访问
这是一种很直观但并 不 被官方推荐的方法，原因在于Seurat在不断更新，也许某一个数据槽在新版本Seurat中就会改叫其他的名字；另外，这种方式过于底层，容易破坏对象的完整性和一致性。如果直接修改槽内容，可能导致不可预期的行为。
不过，我们还是在这里介绍一下访问方法。如下：
1
2
3
4
|
metadata = pbmc@meta.data
class(metadata)
metadata[1:10,1:3]
|
输出：
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
| > metadata = pbmc@meta.data
> class(metadata)
[1] "data.frame"
> metadata[1:10,1:3]
orig.ident nCount_RNA nFeature_RNA
AAACATACAACCAC-1 pbmc3k 2419 779
AAACATTGAGCTAC-1 pbmc3k 4903 1352
AAACATTGATCAGC-1 pbmc3k 3147 1129
AAACCGTGCTTCCG-1 pbmc3k 2639 960
AAACCGTGTATGCG-1 pbmc3k 980 521
AAACGCACTGGTAC-1 pbmc3k 2163 781
AAACGCTGACCAGT-1 pbmc3k 2175 782
AAACGCTGGTTCTT-1 pbmc3k 2260 790
AAACGCTGTAGCCA-1 pbmc3k 1275 532
AAACGCTGTTTCTG-1 pbmc3k 1103 550
>
|
如果我们想要看表达矩阵的信息，则可以像下面这样（虽然比较麻烦）
1
2
3
4
5
|
RNA = pbmc@assays$RNA
counts = RNA@layers$counts
counts [1:15,1:15]
|
输出：
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
| > counts [1:15,1:15]
15 x 15 sparse Matrix of class "dgCMatrix"
[1,] . . . . . . . . . . . . . . .
[2,] . . . . . . . . . . . . . . .
[3,] . . . . . . . . . . . . . . .
[4,] . . . . . . . . . . . . . . .
[5,] . . . . . . . . . . . . . . .
[6,] . . . . . . . . . . . 1 . . .
[7,] . . . . . . . . . . . . . . .
[8,] . . . . . . . . . . . . . . .
[9,] . . . . . . . . . . . . . . .
[10,] . . . . . . . . . . . . . . .
[11,] . . . . . . . . . . . . . . .
[12,] . . 1 9 . 1 . . . 3 . . 1 5 .
[13,] . . . . . . . . . . . . . . .
[14,] . . . . . . . . . . . . . . .
[15,] . 2 . . . . . . . . . . . . .
>
|
3.2 使用GetAssayData()或LayerData()更安全的访问数据
GetAssayData() 和 LayerData() 是Seurat对数据访问接口的封装。官方也更推荐用这个进行访问。
二者其实是同一个函数，只不过一个是在Seurat V4当中的名称（GetAssayData），一个是在Seurat V5当中的名称（ LayerData ），使用前需要检查Seurat的版本，然后根据版本使用正确的函数名。
笔者电脑上的Seurat 版本为V5，因此下面的演示以LayerData()为例。下面的代码块展示了如何从一个 seurat_obj 中获得原始表达矩阵和归一化的表达矩阵。
1
2
| counts <- LayerData(seurat_obj, assay = "RNA", layer = "counts")
normalized_data <- LayerData(seurat_obj, assay = "SCT", layer = "data")
|
例如，上述pbmc的数据：（注意，SCTransformer矫正前，已经根据基因表达量过滤了一部分细胞，因此原始表达矩阵counts和归一化以后的表达矩阵normalized_data的行和列并非完全一一对应。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
| > counts [1:15,1:15]
15 x 15 sparse Matrix of class "dgCMatrix"
[[ suppressing 15 column names ‘AAACATACAACCAC-1’, ‘AAACATTGAGCTAC-1’, ‘AAACATTGATCAGC-1’ ... ]]
AL627309.1 . . . . . . . . . . . . . . .
AP006222.2 . . . . . . . . . . . . . . .
RP11-206L10.2 . . . . . . . . . . . . . . .
RP11-206L10.9 . . . . . . . . . . . . . . .
LINC00115 . . . . . . . . . . . . . . .
NOC2L . . . . . . . . . . . 1 . . .
KLHL17 . . . . . . . . . . . . . . .
PLEKHN1 . . . . . . . . . . . . . . .
RP11-54O7.17 . . . . . . . . . . . . . . .
HES4 . . . . . . . . . . . . . . .
RP11-54O7.11 . . . . . . . . . . . . . . .
ISG15 . . 1 9 . 1 . . . 3 . . 1 5 .
AGRN . . . . . . . . . . . . . . .
C1orf159 . . . . . . . . . . . . . . .
TNFRSF18 . 2 . . . . . . . . . . . . .
> normalized_data [1:15,1:5]
15 x 5 sparse Matrix of class "dgCMatrix"
AAACATACAACCAC-1 AAACATTGAGCTAC-1 AAACATTGATCAGC-1 AAACCGTGCTTCCG-1 AAACCGTGTATGCG-1
AL627309.1 . . . . .
RP11-206L10.2 . . . . .
LINC00115 . . . . .
NOC2L . . . . .
KLHL17 . . . . .
PLEKHN1 . . . . .
HES4 . . . . .
ISG15 . . 0.6931472 2.197225 .
AGRN . . . . .
C1orf159 . . . . .
TNFRSF18 . 0.6931472 . . .
TNFRSF4 . . . . .
SDF4 . . 0.6931472 . .
B3GALT6 . . . . .
UBE2J2 . . . . .
>
|
3.3 使用 [[ 访问嵌套列表中的元素
- 作用：用于访问 Reductions、Graphs 等嵌套列表中的具体对象。
1
2
| seurat_obj[["reductions"]][["pca"]]
seurat_obj[["graphs"]][["SNN"]]
|
- 适用场景：当需要直接操作嵌套列表时使用（但推荐用封装函数如
Embeddings()）。
- 注意：这种方式较为灵活，但也容易出错，建议结合封装函数使用。
四、一点补充说明：Seurat对象中的四个表达量矩阵
这一部分是我在工作中的发现。
跟着Seurat官方教程跑完了整个pipeline（包括归一化、聚类、分群），检查数据时发现数据当中存在多个表达矩阵（如下表第一列）。
查询各种资料以及文档以后，了解到了他们的区别如下：
| 访问方法 |
数据说明 |
GetAssayData(pbmc, assay = "RNA", slot = "counts") |
counts（原始UMI） |
GetAssayData(pbmc, assay = "SCT", slot = "counts") |
counts（过滤后） |
GetAssayData(pbmc, assay = "SCT", slot = "data") |
SCTransform归一化数据data |
GetAssayData(pbmc, assay = "SCT", slot = "scale.data") |
经过z-score标准化之后的SCTransform归一化数据data |
如果要访问原始矩阵，推荐使用 GetAssayData(pbmc, assay = "RNA", slot = "counts") 。如果要查询SCTransform归一化并经过z-score标准化以后的数据，则需要使用 GetAssayData(pbmc, assay = "SCT", slot = "scale.data") 。
以上。