Skills enhancement

10XscATAC的结果说明（4）之HDF5

方老师小课堂

前言

上一期我们介绍了10XscATAC数据分析的Matrices结果文件，本期接上一期的结果，介绍存储ATAC数据结果的矩阵文件的另外一类格式HDF5文件，也就是我们在数据分析中经常能见到的h5文件。

何为h5文件

H5（HDF5）全称是 Hierarchical Data Format 5，是一种分层数据格式，专门用于存储和组织大规模、复杂的科学数据。你可以把它理解为一个 “数据容器”—— 就像一个文件夹，里面可以存放不同类型的数据（数组、字符串、表格等），还能给数据加 “标签 / 属性”，方便管理和查找。

核心特点：

● 支持海量数据存储（TB/PB 级），读写效率高；

● 跨平台、跨语言（Python/R 等都支持）；

● 分层结构，数据组织清晰；

● 可以压缩存储，节省空间。

随着生物数据越来越庞大，h5文件的应用频率越来越高，因此，大家对其需要了解如何使用读取存储它。

h5文件数据结构

上面提到了H5文件格式是分层结构，在cellrangerATAC输出的filtered_peak_bc_matrix.h5文件的格式究竟是什么样子？每个层级存储的是什么？我们最关心的矩阵是在哪一层呢？下面我们来具体看下它的数据结构。

该文件第一层是名叫matrix，存放矩阵信息。然后在matrix下的第二层有6个。分别对应的是：

barcodes：常见三联文件的barcodes.tsv.gz结果；

data：常见三联文件的matrix.mtx.gz的结果；

features：常见三联文件的features.tsv.gz结果，该结果下方还有第三层级数据结构，包含6个。

features/_all_tag_keys：除id、name、feature_type以外的对features的额外描述属性。

features/derivation：features的额外描述属性。

features/feature_type：对features的定义，Peaks还是motif还是什么种类。

features/genome：对features所属的基因组种类描述。

features/id：features的各个id，比如peaks就是chr1:1000-2000这样式，比如motif就是SPI1_HUMAN.MA0080.4的名称。

features/name：features的各个名称，基本同id，存在版本差异。这个就和我们常见对基因在不同数据库的名称不一样相似。

indices：data那个数据存放的对应的行索引index位置。

indptr：data那个数据关联indices的每列起始位置的索引index位置。

shape：数据的维度，也就是行列数。

第一种我称为原生态方法，在python中，我们可以使用scipy和tables工具对h5文件进行读取。

import collectionsimport scipy.sparse as sp_sparseimport tablesFeatureBCMatrix = collections.namedtuple('FeatureBCMatrix', ['ids', 'names', 'barcodes', 'matrix'])def get_matrix_from_h5(filename, genome):    with tables.open_file(filename, 'r') as f:        try:            group = f.get_node(f.root, 'matrix')        except tables.NoSuchNodeError:print"Matrix group does not exist in this file."return None        feature_group = getattr(group, 'features').read()        ids = getattr(feature_group, 'id').read()        names = getattr(feature_group, 'name').read()        barcodes = getattr(group, 'barcodes').read()        data = getattr(group, 'data').read()        indices = getattr(group, 'indices').read()        indptr = getattr(group, 'indptr').read()        shape = getattr(group, 'shape').read()        matrix = sp_sparse.csc_matrix((data, indices, indptr), shape=shape)return FeatureBCMatrix(ids, names, barcodes, matrix)filtered_matrix_h5 = "10k_pbmc_ATACv2_nextgem_Chromium_Controller_filtered_peak_bc_matrix.h5"peak_bc_matrix = get_matrix_from_h5(filtered_matrix_h5)matrix = peak_bc_matrix.m

第二种方法，在python中，使用scanpy的read_10x_h5函数进行， 也是我们后续在python中常用的方法。

import pandas as pdimport scanpy as scadata = sc.read_10x_h5("10k_pbmc_ATACv2_nextgem_Chromium_Controller_filtered_peak_bc_matrix.h5")

第三种方法，在R语言中，我们可以使用seurat包的Read10X_h5函数进行读取， 也是我们后续在R中常用的方法。

library(Seurat)counts <- Read10X_h5(filename = "10k_pbmc_ATACv2_nextgem_Chromium_Controller_filtered_peak_bc_matrix.h5")

为什么我们需要对各类结果文件进行解读，主要目的是为了让我们足够了解数据的结构，从而，在面对形形色色的不同格式不同类型数据的读取处理中能够不至于匆匆忙忙，反而更加游刃有余。数据读取的前提建立在你对数据结构的足够了解的基础上，对于后续处理的结构范式的清晰认知，因为后续规定的数据对象的输入也是一个规定好的数据结构，采用不同方法读取后，后续处理，结构转换，进行后续标准流程分析也是有所依赖的。这也是我们在生信分析中常说的数据接口的转换，用通俗的话讲，举个例子，我们可能会遇到scanpy分析结果转R的seurat的机械能后续的，那么这就意味着python的anndata数据对象和R的seurat数据对象的数据转换，这里面就涉及了输出h5文件结构，然后R读取h5结构文件的事情。

目前笔记介绍的相对是比较标准规范的类型。后续我们有机会给大家介绍python的anndata数据对象和R的seurat数据对象的数据相互转换。大家多多关注，方老师会持续分享各类生信小技巧，小知识给大家。谢谢！