peak_detection.Rmd

---
title: "Chromatographic peak detection of the CHRIS untargeted metabolomics data"
author: "Mar Garcia-Aloy, 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.Rmd")$mtime`<br />
**Compiled**: `r date()`

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

## Define paths:
filename <- "peak_detection"
## Path to save the images; remove all old images.
IMAGE_PATH <- paste0("images/", filename, "/")
dir.create(IMAGE_PATH, recursive = TRUE, showWarnings = FALSE)
## Path to store RData files
RDATA_PATH <- paste0("data/RData/", filename, "/")
dir.create(RDATA_PATH, recursive = TRUE, showWarnings = FALSE)

exclusion <- c("041016_BLANK_POS_1",
               "041016_BLANK_POS_2-041016_BLANK_POS_1",
               "041016_BLANK_POS_2",
               "041016_BLANK_POS_1 (3)",
               "041016_BLANK_POS_1 (4)",
               "06102032_BLANK_POS_b",
               "251016_BLANK_POS_1",
               "251016_BLANK_POS_2",
               "20161025_BLANK_POS_a",
               "20161025_BLANK_POS_b",
               "211116_BLANK_POS_1",
               "211116_BLANK_POS_2",
               "051216_BLANK_POS_1",
               "051216_BLANK_POS_2",
               "20161205_BLANK_POS_a",
               "20161205_BLANK_POS_b",
               "20161205_BLANK_POS_a_R",
               "20161205_BLANK_POS_a_R (2)",
               "20161205_BLANK_POS_b_R",
               "20161205_BLANK_POS_a_R1",
               "090117_BLANK_POS_1",
               "090117_BLANK_POS_2",
               "20170112_BLANK_POS_a",
               "20170112_BLANK_POS_b",
               "120117_BLANK_POS_1",
               "20170209_BLANK_POS_a",
               "20170406_BLANK_POS_a",
               "20170731_BLANK_POS (2)",

               "06102016_0010060146_POS",
               "20170615_0010282264_POS_R",

               "20161116_POOL_POS_2",
               "20170405_0010007645_POS",

               "20170209_BLANK_NEG_a",
               "20170209_BLANK_NEG_b",
               "20170220_BLANK_NEG_a",
               "20170220_BLANK_NEG_a (2)",
               "20170224_BLANK_NEG_a",
               "20170407_BLANK_NEG_a",
               "20170731_BLANK_NEG",
               "20170109_BLANK_NEG_a",
               "20170111_BLANK_NEG_a",
               "20170206_BLANK_NEG_a",
               "20170207_BLANK_NEG_a",

               "20161024_0010302963_NEG",
               "20161201_POOL_NEG_1.mzML",
               "20161201_POOL_NEG_2.mzML")

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

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 the peak detection of the untargeted metabolomics
data of the CHRIS population study, both for positive and negative polarity.

We first load and define the phenodata table with all samples/mzML files which
are part of the present data set.

```{r libraries, message = FALSE}
library(readxl)
library(RColorBrewer)
library(xcms)

injections <- as.data.frame(read_xlsx("data/chris-files-annotated.xlsx"))
injections$id <- gsub(".mzML", "", injections$mzML_file)
injections <- injections[order(injections$timestamp), ]
injections$batch <- injections$run_id
injections <- injections[!injections$sample_name %in% exclusion, ]
injections <- injections[!injections$type %in% c("ACN_IS", "EQ_MIX"), ]

## Define the injection index
runids <- split(injections$run_id,
                f = factor(injections$run_id,
                           levels = unique(injections$run_id)))
injections$inj_idx <- unlist(lapply(lengths(runids), seq), use.names = FALSE)

## split by polarity
pd_pos <- injections[injections$polarity == "POS", ]
pd_neg <- injections[injections$polarity == "NEG", ]
```

Performing first the analysis for the positive polarity data. Negative polarity
data will be processed afterwards using the same settings. Below we load the
positive polarity data, subset to the defined retention time range and remove
unnecessary spectrum metadata.

```{r pos-load, eval = !file.exists(paste0(RDATA_PATH, "data_pos.RData"))}
data_pos <- readMSData(files = paste0(MZML_PATH, pd_pos$mzML_file),
                       pdata = as(AnnotatedDataFrame(pd_pos),
                                  "NAnnotatedDataFrame"),
                       mode = "onDisk")
