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

## ----installation, eval=FALSE-------------------------------------------------
# if (!requireNamespace("BiocManager", quietly = TRUE)) {
#   install.packages("BiocManager")
# }
# 
# BiocManager::install("LIPIDIFy")

## ----load_package-------------------------------------------------------------
library(LIPIDIFy)

## ----launch-app, eval=FALSE---------------------------------------------------
# launch_lipidomics_app()

## ----load-data----------------------------------------------------------------
file_path <- system.file("extdata", "ST001359_lipidomics.csv", package = "LIPIDIFy")

raw_data <- load_lipidomics_data(
  file_path,
  metadata_columns = c("Sample Name", "Sample Group")
)

cat("Samples  :", nrow(raw_data$numeric_data), "\n")
cat("Lipids   :", ncol(raw_data$numeric_data), "\n")
cat("Groups   :", paste(unique(raw_data$metadata$`Sample Group`), collapse = ", "), "\n")

## ----peek-data----------------------------------------------------------------
head(raw_data$metadata, 6)
raw_data$numeric_data[, 1:4]

## ----classify-----------------------------------------------------------------
lipid_names <- colnames(raw_data$numeric_data)
classification <- classify_lipids(lipid_names)

head(classification, 8)

## ----classify-summary---------------------------------------------------------
cat("Lipid groups:\n")
print(sort(table(classification$LipidGroup), decreasing = TRUE))

cat("\nFatty-acid saturation:\n")
print(table(classification$Saturation))

## ----raw-boxplot, fig.cap="Per-sample intensity distributions before normalization. Boxes are coloured by group."----
raw_list <- list(numeric_data = raw_data$numeric_data, metadata = raw_data$metadata)

visualize_raw_data_improved(
  raw_list,
  plot_type    = "boxplot",
  view_mode    = "sample",
  metadata     = raw_data$metadata,
  group_column = "Sample Group"
)

## ----raw-density, fig.cap="Intensity density curves before normalization, one curve per sample."----
visualize_raw_data_improved(
  raw_list,
  plot_type    = "density",
  view_mode    = "sample",
  metadata     = raw_data$metadata,
  group_column = "Sample Group"
)

## ----norm-methods-------------------------------------------------------------
get_normalization_methods()

## ----norm-vsn, eval = requireNamespace("vsn", quietly = TRUE)-----------------
synthetic_obj_vsn <- load_lipidomics_data_from_df(generate_example_data())

vsn_norm <- apply_normalizations(synthetic_obj_vsn$numeric_data, "VSN")
dim(vsn_norm) # samples x lipids - orientation and dimnames preserved

# Distinct from Log2Median, which is a different method:
l2m_norm <- apply_normalizations(synthetic_obj_vsn$numeric_data, "Log2Median")
isTRUE(all.equal(vsn_norm, l2m_norm))

## ----norm-compare, fig.cap="Side-by-side comparison of two normalization pipelines applied to the same raw, synthetic data.", fig.height=9----
synthetic_df <- generate_example_data()
synthetic_obj <- load_lipidomics_data_from_df(synthetic_df)

pipeline1 <- apply_normalizations(synthetic_obj$numeric_data, c("TIC", "Log2"))
pipeline2 <- apply_normalizations(synthetic_obj$numeric_data, c("PQN", "Log2"))

p1 <- create_pipeline_plot(pipeline1,
  title     = "Pipeline 1: TIC + Log2",
  metadata  = synthetic_obj$metadata,
  plot_type = "boxplot"
)

p2 <- create_pipeline_plot(pipeline2,
  title     = "Pipeline 2: PQN + Log2",
  metadata  = synthetic_obj$metadata,
  plot_type = "boxplot"
)

gridExtra::grid.arrange(p1, p2, ncol = 1)

## ----normalize----------------------------------------------------------------
norm_data <- apply_normalizations(raw_data$numeric_data, methods = "Median")

dim(norm_data) # samples x lipids - same dimensions as raw data

