单细胞小记2:将h5ad格式转化为seurat格式

h5ad 是用于存储单细胞 RNA 测序(scRNA-seq)数据的文件格式,常用于 Python 中的 Scanpy 软件包。它基于 HDF5 文件格式,支持大规模数据的高效存储和处理,特别适合处理复杂的单细胞数据集
ꔷ h5ad 文件格式的关键特点:
1)高效存储:与 HDF5 文件类似,h5ad 文件可以存储大规模稀疏矩阵(如基因表达矩阵)并进行压缩,从而减少存储空间的占用。
2)多维数据支持:能够存储单细胞 RNA-seq 数据的多个层次,包括:
表达矩阵:细胞与基因的稀疏矩阵。
元数据:包括细胞注释、基因注释等信息。
降维数据:如 PCA、UMAP 等嵌入信息。
图结构数据:如邻接图、基于细胞间关系的图数据。
3)跨平台兼容:虽然 h5ad 主要用于 Python 中的 Scanpy,但它可以通过格式转换与其他工具(如 R 中的 Seurat)互通。
4)随机访问:由于 h5ad 文件是基于 HDF5 的,它允许在不加载整个文件的情况下访问文件中的特定部分,从而提高了处理大数据集时的效率

h5Seurat 文件格式专门用于存储和分析多模态单细胞和空间分辨表达实验数据,例如来自 CITE-seq 或 10X Visium 技术的数据。它包含所有分子信息和相关的元数据,包括最近邻图、维度缩减信息、空间坐标和图像数据以及聚类标签。同时还支持 h5Seurat 和 AnnData 对象之间的快速和磁盘上的转换,旨在增强 Seurat 和 Scanpy 之间的互操作性
ꔷ h5Seurat 文件格式的关键特点:
1)高效存储:h5Seurat 使用 HDF5 文件格式,因此能够压缩和存储大规模的矩阵和元数据,节省硬盘空间。
2)跨平台兼容性:可以在不同的机器或系统之间轻松传输,确保兼容性。
3)存储多种数据类型:能够存储稀疏矩阵、嵌入信息(如 PCA、UMAP)、注释数据等多个层次的信息。
4)转换功能:h5Seurat 文件可以转换为 AnnData 格式(用于 Scanpy 分析),这使得 Seurat 和 Scanpy 用户可以相互共享数据。

相似点

1)基于 HDF5:二者都基于 HDF5 格式,HDF5 是一种高效的文件格式,支持大量数据的存储和压缩,尤其适合稀疏矩阵数据的保存。
2)高效存储:h5Seurat 和 h5ad 都支持大规模单细胞数据的存储,包括表达矩阵、元数据、降维结果等。
3)灵活性:它们都允许随机访问,即可以在不加载整个文件的情况下,提取特定部分的数据,如元数据或基因表达矩阵,尤其在处理大数据集时,这非常有用。
4)可转换性:通过工具(如 SeuratDisk 包),h5Seurat 和 h5ad 文件可以互相转换,支持跨平台的分析工作流。

区别

1)面向的分析工具:h5Seurat主要用于 R 中的 Seura,它存储的内容与 Seurat 的分析流程相对应,包括 Seurat 对象的特定结构,比如 assays、reductions(降维结果)、graphs(邻接图结构)、meta.data(元数据)等。h5ad主要用于 Python 中的 Scanpy,存储的内容与 AnnData 对象相关。Scanpy 的 AnnData 对象将表达矩阵、元数据和其他分析结果组合在一起,形成一个单独的对象,具有 Scanpy 特有的数据存储结构。
2)存储结构:h5Seurat是对 Seurat 对象的直接保存,存储了多个 assays(如 RNA、ADT 等),还有 Seurat 特有的嵌入、图形和簇注释等多层数据结构。
h5ad存储了 Scanpy 中的 AnnData 对象,它结构较为统一,主要包含 X(表达矩阵),obs(细胞元数据),var(基因元数据),以及 .obsm、.varm 中存储的降维结果等信息。
3)元数据格式:h5Seurat存储的元数据较为复杂,可以包含多个层次的注释(例如细胞、基因的多个分类注释)。h5ad则对应Scanpy 中的 AnnData 则结构更简单,元数据通常存储在 obs(细胞元数据)和 var(基因元数据)两个数据框中。
4)扩展支持:h5Seurat更加直接地支持 Seurat 中特有的功能,比如基于 RNA velocity 的分析结果、多个不同的 assays、图和网络信息等。h5ad与 AnnData 的设计紧密相关,主要用于与 Scanpy 兼容的分析任务和存储类型,如邻接图、基因表达图谱等。