data_pos <- filterRt(data_pos, c(0, rt_cut))
data_pos <- selectFeatureData(
    data_pos, fcol = c("fileIdx", "spIdx", "acquisitionNum", "retentionTime",
                       "msLevel", "precursorScanNum", "polarity", "centroided",
                       "seqNum")
)
```

Before proceeding to the peak detection we extract base peak and total ion
chromatograms.

```{r bpc-pos, eval = !file.exists(paste0(RDATA_PATH, "bpc_pos.RData"))}
bpc_pos <- chromatogram(data_pos, aggregationFun = "max")
save(bpc_pos, file = paste0(RDATA_PATH, "bpc_pos.RData"))
rm(bpc_pos)
gc()
tic_pos <- chromatogram(data_pos, aggregationFun = "sum")
save(tic_pos, file = paste0(RDATA_PATH, "tic_pos.RData"))
rm(tic_pos)
gc()
```


# Peak detection and refinement

We perform a centWave-based peak detection followed by a *peak refinement* step
that aims at reducing the peak detection artifacts, such as split peaks or
overlapping peaks. Below we define the parameters for the peak detection and
subsequently perform the analysis.

```{r pos-peak-detection, warning = FALSE, message = FALSE, eval = !file.exists(paste0(RDATA_PATH, "data_pos.RData"))}
cwp <- CentWaveParam(
    peakwidth = c(2, 20),
    ppm = 50,
    snthresh = 5,
    mzdiff = 0.001,
    prefilter = c(3, 500),
    noise = 100,
    integrate = 2)

register(MulticoreParam(ncores - 1L))
tm <- system.time(
    data_pos <- findChromPeaks(data_pos, param = cwp)
)
save(data_pos, file = paste0(RDATA_PATH, "data_pos.RData"))
```

```{r pos-peak-detection-load, echo = FALSE, warning = FALSE, message = FALSE}
load(paste0(RDATA_PATH, "data_pos.RData"))
if (!any(ls() == "tm"))
    tm <- c(elapsed = 1)
```

Next we refine the identified chromatographic peaks. Processing took
`r tm["elapsed"] / 60` minutes.

```{r pos-peak-refinement, warning = FALSE, message = FALSE, eval = !file.exists(paste0(RDATA_PATH, "data_pos_ref.RData"))}
mnp <- MergeNeighboringPeaksParam(
    expandRt = 2,
    expandMz = 0.001,
    ppm = 10,
    minProp = 0.66)

register(MulticoreParam(workers = ceiling(ncores / 2),
                        timeout = 86400L))
tm <- system.time(
    data_pos <- refineChromPeaks(data_pos, param = mnp, BPPARAM = bpparam())
)
save(data_pos, file = paste0(RDATA_PATH, "data_pos_ref.RData"))
```

```{r pos-peak-refinement-load, echo = FALSE, warning = FALSE, message = FALSE}
load(paste0(RDATA_PATH, "data_pos_ref.RData"))
if (!any(ls() == "tm"))
    tm <- c(elapsed = 1)
```

Processing took `r tm["elapsed"] / 60` minutes.

```{r neg-polarity, echo = FALSE}
if (!file.exists(paste0(RDATA_PATH, "data_neg.RData"))) {
    data_neg <- readMSData(files = paste0(MZML_PATH, pd_neg$mzML_file),
                           pdata = as(AnnotatedDataFrame(pd_neg),
                                      "NAnnotatedDataFrame"),
                           mode = "onDisk")
    data_neg <- filterRt(data_neg, c(0, rt_cut))
    data_neg <- selectFeatureData(
        data_neg, fcol = c("fileIdx", "spIdx", "acquisitionNum",
                           "retentionTime", "msLevel", "precursorScanNum",
                           "polarity", "centroided", "seqNum")
    )
    register(MulticoreParam(ncores - 1L))
    if (!file.exists(paste0(RDATA_PATH, "bpc_neg.RData"))) {
        bpc_neg <- chromatogram(data_neg, aggregationFun = "max")
        save(bpc_neg, file = paste0(RDATA_PATH, "bpc_neg.RData"))
        rm(bpc_neg)
        gc()
    }
    if (!file.exists(paste0(RDATA_PATH, "tic_neg.RData"))) {
        tic_neg <- chromatogram(data_neg, aggregationFun = "sum")
        save(tic_neg, file = paste0(RDATA_PATH, "tic_neg.RData"))
        rm(tic_neg)
        gc()
    }
    cwp <- CentWaveParam(
        peakwidth = c(2, 20),
        ppm = 50,
        snthresh = 5,
        mzdiff = 0.001,
        prefilter = c(3, 500),
        noise = 100,
        integrate = 2)
    system.time(
        data_neg <- findChromPeaks(data_neg, param = cwp)
    )
    save(data_neg, file = paste0(RDATA_PATH, "data_neg.RData"))
} else load(paste0(RDATA_PATH, "data_neg.RData"))

