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

ensure_pkgs_available <- function(pkgs) {
  missing <- pkgs[!vapply(pkgs, requireNamespace, logical(1), quietly = TRUE)]
  if (length(missing)) {
    stop("Missing required packages: ", paste(missing, collapse = ", "))
  }
}

# This vignette requires all dependencies to already be installed.
ensure_pkgs_available(c("VISTA", "ggplot2", "airway", "org.Hs.eg.db"))

## ----install-packages, eval=FALSE---------------------------------------------
# if (!requireNamespace("BiocManager", quietly = TRUE)) {
#   install.packages("BiocManager")
# }
# 
# BiocManager::install(c("VISTA", "airway", "org.Hs.eg.db"))
# install.packages("ggplot2")

## ----load-packages------------------------------------------------------------
# Load required packages
library(VISTA)
library(ggplot2)         # For plotting functions
library(airway)          # Dataset
library(org.Hs.eg.db)    # Human gene annotations
library(magrittr)    # %>% 

## ----load-data----------------------------------------------------------------
# Load the SummarizedExperiment object
data("airway", package = "airway")

# Examine the structure
airway

## ----extract-data-------------------------------------------------------------
# Extract count matrix
counts_matrix <- assay(airway, "counts")

# Preview counts (first 5 genes, first 4 samples)
counts_matrix[1:5, 1:4]

# Extract sample metadata
sample_metadata <- as.data.frame(colData(airway))
sample_metadata

# Simplify treatment labels for clarity
sample_metadata$treatment <- ifelse(
  sample_metadata$dex == "trt",
  "Dexamethasone",
  "Untreated"
)

## ----prepare-inputs-----------------------------------------------------------
prepared_counts <- read_vista_counts(
  counts_matrix,
  format = "matrix"
)

prepared_samples <- read_vista_metadata(
  sample_metadata %>%
    tibble::as_tibble() %>%
    dplyr::rename("sample_names" = "Run")
)

matched_inputs <- match_vista_inputs(prepared_counts, prepared_samples)

# Preview the create_vista-ready count table
matched_inputs$counts[1:5, 1:5]

## ----prepare-samples----------------------------------------------------------
matched_inputs$sample_info[, c("sample_names", "cell", "treatment", "dex")]

## ----create-vista-------------------------------------------------------------
# Create VISTA object with DESeq2 backend
vista <- create_vista(
  counts = matched_inputs$counts,
  sample_info = matched_inputs$sample_info,
  column_geneid = matched_inputs$column_geneid,
  group_column = "treatment",
  group_numerator = "Dexamethasone",
  group_denominator = "Untreated",
  method = "deseq2",
  min_counts = 10,
  min_replicates = 2,
  log2fc_cutoff = 1.0,
  pval_cutoff = 0.05,
  p_value_type = "padj"
)

# Examine the VISTA object
vista

## ----validate-vista-----------------------------------------------------------
validate_vista(vista, level = "full")

## ----as-vista-example, eval=FALSE---------------------------------------------
# se <- SummarizedExperiment::SummarizedExperiment(
#   assays = list(norm_counts = norm_counts(vista)),
#   colData = S4Vectors::DataFrame(sample_info(vista), row.names = sample_info(vista)$sample_names),
#   rowData = S4Vectors::DataFrame(row_data(vista), row.names = rownames(norm_counts(vista)))
# )
# vista2 <- as_vista(se, group_column = "treatment")
# validate_vista(vista2, level = "full")

## ----create-vista-edger, eval=FALSE-------------------------------------------
# # Create VISTA object with edgeR backend
# vista_edger <- create_vista(
#   counts = count_data,
#   sample_info = sample_info,
#   column_geneid = "gene_id",
#   group_column = "treatment",
#   group_numerator = "Dexamethasone",
#   group_denominator = "Untreated",
#   method = "edger",  # Use edgeR instead of DESeq2
#   min_counts = 10,
#   min_replicates = 2,
#   log2fc_cutoff = 1.0,
#   pval_cutoff = 0.05,
#   p_value_type = "padj"
# )

## ----create-vista-limma, eval=FALSE-------------------------------------------
# # Create VISTA object with limma-voom backend
# vista_limma <- create_vista(
#   counts = count_data,
#   sample_info = sample_info,
#   column_geneid = "gene_id",
#   group_column = "treatment",
#   group_numerator = "Dexamethasone",
#   group_denominator = "Untreated",
#   method = "limma",
#   min_counts = 10,
#   min_replicates = 2,
#   log2fc_cutoff = 1.0,
#   pval_cutoff = 0.05,
#   p_value_type = "padj"
# )

