sts_untargeted_neg.Rmd

---
title: "Differential abundance analysis, untargeted approach"
subtitle: "Negative polarity"
author: "Christa Malfertheiner"
date: "29 October 2021"
output:
  BiocStyle::html_document:
    toc: true
    number_sections: false
    toc_float: true
bibliography: references.bib
csl: biomed-central.csl
---

```{r setup, echo = FALSE, results = "asis", warning = FALSE}
library(BiocStyle)
BiocStyle::markdown()
knitr::opts_chunk$set(echo = TRUE, message = FALSE, dev = c("png", "pdf"))
```

```{r parameters, echo = FALSE, warning = FALSE}
## Set general parameters
polarity <- "NEG" # specify "POS" or "NEG"
p.cut <- 0.05     # cut-off for significance.
m.cut <- 0.7      # cut-off for log2 fold change
set.seed(123)
## Setting golden ratio to save images
phi <- (1+sqrt(5))/2
FILE_NAME <- "sts_untargeted_neg"
## Define paths:
IMAGE_PATH <- paste0("images/", FILE_NAME, "/")
if (dir.exists(IMAGE_PATH)) unlink(IMAGE_PATH, recursive = TRUE, force = TRUE)
dir.create(IMAGE_PATH, recursive = TRUE, showWarnings = FALSE)
RDATA_PATH <- paste0("data/RData/", FILE_NAME, "/")
dir.create(RDATA_PATH, recursive = TRUE, showWarnings = FALSE)
RESULT_PATH <- paste0("data/results/", FILE_NAME, "/")
dir.create(RESULT_PATH, recursive = TRUE, showWarnings = FALSE)
```

# Introduction

In this document we perform the differential abundance analysis of the features
from short term stability data collected with mitra tips, with the aim of
finding significant storage-related features. This task is performed by
hypothesis testing, where we try to identify which metabolites have the most
different concentrations.
We follow an untargeted approach, the analysis comprises feature pre-filtering,
exploratory analysis and differential abundance analysis.


# Data import

First, we load the required packages and the data, after pre-processing and
normalization. The end result of these steps is a `SummarizedExperiment` that
contains aligned data, where features are grouped (after correspondence), and
that have undergone gap filling, normalization by the median, linear fitting and
per-feature between-batch normalization to remove any unwanted variability.
The `SummarizedExperiment` lets us store all the information regarding the
normalization steps in the form of `assays`, which we are still able to access
to proceed with the analysis.

```{r load-data, echo = FALSE, warning = FALSE}
library(xcms)
library(limma)
library(pheatmap)
library(writexl)
library(SummarizedExperiment)
library(RColorBrewer)
library(MsFeatures)
library(CompMetaboTools)
library(pander)
load("data/RData/sts_normalization_neg/res_neg.RData")
```

It is important now to remove the `QC` samples from the data set, because the
analysis has to be performed only on study samples; the `QC` samples, though
are still required to evaluate the goodness of the detected features, therefore
they will be stored in a separate `SummarizedExperiment` object that can be
accessed when needed.

We assign the colours as seen before.

```{r split-qc, echo = TRUE}
res_qc <- res_neg[, res_neg$storage == "QC"]
res_neg <- res_neg[, res_neg$storage != "QC"]
res_neg$storage <- factor(as.character(res_neg$storage))
## Define colors for the groups.
col_storage <- brewer.pal(12, name = "Paired")[c(6, 2, 10, 1, 9)]
names(col_storage) <- c("RT",           # red
                        "BAG",          # blue
                        "VACUUM",       # purple
                        "4C_BAG",       # light blue
                        "4C_VACUUM")
col_time <- brewer.pal(9, name = "OrRd")[c(1, 3, 5, 6, 7, 8, 9)]
names(col_time) <- c("2h",
                     "6h",
                     "1d",
                     "2d",
                     "3d",
                     "1w",
                     "2w")
## Setting golden ratio to save images
phi <- (1+sqrt(5))/2
```

The samples used in this analysis are listed below.

```{r, echo = FALSE, results = "asis"}
tab <- colData(res_neg)[, c("storage", "time")]
pandoc.table(as.data.frame(tab), style = "rmarkdown",
             caption = "Samples used in this analysis")
```

# Feature pre-filtering

Feature pre-filtering is an important step of data analysis that aims to
reduce as much as possible the random error that occurs during the measurement
of an analyte: during this process, features with high noise and features that
were detected in a low number of samples are removed. As a side effect, by
reducing the number of features that are being tested later, the pre-filtering
reduces also the loss of power by the subsequent adjustment for multiple
hypothesis testing.

This step is fundamental, though one must be careful not to pre-filter for a
characteristic that will be tested later: the pre-filtering must be, as a matter
of fact, independent from later analyses.

