---
title: "OmicsLake Comprehensive Guide"
author: "OmicsLake Development Team"
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
vignette: >
    %\VignetteIndexEntry{OmicsLake Comprehensive Guide}
    %\VignetteEngine{knitr::rmarkdown}
    %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(
    collapse = TRUE,
    comment = "#>",
    eval = isTRUE(getOption("omicslake.vignette.eval", FALSE)),
    fig.width = 7,
    fig.height = 5
)
```

# Introduction

OmicsLake is an R package that provides a comprehensive version management and
dependency tracking system for bioinformatics data analysis.  This vignette
introduces all of the major features with examples.

## Before you run all chunks

This vignette is configured for non-evaluated chunks during package checks.
To execute all chunks locally:

```r
options(omicslake.vignette.eval = TRUE)
```

If you prefer a ready-to-run script for a canonical RNA-seq style flow, use:

```bash
bash tools/run_demo_count_edger_limma_voom_ora.sh
```

For representative layer-by-layer use cases
(SE/SCE/MAE/Spectra/QFeatures/MsExperiment), see:

- `vignettes/omicslake_layer_use_cases_EN.Rmd`

## Main package features

- **Project management**: initialize and manage projects based on DuckDB
- **Data storage**: save tables and R objects
- **Version management**: tagging, history management, version comparison
- **Dependency tracking**: automatically track dependencies in data analysis flows
- **Visualization**: visualize dependency graphs
- **Bioconductor integration**: support for *SE/SCE/MAE/Spectra/QFeatures/MsExperiment* workflows

# 1. Basic operations

## 1.1 Project initialization

```{r init}
library(OmicsLake)

# Initialize a project
ol_init("de_analysis")
```

`ol_init()` creates a project with the specified name and initializes a DuckDB
database.

## 1.2 Writing and reading tables

```{r tables}
# Create sample data
set.seed(123)
raw_counts <- data.frame(
    gene_id = paste0("GENE", 1:100),
    sample1 = rpois(100, 100),
    sample2 = rpois(100, 100),
    sample3 = rpois(100, 120),
    sample4 = rpois(100, 120)
)

# Save as a table
ol_write("raw_counts", raw_counts, mode = "create")

# Read the table
counts <- ol_read("raw_counts")
head(counts)

# List tables
ol_list_tables()
```

**How to interpret the output**: `ol_list_tables()` returns the names of all
tables stored in the current project.

## 1.3 Advanced analysis with SQL queries (`ol_query`)

Using `ol_query()` makes it possible to perform advanced data analysis using
all of DuckDB’s functionality, such as JOINs, aggregations, and window
functions.

```{r query_basic}
# Basic query
result <- ol_query("SELECT * FROM raw_counts WHERE sample1 > 100")
head(result)