## ----advanced-design-consensus, eval=FALSE------------------------------------
# # Covariate-adjusted model (additive design)
# vista_cov <- create_vista(
#   counts = count_data,
#   sample_info = sample_info,
#   column_geneid = "gene_id",
#   group_column = "treatment",
#   group_numerator = "Dexamethasone",
#   group_denominator = "Untreated",
#   method = "deseq2",
#   covariates = c("cell")
# )
# 
# # Equivalent explicit model formula
# vista_formula <- create_vista(
#   counts = count_data,
#   sample_info = sample_info,
#   column_geneid = "gene_id",
#   group_column = "treatment",
#   group_numerator = "Dexamethasone",
#   group_denominator = "Untreated",
#   method = "deseq2",
#   design_formula = "~ cell + treatment"
# )
# 
# # Run both DESeq2 and edgeR and keep consensus as active source
# vista_both <- create_vista(
#   counts = count_data,
#   sample_info = sample_info,
#   column_geneid = "gene_id",
#   group_column = "treatment",
#   group_numerator = "Dexamethasone",
#   group_denominator = "Untreated",
#   method = "both",
#   consensus_mode = "intersection",  # or "union"
#   result_source = "consensus"       # or "deseq2"/"edger"
# )
# 
# # Access source-specific outputs
# comparisons(vista_both, source = "consensus")
# comparisons(vista_both, source = "deseq2")
# comparisons(vista_both, source = "edger")
# 
# # Switch the active source used by plotting functions
# vista_both <- set_de_source(vista_both, "edger")

## ----add-annotations----------------------------------------------------------
vista <- set_rowdata(
  vista,
  orgdb = org.Hs.eg.db,
  columns = c("SYMBOL", "GENENAME", "ENTREZID"),
  keytype = "ENSEMBL"
)

# View updated gene annotations
head(rowData(vista))

## ----access-results-----------------------------------------------------------
# Get comparison names
comp_names <- names(comparisons(vista))
comp_names

# Get DE results for the comparison
de_results <- comparisons(vista)[[1]]
head(de_results)

# Get DEG summary
deg_summary(vista)

# Get analysis cutoffs
cutoffs(vista)

## ----count-degs---------------------------------------------------------------
# Extract upregulated genes
up_genes <- get_genes_by_regulation(
  vista,
  sample_comparisons = comp_names[1],
  regulation = "Up",
  #top_n = 50
)

# Extract downregulated genes
down_genes <- get_genes_by_regulation(
  vista,
  sample_comparisons = comp_names[1],
  regulation = "Down",
  #top_n = 50
)

# Summary
cat("Upregulated genes:", length(up_genes[[1]]), "\n")
cat("Downregulated genes:", length(down_genes[[1]]), "\n")

## ----corr-heatmap-basic, fig.width=7, fig.height=6----------------------------
  get_corr_heatmap(vista, label_size = 18,base_size = 18, viridis_direction = -1)

## ----corr-heatmap-colors, fig.width=7, fig.height=6---------------------------
# Reverse viridis color direction
get_corr_heatmap(
  vista,
  viridis_direction = -1,
  viridis_option = "plasma",
  label_size = 18,
  base_size = 18
)

## ----corr-heatmap-values, fig.width=7, fig.height=6---------------------------
# Display correlation coefficients
get_corr_heatmap(
  vista,
  viridis_direction = -1,
  show_corr_values = TRUE,
  col_corr_values = 'white',
  viridis_option = "mako",
  label_size = 14,
  base_size = 14
) 

## ----pca-basic, fig.width=8, fig.height=6-------------------------------------
get_pca_plot(
  vista,
  label = TRUE,label_size = 5
)

## ----pca-shape, fig.width=8, fig.height=6-------------------------------------
# Shape points by cell line
get_pca_plot(
  vista,
  label = TRUE,
  label_size = 5,
  shape_by = "cell"
)

## ----pca-top-genes, fig.width=8, fig.height=6---------------------------------
# Use top 500 most variable genes
get_pca_plot(
  vista,
  top_n_genes = 500,
  show_clusters = TRUE,
  shape_by = "cell"
)

## ----pca-circle-size, fig.width=8, fig.height=6-------------------------------
# Larger points for better visibility
get_pca_plot(
  vista,
  label = TRUE,
  point_size = 15,label_size = 5
)

## ----pca-no-labels, fig.width=8, fig.height=6---------------------------------
# Clean plot without sample labels
get_pca_plot(
  vista,
  label = FALSE,
  point_size = 12
)

