peak_detection_qa.Rmd

---
title: "Quality assessment of the CHRIS untargeted peak detection"
author: "Johannes Rainer"
output:
  BiocStyle::html_document:
    toc: true
    number_sections: false
    toc_float: true
bibliography: references.bib
csl: biomed-central.csl
---

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

**Modified**: `r file.info("peak_detection_qa.Rmd")$mtime`<br />
**Compiled**: `r date()`

```{r settings, echo = FALSE}
## Set general options
options(useFancyQuotes = FALSE)
set.seed(123)

## Define paths:
FILE_NAME <- "peak_detection_qa"
## Path to save the images; remove all old images.
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)
## Path to store RData files
RDATA_PATH <- paste0("data/RData/", FILE_NAME, "/")
dir.create(RDATA_PATH, recursive = TRUE, showWarnings = FALSE)

## Get the number of cpus allocated or fall back to 3
ncores <- as.integer(Sys.getenv("SLURM_JOB_CPUS_PER_NODE", 4))

rt_cut <- 340

MZML_PATH <- "/data/massspec/mzML/"
if (!dir.exists(MZML_PATH))
    stop("Can not find the directory with the mzML files: ", MZML_PATH)
```

# Introduction

In this document we perform a quality assessment of the peak detection step of
the CHRIS untargeted metabolomics data (defined and performed in
[peak_detection.Rmd](peak_detection.Rmd). The quality assessment is performed on
the QC samples from each batch of the data set, which represent repeated
measurement of the same original samples (the *CHRIS POOL*). Also, retention
time drifts and intensity differences are evaluated on chromatographic peaks
representing signal from internal standards or a set of known compounds (the
lab-specific standards).

# Data import

We next load all required libraries and the data. The full data set is reduced
to only QC samples.

```{r libraries}
library(xcms)
library(Rdisop)
library(pander)
library(RColorBrewer)
library(MetaboCoreUtils)
library(CompMetaboTools)

register(bpstart(MulticoreParam(ncores - 1L)))

load("data/RData/peak_detection/data_pos_ref.RData")
blnk_pos <- filterFile(data_pos, file = which(data_pos$type == "BLANK"))
qc_pos <- filterFile(data_pos, file = which(data_pos$type == "QC"))
rm(data_pos)
load("data/RData/peak_detection/data_neg_ref.RData")
blnk_neg <- filterFile(data_neg, file = which(data_neg$type == "BLANK"))
qc_neg <- filterFile(data_neg, file = which(data_neg$type == "QC"))
rm(data_neg)

## Load the base peak chromatogram data
load("data/RData/peak_detection/bpc_pos.RData")
bpc_qc_pos <- bpc_pos[, which(bpc_pos$type == "QC")]
rm(bpc_pos)
load("data/RData/peak_detection/bpc_neg.RData")
bpc_qc_neg <- bpc_neg[, which(bpc_neg$type == "QC")]
rm(bpc_neg)

## At last define colors for each batch.
btches <- union(unique(qc_pos$batch), unique(qc_neg$batch))
## Ordering can be tricky because we have 12 and 104, so 104 would be
## ordered before 12.
btch_mod <- sapply(btches, function(x) {
    tmp <- strsplit(x, "_")[[1L]]
    paste0(tmp[1L], "_",
           ifelse(nchar(tmp[2L]) == 2L, yes = "0", no = ""),
           tmp[2L])
})
btches <- btches[order(btch_mod)]
col_batch <- rainbow(length(btches))
names(col_batch) <- btches
```


# Global quality assessment

For the global quality assessment we evaluate and compare the signal in QC
samples across all batches. At first we compare the average number of identified
chromatographic peaks in QC samples per batch.

```{r peak-counts, echo = FALSE, fig.path = IMAGE_PATH, fig.width = 7, fig.height = 7, fig.cap = "Average number of detected peaks in QC samples per batch. Batches are ordered by measurement time."}
label_year <- function(x, ypos = 10000) {
    yrs <- vapply(strsplit(names(x), "_"),
                  function(x) as.integer(x[1]), integer(1))
    year_change <- c(1, which(diff(yrs) > 0))
    abline(v = year_change, lty = 2)
    text(x = year_change, y = rep(ypos, length(year_change)),
         labels = yrs[year_change + 1], pos = 4)
}

pks_pos <- table(factor(as.integer(chromPeaks(qc_pos)[, "sample"]),
                        levels = seq_along(fileNames(qc_pos))))
pks_pos <- split(pks_pos, f = qc_pos$batch)
pks_pos <- vapply(pks_pos, mean, numeric(1))
## negative
pks_neg <- table(factor(as.integer(chromPeaks(qc_neg)[, "sample"]),
                        levels = seq_along(fileNames(qc_neg))))
pks_neg <- split(pks_neg, f = qc_neg$batch)
pks_neg <- vapply(pks_neg, mean, numeric(1))
par(mfrow = c(2, 1), mar = c(1, 4.5, 1, 1))
plot(y = pks_pos, x = seq_along(pks_pos), bg = col_batch[names(pks_pos)],
     main = "positive polarity", xaxt = "n", ylab = "average peak count",
     pch = 21, xlab = "")
grid()
label_year(pks_pos, ypos = 10000)
plot(y = pks_neg, x = seq_along(pks_neg), bg = col_batch[names(pks_neg)],
     main = "negative polarity", xaxt = "n", ylab = "average peak count",
     pch = 21, xlab = "")
grid()
label_year(pks_neg, ypos = 9000)
```

The average number of identified peaks in QC samples varies considerably between
batches (over time). The highest number was detected in the last four batches
from 2017. The number of peaks per individual QC sample matches well with the
average number per batch (data not shown here). Note that this is related to the
average intensity across batches which we plot below.

```{r peak-count-individual-qc-samples, echo = FALSE}
## Peak counts for individual QC samples.
png("peak-count-qc-samples.png", width = 18, height = 8, units = "cm",
    pointsize = 4, res = 300)
year_change <- c(1, which(diff(as.numeric(qc_pos$year)) > 0) + 1)
intos <- split(log2(chromPeaks(qc_pos)[, "into"]),
               chromPeaks(qc_pos)[, "sample"])
par(mfrow = c(2, 1), mar = c(2, 4.3, 1.5, 0.5))
tmp <- barplot(lengths(intos), col = col_batch[qc_pos$batch],
               ylab = "peak count", main = "QC samples, positive polarity",
               xaxt = "n", border = "#00000040")
abline(v = tmp[year_change, 1], lty = 2)
text(x = tmp[year_change, ], y = rep(10000, length(year_change)),
     labels = qc_pos$year[year_change], pos = 4)
## Negative polarity
year_change <- c(1, which(diff(as.numeric(qc_neg$year)) > 0) + 1)
intos <- split(log2(chromPeaks(qc_neg)[, "into"]),
               chromPeaks(qc_neg)[, "sample"])
tmp <- barplot(lengths(intos), col = col_batch[qc_neg$batch],
               ylab = "peak count", main = "QC samples, negative polarity",
               xaxt = "n", border = "#00000040")
abline(v = tmp[year_change, 1], lty = 2)
text(x = tmp[year_change, ], y = rep(9000, length(year_change)),
     labels = qc_neg$year[year_change], pos = 4)
dev.off()
```

```{r bpc-boxplots-qc-samples, fig.path = IMAGE_PATH, echo = FALSE, fig.width = 12, fig.height = 8, fig.cap = "Per-sample base peak chromatogram (BPC) intensity distribution."}
## Generate the BPC signal distribution.
ints <- lapply(bpc_qc_pos, intensity)
ints <- lapply(ints, log2)
par(mar = c(2, 4.3, 1.5, 0.5), mfrow = c(2, 1))
tmp <- boxplot(ints, varwidth = TRUE, col = col_batch[bpc_qc_pos$batch],
               outline = FALSE, main = "QC samples",
               ylab = expression(log[2]~base~peak~signal),
               border = "#00000040")
year_change <- c(1, which(diff(as.numeric(bpc_qc_pos$year)) > 0) + 1)
abline(v = year_change, lty = 2)
text(x = year_change, y = rep(19, length(year_change)),
     labels = bpc_qc_pos$year[year_change], pos = 4)
## Neg polarity
ints <- lapply(bpc_qc_neg, intensity)
ints <- lapply(ints, log2)
tmp <- boxplot(ints, varwidth = TRUE, col = col_batch[bpc_qc_neg$batch],
               outline = FALSE, main = "QC samples, negative polarity",
               ylab = expression(log[2]~base~peak~signal),
               border = "#00000040")
year_change <- c(1, which(diff(as.numeric(bpc_qc_neg$year)) > 0) + 1)
abline(v = year_change, lty = 2)
text(x = year_change, y = rep(18, length(year_change)),
     labels = bpc_qc_neg$year[year_change], pos = 4)
```

Positive and negative mode data was not generated at the same (or consecutive)
time which explains the differences in drifts observed for the base peak signal
in positive and negative data (see Figure above).

At last we create also the BPC for positive and negative data.

```{r bpc-plot-qc-samples, fig.path = IMAGE_PATH, echo = FALSE, fig.width = 12, fig.height = 8, fig.cap = "Base peak chromatograms."}
par(mar = c(4.5, 4.3, 1.5, 0.5), mfrow = c(2, 1))
plot(bpc_qc_pos, col = paste0(col_batch[bpc_qc_pos$batch], 20),
     main = "BPC, QC samples, positive polarity")
plot(bpc_qc_neg, col = paste0(col_batch[bpc_qc_neg$batch], 20),
     main = "BPC, QC samples, negative polarity")
```

The BPC look more or less similar, but some drifts are observable.

Next we specifically evaluate the signal from internal standards and lab
standards.


## Internal standards


Next we evaluate the signal of peaks for internal standards in the data set. To
this end we first load the definition of internal standards (i.e. their expected
m/z and retention time), identify chromatographic peaks in the data set
potentially matching these and plot this data (code not shown).

```{r is, echo = FALSE, eval = TRUE}
is_info <- read.table(
    "https://raw.githubusercontent.com/EuracBiomedicalResearch/lcms-standards/master/data/internal_standards.txt",
    sep = "\t", header = TRUE, as.is = TRUE)
rownames(is_info) <- is_info$name
is_info <- is_info[order(is_info$name), ]
is_info$mzneut = NA
is_info$mz_ion_pos = NA
is_info$mz_ion_neg = NA
for (i in seq(nrow(is_info))) {
    if (grepl("C", is_info$formula[i])){
        is_info$mzneut[i] <- getMolecule(
            as.character(is_info$formula[i]))$exactmass
    } else
        is_info$mzneut[i] = as.numeric(is_info$formula[i])
    ## Calculate also the m/z
    if (!is.na(is_info$POS[i]))
        is_info$mz_ion_pos[i] <-
            mass2mz(is_info$mzneut[i], adduct = is_info$POS[i])[1, 1]
    if (!is.na(is_info$NEG[i]))
        is_info$mz_ion_neg[i] <-
            mass2mz(is_info$mzneut[i], adduct = is_info$NEG[i])[1, 1]
}
is_pos <- is_info[!is.na(is_info$mz_ion_pos), ]
is_pos$mz_ion <- is_pos$mz_ion_pos
is_neg <- is_info[!is.na(is_info$mz_ion_neg), ]
is_neg$mz_ion <- is_neg$mz_ion_neg

plot_chrs <- function(x, dir, prefix = "IS-", suffix = "_POS", rect = TRUE) {
    sample_colors <- col_batch[x$batch]
    for (i in seq_len(nrow(x))) {
        cmpname <- fData(x)$name[i]
        filename <- paste0(dr, prefix, gsub("%", "p", cmpname, fixed = TRUE),
                           suffix, ".png")
        png(file = filename, width = 16, height = 8, units = "cm",
            res = 300, pointsize = 4)
        chr <- x[i, ]
        pks <- chromPeaks(chr)
        plot(chr, col = "#00000040",
             main = paste0(cmpname, ": ",
                           format(mz(chr)[1], digits = 6), "-",
                           format(mz(chr)[2], digits = 6)),
             peakCol = paste0(sample_colors[pks[, "column"]], 50),
             peakBg = paste0(sample_colors[pks[, "column"]], 10))
        ## This is just to indicate identified peaks with a rectangle.
        if (rect)
            xcms:::.add_chromatogram_peaks(
                       chr, pks,
                       col = paste0(sample_colors[pks[, "column"]], 10),
                       bg = NA,
                       type = "rectangle")
        abline(v = fData(x)$RT[i])
        dev.off()
    }
}

is_pos_mz_rt <- chromPeakArea(qc_pos, mz = is_pos$mz_ion, rt = is_pos$RT,
                              diffRt = 40, ppm = 10)
## is_pos_mz_rt <- peak_regions_for_mz_rt(
##     qc_pos, mz = is_pos$mz_ion, rt = is_pos$RT)
is_pos <- is_pos[!is.na(is_pos_mz_rt[, "mzmin"]), ]
is_pos_mz_rt <- is_pos_mz_rt[!is.na(is_pos_mz_rt[, "mzmin"]), ]
## Extract the ion chromatogram for each standard
rtr <- is_pos_mz_rt[, c("rtmin", "rtmax")]
rtr[, 1] <- rtr[, 1] - 10
rtr[, 2] <- rtr[, 2] + 10
is_pos_chr <- chromatogram(qc_pos, mz = is_pos_mz_rt[, c("mzmin", "mzmax")],
                           rt = rtr, aggregationFun = "max")
fData(is_pos_chr) <- cbind(
    fData(is_pos_chr), is_pos[, c("name", "abbreviation", "formula",
                                  "POS", "RT")])

## Plot the extracted ion chromatograms
dr <- paste0(IMAGE_PATH, "peakdetection/")
dir.create(dr, recursive = TRUE, showWarnings = FALSE)
plot_chrs(is_pos_chr, dr)

## Define the rt window for standards that look ~ OK.
is_pos$rtmin <- NA_real_
is_pos$rtmax <- NA_real_
is_pos$rtmin[is_pos$name == "Glycine (13C2, 99%; 15N, 99%)"] <- 165
is_pos$rtmax[is_pos$name == "Glycine (13C2, 99%; 15N, 99%)"] <- 183
is_pos$rtmin[is_pos$name == "L-Alanine (13C3, 99%; 15N, 99%)"] <- 160
is_pos$rtmax[is_pos$name == "L-Alanine (13C3, 99%; 15N, 99%)"] <- 180
is_pos$rtmin[is_pos$name == "L-Arginine HCl (13C6, 99%; 15N4, 99%)"] <- 165
is_pos$rtmax[is_pos$name == "L-Arginine HCl (13C6, 99%; 15N4, 99%)"] <- 220
is_pos$rtmin[is_pos$name == "L-Aspartic acid (13C4, 99%; 15N, 99%)"] <- 170
is_pos$rtmax[is_pos$name == "L-Aspartic acid (13C4, 99%; 15N, 99%)"] <- 210
is_pos$rtmin[is_pos$name == "L-Cystine (13C6, 99%; 15N2, 99%)"] <- 205
is_pos$rtmax[is_pos$name == "L-Cystine (13C6, 99%; 15N2, 99%)"] <- 230
is_pos$rtmin[is_pos$name == "L-Glutamic acid (13C5, 99%; 15N, 99%)"] <- 158
is_pos$rtmax[is_pos$name == "L-Glutamic acid (13C5, 99%; 15N, 99%)"] <- 200
is_pos$rtmin[is_pos$name == "L-Histidine HCl H2O (13C6; 15N3, 99%)"] <- 165
is_pos$rtmax[is_pos$name == "L-Histidine HCl H2O (13C6; 15N3, 99%)"] <- 240
is_pos$rtmin[is_pos$name == "L-Lysine 2HCl (13C6, 99%; 15N2, 99%)"] <- 172
is_pos$rtmax[is_pos$name == "L-Lysine 2HCl (13C6, 99%; 15N2, 99%)"] <- 215
is_pos$rtmin[is_pos$name == "L-Methionine (13C5, 99%; 15N, 99%)"] <- 150
is_pos$rtmax[is_pos$name == "L-Methionine (13C5, 99%; 15N, 99%)"] <- 165
is_pos$rtmin[is_pos$name == "L-Phenylalanine (13C9, 99%; 15N, 99%)"] <- 135
is_pos$rtmax[is_pos$name == "L-Phenylalanine (13C9, 99%; 15N, 99%)"] <- 160
is_pos$rtmin[is_pos$name == "L-Proline (13C5, 99%; 15N, 99%)"] <- 158
is_pos$rtmax[is_pos$name == "L-Proline (13C5, 99%; 15N, 99%)"] <- 185
is_pos$rtmin[is_pos$name == "L-Serine (13C3, 99%; 15N, 99%)"] <- 175
is_pos$rtmax[is_pos$name == "L-Serine (13C3, 99%; 15N, 99%)"] <- 195
is_pos$rtmin[is_pos$name == "L-Threonine (13C4, 99%; 15N, 99%)"] <- 165
is_pos$rtmax[is_pos$name == "L-Threonine (13C4, 99%; 15N, 99%)"] <- 190
is_pos$rtmin[is_pos$name == "L-Tyrosine (13C9, 99%; 15N, 99%)"] <- 155
is_pos$rtmax[is_pos$name == "L-Tyrosine (13C9, 99%; 15N, 99%)"] <- 185
is_pos <- is_pos[!is.na(is_pos$rtmin), ]

## Re-create the m/z and rt range matrix
is_pos_mz_rt <- chromPeakArea(qc_pos, mz = is_pos$mz_ion, ppm = 10,
                              rt = is_pos[, c("rtmin", "rtmax")], diffRt = 0)
is_pos_chr <- chromatogram(qc_pos, mz = is_pos_mz_rt[, c("mzmin", "mzmax")],
                           rt = is_pos_mz_rt[, c("rtmin", "rtmax")],
                           aggregationFun = "max", include = "apex_within")
fData(is_pos_chr) <- cbind(
    fData(is_pos_chr), is_pos[, c("name", "abbreviation", "formula",
                                  "POS", "RT")])

## For each IS (row), in each sample (column) keep the chrom peak with the
## highest peak.
is_pos_chr <- filterChromPeaks(is_pos_chr, method = "keepTop",
                               n = 1L, order = "maxo")
save(is_pos_chr, file = paste0(RDATA_PATH, "is_pos_chr.RData"))
## Plot again.
plot_chrs(is_pos_chr, dr, prefix = "IS-", suffix = "-refined_POS", rect = FALSE)


## --------------------------------------------------
## Repeat for negative polarity.
is_neg_mz_rt <- chromPeakArea(qc_neg, mz = is_neg$mz_ion, rt = is_neg$RT,
                              diffRt = 40, ppm = 10)
is_neg <- is_neg[!is.na(is_neg_mz_rt[, "mzmin"]), ]
is_neg_mz_rt <- is_neg_mz_rt[!is.na(is_neg_mz_rt[, "mzmin"]), ]
## Extract the ion chromatogram for each standard
rtr <- is_neg_mz_rt[, c("rtmin", "rtmax")]
rtr[, 1] <- rtr[, 1] - 10
rtr[, 2] <- rtr[, 2] + 10
is_neg_chr <- chromatogram(qc_neg, mz = is_neg_mz_rt[, c("mzmin", "mzmax")],
                           rt = rtr, aggregationFun = "max")
fData(is_neg_chr) <- cbind(
    fData(is_neg_chr), is_neg[, c("name", "abbreviation", "formula",
                                  "NEG", "RT")])

## Plot the extracted ion chromatograms
plot_chrs(is_neg_chr, dr, prefix = "IS-", suffix = "_NEG")

## Define the rt window for standards that look ~ OK.
is_neg$rtmin <- NA_real_
is_neg$rtmax <- NA_real_
is_neg$rtmin[is_neg$name == "Creatinine (N-Methhyl_D3)"] <- 120
is_neg$rtmax[is_neg$name == "Creatinine (N-Methhyl_D3)"] <- 145
is_neg$rtmin[is_neg$name == "L-Alanine (13C3, 99%; 15N, 99%)"] <- 175
is_neg$rtmax[is_neg$name == "L-Alanine (13C3, 99%; 15N, 99%)"] <- 187
is_neg$rtmin[is_neg$name == "L-Arginine HCl (13C6, 99%; 15N4, 99%)"] <- 175
is_neg$rtmax[is_neg$name == "L-Arginine HCl (13C6, 99%; 15N4, 99%)"] <- 200
is_neg$rtmin[is_neg$name == "L-Aspartic acid (13C4, 99%; 15N, 99%)"] <- 170
is_neg$rtmax[is_neg$name == "L-Aspartic acid (13C4, 99%; 15N, 99%)"] <- 200
is_neg$rtmin[is_neg$name == "L-Cystine (13C6, 99%; 15N2, 99%)"] <- 202
is_neg$rtmax[is_neg$name == "L-Cystine (13C6, 99%; 15N2, 99%)"] <- 230
is_neg$rtmin[is_neg$name == "L-Glutamic acid (13C5, 99%; 15N, 99%)"] <- 165
is_neg$rtmax[is_neg$name == "L-Glutamic acid (13C5, 99%; 15N, 99%)"] <- 200
is_neg$rtmin[is_neg$name == "L-Histidine HCl H2O (13C6; 15N3, 99%)"] <- 175
is_neg$rtmax[is_neg$name == "L-Histidine HCl H2O (13C6; 15N3, 99%)"] <- 205
is_neg$rtmin[is_neg$name == "L-Lysine 2HCl (13C6, 99%; 15N2, 99%)"] <- 175
is_neg$rtmax[is_neg$name == "L-Lysine 2HCl (13C6, 99%; 15N2, 99%)"] <- 205
is_neg$rtmin[is_neg$name == "L-Methionine (13C5, 99%; 15N, 99%)"] <- 155
is_neg$rtmax[is_neg$name == "L-Methionine (13C5, 99%; 15N, 99%)"] <- 170
is_neg$rtmin[is_neg$name == "L-Phenylalanine (13C9, 99%; 15N, 99%)"] <- 140
is_neg$rtmax[is_neg$name == "L-Phenylalanine (13C9, 99%; 15N, 99%)"] <- 161
is_neg$rtmin[is_neg$name == "L-Serine (13C3, 99%; 15N, 99%)"] <- 172
is_neg$rtmax[is_neg$name == "L-Serine (13C3, 99%; 15N, 99%)"] <- 190
is_neg$rtmin[is_neg$name == "L-Threonine (13C4, 99%; 15N, 99%)"] <- 163
is_neg$rtmax[is_neg$name == "L-Threonine (13C4, 99%; 15N, 99%)"] <- 183
is_neg$rtmin[is_neg$name == "L-Tyrosine (13C9, 99%; 15N, 99%)"] <- 160
is_neg$rtmax[is_neg$name == "L-Tyrosine (13C9, 99%; 15N, 99%)"] <- 175
is_neg$rtmin[is_neg$name == "L-Valine (13C5, 99%; 15N, 99%)"] <- 160
is_neg$rtmax[is_neg$name == "L-Valine (13C5, 99%; 15N, 99%)"] <- 175
is_neg <- is_neg[!is.na(is_neg$rtmin), ]

## Re-create the m/z and rt range matrix
is_neg_mz_rt <- chromPeakArea(qc_neg, mz = is_neg$mz_ion, ppm = 10,
                              rt = is_neg[, c("rtmin", "rtmax")], diffRt = 0)
is_neg_chr <- chromatogram(
    qc_neg, mz = is_neg_mz_rt[, c("mzmin", "mzmax")],
    rt = is_neg_mz_rt[, c("rtmin", "rtmax")],
    aggregationFun = "max", include = "apex_within")
fData(is_neg_chr) <- cbind(
    fData(is_neg_chr), is_neg[, c("name", "abbreviation", "formula",
                                  "NEG", "RT")])

## For each IS (row), in each sample (column) keep the chrom peak with the
## highest peak.
is_neg_chr <- filterChromPeaks(is_neg_chr, method = "keepTop",
                               n = 1L, order = "maxo")
save(is_neg_chr, file = paste0(RDATA_PATH, "is_neg_chr.RData"))
## Plot again.
plot_chrs(is_neg_chr, dr, prefix = "IS-", suffix = "-refined_NEG", rect = FALSE)

```

Next we quantify the signal for the internal standards for each sample. For
samples in which more than one chromatographic peak is present in the retention
time region, the one with the higher signal is used. We're next evaluating the
retention time drifts and signal/abundance differences for the such defined
internal standard data sets (positive and negative polarity).

### Evaluating retention time drifts

In this section we evaluate drifts in retention times for the internal
standards. We calculate first for each standard the average retention time and
its standard deviation within each batch and calculate next the standard
deviation of the average retention times across batches.

```{r rt-sd-mean}
load(paste0(RDATA_PATH, "is_pos_chr.RData"))
load(paste0(RDATA_PATH, "is_neg_chr.RData"))

#' Per standard and batch summaries (mean and sd for into and rt
batch_summary <- function(x) {
    tmp <- chromPeaks(x)
    tmp <- data.frame(pData(x)[tmp[, "sample"],
                               c("batch", "year", "month", "day")], tmp)
    ## split by row/standard
    tmpl <- split(tmp, tmp$row)
    res <- lapply(tmpl, function(z) {
        btch <- split(z, factor(z$batch,
                                levels = unique(z$batch, fromLast = TRUE)))
        do.call(rbind, lapply(btch, function(vals) {
            data.frame(batch = vals$batch[1L],
                       year = vals$year[1L],
                       month = vals$month[1L],
                       day = vals$day[1L],
                       row = vals$row[1L],
                       rt_mean = mean(vals$rt, na.rm = TRUE),
                       rt_sd = sd(vals$rt, na.rm = TRUE),
                       into_mean = mean(vals$into, na.rm = TRUE),
                       into_sd = sd(vals$into, na.rm = TRUE),
                       into_mean_log2 = mean(log2(vals$into), na.rm = TRUE),
                       into_sd_log2 = sd(log2(vals$into), na.rm = TRUE))
        }))
    })
    res <- do.call(rbind, res)
    res <- cbind(fData(x)[res$row, ], res)
    rownames(res) <- NULL
    res
}

batch_sum_pos <- batch_summary(is_pos_chr)
batch_sum_neg <- batch_summary(is_neg_chr)
```

```{r is-rt-within-between-batch, fig.path = IMAGE_PATH, fig.width = 8, fig.height = 4, fig.cap = "Evaluation of within and between batch retention time shifts. Each point represents one internal standard. X-axis: mean within-batch standard deviation of retention times. Y-axis: standard deviation of within-batch mean retention times.", echo = FALSE}
par(mfrow = c(1, 2))
tmp <- split(batch_sum_pos, batch_sum_pos$name)
## Calculate mean across within batch sd and
tmp_pos <- do.call(rbind, lapply(tmp, function(z) {
    c(mean_sd = mean(z$rt_sd, na.rm = TRUE),
      sd_mean = sd(z$rt_mean, na.rm = TRUE))
}))
plot(tmp_pos[, 1], tmp_pos[, 2], main = "positive polarity",
     xlab = "average sd of RT", ylab = "sd of average RT")
tmp <- split(batch_sum_neg, batch_sum_neg$name)
## Calculate mean across within batch sd and
tmp_neg <- do.call(rbind, lapply(tmp, function(z) {
    c(mean_sd = mean(z$rt_sd, na.rm = TRUE),
      sd_mean = sd(z$rt_mean, na.rm = TRUE))
}))
plot(tmp_neg[, 1], tmp_neg[, 2], main = "negative polarity",
     xlab = "average sd of RT", ylab = "sd of average RT")
```

As expected, the within-batch variation of retention times is on average much
smaller than the variation of the per-batch average retention time across the
data set. Three of the internal standards show exceptionally large variation
across batches. These are `r paste0(rownames(tmp_pos[order(tmp_pos[,
2]),])[1:3], collapse = ", ")` for positive and `r
paste0(rownames(tmp_neg[order(tmp_neg[, 2]),])[1:3], collapse = ", ")` for
negative polarity.

We next evaluate the retention time shifts for the individual standards. We base
this on the average retention time of each standard in each batch (since we've
shown above that the retention time shifts within each batch are much smaller
than between batches).

```{r is-rt-per-analyte-boxplot, fig.path = IMAGE_PATH, echo = FALSE, fig.width = 10, fig.height = 8, fig.cap = "Variation of internal standards' retention times across batches."}
batch_sum_pos_analyte <- split(batch_sum_pos, batch_sum_pos$name)
batch_sum_neg_analyte <- split(batch_sum_neg, batch_sum_neg$name)

rt_pos <- lapply(batch_sum_pos_analyte, function(z) z$rt_mean)
rt_pos <- rt_pos[order(vapply(rt_pos, mean, numeric(1)))]
rt_neg <- lapply(batch_sum_neg_analyte, function(z) z$rt_mean)
rt_neg <- rt_neg[order(vapply(rt_neg, mean, numeric(1)))]

par(mfrow = c(1, 2), mar = c(9, 4, 1, 0.5))
boxplot(rt_pos, las = 2, ylab = "RT", main = "positive polarity")
grid()
boxplot(rt_neg, las = 2, ylab = "RT", main = "negative polarity")
grid()
```

The variation in retention times seems to be larger for analytes eluting between
185 and 195 seconds. Next we investigate whether retention time shift are
*constant* per batch, i.e. whether all standards within a batch show the same
retention time shift.

```{r is-rt-per-analyte-lines, fig.path = IMAGE_PATH, echo = FALSE, fig.width = 10, fig.height = 12, fig.cap = "Per individual internal standard (within-batch average) retention time across batches."}
par(mfrow = c(2, 1), mar = c(7, 4.5, 1, 0.5))
btches <- unique(is_pos_chr$batch)
yl <- range(lapply(rt_pos, range))
xl <- c(1, length(btches))
plot(3, 3, pch = NA, xlim = xl, ylim = yl, xlab = "", xaxt = "n",
     ylab = "RT", main = "positive polarity")
grid()
tmp <- lapply(batch_sum_pos_analyte,
              function(z) points(x = match(z$batch, btches), y = z$rt_mean,
                                 type = "l", col = "#00000060"))
axis(side = 1, at = seq_along(btches),
     label = btches, las = 2)

btches <- unique(is_neg_chr$batch)
yl <- range(lapply(rt_neg, range))
xl <- c(1, length(btches))
plot(3, 3, pch = NA, xlim = xl, ylim = yl, xlab = "", xaxt = "n",
     ylab = "RT", main = "negative polarity")
grid()
tmp <- lapply(batch_sum_neg_analyte,
              function(z) points(x = match(z$batch, btches), y = z$rt_mean,
                                 type = "l", col = "#00000060"))
axis(side = 1, at = seq_along(btches),
     label = btches, las = 2)

```

The lines representing the retention times of individual internal standards show
about the same pattern across batches. Thus retention times of individual
standards appear to be affected by a similar shift within each batch. It should
hence be possible to align the samples without too large problems.


### Evaluating signal abundances

We next evaluate the variances of signal abundances within and between batches.

```{r is-into-within-between-batch, fig.path = IMAGE_PATH, fig.width = 8, fig.height = 4, fig.cap = "Evaluation of within and between batch intensity differences. Each point represents one internal standard. X-axis: mean within-batch standard deviation of intensities. Y-axis: standard deviation of within-batch mean intensities.", echo = FALSE}
par(mfrow = c(1, 2))
tmp <- split(batch_sum_pos, batch_sum_pos$name)
## Calculate mean across within batch sd and
tmp_pos <- do.call(rbind, lapply(tmp, function(z) {
    c(mean_sd = mean(z$into_sd_log2, na.rm = TRUE),
      sd_mean = sd(z$into_mean_log2, na.rm = TRUE))
}))
plot(tmp_pos[, 1], tmp_pos[, 2], main = "positive polarity",
     xlab = "average sd of intensity", ylab = "sd of average intensity")
tmp <- split(batch_sum_neg, batch_sum_neg$name)
## Calculate mean across within batch sd and
tmp_neg <- do.call(rbind, lapply(tmp, function(z) {
    c(mean_sd = mean(z$into_sd_log2, na.rm = TRUE),
      sd_mean = sd(z$into_mean_log2, na.rm = TRUE))
}))
plot(tmp_neg[, 1], tmp_neg[, 2], main = "negative polarity",
     xlab = "average sd of intensity", ylab = "sd of average intensity")
```

As expected, also the within-batch variation of intensities is on average much
smaller than the variation of the per-batch average intensities across the
data set. We consider the per-batch average intensities for each internal
standard for the further analyses.

We next evaluate the distribution of intensities for each internal standard
across the whole data set to evaluate whether some standards show a larger
variation than others.

```{r}
tmp <- split(batch_sum_pos$into_mean_log2, batch_sum_pos$name)
tmp <- tmp[order(vapply(tmp, mean, numeric(1)))]
par(mfrow = c(1, 2), mar = c(7, 4, 1, 0.5))
boxplot(tmp, las = 2, main = "positive polarity",
        ylab = expression(log[2]~average~into))
grid(nx = NA, ny = NULL)
tmp <- split(batch_sum_neg$into_mean_log2, batch_sum_neg$name)
tmp <- tmp[order(vapply(tmp, mean, numeric(1)))]
boxplot(tmp, las = 2, main = "negative polarity",
        ylab = expression(log[2]~average~into))
grid(nx = NA, ny = NULL)
```

While for positive polarity all internal standards seem to have about the same
variance, some standards exhibit a larger variance than others for the negative
polarity data. For both polarities it seems however that the variance is
independent of the absolute intensity (i.e. standards with lower intensity don't
have larger variance).

Next we evaluate whether the signal intensities of internal standards per batch
are related to the absolute signals measured for that batch. To this end we
compare the average intensity of internal standards per batch to the average
total ion current per sample as well as to the average intensity of the top 200
peaks with highest abundances per batch. Below we thus identify the 200 peaks
with highest intensities per sample and calculate per batch averages of these
(first taking the average across the 200 peaks per sample and then calculating
the average across all samples per batch).

```{r top200-peaks}
peaks_200_pos <- split(chromPeaks(qc_pos)[, "into"],
                       chromPeaks(qc_pos)[, "sample"]) |>
vapply(function(x) mean(sort(x, decreasing = TRUE)[seq_len(200)]), numeric(1))

peaks_200_neg <- split(chromPeaks(qc_neg)[, "into"],
                       chromPeaks(qc_neg)[, "sample"]) |>
vapply(function(x) mean(sort(x, decreasing = TRUE)[seq_len(200)]), numeric(1))

peaks_200_batch_pos <- vapply(split(peaks_200_pos, qc_pos$batch),
                              mean, numeric(1))
peaks_200_batch_neg <- vapply(split(peaks_200_neg, qc_neg$batch),
                              mean, numeric(1))
```

In addition we calculate the average total ion current of all files within a
batch.

```{r tic-signal, eval = FALSE}
load(paste0(RDATA_PATH, "tic_qc_pos.RData"))
load(paste0(RDATA_PATH, "tic_qc_neg.RData"))


tic_pos <- split(spectrapply(qc_pos, function(z) sum(intensity(z))),
                 qc_pos$batch) |>
vapply(mean, numeric(1))

tic_neg <- split(spectrapply(qc_neg, function(z) sum(intensity(z))),
                 qc_neg$batch) |>
vapply(mean, numeric(1))
```

```{r is-relationship-into-signal-batch, echo = FALSE, fig.path = IMAGE_PATH, fig.width = 8, fig.height = 8, fit.cap = "Relationship between per-batch average intensities of all internal standards and per-batch average intensity of the top 200 chromatographic peaks with highest abundance (left) and per-batch TIC (right)."}
## Average intensity across batch
par(mfrow = c(2, 2))
is_avg_batch_pos <- vapply(split(batch_sum_pos$into_mean_log2,
                                 batch_sum_pos$batch),
                           mean, numeric(1))
plot(is_avg_batch_pos[names(peaks_200_batch_pos)], log2(peaks_200_batch_pos),
     main = "positive polarity",
     xlab = expression(log[2]~signal[top200]),
     ylab = expression(log[2]~signal[internal~standards]),
     bg = paste0(col_batch[names(peaks_200_batch_pos)], 60), pch = 21)
grid()
L <- lm(log2(peaks_200_batch_pos) ~ is_avg_batch_pos)
abline(L)
legend("topleft", legend = paste0("R2 = ", round(summary(L)$r.squared, 2)))
plot(3, 3) # internal standards vs tic
## NEG
is_avg_batch_neg <- vapply(split(batch_sum_neg$into_mean_log2,
                                 batch_sum_neg$batch),
                           mean, numeric(1))
plot(is_avg_batch_neg[names(peaks_200_batch_neg)], log2(peaks_200_batch_neg),
     main = "negative polarity",
     xlab = expression(log[2]~signal[top200]),
     ylab = expression(log[2]~signal[internal~standards]),
     bg = paste0(col_batch[names(peaks_200_batch_neg)], 60), pch = 21)
grid()
L <- lm(log2(peaks_200_batch_neg) ~ is_avg_batch_neg)
abline(L)
legend("topleft", legend = paste0("R2 = ", round(summary(L)$r.squared, 2)))
plot(3, 3) # internal standards vs tic
```

A high correlation, especially for the positive polarity data, can be seen
between the average intensities of internal standards and the average abundance
of the top 200 chromatographic peaks with highest signal. The same relationship
is present also if individual samples are considered (no average across samples
within batch). This suggests that variability of internal standard measurements
between samples and batches is related to the total signal.

```{r is-relationship-into-signal-individual, echo = FALSE, fig.path = IMAGE_PATH, fig.width = 8, fig.height = 8, fit.cap = "Relationship between per-sample average intensities of all internal standards and per-sample average intensity of the top 200 chromatographic peaks with highest abundance."}
is_avg_pos <- split(log2(chromPeaks(is_pos_chr)[, "into"]),
                    chromPeaks(is_pos_chr)[, "sample"]) |>
vapply(mean, numeric(1))
is_avg_neg <- split(log2(chromPeaks(is_neg_chr)[, "into"]),
                    chromPeaks(is_neg_chr)[, "sample"]) |>
vapply(mean, numeric(1))

par(mfrow = c(1, 2))
plot(is_avg_pos, log2(peaks_200_pos),
     main = "positive polarity",
     xlab = expression(log[2]~signal[top200]),
     ylab = expression(log[2]~signal[internal~standards]),
     bg = paste0(col_batch[is_pos_chr$batch], 60), pch = 21)
grid()
L <- lm(log2(peaks_200_pos) ~ is_avg_pos)
abline(L)
legend("topleft", legend = paste0("R2 = ", round(summary(L)$r.squared, 2)))
plot(is_avg_neg, log2(peaks_200_neg),
     main = "negative polarity",
     xlab = expression(log[2]~signal[top200]),
     ylab = expression(log[2]~signal[internal~standards]),
     bg = paste0(col_batch[is_neg_chr$batch], 60), pch = 21)
grid()
L <- lm(log2(peaks_200_neg) ~ is_avg_neg)
abline(L)
legend("topleft", legend = paste0("R2 = ", round(summary(L)$r.squared, 2)))

```

In addition we create and export such correlation plots separately for each
internal standard (not displayed in this report). For some internal standards
a strong relationship is visible while for others it is not.

```{r is-relationship-into-signal-per-analyte}

nms <- union(batch_sum_pos$name, batch_sum_neg$name)
for (is in nms) {
        filename <- paste0(IMAGE_PATH, "IS_", gsub("%", "p", is, fixed = TRUE),
                           "_signal_relationship.png")
        png(file = filename, width = 16, height = 8, units = "cm",
            res = 300, pointsize = 4)
        par(mfrow = c(1, 2))
        keep <- batch_sum_pos$name == is
        if (any(keep)) {
            is_signal <- batch_sum_pos[batch_sum_pos$name == is, ]
            plot(is_signal$into_mean_log2,
                 log2(peaks_200_batch_pos[is_signal$batch]),
                 bg = paste0(col_batch[is_signal$batch], 80), pch = 21,
                 xlab = expression(log[2]~signal[top200]),
                 ylab = expression(log[2]~signal[internal~standards]),
                 main = is)
        } else plot(3, 3, pch = NA,
                    xlab = expression(log[2]~signal[top200]),
                    ylab = expression(log[2]~signal[internal~standards]),
                    main = is)
        grid()
        keep <- batch_sum_neg$name == is
        if (any(keep)) {
            is_signal <- batch_sum_neg[batch_sum_neg$name == is, ]
            plot(is_signal$into_mean_log2,
                 log2(peaks_200_batch_neg[is_signal$batch]),
                 bg = paste0(col_batch[is_signal$batch], 80), pch = 21,
                 xlab = expression(log[2]~signal[top200]),
                 ylab = expression(log[2]~signal[internal~standards]),
                 main = is)
        } else plot(3, 3, pch = NA,
                    xlab = expression(log[2]~signal[top200]),
                    ylab = expression(log[2]~signal[internal~standards]),
                    main = is)
        grid()
        dev.off()
}

```

TODO Summary.
So, summarizing.

- retention time shifts within batches are smaller than between
  batches. Alignment should ideally be performed on the QC samples.
- extent of retention time shifts dependent on the rt: ranges with smaller
  shifts and ranges with larger shifts (ion suppression/column overload?).
- variance in intensity levels seems to be related (for the largest part/most
  internal standards) to absolute signal per sample. Also here, variance within
  batch is lower than between batches.


## Lab standards

We repeat the analysis performed on the internal standards above also on the
*lab standards*, i.e. the set of standards that were used to optimize and
validate the LC-MS setup for the CHRIS study.

```{r std}
std_info <- read.table(
    "https://raw.githubusercontent.com/EuracBiomedicalResearch/lcms-standards/master/data/standards_dilution.txt",
    sep = "\t", header = TRUE, as.is = TRUE)
std_info <- unique(std_info[, -1])
rownames(std_info) <- std_info$name
std_info <- std_info[order(std_info$name), ]
std_info$mzneut = NA
std_info$mz_ion_pos = NA
std_info$mz_ion_neg = NA
for (i in seq(nrow(std_info))) {
    if (grepl("C", std_info$formula[i])){
        std_info$mzneut[i] <- getMolecule(
            as.character(std_info$formula[i]))$exactmass
    } else
        std_info$mzneut[i] = as.numeric(std_info$formula[i])
    ## Calculate also the m/z
    if (!is.na(std_info$POS[i]))
        std_info$mz_ion_pos[i] <-
            mass2mz(std_info$mzneut[i], adduct = std_info$POS[i])[1, 1]
    if (!is.na(std_info$NEG[i]))
        std_info$mz_ion_neg[i] <-
            mass2mz(std_info$mzneut[i], adduct = std_info$NEG[i])[1, 1]
}
std_pos <- std_info[!is.na(std_info$mz_ion_pos), ]
std_pos$mz_ion <- std_pos$mz_ion_pos
std_neg <- std_info[!is.na(std_info$mz_ion_neg), ]
std_neg$mz_ion <- std_neg$mz_ion_neg

std_pos_mz_rt <- chromPeakArea(qc_pos, mz = std_pos$mz_ion, rt = std_pos$RT,
                               diffRt = 40, ppm = 10)
std_pos <- std_pos[!is.na(std_pos_mz_rt[, "mzmin"]), ]
std_pos_mz_rt <- std_pos_mz_rt[!is.na(std_pos_mz_rt[, "mzmin"]), ]
## Extract the ion chromatogram for each standard
rtr <- std_pos_mz_rt[, c("rtmin", "rtmax")]
rtr[, 1] <- rtr[, 1] - 10
rtr[, 2] <- rtr[, 2] + 10
std_pos_chr <- chromatogram(qc_pos, mz = std_pos_mz_rt[, c("mzmin", "mzmax")],
                            rt = rtr, aggregationFun = "max")
fData(std_pos_chr) <- cbind(
    fData(std_pos_chr), std_pos[, c("name", "abbreviation", "formula",
                                    "POS", "RT")])

## Plot the extracted ion chromatograms
dr <- paste0(IMAGE_PATH, "peakdetection/")
dir.create(dr, recursive = TRUE, showWarnings = FALSE)
plot_chrs(std_pos_chr, dr, prefix = "STD-")

## Define the rt window for standards that look ~ OK.
std_pos$rtmin <- NA_real_
std_pos$rtmax <- NA_real_
## Note: just assigning ranges for standards with a signal in MOST samples!
## Those with a * have a quite broad m/z range and should be rechecked
std_pos[std_pos$name == "1-Methylhistidine", c("rtmin", "rtmax")] <- c(165, 210)
std_pos[std_pos$name == "3-Methylhistidine", c("rtmin", "rtmax")] <- c(165, 210)
std_pos[std_pos$name == "8-Oxo-2-Deoxyguanosine", c("rtmin", "rtmax")] <- c(152, 170) # *
std_pos[std_pos$name == "Acetylalanine", c("rtmin", "rtmax")] <- c(33, 45)
std_pos[std_pos$name == "Acetylhistidine", c("rtmin", "rtmax")] <- c(171, 193)
std_pos[std_pos$name == "Adenine", c("rtmin", "rtmax")] <- c(165, 185)
std_pos[std_pos$name == "ADMA", c("rtmin", "rtmax")] <- c(163, 200)
std_pos[std_pos$name == "Alanine", c("rtmin", "rtmax")] <- c(158, 180)
std_pos[std_pos$name == "Arginine", c("rtmin", "rtmax")] <- c(165, 210) # *
std_pos[std_pos$name == "Asparagine", c("rtmin", "rtmax")] <- c(175, 195) # *
std_pos[std_pos$name == "Betaine", c("rtmin", "rtmax")] <- c(150, 170)
std_pos[std_pos$name == "Caffeine", c("rtmin", "rtmax")] <- c(32, 38) # *
std_pos[std_pos$name == "cGMP", c("rtmin", "rtmax")] <- c(170, 190)
std_pos[std_pos$name == "Citrulline", c("rtmin", "rtmax")] <- c(170, 195) # *
std_pos[std_pos$name == "Creatine", c("rtmin", "rtmax")] <- c(158, 185)   # *
std_pos[std_pos$name == "Cystine", c("rtmin", "rtmax")] <- c(200, 220)
std_pos[std_pos$name == "Dimethylglycine", c("rtmin", "rtmax")] <- c(158, 180)
std_pos[std_pos$name == "Galactitol", c("rtmin", "rtmax")] <- c(155, 169) # *
std_pos[std_pos$name == "Glucuronic Acid", c("rtmin", "rtmax")] <- c(140, 153)
std_pos[std_pos$name == "Glutamine", c("rtmin", "rtmax")] <- c(165, 190) # *
std_pos[std_pos$name == "Glyceraldehyde 2-phosphate", c("rtmin", "rtmax")] <- c(180, 230)
std_pos[std_pos$name == "Glycero-phosphocholine", c("rtmin", "rtmax")] <- c(185, 200)
std_pos[std_pos$name == "Glycine", c("rtmin", "rtmax")] <- c(164, 185)
std_pos[std_pos$name == "Guanine", c("rtmin", "rtmax")] <- c(156, 170)
std_pos[std_pos$name == "Guanosine", c("rtmin", "rtmax")] <- c(157, 167)
std_pos[std_pos$name == "Histamine", c("rtmin", "rtmax")] <- c(175, 105)
std_pos[std_pos$name == "Histidine", c("rtmin", "rtmax")] <- c(170, 210)
std_pos[std_pos$name == "Hypoxanthine", c("rtmin", "rtmax")] <- c(100, 135)
std_pos[std_pos$name == "Indoleacetic acid", c("rtmin", "rtmax")] <- c(32, 38)
std_pos[std_pos$name == "Indolelactic acid", c("rtmin", "rtmax")] <- c(32, 45)
std_pos[std_pos$name == "Inosine", c("rtmin", "rtmax")] <- c(135, 155) # *
std_pos[std_pos$name == "L-Aspartic Acid", c("rtmin", "rtmax")] <- c(172, 195)
std_pos[std_pos$name == "L-Glutamic Acid", c("rtmin", "rtmax")] <- c(165, 185)
std_pos[std_pos$name == "Lysine", c("rtmin", "rtmax")] <- c(173, 210)
std_pos[std_pos$name == "Methionine", c("rtmin", "rtmax")] <- c(155, 169)
std_pos[std_pos$name == "Methioninesulfoxide", c("rtmin", "rtmax")] <- c(174, 190)
std_pos[std_pos$name == "Ornithine", c("rtmin", "rtmax")] <- c(175, 210)
std_pos[std_pos$name == "Proline", c("rtmin", "rtmax")] <- c(158, 180)
std_pos[std_pos$name == "Putrescine", c("rtmin", "rtmax")] <- c(150, 173)
std_pos[std_pos$name == "Salicylic acid", c("rtmin", "rtmax")] <- c(26, 36)
std_pos[std_pos$name == "SDMA", c("rtmin", "rtmax")] <- c(163, 200)
std_pos[std_pos$name == "Serine", c("rtmin", "rtmax")] <- c(171, 192)
std_pos[std_pos$name == "Sphingosine", c("rtmin", "rtmax")] <- c(28, 38)
std_pos[std_pos$name == "Taurine", c("rtmin", "rtmax")] <- c(165, 178)
std_pos[std_pos$name == "Threonine", c("rtmin", "rtmax")] <- c(164, 185)
std_pos[std_pos$name == "Tryptophan", c("rtmin", "rtmax")] <- c(140, 160) # *
std_pos[std_pos$name == "Tyrosine", c("rtmin", "rtmax")] <- c(159, 175)
std_pos[std_pos$name == "Valine", c("rtmin", "rtmax")] <- c(150, 172)
std_pos <- std_pos[!is.na(std_pos$rtmin), ]

## Re-create the m/z and rt range matrix
std_pos_mz_rt <- chromPeakArea(qc_pos, mz = std_pos$mz_ion, ppm = 10,
                               rt = std_pos[, c("rtmin", "rtmax")], diffRt = 0)
std_pos_chr <- chromatogram(qc_pos, mz = std_pos_mz_rt[, c("mzmin", "mzmax")],
                            rt = std_pos_mz_rt[, c("rtmin", "rtmax")],
                            aggregationFun = "max", include = "apex_within")
fData(std_pos_chr) <- cbind(
    fData(std_pos_chr), std_pos[, c("name", "abbreviation", "formula",
                                    "POS", "RT")])

## For each IS (row), in each sample (column) keep the chrom peak with the
## highest peak.
std_pos_chr <- filterChromPeaks(std_pos_chr, method = "keepTop",
                                n = 1L, order = "maxo")
save(std_pos_chr, file = paste0(RDATA_PATH, "std_pos_chr.RData"))
## Plot again.
plot_chrs(std_pos_chr, dr, prefix = "STD-",
          suffix = "-refined_POS", rect = FALSE)


## --------------------------------------------------
## Repeat for negative polarity.
is_neg_mz_rt <- chromPeakArea(qc_neg, mz = is_neg$mz_ion, rt = is_neg$RT,
                              diffRt = 40, ppm = 10)
is_neg <- is_neg[!is.na(is_neg_mz_rt[, "mzmin"]), ]
is_neg_mz_rt <- is_neg_mz_rt[!is.na(is_neg_mz_rt[, "mzmin"]), ]
## Extract the ion chromatogram for each standard
rtr <- is_neg_mz_rt[, c("rtmin", "rtmax")]
rtr[, 1] <- rtr[, 1] - 10
rtr[, 2] <- rtr[, 2] + 10
is_neg_chr <- chromatogram(qc_neg, mz = is_neg_mz_rt[, c("mzmin", "mzmax")],
                           rt = rtr, aggregationFun = "max")
fData(is_neg_chr) <- cbind(
    fData(is_neg_chr), is_neg[, c("name", "abbreviation", "formula",
                                  "NEG", "RT")])

## Plot the extracted ion chromatograms
plot_chrs(is_neg_chr, dr, prefix = "IS-", suffix = "_NEG")

## Define the rt window for standards that look ~ OK.
is_neg$rtmin <- NA_real_
is_neg$rtmax <- NA_real_
is_neg$rtmin[is_neg$name == "Creatinine (N-Methhyl_D3)"] <- 120
is_neg$rtmax[is_neg$name == "Creatinine (N-Methhyl_D3)"] <- 145
is_neg$rtmin[is_neg$name == "L-Alanine (13C3, 99%; 15N, 99%)"] <- 175
is_neg$rtmax[is_neg$name == "L-Alanine (13C3, 99%; 15N, 99%)"] <- 187
is_neg$rtmin[is_neg$name == "L-Arginine HCl (13C6, 99%; 15N4, 99%)"] <- 175
is_neg$rtmax[is_neg$name == "L-Arginine HCl (13C6, 99%; 15N4, 99%)"] <- 200
is_neg$rtmin[is_neg$name == "L-Aspartic acid (13C4, 99%; 15N, 99%)"] <- 170
is_neg$rtmax[is_neg$name == "L-Aspartic acid (13C4, 99%; 15N, 99%)"] <- 200
is_neg$rtmin[is_neg$name == "L-Cystine (13C6, 99%; 15N2, 99%)"] <- 202
is_neg$rtmax[is_neg$name == "L-Cystine (13C6, 99%; 15N2, 99%)"] <- 230
is_neg$rtmin[is_neg$name == "L-Glutamic acid (13C5, 99%; 15N, 99%)"] <- 165
is_neg$rtmax[is_neg$name == "L-Glutamic acid (13C5, 99%; 15N, 99%)"] <- 200
is_neg$rtmin[is_neg$name == "L-Histidine HCl H2O (13C6; 15N3, 99%)"] <- 175
is_neg$rtmax[is_neg$name == "L-Histidine HCl H2O (13C6; 15N3, 99%)"] <- 205
is_neg$rtmin[is_neg$name == "L-Lysine 2HCl (13C6, 99%; 15N2, 99%)"] <- 175
is_neg$rtmax[is_neg$name == "L-Lysine 2HCl (13C6, 99%; 15N2, 99%)"] <- 205
is_neg$rtmin[is_neg$name == "L-Methionine (13C5, 99%; 15N, 99%)"] <- 155
is_neg$rtmax[is_neg$name == "L-Methionine (13C5, 99%; 15N, 99%)"] <- 170
is_neg$rtmin[is_neg$name == "L-Phenylalanine (13C9, 99%; 15N, 99%)"] <- 140
is_neg$rtmax[is_neg$name == "L-Phenylalanine (13C9, 99%; 15N, 99%)"] <- 161
is_neg$rtmin[is_neg$name == "L-Serine (13C3, 99%; 15N, 99%)"] <- 172
is_neg$rtmax[is_neg$name == "L-Serine (13C3, 99%; 15N, 99%)"] <- 190
is_neg$rtmin[is_neg$name == "L-Threonine (13C4, 99%; 15N, 99%)"] <- 163
is_neg$rtmax[is_neg$name == "L-Threonine (13C4, 99%; 15N, 99%)"] <- 183
is_neg$rtmin[is_neg$name == "L-Tyrosine (13C9, 99%; 15N, 99%)"] <- 160
is_neg$rtmax[is_neg$name == "L-Tyrosine (13C9, 99%; 15N, 99%)"] <- 175
is_neg$rtmin[is_neg$name == "L-Valine (13C5, 99%; 15N, 99%)"] <- 160
is_neg$rtmax[is_neg$name == "L-Valine (13C5, 99%; 15N, 99%)"] <- 175
is_neg <- is_neg[!is.na(is_neg$rtmin), ]

## Re-create the m/z and rt range matrix
is_neg_mz_rt <- chromPeakArea(qc_neg, mz = is_neg$mz_ion, ppm = 10,
                              rt = is_neg[, c("rtmin", "rtmax")], diffRt = 0)
is_neg_chr <- chromatogram(
    qc_neg, mz = is_neg_mz_rt[, c("mzmin", "mzmax")],
    rt = is_neg_mz_rt[, c("rtmin", "rtmax")],
    aggregationFun = "max", include = "apex_within")
fData(is_neg_chr) <- cbind(
    fData(is_neg_chr), is_neg[, c("name", "abbreviation", "formula",
                                  "NEG", "RT")])

## For each IS (row), in each sample (column) keep the chrom peak with the
## highest peak.
is_neg_chr <- filterChromPeaks(is_neg_chr, method = "keepTop",
                               n = 1L, order = "maxo")
save(is_neg_chr, file = paste0(RDATA_PATH, "is_neg_chr.RData"))
## Plot again.
plot_chrs(is_neg_chr, dr, prefix = "IS-", suffix = "-refined_NEG", rect = FALSE)

```


## Evaluating the *dilution QC samples*

- check how the QC dilution look.

# Per-batch quality assessment

In this section we evaluate the signal within each batch, comparing the data
from study samples against those of the QC samples.

# Session information

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