## ----setup, include=FALSE-----------------------------------------------------
knitr::opts_chunk$set(collapse = TRUE, comment = "#>",
                      error = FALSE, warning = FALSE, message = FALSE)
has_airway <- requireNamespace("airway", quietly = TRUE)

## ----install, eval=FALSE------------------------------------------------------
# if (!require("BiocManager"))
#     install.packages("BiocManager")
# BiocManager::install("BiocDuckDB")

## ----load---------------------------------------------------------------------
library(BiocDuckDB)
library(SummarizedExperiment)

## ----airway-check, echo=FALSE, results='asis'---------------------------------
if (!has_airway) {
    cat("_The `airway` package is not installed; the code below is shown but",
        "not evaluated._\n")
}

## ----roundtrip, eval=has_airway-----------------------------------------------
data(airway, package = "airway")

se_path <- file.path(tempdir(), "airway_se")
writeParquet(airway, se_path)

airway_ddb <- readParquet(se_path)
class(airway_ddb)
dim(airway_ddb)

## ----layout, eval=has_airway--------------------------------------------------
head(list.files(se_path, recursive = TRUE), 10)

## ----ingest, eval=FALSE-------------------------------------------------------
# resources <- c(
#     writeParquet(features, file.path(dir, "features"), dimension = "feature"),
#     writeParquet(assays,   dir,                        dimension = "crossed"))
# writeDatapackage("summarized_experiment", resources, dir)

## ----attach, eval=FALSE-------------------------------------------------------
# mat <- DuckDBMatrix(coo_path, datacol = "value",
#                     keycols = c("__feature__", "__sample__"))

## ----operations, eval=has_airway----------------------------------------------
assay(airway_ddb, "counts")
colData(airway_ddb)[, 1:3]

## subset on disk, then summarize
sub <- airway_ddb[1:1000, 1:4]
colSums(assay(sub, "counts"))

## ----memory, eval=has_airway--------------------------------------------------
c(in_memory = format(object.size(airway), units = "MB"),
  duckdb    = format(object.size(airway_ddb), units = "MB"))

## ----ranges, eval=has_airway--------------------------------------------------
rr <- rowRanges(airway_ddb)
class(rr)
elementNROWS(rr)[1:8]      # exons per gene, queried from disk

## ----sce----------------------------------------------------------------------
library(SingleCellExperiment)
library(Matrix)

set.seed(1L)
counts <- Matrix(rpois(2000 * 500, lambda = 0.3), nrow = 2000, ncol = 500,
                 sparse = TRUE)
rownames(counts) <- paste0("Gene", seq_len(nrow(counts)))
colnames(counts) <- paste0("Cell", seq_len(ncol(counts)))
sce <- SingleCellExperiment(assays = list(counts = counts))

sce_path <- file.path(tempdir(), "demo_sce")
writeParquet(sce, sce_path)
sce_ddb <- readParquet(sce_path)

## QC filtering stays on disk (SQL), then summarize
totals <- colSums(assay(sce_ddb, "counts"))
keep <- sce_ddb[, totals > median(totals)]
dim(keep)

## ----pattern, eval=FALSE------------------------------------------------------
# library(scran)
# library(scater)
# 
# sce_ddb <- readParquet("data/raw_counts")          # DuckDB-backed, low RAM
# 
# ## 1. filter on disk (SQL-optimized)
# keep <- colSums(assay(sce_ddb, "counts")) > 1000
# sce_ddb <- sce_ddb[, keep]
# detected <- rowSums(assay(sce_ddb, "counts") > 0)
# sce_ddb <- sce_ddb[detected >= 10, ]
# 
# ## 2. realize the filtered subset into memory
# sce <- as(sce_ddb, "SingleCellExperiment")
# 
# ## 3. standard in-memory pipeline
# sce <- logNormCounts(sce)
# dec <- modelGeneVar(sce)
# sce <- runPCA(sce[getTopHVGs(dec, n = 2000), ])
# 
# ## 4. persist results back to Parquet
# writeParquet(sce, "data/processed")

## ----mae, eval=FALSE----------------------------------------------------------
# library(MultiAssayExperiment)
# 
# writeParquet(rna_sce,     file.path(mae_path, "rna"))
# writeParquet(protein_se,  file.path(mae_path, "protein"))
# 
# mae <- MultiAssayExperiment(experiments = list(
#     rna     = readParquet(file.path(mae_path, "rna")),
#     protein = readParquet(file.path(mae_path, "protein"))))

## ----sessioninfo--------------------------------------------------------------
sessionInfo()