chromatographR

Workflow

Loading data

chromatographR now includes an option to directly convert files from ChemStation or MassHunter (Bovee 2022a, 2022b). Alternatively, the data from the spectrometer can be converted into an R-friendly format like a csv file. The load_chroms function provides a simple interface to read chromatograms into R from any of these formats.

> # single folder
> load_chroms(paths = path, format.in="csv")
> 
> # multiple folders
> path = 'foo'
> folders <- list.files(path = path, pattern = "EXPORT3D")
> dat <- load_chroms(folders)

Example data

We have included some example data consisting of root extracts from tall goldenrod (Solidago altissima). Roots were extracted in 90% methanol and run on an Agilent 1100 HPLC coupled with a DAD detector, according to a previously described method (Uesugi and Kessler 2013). The dataset is called (short for Solidago altissima).

> data(Sa)

Pre-processing data

Data from liquid chromatography often suffer from a variety of non-informative artifacts, such as noise or a drifting baseline. In addition, the data produced by the instrument may have a higher resolution or wider range (along either the time or spectral dimensions) than we require. Fortunately, most of these issues can be remedied fairly easily. For example, smoothing can reduce noise in the spectral direction while baseline subtraction can help correct a drifting baseline. Interpolation of wavelengths and retention times can be used to reduce the dimensionality of the data, facilitating comparison across samples and reducing the time and computational load required by downstream analyses. All of these functions (smoothing, baseline correction, and interpolation) are available through the preprocess function and are enabled by default.

To select a narrower range of times and/or wavelengths, arguments can be provided to dim1 and dim2 arguments. The baseline_cor function from the ptw package (Bloemberg et al. 2010) takes arguments p (an asymmetry parameter) and lambda (a smoothing parameter). You can read more about these in the documentation for ptw::asysm. You may want to experiment with these parameters before choosing values to use on your whole dataset.

> i=2 # chromatogram number in list of data
> tpoints <- as.numeric(rownames(Sa[[i]]))
> lambda='200.00000'
> 
> matplot(x=tpoints, y=Sa[[i]][,lambda],
+        type='l', ylab='Abs (mAU)', xlab='Time (min)')
> matplot(x=tpoints, y = ptw::baseline.corr(Sa[[i]][,lambda],p=.001,lambda=1e5),
+        type='l', add = T, col='blue', lty = 3)

> new.ts <- seq(10,18.66,by=.01) # choose time-points
> new.lambdas <- seq(200, 318, by = 2) # choose wavelengths

> dat.pr <- preprocess(Sa, dim1=new.ts, dim2=new.lambdas,
+            parallel=F, p=.001, lambda=1e5)

Dealing with retention time shifts

In many cases, liquid chromatography can suffer from retention time shifts (e.g. due to temperature fluctuations, column degradation, or subtle changes in mobile-phase composition), which can make it very difficult to compare peaks across samples. Luckily, a number of “time-warping” algorithms have been developed for correcting these kinds of shifts. In chromatographR, parametric time warping (ptw) (Eilers 2004; Bloemberg et al. 2010) and variable penalty dynamic time warping (VPdtw) (Clifford et al. 2009; Clifford and Stone 2012) are available for correcting retention time shifts through the correct_rt function. Both warping functions aim to produce a better alignment of retention times by “warping” the time-axis of each supplied chromatogram to match a reference chromatogram. (The reference chromatogram can be determined algorithmically or selected manually by setting the reference argument).

Parametric time warping

The ptw option can take a single wavelength or a list of wavelengths provided by the user using the lambdas argument. For each chromatogram, ptw then produces a “global” warping function across all the wavelengths supplied by the user. The code block below creates warping models for the samples in the provided list of data matrices. The same function is then used to warp each chromatogram according to the corresponding model, by setting the models parameter. Depending on the variety of your samples and the severity of the retention time shifts, it may take some experimentation with the warping parameters to get satisfactory results. Sometimes “less” can actually be “more” here – for example, wavelengths with fewer peaks may sometimes yield better warping models. (Also see the documentation for ptw for more guidance on warp function optimization).

> warping.models <- correct_rt(dat.pr, what = "models", lambdas=c("210"), scale=TRUE)
> warp <- correct_rt(chrom_list=dat.pr, models=warping.models, what="corrected.values")

We can then use the following code snippet to compare the alignment of warped (top panel) and unwarped (bottom panel) chromatograms.

> par(mfrow=c(2,1))
> lambdas=c('210','260')
> plot.new()
> ts <- as.numeric(rownames(warp[[i]]))
> plot.window(xlim=c(head(ts,1), tail(ts,1)),ylim=c(0,1000))
> for (i in 1:length(warp)){
+   matplot(ts, warp[[i]][,lambdas],type='l',add=T)
+ }
> legend("topright", legend="ptw", bty = "n")
> 
> plot.new()
> ts <- as.numeric(rownames(dat.pr[[i]]))
> plot.window(xlim=c(head(ts,1),tail(ts,1)),ylim=c(0,1000))
> for (i in 1:length(dat.pr)){
+   matplot(ts, dat.pr[[i]][,lambdas],type='l',add=T)
+ }
> legend("topright", legend="raw", bty = "n")

