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

## ----results = 'markup'-------------------------------------------------------
library(notame)
library(notameViz)
library(notameStats)
library(dplyr)
ppath <- tempdir()
init_log(log_file = file.path(ppath, "log.txt"))

## -----------------------------------------------------------------------------
data <- import_from_excel(
  file = system.file("extdata", "toy_notame_set.xlsx", package = "notame"), 
  sheet = 1, corner_row = 11, corner_column = "H", 
  split_by = c("Column", "Ion_mode"))
  
names(assays(data)) <- "abundances"

## -----------------------------------------------------------------------------
modes <- fix_object(data, split_data = TRUE, assay.type = "abundances")

## ----results = 'markup'-------------------------------------------------------
names(modes)
sapply(modes, class)

## -----------------------------------------------------------------------------
# Initialize empty list for processed objects
processed <- list()
for (i in seq_along(modes)) {
  name <- names(modes)[i]
  mode <- modes[[i]]
  # Set all zero abundances to NA
  mode <- mark_nas(mode, value = 0)
  # Flag features with low detection rate
  mode <- flag_detection(mode, group = "Group")
  # Visualize data before drift correction
  save_QC_plots(mode, prefix = paste0(ppath, "figures/", name, "_ORIG"),
                 perplexity = 5, group = "Group", time = "Time", 
                 id = "Subject_ID", color = "Group")
  # Correct drift
  corrected <- correct_drift(mode)
  # Visualize data after drift correction
  save_QC_plots(corrected, prefix = paste0(ppath, "figures/", name, "_DRIFT"),
                perplexity = 5, group = "Group", time = "Time", 
                id = "Subject_ID", color = "Group")
  # Flag low-quality features
  corrected <- corrected %>% 
    assess_quality() %>% 
    flag_quality()
  # Visualize data after removal of low-quality features
  save_QC_plots(corrected, prefix = paste0(ppath, "figures/", name, "_CLEAN"),
                perplexity = 5, group = "Group", time = "Time",
                id = "Subject_ID", color = "Group")
  # Save result of iteration
  processed[[i]] <- corrected
}

## ----results = 'markup'-------------------------------------------------------
rowData(processed[[1]])$DC_note

## -----------------------------------------------------------------------------
merged <- merge_notame_sets(processed)
save_QC_plots(merged, prefix = paste0(ppath, "figures/_FULL"),
              group = "Group", time = "Time",
              id = "Subject_ID", color = "Group")

## -----------------------------------------------------------------------------
merged_no_qc <- drop_qcs(merged)
save_QC_plots(merged_no_qc, prefix = paste0(ppath, "figures/FULL_NO_QC"),
              group = "Group", time = "Time",
              id = "Subject_ID", color = "Group")

## -----------------------------------------------------------------------------
set.seed(2024)
imputed <- impute_rf(merged_no_qc)
save_QC_plots(imputed, prefix = paste0(ppath, "figures/FULL_IMPUTED"),
              group = "Group", time = "Time",
              id = "Subject_ID", color = "Group")

## -----------------------------------------------------------------------------
base <- impute_rf(imputed, all_features = TRUE)

## -----------------------------------------------------------------------------
save(base, file = paste0(ppath, "full_data.RData"))

## -----------------------------------------------------------------------------
assay(base, "log") <- log(assay(base))
lm_results <- perform_lm(base, assay.type = "log", 
  formula_char = "Feature ~ Group")
base <- join_rowData(base, lm_results)

## -----------------------------------------------------------------------------
write_to_excel(base, file = paste0(ppath, "results.xlsx"))
finish_log()

## ----echo = FALSE, results = 'markup'-----------------------------------------
sessionInfo()