前言

  今天想跟大家分享一下单细胞分析中的一个需求，就是如何将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('GSM4475049_C52_filtered_feature_bc_matrix.h5','/') #这里的‘/’可以替换为‘matrix’，这样数据的结构就会少一层
> str(hd5)
List of 1
 $ matrix:List of 6
  ..$ barcodes: chr [1:10366(1d)] "AAACCTGAGATATACG-1" "AAACCTGAGATCCGAG-1" "AAACCTGAGCACAGGT-1" "AAACCTGAGCGTGAAC-1" ...
  ..$ data    : int [1:9313106(1d)] 3 1 2 3 8 10 1 8 2 2 ...
  ..$ features:List of 5
  .. ..$ _all_tag_keys: chr [1(1d)] "genome"
  .. ..$ feature_type : chr [1:33538(1d)] "Gene Expression" "Gene Expression" "Gene Expression" "Gene Expression" ...
  .. ..$ genome       : chr [1:33538(1d)] "GRCh38" "GRCh38" "GRCh38" "GRCh38" ...
  .. ..$ id           : chr [1:33538(1d)] "ENSG00000243485" "ENSG00000237613" "ENSG00000186092" "ENSG00000238009" ...
  .. ..$ name         : chr [1:33538(1d)] "MIR1302-2HG" "FAM138A" "OR4F5" "AL627309.1" ...
  ..$ indices : int [1:9313106(1d)] 33508 33507 33506 33505 33503 33502 33499 33498 33497 33496 ...
  ..$ indptr  : int [1:10367(1d)] 0 786 1076 1606 2634 3562 4575 5453 5878 6930 ...
  ..$ shape   : int [1:2(1d)] 33538 10366

解释一下数据里面的变量：
  barcodes：细胞的barcode
  data：就是每个基因在每个细胞中的表达量
  features：记录了5个基因方面的属性
  indices：记录的是features中的基因位置
  indptr：可以计算出细胞中检测到哪些基因
  shape：概况描述，如这里是33538基因，10366细胞
下面来展示一下具体怎么转化，示例代码如下：

>library(Matrix)

>barcodes <- hd5$matrix$barcodes
>barcodes <- gsub('-1','',barcodes) #为了去除barcode中的‘-1’
>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), giveCsparse = FALSE)
#查看一下新建的表达矩阵的结构
> str(mat)
Formal class 'dgTMatrix' [package "Matrix"] with 6 slots
  ..@ i       : int [1:9313106] 33508 33507 33506 33505 33503 33502 33499 33498 33497 33496 ...
  ..@ j       : int [1:9313106] 0 0 0 0 0 0 0 0 0 0 ...
  ..@ Dim     : int [1:2] 33538 10366
  ..@ Dimnames:List of 2
  .. ..$ : chr [1:33538(1d)] "ENSG00000243485" "ENSG00000237613" "ENSG00000186092" "ENSG00000238009" ...
  .. ..$ : chr [1:10366(1d)] "AAACCTGAGATATACG" "AAACCTGAGATCCGAG" "AAACCTGAGCACAGGT" "AAACCTGAGCGTGAAC" ...
  ..@ x       : num [1:9313106] 3 1 2 3 8 10 1 8 2 2 ...
  ..@ factors : list()
> mat[500:510,1:30]
11 x 30 sparse Matrix of class "dgTMatrix"
   [[ suppressing 30 column names ‘AAACCTGAGATATACG’, ‘AAACCTGAGATCCGAG’, ‘AAACCTGAGCACAGGT’ ... ]]
ENSG00000117305 . . . 1 . . 1 . . . . . . . . . . . . . . . . . . . . 1 . .
ENSG00000179163 . . . . 3 1 . . . . . . . 1 . . . . . . . 1 . . . . 1 . . .
ENSG00000188822 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
ENSG00000232557 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
ENSG00000189266 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
ENSG00000188529 . . . . . . . . . . . . . . . . . . . . 1 . . . . . . . . .

  现在是不是觉得很简单，有点意犹未竟的感觉。以后在遇到只有h5表达量文件的时候，你就可以信手拈来不必惊慌了。

最后

  emm，今天就分享到这里，希望可以帮到还不会的朋友们，你是不是也学会了，赶快操作起来吧！