使用示例:

1 h5ad格式转换为h5seurat格式

library(Seurat)
library(SeuratDisk)

# 将h5ad格式文件转换为h5seurat格式文件,同时指定使用的assay为"RNA"
Convert("GSE174188_CLUES1_adjusted.h5ad" , "h5seurat" , overwrite = TRUE , assay = "RNA")
## Warning: Unknown file type: h5ad
## Creating h5Seurat file for version 3.1.5.9900
## Adding X as scale.data
## Adding raw/X as data
## Adding raw/X as counts
## Adding meta.features from raw/var
## Merging gene_ids from scaled feature-level metadata
## Error in dfgroup$obj_copy_to(dst_loc = dfgroup, dst_name = tname, src_name = i) : 
##   HDF5-API Errors:
##     error #000: ../../../src/H5Ocopy.c in H5Ocopy(): line 233: unable to copy object
##         class: HDF5
##         major: Object header
##         minor: Unable to copy object
## 
##     error #001: ../../../src/H5Ocopy.c in H5O__copy(): line 317: unable to copy object
##         class: HDF5
##         major: Object header
##         minor: Unable to copy object
## 
##     error #002: ../../../src/H5Ocopy.c in H5O__copy_obj(): line 1221: unable to copy object
##         class: HDF5
##         major: Object header
##         minor: Unable to copy object
## 
##     error #003: ../../../src/H5Ocopy.c in H5O__copy_header(): line 1165: unable to copy object
##         class: HDF5
##         major: Object header
##         minor: Unable to copy object
## 
##     error #004: ../../../src/H5Ocopy.c in H5O__copy_header_real(): line 890: unable to re-tag metadata entries
##         class: HDF5
##         major: Object cache
##         minor: Unable to tag metadata in the cache
## 
##     error #005: ../../../src/H5AC.c in H5AC_retag_copied
## Error in private$closeFun(id) : HDF5-API Errors:
##     error #000: ../../../src/H5F.c in H5Fclose(): line 668: closing file ID failed
##         class: HDF5
##         major: File accessibilty
##         minor: Unable to close file
## 
##     error #001: ../../../src/H5Fint.c in H5F__close(): line 2047: decrementing file ID failed
##         class: HDF5
##         major: File accessibilty
##         minor: Unable to decrement reference count
## 
##     error #002: ../../../src/H5I.c in H5I_dec_app_ref(): line 1309: can't decrement ID ref count
##         class: HDF5
##         major: Object atom
##         minor: Unable to decrement reference count
## 
##     error #003: ../../../src/H5Fint.c in H5F__close_cb(): line 2105: can't close file
##         class: HDF5
##         major: File accessibilty
##         minor: Unable to close file
## 
##     error #004: ../../../src/H5Fint.c in H5F_try_close(): line 2276: problems closing file
##         class: HDF5
##         major: File accessibilty
##         minor: Unable to close file
## 
##     error #005: ../../../src/H5Fint.c in H5F__des

\color{green}{错误解析及解决方案:}

1.1 检查 h5ad 文件完整性(python)
>>> import scanpy as sc

# 尝试读取 h5ad 文件
>>> adata = sc.read_h5ad("GSE174188_CLUES1_adjusted.h5ad")
>>> print(adata)
## AnnData object with n_obs × n_vars = 1263676 × 1999
##     obs: 'batch_cov', 'ind_cov', 'Processing_Cohort', 'louvain', 'cg_cov', 'ct_cov', 'L3', 'ind_cov_batch_cov', 'Age', 'Sex', 'pop_cov', 'Status', 'SLE_status'
##     var: 'gene_ids'
##     uns: 'neighbors', 'pca', 'rank_genes_groups', 'umap'
##     obsm: 'X_pca', 'X_umap'
##     varm: 'PCs'
##     obsp: 'connectivities', 'distances'
1.2. 检查 h5ad 文件的元数据
>>> print(adata.obs.head())  # 查看前几行元数据
##                                                              batch_cov  ... SLE_status
## CAAGGCCAGTATCGAA-1-1-0-0-0-0-0-0-0-0-0-0-0-0-0-...  dmx_YS-JY-22_pool6  ...    Healthy
## CTAACTTCAATGAATG-1-1-0-0-0-0-0-0-0-0-0-0-0-0-0-...  dmx_YS-JY-22_pool6  ...        SLE
## AAGTCTGGTCTACCTC-1-1-0-0-0-0-0-0-0-0-0-0-0-0-0-...       dmx_AbFlare-3  ...        SLE
## GGCTCGATCGTTGACA-1-1-0-0-0-0-0-0-0-0-0-0-0-0-0-...  dmx_YS-JY-20_pool3  ...        SLE
## ACACCGGCACACAGAG-1-1-0-0-0-0-0-0-0-0-0-0-0-0-0-...           dmx_YE110  ...        SLE
## 
## [5 rows x 13 columns]