## ----mds-basic, fig.width=8, fig.height=6-------------------------------------
get_mds_plot(
  vista,
  label = TRUE
)

## ----mds-top-genes, fig.width=8, fig.height=6---------------------------------
get_mds_plot(
  vista,
  top_n_genes = 500,
  label = TRUE
)

## ----mds-shapes, fig.width=8, fig.height=6------------------------------------
# Shape points by cell line
get_mds_plot(
  vista,
  label = TRUE,
  shape_by = "cell"
)

## ----umap-basic, fig.width=8, fig.height=6, eval=requireNamespace("uwot", quietly=TRUE)----
get_umap_plot(
  vista,
  label = TRUE
)

## ----umap-color-by-cell, fig.width=8, fig.height=6, eval=requireNamespace("uwot", quietly=TRUE)----
get_umap_plot(
  vista,
  color_by = "cell",
  shape_by = "treatment",
  label = TRUE
)

## ----deg-count-basic, fig.width=7, fig.height=5-------------------------------
get_deg_count_barplot(vista)

## ----deg-count-facet, fig.width=8, fig.height=5-------------------------------
get_deg_count_barplot(
  vista,
  facet_by = "regulation"
)

## ----volcano-basic, fig.width=8, fig.height=7, eval=requireNamespace("EnhancedVolcano", quietly=TRUE)----
get_volcano_plot(
  vista,
  sample_comparison = comp_names[1]
)

## ----volcano-custom, fig.width=8, fig.height=7, eval=requireNamespace("EnhancedVolcano", quietly=TRUE)----
get_volcano_plot(
  vista,
  sample_comparison = comp_names[1],
  log2fc_cutoff = 1.5,
  pval_cutoff = 0.01,
  label_size = 3,
  display_id = "SYMBOL"
)

## ----volcano-colors, fig.width=8, fig.height=7, eval=requireNamespace("EnhancedVolcano", quietly=TRUE)----
# Custom colors for up/down regulated genes
get_volcano_plot(
  vista,
  sample_comparison = comp_names[1],
  col_up = "red",
  col_down = "blue",
  display_id = "SYMBOL"
)

## ----ma-basic, fig.width=8, fig.height=6--------------------------------------
get_ma_plot(
  vista,
  sample_comparison = comp_names[1]
)

## ----ma-labeled, fig.width=8, fig.height=6------------------------------------
get_ma_plot(
  vista,
  sample_comparison = comp_names[1],
  label_n = 10,
  point_size = 2,
  display_id = "SYMBOL"
)

## ----ma-custom, fig.width=8, fig.height=6-------------------------------------
get_ma_plot(
  vista,
  sample_comparison = comp_names[1],
  label_n = 5,
  display_id = "SYMBOL",
  point_size = 1.5,
  alpha = 0.8
)

## ----prep-gene-sets-----------------------------------------------------------
# Get top 50 DEGs by adjusted p-value
de_table <- comparisons(vista)[[1]]

top_degs <- get_genes_by_regulation(vista, names(comparisons(vista))[1],regulation = "Both",top_n = 50)[[1]]  # top 50 by abs fold change 

top_up <- get_genes_by_regulation(vista, names(comparisons(vista))[1],regulation = "Up",top_n = 50)[[1]]  

top_down <- get_genes_by_regulation(vista, names(comparisons(vista))[1],regulation = "Down",top_n = 50)[[1]]  


cat("Top upregulated genes:\n")
print(top_up)
cat("\nTop downregulated genes:\n")
print(top_down)

## ----heatmap-basic, fig.width=8, fig.height=9, eval=requireNamespace("ComplexHeatmap", quietly=TRUE)----
get_expression_heatmap(vista)

## ----heatmap-explicit, fig.width=8, fig.height=9, eval=requireNamespace("ComplexHeatmap", quietly=TRUE)----
get_expression_heatmap(
  vista,
  sample_group = levels(colData(vista)$treatment),
  genes = top_degs,
  display_id = "SYMBOL",
  show_row_names = TRUE
)

## ----heatmap-kmeans, fig.width=8, fig.height=9, eval=requireNamespace("ComplexHeatmap", quietly=TRUE)----
get_expression_heatmap(
  vista,
  sample_group = levels(colData(vista)$treatment),
  genes = top_degs,
  display_id = "SYMBOL",
  kmeans_k = 3,
  show_row_names = TRUE
)

