Annotating LC-MS data with PredRet in R

In this tutorial I will show you how you can use the PredRet database to annotate your LC-MS metabolomics data directly in R. The annotation function I will be using is general-purpose, though, so you can use it to annotate any data where you have a list of compounds with known m/z’s and retention times (RTs).

In short PredRet is a user-driven database of compound retention times. The purpose of PredRet is to be able to predict the RT of a compound in one (your!) chromatographic system if it has been experimentally determined in another chromatographic system by someone, somewhere in the world. You can download the paper here and visit the project’s home page at PredRet.org.

Pulling data from PredRet

So lets get going.
First we will download the PredRet database with the PredRetR package. We will get both the experimental and the predicted values for the chromatographic system “LIFE_old”.

library(PredRetR)

predret_db <-    PredRet_get_db(exp_pred = c("exp","pred"),
                                systems = "LIFE_old"
                               )

Then lets take a look at the structure of the data.frame we have retrieved. If the recorded_rt column has data we have an experimentally determined RT, if the predicted_rt column has data we have a RT predicted with the PredRet systems. The ci_lower and ci_upper columns show the prediction interval for the predictions.

dim(predret_db)
head(predret_db)

## [1] 573  12

We have the RT directly in the database above but we do not have the m/z. Since we have the InChI the easiest way to get the m/z is to extract the molecular formula and then use the Rdisop package to get the mass.

library(Rdisop)
library(stringr)

formulas   <- str_split_fixed(predret_db[,"inchi"],"/",3)[,2]
masses     <- sapply(formulas, function(x) getMass(getMolecule(x))   )

predret_db <- cbind.data.frame(predret_db, mz_pos = masses + 1.0073)

We can then split the database in two. One for the experimental RTs and one for the predicted. Here I have used some pipe/dplyr style code and even a forward assign. If you are not familiar with this type of R code I urge you to look into it. It really makes code much more readable.

library(dplyr)
predret_db %>% filter(!is.na(recorded_rt)) -> predret_db_exp
predret_db %>% filter(is.na(recorded_rt))  -> predret_db_pred

Annotating a dataset

We now have the database ready for annotation.
So we can load a dataset/peaklist. This peaklist was previously created with XCMS and fragments/adducts annotated with CAMERA. But again any peaklist will do.

library(readxl)

data <- read_excel("../material/data/peaklist_pos.xlsx") %>% 
        replace(is.na(.), "") %>% 
        as.data.frame

Lets take a look at the interesting columns we have in the dataset.

info_cols <- c("mz","rt","isotopes","adduct","pcgroup")
dim(data)

## [1] 4232  470

head(    data[,info_cols]    )

##          mz       rt   isotopes              adduct pcgroup
## 1  62.06013 273.8827            [M+H-C5H8O]+ 145.11       1
## 2  95.08600 271.5567                                      1
## 3  98.09642 268.1695                                      1
## 4 104.10687 271.2400   [12][M]+ [M+H-COCH2]+ 145.11       1
## 5 105.11030 271.2415 [12][M+1]+                           1
## 6 114.09119 268.9771            [M+H-CH3OH]+ 145.11       1

Now we can use db.comp.assign from my chemhelper package to annotate the dataset.
We would probably want one column in our peaklist for the RTs we have determined experimentally and one for the predicted RTs. So we do the annotation twice.
The first two arguments to the function are the m/z and RT of the dataset. The next three are the name, m/z and RT of the database of known compounds. Lastly we give the tolerance for the m/z and RT for a database match.

library(chemhelper)

anno_exp <- db.comp.assign(mz           = data[,"mz"],
                           rt           = data[,"rt"],
                           comp_name_db = predret_db_exp[,"name"],
                           mz_db        = predret_db_exp[,"mz_pos"],
                           rt_db        = predret_db_exp[,"recorded_rt"] *60, # *60 since XCMS works in seconds but PredRet is in minutes.
                           mzabs        = 0.01,
                           ppm          = 15,
                           ret_tol      = 10
                           )

##                                         [,1]
## Unique hits                               86
## Non-unique hits                           10
## Compounds not found                      207
## Markers had multiple compounds assigned   50

anno_pred <- db.comp.assign(mz          = data[,"mz"],
                           rt           = data[,"rt"],
                           comp_name_db = predret_db_pred[,"name"],
                           mz_db        = predret_db_pred[,"mz_pos"],
                           rt_db        = predret_db_pred[,"predicted_rt"] *60,
                           mzabs        = 0.01,
                           ppm          = 15,
                           ret_tol      = 10
                           )

##                                         [,1]
## Unique hits                               34
## Non-unique hits                            2
## Compounds not found                      234
## Markers had multiple compounds assigned    5

Now lets put the annotations together with our dataset.

data_annotated <- cbind.data.frame(data,anno_exp,anno_pred)

The result

Then lets take a look at one of the feature groups (from CAMERA) where we got an annotation. In this example we have a feature that was annotated as tryptophan. The “OR” is because tryptophan was in the database several times. If multiple compounds would fit the m/z and RT they would be written with “OR” between them.
There is also a fragment annotated as adipic acid and 2-methylglutaric acid using a hit from the predicted RTs. In this case we have almost perfect CAMERA annotation suggesting the pseudo-molecular ion is the m/z = 205.0979 feature and the compound is very likely tryptophan.

data_annotated %>% filter(pcgroup==16) %>% select(one_of(info_cols,"anno_exp","anno_pred"))

Lets take a look at another feature. This time we have experimental RTs that say the feature is either 1,7-dimethylxanthine OR theobromine but a prediction from the PredRet database also suggest that the feature could be theophylline as well.

data_annotated %>% filter(pcgroup==400) %>% select(one_of(info_cols,"anno_exp","anno_pred"))

That’s it!