## ----setup, include=FALSE-----------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment  = "#>",
  message  = FALSE,
  warning  = FALSE
)

## ----load---------------------------------------------------------------------
library(CONCERTDR)

## ----load-data----------------------------------------------------------------
sig_file  <- system.file("extdata", "example_signature.txt",
                          package = "CONCERTDR")
ref_file  <- system.file("extdata", "example_reference_df.csv",
                          package = "CONCERTDR")
sig_info_file  <- system.file("extdata", "example_siginfo.txt",
                               package = "CONCERTDR")
comp_info_file <- system.file("extdata", "example_compoundinfo.txt",
                               package = "CONCERTDR")

# Load the reference matrix
reference_df <- read.csv(ref_file, row.names = 1, check.names = FALSE)
reference_df$gene_symbol <- rownames(reference_df)

cat("Reference matrix:", nrow(reference_df), "genes ×",
    ncol(reference_df) - 1, "signatures\n")

## ----view-signature-----------------------------------------------------------
signature <- read.delim(sig_file)
cat("Up-regulated genes:  ", sum(signature$log2FC > 0), "\n")
cat("Down-regulated genes:", sum(signature$log2FC < 0), "\n")
head(signature[order(-abs(signature$log2FC)), ], 10)

## ----view-siginfo-------------------------------------------------------------
siginfo <- read.delim(sig_info_file, stringsAsFactors = FALSE)
cat("Total signatures:", nrow(siginfo), "\n")
table(Type = siginfo$pert_type, Cell = siginfo$cell_iname)

## ----filter-compounds---------------------------------------------------------
cp_siginfo <- subset_siginfo_beta(
  siginfo_file = sig_info_file,
  interactive  = FALSE,
  filters      = list(pert_type = "trt_cp"),
  verbose      = FALSE,
  show_preview = FALSE
)
cat("Compound signatures:", nrow(cp_siginfo), "\n")
table(cp_siginfo$cell_iname)

## ----subset-reference-cp------------------------------------------------------
cp_sig_ids <- cp_siginfo$sig_id
cp_ref     <- reference_df[, c("gene_symbol",
                                intersect(cp_sig_ids, colnames(reference_df))),
                            drop = FALSE]
cat("Compound reference:", nrow(cp_ref), "genes ×",
    ncol(cp_ref) - 1, "signatures\n")

## ----score-compounds----------------------------------------------------------
results_cp <- process_signature_with_df(
  signature_file = sig_file,
  reference_df   = cp_ref,
  methods        = c("xsum", "ks"),
  topN           = 30,
  permutations   = 100,
  save_files     = FALSE
)
print(results_cp)

## ----top-compounds------------------------------------------------------------
xsum_cp <- results_cp$results$xsum
top_cp  <- head(xsum_cp[order(xsum_cp$Score), ], 20)
top_cp[, c("compound", "Score", "pValue", "pAdjValue")]

## ----plot-compound-scores, fig.width=8, fig.height=4--------------------------
plot(results_cp, method = "xsum", plot_type = "scores", top_n = 20)

## ----plot-compound-volcano, fig.width=6, fig.height=4-------------------------
plot(results_cp, method = "xsum", plot_type = "volcano")

## ----annotate-compounds-------------------------------------------------------
compinfo <- read.delim(comp_info_file, stringsAsFactors = FALSE)

views_cp <- annotate_drug_results(
  results_df     = xsum_cp,
  sig_info_file  = cp_siginfo,
  comp_info_file = compinfo,
  write_outputs  = FALSE,
  verbose        = FALSE
)

# Drug-level summary: one row per compound, best score across all contexts
head(views_cp$wetlab_drug_view[,
  c("perturbation_name", "Score", "effect_direction", "moa", "target")
], 15)

## ----drug-context-summary-----------------------------------------------------
# How many cell-line contexts did each top drug appear in?
head(views_cp$drug_context_summary[,
  c("perturbation_name", "best_score", "n_contexts", "n_cell_lines", "moa_status")
], 10)

## ----extract-zscores-cp-------------------------------------------------------
z_cp <- extract_signature_zscores(
  results_df      = views_cp$tech_view_all,
  signature_file  = sig_file,
  reference_df    = cp_ref,
  max_genes       = 50,           # top-50 up-regulated + top-50 down-regulated
  split_direction = TRUE,
  max_perts       = 50,
  verbose         = FALSE
)
cat("Heatmap matrix:", nrow(z_cp$z_plot), "perturbations ×",
    ncol(z_cp$z_plot), "genes\n")

## ----heatmap-compounds, fig.width=14, fig.height=10---------------------------
if (requireNamespace("ComplexHeatmap", quietly = TRUE) &&
    requireNamespace("circlize",       quietly = TRUE)) {
  plot_signature_direction_tile_barcode(
    precomputed         = z_cp,
    cluster_rows        = TRUE,
    show_row_dendrogram = TRUE,
    cluster_cols        = FALSE,
    verbose             = FALSE
  )
} else {
  message("Install ComplexHeatmap and circlize for heatmap rendering:\n",
          "  BiocManager::install('ComplexHeatmap')\n",
          "  install.packages('circlize')")
}

## ----heatmap-compounds-split, fig.width=16, fig.height=10---------------------
if (requireNamespace("ComplexHeatmap", quietly = TRUE) &&
    requireNamespace("circlize",       quietly = TRUE)) {
  plot_signature_direction_tile_barcode(
    precomputed         = z_cp,
    split_direction     = TRUE,   # up-regulated panel (left) | down-regulated panel (right)
    gap_width           = 6,      # mm gap between the two panels
    cluster_rows        = TRUE,
    show_row_dendrogram = TRUE,
    cluster_cols        = FALSE,
    verbose             = FALSE
  )
}

