前言
今天想跟大家分享一下单细胞分析中的一个需求，就是如何将h5的表达量文件转化为常见的表达矩阵。大部分情况下，我们是不需要用到这个需求，因为通常我们见到的表达量文件是下面这三个文件：barcodes.tsv.gz、features.tsv.gz、matrix.mtx.gz。大部分分析单细胞的软件都可以读取这个文件做分析，如果有些软件不行的话，你也可以用单细胞的三大R包Seurat、Monocle、Scater来读取转化需要的格式，曲线救国也是可以的。
假如现在你只有h5文件怎么办？你可能会说既然都能找到h5文件肯定也可以找到常用的表达量文件，当然了我也不会反对这种可能的。但是你有没有想过用h5文件来直接转化为表达矩阵呢！那么下面我就跟大家分享一下具体如何操作。
首先先来介绍一下h5文件是什么，我估计有不少小伙伴都不怎么了解h5文件，下面来科普一下h5的基本格式介绍：H5文件是层次数据格式第5代的版本(Hierarchical Data Format，HDF5)，h5格式主要有两个主要概念：dataset和group。 Dataset是类似于数组的数据集，而group是类似文件夹一样的容器，存放dataset和其他group。为了有直观地印象，可以把H5文件想想成一个文件系统，一个h5的最顶层是‘/’，然后是不同的组‘group’，每个组下面是数据‘dataset’，格式类似如下：
+-- /
| +-- group1
| | +-- dataset1
| | | +-- attribute1
| | | +-- attribute2
| | | +-- ...
| | |
| | +-- dataset2
| | | +-- attribute1
| | | +-- attribute2
| | | +-- ...
现在脑海中有了直观的印象，那么现在我们打开表达量的h5文件了，R包rhdf5、hdf5r都可以用来打开文件，这里我使用rhdf5包，示例代码如下：
library(rhdf5)
> hd5 <- h5read('sample_bc_matrix.h5','/') #这里的‘/’可以替换为‘matrix’，这样数据的结构就会少一层
> str(hd5)
List of 1
$ matrix:List of 6
..$ barcodes: chr [1:23130(1d)] "ccRCC_1_rep1@AAATGCCGTCTCAACA-1" "ccRCC_1_rep1@AACACGTAGATCACGG-1" "ccRCC_1_rep1@AACTTTCAGTCATCCA-1" "ccRCC_1_rep1@ACACCCTGTCGACTAT-1" ...
..$ data : num [1:58298590(1d)] 1.54 1.54 1.54 1.54 2.13 ...
..$ features:List of 5
.. ..$ _all_tag_keys: chr [1(1d)] "genome"
.. ..$ feature_type : chr [1:18104(1d)] "Gene" "Gene" "Gene" "Gene" ...
.. ..$ genome : chr [1:18104(1d)] "GRCh38" "GRCh38" "GRCh38" "GRCh38" ...
.. ..$ id : chr [1:18104(1d)] "FO538757.1" "AP006222.1" "AL732372.2" "AL669831.5" ...
.. ..$ name : chr [1:18104(1d)] "FO538757.1" "AP006222.1" "AL732372.2" "AL669831.5" ...
..$ indices : int [1:58298590(1d)] 19 31 37 58 81 98 109 237 248 264 ...
..$ indptr : int [1:23131(1d)] 0 968 2263 3364 3919 4920 6049 7129 8248 9509 ...
..$ shape : int [1:2(1d)] 18104 23130
解释一下数据里面的变量：
barcodes：细胞的barcode
data：就是每个基因在每个细胞中的表达量
features：记录了5个基因方面的属性
indices：记录的是features中的基因位置
indptr：可以计算出细胞中检测到哪些基因
shape：概况描述，如这里是18104 基因，23130细胞
下面来展示一下具体怎么转化，示例代码如下：
>library(Matrix)
>barcodes <- hd5$matrix$barcodes
>counts <- hd5$matrix$data
>gene_id <- hd5$matrix$features$id
>indices <- hd5$matrix$indices
>indptr <- hd5$matrix$indptr
>shape <- hd5$matrix$shape
#用Matrix包的函数轻松转化为表达矩阵
>mat <- sparseMatrix(i = indices[], p = indptr[], x = as.numeric(x = counts[]), dims = shape[], dimnames = list(gene_id,barcodes))
#查看一下新建的表达矩阵的结构
>str(mat)
Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
..@ i : int [1:58298590] 19 31 37 58 81 98 109 237 248 264 ...
..@ p : int [1:23131] 0 968 2263 3364 3919 4920 6049 7129 8248 9509 ...
..@ Dim : int [1:2] 18104 23130
..@ Dimnames:List of 2
.. ..$ : chr [1:18104(1d)] "FO538757.1" "AP006222.1" "AL732372.2" "AL669831.5" ...
.. ..$ : chr [1:23130(1d)] "ccRCC_1_rep1@AAATGCCGTCTCAACA-1" "ccRCC_1_rep1@AACACGTAGATCACGG-1" "ccRCC_1_rep1@AACTTTCAGTCATCCA-1" "ccRCC_1_rep1@ACACCCTGTCGACTAT-1" ...
..@ x : num [1:58298590] 1.54 1.54 1.54 1.54 2.13 ...
..@ factors : list()
>mat[19:26,58:59]
8 x 2 sparse Matrix of class "dgCMatrix"
ccRCC_1_rep2@CAGATCAAGCAATATG-1 ccRCC_1_rep2@CAGCGACCACAGACTT-1
TNFRSF4 1.2251890 .
SDF4 0.9567439 .
B3GALT6 . .
C1QTNF12 . .
AL162741.1 . .
UBE2J2 0.5886769 .
SCNN1D . .
ACAP3 . 1.58509
现在是不是觉得很简单，有点意犹未竟的感觉。以后在遇到只有h5表达量文件的时候，你就可以信手拈来不必惊慌了。
最后
emm，今天就分享到这里，希望可以帮到还不会的朋友们，你是不是也学会了，赶快操作起来吧！
