Learning Single-Cell Analysis with Theis Lab: Format Conversion — Understanding the Underlying Logic (Part 1)
Single-Cell Best Practices
Disk-Based Format Conversion
SECTION PREVIEW
1. Why cross-ecosystem interoperability is necessary, what challenges it poses, and how to achieve it
2. Understanding disk-based data interoperability methods, including HDF5 and related formats
3. Building practical workflow strategies for cross-ecosystem analysis
This section corresponds to chapters 5.1 Motivation, 5.2 Nomenclature, and 5.3 Disk-based Interoperability of the original book.
01
Why Format Conversion?
If you're new to single-cell analysis, you've probably asked yourself: "Which ecosystem is best? Where should I even start?" In previous articles, we briefly introduced the origins and characteristics of the three main ecosystems — Bioconductor, Seurat, and Scanpy. Because many tools only support a single ecosystem, or perform best within a particular one (for example, scDblFinder is built on SingleCellExperiment, while scVelo is Python-based), researchers often find themselves needing to switch between formats.
First, let's be clear: all three ecosystems are capable of performing standard single-cell analyses. However, for advanced analyses, capable researchers are best served by having familiarity with all three ecosystems and the flexibility to switch between them. Only then can you take full advantage of cutting-edge tools without being constrained by format incompatibilities.
Take the two most popular ecosystems today — Seurat and Scanpy — as an example. Switching between them isn't just a matter of file format; it also requires familiarity with both R and Python and an understanding of their respective strengths. R excels at complex statistical modeling, while virtually all machine learning libraries are developed in Python. Knowing these differences helps researchers navigate tutorials and documentation more effectively and choose the right tool for the job.
The good news: in recent years, various tools have made cross-ecosystem data exchange quite straightforward! There are two main approaches: disk-based interoperability and in-memory interoperability. Before we dive into either, let's first establish a notation convention for R and Python packages/libraries and functions:
{package}— R packagepackage::function()— a function from an R packagepackage— Python librarypackage.function()— a function from a Python library
02
Disk-Based Interoperability
The first approach to data conversion between programming languages is disk-based interoperability. In simple terms, this means writing a file to disk in one language and reading it in another. In many situations, this approach is simpler, more reliable, and easier to scale — though it requires more storage space and sacrifices some interactivity. It is especially well-suited for end-to-end analysis pipelines (such as those defined with Nextflow or Snakemake), where data is written and read between predefined steps. However, for exploratory data analysis or iterative experimentation (e.g., side-by-side comparisons, parameter tuning), which demand higher interactivity, disk-based approaches are less convenient.
Plain Text Formats
Before discussing formats specifically designed for single-cell data, it's worth noting that common plain text formats — CSV, TSV, and JSON — can satisfy some basic data transfer needs. For small datasets (hundreds to thousands of cells) or when only metadata is needed, plain text files can do the job. Their biggest advantage is near-universal compatibility with any programming language, with no need for single-cell-specific packages to read or write them. However, as data complexity grows, these formats quickly become impractical. As we explored in a previous article, an AnnData object has many slots — storing it as plain text would require a separate file for each slot, and every read would require reconstructing the entire AnnData object. Clearly not ideal.
Here's an example using CSV format:
# R
# Extract cell metadata from a Seurat object and save as CSV
cell_metadata <- [email protected]
write.csv(cell_metadata, "cell_metadata.csv")
# Read back the CSV file
metadata_read <- read.csv("cell_metadata.csv", row.names = 1)
# Python
import pandas as pd
# Extract cell metadata from an AnnData object and save as CSV
cell_metadata = adata.obs
cell_metadata.to_csv("cell_metadata.csv")
# Read back the CSV file
metadata_read = pd.read_csv("cell_metadata.csv", index_col=0)
HDF5 and Derivative Formats
HDF5 (Hierarchical Data Format version 5) is the dominant format for single-cell data exchange today. It is an open-source file format designed to store large, complex, and heterogeneous datasets.
HDF5 has a filesystem-like directory structure, allowing multiple types of data to be stored within a single file. Despite its great flexibility, using it correctly requires knowing where and how data is stored. To address this, standardized specifications have been developed for single-cell researchers.
1. H5AD
H5AD is the HDF5 on-disk representation of AnnData objects in the scverse ecosystem and has become a widely adopted format for sharing single-cell datasets.
Let's first create a random AnnData object for demonstration:
# Python
# Generate random count data
counts = csr_matrix(
numpy.random.default_generator().poisson(1, size=(100, 2000)), dtype=numpy.float32
)
adata = anndata.AnnData(counts)
# Generate cell/gene names
adata.obs_names = [f"Cell_{i:d}" for i in range(adata.n_obs)]
adata.var_names = [f"Gene_{i:d}" for i in range(adata.n_vars)]
# Basic preprocessing
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'
Reading and writing H5AD in Python is straightforward, since it is the native format of the scverse ecosystem:
# Write to H5AD
adata.write_h5ad("example.h5ad")
# Read from H5AD
adata = scanpy.read_h5ad("example.h5ad")
Several R packages support reading and writing H5AD files. However, since H5AD is a Python-native format, some packages merely provide a wrapper interface and internally rely on the Python anndata library for the actual I/O.
{zellkonverter} (R package, H5AD ⇔ SingleCellExperiment)
{zellkonverter} uses the {basilisk} package to manage a Python environment under the hood, which means users need no Python knowledge to read and write H5AD files. The first time you call any of its functions, it creates a dedicated conda environment and internally invokes Python for the file operations.
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
Setting verbose = TRUE lets us trace exactly how {zellkonverter} reads the data through Python and maps each slot of the AnnData to its SingleCellExperiment counterpart. Let's inspect the converted object:
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):
The converted object is fully compatible with downstream Bioconductor analysis.
Now let's look at writing back to 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]
The written file can be read back normally with scanpy in Python:
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 package, H5AD ⇔ H5Seurat ⇔ Seurat object)
Converting between H5AD and Seurat objects using {SeuratDisk} requires an intermediate file format: H5Seurat. This is a custom HDF5 specification developed by the Seurat team that mirrors Seurat's internal structure — analogous to the relationship between H5AD and AnnData.
⚠️ Note: SeuratDisk does not currently support the V5 Assay format ("Assay5"). You must first downgrade to V3 for the conversion to work properly:
obj[["RNA"]] <- as(obj[["RNA"]], "Assay")
# Reading: first convert to H5Seurat format, then load
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 for 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
Because the Seurat object structure differs substantially from that of SingleCellExperiment and AnnData, the conversion is more involved than a simple slot remapping.
{sceasy} (R package, SingleCellExperiment / AnnData / Loom ⇔ Seurat)
{sceasy} also leverages Python for one-step format conversion, but unlike {zellkonverter}, it requires users to set up their own conda environment beforehand.
⚠️ Note: A conda environment must be configured before installing sceasy. The setup steps are:
conda install -c bioconda anndata loompyYou also need to install the following R packages:
if (!requireNamespace("BiocManager", quietly = TRUE)) install.packages("BiocManager") BiocManager::install(c("LoomExperiment", "SingleCellExperiment")) install.packages("reticulate")
The example below demonstrates conversion between AnnData and Seurat. For other format pairs, simply change the from and to arguments accordingly.
# R
# Installation:
devtools::install_github("cellgeni/sceasy")
library(sceasy)
library(reticulate)
use_condaenv('EnvironmentName')
# AnnData to Seurat
sceasy_seurat <- sceasy::convertFormat(h5ad_file, from = "anndata", to = "seurat")
sceasy_seurat
# Seurat to 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 package, for reading H5AD only)
There is also an R package called {anndata}, which provides an interface to Python's AnnData object for reading data. However, this package does not convert AnnData into native R objects, and most downstream analysis tools do not accept this format as input — so further conversion is still required in practice.
Loom
The Loom file format is another HDF5-based specification for omics data. Unlike H5AD, it is not tied to any particular analysis ecosystem, though its structure is similar. Loom files are commonly encountered in RNA velocity analysis output from velocyto. Both R and Python have implementations of the Loom format (e.g., Python's loompy library and the {sceasy} package mentioned above), but the high-level interfaces provided by core ecosystem packages are generally more convenient to use. Since Loom is rarely used in general omics analysis pipelines, we won't go into further detail here.
RDS
The RDS format is an R-specific binary format for serializing arbitrary R objects to disk. It is commonly used to share R analysis results, particularly SingleCellExperiment and Seurat objects.
While RDS is convenient within the R ecosystem, it is not recommended for public data sharing or cross-collaborator exchange due to its complete lack of interoperability with other ecosystems. The HDF5-based formats described above are preferred for such use cases.
⚠️ Note: As mentioned in a recent post, RDS can also be slow to read/write and consume substantial memory and disk space. For reading and writing single-cell data in R, consider using {qs} or {qs2} as a faster alternative.
Emerging Disk Formats
While HDF5 and its derivatives remain the current industry standard, several emerging formats — such as Zarr and TileDB — offer unique advantages for storing very large datasets or multimodal data. It is foreseeable that these formats will develop community specifications and gain broader adoption (anndata already supports Zarr files), but they are beyond the scope of this article. Interested readers are encouraged to explore them further.
03
Closing Thoughts
In this installment, we covered disk-based data conversion methods between R and Python, and walked through several R packages that facilitate these conversions. However, as anyone with hands-on experience knows, these format conversion packages often suffer from a lack of long-term maintenance and struggle to keep pace with updates to their respective core ecosystems (looking at you, H5Seurat). That's why it's important to be familiar with multiple conversion strategies. In the next installment, we'll discuss an alternative approach: in-memory interoperability — stay tuned!
Originally posted on WeChat Official Account: bioinfo_studyroom