The first step of pre-filtering consists of removing features with high
technical variance; several methods have been developed to determine which
signals must be removed, the most common of which relies on the **relative
standard deviation (RSD)**, which is defined as the ratio between the standard
deviation and the mean:

$RSD = \dfrac{s_{i,qc}}{\bar{m}_{i,qc}}$

This value is calculated for each feature found in the pooled QC samples: when
this is higher than 30%, the feature is removed from the data set
[@broadhurstGuidelinesConsiderationsUse2018].

Another common approach is based on the **dispersion ratio (D-ratio)**: this is
defined as the ratio between the sample standard deviation for the pooled QC
samples and the sample standard deviation for the study samples (the former
expected to represent technical variance, the latter a combination of technical
and biological variance):

$D-ratio = \dfrac{s_{i,qc}}{s_{i,sample}}$

The interpretation of this value goes as follows: when the D-ratio is 0%, there
is no technical variance in the observed measurements, whereas a D-ratio of
100% represents only noise and no biological variance detected. A common cut-off
for the D-ratio is 0.5, aiming at keeping features whose variation in study
samples is twice as large as the one in QC samples
[@broadhurstGuidelinesConsiderationsUse2018].

For the present data set calculate the D-ratio **separately** for each source
(sample matrix) because the variance between sample matrices is expected to be
very large. For each feature the mean D-ratio across the 3 sample matrices is
used for the filtering.

```{r filter-rsd}
rsds <- rowRsd(assay(res_qc, "normalized_filled"))
dratios <- apply(
    log2(assay(res_qc, "normalized_filled")), 1, sd, na.rm = TRUE) /
    apply(log2(assay(res_neg, "normalized_filled")), 1, sd, na.rm = TRUE)
```

The distribution of RSD values and D-ratio is shown in the plot below:

```{r filter-rsd-plot, fig.cap = "Distribution of RSD values and D-ratios in the data set. The dashed vertical red line represents the cut-off value for the RSD and D-ratio, respectively.", echo = FALSE}
par(mfrow = c(1, 2), mar = c(5, 5, 5, 1))
plot(density(rsds, na.rm = TRUE), xlab = "RSD",
     main = "Distribution of RSD values",
     cex.main = 1.5, cex.lab = 1.5, cex.axis = 1.3)
abline(v = 0.3, col = "red", lty = 2)
plot(density(dratios, na.rm = TRUE), xlab = "D-ratio",
     main = "Distribution of D-ratios",
     cex.main = 1.5, cex.lab = 1.5, cex.axis = 1.3)
abline(v = 0.5, col = "red", lty = 2)
```

The plot below directly compares the RSD and D-ratio for each feature.

```{r filter-rsd-vs-dratio-plot, fig.path = IMAGE_PATH, fig.width = 5, fig.height = 5, fig.cap = "Direct comparison of RSD and D-ratios.", echo = FALSE}
plot(log2(rsds), log2(dratios), xlab = expression(log[2]~RSD),
     ylab = expression(log[2]~D-ratio), pch = 16, col = "#00000040")
abline(v = log2(0.3), col = "red", lty = 2)
abline(h = log2(0.5), col = "red", lty = 2)
```

The plot shows a correlation between RSD and D-ratios, though the two methods
are not interchangeable. Below we pre-filter the data using the D-ratio.

```{r do-filter}
res <- res_neg[which(dratios < 0.5), ]
```

This reduced the data set from `r length(dratios)` to `r nrow(res)` features.
Next, we discard the features that have not been identified in at least half of
the samples in any of the sample groups.

```{r filter-proportion}
keep <- moreAreValidThan(assay(res, "raw"), f = res$storage_time, prop = 0.5)
res <- res[keep, ]
```

The data set has been reduced from `r length(rsds)` to `r nrow(res)` features.


# Exploratory analysis: PCA

Next, we perform a PCA analysis: this allows us to gather information about any
possible similarities among the samples, based on the measured metabolite
intensities.

```{r pca-all}
pc <- prcomp(t(log2(assay(res, "normalized_filled_imputed"))),
             center = TRUE, scale. = FALSE)
```

```{r pca-plot, fig.path = IMAGE_PATH, fig.cap = "PCA of the samples based on feature intensities.", fig.width = 7 * phi, fig.height = 7, echo = FALSE}
par(mfrow = c(1, 2))
plot_pca(pc, col = "#00000080",
         bg = paste0(col_storage[as.character(res$storage)], 90),
         pc_x = 1, pc_y = 2, pch = 21)
plot_pca(pc, col = "#00000080",
         bg = paste0(col_storage[as.character(res$storage)], 90),
         pc_x = 3, pc_y = 4, pch = 21)
legend("bottomright", col = col_storage, legend = names(col_storage),
       title = "storage", pch = 16, ncol = 2)
```