## ----heatmap-annotated, fig.width=8, fig.height=9, eval=requireNamespace("ComplexHeatmap", quietly=TRUE)----
get_expression_heatmap(
  vista,
  sample_group = levels(colData(vista)$treatment),
  genes = top_degs,
  show_row_names = TRUE,
  display_id = "SYMBOL",
  kmeans_k = 3,
  cluster_row_slice = FALSE,
  summarise_replicates = FALSE,
  annotate_columns = TRUE
)

## ----heatmap-annotated-multi, fig.width=8, fig.height=9, eval=requireNamespace("ComplexHeatmap", quietly=TRUE)----
# Use multiple sample-level columns in top annotation.
# By default, columns are split by the first annotation column.


get_expression_heatmap(
  vista,
  sample_group = levels(colData(vista)$treatment),
  genes = top_degs,
  show_row_names = FALSE,
  display_id = "SYMBOL",
  summarise_replicates = FALSE,
  annotate_columns = c("treatment", "cell"),
  cluster_by = "cell"
)

## ----heatmap-summarized, fig.width=7, fig.height=9, eval=requireNamespace("ComplexHeatmap", quietly=TRUE)----
get_expression_heatmap(
  vista,
  sample_group = levels(colData(vista)$treatment),
  genes = top_degs,
  display_id = "SYMBOL",
  summarise_replicates = FALSE,
  show_row_names = FALSE,
  annotate_columns = TRUE,
  kmeans_k = 2
)

## ----barplot-basic, fig.width=9, fig.height=6---------------------------------
get_expression_barplot(
  vista,
  genes = top_up[1:4],
  display_id = "SYMBOL",
)+theme_minimal(base_size = 16)

## ----barplot-log-stats, fig.width=9, fig.height=6-----------------------------
# Add statistical comparisons between groups
get_expression_barplot(
  vista,
  genes = top_up[1:4],
  display_id = "SYMBOL",
  log_transform = TRUE,
  stats_group = TRUE  # Enable statistical annotations
) + theme_minimal(base_size = 16)

## ----barplot-per-sample, fig.width=10, fig.height=6---------------------------
get_expression_barplot(
  vista,
  genes = top_up[1:2],
  display_id = "SYMBOL",
  by = "sample",
  sample_order = "group"
)+ theme(text = element_text(size = 16))

## ----barplot-comparison, fig.width=10, fig.height=6---------------------------
# Compare expression of both up- and down-regulated genes
selected_genes <- c(top_up[1:3], top_down[1:3])
get_expression_barplot(
  vista,
  genes = selected_genes,
  display_id = "SYMBOL",
  log_transform = TRUE,
  stats_group = TRUE
)+ theme(text = element_text(size = 16))

## ----boxplot-basic, fig.width=9, fig.height=6---------------------------------
get_expression_boxplot(
  vista,
  genes = top_up[1:4],
  display_id = "SYMBOL"
)+ theme(text = element_text(size = 16))

## ----boxplot-no-facet, fig.width=9, fig.height=6------------------------------
# All genes overlaid on same plot
get_expression_boxplot(
  vista,
  genes = top_up[1:3],
  display_id = "SYMBOL",
  facet_by = "none"
) + theme(text = element_text(size = 16))

## ----boxplot-facets, fig.width=9, fig.height=6--------------------------------
# Each gene in separate panel - must specify facet_by = "gene"
get_expression_boxplot(
  vista,
  genes = top_up[1:3],
  display_id = "SYMBOL",
  facet_by = "gene",
  facet_scales = "free_y"
) + theme(text = element_text(size = 16))

## ----boxplot-facets-stats, fig.width=10, fig.height=7-------------------------
# Each gene in separate panel WITH statistical comparisons
get_expression_boxplot(
  vista,
  genes = top_up[1:4],
  display_id = "SYMBOL",
  log_transform = TRUE,
  facet_by = "gene",
  facet_scales = "free_y",
  stats_group = TRUE,  # Add statistics to each gene panel
  p.label = "p.signif"
)+ theme(text = element_text(size = 16))

## ----boxplot-pooled, fig.width=8, fig.height=6--------------------------------
# Pool all genes together for group comparison with statistical test
get_expression_boxplot(
  vista,
  genes = top_up[1:5],
  display_id = "SYMBOL",
  log_transform = TRUE,
  pool_genes = TRUE,
  facet_by = "none",
  stats_group = TRUE,  # Required for statistical annotations
  p.label = "p.signif"
)+ theme(text = element_text(size = 16))

## ----boxplot-pvalues, fig.width=9, fig.height=6-------------------------------
# Show statistical comparisons between treatment groups
get_expression_boxplot(
  vista,
  genes = top_up[1:4],
  display_id = "SYMBOL",
  log_transform = TRUE,
  stats_group = TRUE,  # Enable statistical annotations
  p.label = "p.signif"
)+ theme(text = element_text(size = 16))

