XCMS workshop 2017

Why are we here and why do we need to do this?

Simply to do. Difficult to do well. Impossible to do perfectly.

• R package
• First published in 2006[1] by Colin A. Smith et al. at The Scripps Research Institute (Patti Lab).
• Open source: https://github.com/sneumann/xcms
• The 'X' in the XCMS denotes that it can be used for an chromatography (GC/LC)
• Maintained for many years primarily by Steffen Neumann and his group
• New major version being developed by Johannes Rainer (EURAC Research , Bolzano, Italy)

Is this a peak or just noise?

First we locate the files. They need to be in an open format:

• mzML (newest)
• mzData
• mzXML (most widely supported)
• netCDF (obsolete, last resort)

Lets make a list of the files we want to pre-process.

files <- list.files("_data", recursive = TRUE, full.names = TRUE, pattern=".mzXML")
head(files,5)

[1] "_data/A/POOL_1_A_1.mzXML" "_data/A/POOL_1_A_2.mzXML"
[3] "_data/A/POOL_1_A_3.mzXML" "_data/A/POOL_2_A_1.mzXML"
[5] "_data/A/POOL_2_A_2.mzXML"

Blanks interfer with peak alignment since there are so few peaks.
So, better remove them before starting. Same goes for compound mixtures.

Now we can do peak picking. We need to load some packages first.

suppressMessages(library(xcms)) # suppress avoids a flood of messages
suppressMessages(library(BiocParallel))

Now lets do the peak picking.
There are many parameters to set that we will explain in a sec.
There are a few more settings available but these are the most important.

xset <- xcmsSet(files, 
                BPPARAM  = SerialParam(),
                method = 'centWave',
                prefilter = c(3,1E3),
                ppm = 30,
                snthr = 1E2,
                profparam = list(step=0.01),
                peakwidth = c(0.05*60,0.20*60),
                verbose.columns = TRUE,
                fitgauss = TRUE
                )

By default the newest XCMS uses all available cores.
Before we set it to use a single core (SerialParam) just because it works better for generating this presentation.

To control the number of cores you could for example do:

library(parallel)

and set:

BPPARAM  = SnowParam(detectCores()-1, progressbar = TRUE)

• SnowParam() on windows
• MulticoreParam() on linux

Here we will go through peakpicking with the centwave algorithm[2].

• First it finds regions of interest (ROIs).
  Those are where there could be a peak: Stable mass for a certain time.
• Next the peak shape is evaluated.
• The peak is expanded beyond the ROI if needed.

Figure from Tautenhahn, Böttcher and Neumann, 2008 (centWave paper)

You can also use a 2D plot to try to find short and long peaks. Here I plotted with MzMine (remember to use the continuous/profile mode toggle otherwise it looks very wrong).

ppm: Relative deviation in the m/z dimension
    In the centWave context it is the maximum allowed deviation between scans when locating regions of interest (ROIs).

Do not use the vendors number. Like the mileage of a car these numbers are far from real world scenarios.
We need a range that is true also for the tails of the peaks.

What we choose before was:

ppm = 30

A 2D plot is a reasonable way to look at this.

Check the peak bounderies with an EIC.

Check what are actually reasonably sized peaks in the spectrum.

((189.074-189.0690)/189.0690)*1e6

[1] 26.44537

First we need to clear up the confusion about the difference between profile mode data and the profile matrix.

Profile mode: The data is continuous in the m/z dimension.
As opposed to centroid mode data where you have discrete (sticks) in the m/z dimension.

Profile matrix: A rt (or scan rather) x m/z matrix of m/z slices.
XCMS uses this for some procedures (fillPeaks, but not peakpicking) instead of the full raw data. Think of it as binned data with less m/z resolution.

• The profparam parameter says how fine the binning in the Profile matrix is going to be.
• It cannot be set too low as slices will be combined as needed. Only memory usage will increase.
• For high dimensional data the default setting can cause problems. So setting something like this will set it to 0.01 Da slices.

profparam = list(step=0.01)

Signal to noise ratio. It really depends on your data. And it is defined differently for centWave and matchedfilter.
It appears to be very unstable even for very similar peaks. So, I suggest to set it low and filter on intensity afterwards if needed.

For centWave I would start with:

snthr = 1E2

For matchedfilter I would start with:

snthr = 5

If low peaks that you would like to include are missing then lower the value.

Is this peak the same as that peak?

This step tries to group/match features between samples.
This is probably the most important step!

It is a 3 step procedure:

1. Features are grouped according to RT. Loose parameters for matching is used.
2. The matched features are used to model RT shifts. –> retention time correction/alignment.
3. Using the corrected RTs features are matched again with stricter criteria.

