实验记录5：.h5文件解析为.tsv和.mtx格式

2019-01-16 本文已影响0人 MC学公卫

这篇主要针对Seurat包装不上hdf5r包的bug而无法使用Read10X_h5这个函数的读者。

在使用Seurat包读取10Xgenomics测序数据时，我发现github上有很多人面临着与我同样的问题：hdf5r包无法安装，因此在调用Read10X_h5函数时会出现以下报错信息：

Error in Read10X_h5("/bone_marrow/ica_bone_marrow_h5.h5") : 
  Please install the 'hdf5r' package

在用包的时候，若数据格式不是.h5则用不到这个函数。但是偏巧，当你需要用Seurat分析这个数据，而数据格式正好是.h5怎么办呢？

因此我们需要先将.h5解析为Read10X这个函数能够读取的格式。

目标：将.h5文件转化为10X outs中的barcodes.tsv, genes.tsv, matrix.mtx三个文件

文件	格式
barcodes.tsv	1列，为barcode名
genes.tsv	2列，第1列为ENS编号，第2列为基因名
matrix.mtx	3列，第1列为基因编号，第2列为细胞编号，第3列为对应的reads数

梗概：

这是大文件的文本处理，要注意内存和磁盘空间的控制
主要是用R的hdf5r包处理，Seurat的依赖包
tsv文件与csv文件类似，不同点在于分隔符是\t而不是逗号
mtx格式是一个矩阵

Step1：需要一台能装上hdf5r包的设备

本人MacOS系统
R version 3.4.4 (2018-03-15) -- "Someone to Lean On"
Copyright (C) 2018 The R Foundation for Statistical Computing
Platform: x86_64-apple-darwin15.6.0 (64-bit)
能够成功安装hdf5r包。

Step2：查看数据

setwd("/data/immune_census/bone_marrow/data")
library(hdf5r)
bm <- H5File$new("./ica_bone_marrow_h5.h5", mode = "r")
bm
Class: H5File
Filename: /immune_census/bone_marrow/data/ica_bone_marrow_h5.h5
Access type: H5F_ACC_RDONLY
Attributes: TITLE, CLASS, VERSION, FILTERS, PYTABLES_FORMAT_VERSION
Listing:
   name  obj_type dataset.dims dataset.type_class
 GRCh38 H5I_GROUP         <NA>               <NA>

bm$ls()
name     link.type  obj_type num_attrs group.nlinks group.mounted dataset.rank
1 GRCh38 H5L_TYPE_HARD H5I_GROUP         4            7             0           NA
  dataset.dims dataset.maxdims dataset.type_class dataset.space_class committed_type
1         <NA>            <NA>               <NA>                <NA>           <NA>

将GRCh38导出来查看后发现共7组信息。

GRCh38 <- bm[["GRCh38"]]
GRCh38$ls()
        name     link.type    obj_type num_attrs group.nlinks group.mounted dataset.rank dataset.dims
1   barcodes H5L_TYPE_HARD H5I_DATASET         3           NA            NA            1       378000
2       data H5L_TYPE_HARD H5I_DATASET         3           NA            NA            1    303795547
3 gene_names H5L_TYPE_HARD H5I_DATASET         3           NA            NA            1        33694
4      genes H5L_TYPE_HARD H5I_DATASET         3           NA            NA            1        33694
5    indices H5L_TYPE_HARD H5I_DATASET         3           NA            NA            1    303795547
6     indptr H5L_TYPE_HARD H5I_DATASET         3           NA            NA            1       378001
7      shape H5L_TYPE_HARD H5I_DATASET         3           NA            NA            1            2
  dataset.maxdims dataset.type_class dataset.space_class committed_type
1          378000         H5T_STRING          H5S_SIMPLE           <NA>
2       303795547        H5T_INTEGER          H5S_SIMPLE           <NA>
3           33694         H5T_STRING          H5S_SIMPLE           <NA>
4           33694         H5T_STRING          H5S_SIMPLE           <NA>
5       303795547        H5T_INTEGER          H5S_SIMPLE           <NA>
6          378001        H5T_INTEGER          H5S_SIMPLE           <NA>
7           16384        H5T_INTEGER          H5S_SIMPLE           <NA>