## ----filter-gene-perts--------------------------------------------------------
ge_siginfo <- subset_siginfo_beta(
  siginfo_file = sig_info_file,
  interactive  = FALSE,
  filters      = list(pert_type = c("trt_oe", "trt_xpr")),
  verbose      = FALSE,
  show_preview = FALSE
)
cat("Gene perturbation signatures:", nrow(ge_siginfo), "\n")
table(Type = ge_siginfo$pert_type, Cell = ge_siginfo$cell_iname)

## ----subset-reference-ge------------------------------------------------------
ge_sig_ids <- ge_siginfo$sig_id
ge_ref     <- reference_df[, c("gene_symbol",
                                intersect(ge_sig_ids, colnames(reference_df))),
                            drop = FALSE]
cat("Gene perturbation reference:", nrow(ge_ref), "genes ×",
    ncol(ge_ref) - 1, "signatures\n")

## ----score-gene-perts---------------------------------------------------------
results_ge <- process_signature_with_df(
  signature_file = sig_file,
  reference_df   = ge_ref,
  methods        = c("xsum", "ks"),
  topN           = 50,
  permutations   = 100,
  save_files     = FALSE
)
print(results_ge)

## ----top-gene-hits------------------------------------------------------------
xsum_ge <- results_ge$results$xsum
top_ge  <- head(xsum_ge[order(xsum_ge$Score), ], 20)
top_ge[, c("compound", "Score", "pValue", "pAdjValue")]

## ----plot-gene-scores, fig.width=8, fig.height=4------------------------------
plot(results_ge, method = "xsum", plot_type = "scores", top_n = 20)

## ----annotate-gene-perts------------------------------------------------------
views_ge <- annotate_drug_results(
  results_df     = xsum_ge,
  sig_info_file  = ge_siginfo,
  comp_info_file = compinfo,   # compinfo may not cover gene perts; that is expected
  write_outputs  = FALSE,
  verbose        = FALSE
)

# Gene-focused view: perturbation name, mode (OE vs XPR), score
head(views_ge$wetlab_gene_view[,
  c("perturbation_name", "mode", "Score", "effect_direction", "cell_line")
], 15)

## ----extract-zscores-ge-------------------------------------------------------
z_ge <- extract_signature_zscores(
  results_df      = views_ge$tech_view_all,
  signature_file  = sig_file,
  reference_df    = ge_ref,
  max_genes       = 50,
  split_direction = TRUE,
  max_perts       = 50,
  verbose         = FALSE
)

## ----heatmap-gene-perts, fig.width=14, fig.height=10--------------------------
if (requireNamespace("ComplexHeatmap", quietly = TRUE) &&
    requireNamespace("circlize",       quietly = TRUE)) {
  plot_signature_direction_tile_barcode(
    precomputed         = z_ge,
    cluster_rows        = TRUE,
    show_row_dendrogram = TRUE,
    cluster_cols        = FALSE,
    verbose             = FALSE
  )
}

## ----heatmap-gene-perts-split, fig.width=16, fig.height=10--------------------
if (requireNamespace("ComplexHeatmap", quietly = TRUE) &&
    requireNamespace("circlize",       quietly = TRUE)) {
  plot_signature_direction_tile_barcode(
    precomputed     = z_ge,
    split_direction = TRUE,
    gap_width       = 6,
    cluster_rows    = TRUE,
    cluster_cols    = FALSE,
    verbose         = FALSE
  )
}

## ----compare-use-cases--------------------------------------------------------
# Top 10 compounds vs top 10 gene perturbations by XSum score
top10_cp <- head(xsum_cp[order(xsum_cp$Score), "Score", drop = FALSE], 10)
top10_ge <- head(xsum_ge[order(xsum_ge$Score), "Score", drop = FALSE], 10)

cat("Score range — compounds:         ",
    round(range(top10_cp$Score), 3), "\n")
cat("Score range — gene perturbations:",
    round(range(top10_ge$Score), 3), "\n")

## ----full-cmap-setup, eval=FALSE----------------------------------------------
# data_dir <- "/path/to/cmap_data"
# 
# options(
#   CONCERTDR.data_dir          = data_dir,
#   CONCERTDR.gctx_file         = file.path(data_dir, "level5_beta_trt_cp_n720216x12328.gctx"),
#   CONCERTDR.siginfo_file      = file.path(data_dir, "siginfo_beta.txt"),
#   CONCERTDR.geneinfo_file     = file.path(data_dir, "geneinfo_beta.txt"),
#   CONCERTDR.compoundinfo_file = file.path(data_dir, "compoundinfo_beta.txt")
# )
# 
# # Step 1: filter siginfo
# cp_siginfo <- subset_siginfo_beta(
#   getOption("CONCERTDR.siginfo_file"),
#   interactive = FALSE,
#   filters = list(
#     pert_type  = "trt_cp",
#     pert_itime = c("6 h", "24 h"),
#     cell_iname = c("K562", "HL60", "THP1", "JURKAT", "U937", "NB4")
#   )
# )
# 
# # Step 2: build reference matrix from GCTX
# reference_df <- extract_cmap_data_from_siginfo(
#   siginfo_file  = cp_siginfo,
#   geneinfo_file = getOption("CONCERTDR.geneinfo_file"),
#   gctx_file     = getOption("CONCERTDR.gctx_file"),
#   filter_quality = FALSE,
#   landmark       = TRUE
# )
# 
# # Steps 3–7: identical to the vignette above — just replace cp_ref with
# # reference_df and use larger topN / permutations for real-scale analysis:
# #   topN = 400, permutations = 1000

## ----session-info-------------------------------------------------------------
sessionInfo()