>>> print(adata.var.head())  # 查看前几行基因元数据
##                  gene_ids
## HES4      ENSG00000188290
## ISG15     ENSG00000187608
## TNFRSF18  ENSG00000186891
## TNFRSF4   ENSG00000186827
## MIB2      ENSG00000197530
1.3. 使用不同的 assy 参数(确保指定的 assay 名称是正确的)
>>> print(adata.layers.keys())  # 查看所有的 layers (assays)
KeysView(Layers with keys: )

★ 这里表明在所示范的 h5ad 文件中没有定义任何 layers。可能是因为该文件中只包含了基本的表达矩阵而没有额外的 assays 或 layers,需进一步检查

1.3-1. 检查表达矩阵

>>> print(adata.X)  # 查看主要表达矩阵
## [[-0.21663427 -0.8426997  -0.11158545 ... -0.2582464   2.0164974
  -0.06928813]
##  [-0.1984313  -0.744276   -0.11177281 ... -0.25383508 -0.3862635
  -0.0689166 ]
##  [-0.19595543  0.4312931  -0.14541313 ... -0.3167432  -0.45906577
  -0.09076159]
##  ...
##  [-0.18843444  1.6739295  -0.14349106 ... -0.25323954 -0.38606548
  -0.07936437]
##  [-0.1778934  -0.59904104  6.167296   ... -0.42856053 -0.45448688
  -0.07272661]
##  [-0.10157593 -0.48876476 -0.1120045  ... -0.31424302 -0.34174454
  -0.05100624]]

>>> print(adata.obs.head())  # 查看元数据
##                                                              batch_cov  ... SLE_status
## CAAGGCCAGTATCGAA-1-1-0-0-0-0-0-0-0-0-0-0-0-0-0-...  dmx_YS-JY-22_pool6  ...    Healthy
## CTAACTTCAATGAATG-1-1-0-0-0-0-0-0-0-0-0-0-0-0-0-...  dmx_YS-JY-22_pool6  ...        SLE
## AAGTCTGGTCTACCTC-1-1-0-0-0-0-0-0-0-0-0-0-0-0-0-...       dmx_AbFlare-3  ...        SLE
## GGCTCGATCGTTGACA-1-1-0-0-0-0-0-0-0-0-0-0-0-0-0-...  dmx_YS-JY-20_pool3  ...        SLE
## ACACCGGCACACAGAG-1-1-0-0-0-0-0-0-0-0-0-0-0-0-0-...           dmx_YE110  ...        SLE
## 
## [5 rows x 13 columns]

>>> print(adata.var.head())  # 查看基因元数据
##                  gene_ids
## HES4      ENSG00000188290
## ISG15     ENSG00000187608
## TNFRSF18  ENSG00000186891
## TNFRSF4   ENSG00000186827
## MIB2      ENSG00000197530

1.3-2. 使用默认 assay(在没有层(layers)的情况下,通常会直接使用 adata.X 作为默认的表达矩阵,不指定 assay):

Convert("GSE174188_CLUES1_adjusted.h5ad", dest = "h5seurat")
## Warning: Unknown file type: h5ad
## Warning: 'assay' not set, setting to 'RNA'
## Creating h5Seurat file for version 3.1.5.990
## Adding X as scale.data
## Adding raw/X as data
## Adding raw/X as counts
## Adding meta.features from raw/var
## Merging gene_ids from scaled feature-level metadata
## Adding X_pca as cell embeddings for pca
## Adding X_umap as cell embeddings for umap
## Adding PCs as feature loadings fpr pca
## Adding miscellaneous information for pca
## Adding standard deviations for pca
## Adding miscellaneous information for umap
## Adding rank_genes_groups to miscellaneous data

