【생물정보학】 HDF 파일과 filtered_feature_bc_matrix.h5

HDF 파일과 filtered_feature_bc_matrix.h5

 

추천글 : 【생물정보학】 생물정보학 분석 목차

---

1. HDF 파일 [본문]

2. filtered_feature_bc_matrix.h5 [본문]

---

 

1. HDF 파일(hierarchical data format files) [목차]

⑴ 개요

① 대용량 파일을 용량 제한 없이 쉽고 빠르게 입출력하기 위해 사용됨

② environment와 무관하게 파일 교환되는 방식 중 하나 

③ aerospace, physics, engineering, finance, academic research, genomics, astronomy, electronics, medicine 등에서 사용

④ HDF4 파일 (.h4)과 HDF5 파일 (.h5)이 많이 쓰임

⑵ Python을 활용한 HDF 파일 생성, 읽기, 그룹화 등

① numpy를 이용한 쓰기  

 

import numpy as np
import h5py

A = np.random.randint(100, size=(4,4))
print(A)
# [[93 70 70 57]
#  [99 74 27 40]
#  [37  4 71 30]
#  [58 89 86 73]]

B = np.random.randint(100, size=(5,3,3))
print(B)
# [[[49  3 33]
#  [ 4  8 32]
#  [18  8 36]]
# 
# [[90 84 26]
#  [ 6 91 98]
#  [63 55 16]]
# 
# [[37 75 57]
#  [17 68 95]
#  [39 63 96]]
#
# [[19 77 25]
#  [ 4 27 28]
#  [47 34 92]]
#
# [[49  5 33]
#  [85 89 76]
#  [17  5 41]]]

### create a hdf file ###
f1 = h5py.File("myData.hdf5", "w")

### save data in the hdf file ###
dset1 = f1.create_dataset("dataset_01", (4,4), dtype='i', data=A)
dset2 = f1.create_dataset("dataset_02", (5,3,3), dtype = 'i', data=B)

### add metadata ###
dset1.attrs['scale'] = 0.01 # add metadata
dset1.attrs['offset'] = 15 # add metadata

### close the file ###
f1.close

 

② numpy를 이용한 읽기 

 

f2 = h5py.File("myData.hdf5", 'r')
print(list(f2.keys()))
# ['dataset_01', 'dataset_02']

dset1 = f2['dataset_01']
data = dset1[:]
print(data)
# [[93 70 70 57]
#  [99 74 27 40]
#  [37  4 71 30]
#  [58 89 86 73]]

print(list(dset1.attrs.keys()))
# ['offset', 'scale']

print(dset1.attrs['scale'])
# 0.01

 

③ pandas를 이용한 쓰기

 

import pandas as pd
import numpy as np

data = np.arange(1,13)
data = data.reshape(3,4)
print(data)
# [[ 1  2  3  4]
#  [ 5  6  7  8]
#  [ 9 10 11 12]]

columns = ['a','b','c','d']
index = ['A','B','C']

df = pd.DataFrame(data = data, index = index, columns = columns)
print(df)
#    a   b   c   d
# A  1   2   3   4
# B  5   6   7   8
# C  9  10  11  12

### store the data using HDFStore ###
store = pd.HDFStore('myData.hdf5')
store.put('dataset_01', df)
metadata = {'scale':0.1,'offset':15}
store.get_storer('dataset_01').attrs.metadata = metadata
store.close()

 

○ 데이터 쓰기는 2D 데이터만 가능한 듯  

○ 3D 데이터를 쓰려고 하면 다음과 같은 오류가 발생 

 

ValueError: Must pass 2-d input. shape=(512, 512, 3)

 

④ pandas를 이용한 읽기 

 

import pandas as pd

with pd.HDFStore('myData.hdf5') as store:
    data = store['dataset_01']
    metadata = store.get_storer('dataset_01').attrs.metadata
    
print(data)
#    a   b   c   d
# A  1   2   3   4
# B  5   6   7   8
# C  9  10  11  12

print(metadata)
# {'scale': 0.1, 'offset': 15}

 

⑤ 레퍼런스 (ref1, ref2)

 

h5py: reading and writing HDF5 files in Python

How to save a large dataset in a hdf5 file using python ? (Quick Guide)

 

⑶ HDF 파일 뷰어

① Windows용 소프트웨어 : GDAL, Golden Software Surfer, Safe Software FME Desktop

② Mac OS용 소프트웨어 : NCSA HDFView, Basic ENVISAT Atmospheric Toolbox, WaveMetrics IGOR Pro, Netron

③ Linux용 소프트웨어 : GDAL, NCSA HDFView, Basic ENVISAT Atmospheric Toolbox  

④ 기타 : HDFView, h5dump command (Linux)

 

 

2. filtered_feature_bc_matrix.h5 [목차]

⑴ HDF5 file hierarchy 

 

