In PaoloBnn/Retip: Retention Time Prediction for Metabolomics

knitr::opts_chunk$set(echo = TRUE)

Introduction

Retip is a tool for predicting Retention Time (RT) for small molecules in a high pressure liquid chromatography (HPLC) Mass Spectrometry analysis, available as both an R package and a Python package. Retention time calculation can be useful in identifying unknowns and removing false positive annotations. The R package uses six different machine learning algorithms to built a stable, accurate and fast RT prediction model:

Random Forest: a decision tree algorithms.
BRNN: Bayesian Regularized Neural Network.
XGBoost: an extreme Gradient Boosting for tree algorithms.
lightGBM: a gradient boosting framework that uses tree based learning algorithms.
Keras: a high-level neural networks API for Tensorflow.
H2O autoML: an automatic machine learning tool.

Retip also includes useful biochemical databases like: HMDB, KNApSAcK, ChEBI, DrugBank, SMPDB, YMDB, T3DB, FooDB, NANPDB, STOFF, BMDB, LipidMAPS, Urine, Saliva, Feces, ECMDB, CSF, Serum, PubChem.1, PlantCyc, UNPD, BLEXP, NPA and COCONUT.

Get started

To use Retip, you need to prepare a compound retention time library in the format below. The input file needs compound Name, InChIKey, SMILES code and experimental retention time information for each compound. The input must be an MS Excel file. Retip will use this input file to build a the model and will predict retention times for other biochemical databases or an input query list of compounds. It is suggested that the file has at least 300 compounds to build a good retention time prediction model.

Input file view:

| Name | InChIKey | SMILES | RT | |----------------------------------------------|-----------------------------|------------------------------------------------------|----------| | (2-oxo-2,3-dihydro-1H-indol-3-yl)acetic acid | ILGMGHZPXRDCCS-UHFFFAOYSA-N | C1=CC=C2C(=C1)C(C(=O)N2)CC(=O)O | 2.019083 | | 1,1-Dimethyl-4-phenylpiperazinium | WVHNHHJEHFWYHH-UHFFFAOYSA-N | CC1C(NC(CN1)C2=CC=CC=C2)C | 2.60795 | | 1,2-Cyclohexanediol | PFURGBBHAOXLIO-OLQVQODUSA-N | C1CCC\@\@HO | 4.87655 | | 1,2-Cyclohexanedione | OILAIQUEIWYQPH-UHFFFAOYSA-N | C1CCC(=O)C(=O)C1 | 5.772267 | | 1,3,7-Trimethyluric acid | BYXCFUMGEBZDDI-UHFFFAOYSA-N | CN1C2=C(NC1=O)N(C(=O)N(C2=O)C)C | 1.827733 | | 1,3 Cyclohexanedione | HJSLFCCWAKVHIW-UHFFFAOYSA-N | C1CC(=O)CC(=O)C1 | 1.473133 | | 1,4-Cyclohexanedicarboxylic acid | PXGZQGDTEZPERC-UHFFFAOYSA-N | C1CC(CCC1C(=O)O)C(=O)O | 1.560217 |

In case you want to avoid building a library from scratch, you can utilize publicly available libraries from Riken, Japan PlaSMA and the West Coast Metabolomics Center UC Davis, USA MoNA. Experimentation details for these two libraries are available here (link). We have already calculated retention time for several databases for these two experimental conditions.

Packages import

Set keras_installed = TRUE if you have installed keras.

setwd(dirname(rstudioapi::getActiveDocumentContext()$path))
library(Retip)
keras_installed <- FALSE
if (!keras_installed) {
  keras::install_keras()
}

Set up Retip

First of all you have execute function prepare.wizard that is needed to activate the parallel computing to speed up models functions.

Then you have to import your custom excel file or use FiehnLab Hilic or Riken Plasma included library.

#>Starts parallel computing
prep.wizard()