## ----violin-basic, fig.width=9, fig.height=6----------------------------------
get_expression_violinplot(
  vista,
  genes = top_up[1:4],
  display_id = "SYMBOL"
)+ theme(text = element_text(size = 16))

## ----violin-log, fig.width=9, fig.height=6------------------------------------
get_expression_violinplot(
  vista,
  genes = top_up[1:4],,
  display_id = "SYMBOL",
  value_transform = "none"
)+ theme(text = element_text(size = 16))

## ----violin-zscore, fig.width=8, fig.height=6---------------------------------
get_expression_violinplot(
  vista,
  genes = top_up[1:4],
  value_transform = "zscore"
)+ theme(text = element_text(size = 16))

## ----density-plot, fig.width=9, fig.height=6----------------------------------
get_expression_density(
  vista,
  genes = top_up[1:50],
  log_transform = TRUE
)+ theme(text = element_text(size = 16))

## ----scatter-plot, fig.width=8, fig.height=7----------------------------------
# Compare two samples
samples <- colnames(vista)
get_expression_scatter(
  vista,
  sample_x = samples[1],
  sample_y = samples[2],
  log_transform = TRUE,
  label_n = 50,display_id = "SYMBOL",label_size = 4
)+ theme(text = element_text(size = 16))

## ----line-plot, fig.width=9, fig.height=6-------------------------------------
get_expression_lineplot(
  vista,
  genes = top_up[1:3],
  log_transform = TRUE,display_id = "SYMBOL",
  by = "sample",facet_by = "none",
  group_column = "treatment",sample_group = c("Untreated","Dexamethasone")
)

## ----lollipop-plot, fig.width=9, fig.height=6---------------------------------
get_expression_lollipop(
  vista,
  genes = top_up[1:4],
  display_id = "SYMBOL",
  log_transform = TRUE
)

## ----lollipop-per-sample, fig.width=10, fig.height=6--------------------------
get_expression_lollipop(
  vista,
  genes = top_up[1:2],
  display_id = "SYMBOL",
  by = "sample",
  sample_order = "expression"
)

## ----joyplot, fig.width=8, fig.height=7, eval=requireNamespace("ggridges", quietly=TRUE)----
# Ridges by treatment group - shows distribution for each group
get_expression_joyplot(
  vista,
  genes = top_up[1:5],
  log_transform = TRUE,
  y_by = "group",      # Each treatment group gets a ridge
  color_by = "group"   # Color by treatment group
)

## ----joyplot-sample, fig.width=8, fig.height=9, eval=requireNamespace("ggridges", quietly=TRUE)----
# Ridges by individual sample - shows distribution for each sample
get_expression_joyplot(
  vista,
  genes = top_up[1:3],
  log_transform = TRUE,
  y_by = "sample",     # Each sample gets a ridge
  color_by = "group"   # Color by treatment group
)

## ----raincloud-expression-message, eval=!requireNamespace("ggrain", quietly=TRUE), echo=FALSE----
# cat("Package 'ggrain' is not installed; expression raincloud plot example is skipped.\n")

## ----raincloud-expression, fig.width=9, fig.height=6, eval=requireNamespace("ggrain", quietly=TRUE)----
get_expression_raincloud(
  vista,
  genes = top_up,
  value_transform = "log2",
  summarise = TRUE,
  facet_by = "none",
  id.long.var = "gene",
  stats_group = TRUE
)

## ----msigdb-up, eval=requireNamespace("msigdbr", quietly=TRUE)----------------
msig_up <- get_msigdb_enrichment(
  vista,
  sample_comparison = comp_names[1],
  regulation = "Up",
  msigdb_category = "H",  # Hallmark gene sets
  species = "Homo sapiens",
  from_type = "ENSEMBL"
)

# View top enriched pathways
if (!is.null(msig_up$enrich) && nrow(msig_up$enrich@result) > 0) {
  head(msig_up$enrich@result[, c("Description", "pvalue", "p.adjust", "Count")])
}

## ----msigdb-down, eval=requireNamespace("msigdbr", quietly=TRUE)--------------
msig_down <- get_msigdb_enrichment(
  vista,
  sample_comparison = comp_names[1],
  regulation = "Down",
  msigdb_category = "H",
  species = "Homo sapiens",
  from_type = "ENSEMBL"
)

if (!is.null(msig_down$enrich) && nrow(msig_down$enrich@result) > 0) {
  head(msig_down$enrich@result[, c("Description", "pvalue", "p.adjust", "Count")])
}

