# TidyMass LC-MS workflow template. Reviewed 2026-09-23.
# Replace paths/groups and review every analysis choice before execution.
# This script does not contain study data or claim validated instrument parameters.
library(tidymass)
library(dplyr)

# 1. Inputs and raw processing -------------------------------------------------
raw_dir <- "analysis/MS1"                 # mzML/mzXML in QC/Subject subdirectories
sample_csv <- "analysis/sample_info.csv"
output_dir <- "analysis/results"
control_group <- "Control"
case_group <- "Treatment"
dir.create(output_dir, recursive = TRUE, showWarnings = FALSE)
stopifnot(dir.exists(raw_dir), file.exists(sample_csv))
writeLines(capture.output(sessionInfo()), file.path(output_dir, "sessionInfo.txt"))
set.seed(20260923)
massprocesser::process_data(
  path = raw_dir, polarity = "positive",
  ppm = 15, peakwidth = c(5, 30), snthresh = 10,
  threads = 2, min_fraction = 0.5, fill_peaks = FALSE,
  output_tic = TRUE, output_bpc = TRUE
)
result_file <- file.path(raw_dir, "Result", "object")
stopifnot(file.exists(result_file))
saved <- new.env()
load(result_file, envir = saved)
objects <- mget(ls(saved), envir = saved)
keep <- vapply(objects, function(x) inherits(x, "mass_dataset"), logical(1))
stopifnot(sum(keep) == 1L)
object <- objects[[which(keep)]]

# 2. Match metadata and inspect ------------------------------------------------
si <- read.csv(sample_csv, check.names = FALSE, stringsAsFactors = FALSE,
               colClasses = "character")
stopifnot(all(c("sample_id", "group", "class") %in% names(si)),
          !anyNA(si$sample_id), !anyDuplicated(si$sample_id))
ids <- colnames(massdataset::extract_expression_data(object))
stopifnot(setequal(ids, si$sample_id))
si <- si[match(ids, si$sample_id), , drop = FALSE]
if ("injection.order" %in% names(si)) {
  si$injection.order <- as.numeric(si$injection.order)
  stopifnot(!anyNA(si$injection.order))
}
# Reconstruct with the matched metadata before attaching MS2/annotations.
object <- massdataset::create_mass_dataset(
  expression_data = massdataset::extract_expression_data(object),
  sample_info = si,
  variable_info = massdataset::extract_variable_info(object)
)
saveRDS(object, file.path(output_dir, "01-imported.rds"))
massdataset::show_sample_missing_values(object)
massdataset::show_variable_missing_values(object)
# Optional report, with a working report-rendering environment:
# massqc::massqc_report(object, path = file.path(output_dir, "qc_before"), type = "html")

# 3. Filter, impute and normalize ---------------------------------------------
# Illustrative criteria; revise and document for your study.
expr <- massdataset::extract_expression_data(object)
ctrl <- si$sample_id[si$class == "Subject" & si$group == control_group]
case <- si$sample_id[si$class == "Subject" & si$group == case_group]
qc <- si$sample_id[si$class == "QC"]
stopifnot(length(ctrl) >= 3L, length(case) >= 3L)
keep_feature <- rowMeans(is.na(expr[, ctrl, drop = FALSE])) <= 0.5 |
                rowMeans(is.na(expr[, case, drop = FALSE])) <= 0.5
if (length(qc) > 0L) {
  keep_feature <- keep_feature & rowMeans(is.na(expr[, qc, drop = FALSE])) <= 0.2
}
write.csv(data.frame(variable_id = rownames(expr), keep = keep_feature),
          file.path(output_dir, "filter_decisions.csv"), row.names = FALSE)
# Exclude blanks from the illustrative normalized matrix.
stopifnot(any(keep_feature))
object <- object[which(keep_feature), si$sample_id[si$class %in% c("Subject", "QC")]]
stopifnot(nrow(object) > 1L)
cleaned <- masscleaner::impute_mv(object, method = "median")
normalized <- masscleaner::normalize_data(cleaned, method = "median")
saveRDS(normalized, file.path(output_dir, "02-normalized.rds"))
# For QC-based SVR/LOESS, first verify quantitative QC coverage and metadata.
# Review outliers/batches with acquisition records; do not remove PCA outliers blindly.