# import excel file for training and testing data
rp2 <- readxl::read_excel("Plasma_positive.xlsx", sheet = "lib_2",
                          col_types = c("text", "text", "text", "numeric"))

#> or use HILIC database included in Retip
#> HILIC <- HILIC

Compute Chemical Descriptors with CDK

Now it is time to compute chemical descriptors with CDK a JAVA based open source project aimed at cheminformatics. Several core people around the Steinbeck group developed the code. According to Ohloh the development took 136 Person, several years and translated into several million dollars in development cost. For you will be simple as type getCD. Depending on your library length and your PC it will take few minutes to several hours. There is a counter to visualize the progress. It’s possible that the return table have less compounds respect the one you have start with. This is because some of yours smile code is not well formatted or for some compounds CDK is not able to calculate chemical descriptors. It is a pity but in real life happens that you have to leave behind even your favorite compound. If happens with most of compounds check SMILES, or better use Pubchem interchange service to get good working SMILES for your library.

Remember: clean your SMILES from salts and metal containing compounds before starting Retip to avoid cavities. You can use for this purpose ChemAxon standardizer, for example.

descs <- getCD(rp2)

Clean dataset

Some machine learning model did not like NA value, or low variance columns. So they are eliminated to improve model performance. This is done in Retip function proc.data.

#> Clean dataset from NA and low variance value
db_rt <- proc.data(descs)

Plot ChemSpace

You have now created your library with chemical descriptors. It is time to visualize your data and answer a two crucial question:

What do you want to predict?
Can I use my library to predict the whole Human Metabolome Database? Or PlantCyc?

The answer is not easy and pass trough another question: is the chemical space of my library large enough? So we have created function chem.space to plot your indulged data into chemical reality. You can chose a target between several database (HMDB, KNAPSACK, CHEBI, DRUGBANK, SMPDB, YMDB, T3DB, FOODB, NANPDB, STOFF, BMDB, LIPIDMAPS, URINE, SALIVA, FECES, ECMDB, CSF, SERUM, PUBCHEM.1, PLANTCYC, UNPD, BLEXP, NPA and COCONUT). In red you will see your library, in gray the chosen database.

chem.space(db_rt, target = "HMDB", title = "HMDB - RIKEN PLASMA")

Center and scale - Optional - Warning

Centering and scaling can be useful for Keras model, but can decrease performance of others models. Actually is a very delicate step that you have to pay a lot of attention. Set cesc=TRUE if you want to perform this step. In that case, cesc input will be included when you predict your target database. It is necessary to avoid obtaining a completely wrong prediction.

Pay attention!

#> Optional center and scale data
#> this option improves Keras model
cesc <- FALSE
if (cesc) {
  # Build a model to use for center and scale a dataframe
  preproc <- cesc(db_rt)
  # use the above created model for center and scale dataframe
  db_rt <- predict(preproc, db_rt)
}

#> IMPORTANT : if you use this option remember to set cesc variable into
#> rt.spell and getRT.smile functions

Training and Testing set-up

A very important step is to split in training and testing your data frame.

Guide lines to build effective model explicitly specify: “QSAR/QSPR model have to be properly validated using data which was not in the training set.”

This is not new, is known from the Ancient Rome “divide et impera” (divide and conquer).

We can use caret function createDataPartition as follow:

#> Split in training and testing using caret::createDataPartition
set.seed(101)
in_training <- caret::createDataPartition(db_rt$XLogP, p = .8, list = FALSE)
training <- db_rt[in_training, ]
testing <- db_rt[-in_training, ]

Building QSAR models

Now is time to put your hands on Retip ALMA mater: Advanced Learning Machine Algorithms.

Our suggestion is to compute all of them. Each function have an integrated tuning function to find out best parameters on your own library. Do not worry too much about over fitting because we have done 10 times fold cross validation. Each model have his advantage or disadvantage, it is simply depends on your data the one that fit better. Your models will be saved after building them so you can directly load them the next time that you open your project. Set buildModels=FALSE to load your previously built models.