## ----enrichment-dotplot, fig.width=9, fig.height=7, eval=requireNamespace("msigdbr", quietly=TRUE)----
# VISTA's wrapper function
if (!is.null(msig_up$enrich) && nrow(msig_up$enrich@result) > 0) {
  get_enrichment_plot(msig_up$enrich, top_n = 10)
}

## ----enrichment-barplot, fig.width=9, fig.height=6, eval=requireNamespace("msigdbr", quietly=TRUE) && requireNamespace("enrichplot", quietly=TRUE)----
# Use generic barplot with enrichResult method
if (!is.null(msig_up$enrich) && nrow(msig_up$enrich@result) > 0) {
  barplot(msig_up$enrich, showCategory = 10)
}

## ----enrichment-dotplot-native, fig.width=9, fig.height=7, eval=requireNamespace("msigdbr", quietly=TRUE) && requireNamespace("enrichplot", quietly=TRUE)----
# Customized dotplot with more categories
if (!is.null(msig_up$enrich) && nrow(msig_up$enrich@result) > 0) {
  enrichplot::dotplot(msig_up$enrich, showCategory = 20, font.size = 12)
}

## ----enrichment-network, fig.width=10, fig.height=9, eval=requireNamespace("msigdbr", quietly=TRUE) && requireNamespace("enrichplot", quietly=TRUE)----
# Gene-concept network showing gene-pathway relationships
if (!is.null(msig_up$enrich) && nrow(msig_up$enrich@result) > 0) {
  enrichplot::cnetplot(msig_up$enrich, showCategory = 5)
}

## ----enrichment-chord-pathway, fig.width=8, fig.height=8, eval=requireNamespace("msigdbr", quietly=TRUE) && requireNamespace("circlize", quietly=TRUE)----
# Pathway-coloured chord diagram (no VISTA object needed)
if (!is.null(msig_up$enrich) && nrow(msig_up$enrich@result) > 0) {
  get_enrichment_chord(
    msig_up,
    top_n    = 6,
    color_by = "pathway",
    title    = "Gene-Pathway Membership (Up-regulated, Hallmark)"
  )
}

## ----enrichment-chord-fc, fig.width=8, fig.height=8, eval=requireNamespace("msigdbr", quietly=TRUE) && requireNamespace("circlize", quietly=TRUE)----
# Fold-change coloured chords
if (!is.null(msig_up$enrich) && nrow(msig_up$enrich@result) > 0) {
  get_enrichment_chord(
    msig_up,
    vista  = vista,
    sample_comparison = names(comparisons(vista))[1],
    top_n      = 6,
    color_by   = "foldchange",
    display_id = "SYMBOL",
    title      = "Gene-Pathway Chord (coloured by log2FC)"
  )
}

## ----enrichment-chord-hub, fig.width=8, fig.height=8, eval=requireNamespace("msigdbr", quietly=TRUE) && requireNamespace("circlize", quietly=TRUE)----
# Show only hub genes shared across 2+ pathways
if (!is.null(msig_up$enrich) && nrow(msig_up$enrich@result) > 0) {
  chord_result <- get_enrichment_chord(
    msig_up,
    vista    = vista,
    sample_comparison   = names(comparisons(vista))[1],
    top_n        = 8,
    min_pathways = 2,
    color_by     = "regulation",
    display_id   = "SYMBOL",
    title        = "Hub Genes Bridging Multiple Pathways"
  )

  # Inspect hub genes returned invisibly
  if (length(chord_result$hub_genes)) {
    cat("Hub genes:", paste(head(chord_result$hub_genes, 10), collapse = ", "), "\n")
  }
}

## ----pathway-genes, eval=requireNamespace("msigdbr", quietly=TRUE)------------
if (!is.null(msig_up$enrich) && nrow(msig_up$enrich@result) > 0) {
  pathway_gene_list <- get_pathway_genes(
    msig_up$enrich,
    top_n = 3,
    return_type = "list"
  )

  # Preview first few genes per pathway
  lapply(pathway_gene_list, head, n = 5)
}

## ----pathway-heatmap, fig.width=8, fig.height=10, eval=requireNamespace("msigdbr", quietly=TRUE) && requireNamespace("ComplexHeatmap", quietly=TRUE)----
if (!is.null(msig_up$enrich) && nrow(msig_up$enrich@result) > 0) {
  get_pathway_heatmap(
    x = vista,
    enrichment = msig_up$enrich,
    sample_group = c("Untreated", "Dexamethasone"),
    top_n = 2,
    gene_mode = "union",
    max_genes = 60,
    value_transform = "zscore",
    display_id = "SYMBOL",
    annotate_columns = TRUE,
    summarise_replicates = FALSE,
    show_row_names = FALSE
  )
}

