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


## -----------------------------------------------------------------------------
#| 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


## -----------------------------------------------------------------------------
#| 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)


## -----------------------------------------------------------------------------
#| 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


## -----------------------------------------------------------------------------
#| 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


## -----------------------------------------------------------------------------
#| 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


## -----------------------------------------------------------------------------
#| 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


## -----------------------------------------------------------------------------
#| 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


## -----------------------------------------------------------------------------
#| label: show-gr
gr


## -----------------------------------------------------------------------------
#| 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)


## -----------------------------------------------------------------------------
#| 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


## -----------------------------------------------------------------------------
#| 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


## -----------------------------------------------------------------------------
#| 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


## -----------------------------------------------------------------------------
#| 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)


## -----------------------------------------------------------------------------
#| label: session-info
sessionInfo()