PC3 and 4 separate room temperature samples from all other samples.

```{r, echo = FALSE}
png(paste0(IMAGE_PATH, "FIG-1-PCA-untargeted-neg.png"), width = 12, height = 6,
    units = "cm", res = 600, pointsize = 5)
par(mfrow = c(1, 2), mar = c(4.3, 4.5, 0.5, 0.5), cex.lab = 1.5, cex.axis = 1.5)
plot_pca(pc, col = "#00000080",
         bg = paste0(col_storage[as.character(res$storage)], "ce"),
         pc_x = 1, pc_y = 2, pch = 21, ylim = c(-200, 330), cex = 1.5)
legend("top", horiz = FALSE, legend = names(col_storage),
       pt.bg = col_storage, pch = 21, ncol = 3, col = "#00000080", cex = 1.1)
pchs <- rep(21, 5)
names(pchs) <- c("RT", "BAG", "4C_BAG", "VACUUM", "4C_VACUUM")
plot_pca(pc, col = "#00000080",
         bg = paste0(col_time[as.character(res$time)], "ce"),
         pc_x = 1, pc_y = 2, pch = pchs[as.character(res$storage)],
         ylim = c(-220, 330), cex = 1.5)
legend("top", horiz = TRUE, legend = names(col_time),
       pt.bg = col_time, pch = 21, cex = 1.1, col = "#00000080")
dev.off()

```


# Differential abundance analysis

In this section, we perform a differential abundance analysis to identify
features that have significantly different abundances between male and female
samples. The analysis is based on feature-wise multiple linear regression:
the aim of such analysis is to find the relationship between the independent
variables (age and sex) and the response variable (signal intensity).
In short, multiple linear regression is a form of linear regression that is used
when there are two or more predictors.
Multiple linear regression is preferred over separate simple linear regression
in order to avoid wrong predictions: this could happen because the input
variables may be correlated, which could lead to unsatisfactory results. The
formula for multiple regression model is:

$Y = \beta_0 + \beta_1x_{i1} + \beta_2x_{i2} \ldots, + \beta_px_{ip} + \epsilon$

where: \
- $Y =$ predicted value (dependent variable) \
- $\beta_0 =$ y intercept, constant term
- $\beta_1, \beta_2, \ldots, \beta_p =$ regression coefficients \
- $x_i =$ independent variables \
- $\epsilon =$ residuals.

The `limma` package contains the `lmFit` function, which calculates the linear
model that best describes the data. The results are stored in a `MArrayLM`
(Microarray Linear Model Fit) object.

This model, though, is not enough to accept whether or not there is actually a
relationship among the response and the independent variables, therefore we must
perform a hypothesis test: we define the *null hypothesis* as there not being
any differences in the abundances of metabolites between the different storage
methods. The *alternative hypothesis* is therefore defined when there are
differences in the intensities of the signals coming from the metabolites in
the different experimental groups.
To accept or reject the alternative hypothesis, it is necessary to calculate the
p-value: the function that allows us to add the results to the `MArrayLM` object
created before is `eBayes`. The `eBayes` function computes several statistics,
including the moderated t-test, which is defined as follows:

$\dfrac{d}{s + s_0}$

where \
- $d =$ difference in two group means ($m_1 - m_2$) \
- $s =$ pooled standard deviation \
- $s_0 =$ small constant (it depends on the variance within the group).

The constant is added to the denominator in order to avoid a division by an
extremely low number, which would of course increase the result of the statistic
falsely inducing us into rejecting the null hypothesis, thus considering the
difference as significant, when it is not. When performing multiple hypothesis
testing, though, there is a high chance of rejecting the null hypothesis when it
is true (type I error), thus a method to control the False Discovery Rate is
required: in this case we opted for Benjamini-Hochberg correction. In
conclusion, the alternative hypothesis is rejected when the adjusted p-value is
smaller than the confidence threshold that was set at the beginning of this
document to `r p.cut`. This means we accept `r p.cut * 100`% false positives
among the features called *significant*.