1.3-3. 简化数据 (创建一个新的 AnnData 对象,只包含主要的表达数据):

>>> adata_simple = sc.AnnData(X=adata.X, obs=adata.obs, var=adata.var)
>>> adata_simple.write("simple_data.h5ad")
Convert("simple_data.h5ad", dest = "h5seurat")
## Warning: Unknown file type: h5ad
## Warning: 'assay' not set, setting to 'RNA'
## Creating h5Seurat file for version 3.1.5.9900
## Adding X as scale.data
## Adding X as data
## Adding X as counts
## Adding meta.features from var

2 使用LoadH5Seurat()函数加载h5seurat格式文件,并创建Seurat对象

seurat_obj <- LoadH5Seurat("GSE174188_CLUES1_adjusted.h5seurat")
## Validating h5Seurat file
## Warning: Feature names cannot have underscores ('_'), replacing with dashes ('-')
## Initializing RNA with data
## Adding counts for RNA(

## Adding scale.data for RNA
## Adding feature-level metadata for RNA
## Adding reduction pca
## Adding cell embeddings for pca
## Adding feature loadings for pca
## Adding miscellaneous information for pca
## Adding reduction umap
## Adding cell embeddings for umap
## Adding miscellaneous information for umap
## Adding command information
## Adding cell-level metadata
## Error: Too many values for levels provided

\color{green}{错误解析及解决方案:}

2.1 只加载表达矩阵部分,先跳过元数据(meta.data)和其他辅助信息(misc),然后再手动添加
scRNA <- LoadH5Seurat("GSE174188_CLUES1_adjusted.h5seurat", meta.data = FALSE, misc = FALSE)
## Validating h5Seurat file
## Warning: Feature names cannot have underscores ('_'), replacing with dashes ('-')
## Initializing RNA with data
## Adding counts for RNA
## Adding scale.data for RNA
## Adding feature-level metadata for RNA
## Adding reduction pca
## Adding cell embeddings for pca
## Adding feature loadings for pca
## Adding miscellaneous information for pca
## Adding reduction umap
## Adding cell embeddings for umap
## Adding miscellaneous information for umap
## Adding command information
## Adding tool-specific results


# 读取元数据
meta_data <- h5read("GSE174188_CLUES1_adjusted.h5seurat", "meta.data")

str(meta_data)
## List of 13
##  $ Age              :List of 2
##   ..$ levels: chr [1:55(1d)] "20.0" "21.0" "22.0" "23.0" ...
##   ..$ values: int [1:1263676(1d)] 9 26 15 51 9 33 45 14 26 11 ...
##  $ L3               :List of 2
##   ..$ levels: chr [1:2(1d)] "0.0" "1.0"
##   ..$ values: int [1:1263676(1d)] 2 2 1 2 1 2 2 1 1 1 ...
##  $ Processing_Cohort:List of 2
##   ..$ levels: chr [1:4(1d)] "1.0" "2.0" "3.0" "4.0"
##   ..$ values: int [1:1263676(1d)] 4 4 3 4 2 2 4 2 2 2 ...
##  $ SLE_status       :List of 2
##   ..$ levels: chr [1:2(1d)] "Healthy" "SLE"
##   ..$ values: int [1:1263676(1d)] 1 2 2 2 2 2 2 1 2 1 ...
##  $ Sex              :List of 2
##   ..$ levels: chr [1:2(1d)] "Female" "Male"
##   ..$ values: int [1:1263676(1d)] 1 1 1 1 1 1 1 1 2 1 ...
##  $ Status           :List of 2
##   ..$ levels: chr [1:4(1d)] "Flare" "Healthy" "Managed" "Treated"
##   ..$ values: int [1:1263676(1d)] 2 3 1 3 3 3 3 2 3 2 ...
##  $ batch_cov        :List of 2
##   ..$ levels: chr [1:23(1d)] "dmx_AbFlare-3" "dmx_AbFlare-4" "dmx_YE110" "dmx_YE_7-13" ...
##   ..$ values: int [1:1263676(1d)] 18 18 1 13 3 8 17 4 8 6 ...
##  $ cg_cov           :List of 2
##   ..$ levels: chr [1:11(1d)] "B" "NK" "PB" "Progen" ...
##   ..$ values: int [1:1263676(1d)] 6 9 9 1 6 7 6 6 6 1 ...
##  $ ct_cov           :List of 2
##   ..$ levels: chr [1:14(1d)] "T4_naive" "T8_naive" "T4_em" "T4_reg" ...
##   ..$ values: int [1:1263676(1d)] 1 0 0 10 0 5 1 1 0 10 ...
##  $ ind_cov          :List of 2
##   ..$ levels: chr [1:261(1d)] "1004_1004" "1014_1014" "1015_1015" "1019_1019" ...
##   ..$ values: int [1:1263676(1d)] 195 24 152 18 76 60 59 221 66 227 ...
##  $ ind_cov_batch_cov:List of 2
##   ..$ levels: chr [1:355(1d)] "1004_1004:dmx_YE_7-13" "1014_1014:dmx_YE_7-13" "1014_1014:dmx_YS-JY-22_pool5" "1015_1015:dmx_YS-JY-20_pool4" ...
##   ..$ values: int [1:1263676(1d)] 245 31 194 22 94 75 74 292 82 310 ...
##  $ louvain          :List of 2
##   ..$ levels: chr [1:24(1d)] "0" "1" "2" "3" ...
##   ..$ values: int [1:1263676(1d)] 2 8 1 4 3 5 16 2 2 4 ...
##  $ pop_cov          :List of 2
##   ..$ levels: chr [1:4(1d)] "African American" "Asian" "European" "Hispanic"
##   ..$ values: int [1:1263676(1d)] 2 3 3 3 2 2 2 3 3 3 ...
2.2 创建新的 meta.data 列表
new_meta_data <- lapply(meta_data, function(x) {
  # 用 levels 替换 values
  return(x$levels[x$values])
})

