scRNAseq:h5文件转化为matrix表达矩阵

前言

  今天想跟大家分享一下单细胞分析中的一个需求,就是如何将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,今天就分享到这里,希望可以帮到还不会的朋友们,你是不是也学会了,赶快操作起来吧!

最后编辑于
©著作权归作者所有,转载或内容合作请联系作者
平台声明:文章内容(如有图片或视频亦包括在内)由作者上传并发布,文章内容仅代表作者本人观点,简书系信息发布平台,仅提供信息存储服务。

推荐阅读更多精彩内容