Seurat 与 Cellranger 之间互通的二三事

最近，在做单细胞测序的分析，出现了这么一个需求：Cellranger 中没有像 Seurat 一样进行单细胞数据中常见的几类质控，比如 nGene，nUMI, percent of mitochondria genes 等，因此对于 cellranger 得到的矩阵先要经过这类质控，再进行 cellranger analyze 的后续分析。本来 Cellranger 自己有进行后续的 Rkit 包进行分析，但由于第三方包的发展，这个包已经被废弃了。鉴于 Cellranger 主页自己推荐使用 Seurat 进行后续分析，我们也就选择了这个软件进行质控及后续分析。

问题来了

在 Seurat 中，对于 Cellranger 数据的导入，它提供了两个函数 Read10X 和 Read10X_h5。前者用于读取 Cellranger 生成的 3 个文件：barcodes.tsv，features.tsv，matrix.mtx，后者则用于读取生成的 xxxxxx.h5文件。本来说前一个函数读入来其实就很方便了，但是由于 Cellranger 的 analyze 子程序的输入文件是一个 .h5 文件，所以直接读入 xxxxxx.h5 文件似乎更好一些。然而，第一个问题就来了。Read10X_h5 函数读取时，会报这样一个错误：

1

Data has no data["matrix/gene_name"]

查看源码发现，gene_name 这个数据在 Read10X_h5 中是直接从 data[["matrix"]] 中读取的，但是，如果你仔细看过的 cellranger 官网上有关 h5 数据的介绍的话，它是没有 gene_name 这一个数据项的，而是保存在 data[["matrix/features/name"]] 中，不知道 Seurat 的这个函数是只兼容到 cellranger 3.0 之前的结果还是怎么的。既然如此，我只好自己加上了。比如这样：

1
2
3
4

require(hdf5r)
infile <- H5File$new(filename, "r+")
infile[["matrix/gene_name"]] <- infile[["matrix/features/name"]][]
infile$close_all()

到这样一切顺利，可以使用 Read10X_h5 了。紧接着进行那些质控什么的，筛选一下数据啊，生成一些图片啊，高大上，很令人兴奋有没有？然后赶紧保存为 .h5 格式的文件给 cellranger 使用，可是这一切并没有那么顺利，第二个问题这么快就出现了。

查了半天文档，WTF，Seurat 并没有给函数把它的对象导出去啊，更别提什么导出为 h5 格式的文件。唉，不对，应该有的。仔细查查罢。Bingo，有个 Convert 函数，它可以导出为 loom 格式的文件，号称是格式更为严格的 hdf5 类型。然而这似乎并没有什么卵用。因为 Seurat 连读入 cellranger 的 h5 文件都没兼容好，就这个转为 loom 文件能兼容好？索性放弃之，不能乱踩坑的（666）。那就自己来吧。

基于 Seurat 对象构建 .h5 文件

这里就不赘述其中的踩坑过程了，反正就是写 Bug 然后 Debug 的过程，我只把关键部分放在这里。前面我们提到，我们看了 Read10X_h5 的源码对不对，这样就肯定知道 .h5 文件是怎么变为 Seurat 对象的。首先是读取文件变成原始数据，然后把原始数据传给 CreateSeuratObject 构建 Seurat 对象。所以关键的部分就是传给这个函数的原始数据了。我们可以看一下这个原始数据的结构。

请注意，这个原始数据是个 S4 object of class dgTMatrix，也就是一个稀疏矩阵的一种表示法，关于这个，请看另外一篇笔记稀疏矩阵的表示法。这里需要关注的点就是：

• i ：读入 .h5 文件的 “matrix/indices” 数据集生成，即生成的稀疏矩阵表示法中非 0 值的 row index
• J ：读入 .h5 文件后，生成的稀疏矩阵表示法中非 0 值的 colname index
• Dim：原始矩阵的维度
• x：原始矩阵中的所有非零值的集合，值的顺序是原始矩阵中从左到右，从上到下排列
• Dimnames：原始矩阵的行列名，也就是 barcodes 和 gene names

但是我们可以看一下 cellranger 需要的数据内容：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17

(root)
└── matrix [HDF5 group]
    ├── barcodes
    ├── data
    ├── indices
    ├── indptr
    ├── data
    ├── shape
    └── features [HDF5 group]
        ├─ _all_tag_keys
        ├─ feature_type
        ├─ genome
        ├─ id
        ├─ name
        ├─ pattern [Feature Barcoding only]
        ├─ read [Feature Barcoding only]
        └─ sequence [Feature Barcoding only]

最主要的就是这个：

我们可以得到：data == x， barcodes == Dimnames[2]，indices == i，shape == Dim。然后只剩下 indptr 了。在开始的时候，并不知道怎么去获得 indptr 这个参数，总不能自己去写代码搞吧。然后仔细读 cellranger 另一个有关数据格式的页面中发现，cellranger 使用的是 Compressed sparse column format。基于此，就好办一些了，使用下面这个函数把原来的 dgTMatrix 转换为 dgCMatrix。

1

scrna.raw2 <- as(scrna.raw, "dgCMatrix")

转换后的 scrna.raw2 就不包含前面所说的 j 数据，而是多了一个 p 数据，它就是 indptr.

得到这些数据我们就可以利用它们进行生成 .h5 文件了。但是，经过我的踩坑，不推荐自己从头创建出新的 .h5 文件，而是基于 cellranger 的 .h5 文件进行修改。这样做的原因是，cellranger.h5 文件的数据集的 dataType 同 hdf5r 默认创建的 dataType 在具体的数据格式参数上有区别，但是 hdf5r 并没有提供 API 进行修改，你要修改的话，只能在更底层的层面进行（hdf5r 只是提供了接口来调用其他语言的包来交互 hdf5 文件）。如此做太费劲，不值当。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24

require(hdf5r)
cellranger.h5 <- "/path/to/cellranger.h5"
filename <- "my_seurat.h5"
file.copy(cellranger.h5, filename)
newinfile <- H5File$new(filename, "r+")	# 修改模式

mtx <- newinfile[["matrix"]]
## 要覆盖原值，先删去原始值
for (property in c("data", "indices", "indptr", "shape", "barcodes")) {
	mtx$link_delete(property)
}
mtx[["data"]] <- scrna.raw2@x
mtx[["indices"]] <- scrna.raw2@i
mtx[["indptr"]] <- scrna.raw2@p
mtx[["shape"]] <- scrna.raw2@Dim
mtx[["barcodes"]] <- scrna.raw2@Dimnames[2]
mtx[["gene_names"]] <- scrna.raw2@Dimnames[1]  # 修复 Read10X_h5 的 Bug

## 这里的操作方式同前面一样
## 只是各项的值，你需要根据前面的 barcodes、gene_names 的内容和长度来相对应
features <- mtx[["features"]]
## 省略了修改 features 的好多代码

newinfile.close_all()