# 将新列表转换为数据框
new_meta_df <- as.data.frame(new_meta_data, stringsAsFactors = FALSE)

## Error in (function (..., row.names = NULL, check.rows = FALSE, check.names = TRUE,  : 
##   arguments imply differing number of rows: 1263676, 793873
h5ad-1

★ 发现"ct_cov"在转换过程中出现问题,进一步检查,发现“levels”(14个)和“values”(15个)不对应

meta_data[["ct_cov"]][["levels"]]
##  [1] "T4_naive"    "T8_naive"    "T4_em"       "T4_reg"      "CytoT_GZMH+" "CytoT_GZMK+" "T_mait"      "NK_bright"   "NK_dim"     
## [10] "B_naive"     "B_mem"       "B_plasma"    "B_atypical"  "Progen"     

table(meta_data[["ct_cov"]][["values"]])
## 
##      0      1      2      3      4      5      6      7      8      9     10     11     12     13     14 
## 469803 207629  72202  86979  31610  82466  56218  14284  13352  74995  95726  40804   1201   4205  12202 
2.3 根据新的 meta.data 列表,重新赋值(“ct_cov”项选择后续自己重新注释)
scRNA[["Age"]] <- new_meta_data[["Age"]]
scRNA[["L3"]] <- new_meta_data[["L3"]]
scRNA[["Processing_Cohort"]] <- new_meta_data[["Processing_Cohort"]]
scRNA[["SLE_status"]] <- new_meta_data[["SLE_status"]]
scRNA[["Sex"]] <- new_meta_data[["Sex"]]
scRNA[["Status"]] <- new_meta_data[["Status"]]
scRNA[["batch_cov"]] <- new_meta_data[["batch_cov"]]
scRNA[["cg_cov"]] <- new_meta_data[["cg_cov"]]
scRNA[["ind_cov"]] <- new_meta_data[["ind_cov"]]
scRNA[["ind_cov_batch_cov"]] <- new_meta_data[["ind_cov_batch_cov"]]
scRNA[["louvain"]] <- new_meta_data[["louvain"]]
scRNA[["pop_cov"]] <- new_meta_data[["pop_cov"]]
最后编辑于
©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 217,907评论 6 506
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 92,987评论 3 395
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 164,298评论 0 354
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 58,586评论 1 293
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 67,633评论 6 392
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 51,488评论 1 302
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 40,275评论 3 418
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 39,176评论 0 276
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 45,619评论 1 314
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 37,819评论 3 336
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 39,932评论 1 348
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 35,655评论 5 346
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 41,265评论 3 329
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 31,871评论 0 22
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 32,994评论 1 269
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 48,095评论 3 370
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 44,884评论 2 354

推荐阅读更多精彩内容