Have fun!

build_models <- TRUE

if (build_models) {
  rf  <- fit.rf(training)
  saveRDS(rf, "rf_model.rds")

  brnn <- fit.brnn(training)
  saveRDS(brnn, "brnn_model.rds")

  keras <- fit.keras(training, testing)
  save_model_hdf5(keras, filepath = "keras_model.h5")

  lightgbm <- fit.lightgbm(training, testing)
  saveRDS(lightgbm, "lightgbm_model.rds")

  xgb <- fit.xgboost(training)
  saveRDS(xgb, "xgb_model.rds")

  aml <- fit.automl.h2o(training)
  # It saves by default the best model
  h2o::h2o.saveModel(aml@leader, "automl_h2o_model")
} else {
  xgb <- readRDS("xgb_model.rds")
  rf <- readRDS("rf_model.rds")
  brnn <- readRDS("brnn_model.rds")
  keras <- keras::load_model_hdf5("keras_model.h5")
  lightgbm <- readRDS("lightgbm_model.rds")
  h2o::h2o.init(nthreads = -1, strict_version_check=FALSE)
  aml <- h2o.loadModel("automl_h2o_model/###") # replace ### with the name of your saved model
}

Testing and plot models

Congratulations you have built your model!

Let’s see how much accurate it is. First of all check the scores with function get.score. Returns:

RMSE: root-mean-square error (the most important one)
R2: coefficient of determination
MAE: Mean Absolute Error
95% ± min: is the 95% confidence in minutes

The best model is always the first based on lower RMSE value. You can check only one model or all of them (maximum 6) at same time.

#> first you have to put the testing dataframe and then the name of the models
#> you have computed
stat <- get.score(testing, xgb, rf, brnn, keras, lightgbm, aml)
stat

You can also visualize model performance using function p.model using testing dataframe as input. You can also use training data or whole data to see how it change, even if the only really important is the testing one.

XGBoost

#> first you have to put the testing daframe and then the name of the models
#> you want to visualize.
#> Last value is the title of the graphics.
p.model(testing, m = xgb, title = "XGBoost - RIKEN PLASMA")

Random forest

#> first you have to put the testing daframe and then the name of the models
#> you want to visualize.
#> Last value is the title of the graphics.
p.model(testing, m = rf, title = "Random forest - RIKEN PLASMA")

BRNN

#> first you have to put the testing daframe and then the name of the models
#> you want to visualize.
#> Last value is the title of the graphics.
p.model(testing, m = brnn, title = "BRNN - RIKEN PLASMA")

LightGBM

#> first you have to put the testing daframe and then the name of the models
#> you want to visualize.
#> Last value is the title of the graphics.
p.model(testing, m = lightgbm, title = "LightGBM - RIKEN PLASMA")

Keras

#> first you have to put the testing daframe and then the name of the models
#> you want to visualize.
#> Last value is the title of the graphics.
p.model(testing, m = keras, title = "Keras - RIKEN PLASMA")

H2O autoML

#> first you have to put the testing daframe and then the name of the models
#> you want to visualize.
#> Last value is the title of the graphics.
p.model(testing, m = aml, title = "H2O autoML - RIKEN PLASMA")

For the H2O autoML model, the statistics and the figures are obtained from the best model (aml@leader). If you want to check one of the other models, you have to select it using the model ID. The same applies when saving your H2O AutoML model. You will find an example of how to get a different model in the following command block.

#> First, get the leaderboard of your model
# aml_lb <- as.data.frame(h2o::h2o.get_leaderboard(object = aml, extra_columns = "ALL"))
#> Second, decide the model you want to inspect and select
#> The number corresponds to the index
# model_id <- aml_lb[3, "model_id"]
# selected_model <- h2o::h2o.getModel(model_id)

Feature importance