```{r analysis}
## Factor sample source, sex and age
storage <- factor(res$storage)
time <- factor(res$time)
storetime <- factor(res$storage_time)
## Fit the data to the desired design
dsgn <- model.matrix(~ 0 + storetime)
fit <- lmFit(log2(assay(res, "normalized_filled_imputed")), design = dsgn)
## Fit the actual contrasts of interest
contr_mat <- makeContrasts(
  RT6hvsRT2h = storetimeRT_6h - storetimeRT_2h,
  RT1dvsRT2h = storetimeRT_1d - storetimeRT_2h,
  RT2dvsRT2h = storetimeRT_2d - storetimeRT_2h,
  RT3dvsRT2h = storetimeRT_3d - storetimeRT_2h,
  RT1wvsRT2h = storetimeRT_1w - storetimeRT_2h,
  RT2wvsRT2h = storetimeRT_2w - storetimeRT_2h,
  BAG6hvsRT2h = storetimeBAG_6h - storetimeRT_2h,
  BAG1dvsRT2h = storetimeBAG_1d - storetimeRT_2h,
  BAG2dvsRT2h = storetimeBAG_2d - storetimeRT_2h,
  BAG3dvsRT2h = storetimeBAG_3d - storetimeRT_2h,
  BAG1wvsRT2h = storetimeBAG_1w - storetimeRT_2h,
  BAG2wvsRT2h = storetimeBAG_2w - storetimeRT_2h,
  VACUUM6hvsRT2h = storetimeVACUUM_6h - storetimeRT_2h,
  VACUUM1dvsRT2h = storetimeVACUUM_1d - storetimeRT_2h,
  VACUUM2dvsRT2h = storetimeVACUUM_2d - storetimeRT_2h,
  VACUUM3dvsRT2h = storetimeVACUUM_3d - storetimeRT_2h,
  VACUUM1wvsRT2h = storetimeVACUUM_1w - storetimeRT_2h,
  VACUUM2wvsRT2h = storetimeVACUUM_2w - storetimeRT_2h,
  `4CBAG6hvsRT2h` = storetime4C_BAG_6h - storetimeRT_2h,
  `4CBAG1dvsRT2h` = storetime4C_BAG_1d - storetimeRT_2h,
  `4CBAG2dvsRT2h` = storetime4C_BAG_2d- storetimeRT_2h,
  `4CBAG3dvsRT2h` = storetime4C_BAG_3d - storetimeRT_2h,
  `4CBAG1wvsRT2h` = storetime4C_BAG_1w- storetimeRT_2h,
  `4CBAG2wvsRT2h` = storetime4C_BAG_2w - storetimeRT_2h,
  `4CVACUUM6hvsRT2h` = storetime4C_VACUUM_6h - storetimeRT_2h,
  `4CVACUUM1dvsRT2h` = storetime4C_VACUUM_1d- storetimeRT_2h,
  `4CVACUUM2dvsRT2h` = storetime4C_VACUUM_2d - storetimeRT_2h,
  `4CVACUUM3dvsRT2h` = storetime4C_VACUUM_3d - storetimeRT_2h,
  `4CVACUUM1wvsRT2h` = storetime4C_VACUUM_1w - storetimeRT_2h,
  `4CVACUUM2wvsRT2h` = storetime4C_VACUUM_2w - storetimeRT_2h,
  levels = dsgn)
fit <- contrasts.fit(fit, contrasts = contr_mat)
fit <- eBayes(fit)
adjp <- apply(fit$p.value, 2, p.adjust, method = "BH")
tmp <- data.frame(
    coef = fit$coefficient,
    pvalue = fit$p.value,
    adjp = adjp,
    significant = adjp < p.cut & abs(fit$coefficient) > m.cut,
    check.names = FALSE
    )
avgs <- lapply(unique(res$storage_time), function(z) {
    rowMeans(log2(assay(
        res, "normalized_filled_imputed")[, res$storage_time == z]))
})
avgs <- do.call(cbind, avgs)
colnames(avgs) <- paste0("avg.", sub("_", "", unique(res$storage_time)))
rowData(res) <- cbind(rowData(res), tmp, avgs)
```

A table with the number of significant metabolites is shown below.

```{r table-sig, echo = FALSE, results = "asis"}
tab <- colSums(as.matrix(rowData(res)[, grep("significant",
                                             colnames(rowData(res)))]))
tab <- data.frame(comparison = sub("significant.", "", names(tab)), count = tab)
tab$storage <- c(rep("RT", 6), rep("BAG", 6), rep("VACUUM", 6),
                 rep("4C_BAG", 6), rep("4C_VACUUM", 6))
tab$time <- rep(c("6h", "1d", "2d", "3d", "1w", "2w"), 5)
rownames(tab) <- NULL
pandoc.table(tab[, c("comparison", "storage", "time", "count")],
             style = "rmarkdown",
             caption = paste0("Number of significant features of the in",
                              " total", nrow(res), "analyzed features."))
```

The number of significant features is then shown in a barplot:

```{r sig-features-barplot, echo = FALSE, fig.path = IMAGE_PATH, fig.width = 10 * phi, fig.cap = "Amount of significant features per storage condition over time"}
tab_time <- tab
tab_time$time <- factor(tab_time$time,
                        levels = c("6h", "1d", "2d", "3d", "1w", "2w"))
tab_time <- tab_time[order(as.integer(tab_time$time)), ]
x <- barplot(tab_time$count, space = c(rep(0, 5), 0.2, rep(0, 4), 0.2,
                                       rep(0, 4), 0.2, rep(0, 4), 0.2,
                                       rep(0, 4), 0.2, rep(0, 4)),
             col = col_storage[tab_time$storage], ylab = "count",
             main = "Number of significant features")
legend("topleft", pch = 22, col = "black", pt.bg = col_storage,
       legend = names(col_storage))
grid(nx = NA, ny = NULL)
mtext(at = vapply(split(x, tab_time$time), mean, numeric(1)),
      text = levels(tab_time$time), side = 1, cex = par("cex.axis"),
      line = 1.5)
```

Like already observed in positive mode,  the samples stored at room temperature
show the smallest changes in feature abundance at nearly all time points
measured in negative mode, too. Also, the amount of features with a significant
difference in abundance drops after one week as seen in positive mode.

We then calculated the percentage of the metabolome changing over storage time:

```{r sig-features-barplot-perc, echo = FALSE, fig.path = IMAGE_PATH, fig.width = 10 * phi, fig.cap = "Percentage of significant features per storage condition over time"}
col_storage <- col_storage[names(col_storage) != "QC"]
tab_time$percentage <- tab_time$count / nrow(res) * 100
x <- barplot(tab_time$percentage, space = c(rep(0, 5), 0.2, rep(0, 4), 0.2,
                                            rep(0, 4), 0.2, rep(0, 4), 0.2,
                                            rep(0, 4), 0.2, rep(0, 4)),
             col = col_storage[as.character(tab_time$storage)], ylab = "%",
             main = "% of features being significant")
legend("topleft", pch = 22, col = "black", pt.bg = col_storage,
       legend = names(col_storage))
grid(nx = NA, ny = NULL)
mtext(at = vapply(split(x, tab_time$time), mean, numeric(1)),
      text = levels(tab_time$time), side = 1, cex = par("cex.axis"),
      line = 1.5)
```

We plot the same information as lines instead of bar plots.

```{r sig-features-lines, echo = FALSE, fig.path = IMAGE_PATH, fig.width = 10 * phi, fig.cap = "Amount of significant features per storage condition over time"}
tab_time$days <- 0.25
tab_time$days[tab_time$time == "1d"] <- 1
tab_time$days[tab_time$time == "2d"] <- 2
tab_time$days[tab_time$time == "3d"] <- 3
tab_time$days[tab_time$time == "1w"] <- 7
tab_time$days[tab_time$time == "2w"] <- 14

tab_condition <- split(tab_time, tab_time$storage)
plot(NA, NA, main = "Number of significant features", xlim = c(0, 14),
     ylim = c(0, max(tab_time$count)), ylab = "count", xlab = "days")
for (i in seq_along(tab_condition)) {
    points(x = tab_condition[[i]]$days,
           y = tab_condition[[i]]$count,
           type = "b", pch = 21,
           bg = paste0(col_storage[names(tab_condition)[i]], 40),
           col = paste0(col_storage[names(tab_condition)[i]], 80))
}
grid()
legend("bottomright", pch = 21, col = col_storage,
       pt.bg = paste0(col_storage, 40),
       legend = names(col_storage), bg = "white")
```

And the same plot using percentages instead of absolute counts.

```{r sig-features-lines-perc, echo = FALSE, fig.path = IMAGE_PATH, fig.width = 10 * phi, fig.cap = "Percentage of significant features per storage condition over time"}
plot(NA, NA, main = "Percentage of significant features", xlim = c(0, 14),
     ylim = c(0, max(tab_time$percentage)), ylab = "%", xlab = "days")
for (i in seq_along(tab_condition)) {
    points(x = tab_condition[[i]]$days,
           y = tab_condition[[i]]$percentage,
           type = "b", pch = 21,
           bg = paste0(col_storage[names(tab_condition)[i]], 40),
           col = paste0(col_storage[names(tab_condition)[i]], 80))
}
grid()
legend("bottomright", pch = 21, col = col_storage,
       pt.bg = paste0(col_storage, 40),
       legend = names(col_storage), bg = "white")
```

