Package: xcms
Authors: Johannes Rainer
Modified: 2025-02-14 13:40:33.337701
Compiled: Mon Mar 17 20:14:09 2025

1 Introduction

The xcms package provides the functionality to perform the preprocessing of LC-MS, GC-MS or LC-MS/MS data in which raw signals from mzML, mzXML or CDF files are processed into feature abundances. This preprocessing includes chromatographic peak detection, sample alignment and correspondence analysis.

The first version of the package was already published in 2006 [1] and has since been updated and modernized in several rounds to better integrate it with other R-based packages for the analysis of untargeted metabolomics data. This includes version 3 of xcms that used the MSnbase package for MS data representation [2]. The most recent update (xcms version 4) enables in addition preprocessing of MS data represented by the modern MsExperiment and Spectra packages which provides an even better integration with the RforMassSpectrometry R package ecosystem simplifying e.g. also compound annotation [3].

This document describes data import, exploration and preprocessing of a simple test LC-MS data set with the xcms package version >= 4. The same functions can be applied to the older MSnbase-based workflows (xcms version 3). Additional documents and tutorials covering also other topics of untargeted metabolomics analysis are listed at the end of this document. There is also a xcms tutorial available with more examples and details. To get a complete overview of LCMS-MS analysis, an end-to-end workflow Metabonaut website, which integrate the xcms preprocessing steps with the downstream analysis, is available.

2 Preprocessing of LC-MS data

2.1 Data import

xcms supports analysis of any LC-MS(/MS) data that can be imported with the Spectra package. Such data will typically be provided in (AIA/ANDI) NetCDF, mzXML and mzML format but can, through dedicated extensions to the Spectra package, also be imported from other sources, e.g. also directly from raw data files in manufacturer’s formats.

For demonstration purpose we will analyze in this document a small subset of the data from [4] in which the metabolic consequences of the knock-out of the fatty acid amide hydrolase (FAAH) gene in mice was investigated. The raw data files (in NetCDF format) are provided through the faahKO data package. The data set consists of samples from the spinal cords of 6 knock-out and 6 wild-type mice. Each file contains data in centroid mode acquired in positive ion polarity from 200-600 m/z and 2500-4500 seconds. To speed-up processing of this vignette we will restrict the analysis to only 8 files.

Below we load all required packages, locate the raw CDF files within the faahKO package and build a phenodata data.frame describing the experimental setup. Generally, such data frames should contain all relevant experimental variables and sample descriptions (including also the names of the raw data files) and will be imported into R using either the read.table() function (if the file is in csv or tabulator delimited text file format) or also using functions from the readxl R package if it is in Excel file format.

library(xcms)
library(faahKO)
library(RColorBrewer)
library(pander)
library(pheatmap)
library(MsExperiment)

## Get the full path to the CDF files
cdfs <- dir(system.file("cdf", package = "faahKO"), full.names = TRUE,
            recursive = TRUE)[c(1, 2, 5, 6, 7, 8, 11, 12)]
## Create a phenodata data.frame
pd <- data.frame(sample_name = sub(basename(cdfs), pattern = ".CDF",
                                   replacement = "", fixed = TRUE),
                 sample_group = c(rep("KO", 4), rep("WT", 4)),
                 stringsAsFactors = FALSE)

We next load our data using the readMsExperiment function from the MsExperiment package.

faahko <- readMsExperiment(spectraFiles = cdfs, sampleData = pd)
faahko
## Object of class MsExperiment 
##  Spectra: MS1 (10224) 
##  Experiment data: 8 sample(s)
##  Sample data links:
##   - spectra: 8 sample(s) to 10224 element(s).

The MS spectra data from our experiment is now available as a Spectra object within faahko. Note that this MsExperiment container could in addition to spectra data also contain other types of data or also references to other files. See the vignette from the MsExperiment for more details. Also, when loading data from mzML, mzXML or CDF files, by default only general spectra data is loaded into memory while the actual peaks data, i.e. the m/z and intensity values are only retrieved on-the-fly from the raw files when needed (this is similar to the MSnbase on-disk mode described in [2]). This guarantees a low memory footprint hence allowing to analyze also large experiments without the need of high performance computing environments. Note that also different alternative backends (and hence data representations) could be used for the Spectra object within faahko with eventually even lower memory footprint, or higher performance. See the package vignette from the Spectra package or the SpectraTutorials tutorial for more details on Spectra backends and how to change between them.

2.2 Initial data inspection

The MsExperiment object is a simple and flexible container for MS experiments. The raw MS data is stored as a Spectra object that can be accessed through the spectra() function.

spectra(faahko)
## MSn data (Spectra) with 10224 spectra in a MsBackendMzR backend:
##         msLevel     rtime scanIndex
##       <integer> <numeric> <integer>
## 1             1   2501.38         1
## 2             1   2502.94         2
## 3             1   2504.51         3
## 4             1   2506.07         4
## 5             1   2507.64         5
## ...         ...       ...       ...
## 10220         1   4493.56      1274
## 10221         1   4495.13      1275
## 10222         1   4496.69      1276
## 10223         1   4498.26      1277
## 10224         1   4499.82      1278
##  ... 33 more variables/columns.
## 
## file(s):
## ko15.CDF
## ko16.CDF
## ko21.CDF
##  ... 5 more files

All spectra are organized sequentially (i.e., not by file) but the fromFile() function can be used to get for each spectrum the information to which of the data files it belongs. Below we simply count the number of spectra per file.

table(fromFile(faahko))
## 
##    1    2    3    4    5    6    7    8 
## 1278 1278 1278 1278 1278 1278 1278 1278