## ----go-bp, eval=requireNamespace("clusterProfiler", quietly=TRUE)------------
go_bp <- get_go_enrichment(
  vista,
  sample_comparison = comp_names[1],
  regulation = "Up",
  ont = "BP",  # Biological Process
  species = "Homo sapiens",
  from_type = "ENSEMBL"
)

if (!is.null(go_bp$enrich) && nrow(go_bp$enrich@result) > 0) {
  head(go_bp$enrich@result[, c("Description", "pvalue", "p.adjust", "Count")], n = 10)
}

## ----go-plot, fig.width=9, fig.height=7, eval=requireNamespace("clusterProfiler", quietly=TRUE)----
if (!is.null(go_bp$enrich) && nrow(go_bp$enrich@result) > 0) {
  get_enrichment_plot(go_bp$enrich, top_n = 20)
}

## ----gsea-hallmark, eval=requireNamespace("msigdbr", quietly=TRUE) && requireNamespace("clusterProfiler", quietly=TRUE)----
# Run GSEA using VISTA's native function
gsea_results <- get_gsea(
  vista,
  sample_comparison = comp_names[1],
  set_type = "msigdb",
  from_type = "ENSEMBL",
  species = "Homo sapiens",
  msigdb_category = "H",  # Hallmark gene sets
  pvalueCutoff = 0.05,
  pAdjustMethod = "BH"
)

# Show results
if (!is.null(gsea_results$enrich) && nrow(gsea_results$enrich@result) > 0) {
  head(gsea_results$enrich@result[, c("Description", "NES", "pvalue", "p.adjust")], n = 10)
}

## ----gsea-go, eval=requireNamespace("clusterProfiler", quietly=TRUE)----------
# Run GSEA with GO terms
gsea_go <- get_gsea(
  vista,
  sample_comparison = comp_names[1],
  set_type = "go",
  from_type = "ENSEMBL",
  orgdb = org.Hs.eg.db,
  ont = "BP",  # Biological Process
  pvalueCutoff = 0.05
)

# Show results
if (!is.null(gsea_go$enrich) && nrow(gsea_go$enrich@result) > 0) {
  head(gsea_go$enrich@result[, c("Description", "NES", "pvalue", "p.adjust")], n = 10)
}

## ----gsea-dotplot, fig.width=9, fig.height=7, eval=requireNamespace("msigdbr", quietly=TRUE) && requireNamespace("clusterProfiler", quietly=TRUE)----
# Show all significant pathways using VISTA's visualization
if (!is.null(gsea_results$enrich) && nrow(gsea_results$enrich@result) > 0) {
  get_enrichment_plot(gsea_results$enrich, top_n = 15)
}

## ----gsea-pathway, fig.width=10, fig.height=6, eval=requireNamespace("msigdbr", quietly=TRUE) && requireNamespace("enrichplot", quietly=TRUE)----
# Show enrichment plot for the top pathway
if (!is.null(gsea_results$enrich) && nrow(gsea_results$enrich@result) > 0) {
  # Create GSEA enrichment plot with running enrichment score
  enrichplot::gseaplot2(
    gsea_results$enrich,
    geneSetID = 1,  # Top pathway
    title = gsea_results$enrich@result$Description[1],
    pvalue_table = TRUE
  )
}

## ----gsea-multi, fig.width=10, fig.height=8, eval=requireNamespace("msigdbr", quietly=TRUE) && requireNamespace("enrichplot", quietly=TRUE)----
# Show top 3 pathways together
if (!is.null(gsea_results$enrich) && nrow(gsea_results$enrich@result) > 0) {
  enrichplot::gseaplot2(
    gsea_results$enrich,
    geneSetID = 1:3,  # Top 3 pathways
    pvalue_table = TRUE,
    ES_geom = "line"
  )
}

## ----gsea-go-plot, fig.width=9, fig.height=7, eval=requireNamespace("clusterProfiler", quietly=TRUE)----
# Visualize GO GSEA results
if (!is.null(gsea_go$enrich) && nrow(gsea_go$enrich@result) > 0) {
  get_enrichment_plot(gsea_go$enrich, top_n = 15)
}