将所有的信息导出来

### 1细胞barcode：
barcodes <- GRCh38[["barcodes"]]
barcodes[1:5]
[1] "MantonBM1_HiSeq_1-AAACCTGAGCAGGTCA-1" "MantonBM1_HiSeq_1-AAACCTGCACACTGCG-1"
[3] "MantonBM1_HiSeq_1-AAACCTGCACCGGAAA-1" "MantonBM1_HiSeq_1-AAACCTGCATAGACTC-1"
[5] "MantonBM1_HiSeq_1-AAACCTGCATCGATGT-1"

### 2基因ENS编号
genes <- GRCh38[["genes"]]
genes[1:10]
 [1] "ENSG00000243485" "ENSG00000237613" "ENSG00000186092" "ENSG00000238009" "ENSG00000239945"
 [6] "ENSG00000239906" "ENSG00000241599" "ENSG00000279928" "ENSG00000279457" "ENSG00000228463"

### 3基因名
gene_names <- GRCh38[["gene_names"]]
gene_names[1:10]
 [1] "RP11-34P13.3"  "FAM138A"       "OR4F5"         "RP11-34P13.7"  "RP11-34P13.8" 
 [6] "RP11-34P13.14" "RP11-34P13.9"  "FO538757.3"    "FO538757.2"    "AP006222.2"   

### 4 基因索引号
indices <- GRCh38[["indices"]]
indices[1:10]
[1] 33664 33662 33661 33659 33658 33657 33655 33654 33653 33652

### 5 reads数：
data <- GRCh38[["data"]]
data[1:10]
[1]  5  5 14 12 35  9 15 26 17  5

### 6 
indptr <- GRCh38[["indptr"]]
indptr[1:10]
[1]    0  785 1466 2010 2942 4605 5569 5678 6374 7096

### 7 矩阵shape（33694个基因，378000个细胞）
shape <- GRCh38[["shape"]]
shape[1:2]
[1]  33694 378000

Step3：转化为R可操作的dataframe

现在仍然是h5的compounds数据类型，要把它转化为R可操作的dataframe

### 读取数据                                   ###数据个数（用length查看）
barcodes <- bm$open("/GRCh38/barcodes")       #378000
genes <- bm$open("/GRCh38/genes")             #33694
gene_names <- bm$open("/GRCh38/gene_names")   #33694
indices <- bm$open("/GRCh38/indices")         #303795547
data <- bm$open("/GRCh38/data")               #303795547
indptr <- bm$open("/GRCh38/indptr")           #378001
shape <- bm$open("/GRCh38/shape")             #[1]  33694 378000

####### 文件1 barcodes.tsv 的生成 ####### 
barcodes <- barcodes[]
write.csv(barcodes, file = "barcodes.csv") #再手动重命名为.tsv格式

####### 文件2 genes.tsv 的生成 ####### 
genes <- genes[]
gene_names <- gene_names[]
df <- data.frame(gene_names, row.names = c(genes))
write.csv(df, file = "genes.csv") #再手动重命名为.tsv格式

文件3 matrix.mtx 比较特殊，列名是基因的索引号，表示该基因在该细胞中的reads数。第二列没有现成的变量可以一次性填进去，需要进行一些转化和运算：

第一列是indices，作为row.names;
第三列是data，也就是reads数;
第二列是细胞的编号，如1，1，1，2，2，2，2，2，2，3，4，……，33694。这些数字表示即第1个细胞至第33694个细胞。如第一行为33688 1 2 则表示33688号基因在细胞1中的reads数为2。

indptr是一个递增数列，a2-a1 = 细胞cell1的数量。因此也就是说数字1要重复a2-a1次；数字2要重复a3-a2次，以此类推。

indices <- indices[]
indptr <- indptr[]
data <- data[]

####### 文件3 matrix.csv 的生成 ####### 

# 创建两个相差一个元素的数列，相减得到细胞数量列表
indptr1 <- indptr[2:length(indptr)]
indptr2 <- indptr[1:length(indptr)-1]
cell_numb <- indptr1 - indptr2

column2_of_file3 <- rep(c(1:length(indptr1)),c(cell_numb)) # rep函数万岁