```{r, echo = FALSE}

png(paste0(IMAGE_PATH, "FIG-1-perc-untargeted-neg.png"), width = 8,
    height = 4, units = "cm", res = 600, pointsize = 5)
par(mar = c(4.3, 4.5, 1.5, 0.5), cex.axis = 1.5, cex.lab = 1.5)
plot(NA, NA, main = "Percentage of significant features", xlim = c(0, 14),
     ylim = c(0, max(tab_time$percentage)), ylab = "%", xlab = "days")
for (i in seq_along(tab_condition)) {
    points(x = tab_condition[[i]]$days,
           y = tab_condition[[i]]$percentage,
           type = "b", pch = 21,
           bg = paste0(col_storage[names(tab_condition)[i]], 80),
           col = paste0(col_storage[names(tab_condition)[i]], "ce"),
           cex = 1.4)
}
grid()
legend("bottomright", pch = 21, col = col_storage,
       pt.bg = paste0(col_storage, 40),
       legend = names(col_storage), bg = "white")

dev.off()

```


In addition, we create a PCA plot after averaging the replicates.

```{r average-replicates, echo = FALSE}
#' Average
averageSE <- function(x, column = character(), mainAssay = character()) {
    if (!column %in% colnames(colData(x)))
        stop("Column '", "' not found in 'colData' of 'x'")
    f <- factor(colData(x)[, column], levels = unique(colData(x)[, column]))
    ## new colData: take the first element for each replicate.
    cd <- colData(x)[match(levels(f), f), ]
    rownames(cd) <- cd[, column]
    ## loop over the assays and average them.
    a <- lapply(assays(x), function(z) {
        z <- split.data.frame(t(z), f = f)
        z <- do.call(cbind, lapply(z, colMeans, na.rm = TRUE))
        z[is.na(z)] <- NA
        z
    })
    if (length(mainAssay)) {
        tmp <- split.data.frame(t(assay(x, mainAssay)), f = f)
        tmp <- do.call(cbind, lapply(tmp, function(y) {
            apply(y, MARGIN = 2, FUN = sd, na.rm = TRUE)
        }))
        tmp[is.na(tmp)] <- NA
        a[[paste0(mainAssay, "_sd")]] <- tmp
    }
    SummarizedExperiment(assays = a, rowData = rowData(x),
                         colData = cd, metadata = metadata(x))
}
## Average technical replicates:
res_avg <- averageSE(res, column = "storage_time")
```

```{r standards-pca-avg, echo = FALSE}
pc <- prcomp(t(log2(assay(res_avg, "normalized_filled_imputed"))),
                 center = TRUE, scale. = FALSE)
```

```{r standards-pca-plot-avg, fig.path = IMAGE_PATH, fig.cap = "PCA of the samples based on intensities of known compounds.", fig.width = 7 * phi, fig.height = 7, echo = FALSE}
par(mfrow = c(1, 2))
plot_pca(pc, col = "#00000080",
         bg = paste0(col_storage[as.character(res_avg$storage)], 90),
         pc_x = 1, pc_y = 2, pch = 21)
legend("topright", col = col_storage, legend = names(col_storage),
       title = "storage", pch = 16, ncol = 2)
plot_pca(pc, col = "#00000080",
         bg = paste0(col_storage[as.character(res_avg$storage)], 90),
         pc_x = 3, pc_y = 4, pch = 21)
```


```{r standards-pca-plot-time-avg, fig.path = IMAGE_PATH, fig.cap = "PCA of the samples based on intensities of known compounds.", fig.width = 7 * phi, fig.height = 7, echo = FALSE}
par(mfrow = c(1, 2))
pchs <- 21:25
names(pchs) <- c("RT", "BAG", "4C_BAG", "VACUUM", "4C_VACUUM")
plot_pca(pc, col = "#00000080",
         bg = paste0(col_time[as.character(res_avg$time)], 90),
         pc_x = 1, pc_y = 2, pch = pchs[as.character(res_avg$storage)])
legend("topright", pch = pchs, legend = names(pchs),
       title = "storage", ncol = 2)
plot_pca(pc, col = "#00000080",
         bg = paste0(col_time[as.character(res_avg$time)], 90),
         pc_x = 3, pc_y = 4, pch = pchs[as.character(res_avg$storage)])
```

PC1 seems to represent the time.

```{r, echo = FALSE}

png(paste0(IMAGE_PATH, "FIG-1-PCA-avg-untargeted-neg.png"), width = 12,
    height = 6, units = "cm", res = 600, pointsize = 5)
par(mfrow = c(1, 2), mar = c(4.3, 4.5, 0.5, 0.5), cex.lab = 1.5, cex.axis = 1.5)
plot_pca(pc, col = "#00000080",
         bg = paste0(col_storage[as.character(res_avg$storage)], "ce"),
         pc_x = 1, pc_y = 2, pch = 21, ylim = c(-200, 240), cex = 1.5)
legend("top", horiz = FALSE, legend = names(col_storage), cex = 1.1,
       pt.bg = col_storage, pch = 21, ncol = 3, col = "#00000080")
pchs <- rep(21, 5)
names(pchs) <- c("RT", "BAG", "4C_BAG", "VACUUM", "4C_VACUUM")
plot_pca(pc, col = "#00000080",
         bg = paste0(col_time[as.character(res_avg$time)], "ce"),
         pc_x = 1, pc_y = 2, pch = pchs[as.character(res_avg$storage)],
         ylim = c(-200, 240), cex = 1.5)
legend("top", horiz = TRUE, legend = names(col_time), cex = 1.1,
       pt.bg = col_time, pch = 21, col = "#00000080")
dev.off()

```