## ----norm-boxplot, fig.cap="Per-sample distributions after median-centring."----
norm_list <- list(numeric_data = norm_data, metadata = raw_data$metadata)

visualize_raw_data_improved(
  norm_list,
  plot_type    = "boxplot",
  view_mode    = "sample",
  metadata     = raw_data$metadata,
  group_column = "Sample Group"
)

## ----impute-------------------------------------------------------------------
cat("Missing values in this dataset:", sum(is.na(norm_data)), "\n")
cat("Imputation methods available:\n")
print(get_imputation_methods())

## ----impute-apply-------------------------------------------------------------
# For demonstration: introduce 5% artificial missing values
demo_norm <- norm_data
set.seed(42)
demo_norm[sample(length(demo_norm), size = round(0.05 * length(demo_norm)))] <- NA
cat(
  "Missingness introduced:", sum(is.na(demo_norm)),
  sprintf("(%.1f%%)\n", 100 * mean(is.na(demo_norm)))
)

imputed_norm <- impute_missing_values(demo_norm, method = "half_min")
cat("Missing after imputation:", sum(is.na(imputed_norm)), "\n")

## ----batch-correction, eval=FALSE---------------------------------------------
# # Simulate batch labels (real data would have this in the metadata file)
# raw_data$metadata$Batch <- rep(c("Run1", "Run2"), each = 3)
# 
# norm_corrected <- correct_batch_effects(
#   data_matrix  = norm_data,
#   metadata     = raw_data$metadata,
#   batch_column = "Batch",
#   group_column = "Sample Group",
#   method       = "limma"
# )
# 
# # Verify: PCA should now show clustering by biology, not by batch
# pca_corrected <- perform_pca(norm_corrected, raw_data$metadata, "Sample Group")
# create_pca_plot_with_ellipses(
#   pca_data           = pca_corrected$pca_data,
#   variance_explained = pca_corrected$variance_explained,
#   ellipse_type       = "confidence",
#   title              = "PCA after Batch Correction"
# )

## ----pca, fig.cap="PCA of the normalized ST001359 data."----------------------
pca_res <- perform_pca(norm_data, raw_data$metadata, "Sample Group")

create_pca_plot_with_ellipses(
  pca_data           = pca_res$pca_data,
  variance_explained = pca_res$variance_explained,
  ellipse_type       = "none",
  show_sample_labels = TRUE,
  title              = "PCA - Control vs. SDC Inhibition"
)

## ----contrasts----------------------------------------------------------------
groups <- unique(raw_data$metadata$`Sample Group`)
contrasts <- create_default_contrasts(groups)
print(contrasts)

## ----diff-analysis------------------------------------------------------------
diff_res <- perform_differential_analysis(
  data_matrix    = norm_data,
  metadata       = raw_data$metadata,
  group_column   = "Sample Group",
  contrasts_list = contrasts,
  method         = "limma"
)

## ----diff-summary-------------------------------------------------------------
cat(sprintf("%-20s  %5s  %5s  %4s  %5s\n", "Contrast", "Total", "Sig", "Up", "Down"))
cat(strrep("-", 48), "\n")

for (nm in names(diff_res$results)) {
  res <- diff_res$results[[nm]]
  n_sig <- sum(res$adj.P.Val < 0.05, na.rm = TRUE)
  n_up <- sum(res$adj.P.Val < 0.05 & res$logFC > 0, na.rm = TRUE)
  n_down <- sum(res$adj.P.Val < 0.05 & res$logFC < 0, na.rm = TRUE)
  cat(sprintf("%-20s  %5d  %5d  %4d  %5d\n", nm, nrow(res), n_sig, n_up, n_down))
}

## ----top-results--------------------------------------------------------------
first_contrast <- names(diff_res$results)[1]
res_df <- diff_res$results[[first_contrast]]
res_df <- res_df[order(res_df$adj.P.Val), ]

