跟着Theis Lab学单细胞: 格式转换，你需要知道底层逻辑（上）
跟着Theis Lab学单细胞: 格式转换，你需要知道底层逻辑（上）
单细胞最佳实践
基于磁盘的格式转换
本章导读 SECTION PREVIEW
1. 跨生态系统的必要性、难点和实现方式
2. 了解基于磁盘的数据互操作性方法，如HDF5相关格式
3. 建立跨生态系统分析的实用工作流程策略
本节对应原书 5.1 Motivation， 5.2 Nomenclature，5.3 Disk-based Interoperability
01
格式转换的必要性
对于单细胞分析的初学者，大家很可能都疑惑过：“哪个生态系统最好？我应该从哪个学起？”在前文中，我们也简要的介绍过每个生态系统（Bioconductor， Seurat和Scanpy）的起源和特点。由于许多工具都只支持一种生态系统或者在某种生态系统中表现更优（比如scDblFinder是基于SingleCellExperiment，scVelo基于python），因此我们可能不得不面对在不同格式间切换的问题。
首先我们应该明晰，所有生态系统都可以完成基础的标准分析。然而一旦上升到高级分析，有能力的研究者最好对所有三种生态系统都有了解并且能自如切换，这样才能尽可能地使用最先进、最好的分析工具，而不被格式桎梏。
以目前主流的Seurat和scanpy为例，在这两者间转换不仅涉及格式的区别，更需要研究者熟悉R和Python两种语言及其区别。比如，R对复杂统计学模型有强大的支持，而跟机器学习有关的库几乎都基于Python开发。了解这些区别有助于研究者理解这些方法的教程和文档，从而选择更合适的工具。
好消息是：近年来各种工具已经让不同生态系统之间的数据交换变得相对简单了！其实现方式主要有两种模式：基于磁盘的互操作性和内存中的互操作性。在我们详细介绍这两种模式之前，先区分一下R和Python的包/库和函数：
{package} - R包
package::function() - R包的函数
package - Python库
package.function() - A function in a Python package
02
基于磁盘的互操作性
实现编程语言间数据转换的第一种方式是基于磁盘的互操作性（disk-based interoperability）。简单来说，也就是用一种语言将文件写入到磁盘，再用另一种语言读取的方式。很多情况下这种方式更简单可靠、易于扩展，但对存储空间有较高要求，也牺牲了一定的交互性。这种方式尤其适合完整的分析流程（比如通过nextflow或者snakemake定义的流程），在预先构建好的步骤之间储存/读取文件。然而，对于数据探索（data exploration）和试验各种分析方式（如横向比较、调参等），对交互性要求较高时，基于磁盘的互操作方式就没那么方便了。
简单文本格式
在讨论那些为了单细胞专门设计的格式前，我们可以注意到其实常见的简单文本格式，如CSV、TSV和JSON等，就可以满足一些数据传输的要求。比如，数据较小（几百-几千细胞）、只需要元数据（metadata）的情况下，文本文件就够用了。使用文本文件的好处是几乎与所有编程语言兼容，也不需要单细胞特定的包就能读写。然而，随着数据复杂度增加，这些简单格式很快就会变得不实用。在之前的文章中我们研究过anndata的数据格式，如果使用文本格式，那么每一个slot都需要存为一个单独的文件，而每次读取都要重新构建anndata，显然不太高效。
以CSV格式为例：
# R
# 从Seurat对象提取细胞元数据并保存为CSV
cell_metadata <- [email protected]
write.csv(cell_metadata,"cell_metadata.csv")
# 读取CSV文件
metadata_read <- read.csv("cell_metadata.csv", row.names=1)
# Python
import pandas as pd
# 从AnnData提取细胞元数据并保存为CSV
cell_metadata = adata.obs
cell_metadata.to_csv("cell_metadata.csv")
# 读取CSV文件
metadata_read = pd.read_csv("cell_metadata.csv", index_col=0)
HDF5及衍生格式
HDF5（Hierarchical Data Format version 5）是目前单细胞数据交换的主要格式。它是一种开源文件格式，专为存储大型、复杂和异构数据而设计。
HDF5格式像文件系统一样有目录结构，可以在单个文件中存储多种类型的数据。尽管它非常灵活，但要正确使用，我们需要知道数据存储在哪里以及如何存储。因此，针对单细胞数据研究者开发了一些标准。
1. H5AD
H5AD是scverse生态系统中AnnData对象的HDF5磁盘表示形式，已成为共享单细胞数据集的常用格式。
先创建一个随机的AnnData数据用于演示：
# Python
# 随机生成count数据
counts = csr_matrix(
numpy.random.default_generator().poisson(1, size=(100, 2000)), dtype=numpy.float32
)
adata = anndata.AnnData(counts)
# 生成细胞/基因名
adata.obs_names = [f"Cell_{i:d}"for i inrange(adata.n_obs)]
adata.var_names = [f"Gene_{i:d}"for i inrange(adata.n_vars)]
# 进行一些基础预处理
scanpy.pp.calculate_qc_metrics(adata, inplace=True)
adata.layers["counts"] = adata.X.copy()
scanpy.pp.normalize_total(adata, inplace=True)
scanpy.pp.log1p(adata)
scanpy.pp.highly_variable_genes(adata, inplace=True)
scanpy.tl.pca(adata)
scanpy.pp.neighbors(adata)
scanpy.tl.umap(adata)
adata
AnnData object with n_obs × n_vars = 100 × 2000
obs: 'n_genes_by_counts', 'log1p_n_genes_by_counts', 'total_counts', 'log1p_total_counts', 'pct_counts_in_top_50_genes', 'pct_counts_in_top_100_genes', 'pct_counts_in_top_200_genes', 'pct_counts_in_top_500_genes'
var: 'n_cells_by_counts', 'mean_counts', 'log1p_mean_counts', 'pct_dropout_by_counts', 'total_counts', 'log1p_total_counts', 'highly_variable', 'means', 'dispersions', 'dispersions_norm'
uns: 'log1p', 'hvg', 'pca', 'neighbors', 'umap'
obsm: 'X_pca', 'X_umap'
varm: 'PCs'
layers: 'counts'
obsp: 'distances', 'connectivities'
在Python中读写H5AD非常直接，因为这本就是scverse原生格式：
# 写入H5AD文件
adata.write_h5ad("example.h5ad")
# 读取H5AD文件
adata = scanpy.read_h5ad("example.h5ad")
R语言中有好几个读写H5AD文件的包。然而，由于H5AD是Python的原生格式，部分包只提供一个接口，实际上仍然依靠包装anndata库进行读写。
{Zellkonverter} (R包，H5AD ⇔ SingleCellExperiment)
{zellkonverter}使用{basilisk}包管理Python环境，这意味着用户不需要任何Python知识就能读写H5AD文件。当你首次运行这个包的函数时，它会创建一个特殊的conda环境，然后在内部调用python进行读取。
sce <- zellkonverter::readH5AD(h5ad_file, verbose =TRUE)
ℹ Using the Python reader
ℹ Using anndata version 0.8.0
✔ Read /.../luke.zappia/Downloads/example.h5ad [113ms]
✔ uns$hvg$flavor converted [17ms]
✔ uns$hvg converted [50ms]
✔ uns$log1p converted [25ms]
✔ uns$neighbors converted [18ms]
✔ uns$pca$params$use_highly_variable converted [16ms]
✔ uns$pca$params$zero_center converted [16ms]
✔ uns$pca$params converted [80ms]
✔ uns$pca$variance converted [17ms]
✔ uns$pca$variance_ratio converted [16ms]
✔ uns$pca converted [184ms]
✔ uns$umap$params$a converted [16ms]
✔ uns$umap$params$b converted [16ms]
✔ uns$umap$params converted [80ms]
✔ uns$umap converted [112ms]
✔ uns converted [490ms]
✔ Converting uns to metadata ... done
✔ X matrix converted to assay [29ms]
✔ layers$counts converted [27ms]
✔ Converting layers to assays ... done
✔ var converted to rowData [25ms]
✔ obs converted to colData [24ms]
✔ varm$PCs converted [18ms]
✔ varm converted [47ms]
✔ Converting varm to rowData$varm ... done
✔ obsm$X_pca converted [15ms]
✔ obsm$X_umap converted [16ms]
✔ obsm converted [80ms]
✔ Converting obsm to reducedDims ... done
ℹ varp is empty and was skipped
✔ obsp$connectivities converted [22ms]
✔ obsp$distances converted [23ms]
✔ obsp converted [92ms]
✔ Converting obsp to colPairs ... done
✔ SingleCellExperiment constructed [164ms]
ℹ Skipping conversion of raw
✔ Converting AnnData to SingleCellExperiment ... done
设置verbose=TRUE可以让我们看到{zellkonverter}是如何用python读取数据，并将AnnData的每个部分转换为SingleCellExperiment的。我们再来看看转换完的数据：
sce
class: SingleCellExperiment
dim: 2000 100
metadata(5): hvg log1p neighbors pca umap
assays(2): X counts
rownames(2000): Gene_0 Gene_1 ... Gene_1998 Gene_1999
rowData names(11): n_cells_by_counts mean_counts ... dispersions_norm
varm
colnames(100): Cell_0 Cell_1 ... Cell_98 Cell_99
colData names(8): n_genes_by_counts log1p_n_genes_by_counts ...
pct_counts_in_top_200_genes pct_counts_in_top_500_genes
reducedDimNames(2): X_pca X_umap
mainExpName: NULL
altExpNames(0):
转换完的数据可以正常进行Bioconductor生态的分析。
再来看看写入为H5AD：
zellkonverter_h5ad_file <- tempfile(fileext =".h5ad")
zellkonverter::writeH5AD(sce, zellkonverter_h5ad_file, verbose =TRUE)
ℹ Using anndata version 0.8.0
ℹ Using the 'X' assay as the X matrix
✔ Selected X matrix [29ms]
✔ assays$X converted to X matrix [50ms]
✔ additional assays converted to layers [30ms]
✔ rowData$varm converted to varm [28ms]
✔ reducedDims converted to obsm [68ms]
✔ metadata converted to uns [24ms]
ℹ rowPairs is empty and was skipped
✔ Converting AnnData to SingleCellExperiment ... done
✔ Wrote '/.../.../rj/.../T/.../file102cfa97cc51.h5ad ' [133ms]
写入的文件在Python中用scanpy正常读取：
scanpy.read_h5ad(zellkonverter_h5ad_file)
AnnData object with n_obs × n_vars = 100 × 2000
obs: 'n_genes_by_counts', 'log1p_n_genes_by_counts', 'total_counts', 'log1p_total_counts', 'pct_counts_in_top_50_genes', 'pct_counts_in_top_100_genes', 'pct_counts_in_top_200_genes', 'pct_counts_in_top_500_genes'
var: 'n_cells_by_counts', 'mean_counts', 'log1p_mean_counts', 'pct_dropout_by_counts', 'total_counts', 'log1p_total_counts', 'highly_variable', 'means', 'dispersions', 'dispersions_norm'
uns: 'X_name', 'hvg', 'log1p', 'neighbors', 'pca', 'umap'
obsm: 'X_pca', 'X_umap'
varm: 'PCs'
layers: 'counts'
obsp: 'connectivities', 'distances'
{SeuratDisk} (R包，H5AD ⇔ H5Seurat ⇔ Seurat对象)
使用{SeuratDisk}包在H5AD和Seurat对象中转换需要经过一个中间文件：H5Seurat。这是Seurat的开发者们根据Seurat的格式制作的一个自定义H5格式（与H5AD和AnnData的关系类似）。
⚠️ 注：SeuratDisk暂不支持V5 Assay ("Assay5")，需要转换为V3才能正常运行。
obj[["RNA"]]<- as(obj[["RNA"]],"Assay")
# 读取：先转换为H5Seurat格式，再加载
message("Converting H5AD to H5Seurat...")
SeuratDisk::Convert(h5ad_file, dest ="h5seurat", overwrite =TRUE)
message("Reading H5Seurat...")
h5seurat_file <- gsub(".h5ad",".h5seurat", h5ad_file)
seurat <- SeuratDisk::LoadH5Seurat(h5seurat_file, assays ="RNA")
message("Read Seurat object:")
seurat
Converting H5AD to H5Seurat...
The legacy packages maptools, rgdal, and rgeos, underpinning the sp package,
which was just loaded, will retire in October 2023.
Please refer to R-spatial evolution reports for details, especially
https://r-spatial.org/r/2023/05/15/evolution4.html.
It may be desirable to make the sf package available;
package maintainers should consider adding sf to Suggests:.
The sp package is now running under evolution status 2
(status 2 uses the sf package in place of rgdal)
Registered S3 method overwritten by 'SeuratDisk':
method from
as.sparse.H5Group Seurat
Warning: Unknown file type: h5ad
Warning: 'assay' not set, setting to 'RNA'
Creating h5Seurat file for version 3.1.5.9900
Adding X as data
Adding X as counts
Adding meta.features from var
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 hvg to miscellaneous data
Adding log1p to miscellaneous data
Adding layer counts as data in assay counts
Adding layer counts as counts in assay counts
Reading H5Seurat...
Validating h5Seurat file
Warning: Feature names cannot have underscores ('_'), replacing with dashes ('-')
Initializing RNA with data
Adding counts 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
Read Seurat object:
An object of class Seurat
2000 features across 100 samples within 1 assay
Active assay: RNA (2000 features, 0 variable features)
2 dimensional reductions calculated: pca, umap
由于Seurat对象跟SingleCellExperiment和AnnData的结构有很大的区别，因此转换会更加复杂。
{sceasy} (R包，SingleCellExperiment，AnnData，Loom ⇔ Seurat)
sceasy也通过Python实现一步到位的转换，不过与zellkonverter不同，sceasy需要用户自行设置conda环境。
⚠️ 注：sceasy安装前需要事先设置好conda环境。具体操作如下：
conda install -c bioconda anndata loompy还要安装以下R包：
if(!requireNamespace("BiocManager", quietly =TRUE)) > install.packages("BiocManager") > > BiocManager::install(c("LoomExperiment","SingleCellExperiment")) > install.packages('reticulate')
这里展示AnnData和Seurat对象之间的转换。其他格式，只需要更改“from” 和“to”的格式就可以了。
# R
# 安装：
devtools::install_github("cellgeni/sceasy")
library(sceasy)
library(reticulate)
use_condaenv('EnvironmentName')
# anndata到seurat
sceasy_seurat <- sceasy::convertFormat(h5ad_file, from="anndata", to="seurat")
sceasy_seurat
# seurat到anndata
sceasy::convertFormat(seurat_object, from="seurat", to="anndata",
outFile='filename.h5ad')
Warning: Feature names cannot have underscores ('_'), replacing with dashes ('-')
X -> counts
An object of class Seurat
2000 features across 100 samples within 1 assay
Active assay: RNA (2000 features, 0 variable features)
2 dimensional reductions calculated: pca, umap
{anndata} (R包，仅用于读取H5AD)
R语言也有一个{anndata}包，用于提供与python对象AnnData的交互界面，以读取对象中的数据。但是，这个包并不将AnnData转换为原生的R对象，而且大部分分析工具也不接受这种格式作为输入，所以多数情况下仍然需要进一步转换。
Loom
Loom文件格式是另一种用于组学数据的HDF5规范。与H5AD不同，它不与特定分析生态系统绑定，但结构相似。Loom格式在velocyto（RNA速率分析）生成的数据中很常见。R和Python都有实现Loom格式的包（如python的loompy库，和上文提到的{sceasy}），但通常使用核心生态系统包提供的高级接口更方便。由于在一般组学数据分析中使用较少，此处我们不做赘述。
RDS
RDS格式是R特有的二进制格式，可用于存储/加载任意R对象。RDS文件常用于共享R分析结果，尤其是SingleCellExperiment和Seurat对象。
虽然RDS在R生态系统内部使用方便，但不建议将其用于公开共享数据或与合作者交流，因为它缺乏与其他生态系统的互操作性。相反，我们推荐使用上述提到的HDF5格式。
⚠️注：在前几天的推文中，我们也提到rds存在内存和磁盘占用大、读写耗时长的问题。在R中读写单细胞数据，可以考虑使用{qs}/{qs2}。
新兴的磁盘格式
虽然目前的行业标准是HDF5和衍生格式，有一些新兴的格式，比如Zarr、TileDB，在储存超大数据或多模态数据时也有独特优势。可以预见，未来这些格式也会开发出规范并被社区采用（anndata已经提供对Zarr文件的支持），但我们在此不展开讨论。有兴趣的读者可以自行了解。
03
后记
本期我们介绍了R和Python之间的基于磁盘的数据转换方法，并介绍了几个R中的转换包。然而，有相关经验的小伙伴应该知道，这些数据转换包常常缺乏长期维护、跟不上主要生态系统的更新（点名h5Seurat）。因此，知道多种转换方法很有必要。下一期我们将讨论另一种策略：基于内存的转换，敬请期待！本章英文原文及完整代码，可以在后台回复scbp05获取。
----------
单细胞领域大牛编写的最强硬核单细胞教程来啦！想要追更学习笔记，请给我们加一个小✨，这样最新的推送才会显示在主页：