At last we create also a heatmap of the coefficients.

```{r standards-heatmap-coefficients, echo = FALSE, fig.path = IMAGE_PATH, fig.height = 8 * phi, fig.width = 8, fig.cap = "Heatmap of coefficients from the comparisons of each storage timepoint against 2h RT. For better visibility the color bar was cut at a value of 5."}
coefs <- as.matrix(rowData(res)[, grep("coef", colnames(rowData(res)))])
anns <- data.frame(storage = rep(c("RT", "BAG", "VACUUM",
                                   "4C_BAG", "4C_VACUUM"), each = 6),
                   time = rep(c("6h", "1d", "2d", "3d", "1w", "2w"), 5))
rownames(anns) <- colnames(coefs)
rownames(coefs) <- rowData(res)$name
brks <- seq(-5, 5, length.out = 101)
pheatmap(coefs, show_colname = FALSE, annotation_col = anns, breaks = brks,
         annotation_colors = list(storage = col_storage, time = col_time))
```

# Evaluation of variance between replicated measurements

The rather surprising decrease of the number of significant metabolites for
later time points seen in the barplots above might eventually be caused by an
increased variance at these time points which then affect the statistical
test. To evaluate this we compare next the total number of feature with an
absolute coefficient larger than a certain cut-off and then also the variance
between the 3 replicated measurements for each time points.

The table below lists the number of features with an more than 2-fold (absolute)
difference in abundance compared to the initial time point.

```{r , echo = FALSE, results = "asis"}
## Get the coefficients for all comparisons.
tmp <- as.matrix(rowData(res)[, grep("coef", colnames(rowData(res)))])
## Count the number of features with abs coef > 1
tab_fc <- apply(abs(tmp), MARGIN = 2, function(z) sum(z > 1, na.rm = TRUE))
tab_fc <- data.frame(comparison = sub("coef.", "", names(tab_fc)),
                     count = tab_fc)
tab_fc$storage <- c(rep("RT", 6), rep("BAG", 6), rep("VACUUM", 6),
                    rep("4C_BAG", 6), rep("4C_VACUUM", 6))
tab_fc$time <- rep(c("6h", "1d", "2d", "3d", "1w", "2w"), 5)
rownames(tab_fc) <- NULL
pandoc.table(tab_fc[, c("comparison", "storage", "time", "count")],
             style = "rmarkdown",
             caption = paste0("Number of features of the in",
                              " total", nrow(res), "analyzed features ",
                              "with a more that 2-fold difference in abundance"))
```

```{r fold-change-features-barplot, echo = FALSE, fig.path = IMAGE_PATH, fig.width = 10 * phi, fig.cap = "Amount of features with a more than two-fold difference in abundance per storage condition over time."}
tab_time <- tab_fc
tab_time$time <- factor(tab_time$time,
                        levels = c("6h", "1d", "2d", "3d", "1w", "2w"))
tab_time <- tab_time[order(as.integer(tab_time$time)), ]
x <- barplot(tab_time$count, space = c(rep(0, 5), 0.2, rep(0, 4), 0.2,
                                       rep(0, 4), 0.2, rep(0, 4), 0.2,
                                       rep(0, 4), 0.2, rep(0, 4)),
             col = col_storage[as.character(tab_time$storage)], ylab = "count",
             main = "Number of significant features")
legend("topleft", pch = 22, col = "black", pt.bg = col_storage,
       legend = names(col_storage))
grid(nx = NA, ny = NULL)
mtext(at = vapply(split(x, tab_time$time), mean, numeric(1)),
      text = levels(tab_time$time), side = 1, cex = par("cex.axis"),
      line = 1.5)
```

We next evaluate the average abundance of the replicates as well as their
standard deviation over time and condition.

```{r}
res$storage_time <- droplevels(res$storage_time)
sds <- lapply(levels(res$storage_time), function(z) {
    tmp <- res[, res$storage_time == z]
    apply(log2(assay(tmp, "normalized_filled_imputed")),
          MARGIN = 1, sd, na.rm = TRUE)
})
names(sds) <- levels(res$storage_time)
## Mean abundance
means <- lapply(levels(res$storage_time), function(z) {
    tmp <- res[, res$storage_time == z]
    apply(log2(assay(tmp, "normalized_filled_imputed")),
          MARGIN = 1, mean, na.rm = TRUE)
})
names(means) <- levels(res$storage_time)
## Median abundance
medians <- lapply(levels(res$storage_time), function(z) {
    tmp <- res[, res$storage_time == z]
    apply(log2(assay(tmp, "normalized_filled_imputed")),
          MARGIN = 1, median, na.rm = TRUE)
})
names(medians) <- levels(res$storage_time)
```