(root)
└── matrix [HDF5 group]
    ├── barcodes
    ├── data
    ├── indices
    ├── indptr
    ├── shape
    └── features [HDF5 group]
        ├─ _all_tag_keys
        ├─ target_sets [for Targeted Gene Expression]
        │   └─ [target set name]
        ├─ feature_type
        ├─ genome
        ├─ id
        ├─ name
        ├─ pattern [Feature Barcode only]
        ├─ read [Feature Barcode only]
        └─ sequence [Feature Barcode only]

 

⑵ 문제점 1. (ST) 대부분의 공간전사체 파이프라인은 .h5 파일이 있다는 가정 하에 진행됨

① 개요  

○ scanpy.read_visium (python), Load10X_Spatial (R) 등의 경우 .h5 파일만을 입력으로 함

○ 어떤 데이터베이스의 경우 .h5 파일을 공개하지 않고, 그에 파생적인 다음 파일들만을 제공

○ barcodes.tsv.gz 

○ features.tsv.gz

○ matrix.mtx.gz 

○ 위 세 가지 파일만을 가지고 .h5 파일을 재구성(creation or restoration)해야 할 필요가 있음

② 해결방법 1. R sceasy::convertFormat : 실패 

 

library(Seurat)
library(reticulate)
library(sceasy)

pbmc.data <- Seurat::Read10X(data.dir = 
    "./Data/CID44971_spatial/CID44971_raw_feature_bc_matrix/")
pbmc <- Seurat::CreateSeuratObject(counts = pbmc.data)
sceasy::convertFormat(pbmc, from="seurat", to="anndata", 
    outFile='./Data/CID44971_spatial/CID44971_raw_feature_bc_matrix/filtered_feature_bc_matrix.h5')

 

○ 원인 : .h5 파일을 생성할 때 v2 또는 v3에 따라 genome 정보가 추가되는지가 달라져서 위 문제가 발생함

③ 해결방법 2. tissue_dir 디렉토리 내 matrix.mtx, barcodes.tsv, features.tsv 및 spatial 폴더가 있을 때, Visium 데이터를 R에서 읽는 코드 : 성공 

 

library(Seurat)

b_data <-ReadMtx('./matrix.mtx.gz',
               './barcodes.tsv.gz',
               './features.tsv.gz',
               feature.column=1)
b_data = CreateSeuratObject(b_data, assay='Spatial')
b_image = Read10X_Image(paste('./spatial/',sep=''))
b_image_ = b_image[Cells(b_data)]
DefaultAssay(object = b_image_) <- 'Spatial'
b_data[["slice1"]]=b_image_
# you can change "slice1"

b_data <- SCTransform(b_data, assay = "Spatial", verbose = FALSE, variable.features.n = 1000)
# Some errors might occur when running SpatialFeaturePlot without SCTransform

SpatialFeaturePlot(b_data, rownames(b_data)[1])
# check whether it works properly

 

④ 해결방법 3. tissue_dir 디렉토리 내 matrix.mtx, barcodes.tsv, features.tsv 및 spatial 폴더가 있을 때, Visium 데이터를 Python에서 읽는 코드 : 성공 

 

import scanpy as sc
import pandas as pd
import matplotlib.pyplot as plt
from matplotlib.image import imread
import numpy as np
import cv2
import math

from pathlib import Path, PurePath
from typing import Union, Dict, Optional, Tuple, BinaryIO
import json

def to_spatial(adata, path, library_id = 'TNBC1'):
    
    adata.uns["spatial"][library_id] = dict()
    path = Path(path)
    
    files = dict(
        tissue_positions_file=path / 'spatial/tissue_positions_list.csv',
        scalefactors_json_file=path /  'spatial/scalefactors_json.json',
        hires_image=path / 'spatial/tissue_hires_image.png',
        lowres_image=path / 'spatial/tissue_lowres_image.png',
    )

    # check if files exists, continue if images are missing
    for f in files.values():
        if not f.exists():
            if any(x in str(f) for x in ["hires_image", "lowres_image"]):
                logg.warning(
                    f"You seem to be missing an image file.\n"
                    f"Could not find '{f}'."
                )
            else:
                raise OSError(f"Could not find '{f}'")

    adata.uns["spatial"][library_id]['images'] = dict()
    for res in ['hires', 'lowres']:
        try:
            adata.uns["spatial"][library_id]['images'][res] = imread(
                str(files[f'{res}_image'])
            )
        except Exception:
            raise OSError(f"Could not find '{res}_image'")

    # read json scalefactors
    adata.uns["spatial"][library_id]['scalefactors'] = json.loads(
        files['scalefactors_json_file'].read_bytes()
    )

    #adata.uns["spatial"][library_id]["metadata"] = {
    #    k: (str(attrs[k], "utf-8") if isinstance(attrs[k], bytes) else attrs[k])
    #    for k in ("chemistry_description", "software_version")
    #    if k in attrs
    #}

    # read coordinates
    positions = pd.read_csv(files['tissue_positions_file'], header=None)
    positions.columns = [
        'barcode',
        'in_tissue',
        'array_row',
        'array_col',
        'pxl_col_in_fullres',
        'pxl_row_in_fullres',
    ]
    positions.index = positions['barcode']

    adata.obs = adata.obs.join(positions, how="left")
    
    adata.obsm['spatial'] = adata.obs[
        ['pxl_row_in_fullres', 'pxl_col_in_fullres']
    ].to_numpy()
    
    adata.obs.drop(
        columns=['barcode', 'pxl_row_in_fullres', 'pxl_col_in_fullres'],
        inplace=True,
    )

    return adata