df <- data.frame(indices, column2_of_file3, data)
save(df, file = "df_marrow_bone.RData") # 10分钟后运行结束，文件大小672Mb

这里仅仅只是把变量保存，由于磁盘限制没有写入csv中。

这里是利用另一个数据的处理过程，和之前的过程一样，不展开说明：

setwd("/Users/shinianyike/Desktop/zll/data/immune_census/cord_blood/data")
library(hdf5r)
bm <- H5File$new("./ica_cord_blood_h5.h5", mode = "r")
GRCh38 <- bm[["GRCh38"]]

### 读取数据                                   ###数据个数
barcodes <- bm$open("/GRCh38/barcodes")       #384000
genes <- bm$open("/GRCh38/genes")             #33694
gene_names <- bm$open("/GRCh38/gene_names")   #33694
indices <- bm$open("/GRCh38/indices")         #260473471
data <- bm$open("/GRCh38/data")               #260473471
indptr <- bm$open("/GRCh38/indptr")           #384001
shape <- bm$open("/GRCh38/shape")             #[1]  33694 384000

#### 将.h5格式转化为R可操作的list
#### 注意data <- GRCh38[["data"]] 这种读取方式仍然得到的是.h5的格式
barcodes <- barcodes[]
genes <- genes[]
gene_names <- gene_names[]
indices <- indices[]
indptr <- indptr[]
data <- data[]

####### 文件1 barcodes.tsv 的生成 ####### 
write.csv(barcodes, file = "barcodes.csv", row.names = FALSE, col.names = FALSE) #再手动将文件重命名为.tsv结尾

####### 文件2 genes.tsv 的生成 ####### 
df <- data.frame(gene_names, row.names = c(genes))
write.csv(df, file = "genes.csv", row.names = FALSE, col.names = FALSE) #再手动将文件重命名为.tsv结尾

####### 文件3 matrix.tsv 的生成 ####### 

# 创建两个相差一个元素的数列，相减得到细胞数量列表
indptr1 <- indptr[2:length(indptr)]
indptr2 <- indptr[1:length(indptr)-1]
cell_numb <- indptr1 - indptr2

column2_of_file3 <- rep(c(1:length(indptr1)),c(cell_numb)) # rep函数万岁

df <- data.frame(indices, column2_of_file3, data)
save(df, file = "df_cord_blood.RData")

Step4：将数据转换为mtx格式

上述的文件3仅仅只是生成了一个csv文件，然而与输入的mtx格式是不同的。

mtx文件描述：

引自：http://people.sc.fsu.edu/~jburkardt/data/mm/mm.html

A file in the Matrix Market format comprises four parts:

Header line: contains an identifier, and four text fields;

Comment lines: allow a user to store information and comments;

Size line: specifies the number of rows and columns, and the number of nonzero elements;

Data lines: specify the location of the matrix entries (implicitly or explicitly) and their values.

The header line has the form:
%%MatrixMarket object format field symmetry

The header line must be the first line of the file, and the header line must begin with the string %%MatrixMarket. The four fields that follow that string are

object is usually matrix, and that is the case we will consider here. Another legal value is vector, whose format is similar, but with some obvious simplifications.

format is either coordinate or array;

field is either real, double, complex, integer or pattern.

symmetry is either general (legal for real, complex, integer or pattern fields), symmetric (real, complex, integer or pattern), skew-symmetric (real, complex or integer), or hermitian (complex only).

例子：人外周血10Xgenomics数据比对方法所得RNA-seq数据的mtx文件（下图）以%%MatrixMarket字符串开头，object为matrix，format为coordinate，field为real，symmetry为general。分隔符为空格。
第一行为%%MatrixMarket matrix coordinate real general
第二行为%
第三行为基因总数，细胞总数，以及reads总数（第三列的和）
第三行的前两个数字已经有了，只需要计算一下reads总数。

pbmc.mtx

脾脏数据的第一行则为
%%MatrixMarket matrix coordinate integer general

载入R变量

> setwd("/Users/shinianyike/Desktop/zll/data/immune_census/bone_marrow/")
> load("./data/10X/df_marrow_bone.RData")
> head(df)
  indices column2_of_file3 data