## ----kegg-up, eval=FALSE------------------------------------------------------
# kegg_up <- get_kegg_enrichment(
#   vista,
#   sample_comparison = comp_names[1],
#   regulation = "Up",
#   species = "Homo sapiens",
#   from_type = "ENSEMBL"
# )
# 
# if (!is.null(kegg_up$enrich) && nrow(kegg_up$enrich@result) > 0) {
#   head(kegg_up$enrich@result[, c("Description", "pvalue", "p.adjust", "Count")], n = 10)
# }

## ----kegg-down, eval=FALSE----------------------------------------------------
# kegg_down <- get_kegg_enrichment(
#   vista,
#   sample_comparison = comp_names[1],
#   regulation = "Down",
#   species = "Homo sapiens",
#   from_type = "ENSEMBL"
# )
# 
# if (!is.null(kegg_down$enrich) && nrow(kegg_down$enrich@result) > 0) {
#   head(kegg_down$enrich@result[, c("Description", "pvalue", "p.adjust", "Count")])
# }

## ----kegg-plot, fig.width=9, fig.height=6, eval=FALSE-------------------------
# if (!is.null(kegg_up$enrich) && nrow(kegg_up$enrich@result) > 0) {
#   get_enrichment_plot(kegg_up$enrich, top_n = 15)
# }

## ----foldchange-matrix--------------------------------------------------------
fc_matrix <- get_foldchange_matrix(vista)
head(fc_matrix, n = 10)

## ----fc-barplot-gene, fig.width=9, fig.height=6-------------------------------
get_foldchange_barplot(
  vista,
  genes = top_up[1:3],
  sample_comparisons = comp_names,
  display_id = "SYMBOL",
  facet_by = "gene"
)

## ----fc-lollipop-gene, fig.width=9, fig.height=6------------------------------
get_foldchange_lollipop(
  vista,
  sample_comparison = comp_names[1],
  genes = top_up[1:3],
  display_id = "SYMBOL",
  facet_by = "gene"
)

## ----raincloud-foldchange-message, eval=!requireNamespace("ggrain", quietly=TRUE), echo=FALSE----
# cat("Package 'ggrain' is not installed; fold-change raincloud plot example is skipped.\n")

## ----raincloud-foldchange, fig.width=9, fig.height=6, eval=requireNamespace("ggrain", quietly=TRUE)----
get_foldchange_raincloud(
  vista,
  sample_comparisons = comp_names,
  facet_by = "none",
  id.long.var = "gene_id",
  stats_group = TRUE
)

## ----fc-heatmap-basic, fig.width=7, fig.height=9, eval=requireNamespace("ComplexHeatmap", quietly=TRUE)----
get_foldchange_heatmap(vista)

## ----fc-heatmap-explicit, fig.width=7, fig.height=9, eval=requireNamespace("ComplexHeatmap", quietly=TRUE)----
# Select genes with large fold-changes
fc_genes <- rownames(fc_matrix)[abs(fc_matrix[, 1]) > 2][1:30]

get_foldchange_heatmap(
  vista,
  sample_comparisons = comp_names,
  genes = fc_genes,
  display_id = "SYMBOL"
)

## ----fc-heatmap-names, fig.width=7, fig.height=10, eval=requireNamespace("ComplexHeatmap", quietly=TRUE)----
get_foldchange_heatmap(
  vista,
  sample_comparisons = comp_names,
  genes = fc_genes[1:25],
  show_row_names = TRUE,
  display_id = "SYMBOL"
)

## ----fc-heatmap-custom, fig.width=7, fig.height=8, eval=requireNamespace("ComplexHeatmap", quietly=TRUE)----
# Use top upregulated genes
get_foldchange_heatmap(
  vista,
  sample_comparisons = comp_names,
  genes = top_up,
  show_row_names = TRUE,
  display_id = "SYMBOL"
)

## ----export-results, eval=FALSE-----------------------------------------------
# # Export complete DE table with annotations
# de_annotated <- merge(
#   comparisons(vista)[[1]],
#   as.data.frame(rowData(vista)),
#   by.x = "gene_id",
#   by.y = "row.names"
# )
# 
# # Write to CSV
# write.csv(
#   de_annotated,
#   file = "airway_dexamethasone_vs_untreated_DE_results.csv",
#   row.names = FALSE
# )
# 
# # Export significant genes only
# sig_genes <- de_annotated[de_annotated$regulation %in% c("Up", "Down"), ]
# write.csv(
#   sig_genes,
#   file = "airway_significant_DEGs.csv",
#   row.names = FALSE
# )

## ----save-vista, eval=FALSE---------------------------------------------------
# # Save the complete VISTA object for later use
# saveRDS(vista, file = "airway_vistaect.rds")
# 
# # Load it back
# # vista <- readRDS("airway_vistaect.rds")

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