# Aggregation query
summary <- ol_query("\
    SELECT
    COUNT(*) as total_genes,
    AVG(sample1) as avg_sample1,
    MAX(sample1) as max_sample1
    FROM raw_counts
")
print(summary)
```

### JOINs and subqueries

```{r query_join}
# Join tables (example: join with gene annotation data)
ol_write("gene_annotations", data.frame(
    gene_id = paste0("gene", 1:5),
    name = paste0("GeneA", 1:5),
    chromosome = c("chr1", "chr2", "chr1", "chr3", "chr2")
))

joined <- ol_query("\
    SELECT r.gene_id, r.sample1, a.name, a.chromosome
    FROM raw_counts r
    JOIN gene_annotations a ON r.gene_id = a.gene_id
    WHERE r.sample1 > 50
")
head(joined)
```

### Lazy evaluation and integration with `dplyr`

```{r query_lazy}
# collect=FALSE for lazy evaluation
if (requireNamespace("dplyr", quietly = TRUE)) {
    lazy_tbl <- ol_query("SELECT * FROM raw_counts", collect = FALSE)

    # Process with dplyr
    result <- lazy_tbl %>%
        dplyr::filter(sample1 > 100) %>%
        dplyr::select(gene_id, sample1, sample2) %>%
        dplyr::arrange(dplyr::desc(sample1)) %>%
        dplyr::collect()

    head(result)
}
```

**Note**: Table names can be referred to without the `ol.` prefix (for example,
`genes` is fine, but `ol.genes` can also be used).

## 1.4 Advanced analysis with aggregation and window functions

By using the Phase 2 functions, common aggregation and window function
operations become easier.  These functions leverage DuckDB’s advanced analytics
features to efficiently process gene expression data.

### Calculating gene expression statistics

```{r aggregate_stats}
# Statistics across all samples
overall_stats <- ol_aggregate("raw_counts",
    mean_sample1 = list(func = "avg", col = "sample1"),
    median_sample1 = list(func = "median", col = "sample1"),
    sd_sample1 = list(func = "stddev", col = "sample1")
)
print(overall_stats)

# Statistics by category (e.g., high vs low expression in sample1)
# First add a category for high/low expression
ol_query("\
    CREATE OR REPLACE TABLE raw_counts_categorized AS
    SELECT *,
    CASE WHEN sample1 > 100 THEN 'high' ELSE 'low' END as expr_category
    FROM raw_counts
")

grouped_stats <- ol_aggregate("raw_counts_categorized",
    group_by = "expr_category",
    count = list(func = "count", col = "*"),
    mean_s1 = list(func = "avg", col = "sample1"),
    mean_s2 = list(func = "avg", col = "sample2")
)
print(grouped_stats)
```

### Ranking genes

```{r ranking}
# Rank genes by expression
ranked_genes <- ol_add_rank("raw_counts",
    rank_by = "sample1",
    method = "row_number",
    descending = TRUE,
    as_column = "expression_rank"
)
head(ranked_genes)

# Get the top 10 genes
top_genes <- ol_top_n("raw_counts",
    n = 10,
    order_by = "sample1",
    descending = TRUE
)
print(top_genes)

# Top 5 genes per sample
# Convert to long format first
ol_query("\
    CREATE OR REPLACE TABLE raw_counts_long AS
    SELECT gene_id, 'sample1' as sample, sample1 as expression FROM raw_counts
    UNION ALL
    SELECT gene_id, 'sample2' as sample, sample2 as expression FROM raw_counts
    UNION ALL
    SELECT gene_id, 'sample3' as sample, sample3 as expression FROM raw_counts
    UNION ALL
    SELECT gene_id, 'sample4' as sample, sample4 as expression FROM raw_counts
")

top_per_sample <- ol_top_n("raw_counts_long",
    n = 5,
    order_by = "expression",
    partition_by = "sample",
    descending = TRUE
)
print(top_per_sample)
```

### Moving average and cumulative sum

```{r moving_cumulative}
# Moving average (smoothing)
# Sort by gene ID and compute a 3-gene moving average
smoothed <- ol_moving_avg("raw_counts",
    column = "sample1",
    window_size = 3,
    order_by = "gene_id",
    as_column = "sample1_smoothed"
)
head(smoothed, 10)

# Larger window size
smoothed_5 <- ol_moving_avg("raw_counts",
    column = "sample1",
    window_size = 5,
    order_by = "gene_id",
    as_column = "sample1_ma5"
)
head(smoothed_5, 10)

# Cumulative sum (sorted by gene ID)
cumulative <- ol_cumulative_sum("raw_counts",
    column = "sample1",
    order_by = "gene_id",
    as_column = "sample1_cumsum"
)
head(cumulative, 10)
tail(cumulative)

# Cumulative sum per sample (using the long-format data)
cumulative_per_sample <- ol_cumulative_sum("raw_counts_long",
    column = "expression",
    partition_by = "sample",
    order_by = "gene_id",
    as_column = "cumulative_expr"
)
head(cumulative_per_sample, 20)
```

### Combining with lazy evaluation

```{r lazy_aggregation}
# Aggregation functions also support lazy evaluation
if (requireNamespace("dplyr", quietly = TRUE)) {
    # Compute statistics lazily
    lazy_stats <- ol_aggregate("raw_counts",
        mean_s1 = list(func = "avg", col = "sample1"),
        collect = FALSE
    )

    # Further process with dplyr
    result <- lazy_stats %>%
        dplyr::collect()
    print(result)

    # Combine with ranking
    lazy_ranked <- ol_add_rank("raw_counts",
        rank_by = "sample1",
        as_column = "rank",
        collect = FALSE
    )

    top_ranked <- lazy_ranked %>%
        dplyr::filter(rank <= 10) %>%
        dplyr::select(gene_id, sample1, rank) %>%
        dplyr::arrange(rank) %>%
        dplyr::collect()
    print(top_ranked)
}
```

By using these functions, you can easily perform common analysis patterns
without writing complex SQL.

## 1.5 Saving and loading R objects

```{r objects}
# Compute normalization factors (example)
norm_factors <- list(
    method = "TMM",
    factors = runif(100, 0.8, 1.2)
)

# Save as an R object
ol_save("norm_factors", norm_factors)

# Load the object
loaded_factors <- ol_read_object("norm_factors")
str(loaded_factors)

# List objects
ol_list_objects()
```

**How to interpret the output**: `ol_list_objects()` returns the names of all R
objects that have been saved.

# 2. Tracking dependencies

## 2.1 Saving data with specified dependencies

```{r dependencies}
# Compute normalized counts (specify dependencies via depends_on)
normalized_counts <- raw_counts
for (i in 2:5) {
    normalized_counts[[i]] <- normalized_counts[[i]] * norm_factors$factors
}

# Save with dependencies
ol_write("normalized_counts", normalized_counts,
    depends_on = c("raw_counts", "norm_factors")
)

# Save DE parameters
de_params <- list(
    method = "DESeq2",
    alpha = 0.05,
    lfc_threshold = 1.0
)
ol_save("de_params", de_params)

# Compute DE results (simplified example)
de_results <- data.frame(
    gene_id = paste0("GENE", 1:20),
    log2FC = rnorm(20, 0, 2),
    pvalue = runif(20, 0, 0.1),
    padj = runif(20, 0, 0.1)
)

# Save DE results with dependencies
ol_save("de_results", de_results,
    depends_on = c("normalized_counts", "de_params")
)
```

## 2.2 Checking dependencies

```{r view_dependencies}
# Upstream dependencies (what this object depends on)
upstream <- ol_get_dependencies("de_results", direction = "upstream")
print(upstream)

# Downstream dependencies (what depends on this object)
downstream <- ol_get_dependencies("raw_counts", direction = "downstream")
print(downstream)

# Show the complete lineage tree
lineage <- ol_show_lineage("de_results", direction = "upstream")
print(lineage)
```

**How to interpret the output**:

- `ol_get_dependencies()`: returns the direct dependencies as a `data.frame`
- `ol_show_lineage()`: recursively explores all ancestors/descendants and
- returns the complete lineage tree

## 2.3 Importing and exporting Parquet files

OmicsLake uses DuckDB’s high‑performance Parquet support to efficiently export
tables to Parquet files and import from Parquet files.

### Parquet export

```{r parquet_export}
# Basic export (default is Snappy compression)
ol_export_parquet("genes", "genes.parquet")

# Export with zstd compression (high compression ratio)
ol_export_parquet("genes", "genes_zstd.parquet",
    compression = "zstd",
    compression_level = 3
)

# Specify row group size
ol_export_parquet("genes", "genes_optimized.parquet",
    compression = "zstd",
    row_group_size = 50000
)

# Uncompressed export (fastest)
ol_export_parquet("genes", "genes_uncompressed.parquet",
    compression = "uncompressed"
)
```

### Parquet import

```{r parquet_import}
# Create a table from a Parquet file
ol_import_parquet("genes.parquet", "imported_genes", mode = "create")

# Overwrite an existing table
ol_import_parquet("genes_new.parquet", "genes", mode = "overwrite")

# Append data to an existing table
ol_import_parquet("genes_batch2.parquet", "genes", mode = "append")

# Import multiple Parquet files at once
ol_import_parquet(c("genes_part1.parquet", "genes_part2.parquet"),
    "all_genes",
    mode = "create"
)

# Record dependencies when importing
ol_import_parquet("processed_genes.parquet",
    "final_genes",
    depends_on = "raw_genes",
    mode = "create"
)
```

### Performance tips

**Choosing a compression algorithm:**

- `snappy`: balanced and recommended by default
- `zstd`: highest compression ratio, also good read/write speeds
- `lz4`: fastest compression/decompression
- `brotli`: highest compression ratio but slow
- `uncompressed`: no compression, fastest but large file size

**Row group size:**

- Small values (10,000–50,000): good memory efficiency, advantageous for
- filtering
- Large values (100,000–500,000): higher compression ratio, advantageous for
- scan performance
- Default (100,000): well‑balanced

```{r parquet_performance}
# Example optimization for a large dataset
ol_export_parquet("large_dataset",
    "large_dataset.parquet",
    compression = "zstd",
    compression_level = 1, # fastest zstd
    row_group_size = 250000
)
```


# 3. Version management

## 3.1 Tagging and labeling

```{r tagging}
# Tag an object
ol_tag_object("de_results", "baseline_analysis")

# Label the entire project (tags all tables/objects)
ol_label("experiment_v1")

# List tags
ol_list_tags()

# List project labels
ol_list_labels()
```

## 3.2 Creating multiple versions

```{r multiple_versions}
# Change parameters and re‑run analysis
de_params_v2 <- list(
    method = "DESeq2",
    alpha = 0.01,
    lfc_threshold = 1.5
)
ol_save("de_params", de_params_v2)

# Recompute DE results with new parameters
de_results_v2 <- data.frame(
    gene_id = paste0("GENE", 1:10),
    log2FC = rnorm(10, 0, 2.5),
    pvalue = runif(10, 0, 0.01),
    padj = runif(10, 0, 0.01)
)
ol_save("de_results", de_results_v2,
    depends_on = c("normalized_counts", "de_params")
)

# Tag the new version
ol_tag_object("de_results", "strict_analysis")

# Try a different analysis method
de_params_edger <- list(
    method = "edgeR",
    alpha = 0.05,
    lfc_threshold = 1.0
)
ol_save("de_params", de_params_edger)

de_results_edger <- data.frame(
    gene_id = paste0("GENE", c(1:12, 25:30)),
    log2FC = rnorm(18, 0, 2),
    pvalue = runif(18, 0, 0.1),
    padj = runif(18, 0, 0.1)
)
ol_save("de_results", de_results_edger,
    depends_on = c("normalized_counts", "de_params")
)
ol_tag_object("de_results", "edger_analysis")
```

## 3.3 Listing and comparing versions

```{r version_comparison}
# List all versions
versions <- ol_list_object_versions("de_results")
print(versions)

# Compare versions
comparison <- ol_compare_versions("de_results")
print(comparison)

# Compare specific tags
tag_comparison <- ol_compare_versions("de_results",
    versions = c("baseline_analysis", "strict_analysis")
)
print(tag_comparison)
```

**How to interpret the output**:

- `version_ts`: timestamp of the version
- `tags`: tags applied to that version
- `size_bytes`: size of the object (bytes)
- `dependencies`: dependencies of that version
- `size_change`: size change from the previous version
- `time_since_previous`: elapsed time since the previous version
- `deps_added`: dependencies added
- `deps_removed`: dependencies removed

## 3.4 Loading a specific version

```{r load_versions}
# Load by tag
baseline_results <- ol_read_object(
    "de_results",
    ref = "@tag(baseline_analysis)"
)
strict_results <- ol_read_object("de_results", ref = "@tag(strict_analysis)")

cat("Baseline analysis:", nrow(baseline_results), "genes\n")
cat("Strict analysis:", nrow(strict_results), "genes\n")
```

# 4. Commits and history management

## 4.1 Creating commits

```{r commits}
# Create a commit at a milestone
commit_id <- ol_commit(
    note = "Completed initial DE analysis with three methods",
    params = list(
        methods = c("DESeq2_baseline", "DESeq2_strict", "edgeR"),
        date = as.character(Sys.Date())
    )
)

cat("Commit ID:", commit_id, "\n")
```

## 4.2 Viewing history

```{r history}
# Show commit history
commits <- ol_log_commits(n = 10)
print(commits)

# History of a specific table
table_log <- ol_log("raw_counts")
print(table_log)
```

# 5. Visualizing dependencies

## 5.1 Creating dependency graphs

```{r visualization, fig.width=8, fig.height=6}
# The igraph package is required
if (requireNamespace("igraph", quietly = TRUE)) {
    # Visualize upstream dependencies
    ol_plot_lineage("de_results",
        direction = "upstream",
        layout = "sugiyama",
        main = "DE Results - Upstream Dependencies"
    )

    # Visualize dependencies in both directions
    ol_plot_lineage("normalized_counts",
        direction = "both",
        layout = "tree",
        main = "Normalized Counts - Full Lineage"
    )
}
```

**How to interpret the visualization**:

- Blue nodes: tables
- Orange nodes: R objects
- Red nodes: the specified focal node
- Arrows: direction of dependencies (parent → child)

# 6. Restoring a project

## 6.1 Restoring state using labels

```{r checkout}
# Add a new table
filtered_results <- de_results[de_results$padj < 0.05, ]
ol_write("filtered_de", filtered_results)

# Revert to an earlier label state
ol_checkout("experiment_v1")

# The filtered_de table will no longer exist
current_tables <- ol_list_tables()
cat(
    "Tables after checkout:",
    paste(current_tables$table_name, collapse = ", "),
    "\n"
)
```

`ol_checkout()` restores the state of all tables and objects at the specified
label.

# 7. Advanced features

## 7.1 Filtered data loading

```{r fread}
# Recreate raw_counts (may have been removed by checkout)
ol_write("raw_counts", raw_counts, mode = "overwrite")

# Select specific columns and read
selected <- ol_fread("raw_counts",
    select = c("gene_id", "sample1", "sample2"),
    nrows = 10
)
head(selected)

# Filtering by condition
filtered <- ol_fread("raw_counts",
    filter = "sample1 > 100"
)
head(filtered)
```

## 7.2 Loading with lazy evaluation

```{r lazy_eval}
if (requireNamespace("dplyr", quietly = TRUE)) {
    # collect = FALSE for lazy evaluation
    lazy_tbl <- ol_read("raw_counts", collect = FALSE)

    # Process with dplyr
    result <- lazy_tbl %>%
        dplyr::filter(sample1 > 100) %>%
        dplyr::select(gene_id, sample1) %>%
        dplyr::collect()

    head(result)
}
```

## 7.3 Deleting tables

```{r drop}
# Remove an unnecessary table
ol_drop("filtered_de")
```

# 8. Bioconductor integration

## 8.1 Creating a SummarizedExperiment

```{r bioc_se}
if (requireNamespace("SummarizedExperiment", quietly = TRUE)) {
    # Create long-format data
    long_counts <- data.frame(
        feature = rep(paste0("GENE", 1:100), each = 4),
        sample = rep(paste0("sample", 1:4), times = 100),
        value = rpois(400, 100)
    )
    ol_write("long_counts", long_counts, mode = "overwrite")

    # Load as a SummarizedExperiment
    se <- ol_read_se("long_counts",
        feature_col = "feature",
        sample_col = "sample",
        value_col = "value",
        backing = "memory"
    )

    print(se)
}
```

## 8.2 Creating a MultiAssayExperiment

```{r bioc_mae}
if (requireNamespace("MultiAssayExperiment", quietly = TRUE)) {
    # Prepare multiple assay datasets
    rna_data <- data.frame(
        feature = rep(paste0("GENE", 1:50), each = 4),
        sample = rep(paste0("sample", 1:4), times = 50),
        value = rpois(200, 100)
    )
    ol_write("rna_assay", rna_data, mode = "overwrite")

    protein_data <- data.frame(
        feature = rep(paste0("PROT", 1:30), each = 4),
        sample = rep(paste0("sample", 1:4), times = 30),
        value = rnorm(120, 50, 10)
    )
    ol_write("protein_assay", protein_data, mode = "overwrite")

    # Load as a MultiAssayExperiment
    mae <- ol_read_mae(
        assays = list(
            rna = list(name = "rna_assay"),
            protein = list(name = "protein_assay")
        ),
        backing = "memory"
    )

    print(mae)
}
```

# 9. Practical workflow example

## 9.1 Complete differential expression analysis workflow

```{r complete_workflow}
# 1. Initialize the project
ol_init("rnaseq_project")

# 2. Save raw data
set.seed(456)
raw_data <- data.frame(
    gene_id = paste0("GENE", 1:200),
    control_1 = rpois(200, 100),
    control_2 = rpois(200, 100),
    control_3 = rpois(200, 100),
    treated_1 = rpois(200, 120),
    treated_2 = rpois(200, 120),
    treated_3 = rpois(200, 120)
)
ol_write("raw_expression", raw_data)

# 3. Save QC parameters
qc_params <- list(
    min_counts = 10,
    min_samples = 3,
    method = "standard"
)
ol_save("qc_parameters", qc_params)

# 4. Post-QC data
qc_data <- raw_data[rowSums(raw_data[, -1] > 10) >= 3, ]
ol_write(
    "qc_filtered",
    qc_data,
    depends_on = c("raw_expression", "qc_parameters")
)

# 5. Normalization parameters
norm_params <- list(
    method = "TMM",
    log_transform = TRUE
)
ol_save("norm_parameters", norm_params)

# 6. Normalized data
norm_data <- qc_data
for (i in 2:ncol(norm_data)) {
    norm_data[[i]] <- log2(norm_data[[i]] + 1)
}
ol_write("normalized", norm_data,
    depends_on = c("qc_filtered", "norm_parameters")
)

# 7. Statistical analysis parameters
stats_params <- list(
    test = "t-test",
    alpha = 0.05,
    correction = "BH"
)
ol_save("stats_parameters", stats_params)

# 8. Statistical analysis results
set.seed(789)
de_final <- data.frame(
    gene_id = sample(qc_data$gene_id, 30),
    log2fc = rnorm(30, 1, 0.5),
    pvalue = runif(30, 0, 0.05),
    padj = runif(30, 0, 0.05)
)
ol_save("final_de_results", de_final,
    depends_on = c("normalized", "stats_parameters")
)

# 9. Commit the analysis
ol_commit(
    note = "Complete RNA-seq differential expression analysis",
    params = list(
        samples = 6,
        genes_tested = nrow(qc_data),
        significant = nrow(de_final),
        date = as.character(Sys.Date())
    )
)

# 10. Create a label
ol_label("final_analysis_v1")

# 11. View the full lineage
full_lineage <- ol_show_lineage("final_de_results", direction = "upstream")
print(full_lineage)

# 12. Visualize dependencies
if (requireNamespace("igraph", quietly = TRUE)) {
    ol_plot_lineage("final_de_results",
        direction = "upstream",
        layout = "sugiyama",
        main = "RNA-seq Analysis Pipeline"
    )
}
```

# 10. Conclusion

## List of key functions

### Project management
- `ol_init()`: initialize a project
- `ol_label()`: label the entire project
- `ol_checkout()`: restore to a labeled state

### Data saving and loading
- `ol_write()`: save a table
- `ol_read()`: read a table
- `ol_save()`: save an R object
- `ol_read_object()`: load an R object
- `ol_load()`: alias for `ol_read()`
- `ol_fread()`: read with filtering

### Version management
- `ol_tag()`: tag a table
- `ol_tag_object()`: tag an object
- `ol_list_object_versions()`: list versions
- `ol_compare_versions()`: compare versions

### History management
- `ol_commit()`: create a commit
- `ol_log()`: show table history
- `ol_log_commits()`: show commit history

### Dependencies
- `ol_get_dependencies()`: get dependencies
- `ol_show_lineage()`: display the full lineage tree
- `ol_plot_lineage()`: visualize dependency graphs

### Listing
- `ol_list_tables()`: list tables
- `ol_list_objects()`: list objects
- `ol_list_tags()`: list tags
- `ol_list_labels()`: list labels

### Bioconductor
- `ol_read_se()`: create a *SummarizedExperiment*
- `ol_read_mae()`: create a *MultiAssayExperiment*

### Others
- `ol_drop()`: delete a table

## Best practices

1. **Explicitly specify dependencies**: use the `depends_on` parameter to
1. record the analysis flow.
2. **Create labels at important milestones**: run `ol_label()` at key points to
2. ensure reproducibility.
3. **Leverage tags**: assign clear tags for different analysis methods or
3. parameter sets.
4. **Write detailed commit messages**: record the analysis content and results
4. in `ol_commit()`.
5. **Confirm with visualization**: visualize complex analysis flows with
5. `ol_plot_lineage()`.
6. **Use version comparison**: track the impact of parameter changes with
6. `ol_compare_versions()`.

## Session Information

```{r session_info}
sessionInfo()
```