if (!file.exists(paste0(RDATA_PATH, "data_neg_ref.RData"))) {
    mnp <- MergeNeighboringPeaksParam(
        expandRt = 2,
        expandMz = 0.001,
        ppm = 10,
        minProp = 0.66)
    ## Peak refinement
    register(MulticoreParam(workers = ceiling(ncores / 2),
                            timeout = 86400L))
    tm <- system.time(
        data_neg <- refineChromPeaks(data_neg, param = mnp, BPPARAM = bpparam())
    )
    save(data_neg, file = paste0(RDATA_PATH, "data_neg_ref.RData"))
} else {
    load(paste0(RDATA_PATH, "data_neg_ref.RData"))
}
```


# Initial quality assessment

In this section we perform an initial quality assessment comparing the base peak
chromatograms of QC samples and the number of identified chromatographic peaks
across the whole data set. To this end we first extract the base peak
chromatograms for the positive polarity data (the same code is used to
create the BPC for negative polarity.

```{r bpc-pos-load, echo = FALSE}
load(paste0(RDATA_PATH, "bpc_pos.RData"))
```

Next we plot the base peak chromatogram for all QC samples in the data set.

```{r pos-bpc-qc, fig.path = IMAGE_PATH, caption = "BPC of QC samples. Positive polarity", fig.width = 16, fig.height = 8, echo = FALSE}
## pdf("pos-bpc-qc.pdf", width = 16, height = 8)
plot(bpc_pos[, bpc_pos$type == "QC"], main = "QC samples, positive polarity",
     col = "#00000020")
## dev.off()
```

Because of the large number of samples being plotted it is hard to see
differences between the BPC. As expected from such a large data set, retention
time shifts seem to be present. Next we compare the distribution of base peak
intensities of all QC samples in the data set.

```{r pos-bpc-qc-boxplot, echo = FALSE, fig.cap = "Distribution of base peak intensities of QC samples.", fig.width = 16, fig.height = 8}
## pdf("pos-bpc-qc-boxplot.pdf", width = 16, height = 8)
tmp <- bpc_pos[, bpc_pos$type == "QC"]
par(mar = c(5, 4.5, 1.5, 0.5))
boxplot(lapply(tmp, function(z) log2(intensity(z))), outline = FALSE,
        ylab = expression(log[2]~base~peak~intensity), xaxt = "n",
        main = "QC samples, positive polarity")
batch_idx <- as.integer(factor(tmp$batch, levels = unique(tmp$batch)))
batch_tick <- c(1, which(diff(batch_idx) == 1), ncol(tmp))
axis(side = 1, at = batch_tick, labels = rep("", length(batch_tick)))
axis(side = 1, at = batch_tick[-length(batch_tick)] + diff(batch_tick) / 2,
     labels = unique(tmp$batch), las = 2, tick = 0, line = -0.8)
## dev.off()

```

We also create base peak chromatograms separately for each plate (code and plots
not shown here).

```{r pos-bpc-plate, echo = FALSE}
col_type <- brewer.pal(7, "Set1")[c(2, 3, 5, 4)] # blue, green, orange, purple
names(col_type) <- c("study", "QC", "BLANK", "QC_DIL")
batches <- unique(bpc_pos$batch)
dr <- paste0(IMAGE_PATH, "pos-per-batch/")
dir.create(dr, showWarnings = FALSE, recursive = TRUE)
for (batch in batches) {
    tmp <- bpc_pos[, bpc_pos$batch == batch]
    png(paste0(dr, "BPC-chromatogram-", batch, ".png"), width = 16, height = 8,
        units = "cm", res = 300, pointsize = 5)
    plot(tmp, col = paste0(col_type[as.character(tmp$type)], 60),
         main = paste0(batch, ": ", tmp$year[1], "-",
                       tmp$month[1], "-", tmp$day[1]))
    dev.off()
    png(paste0(dr, "BPC-boxplot-", batch, ".png"), width = 16, height = 8,
        units = "cm", res = 300, pointsize = 5)
    boxplot(lapply(tmp, function(z) log2(intensity(z))), outline = FALSE,
            col = paste0(col_type[as.character(tmp$type)], 60),
            ylab = expression(log[2]~base~peak~intensity), las = 2,
            main = paste0(batch, ": ", tmp$year[1], "-",
                          tmp$month[1], "-", tmp$day[1]))
    dev.off()
}
```

```{r pos-sirocco, echo = FALSE, eval = FALSE}
png("BPLT00000063-BPC-mirror.png", width = 12, height = 12, units = "cm",
    res = 300, pointsize = 5)