You can visualize the feature importance of your model using the p.model.features function. You have to indicate your model (model) and your model type (mdl_type). The accepted inputs for mdl_type are XGBoost (xgb), Random Forest (rf), BRNN (brnn), LightGBM (lightgbm) and H2O autoML (h2o). The Keras model is not supported.

XGBoost

p.model.features(xgb, mdl_type = "xgb")

Random forest

p.model.features(rf, mdl_type = "rf")

BRNN

p.model.features(brnn, mdl_type = "brnn")

LightGBM

p.model.features(lightgbm, mdl_type = "lightgbm")

H2O autoML

It is not possible to plot the stacked ensemble models, but it is possible to plot the best 20 models used to construct the stacked one. For this reason, you have to indicate the number of the model (h2o_number) you want to see. You can find the correspondence between model ID and number running the following command line.

# aml_model_number <- as.data.frame(rownames(h2o::h2o.varimp(aml)))

p.model.features(aml, mdl_type = "h2o", h2o_number = 1)

Retention Time Prediction SPELL

This is the the final step. All the effort made before have been accomplished to compute the prediction on a whole database.

You have 3 options to make a prediction:

for your personal database
for downloaded mona.msp
for a Retip included database

1. Personal database

In that case, you have to import your excel with compound Name, InChIKey and SMILES code, compute CD and then make the prediction.

# import excel file for external validation set
rp_ext <- readxl::read_excel("Plasma_positive.xlsx", sheet = "ext",
                             col_types = c("text", "text", "text", "numeric"))
#> compute Chemical descriptors
rp_ext_desc <- getCD(rp_ext)

Clean data set from NA and low variance values using the Retip function proc.data as you did before with your data.

#> Clean dataset from NA and low variance value
db_rt_ext <- proc.data(rp_ext_desc)

Center and scale your data if you performed this step before.

if (cesc) {
  # Build a model to use for center and scale a dataframe
  preproc_ext <- cesc(db_rt_ext)
  # use the above created model for center and scale dataframe
  db_rt_ext <- predict(preproc_ext, db_rt_ext)
}

XGBoost

#> perform the RT spell
if (cesc) {
  rp_ext_pred_xgb <- Retip::RT.spell(training, rp_ext_desc, model = xgb,
                                     cesc = preproc)
} else {
  rp_ext_pred_xgb <- Retip::RT.spell(training, rp_ext_desc, model = xgb)
}

Random forest

#> perform the RT spell
if (cesc) {
  rp_ext_pred_rf <- Retip::RT.spell(training, rp_ext_desc, model = rf,
                                    cesc = preproc)
} else {
  rp_ext_pred_rf <- Retip::RT.spell(training, rp_ext_desc, model = rf)
}

BRNN

#> perform the RT spell
if (cesc) {
  rp_ext_pred_brnn <- Retip::RT.spell(training, rp_ext_desc, model = brnn,
                                      cesc = preproc)
} else {
  rp_ext_pred_brnn <- Retip::RT.spell(training, rp_ext_desc, model = brnn)
}

LightGBM

#> perform the RT spell
if (cesc) {
  rp_ext_pred_lightgbm <- RT.spell(training, rp_ext_desc, model = lightgbm,
                                   cesc = preproc)
} else {
  rp_ext_pred_lightgbm <- RT.spell(training, rp_ext_desc, model = lightgbm)
}

Keras

#> perform the RT spell
if (cesc) {
  rp_ext_pred_keras <- RT.spell(training, rp_ext_desc, model = keras,
                                cesc = preproc)
} else {
  rp_ext_pred_keras <- RT.spell(training, rp_ext_desc, model = keras)
}

H2O autoML

#> perform the RT spell
if (cesc) {
  rp_ext_pred_aml <- RT.spell(training, rp_ext_desc, model = aml,
                                cesc = preproc)
} else {
  rp_ext_pred_aml <- RT.spell(training, rp_ext_desc, model = aml)
}