def load_adata(tissue_dir):
    adata1 = sc.read_mtx(tissue_dir + 'matrix.mtx') 
    adata1 = adata1.transpose() 
    cellname = pd.read_csv(tissue_dir + 'barcodes.tsv', sep="\t", header=None, index_col=1, names=['idx', 'barcode']) 
    cellname.index = cellname['idx']
    featurename = pd.read_csv(tissue_dir + 'features.tsv', sep='\t', header=None, index_col=1, names=['gene_ids', 'feature_types'])
    adata1.var = featurename
    adata1.obs = cellname
    adata1.uns["spatial"] = dict()
    adata1 = to_spatial(adata1, tissue_dir, library_id='Library_1')
    adata1 = adata1[adata1.obs.in_tissue == 1,:]
    adata1.var_names_make_unique()
    adata1.obs['unnormalized_counts'] = adata1.X.sum(axis=1).A1
    sc.pp.normalize_total(adata1, inplace=True)
    adata1.obs['normalized_counts'] = adata1.X.sum(axis=1).A1
    sc.pp.log1p(adata1)
    sc.pp.highly_variable_genes(adata1, flavor="seurat", n_top_genes=2000)

    return adata1

 

⑤ 트러블슈팅 1. Warning message: "Invalid name supplied, making object name syntactically valid. New object name is X1160920F; see ?make.names for more details on syntax validity" 

○  디렉토리명이 숫자가 아니라 문자로 시작되도록 하면 해결됨

○ 위 코드의 경우 디렉토리명을 구체적으로 언급하지 않았으므로 해당 없음 

⑥ 트러블슈팅 2. Exception: File is missing one or more required datasets. (파이썬)

○ 정해진 문법에 따라 생성된 .h5 파일이 아니면 위 오류 메시지가 출력됨

⑶ 문제점 2. (sc 또는 ST) matrix.mtx.gz, barcodes.tsv, 또는 features.tsv가 깨진 파일인 경우

① matrix.mtx.gz 파일이 깨진 경우 R 환경에서 readMM으로 읽어들인 뒤 writeMM으로 깨진 파일을 복구

② barcodes.tsv 및 features.tsv의 경우 R 환경에서 read.table로 읽어들인 뒤 write.table로 깨진 파일을 복구 

⑷ 문제점 3. (sc 또는 ST) R 혹은 파이썬 객체로부터 raw data를 추출

① 해결방법 1. R Seurat object가 있을 때 matrix.mtx, barcodes.tsv, features.tsv 파일을 추출하는 코드 : 성공

 

library(Seurat)

counts <- seurat_obj@assays$RNA@counts 
barcodes <- rownames(seurat_obj@meta.data)
features <- rownames(counts)

library(Matrix)
writeMM(counts, file = "~/Downloads/matrix.mtx")

write.table(barcodes, file = "~/Downloads/barcodes.tsv", row.names = FALSE, col.names = FALSE, quote = FALSE, sep = "\t")

gene_ids <- features
gene_names <- features
features_df <- data.frame(gene_ids, gene_names)
write.table(features_df, file = "~/Downloads/features.tsv", row.names = FALSE, col.names = FALSE, quote = FALSE, sep = "\t")

 

⑸ 응용 1. h5 파일로부터 각 바코드별 리드수를 알려주는 파이썬 코드

 

def get_mapped_reads(h5_file):
    with h5py.File(h5_file, 'r') as f:
        # Access the matrix group in the HDF5 file
        matrix_group = f['matrix']

        # Read the barcode, gene, and data datasets
        barcodes = np.array(matrix_group['barcodes']).astype(str)
        genes = np.array(matrix_group['features/name']).astype(str)
        data = np.array(matrix_group['data'])

        # Sum the data for each barcode to get the number of mapped reads per barcode
        barcode_read_counts = np.array(matrix_group['indptr'][1:] - matrix_group['indptr'][:-1])

        return barcodes, genes, barcode_read_counts

 

입력: 2022.01.04 09:04

수정: 2023.10.17 18:22