top10 <- utils::head(res_df[, c("logFC", "AveExpr", "adj.P.Val")], 10)
top10$logFC <- round(top10$logFC, 3)
top10$AveExpr <- round(top10$AveExpr, 2)
top10$adj.P.Val <- signif(top10$adj.P.Val, 3)

cat("Top 10 features in:", first_contrast, "\n\n")
print(top10)

## ----volcano, fig.cap="Volcano plot coloured by lipid class."-----------------
create_volcano_plot_labeled(
  results             = diff_res$results[[first_contrast]],
  title               = paste("Volcano:", first_contrast),
  logfc_threshold     = 1,
  pval_threshold      = 0.05,
  top_labels          = 10,
  classification_data = classification,
  color_by            = "LipidGroup"
)

## ----heatmap, fig.cap="Heatmap of significant lipids. Rows (lipids) and columns (samples) are hierarchically clustered. Colour scale is z-scored per lipid.", fig.height=7----
sig_lipids <- rownames(res_df)[
  !is.na(res_df$adj.P.Val) &
    res_df$adj.P.Val < 0.05 &
    abs(res_df$logFC) > 1
]

if (length(sig_lipids) >= 5) {
  top_n <- min(30, length(sig_lipids))
  avail <- intersect(sig_lipids[seq_len(top_n)], colnames(norm_data))
  hm_mat <- t(norm_data[, avail, drop = FALSE])

  create_heatmap_robust(
    data_matrix  = hm_mat,
    metadata     = raw_data$metadata,
    group_column = "Sample Group",
    top_n        = top_n,
    title        = paste("Significant lipids:", first_contrast)
  )
} else {
  message("Fewer than 5 significant lipids at these thresholds.")
}

## ----lipid-expression, fig.cap="Per-sample abundance of the 3 most significant lipids."----
top_lipids <- utils::head(rownames(res_df)[!is.na(res_df$adj.P.Val)], 3)
cat("Selected lipids:\n")
print(top_lipids)

expr_result <- create_lipid_expression_barplot(
  data_matrix     = norm_data,
  metadata        = raw_data$metadata,
  selected_lipids = top_lipids,
  group_column    = "Sample Group",
  data_type       = "normalized"
)

if (inherits(expr_result, "ggplot")) {
  print(expr_result)
} else {
  for (lipid_plot in expr_result) print(lipid_plot)
}

## ----enrichment---------------------------------------------------------------
enrich_res <- perform_enrichment_analysis(
  results_list        = diff_res$results,
  classification_data = classification,
  min_set_size        = 3,
  max_set_size        = 500
)

cat("Categories per contrast:", paste(names(enrich_res[[1]]), collapse = ", "), "\n")

## ----enrichment-table---------------------------------------------------------
grp_enrich <- enrich_res[[first_contrast]][["LipidGroup"]]

if (!is.null(grp_enrich) && nrow(grp_enrich) > 0) {
  grp_enrich <- grp_enrich[order(grp_enrich$pval), ]
  out <- utils::head(grp_enrich[, c("pathway", "NES", "pval", "padj", "size")], 5)
  out$NES <- round(out$NES, 3)
  out$pval <- signif(out$pval, 3)
  out$padj <- signif(out$padj, 3)
  print(out, row.names = FALSE)
} else {
  cat("No enrichment results available.\n")
}

## ----enrichment-plot, fig.cap="Enrichment dot plot for the first contrast. Dot size reflects set size; colour reflects statistical significance."----
if (!is.null(grp_enrich) && nrow(grp_enrich) > 0) {
  create_enrichment_dotplot(
    enrichment_data = grp_enrich,
    title           = paste("Lipid Group Enrichment:", first_contrast),
    max_pathways    = 10
  )
}

## ----summarized-experiment----------------------------------------------------
if (requireNamespace("SummarizedExperiment", quietly = TRUE)) {
  se <- SummarizedExperiment::SummarizedExperiment(
    assays  = list(normalized = t(norm_data)),
    colData = raw_data$metadata
  )
  print(se)
}

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