a <- bpc_pos[, bpc_pos$batch == "2016_09" & bpc_pos$type == "QC"]
b <- bpc_pos[, bpc_pos$batch == "2016_21" & bpc_pos$type == "QC"]
plot(xlim = range(unlist(lapply(a, rtime))), pch = NA, 3, 3,
     ylim = max(unlist(lapply(a, intensity))) * c(-1, 1), xlab = "rtime",
     ylab = "base peak intensity")
abline(0, 1)
grid()
rem <- lapply(a, function(z) points(rtime(z), intensity(z), type = "l",
                                    col = paste0(col_type["QC"], 60)))
rem <- lapply(b, function(z) points(rtime(z), -intensity(z), type = "l",
                                    col = paste0(col_type["QC"], 60)))
text(x = 100, y = 2e6, label = paste0(a$batch[1], ", ", a$plate_name[1]),
     pos = 4)
text(x = 100, y = -2e6, label = paste0(b$batch[1], ", ", b$plate_name[1]),
     pos = 4)
dev.off()

png("BPLT00000067-BPC-mirror.png", width = 12, height = 12, units = "cm",
    res = 300, pointsize = 5)
a <- bpc_pos[, bpc_pos$batch == "2016_10" & bpc_pos$type == "QC"]
b <- bpc_pos[, bpc_pos$batch == "2016_22" & bpc_pos$type == "QC"]
plot(xlim = range(unlist(lapply(a, rtime))), pch = NA, 3, 3,
     ylim = max(unlist(lapply(a, intensity))) * c(-1, 1), xlab = "rtime",
     ylab = "base peak intensity")
abline(0, 1)
grid()
rem <- lapply(a, function(z) points(rtime(z), intensity(z), type = "l",
                                    col = paste0(col_type["QC"], 60)))
rem <- lapply(b, function(z) points(rtime(z), -intensity(z), type = "l",
                                    col = paste0(col_type["QC"], 60)))
text(x = 100, y = 2e6, label = paste0(a$batch[1], ", ", a$plate_name[1]),
     pos = 4)
text(x = 100, y = -2e6, label = paste0(b$batch[1], ", ", b$plate_name[1]),
     pos = 4)
dev.off()

```


```{r pos-peaks-per-file, fig.path = IMAGE_PATH, caption = "Positive polarity data: peak count and signal distribution per sample.", fig.width = 16, fig.height = 10, echo = FALSE}

intos <- split(log2(chromPeaks(data_pos)[, "into"]),
               f = as.integer(chromPeaks(data_pos)[, "sample"]))
intos_mean <- lapply(intos, mean)

## pdf("peak_summary.pdf", width = 12, height = 8)

par(mfrow = c(2, 1), mar = c(0.5, 4.5, 0.5, 0.5))
barplot(lengths(intos), border = paste0(col_type[data_pos$type], 40),
        ylab = "peak count", xaxt = "n", main = "positive polarity")
boxplot(intos, outline = FALSE, medlty = 0, whisklty = 0, staplelty = 0,
        border = paste0(col_type[data_pos$type], 40), ylim = c(10, 15),
        xaxt = "n", ylab = expression(log[2]~intensity))
points(x = seq_along(intos_mean), intos_mean,
       col = paste0(col_type[data_pos$type], 80))

year_month <- paste(data_pos$year, data_pos$month, sep = "-")
is <- table(year_month)
year_month_x <- c(1, cumsum(is))
for (i in 2:length(year_month_x)) {
    arrows(x0 = year_month_x[i - 1], x1 = year_month_x[i],
           y0 = 14.5, y1 = 14.5, code = 3, angle = 90, length = 0.1)
    text(x = mean(year_month_x[c(i - 1, i)]), y = 14.9,
         label = names(year_month_x)[i])
}

## dev.off()

```

While the number of detected peaks is comparable within batch (plate) a large
differences can be observed between some batches. A comparably low number of
peaks was identified in the first 4 batches, while the largest number of
detected peaks are present in the last batches. In contrast to the number of
peaks, the intensity of the detected peaks is comparable across all batches.

- compare BPC of plates with and without phospholipid removal.
- create EICs of standards.

Next we create the same plots for negative polarity.

```{r bpc-neg-load, echo = FALSE}
load(paste0(RDATA_PATH, "bpc_neg.RData"))
```

Next we plot the base peak chromatogram for all QC samples in the data set.

```{r neg-bpc-qc, fig.path = IMAGE_PATH, caption = "BPC of QC samples. Negative polarity", fig.width = 16, fig.height = 8, echo = FALSE}
## pdf("neg-bpc-qc.pdf", width = 16, height = 8)
plot(bpc_neg[, bpc_neg$type == "QC"], main = "QC samples, negative polarity",
     col = "#00000020")
