---
title: "Saving Bioconductor Objects for Sharing"
author:
  - name: Michael Love
    affiliation: Department of Genetics and Department of Biostatistics, University of North Carolina at Chapel Hill
date: today
format:
  html:
    toc: true
    toc-depth: 3
vignette: >
  %\VignetteIndexEntry{Saving Bioconductor Objects for Sharing}
  %\VignetteEngine{quarto::html}
  %\VignetteEncoding{UTF-8}
---

## Introduction

Bioconductor objects are useful for their shared structure within the project, 
and that they enable rich metadata.

How then should one save a Bioconductor object so that a collaborator 
can load and use it? 
Or so that it persists reliably across time? 

The answer depends on several factors:

- Longevity: will the file need to be easily readable in 1 year, 
  or 5-10 years from now?
- Language interoperability: does the recipient use R, Python, or both?
- Object size: is the object small enough to serialize entirely, or does it
  contain large on-disk arrays?
- Is it easy enough to write out in common formats (BED) or plaintext?
- Reproducibility: should the saved form capture provenance and metadata
  alongside data?

This vignette walks through the main options, their trade-offs, and some
recommendations for common scenarios.

## Acknowledgments

The content of this vignette has been informed by discussions on the
[Bioconductor community Zulip](https://chat.bioconductor.org/), including
contributions from Kevin Rue-Albrecht, Johannes Rainer, Lori Shepherd, Jayaram
Kancherla, Aaron Lun, Hervé Pagès, Laurent Gatto, Vince Carey, Sean Davis,
Robert Castelo, Hugo Gruson, and Luke Zappia.

## Version Info

```{r}
#| label: version-info
#| echo: false
#| results: asis
suppressPackageStartupMessages(library(BiocManager))
```

<p>
**R version**: `r R.version.string`
<br />
**Bioconductor version**: `r BiocManager::version()`
<br />
**Package version**: `r packageVersion("savingBiocObjects")`
</p>

## R serialization

### saveRDS and readRDS

The simplest approach is R's built-in binary serialization. `saveRDS()` saves
a single object to an RDS file (the acronym is not defined in the man pages
but likely stands for "R Data Serialization"); `save()` bundles one or more
named objects into an `.RData` (or `.rda`) file.

```{r}
#| label: rds-save-load
#| message: false
library(SummarizedExperiment)

se <- SummarizedExperiment(
  assays = list(counts = matrix(1:12, nrow = 3)),
  colData = DataFrame(
    condition = c("A", "A", "B", "B"), 
    row.names=1:4
    ),
  rowData = DataFrame(
    gene = c("gene1","gene2","gene3"), 
    row.names=1:3
    )
)

# Single object
tmp_rds <- tempfile(fileext = ".rds")
saveRDS(se, file = tmp_rds)
se_from_rds <- readRDS(tmp_rds)
se_from_rds

# Multiple objects in one file
tmp_rda <- tempfile(fileext = ".RData")
save(se, file = tmp_rda)
load(tmp_rda)  # restores 'se' by name into the current environment
```

Prefer `saveRDS()` for most use cases. `save()` / `load()` silently overwrites
any object in the calling environment that shares a name, which is a common
source of confusion. The main remaining use case for `save()` is `.rda` data
files shipped inside R packages (under `data/`).

Advantages:

- Zero setup — works with any R object.
- Compression applied automatically (gzip by default).
- Round-trips perfectly: the loaded object is identical to the saved one.

Disadvantages:

- R-only: Python or other languages cannot read these files without a bridge
  library.
- Large objects are loaded entirely into RAM; there is no lazy / on-disk access.

When to use it: quick sharing between R users on the same project, saving
intermediate objects in a pipeline, anything under ~1 GB.

If your workflow downloads RDS files from a remote URL, consider using
[`BiocFileCache`](https://bioconductor.org/packages/BiocFileCache) to cache
them locally so they are only fetched once.

### Reading RDS files in Python

[BiocPy](https://biocpy.github.io/) is a Python ecosystem that brings
Bioconductor's core data structures to Python, including
[BiocFrame](https://github.com/BiocPy/BiocFrame),
[IRanges](https://github.com/BiocPy/IRanges),
[GenomicRanges](https://github.com/BiocPy/GenomicRanges),
[SummarizedExperiment](https://github.com/BiocPy/SummarizedExperiment),
[SingleCellExperiment](https://github.com/BiocPy/SingleCellExperiment), and
[MultiAssayExperiment](https://github.com/BiocPy/MultiAssayExperiment).

Python users can read RDS files with the
[rds2py](https://github.com/BiocPy/rds2py) package from the BiocPy ecosystem.
Standard R types map to NumPy/SciPy equivalents (e.g. numeric vectors become
`numpy.ndarray`), and Bioconductor classes such as `SummarizedExperiment`,
`SingleCellExperiment`, `GRanges`, and `MultiAssayExperiment` are converted to
their BiocPy counterparts. For unrecognised S4 classes the object falls back
to a dictionary so no data is lost. Support for writing RDS files from Python
is also in development.

### Cross-release stability

There is no guarantee that an S4 object serialized today will work with older
versions of Bioconductor. The most common reasons are that the S4 class is not
defined at all in the earlier version, or that it exists but its definition has
changed — new slots added, slots renamed or removed, or infrastructure moved
between packages.

A concrete example: in Bioconductor 3.22, `Seqinfo` was moved out of
`GenomicRanges` into its own package. A `GRanges` serialized under BioC ≥ 3.22
can still be loaded on a machine running BioC < 3.22, and many basic operations
work — `show()`, `seqnames()`, `ranges()`, `mcols()`, and `[` among them. But
operations that require the `Seqinfo` package (such as `shift()` or `reduce()`)
will fail with a confusing "package not available" error. The practical advice
is to upgrade to BioC ≥ 3.22. The reverse was also true: older `GRanges`
objects loaded into a newer session sometimes required `updateObject()` to
migrate the internal representation:

```{r}
#| label: update-object
#| eval: false
# not evaluated — requires an object saved under an older Bioconductor release
gr <- readRDS("old_granges.rds")
gr <- updateObject(gr, verbose = TRUE)
```

The general lesson is: if you can save your data in a format that is not
tied to a particular version of Bioconductor — or better, not tied to R at
all — you should. For a `GRanges`, for instance, a simple TSV is often
enough:

```{r}
#| label: granges-tsv
#| message: false
library(GenomicRanges)

gr <- GRanges(
  seqnames = "chr1",
  ranges = IRanges(start = c(100, 200, 300), width = 50),
  seqinfo = Seqinfo(seqnames = "chr1", seqlengths = 248956422,
                    isCircular = FALSE, genome = "hg38")
)
names(gr) <- c("peak1", "peak2", "peak3")
gr$score  <- c(500, 800, 300)   # standard BED score column
gr$log2fc <- c(1.2, -0.5, 2.1) # extra metadata column

tmp_tsv <- tempfile(fileext = ".tsv")
write.table(as.data.frame(gr), tmp_tsv, sep = "\t", quote = FALSE)

gr_from_tsv <- makeGRangesFromDataFrame(
  read.table(tmp_tsv, header = TRUE, sep = "\t"),
  keep.extra.columns = TRUE
)
gr_from_tsv
```

This round-trip survives any Bioconductor version and is readable from Python
or the command line. The `alabaster` ecosystem (described below) applies the
same principle more systematically and with better support for complex objects,
using HDF5 and JSON as the underlying storage formats.

## HDF5-backed storage

### Saving with HDF5Array

For large assay matrices (e.g., single-cell count matrices with millions of
cells), it is impractical to hold the entire object in RAM. The
[`HDF5Array`](https://bioconductor.org/packages/HDF5Array) package provides
array classes backed by HDF5 files, enabling lazy loading and out-of-memory
computation.

```{r}
#| label: hdf5-save-load
#| message: false
library(HDF5Array)

tmp_hdf5 <- tempfile()
saveHDF5SummarizedExperiment(se, dir = tmp_hdf5, replace = TRUE)

# Assay data remains on disk until accessed
se_from_hdf5 <- loadHDF5SummarizedExperiment(tmp_hdf5)
se_from_hdf5
```

The saved directory contains an HDF5 file with the assay data and an RDS
file for the non-assay metadata.

Advantages:

- Assay data is stored on disk; only the chunks you access are read into RAM.
- HDF5 is a widely used binary format with readers in Python (`h5py`,
  `anndata`), Julia, C/C++, and more.
- Good for objects with tens of gigabytes of assay data.

Disadvantages:

- The output is a directory, not a single file, which complicates transfer
  (use `tar` or `zip` before sharing).
- The RDS envelope for metadata is still R-specific.
- Write performance can be slower than `saveRDS()` for small objects.

When to use it: large single-cell or spatial datasets where you want
on-disk access; workflows shared between R users who need memory efficiency.

### SummarizedExperiment and AnnData

For single-cell workflows that move between R and Python, the `.h5ad` format
used by scanpy and related tools is often the most convenient path when
collaborators are working in Python. Like `HDF5Array`, `.h5ad` is an
HDF5-based format. The trade-off relative to alabaster is that `.h5ad` is
AnnData-specific rather than a general Bioconductor serialization format.

The recommended starting point is the
[anndataR](https://anndatar.scverse.org/) package, a more recent and complete
implementation that is actively maintained as part of the scverse ecosystem:

```{r}
#| label: anndataR
#| message: false
library(anndataR)
library(SingleCellExperiment)

sce <- as(se, "SingleCellExperiment")

tmp_h5ad <- tempfile(fileext = ".h5ad")
write_h5ad(sce, path = tmp_h5ad)

sce_from_h5ad <- read_h5ad(tmp_h5ad, as = "SingleCellExperiment")
sce_from_h5ad
```

The [zellkonverter](https://bioconductor.org/packages/zellkonverter) package
is an alternative that also converts directly between `SingleCellExperiment`
and `.h5ad`. It remains actively maintained and has advantages in some cases,
though at the cost of managing a Python environment via `basilisk`.

```{r}
#| label: zellkonverter
#| eval: false
# not evaluated — zellkonverter installs a full Python environment via basilisk
# on first use, which takes too long in CI
library(zellkonverter)

writeH5AD(sce, file = "sce.h5ad")

sce_from_h5ad <- readH5AD("sce.h5ad")
sce_from_h5ad
```

## The alabaster ecosystem

### saveObject and readObject

The [`alabaster`](https://bioconductor.org/packages/alabaster.base) family of
packages is part of the broader [ArtifactDB](https://github.com/ArtifactDB)
project, which provides a multi-language system for storing and retrieving
analysis-ready Bioconductor objects. The core idea
is to save objects as directories of standard files (HDF5, JSON, CSV) whose
format is defined by explicit, versioned specifications — meaning the saved
form is readable without R, and can evolve over time without breaking
previously saved objects.

```{r}
#| label: alabaster-save-load
#| message: false
library(alabaster.base)
library(alabaster.se)

tmp_alabaster <- tempfile()
saveObject(se, path = tmp_alabaster)

se_from_alabaster <- readObject(tmp_alabaster)
se_from_alabaster
```

The `alabaster` umbrella package pulls in support for the most common
Bioconductor classes. Individual sub-packages cover specific classes:
`alabaster.se` for `SummarizedExperiment`, `alabaster.sce` for
`SingleCellExperiment`, and so on.

### Validation with takane

A key part of the ArtifactDB design is that saved directories can be
independently validated against the format specification. This is handled by
[takane](https://github.com/ArtifactDB/takane), a C++ library that maintains
separate, versioned specifications for 30+ Bioconductor object types. Calling
`takane::validate()` on a saved directory checks that all files conform to the
expected layout and types, which means a collaborator or downstream tool can
verify the integrity of a saved object without needing to load it into R.
This makes alabaster directories suitable for deposition in data repositories
where format conformance needs to be auditable.

The Python counterpart to `alabaster` is the
[dolomite](https://github.com/ArtifactDB/dolomite-base) family of packages,
which reads and writes the same on-disk format. An object saved with
`alabaster` in R can be read with `dolomite` in Python, and vice versa, with
no conversion step.

Advantages:

- Truly language-agnostic: Python readers exist via the `dolomite` family of
  packages, enabling seamless R ↔ Python interoperability.
- Built on open standards (HDF5, JSON); inspectable without R.
- Versioned format specifications with independent validation via takane.
- A good choice for archives or data portals.

Disadvantages:

- Newer ecosystem; not all Bioconductor classes have `alabaster` support yet.
- Requires installing the relevant `alabaster.*` sub-package for each class.
- Like HDF5Array, the output is a directory.

When to use it: archival storage, data portal submissions, cross-language
workflows, or any situation where you want the saved format to be readable
without R.

## BED format for ranges

### Writing BED files

When the object is a `GRanges` or similar ranges object and the goal is
interoperability with other tools (genome browsers, Python, command-line
utilities), exporting to BED format is often more useful than R-specific
serialization. Our `gr` has range names, a standard BED score column, and
an extra metadata column `log2fc`:

```{r}
#| label: show-gr
gr
```

Both `rtracklayer` and `plyranges` can write BED files:

```{r}
#| label: bed-write
#| message: false
library(rtracklayer)
library(plyranges)

tmp_bed <- tempfile(fileext = ".bed")
export(gr, tmp_bed)

tmp_bed2 <- tempfile(fileext = ".bed")
write_bed(gr, tmp_bed2)
```

When reading back, the standard `score` column is preserved, but `log2fc` is
silently dropped — BED has no mechanism to carry arbitrary metadata columns.
Range names are stored in the BED name field but come back as a `$name`
metadata column rather than as R names on the object; restore them manually:

```{r}
#| label: bed-read
gr_rtracklayer <- import(tmp_bed)
names(gr_rtracklayer) <- gr_rtracklayer$name
gr_rtracklayer$name <- NULL
gr_rtracklayer

gr_plyranges <- read_bed(tmp_bed2)
names(gr_plyranges) <- gr_plyranges$name
gr_plyranges$name <- NULL
gr_plyranges
```

### Saving metadata columns

If preserving all metadata columns is the priority and BED compatibility is
not required, the simplest approach is the TSV round-trip shown earlier:
`write.table(as.data.frame(gr), ...)` followed by
`makeGRangesFromDataFrame(..., keep.extra.columns = TRUE)` restores all mcols
in one step without a sidecar.

When you do need a BED file (e.g. for a genome browser or a tool that expects
BED input), extra mcols can be preserved by writing them to a sidecar file.
Here using plyranges to read the BED back, then reattaching from the sidecar:

```{r}
#| label: bed-mcols-sidecar
tmp_meta <- tempfile(fileext = ".tsv")
write.table(
  data.frame(name = names(gr), log2fc = gr$log2fc),
  tmp_meta, sep = "\t", quote = FALSE, row.names = FALSE
)

gr_restored <- read_bed(tmp_bed2)
meta <- read.table(tmp_meta, header = TRUE, sep = "\t")
gr_restored$log2fc <- meta$log2fc
gr_restored
```

### Saving Seqinfo separately

BED files do not store chromosome lengths or genome build information, so
`Seqinfo` is silently dropped on export. To preserve it, write it out
alongside the BED file and restore it on load:

```{r}
#| label: seqinfo-save-restore
tmp_seqinfo <- tempfile(fileext = ".csv")
write.csv(as.data.frame(seqinfo(gr)), tmp_seqinfo)

df <- read.csv(tmp_seqinfo, row.names = 1)
si <- Seqinfo(
  seqnames   = rownames(df),
  seqlengths = as.integer(df$seqlengths),
  isCircular = as.logical(df$isCircular),
  genome     = as.character(df$genome)
)
si
```

## Propagating object metadata

Object-level metadata stored in `metadata(object)` — things like processing
parameters, provenance notes, or experiment descriptors — is lost in any
format that only encodes the ranges or assay data. This applies whether you
are writing a BED file, an HDF5 matrix, or any other non-R format. Write it
to a JSON sidecar file so it travels with the data:

```{r}
#| label: metadata-json
#| message: false
library(jsonlite)

metadata(se) <- list(
  timestamp = as.POSIXct("2020-01-01 12:00:00", tz = "UTC"),
  pipeline = "v2.1",
  n_samples = 4L
)

tmp_json <- tempfile(fileext = ".json")
writeLines(toJSON(metadata(se), pretty = TRUE, auto_unbox = TRUE), tmp_json)

metadata(se_from_rds) <- fromJSON(tmp_json)
metadata(se_from_rds)
```

`toJSON` handles simple R types (lists, vectors, data frames) well, but
complex objects (S4 instances, environments) need to be simplified or omitted
before serializing.

## Summary and recommendations

| Scenario | Recommended approach |
|---|---|
| Quick sharing between R users, same Bioc release | `saveRDS()` |
| Loading an object from an older Bioc release | `updateObject()` after `readRDS()` |
| Large assay matrices, R-only | `saveHDF5SummarizedExperiment()` |
| SingleCellExperiment ↔ Python AnnData | `anndataR::write_h5ad()` / `read_h5ad()` (or `zellkonverter` for Python-env integration) |
| Cross-language (R + Python), general | `alabaster::saveObject()` |
| Long-term archive / data portal | `alabaster::saveObject()` |

In most new projects we recommend defaulting to `saveRDS()` for convenience
and upgrading to `alabaster` when cross-language access or archival stability
becomes a priority. Be aware that any R-serialized Bioconductor object may
require `updateObject()` when loaded under a different Bioconductor release.

## Session info

```{r}
#| label: session-info
sessionInfo()
```