Information on samples can be retrieved through the sampleData() function.

sampleData(faahko)
## DataFrame with 8 rows and 3 columns
##   sample_name sample_group spectraOrigin
##   <character>  <character>   <character>
## 1        ko15           KO /home/bioc...
## 2        ko16           KO /home/bioc...
## 3        ko21           KO /home/bioc...
## 4        ko22           KO /home/bioc...
## 5        wt15           WT /home/bioc...
## 6        wt16           WT /home/bioc...
## 7        wt21           WT /home/bioc...
## 8        wt22           WT /home/bioc...

Each row in this DataFrame represents one sample (input file). Using [ it is possible to subset a MsExperiment object by sample. Below we subset the faahko to the 3rd sample (file) and access its spectra and sample data.

faahko_3 <- faahko[3]
spectra(faahko_3)
## MSn data (Spectra) with 1278 spectra in a MsBackendMzR backend:
##        msLevel     rtime scanIndex
##      <integer> <numeric> <integer>
## 1            1   2501.38         1
## 2            1   2502.94         2
## 3            1   2504.51         3
## 4            1   2506.07         4
## 5            1   2507.64         5
## ...        ...       ...       ...
## 1274         1   4493.56      1274
## 1275         1   4495.13      1275
## 1276         1   4496.69      1276
## 1277         1   4498.26      1277
## 1278         1   4499.82      1278
##  ... 33 more variables/columns.
## 
## file(s):
## ko21.CDF
sampleData(faahko_3)
## DataFrame with 1 row and 3 columns
##   sample_name sample_group spectraOrigin
##   <character>  <character>   <character>
## 1        ko21           KO /home/bioc...

As a first evaluation of the data we below plot the base peak chromatogram (BPC) for each file in our experiment. We use the chromatogram() method and set the aggregationFun to "max" to return for each spectrum the maximal intensity and hence create the BPC from the raw data. To create a total ion chromatogram we could set aggregationFun to "sum".

## Get the base peak chromatograms. This reads data from the files.
bpis <- chromatogram(faahko, aggregationFun = "max")
## Define colors for the two groups
group_colors <- paste0(brewer.pal(3, "Set1")[1:2], "60")
names(group_colors) <- c("KO", "WT")

## Plot all chromatograms.
plot(bpis, col = group_colors[sampleData(faahko)$sample_group])
Base peak chromatogram.

Figure 1: Base peak chromatogram

The chromatogram() method returned a MChromatograms object that organizes individual Chromatogram objects (which in fact contain the chromatographic data) in a two-dimensional array: columns represent samples and rows (optionally) m/z and/or retention time ranges. Below we extract the chromatogram of the first sample and access its retention time and intensity values.

bpi_1 <- bpis[1, 1]
rtime(bpi_1) |> head()
## [1] 2501.378 2502.943 2504.508 2506.073 2507.638 2509.203
intensity(bpi_1) |> head()
## [1] 43888 43960 43392 42632 42200 42288

From the BPC above it seems that after around 4200 seconds no signal is measured anymore. Thus, we filter below the full data set to a retention time range from 2550 to 4250 seconds using the filterRt() function. Note that at present this will only subset the spectra within the MsExperiment. Subsequently we re-create also the BPC.

faahko <- filterRt(faahko, rt = c(2550, 4250))
## Filter spectra
## creating the BPC on the subsetted data
bpis <- chromatogram(faahko, aggregationFun = "max")

We next create boxplots representing the distribution of the total ion currents per data file. Such plots can be very useful to spot potentially problematic MS runs. To extract this information, we use the tic() function on the Spectra object within faahko and split the values by file using fromFile().

## Get the total ion current by file
tc <- spectra(faahko) |>
    tic() |>
    split(f = fromFile(faahko))
boxplot(tc, col = group_colors[sampleData(faahko)$sample_group],
        ylab = "intensity", main = "Total ion current")
Distribution of total ion currents per file.

Figure 2: Distribution of total ion currents per file

In addition, we can also cluster the samples based on similarity of their base peak chromatograms. Samples would thus be grouped based on similarity of their LC runs. For that we need however to bin the data along the retention time axis, since retention times will generally differ between samples. Below we use the bin() function on the BPC to bin intensities into 2 second wide retention time bins. The clustering is then performed using complete linkage hierarchical clustering on the pairwise correlations of the binned base peak chromatograms.

## Bin the BPC
bpis_bin <- bin(bpis, binSize = 2)

## Calculate correlation on the log2 transformed base peak intensities
cormat <- cor(log2(do.call(cbind, lapply(bpis_bin, intensity))))
colnames(cormat) <- rownames(cormat) <- bpis_bin$sample_name

## Define which phenodata columns should be highlighted in the plot
ann <- data.frame(group = bpis_bin$sample_group)
rownames(ann) <- bpis_bin$sample_name

## Perform the cluster analysis
pheatmap(cormat, annotation = ann,
         annotation_color = list(group = group_colors))
Grouping of samples based on similarity of their base peak chromatogram.

Figure 3: Grouping of samples based on similarity of their base peak chromatogram

The samples cluster in a pairwise manner, with the KO and WT samples for the same sample index having the most similar BPC.

2.3 Chromatographic peak detection

Chromatographic peak detection aims at identifying all signal in each sample created from ions of the same originating compound species. Chromatographic peak detection can be performed in xcms with the findChromPeaks() function and a parameter object which defines and configures the algorithm that should be used (see ?findChromPeaks for a list of supported algorithms). Before running any peak detection it is however strongly suggested to first visually inspect the extracted ion chromatogram of e.g. internal standards or compounds known to be present in the samples in order to evaluate and adapt the settings of the peak detection algorithm since the default settings will not be appropriate for most LC-MS setups.

Below we extract the EIC for one compound using the chromatogram() function by specifying in addition the m/z and retention time range where we would expect the signal for that compound.

## Define the rt and m/z range of the peak area
rtr <- c(2700, 2900)
mzr <- c(334.9, 335.1)
## extract the chromatogram
chr_raw <- chromatogram(faahko, mz = mzr, rt = rtr)
plot(chr_raw, col = group_colors[chr_raw$sample_group])
Extracted ion chromatogram for one peak.

Figure 4: Extracted ion chromatogram for one peak

Note that Chromatogram objects extracted by the chromatogram() method contain an NA value if in a certain scan (i.e. in a spectrum for a specific retention time) no signal was measured in the respective m/z range. This is reflected by the lines not being drawn as continuous lines in the plot above.

The peak above has thus a width of about 50 seconds. We can use this information to define the peakwidth parameter of the centWave peak detection method [5] that we will use for chromatographic peak detection on our data. The peakwidth parameter allows to define the expected lower and upper width (in retention time dimension) of chromatographic peaks. For our data we set it thus to peakwidth = c(20, 80). The second important parameter for centWave is ppm which defines the expected maximum deviation of m/z values of the centroids that should be included into one chromatographic peak (note that this is not the precision of m/z values provided by the MS instrument manufacturer).

For the ppm parameter we extract the full MS data (intensity, retention time and m/z values) corresponding to the above peak. To this end we first filter the raw object by retention time, then by m/z and finally plot the object with type = "XIC" to produce the plot below. We use the pipe (|>) operator to better illustrate the corresponding workflow.

faahko |>
filterRt(rt = rtr) |>
filterMz(mz = mzr) |>
plot(type = "XIC")
Visualization of the raw MS data for one peak. For each plot: upper panel: chromatogram plotting the intensity values against the retention time, lower panel m/z against retention time plot. The individual data points are colored according to the intensity.

Figure 5: Visualization of the raw MS data for one peak
For each plot: upper panel: chromatogram plotting the intensity values against the retention time, lower panel m/z against retention time plot. The individual data points are colored according to the intensity.

In the present data there is actually no variation in the m/z values. Usually the m/z values of the individual centroids (lower panel) in the plots above would scatter around the real m/z value of the compound (in an intensity dependent manner).

The first step of the centWave algorithm defines regions of interest (ROI) that are subsequently screened for the presence of chromatographic peaks. These ROIs are defined based on the difference of m/z values of centroids from consecutive scans (spectra). In detail, centroids from consecutive scans are included into a ROI if the difference between their m/z and the mean m/z of the ROI is smaller than the user defined ppm parameter. A reasonable choice for the ppm could thus be the maximal m/z difference of data points from neighboring scans/spectra that are part of a chromatographic peak for an internal standard of known compound. It is suggested to inspect the ranges of m/z values for several compounds (either internal standards or compounds known to be present in the sample) and define the ppm parameter for centWave according to these. See also this tutorial for additional information and examples on choosing and testing peak detection settings.

Chromatographic peak detection can also be performed on extracted ion chromatograms, which helps testing and tuning peak detection settings. Note however that peak detection on EICs does not involve the first step of centWave described above and will thus not consider the ppm parameter. Also, since less data is available to the algorithms, background signal estimation is performed differently and different settings for snthresh will need to be used (generally a lower snthresh will be used for EICs since the estimated background signal tends to be higher for data subsets than for the full data). Below we perform the peak detection with the findChromPeaks() function on the EIC generated above. The submitted parameter object defines which algorithm will be used and allows to define the settings for this algorithm. We use a CentWaveParam parameter object to use and configure the centWave algorithm with default settings, except for snthresh.

xchr <- findChromPeaks(chr_raw, param = CentWaveParam(snthresh = 2))

We can access the identified chromatographic peaks with the chromPeaks() function.

chromPeaks(xchr)
##        mz mzmin mzmax       rt    rtmin    rtmax        into        intb  maxo
## mzmin 335 334.9 335.1 2781.505 2761.160 2809.674  412134.255  355516.374 16856
## mzmin 335 334.9 335.1 2786.199 2764.290 2812.803 1496244.206 1391821.332 58736
## mzmin 335 334.9 335.1 2734.556 2714.211 2765.855   21579.367   18449.428   899
## mzmin 335 334.9 335.1 2797.154 2775.245 2815.933  159058.782  150289.314  6295
## mzmin 335 334.9 335.1 2784.635 2761.160 2808.109   54947.545   37923.532  2715
## mzmin 335 334.9 335.1 2859.752 2845.668 2878.532   13895.212   13874.868   905
## mzmin 335 334.9 335.1 2897.311 2891.051 2898.876    5457.155    5450.895   995
## mzmin 335 334.9 335.1 2819.064 2808.109 2834.713   19456.914   19438.134  1347
## mzmin 335 334.9 335.1 2789.329 2762.725 2808.109  174473.311   91114.975  8325
## mzmin 335 334.9 335.1 2786.199 2764.290 2812.803  932645.211  569236.246 35856
## mzmin 335 334.9 335.1 2792.461 2768.987 2823.760  876585.527  511324.134 27200
## mzmin 335 334.9 335.1 2789.329 2773.680 2806.544   89582.569   73871.386  5427
##         sn row column
## mzmin   13   1      1
## mzmin   20   1      2
## mzmin    4   1      3
## mzmin   12   1      3
## mzmin    2   1      4
## mzmin  904   1      4
## mzmin  994   1      4
## mzmin 1576   1      4
## mzmin    3   1      5
## mzmin    2   1      6
## mzmin    2   1      7
## mzmin    6   1      8

Parallel to the chromPeaks() matrix there is also a chromPeakData() data frame that allows to add arbitrary annotations to each chromatographic peak, such as e.g. the MS level in which the peak was detected:

chromPeakData(xchr)
## DataFrame with 12 rows and 4 columns
##        ms_level is_filled       row    column
##       <integer> <logical> <integer> <integer>
## mzmin         1     FALSE         1         1
## mzmin         1     FALSE         1         2
## mzmin         1     FALSE         1         3
## mzmin         1     FALSE         1         3
## mzmin         1     FALSE         1         4
## ...         ...       ...       ...       ...
## mzmin         1     FALSE         1         4
## mzmin         1     FALSE         1         5
## mzmin         1     FALSE         1         6
## mzmin         1     FALSE         1         7
## mzmin         1     FALSE         1         8

Below we plot the EIC along with all identified chromatographic peaks using the plot() function on the result object from above. Additional parameters peakCol and peakBg allow to define a foreground and background (fill) color for each identified chromatographic peak in the provided result object (i.e., we need to define one color for each row of chromPeaks(xchr) - column "column" (or "sample" if present) in that peak matrix specifies the sample in which the peak was identified).

## Define a color for each sample
sample_colors <- group_colors[xchr$sample_group]
## Define the background color for each chromatographic peak
bg <- sample_colors[chromPeaks(xchr)[, "column"]]
## Parameter `col` defines the color of each sample/line, `peakBg` of each
## chromatographic peak.
plot(xchr, col = sample_colors, peakBg = bg)
Signal for an example peak. Red and blue colors represent KO and wild type samples, respectively. Peak area of identified chromatographic peaks are highlighted in the sample group color.

Figure 6: Signal for an example peak
Red and blue colors represent KO and wild type samples, respectively. Peak area of identified chromatographic peaks are highlighted in the sample group color.

If we are happy with the settings we can use them for the peak detection on the full data set (i.e. on the MsExperiment object with the data for the whole experiment). Note that below we set the argument prefilter to c(6, 5000) and noise to 5000 to reduce the run time of this vignette. With this setting we consider only ROIs with at least 6 centroids with an intensity larger than 5000 for the centWave chromatographic peak detection.

cwp <- CentWaveParam(peakwidth = c(20, 80), noise = 5000,
                     prefilter = c(6, 5000))
faahko <- findChromPeaks(faahko, param = cwp)

The results of findChromPeaks() on a MsExperiment object are returned as an XcmsExperiment object. This object extends MsExperiment directly (hence providing the same access to all raw data) and contains all xcms preprocessing results. Note also that additional rounds of chromatographic peak detections could be performed and their results being added to existing peak detection results by additional calls to findChromPeaks() on the result object and using parameter add = TRUE.

The chromPeaks function can also here be used to access the results from the chromatographic peak detection. Below we show the first 6 identified chromatographic peaks.

chromPeaks(faahko) |>
    head()
##           mz mzmin mzmax       rt    rtmin    rtmax     into     intb  maxo sn
## CP0001 594.0 594.0 594.0 2601.535 2581.191 2637.529 161042.2 146073.3  7850 11
## CP0002 577.0 577.0 577.0 2604.665 2581.191 2626.574 136105.2 128067.9  6215 11
## CP0003 307.0 307.0 307.0 2618.750 2592.145 2645.354 284782.4 264907.0 16872 20
## CP0004 302.0 302.0 302.0 2617.185 2595.275 2640.659 687146.6 669778.1 30552 43
## CP0005 370.1 370.1 370.1 2673.523 2643.789 2700.127 449284.6 417225.3 25672 17
## CP0006 427.0 427.0 427.0 2675.088 2643.789 2684.478 283334.7 263943.2 11025 13
##        sample
## CP0001      1
## CP0002      1
## CP0003      1
## CP0004      1
## CP0005      1
## CP0006      1

Columns of this chromPeaks() matrix might differ depending on the used peak detection algorithm. Columns that all algorithms have to provide are: "mz", "mzmin", "mzmax", "rt", "rtmin" and "rtmax" that define the m/z and retention time range of the chromatographic peak (i.e. all mass peaks within that area are used to integrate the peak signal) as well as the m/z and retention time of the peak apex. Other mandatory columns are "into" and "maxo" with the absolute integrated peak signal and the maximum peak intensity. Finally, "sample" provides the index of the sample in which the peak was identified.

Additional annotations for each individual peak can be extracted with the chromPeakData() function. This data frame could also be used to add/store arbitrary annotations for each detected peak (that don’t necessarily need to be numeric).

chromPeakData(faahko)
## DataFrame with 2908 rows and 2 columns
##         ms_level is_filled
##        <integer> <logical>
## CP0001         1     FALSE
## CP0002         1     FALSE
## CP0003         1     FALSE
## CP0004         1     FALSE
## CP0005         1     FALSE
## ...          ...       ...
## CP2904         1     FALSE
## CP2905         1     FALSE
## CP2906         1     FALSE
## CP2907         1     FALSE
## CP2908         1     FALSE

2.3.1 Chromatographic peak quality

Based on the publication by Kumler et al. published in 2023 [6], new quality metrics (beta_cor and beta_snr) were added to xcms. beta_cor indicates how bell-shaped the chromatographic peak is and beta_snr is estimating the signal-to-noise ratio using the residuals from this fit. These metrics can be calculated during peak picking by setting verboseBetaColumns = TRUE in the CentWaveParam object, or they can be calculated afterwards by using the chromPeakSummary() function with the XcmsExperiment object and the BetaDistributionParam parameter object as input:

beta_metrics <- chromPeakSummary(faahko, BetaDistributionParam())
head(beta_metrics)
##         beta_cor beta_snr
## CP0001 0.9865868 5.349210
## CP0002 0.9401467 4.928835
## CP0003 0.9846982 5.783004
## CP0004 0.9883246 6.044377
## CP0005 0.6972246 5.293175
## CP0006 0.1358883 4.995241

The result returned by chromPeakSummary() is thus a numeric matrix with the values for these quality estimates, one row for each chromatographic peak. Using summary statistics, one can explore the distribution of these metrics in the data.

summary(beta_metrics)
##     beta_cor          beta_snr    
##  Min.   :-0.8541   Min.   :3.921  
##  1st Qu.: 0.7390   1st Qu.:5.157  
##  Median : 0.9442   Median :5.594  
##  Mean   : 0.7904   Mean   :5.662  
##  3rd Qu.: 0.9771   3rd Qu.:6.097  
##  Max.   : 0.9990   Max.   :7.983  
##  NA's   :6         NA's   :6

Visual inspection gives a better idea of what these metrics represent in terms of peak quality in the data at hand. This information can be used to e.g. filter out peaks that don’t meet a chosen quality metric threshold. In order to plot the detected peaks, their EIC can be extracted with the function chromPeakChromatograms. An example of a peak with a high beta_cor and for a peak with a low beta_cor score is given below.

beta_metrics[c(4, 6), ]
##         beta_cor beta_snr
## CP0004 0.9883246 6.044377
## CP0006 0.1358883 4.995241
eics <- chromPeakChromatograms(
    faahko, peaks = rownames(chromPeaks(faahko))[c(4, 6)])
peak_1 <- eics[1]
peak_2 <- eics[2]
par(mfrow = c(1, 2))
plot(peak_1)
plot(peak_2)
Plots of high and low quality peaks. Left: peak CP0004 with a beta_cor = 0.98, right: peak CP0006 with a beta_cor = 0.13.

Figure 7: Plots of high and low quality peaks
Left: peak CP0004 with a beta_cor = 0.98, right: peak CP0006 with a beta_cor = 0.13.

2.3.2 Refining peak detection

Peak detection will not always work perfectly for all types of peak shapes present in the data set leading to peak detection artifacts, such as (partially or completely) overlapping peaks or artificially split peaks (common issues especially for centWave). xcms provides the refineChromPeaks() function that can be called on peak detection results in order to refine (or clean) peak detection results by either removing identified peaks not passing a certain criteria or by merging artificially split or partially or completely overlapping chromatographic peaks. Different algorithms are available that can again be configured with their respective parameter objects: CleanPeaksParam and FilterIntensityParam allow to remove peaks with their retention time range or intensity being below a specified threshold, respectively. With MergeNeighboringPeaksParam it is possible to merge chromatographic peaks and hence remove many of the above mentioned (centWave) peak detection artifacts. See also ?refineChromPeaks for more information and help on the different methods.

Below we post-process the peak detection results merging peaks that overlap in a 4 second window per file and for which the signal between them is lower than 75% of the smaller peak’s maximal intensity. See the ?MergeNeighboringPeaksParam help page for a detailed description of the settings and the approach.

mpp <- MergeNeighboringPeaksParam(expandRt = 4)
faahko_pp <- refineChromPeaks(faahko, mpp)

An example for a merged peak is given below.

mzr_1 <- 305.1 + c(-0.01, 0.01)
chr_1 <- chromatogram(faahko[1], mz = mzr_1)
## Processing chromatographic peaks
chr_2 <- chromatogram(faahko_pp[1], mz = mzr_1)
## Processing chromatographic peaks
par(mfrow = c(1, 2))
plot(chr_1)
plot(chr_2)
Result from the peak refinement step. Left: data before processing, right: after refinement. The splitted peak was merged into one.

Figure 8: Result from the peak refinement step
Left: data before processing, right: after refinement. The splitted peak was merged into one.

centWave thus detected originally 3 chromatographic peaks in the m/z slice (left panel in the plot above) and peak refinement with MergeNeighboringPeaksParam and the specified settings merged the first two peaks into a single one (right panel in the figure above). Other close peaks, with a lower intensity between them, were however not merged (see below).

mzr_1 <- 496.2 + c(-0.01, 0.01)
chr_1 <- chromatogram(faahko[1], mz = mzr_1)
## Processing chromatographic peaks
chr_2 <- chromatogram(faahko_pp[1], mz = mzr_1)
## Processing chromatographic peaks
par(mfrow = c(1, 2))
plot(chr_1)
plot(chr_2)
Result from the peak refinement step. Left: data before processing, right: after refinement. The peaks were not merged.

Figure 9: Result from the peak refinement step
Left: data before processing, right: after refinement. The peaks were not merged.

It is also possible to perform the peak refinement on extracted ion chromatograms. This could again be used to test and fine-tune the settings for the parameter and to avoid potential problematic behavior. The minProp parameter for example has to be carefully chosen to avoid merging of isomer peaks (like in the example above). With the default minProp = 0.75 only peaks are merged if the signal between the two peaks is higher than 75% of the smaller peak’s maximal intensity. Setting this value too low could eventually result in merging of isomers as shown below.

#' Too low minProp could cause merging of isomers!
res <- refineChromPeaks(chr_1, MergeNeighboringPeaksParam(minProp = 0.05))
chromPeaks(res)
##         mz  mzmin  mzmax       rt    rtmin    rtmax     into intb    maxo  sn
## CPM1 496.2 496.19 496.21 3384.012 3294.809 3412.181 45940118   NA 1128960 177
##      sample row column
## CPM1      1   1      1
plot(res)

Thus, before running such a peak refinement evaluate that isomers present in the data set were not wrongly merged based on the chosen settings.

Before proceeding we next replace the faahko object with the results from the peak refinement step.

faahko <- faahko_pp

Below we use the data from the chromPeaks() matrix to calculate per-file summaries of the peak detection results, such as the number of peaks per file as well as the distribution of the retention time widths.

summary_fun <- function(z)
    c(peak_count = nrow(z), rt = quantile(z[, "rtmax"] - z[, "rtmin"]))

T <- chromPeaks(faahko) |>
    split.data.frame(f = chromPeaks(faahko)[, "sample"]) |>
    lapply(FUN = summary_fun) |>
    do.call(what = rbind)
rownames(T) <- basename(fileNames(faahko))
pandoc.table(
    T,
    caption = paste0("Summary statistics on identified chromatographic",
                     " peaks. Shown are number of identified peaks per",
                     " sample and widths/duration of chromatographic ",
                     "peaks."))
Summary statistics on identified chromatographic peaks. Shown are number of identified peaks per sample and widths/duration of chromatographic peaks.
  peak_count rt.0% rt.25% rt.50% rt.75% rt.100%
ko15.CDF 397 10.95 34.43 42.25 53.21 319.2
ko16.CDF 520 10.95 32.86 42.25 53.21 151.8
ko21.CDF 207 10.95 42.25 50.08 65.73 164.3
ko22.CDF 233 10.95 37.56 46.95 59.47 147.1
wt15.CDF 403 10.95 32.86 42.25 53.21 161.2
wt16.CDF 361 10.95 35.99 45.38 57.9 162.8
wt21.CDF 227 10.95 35.21 48.51 64.16 172.1
wt22.CDF 328 10.95 35.99 45.38 57.9 228.5

While by default chromPeaks() will return all identified chromatographic peaks in a result object it is also possible to extract only chromatographic peaks for a specified m/z and/or rt range:

chromPeaks(faahko, mz = c(334.9, 335.1), rt = c(2700, 2900))
##         mz mzmin mzmax       rt    rtmin    rtmax      into      intb  maxo sn
## CP0038 335   335   335 2781.505 2761.160 2809.674  412134.3  383167.4 16856 23
## CP0494 335   335   335 2786.199 2764.290 2812.803 1496244.2 1461187.3 58736 72
## CP1025 335   335   335 2797.154 2775.245 2815.933  159058.8  149229.6  6295 13
## CP1964 335   335   335 2786.199 2764.290 2812.803  932645.2  915333.8 35856 66
## CP2349 335   335   335 2792.461 2768.987 2823.760  876585.5  848569.1 27200 36
##        sample
## CP0038      1
## CP0494      2
## CP1025      3
## CP1964      6
## CP2349      7

We can also plot the location of the identified chromatographic peaks in the m/z - retention time space for one file using the plotChromPeaks() function. Below we plot this information for the third sample.

plotChromPeaks(faahko, file = 3)
Identified chromatographic peaks in the m/z by retention time space for one sample.

Figure 10: Identified chromatographic peaks in the m/z by retention time space for one sample

As a general overview of the peak detection results we can in addition visualize the number of identified chromatographic peaks per file along the retention time axis. Parameter binSize allows to define the width of the bins in rt dimension in which peaks should be counted. This number of chromatographic peaks within each bin is then shown color-coded in the resulting plot.

plotChromPeakImage(faahko, binSize = 10)
Frequency of identified chromatographic peaks along the retention time axis. The frequency is color coded with higher frequency being represented by yellow-white. Each line shows the peak frequency for one file.

Figure 11: Frequency of identified chromatographic peaks along the retention time axis
The frequency is color coded with higher frequency being represented by yellow-white. Each line shows the peak frequency for one file.

Note that extracting ion chromatorams from an xcms result object will in addition to the chromatographic data also extract any identified chromatographic peaks within that region. This can thus also be used to validate and verify that the used peak detection settings identified e.g. peaks for known compounds or internal standards properly. Below we extract the ion chromatogram for the m/z - rt region above and access the detected peaks in that region using the chromPeaks() function.

chr_ex <- chromatogram(faahko, mz = mzr, rt = rtr)
chromPeaks(chr_ex)
##         mz mzmin mzmax       rt    rtmin    rtmax      into      intb  maxo sn
## CP0038 335   335   335 2781.505 2761.160 2809.674  412134.3  383167.4 16856 23
## CP0494 335   335   335 2786.199 2764.290 2812.803 1496244.2 1461187.3 58736 72
## CP1025 335   335   335 2797.154 2775.245 2815.933  159058.8  149229.6  6295 13
## CP1964 335   335   335 2786.199 2764.290 2812.803  932645.2  915333.8 35856 66
## CP2349 335   335   335 2792.461 2768.987 2823.760  876585.5  848569.1 27200 36
##        sample row column
## CP0038      1   1      1
## CP0494      2   1      2
## CP1025      3   1      3
## CP1964      6   1      6
## CP2349      7   1      7

We can also plot this extracted ion chromatogram which will also visualize all identified chromatographic peaks in that region.

sample_colors <- group_colors[chr_ex$sample_group]
plot(chr_ex, col = group_colors[chr_raw$sample_group], lwd = 2,
     peakBg = sample_colors[chromPeaks(chr_ex)[, "sample"]])
Signal for an example peak. Red and blue colors represent KO and wild type samples, respectively. The signal area of identified chromatographic peaks are filled with a color.

Figure 12: Signal for an example peak
Red and blue colors represent KO and wild type samples, respectively. The signal area of identified chromatographic peaks are filled with a color.

Chromatographic peaks can be visualized also in other ways: by setting peakType = "rectangle" the they are indicated with a rectangle instead of the default highlighting option (peakType = "polygon") used above. As a third alternative it would also possible to just indicate the apex position for each identified chromatographic peak with a single point (peakType = "point"). Below we plot the data again using peakType = "rectangle".

plot(chr_ex, col = sample_colors, peakType = "rectangle",
     peakCol = sample_colors[chromPeaks(chr_ex)[, "sample"]],
     peakBg = NA)
Signal for an example peak. Red and blue colors represent KO and wild type samples, respectively. The rectangles indicate the identified chromatographic peaks per sample.

Figure 13: Signal for an example peak
Red and blue colors represent KO and wild type samples, respectively. The rectangles indicate the identified chromatographic peaks per sample.

Finally we plot also the distribution of peak intensity per sample. This allows to investigate whether systematic differences in peak signals between samples are present.

## Extract a list of per-sample peak intensities (in log2 scale)
ints <- split(log2(chromPeaks(faahko)[, "into"]),
              f = chromPeaks(faahko)[, "sample"])
boxplot(ints, varwidth = TRUE, col = sample_colors,
        ylab = expression(log[2]~intensity), main = "Peak intensities")
grid(nx = NA, ny = NULL)
Peak intensity distribution per sample.

Figure 14: Peak intensity distribution per sample

Over and above the signal of the identified chromatographic peaks is comparable across samples, with the exception of samples 3, 4 and, to some degree, also 7 that show slightly higher average intensities, but also a lower number of detected peaks (indicated by the smaller width of the boxes).

Note that in addition to the above described identification of chromatographic peaks, it is also possible to manually define and add chromatographic peaks with the manualChromPeaks() function (see ?manualChromPeaks help page for more information).

2.4 Alignment

The time at which analytes elute in the chromatography can vary between samples (and even compounds). Such differences were already visible in the BPC and even the extracted ion chromatogram plot in the previous section. The alignment step, also referred to as retention time correction, aims to adjust these differences by shifting signals along the retention time axis and aligning them between different samples within an experiment.

A plethora of alignment algorithms exist (see [7]), with some of them being also implemented in xcms. Alignment of LC-MS data can be performed in xcms using the adjustRtime() method and an algorithm-specific parameter class (see ?adjustRtime for an overview of available methods in xcms).

In the example below we use the obiwarp method [8] to align the samples. We use a binSize = 0.6 which creates warping functions in m/z bins of 0.6. Also here it is advisable to modify and adapt the settings for each experiment.

faahko <- adjustRtime(faahko, param = ObiwarpParam(binSize = 0.6))

Note that adjustRtime(), besides calculating adjusted retention times for each spectrum, adjusts also the retention times of the identified chromatographic peaks in the xcms result object. Adjusted retention times of individual spectra can be extracted from the result object using either the adjustedRtime() function or using rtime() with parameter adjusted = TRUE (the default):

## Extract adjusted retention times
adjustedRtime(faahko) |> head()
## [1] 2551.457 2553.089 2554.720 2556.352 2557.983 2559.615
## Or simply use the rtime method
rtime(faahko) |> head()
## [1] 2551.457 2553.089 2554.720 2556.352 2557.983 2559.615
## Get raw (unadjusted) retention times
rtime(faahko, adjusted = FALSE) |> head()
## [1] 2551.457 2553.022 2554.586 2556.151 2557.716 2559.281

To evaluate the impact of the alignment we plot the BPC on the adjusted data. In addition we plot also the differences between the adjusted and the raw retention times per sample using the plotAdjustedRtime() function. To disable the automatic extraction of all identified chromatographic peaks by the chromatogram() function (which would not make much sense for a BPC) we use chromPeaks = "none" below.

## Get the base peak chromatograms.
bpis_adj <- chromatogram(faahko, aggregationFun = "max", chromPeaks = "none")
par(mfrow = c(3, 1), mar = c(4.5, 4.2, 1, 0.5))
plot(bpis, col = sample_colors)
grid()
plot(bpis_adj, col = sample_colors)
grid()
## Plot also the difference of adjusted to raw retention time.
plotAdjustedRtime(faahko, col = sample_colors)
grid()
Obiwarp aligned data. Base peak chromatogram before (top) and after alignment (middle) and difference between adjusted and raw retention times along the retention time axis (bottom).

Figure 15: Obiwarp aligned data
Base peak chromatogram before (top) and after alignment (middle) and difference between adjusted and raw retention times along the retention time axis (bottom).

Too large differences between adjusted and raw retention times could indicate poorly performing samples or alignment.

At last we evaluate also the impact of the alignment on the test peak.

par(mfrow = c(2, 1))
## Plot the raw data
plot(chr_raw, col = sample_colors)
grid()
## Extract the chromatogram from the adjusted object
chr_adj <- chromatogram(faahko, rt = rtr, mz = mzr)
plot(chr_adj, col = sample_colors, peakType = "none")
grid()
Example extracted ion chromatogram before (top) and after alignment (bottom).

Figure 16: Example extracted ion chromatogram before (top) and after alignment (bottom)

Note: xcms result objects (XcmsExperiment as well as XCMSnExp) contain both the raw and the adjusted retention times for all spectra and subset operation will in many cases drop adjusted retention times. Thus it might sometimes be useful to immediately replace the raw retention times in the data using the applyAdjustedRtime() function.

2.4.1 Subset-based alignment

For some experiments it might be better to perform the alignment based on only a subset of the available samples, e.g. if pooled QC samples were injected at regular intervals or if the experiment contains blanks. All alignment methods in xcms support such a subset-based alignment in which retention time shifts are estimated on only a specified subset of samples followed by an alignment of the whole data set based on the aligned subset.

The subset of samples for such an alignment can be specified with the parameter subset of the PeakGroupsParam or ObiwarpParam object. Parameter subsetAdjust allows to specify the method by which the left-out samples will be adjusted. There are currently two options available:

  • subsetAdjust = "previous": adjust the retention times of a non-subset sample based on the alignment results of the previous subset sample (e.g. a QC sample). If samples are e.g. in the order A1, B1, B2, A2, B3, B4 with A representing QC samples and B study samples, using subset = c(1, 4) and subsetAdjust = "previous" would result in all A samples to be aligned with each other and non-subset samples B1 and B2 being adjusted based on the alignment result of subset samples A1 and B3 and B4 on those of A2.

  • subsetAdjust = "average": adjust retention times of non-subset samples based on an interpolation of the alignment results of the previous and subsequent subset sample. In the example above, B1 would be adjusted based on the average of adjusted retention times between subset (QC) samples A1 and A2. Since there is no subset sample after non-subset samples B3 and B4 these will be adjusted based on the alignment results of A2 alone. Note that a weighted average is used to calculate the adjusted retention time averages, which uses the inverse of the difference of the index of the non-subset sample to the subset samples as weights. Thus, if we have a setup like A1, B1, B2, A2 the adjusted retention times of A1 would get a larger weight than those of A2 in the adjustment of non-subset sample B1 causing it’s adjusted retention times to be closer to those of A1 than to A2. See below for examples.

Important: both cases require a meaningful/correct ordering of the samples within the object (i.e., samples should be ordered by injection index).

The examples below aim to illustrate the effect of these alignment options. We assume samples 1, 4 and 7 in the faahKO data set to be QC pool samples. We thus want to perform the alignment based on these samples and subsequently adjust the retention times of the left-out samples (2, 3, 5, 6 and 8) based on interpolation of the results from the neighboring subset (QC) samples. After initial peak grouping we perform the subset-based alignment with the peak groups method by passing the indices of the QC samples with the subset parameter to the PeakGroupsParam function and set subsetAdjust = "average" to adjust the study samples based on interpolation of the alignment results from neighboring subset/QC samples.

Note that for any subset-alignment all parameters such as minFraction are relative to the subset, not the full experiment!

Below we first remove any previous alignment results with the dropAdjustedRtime() function to allow a fresh alignment using the subset-based option outlined above. In addition to removing adjusted retention times for all spectra, this function will also restore the original retention times for identified chromatographic peaks.

faahko <- dropAdjustedRtime(faahko)

## Define the experimental layout
sampleData(faahko)$sample_type <- "study"
sampleData(faahko)$sample_type[c(1, 4, 7)] <- "QC"

For an alignment with the peak groups method an initial peak grouping (correspondence) analysis is required, because the algorithm estimates retention times shifts between samples using the retention times of hook peaks (i.e. chromatographic peaks present in most/all samples). Here we use the default settings for an peak density method-based correspondence, but it is strongly advised to adapt the parameters for each data set (details in the next section). The definition of the sample groups (i.e. assignment of individual samples to the sample groups in the experiment) is mandatory for the PeakDensityParam. If there are no sample groups in the experiment, sampleGroups should be set to a single value for each file (e.g. rep(1, length(fileNames(faahko))).

## Initial peak grouping. Use sample_type as grouping variable
pdp_subs <- PeakDensityParam(sampleGroups = sampleData(faahko)$sample_type,
                             minFraction = 0.9)
faahko <- groupChromPeaks(faahko, param = pdp_subs)

## Define subset-alignment options and perform the alignment
pgp_subs <- PeakGroupsParam(
    minFraction = 0.85,
    subset = which(sampleData(faahko)$sample_type == "QC"),
    subsetAdjust = "average", span = 0.4)
faahko <- adjustRtime(faahko, param = pgp_subs)

Below we plot the results of the alignment highlighting the subset samples in green. This nicely shows how the interpolation of the subsetAdjust = "average" works: retention times of sample 2 are adjusted based on those from subset sample 1 and 4, giving however more weight to the closer subset sample 1 which results in the adjusted retention times of 2 being more similar to those of sample 1. Sample 3 on the other hand gets adjusted giving more weight to the second subset sample (4).

clrs <- rep("#00000040", 8)
clrs[sampleData(faahko)$sample_type == "QC"] <- c("#00ce0080")
par(mfrow = c(2, 1), mar = c(4, 4.5, 1, 0.5))
plot(chromatogram(faahko, aggregationFun = "max", chromPeaks = "none"),
     col = clrs)
grid()
plotAdjustedRtime(faahko, col = clrs, peakGroupsPch = 1,
                  peakGroupsCol = "#00ce0040")
grid()