options(omicslake.vignette.eval = TRUE)

bash tools/run_demo_count_edger_limma_voom_ora.sh

library(OmicsLake)

# プロジェクトを初期化
ol_init("de_analysis")

# サンプルデータの作成
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)
)

# テーブルとして保存
ol_write("raw_counts", raw_counts, mode = "create")

# テーブルの読み込み
counts <- ol_read("raw_counts")
head(counts)

# テーブル一覧の確認
ol_list_tables()

# 基本的なクエリ
result <- ol_query("SELECT * FROM raw_counts WHERE sample1 > 100")
head(result)

# 集計クエリ
summary <- ol_query("
    SELECT
    COUNT(*) as total_genes,
    AVG(sample1) as avg_sample1,
    MAX(sample1) as max_sample1
    FROM raw_counts
")
print(summary)

# 正規化係数を計算（例）
norm_factors <- list(
    method = "TMM",
    factors = runif(100, 0.8, 1.2)
)

# Rオブジェクトとして保存
ol_save("norm_factors", norm_factors)

# オブジェクトの読み込み
loaded_factors <- ol_read_object("norm_factors")
str(loaded_factors)

# オブジェクト一覧の確認
ol_list_objects()

# 正規化されたカウントを計算（depends_on で依存関係を指定）
normalized_counts <- raw_counts
for (i in 2:5) {
    normalized_counts[[i]] <- normalized_counts[[i]] * norm_factors$factors
}

# 依存関係を指定して保存
ol_write("normalized_counts", normalized_counts,
    depends_on = c("raw_counts", "norm_factors")
)

# DEパラメータを保存
de_params <- list(
    method = "DESeq2",
    alpha = 0.05,
    lfc_threshold = 1.0
)
ol_save("de_params", de_params)

# DE結果を計算（簡略化した例）
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)
)

# DE結果を依存関係付きで保存
ol_save("de_results", de_results,
    depends_on = c("normalized_counts", "de_params")
)

# 上流の依存関係（このオブジェクトが依存しているもの）
upstream <- ol_get_dependencies("de_results", direction = "upstream")
print(upstream)

# 下流の依存関係（このオブジェクトに依存しているもの）
downstream <- ol_get_dependencies("raw_counts", direction = "downstream")
print(downstream)

# 完全な系統樹を表示
lineage <- ol_show_lineage("de_results", direction = "upstream")
print(lineage)

# オブジェクトにタグを付ける
ol_tag_object("de_results", "baseline_analysis")

# プロジェクト全体にラベルを付ける（全テーブル・オブジェクトにタグが付く）
ol_label("experiment_v1")

# タグ一覧の確認
ol_list_tags()

# プロジェクトラベル一覧
ol_list_labels()

# パラメータを変更して再解析
de_params_v2 <- list(
    method = "DESeq2",
    alpha = 0.01,
    lfc_threshold = 1.5
)
ol_save("de_params", de_params_v2)

# 新しいパラメータでDE結果を再計算
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")
)

# 新バージョンにタグ付け
ol_tag_object("de_results", "strict_analysis")

# 別の解析手法を試す
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")

# 全バージョンをリスト表示
versions <- ol_list_object_versions("de_results")
print(versions)

# バージョン間の比較
comparison <- ol_compare_versions("de_results")
print(comparison)

# 特定のタグ間の比較
tag_comparison <- ol_compare_versions("de_results",
    versions = c("baseline_analysis", "strict_analysis")
)
print(tag_comparison)

# タグを指定して読み込み
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")

# 解析の節目でコミットを作成
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")

# コミット履歴を表示
commits <- ol_log_commits(n = 10)
print(commits)

# 特定のテーブルの履歴
table_log <- ol_log("raw_counts")
print(table_log)

# igraphパッケージが必要
if (requireNamespace("igraph", quietly = TRUE)) {
    # 上流の依存関係を可視化
    ol_plot_lineage("de_results",
        direction = "upstream",
        layout = "sugiyama",
        main = "DE Results - Upstream Dependencies"
    )

    # 双方向の依存関係を可視化
    ol_plot_lineage("normalized_counts",
        direction = "both",
        layout = "tree",
        main = "Normalized Counts - Full Lineage"
    )
}

# 新しいテーブルを追加
filtered_results <- de_results[de_results$padj < 0.05, ]
ol_write("filtered_de", filtered_results)

# 以前のラベル状態に戻る
ol_checkout("experiment_v1")

# filtered_deテーブルは存在しなくなる
current_tables <- ol_list_tables()
cat(
    "Tables after checkout:",
    paste(current_tables$table_name, collapse = ", "),
    "\n"
)

# 再度raw_countsを作成（checkoutで消えた可能性があるため）
ol_write("raw_counts", raw_counts, mode = "overwrite")

# 特定の列のみを選択して読み込み
selected <- ol_fread("raw_counts",
    select = c("gene_id", "sample1", "sample2"),
    nrows = 10
)
head(selected)

# 条件でフィルタリング
filtered <- ol_fread("raw_counts",
    filter = "sample1 > 100"
)
head(filtered)

if (requireNamespace("dplyr", quietly = TRUE)) {
    # collect = FALSE で遅延評価
    lazy_tbl <- ol_read("raw_counts", collect = FALSE)

    # dplyrで処理
    result <- lazy_tbl %>%
        dplyr::filter(sample1 > 100) %>%
        dplyr::select(gene_id, sample1) %>%
        dplyr::collect()

    head(result)
}

# 不要なテーブルを削除
ol_drop("filtered_de")

if (requireNamespace("SummarizedExperiment", quietly = TRUE)) {
    # ロングフォーマットのデータを作成
    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")

    # SummarizedExperimentとして読み込み
    se <- ol_read_se("long_counts",
        feature_col = "feature",
        sample_col = "sample",
        value_col = "value",
        backing = "memory"
    )

    print(se)
}

if (requireNamespace("MultiAssayExperiment", quietly = TRUE)) {
    # 複数のアッセイデータを準備
    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")

    # MultiAssayExperimentとして読み込み
    mae <- ol_read_mae(
        assays = list(
            rna = list(name = "rna_assay"),
            protein = list(name = "protein_assay")
        ),
        backing = "memory"
    )

    print(mae)
}

# 1. プロジェクト初期化
ol_init("rnaseq_project")

# 2. 生データの保存
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. QCパラメータの保存
qc_params <- list(
    min_counts = 10,
    min_samples = 3,
    method = "standard"
)
ol_save("qc_parameters", qc_params)

# 4. QC後のデータ
qc_data <- raw_data[rowSums(raw_data[, -1] > 10) >= 3, ]
ol_write(
    "qc_filtered",
    qc_data,
    depends_on = c("raw_expression", "qc_parameters")
)

# 5. 正規化パラメータ
norm_params <- list(
    method = "TMM",
    log_transform = TRUE
)
ol_save("norm_parameters", norm_params)

# 6. 正規化データ
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. 統計解析パラメータ
stats_params <- list(
    test = "t-test",
    alpha = 0.05,
    correction = "BH"
)
ol_save("stats_parameters", stats_params)

# 8. 統計解析結果
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. 解析をコミット
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. ラベルを作成
ol_label("final_analysis_v1")

# 11. 完全な系統を確認
full_lineage <- ol_show_lineage("final_de_results", direction = "upstream")
print(full_lineage)

# 12. 依存関係を可視化
if (requireNamespace("igraph", quietly = TRUE)) {
    ol_plot_lineage("final_de_results",
        direction = "upstream",
        layout = "sugiyama",
        main = "RNA-seq Analysis Pipeline"
    )
}

sessionInfo()