```{r sd-mean-features-boxplot, echo = FALSE, fig.path = IMAGE_PATH, fig.width = 5 * phi, fig.cap = "Distribution of average feature abundances and their standard deviation per time point/condition.", echo = FALSE}
## Define ordering index (6h, 1d, 2d, ...).
idx <- c(2, 8, 14, 20, 26,
         3, 9, 15, 21, 27,
         4, 10, 16, 22, 28,
         5, 11, 17, 23, 29,
         6, 12, 18, 24, 30,
         7, 13, 19, 25, 31)
par(mfrow = c(2, 1), mar = c(0, 4.5, 0, 0))
boxplot(means[idx], xaxt = "n", xlab = "",
        ylab = expression(log[2]~mean~abundance),
        col = col_storage[rep(1:5, 6)])
grid(nx = NA, ny = NULL)
abline(v = c(5.5, 10.5, 15.5, 20.5, 25.5))
par(mar = c(4.5, 4.5, 0, 0))
boxplot(sds[idx], xaxt = "n", xlab = "",
        ylab = expression(SD),
        col = col_storage[rep(1:5, 6)])
grid(nx = NA, ny = NULL)
abline(v = c(5.5, 10.5, 15.5, 20.5, 25.5))
axis(side = 1, at = c(3, 8, 13, 18, 23, 28),
     labels = c("6h", "1d", "2d", "3d", "1w", "2w"))
```

There seems to be no clear systematic increase of the variance or decrease of
absolute signal with time.

# Reduced data sets

Now, we split the data sets into two parts. The first contains all samples
that were stored at room temperature or in plastic bags filled with desiccants,
the second contains all samples stored at rt or in vacuum bags.

```{r split-sets, echo = TRUE}
bags_rt <- subset(
    res_avg, res_avg$storage %in% c("RT", "BAG", "4C_BAG"))
vacuum_rt <- subset(
    res_avg, res_avg$storage %in% c("RT", "VACUUM", "4C_VACUUM"))
```

Then, a PCA plot for each subset is created, beginning with the **BAG** subset:

```{r standards-pca-avg-bag, echo = FALSE}
pc <- prcomp(t(log2(assay(bags_rt, "normalized_filled_imputed"))),
                 center = TRUE, scale. = FALSE)
```

```{r standards-pca-plot-avg-bag, fig.path = IMAGE_PATH, fig.cap = "PCA of the samples based on intensities of known compounds.", fig.width = 7 * phi, fig.height = 7, echo = FALSE}
par(mfrow = c(1, 2))
plot_pca(pc, col = "#00000080",
         bg = paste0(col_storage[as.character(res_avg$storage)], 90),
         pc_x = 1, pc_y = 2, pch = 21)
plot_pca(pc, col = "#00000080",
         bg = paste0(col_storage[as.character(res_avg$storage)], 90),
         pc_x = 3, pc_y = 4, pch = 21)
legend("topleft", col = col_storage, legend = names(col_storage),
       title = "storage method: BAG", pch = 16, ncol = 2)
```

We proceed with the **VACUUM** subset:

```{r standards-pca-avg-vacuum, echo = FALSE}
pc <- prcomp(t(log2(assay(vacuum_rt, "normalized_filled_imputed"))),
                 center = TRUE, scale. = FALSE)
```

```{r standards-pca-plot-avg-vacuum, fig.path = IMAGE_PATH, fig.cap = "PCA of the samples based on intensities of known compounds.", fig.width = 7 * phi, fig.height = 7, echo = FALSE}
par(mfrow = c(1, 2))
plot_pca(pc, col = "#00000080",
         bg = paste0(col_storage[as.character(res_avg$storage)], 90),
         pc_x = 1, pc_y = 2, pch = 21)
plot_pca(pc, col = "#00000080",
         bg = paste0(col_storage[as.character(res_avg$storage)], 90),
         pc_x = 3, pc_y = 4, pch = 21)
legend("topleft", col = col_storage, legend = names(col_storage),
       title = "storage method: VACUUM", pch = 16, ncol = 2)
```

```{r export-to-excel}
library(writexl)

write_xlsx(as.data.frame(rowData(res_neg)),
           "mitra_short_term_stability_untar_neg.xlsx")
```

# Session information

```{r}
sessionInfo()
```