## dev.off()
```

Because of the large number of samples being plotted it is hard to see
differences between the BPC. As expected from such a large data set, retention
time shifts seem to be present. Next we compare the distribution of base peak
intensities of all QC samples in the data set.

```{r neg-bpc-qc-boxplot, echo = FALSE, fig.cap = "Distribution of base peak intensities of QC samples.", fig.width = 16, fig.height = 8}
## pdf("neg-bpc-qc-boxplot.pdf", width = 16, height = 8)
tmp <- bpc_neg[, bpc_neg$type == "QC"]
par(mar = c(5, 4.5, 1.5, 0.5))
boxplot(lapply(tmp, function(z) log2(intensity(z))), outline = FALSE,
        ylab = expression(log[2]~base~peak~intensity), xaxt = "n",
        main = "QC samples, negative polarity")
batch_idx <- as.integer(factor(tmp$batch, levels = unique(tmp$batch)))
batch_tick <- c(1, which(diff(batch_idx) == 1), ncol(tmp))
axis(side = 1, at = batch_tick, labels = rep("", length(batch_tick)))
axis(side = 1, at = batch_tick[-length(batch_tick)] + diff(batch_tick) / 2,
     labels = unique(tmp$batch), las = 2, tick = 0, line = -0.8)
## dev.off()

```

We also create base peak chromatograms separately for each plate (code and plots
not shown here).

```{r neg-bpc-plate, echo = FALSE}
batches <- unique(bpc_neg$batch)
dr <- paste0(IMAGE_PATH, "neg-per-batch/")
dir.create(dr, showWarnings = FALSE, recursive = TRUE)
for (batch in batches) {
    tmp <- bpc_neg[, bpc_neg$batch == batch]
    png(paste0(dr, "BPC-chromatogram-", batch, ".png"), width = 16, height = 8,
        units = "cm", res = 300, pointsize = 5)
    plot(tmp, col = paste0(col_type[as.character(tmp$type)], 60),
         main = paste0(batch, ": ", tmp$year[1], "-",
                       tmp$month[1], "-", tmp$day[1]))
    dev.off()
    png(paste0(dr, "BPC-boxplot-", batch, ".png"), width = 16, height = 8,
        units = "cm", res = 300, pointsize = 5)
    boxplot(lapply(tmp, function(z) log2(intensity(z))), outline = FALSE,
            col = paste0(col_type[as.character(tmp$type)], 60),
            ylab = expression(log[2]~base~peak~intensity), las = 2,
            main = paste0(batch, ": ", tmp$year[1], "-",
                          tmp$month[1], "-", tmp$day[1]))
    dev.off()
}
```

```{r neg-peaks-per-file, fig.path = IMAGE_PATH, caption = "Negative polarity data: peak count and signal distribution per sample.", fig.width = 16, fig.height = 10, echo = FALSE}


intos <- split(log2(chromPeaks(data_neg)[, "into"]),
               f = as.integer(chromPeaks(data_neg)[, "sample"]))
intos_mean <- lapply(intos, mean)

## pdf("peak_summary.pdf", width = 12, height = 8)

par(mfrow = c(2, 1), mar = c(0.5, 4.5, 0.5, 0.5))
barplot(lengths(intos), border = paste0(col_type[data_neg$type], 40),
        ylab = "peak count", xaxt = "n", main = "negative polarity")
boxplot(intos, outline = FALSE, medlty = 0, whisklty = 0, staplelty = 0,
        border = paste0(col_type[data_neg$type], 40), ylim = c(10, 15),
        xaxt = "n", ylab = expression(log[2]~intensity))
points(x = seq_along(intos_mean), intos_mean,
       col = paste0(col_type[data_neg$type], 80))

year_month <- paste(data_neg$year, data_neg$month, sep = "-")
is <- table(year_month)
year_month_x <- c(1, cumsum(is))
for (i in 2:length(year_month_x)) {
    arrows(x0 = year_month_x[i - 1], x1 = year_month_x[i],
           y0 = 14.5, y1 = 14.5, code = 3, angle = 90, length = 0.1)
    text(x = mean(year_month_x[c(i - 1, i)]), y = 14.9,
         label = names(year_month_x)[i])
}

## dev.off()

```


# Session information

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