There are three methods for retention time correction/alignment.

1. LOESS: Fits a LOESS curve between retention times. Works on the peaktable (therefore fast). The default and usually no reason to look further.
2. Obiwarp: Warps the raw data (and hence slow) to a reference sample. Can handle dramatic shifts. Often overfitting. Use if needed with great attention.
3. Linear: When all else fails.

xset_g <- group(xset, 
                bw = 0.2*60, 
                minfrac = 0.5, 
                minsamp = 5, 
                mzwid = 0.01, 
                max = 20, 
                sleep = 0
                )

xset_g

An "xcmsSet" object with 27 samples

Time range: 0.6-720 seconds (0-12 minutes)
Mass range: 60.5133-1418.4605 m/z
Peaks: 97334 (about 3605 per sample)
Peak Groups: 2978 
Sample classes: A, B, C 

Feature detection:
 o Peak picking performed on MS1.
Profile settings: method = bin
                  step = 0.01

Memory usage: 18.4 MB

Example of grouping density.

bw: bandwidth (standard deviation or half width at half maximum) of gaussian smoothing kernel to apply to the peak density chromatogram.

Perhaps a bit obscure definition but it roughly translates to:
The maximum expected RT deviation across samples.

Overlay the BPI and find the peak with the highest deviation.

Here we have zoomed. Compare peak apexes.

Calculate the deviation.

4.90-4.78

[1] 0.12

Since there could be small peaks with larger variation than what we found we bump it up to 0.2 min.

minfrac = 0.5
minsamp = 5

minsamp: minimum number of samples a peak should be found in.

minfrac: minimum fraction of samples a peak should be found in.

These are used per sample group.

Files in different subfolders are considered different groups.
The groups can be manipulated by changing:

head(xset@phenoData)

           class
POOL_1_A_1     A
POOL_1_A_2     A
POOL_1_A_3     A
POOL_2_A_1     A
POOL_2_A_2     A
POOL_2_A_3     A

In this dataset we had 3 groups with 9 samples in each. The first time we do the grouping we are pretty liberal with how often it should be found.

Consider your experimental design and group VERY carefully when you set these parameters!

Example:

1. You put all files on one folder.
2. You have 2 study groups with each 20 samples. So 40 samples in total.

3. You set:

   minfrac = 0.7
   minsamp = 30
   

4. You have now removed everything that is unique to one group!
   So possibly the most important features have been removed.

mzwid = 0.01

How close the m/z need to be the same compound.

max: maximum number of groups to identify in a single m/z slice.
Meaning how many isomers do you have across a chromatogram.

max = 20

If you are processing GC data remember that many fragments will be isomers. So, set this high.

xset_g_r <- retcor(xset_g,
                   method="peakgroups",
                   missing = 3, 
                   extra = 3, 
                   smooth = "loess",
                   span = 0.6, 
                   plottype = "mdevden",
                   col=xset_g@phenoData[,"class"]
                  )

span: degree of smoothing for local polynomial regression fitting.
The higher the less flexible the fitting is.
0 == overfitting
1 == only captures very overall trends

0.4-0.7 is typical in my opinion. 0.2 is default.

We set:

span = 0.6

missing = 3
extra = 3

missing: number of missing samples to allow in retention time correction groups

extra: number of extra peaks to allow in retention time correction correction groups

plottype = "mdevden",
col=xset_g@phenoData[,"class"]

plottype:

• “deviation”: Plots only the RT deviations
• “mdevden”: As “deviation” but adds density plot below

col: We can tell it how to color the points/lines in the deviations plot.

• Can be a factor, numeric or a character vector of colors.
• Useful to color by batch.

Now we do grouping again with the corrected retention times.

Notice that we give more strict bw, minfrac and minsamp since we expect better matching now that we have corrected the retention times.

xset_g_r_g <- group(xset_g_r,
                    bw = 5, 
                    minfrac = 0.75, 
                    minsamp = 7, 
                    mzwid = 0.01, 
                    max = 20, 
                    sleep = 0
                   )

xset_g_r_g

An "xcmsSet" object with 27 samples

Time range: 0.5-720.5 seconds (0-12 minutes)
Mass range: 60.5133-1418.4605 m/z
Peaks: 97334 (about 3605 per sample)
Peak Groups: 1627 
Sample classes: A, B, C 

Feature detection:
 o Peak picking performed on MS1.
Profile settings: method = bin
                  step = 0.01

Memory usage: 20.2 MB

We can redo the retcor to get the plot using corrected RTs. We won't use the returned object.