1   33664                1    5
2   33662                1    5
3   33661                1   14
4   33659                1   12
5   33658                1   35
6   33657                1    9

# reads总数：
sum(df$data)
[1] 1179704196
# 所以mtx文件的第三行为33694 378000 1179704196

总的来说，要做的事情是在df数据框前添加三行：
第一行：%%MatrixMarket matrix coordinate real general
第二行：%
第三行：33694 378000 303795547

由于第一行和第二行是一行字符，并不像数据框那样行列元素个数能一一对应，因此要考虑其他的文本处理方法。

另外，R在第一行添加元素的操作比较繁杂，尤其是大文件。需要创建新空白文件，写入第一行，再把剩下的一行一行写入，缺点在于慢。若一次性写入又会产生内存不够的问题。

但是，利用linux的sed命令就可以轻松解决这个问题。

正式开始操作：
将R变量传上服务器，载入R变量，写入csv文件：

load("./df_marrow_bone.RData")
head(df)
  indices column2_of_file3 data
1   33664                1    5
2   33662                1    5
3   33661                1   14
4   33659                1   12
5   33658                1   35
6   33657                1    9

write.csv(df, file = "matrix.csv", row.names = FALSE, col.names = FALSE)
# 用时15min左右

Warning message:
In write.csv(df, file = "matrix.csv", row.names = FALSE, col.names = FALSE) :
quit()

### 查看文件
head matrix.csv 
"indices","column2_of_file3","data"
33664,1,5
33662,1,5
33661,1,14
33659,1,12
33658,1,35
33657,1,9
33655,1,15
33654,1,26
33653,1,17

### 处理文件格式
sed -i 's/,/ /g' matrix.csv ###将分隔符逗号改为空格
sed -i '1d' matrix.csv      ###将第一行删除
sed -i '1i %%MatrixMarket matrix coordinate real general\
%\
33694 378000 303795547' matrix.csv ###添加前面三行内容

### 检查文件
head matrix.csv 
%%MatrixMarket matrix coordinate real general
%
33694 378000 303795547
33664 1 5
33662 1 5
33661 1 14
33659 1 12
33658 1 35
33657 1 9
33655 1 15

### 更改后缀名
mv ./matrix.csv ./matrix.mtx

wc -l matrix.mtx  ##查看行数
303795550 matrix.mtx ##行数正确

Step6：genes.tsv和barcodes.tsv的处理

sed -i 's/ENSG/\nENSG/g' genes.tsv ##添加空行
sed -i '$a \' genes.tsv 
sed -i '1d' genes.tsv 
wc -l genes.tsv 
33694 genes.tsv

barcodes同理

参考：

hdf5r包的介绍：https://cran.r-project.org/web/packages/hdf5r/vignettes/hdf5r.html#conversion-of-datatypes
mtx文件描述：http://people.sc.fsu.edu/~jburkardt/data/mm/mm.html

https://github.com/cran/Matrix/blob/master/R/HBMM.R

后续读取数据出现报错

bm_data <- Read10X(data.dir ='./10X/')

Error: readMM(): row     values 'i' are not in 1:nr
#6：12

这里说明第一列的数据中不在1：nrow的范围内，检查后发现有3行0的存在（这是什么坑爹数据啊。。。。），只好手动修改。

2月28日

library(Matrix)
head(df)
     indices column2_of_file3 data
[1,]   33664                1    5
[2,]   33662                1    5
[3,]   33661                1   14
[4,]   33659                1   12
[5,]   33658                1   35
[6,]   33657                1    9

which(df[,1]=="0")
[1]   2446541 121241013 212008491
# 这里就是报错来源。
df[2446540,]
#         indices column2_of_file3             data 
#             128             3014                1 
df[121241012,]
#         indices column2_of_file3             data 
#              12           155705                1 
df[212008490,]
#         indices column2_of_file3             data 
#              72           259211                1 
df[2446541,1]=45
df[121241013,1]=8
df[212008491,1]=20


sparseM <- sparseMatrix(i=df[,1],j=df[,2],x=df[,3])
writeMM(obj = sparseM, file="matrix.mtx")