Clearly, the alignment is considerably improved after warping. You can also use the correct_rt function to do a global alignment on multiple wavelengths, by providing a list of wavelengths to the lambdas argument, but this will not always improve the results. While the alignment still isn’t perfect after warping, it is probably good enough to align our peaks and assemble them in the peak table, which is our primary goal.

Variable penalty dynamic time warping

Variable penalty dynamic time warping is another algorithm that can be quite effective for correcting retention time shifts (Clifford et al. 2009; Clifford and Stone 2012). Instead of calculating an explicit warping function, dynamic time warping (DTW) works by repeatedly shifting, expanding or contracting the time axis in small steps until a satisfactory alignment is attained. This process can result in very close matching of peaks, but can also result in severe peak distortion if the time axis is repeatedly expanded and/or contracted (Tomasi, van den Berg, and Andersson 2004). Variable penalty dynamic time warping is a derivative of DTW that constrains expansion and contraction of the time-axis by applying a penalty to these types of risky moves. This can often achieve a very good alignment with minimal peak distortion. VPdtw package must be installed manually, since it is not currently available on CRAN. Unlike , can only take a single wavelength as input (but it returns a global warping on the whole chromatographic matrix). The code snippet below is not run, but shows how one could install VPdtw and use it to warp the example chromatograms.

> install.packages("VPdtw", repos="https://ethanbass.github.io/drat")
> warp <- correct_rt(chrom_list=dat.pr, alg="vpdtw", lambda="210", what="corrected.values")

Peak finding and fitting

The get_peaks function produced a nested list of peaks by looping through the supplied chromatograms at the specified wavelengths, finding peaks, and fitting them to the specified function using non-linear least squares. The current peak-shape options are gaussian and egh (exponential-gaussian hybrid) (Lan and Jorgenson 2001). The area under the curve for each peak is then estimated using trapezoidal approximation. The function returns a nested list of data.frames containing parameters for the peaks identified in each chromatogram.

> pks_gauss <- get_peaks(warp, lambdas = c('210','260'), sd.max=40, fit="gaussian")
> pks_egh <- get_peaks(warp, lambdas = c('210', '260'), sd.max=40, fit="egh")

The plot.peak_list function allows you to visually assess the peak integration accuracy. Below we compare the peaks fitted for the first chromatogram (index=1) using the two algorithms. Usually the egh function performs slightly better for asymmetrical peaks.

> par(mfrow=c(2,1))
> plot(pks_gauss, index=1, lambda='210')
> plot(pks_egh, index=1, lambda='210')

After obtaining a peak_list, the get_peaks function performs complete linkage hierarchical clustering to link peaks across samples. It returns a peak_table object with samples as rows and peaks as columns. The peak_table object also has slots for holding various metadata about peaks, samples, and analyses.

> pk_tab <- get_peaktable(pks_egh, response = "area")
> head(pk_tab$tab[,1:6])

          V1        V2        V3       V4        V5        V6
119 547.5999   0.00000  64.73611 188.9210 1514.4775  31.08266
121 426.2768 109.41736  61.30032   0.0000  800.6312   0.00000
122 839.5934  94.67302 199.57845   0.0000 2595.0200 179.59510
458 620.1729 228.73095 198.89770 530.8515 1399.5000  28.75473

Further analysis and data-visualization

Attaching metadata

To begin analyzing your peak table, you will usually want to attach sample meta-data to the peak_table object. This can be easily accomplished using the attach_metadata function. The peak table can then be normalized (e.g. by dividing out the sample weight) using the normalize_data function.

> path <- system.file("extdata", "Sa_metadata.csv", package = "chromatographR")
> meta <- read.csv(path)
> pk_tab <- attach_metadata(peak_table = pk_tab, metadata = meta, column="vial")
> pk_tab <- normalize_data(peak_table = pk_tab, column="mass")

Data exploration

Mirror plot

The mirror_plot function provides a quick way to visually compare results across treatment groups.

> mirror_plot(pk_tab, lambdas = c("210","260"), var = "trt", legend_size=2)

Plotting spectra

The plot_spectrum function allows you to easily plot or record the spectra associated with a particular peak in your peak table. This is useful for interpreting your results and/or checking for errors in your peak table. For example, you may want to check if the spectra for a particular peak match across different samples, or you may want to compare your spectrum with a known standard. The plot_spectrum function can be used to plot only the spectrum or only the chromatographic trace using the arguments plot_spectrum and plot_trace. By default it will plot the trace and spectrum from the chromatogram with the largest peak in the peak table. Alternatively, you can choose the chromatogram index and wavelength using the chr and lambda arguments.

> oldpar <- par(no.readonly = TRUE)
> par(mfrow=c(2,1))
> peak="V7"
> plot_spectrum(peak, peak_table = pk_tab, chrom_list=warp, verbose=F)

> par(oldpar)

The plot_all_spectra function can be used to visually compare the spectra for a specified peak across all samples.

> peak="V13"
> plot_all_spectra(peak, peak_table=pk_tab, export=F, overlapping=T)

Plot peak table function

The plot.peak_table function provides a simplified interface to various options for plotting data from the peak table. For example, it can be used as a quick interface to the and functions shown above. It can also be used to make a quick box plot to compare results across treatments.

> plot(pk_tab, loc = "V13", box_plot = TRUE, vars = "trt", verbose = FALSE)