xset_g_r_g_r <- retcor(xset_g_r_g,
                       method="peakgroups",
                       missing = 3, 
                       extra = 3, 
                       smooth = "loess",
                       span = 0.6, 
                       plottype = "mdevden",
                       col=xset_g@phenoData[,"class"]
                      )

Is there really nothing here?

• Peaks that were not originally found in some samples are problematic for statistics
• Therefore we re-examine each file again and check the area where peaks were found in other files
• We then integrate (blindly) the area we expect to find a peak
• The integration is done with the slice-size we set with the profparam parameter in xcmsSet

xset_g_r_g_fill <- fillPeaks(xset_g_r_g, 
                             expand.mz=1.5,
                             expand.rt=1.1, 
                             BPPARAM = SerialParam()
                            )

C:/Users/JPZS/Desktop/gits/XCMS_course/_data/A/POOL_1_A_1.mzXML 
method:  bin 
step:  0.01 
C:/Users/JPZS/Desktop/gits/XCMS_course/_data/A/POOL_1_A_2.mzXML 
method:  bin 
step:  0.01 
C:/Users/JPZS/Desktop/gits/XCMS_course/_data/A/POOL_1_A_3.mzXML 
method:  bin 
step:  0.01 
C:/Users/JPZS/Desktop/gits/XCMS_course/_data/A/POOL_2_A_1.mzXML 
method:  bin 
step:  0.01 
C:/Users/JPZS/Desktop/gits/XCMS_course/_data/A/POOL_2_A_2.mzXML 
method:  bin 
step:  0.01 
C:/Users/JPZS/Desktop/gits/XCMS_course/_data/A/POOL_2_A_3.mzXML 
method:  bin 
step:  0.01 
C:/Users/JPZS/Desktop/gits/XCMS_course/_data/A/POOL_3_A_1.mzXML 
method:  bin 
step:  0.01 
C:/Users/JPZS/Desktop/gits/XCMS_course/_data/A/POOL_3_A_2.mzXML 
method:  bin 
step:  0.01 
C:/Users/JPZS/Desktop/gits/XCMS_course/_data/A/POOL_3_A_3.mzXML 
method:  bin 
step:  0.01 
C:/Users/JPZS/Desktop/gits/XCMS_course/_data/B/POOL_1_B_1.mzXML 
method:  bin 
step:  0.01 
C:/Users/JPZS/Desktop/gits/XCMS_course/_data/B/POOL_1_B_2.mzXML 
method:  bin 
step:  0.01 
C:/Users/JPZS/Desktop/gits/XCMS_course/_data/B/POOL_1_B_3.mzXML 
method:  bin 
step:  0.01 
C:/Users/JPZS/Desktop/gits/XCMS_course/_data/B/POOL_2_B_1.mzXML 
method:  bin 
step:  0.01 
C:/Users/JPZS/Desktop/gits/XCMS_course/_data/B/POOL_2_B_2.mzXML 
method:  bin 
step:  0.01 
C:/Users/JPZS/Desktop/gits/XCMS_course/_data/B/POOL_2_B_3.mzXML 
method:  bin 
step:  0.01 
C:/Users/JPZS/Desktop/gits/XCMS_course/_data/B/POOL_3_B_1.mzXML 
method:  bin 
step:  0.01 
C:/Users/JPZS/Desktop/gits/XCMS_course/_data/B/POOL_3_B_2.mzXML 
method:  bin 
step:  0.01 
C:/Users/JPZS/Desktop/gits/XCMS_course/_data/B/POOL_3_B_3.mzXML 
method:  bin 
step:  0.01 
C:/Users/JPZS/Desktop/gits/XCMS_course/_data/C/POOL_1_C_1.mzXML 
method:  bin 
step:  0.01 
C:/Users/JPZS/Desktop/gits/XCMS_course/_data/C/POOL_1_C_2.mzXML 
method:  bin 
step:  0.01 
C:/Users/JPZS/Desktop/gits/XCMS_course/_data/C/POOL_1_C_3.mzXML 
method:  bin 
step:  0.01 
C:/Users/JPZS/Desktop/gits/XCMS_course/_data/C/POOL_2_C_1.mzXML 
method:  bin 
step:  0.01 
C:/Users/JPZS/Desktop/gits/XCMS_course/_data/C/POOL_2_C_2.mzXML 
method:  bin 
step:  0.01 
C:/Users/JPZS/Desktop/gits/XCMS_course/_data/C/POOL_2_C_3.mzXML 
method:  bin 
step:  0.01 
C:/Users/JPZS/Desktop/gits/XCMS_course/_data/C/POOL_3_C_1.mzXML 
method:  bin 
step:  0.01 
C:/Users/JPZS/Desktop/gits/XCMS_course/_data/C/POOL_3_C_2.mzXML 
method:  bin 
step:  0.01 
C:/Users/JPZS/Desktop/gits/XCMS_course/_data/C/POOL_3_C_4.mzXML 
method:  bin 
step:  0.01 

