---
title: "Project example"
author: "Anton Klåvus, Vilhelm Suksi"
date: "`r Sys.Date()`"
output:
  BiocStyle::html_document:
    toc: true
    toc_depth: 2
vignette: >
  %\VignetteIndexEntry{Project example}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
bibliography: references.bib
biblio-style: apalike
---

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

In this project example, core functionality of notame, notameViz and 
notameStats is demonstrated in unison with the associated protocol article 
[@notame].

# Project setup

## Set up path and logging
Let's set up a path and start logging our project. Many functions in the notame 
packages will log information when they have finished. This helps in monitoring 
the analysis process in case something goes wrong and in reviewing the results 
later. The log will also print to console, although we'll exclude console 
output for the rest of the document for brevity.

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

## Read data
The data is mock data with two balanced groups and two balanced time points. 
The time points also correspond to two batches. There are 50 samples with ten 
regularly interspersed QC samples and 80 features divided into four analytical 
modes.

```{r}
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"
```

The function `import_from_excel()` returns a SummarizedExperiment object with 
three main parts and corresponding accessors for streamlined data wrangling.

- `assay`: feature abundances across the samples  
- `colData`: sample information  
- `rowData`: feature information

The example data contains four analytical modes, so we separate them using 
`fix_object()`, which also conveniently checks the object for compatibility 
with all of notame. 

```{r}
modes <- fix_object(data, split_data = TRUE, assay.type = "abundances")
```

Now we have a list of objects, one per mode (LC column x 
ionization mode) or however many unique values were specified in
`rowData(data)$Split`.

```{r, results = 'markup'}
names(modes)
sapply(modes, class)
```

# Tabular data preprocessing

## Preprocessing by mode
Tabular data preprocessing is performed to complete the dataset by way of 
reducing unwanted variation and data preparation dependent on downstream 
methods. Steps performed separately for each mode include marking missing 
values as NA, flagging features with low detection rate [@qcguidelines], drift 
correction [@driftcorrection], and flagging of low-quality features 
[@qcguidelines]. Visualizations are drawn to monitor the process and explore 
the data globally. The `save_QC_plots()` wrapper does not include
feature-wise plots as it may not be feasible to inspect them at this 
stage of the analysis. For example, drift correction visualizations are best 
drawn for a subset of interesting features after feature selection.

```{r}
# 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
}
```

Feature-wise flagging information, quality metrics, and brief drift correction 
notes are included in the feature information after the above steps.

```{r, results = 'markup'}
rowData(processed[[1]])$DC_note
```

## Preprocessing for the complete dataset

Next, it is time to merge the modes together and visualize the complete dataset.

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

Then, QC samples can (and should) be removed, as they are no longer needed, and 
the dataset is visualized anew.

```{r}
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")
```

If there are any missing values in the data, they need to be imputed. Let's use 
random forest imputation to impute these values and visualize the dataset one 
last time before statistical analysis. Seed number should be set before random 
forest imputation to guarantee reproducibility of results.

```{r}
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")
```

By default, the imputation procedure only operates on good quality features, 
i.e. those that have not been flagged. To use flagged features in statistical 
analysis, they should be imputed as well. This can be achieved through a second 
round of imputation, now with all features included. This two-step imputation 
makes sure that low-quality features don't affect the imputation of quality 
features. Imputation could also be performed separately for each mode to reduce 
execution time, especially if the amounts of features still allows for good 
imputation results.

```{r}
base <- impute_rf(imputed, all_features = TRUE)
```

It is a good idea to save the merged and processed data, so experimenting with 
different statistical analyses becomes easier.

```{r}
save(base, file = paste0(ppath, "full_data.RData"))
```

# Feature selection

## Univariate analysis
First, we'll use a linear model to evaluate the features in terms of the 
probability of obtaining the observed abundances given the null hypothesis, 
namely that there is no difference in feature abundance across study groups.

```{r}
assay(base, "log") <- log(assay(base))
lm_results <- perform_lm(base, assay.type = "log", 
  formula_char = "Feature ~ Group")
base <- join_rowData(base, lm_results)
```

Next, univariate and supervised rankings of features could be combined in a 
final ranking. Such a final ranking would combine the qualities of univariate 
analysis and supervised learning. Next, carefully selected features would 
undergo labour-intensive scrutiny for biological meaning, for example by 
identification and pathway analysis.

That concludes our project example. The last thing to do is to write the 
results to an Excel file and finish the logging. Thanks!

```{r}
write_to_excel(base, file = paste0(ppath, "results.xlsx"))
finish_log()
```

# Session information

```{r, echo = FALSE, results = 'markup'}
sessionInfo()
```

# References