# 4. PCA and independent-group statistics ------------------------------------
# Keep scaled values out of fold-change calculation.
expr <- massdataset::extract_expression_data(normalized)
biological <- c(ctrl, case)
expr <- expr[, biological, drop = FALSE]
valid <- apply(expr, 1, function(x) all(is.finite(x)) && stats::sd(x) > 0)
stopifnot(sum(valid) >= 2L)
pca <- stats::prcomp(t(expr[valid, , drop = FALSE]), center = TRUE, scale. = TRUE)
write.csv(pca$x, file.path(output_dir, "pca_scores.csv"))
stats_object <- massstat::mutate_fc(normalized, ctrl, case, mean_median = "mean")
stats_object <- massstat::mutate_p_value(stats_object, ctrl, case,
  method = "t.test", p_adjust_methods = "BH")
feature_stats <- massdataset::extract_variable_info(stats_object)
write.csv(feature_stats, file.path(output_dir, "feature_statistics.csv"), row.names = FALSE)
saveRDS(stats_object, file.path(output_dir, "03-statistics.rds"))

# 5. Attach MS2 and annotate ---------------------------------------------------
# Replace these paths and run this block when these resources are available.
ms2_dir <- "analysis/MS2"
reference_rds <- "analysis/reference_database.rds"
if (dir.exists(ms2_dir) && file.exists(reference_rds)) {
  with_ms2 <- massdataset::mutate_ms2(stats_object, path = ms2_dir,
    column = "rp", polarity = "positive",
    ms1.ms2.match.mz.tol = 15, ms1.ms2.match.rt.tol = 30)
  database <- readRDS(reference_rds)
  stopifnot(inherits(database, "databaseClass"))
  annotated <- metid::annotate_metabolites_mass_dataset(
    object = with_ms2, database = database, polarity = "positive", column = "rp",
    ms1.match.ppm = 15, ms2.match.ppm = 30,
    rt.match.tol = NA, rt.match.weight = 0, candidate.num = 3, threads = 2)
  saveRDS(annotated, file.path(output_dir, "04-annotated.rds"))
  write.csv(massdataset::extract_annotation_table(annotated),
    file.path(output_dir, "annotation_candidates.csv"), row.names = FALSE)
} else {
  message("Annotation not run: provide experimental MS2 and a compatible reference database.")
}

# 6. Pathway analysis after candidate review ----------------------------------
# Prepare reviewed_kegg_ids.csv with one KEGG ID per row in a kegg_id column.
# Select using the prespecified statistical rule and reviewed annotation evidence.
# Do not count multiple features/adducts as independent compounds.
query_csv <- "analysis/reviewed_kegg_ids.csv"
pathway_rds <- "analysis/pathway_database.rds"
if (file.exists(query_csv) && file.exists(pathway_rds)) {
  query <- read.csv(query_csv, colClasses = "character")
  stopifnot("kegg_id" %in% names(query))
  query_ids <- unique(query$kegg_id[!is.na(query$kegg_id) & nzchar(query$kegg_id)])
  stopifnot(length(query_ids) > 0L, all(grepl("^C[0-9]{5}$", query_ids)))
  pathway_database <- readRDS(pathway_rds)
  enrichment <- metpath::enrich_kegg(query_id = query_ids, query_type = "compound",
    id_type = "KEGG", pathway_database = pathway_database,
    p_adjust_method = "BH", method = "hypergeometric", p_cutoff = 0.05, threads = 2)
  saveRDS(enrichment, file.path(output_dir, "05-pathway-enrichment.rds"))
  # Background is defined by the supplied pathway database, not automatically by measured compounds.
} else {
  message("Pathway analysis not run: provide reviewed compound IDs and a pathway database.")
}

# 7. Archive -----------------------------------------------------------------
writeLines(capture.output(sessionInfo()), file.path(output_dir, "sessionInfo.txt"))
# Back up this script, parameter choices, all inputs/reference versions and all results.