2. Downloaded mona.msp

Here, you have to put your msp downloaded from Mona inside your project folder. First, prep.mona function is needed to extract information from msp file and build input table. Then, you have to compute CD using getCD, make the prediction and incorporate again inside msp file. In this way the msp file you have in your project folder is ready to be used for compound identification with MSMS and Retention time prediction.

db_mona <- prep.mona(msp = "MoNA-export-CASMI_2016.msp")

mona <- getCD(db_mona)

XGBoost

if (cesc) {
  mona_rt_xgb <- RT.spell(training, mona, model = xgb, cesc = preproc)
} else {
  mona_rt_xgb <- RT.spell(training, mona, model = xgb)
}

addRT.mona(msp = "MoNA-export-CASMI_2016.msp", mona_rt_xgb)

Random forest

if (cesc) {
  mona_rt_rf <- RT.spell(training, mona, model = rf, cesc = preproc)
} else {
  mona_rt_rf <- RT.spell(training, mona, model = rf)
}

addRT.mona(msp = "MoNA-export-CASMI_2016.msp", mona_rt_rf)

BRNN

if (cesc) {
  mona_rt_brnn <- RT.spell(training, mona, model = brnn,
                           cesc = preproc)
} else {
  mona_rt_brnn <- RT.spell(training, mona, model = brnn)
}

addRT.mona(msp = "MoNA-export-CASMI_2016.msp", mona_rt_brnn)

LightGBM

if (cesc) {
  mona_rt_lightgbm <- RT.spell(training, mona, model = lightgbm,
                               cesc = preproc)
} else {
  mona_rt_lightgbm <- RT.spell(training, mona, model = lightgbm)
}

addRT.mona(msp = "MoNA-export-CASMI_2016.msp", mona_rt_lightgbm)

Keras

if (cesc) {
  mona_rt_keras <- RT.spell(training, mona, model = keras,
                            cesc = preproc)
} else {
  mona_rt_keras <- RT.spell(training, mona, model = keras)
}

addRT.mona(msp = "MoNA-export-CASMI_2016.msp", mona_rt_keras)

H2O autoML

if (cesc) {
  mona_rt_aml <- RT.spell(training, mona, model = aml,
                            cesc = preproc)
} else {
  mona_rt_aml <- RT.spell(training, mona, model = aml)
}

addRT.mona(msp = "MoNA-export-CASMI_2016.msp", mona_rt_aml)

3. Retip included database

This option is the easiest one. We computed more than 300k compounds to save your time and more important save electricity (we love nature). You have just to chose a target database, the type of output and the polarity (pos or neg).

The available outputs are:

MSDIAL: identification with accurate Mass and RTP (RECOMMENDED); get MSDIAL
MSFINDER: identification with accurate Mass, in silico fragmentation MS2 and RTP (RECOMMENDED); get MSFINDER
AGILENT: works with Mass Hunter software
THERMO
WATERS
SCIEX

XGBoost

#> example of all Retip included compounds predicted and exported to MSFINDER
if (cesc) {
  all_pred_xgb <- RT.spell(training, target = "ALL", model = xgb,
                           cesc = preproc)
} else {
  all_pred_xgb <- RT.spell(training, target = "ALL", model = xgb)
}

export_rtp <- RT.export(all_pred_xgb, program = "MSFINDER", pol = "pos")

Random forest

#> example of all Retip included compounds predicted and exported to MSFINDER
if (cesc) {
  all_pred_rf <- RT.spell(training, target = "ALL", model = rf,
                          cesc = preproc)
} else {
  all_pred_rf <- RT.spell(training, target = "ALL", model = rf)
}

export_rtp <- RT.export(all_pred_rf, program = "MSFINDER", pol = "pos")

BRNN