How do I look at the results and how do I know if they are good?

There are two types of peak-tables in XCMS.
First we look at the aligned table that is what you normally want.

peaks_aligned <- peakTable(xset_g_r_g_fill)

kable(peaks_aligned[1:5,])

• npeaks: Number of peaks assigned to the group.
        Can be higher than the number of samples if several original peaks (from peak-picking) were combined.
        If this is a lot higher than the number of samples something is probably wrong.
        If this is a lot lower than the number of samples something is probably wrong.
• A, B, C… (As the groups, or the folders you used, were named): Number of samples where the feature was found (in original peak-picking).
                                   Each file only counted one time.
• Other columns: The intensities/areas in each file.

The other kind of table if the raw peaks as they were detected in each peak.
That includes all peaks that didn't survive the grouping step.

You would usually use this to check if a peak was actually found but removed by grouping.
Also contains extra information for each peak.

peaks <- peaks(xset_g_r_g_fill)

kable(peaks[1:5,])

The sample columns says which sample the peak was found in. Same order as your original input and xset_g_r_g_fill@filepaths.

Troubleshooting

There are many packages that add additional features to XCMS.
Here are some examples:

• CAMERA[4]: Grouping of features and annotation (next talk!)
• X¹³CMS[5]: Unbiased mapping of isotopic fates
• XCMS Online[6]: Online and interactive use of XCMS
• metaXCMS[7]: Finding shared alterations among multiple sample classes
• IPO[8]: Automatic optimization of XCMS parameters
• more exist… See[9]

[1] C. a. Smith, E. J. Want, G. O'Maille, R. Abagyan, et al. “XCMS: Processing mass spectrometry data for metabolite profiling using nonlinear peak alignment, matching, and identification”. In: Analytical Chemistry 78.3 (feb. 2006), pp. 779-787. DOI: 10.1021/ac051437y. .

[2] R. Tautenhahn, C. Böttcher and S. Neumann. “Highly sensitive feature detection for high resolution LC/MS.”. En. In: BMC bioinformatics 9.1 (nov. 2008), p. 504. DOI: 10.1186/1471-2105-9-504. .

[3] C. J. Conley, R. Smith, R. J. O. Torgrip, R. M. Taylor, et al. “Massifquant: Open-source Kalman filter-based XC-MS isotope trace feature detection”. In: Bioinformatics 30.18 (2014), pp. 2636-2643. ISSN: 14602059. DOI: 10.1093/bioinformatics/btu359. .

[4] C. Kuhl, R. Tautenhahn, C. Böttcher, T. R. Larson, et al. “CAMERA: An integrated strategy for compound spectra extraction and annotation of liquid chromatography/mass spectrometry data sets”. In: Analytical Chemistry 84.1 (jan. 2012), pp. 283-289. DOI: 10.1021/ac202450g. eprint: NIHMS150003. .

[5] X. Huang, Y. Chen, K. Cho, I. Nikolskiy, et al. “X13CMS: Global Tracking of Isotopic Labels in Untargeted Metabolomics”. In: Analytical Chemistry 86.3 (2014). PMID: 24397582, pp. 1632-1639. DOI: 10.1021/ac403384n. .

[6] H. Gowda, J. Ivanisevic, C. H. Johnson, M. E. Kurczy, et al. “Interactive XCMS online: Simplifying advanced metabolomic data processing and subsequent statistical analyses”. In: Analytical Chemistry 86.14 (2014), pp. 6931-6939. DOI: 10.1021/ac500734c.

[7] R. Tautenhahn, G. J. Patti, E. Kalisiak, T. Miyamoto, et al. “MetaXCMS: Second-order analysis of untargeted metabolomics data”. In: Analytical Chemistry 83.3 (feb. 2011), pp. 696-700. DOI: 10.1021/ac102980g. .

[8] G. Libiseller, M. Dvorzak, U. Kleb, E. Gander, et al. “IPO: a tool for automated optimization of XCMS parameters”. In: BMC Bioinformatics 16.1 (2015), pp. 1-10. ISSN: 1471-2105. DOI: 10.1186/s12859-015-0562-8. .

[9] N. G. Mahieu, J. L. Genenbacher and G. J. Patti. A roadmap for the XCMS family of software solutions in metabolomics. 2016. DOI: 10.1016/j.cbpa.2015.11.009. .