#> example of all Retip included compounds predicted and exported to MSFINDER
if (cesc) {
  all_pred_brnn <- RT.spell(training, target = "ALL", model = brnn,
                            cesc = preproc)
} else {
  all_pred_brnn <- RT.spell(training, target = "ALL", model = brnn)
}

export_rtp <- RT.export(all_pred_brnn, program = "MSFINDER", pol = "pos")

LightGBM

#> example of all Retip included compounds predicted and exported to MSFINDER
if (cesc) {
  all_pred_lightgbm <- RT.spell(training, target = "ALL", model = lightgbm,
                                cesc = preproc)
} else {
  all_pred_lightgbm <- RT.spell(training, target = "ALL", model = lightgbm)
}

export_rtp <- RT.export(all_pred_lightgbm, program = "MSFINDER", pol = "pos")

Keras

#> example of all Retip included compounds predicted and exported to MSFINDER
if (cesc) {
  all_pred_keras <- RT.spell(training, target = "ALL", model = keras,
                             cesc = preproc)
} else {
  all_pred_keras <- RT.spell(training, target = "ALL", model = keras)
}

export_rtp <- RT.export(all_pred_keras, program = "MSFINDER", pol = "pos")

H2O autoML

#> example of all Retip included compounds predicted and exported to MSFINDER
if (cesc) {
  all_pred_aml <- RT.spell(training, target = "ALL", model = aml,
                             cesc = preproc)
} else {
  all_pred_aml <- RT.spell(training, target = "ALL", model = aml)
}

export_rtp <- RT.export(all_pred_aml, program = "MSFINDER", pol = "pos")

Conclusions remarks

Remember these few tips:

Split your data into training and testing (80/20). If you want to be cool use also an additional external dataset.
Do not cheat yourself doing cherry picking with testing molecules that are out-liner and putting in training data. Leave at it is, is better know the truth than lie to yourself.
As you probably know R use set.seed when you have a random function. This is needed to get the same results when you do it again. There is a set seed before the split training/testing. If you modify you get a slightly different results in your models if the problematic compounds are inside training data. This is not a real cheat because is random driven ;-)
Look at your smiles before import in Retip, if you have cavities in it will not work properly.

Remember Leonardo Da Vinci suggestion:

{css style settings, echo = FALSE} blockquote { padding: 10px 20px; margin: 0 0 20px; font-size: 14px; border-left: 5px solid #eee; }

I have been impressed with the urgency of doing. Knowing is not enough; we must apply. Being willing is not enough; we must do.

Retip born to helps you in identification in metabolomics workflow, it is not an exhibition of knowledge. The success of this app is if helps in your work, in real world metabolomics experiments.

PaoloBnn/Retip documentation built on June 8, 2024, 6:19 a.m.

rdrr.io home R language documentation Run R code online

CRAN packages Bioconductor packages R-Forge packages GitHub packages

Note that we can't provide technical support on individual packages. You should contact the package authors for that.

PaoloBnn/Retip Retention Time Prediction for Metabolomics

In PaoloBnn/Retip: Retention Time Prediction for Metabolomics

Introduction

Get started

Packages import

Set up Retip

Compute Chemical Descriptors with CDK

Clean dataset

Plot ChemSpace

Center and scale - Optional - Warning

Training and Testing set-up

Building QSAR models

Testing and plot models

XGBoost

Random forest

BRNN

LightGBM

Keras

H2O autoML

Feature importance

XGBoost

Random forest

BRNN

LightGBM

H2O autoML

Retention Time Prediction SPELL

1. Personal database

XGBoost

Random forest

BRNN

LightGBM

Keras

H2O autoML

2. Downloaded mona.msp

XGBoost

Random forest

BRNN

LightGBM

Keras

H2O autoML

3. Retip included database

XGBoost

Random forest

BRNN

LightGBM

Keras

H2O autoML

Conclusions remarks

R Package Documentation

Browse R Packages

We want your feedback!

PaoloBnn/Retip
Retention Time Prediction for Metabolomics