For chemoinformatics research, Rcpi calculates various types of molecular descriptors/fingerprints, and computes pairwise similarities derived by fingerprints and maximum common substructure search. Rcpi also provides the searchDrug() function for parallelized molecular similarity search based on these similarity types.
Regression Modeling in QSRR Study of Retention Indices
In Yan et al. (2012), a quantitative structure-retention relationship study was performed for 656 flavor compounds on four stationary phases of different polarities, using constitutional, topological and geometrical molecular descriptors. The gas chromatographic retention indices (RIs) of these compounds were accurately predicted using linear models. Here we choose the molecules and their RIs of one stationary phase (OV101) as our benchmark dataset.
Since it would be rather tedious to implement the complete cross-validation procedures, the R package caret is used here. To run the R code below, users need to install the package and the required predictive modeling packages first. The package provides a unified interface for the modeling tuning task across different statistical machine learning packages. It is particularly helpful in QSAR modeling, for it contains tools for data splitting, pre-processing, feature selection, model tuning and other functionalities.
Just like the last section, we load the Rcpi package, and read the molecules stored in a SMILES file:
library("Rcpi")
RI.smi = system.file(
"vignettedata/FDAMDD.smi", package = "Rcpi")
RI.csv = system.file(
"vignettedata/RI.csv", package = "Rcpi")
x.mol = readMolFromSmi(RI.smi, type = "mol")
x.tab = read.table(RI.csv, sep = "\t", header = TRUE)
y = x.tab$RI
The readMolFromSmi() function is used for reading molecules from SMILES files, for molecules stored in SDF files, use readMolFromSDF() instead.
The CSV file RI.csv contains tabular data for the retention indices, compound name, and odor information of the compounds. Here we only extracted the RI values by calling x.tab$RI.
After the molecules were properly loaded, we calculate several selected molecular descriptors. The corresponding functions for molecular descriptor calculation are all named after extractDrug...() in Rcpi:
# calculate selected molecular descriptors
x = suppressWarnings(cbind(
extractDrugALOGP(x.mol),
extractDrugApol(x.mol),
extractDrugECI(x.mol),
extractDrugTPSA(x.mol),
extractDrugWeight(x.mol),
extractDrugWienerNumbers(x.mol),
extractDrugZagrebIndex(x.mol)))
After the descriptors were calculated, the result x will be a data frame, each row represents one molecule, and each column is one descriptor (predictor). The Rcpi package integrated the molecular descriptors and chemical fingerprints calculated by the rcdk package Steinbeck et al. (2003) and the ChemmineOB package Horan and Girke (2013).
Next, a partial least squares model will be fitted with the pls and the caret package. The cross-validation setting is 5-fold repeated cross-validation (repeat for 10 times).
# regression on training set
library("caret")
library("pls")
# cross-validation settings
ctrl = trainControl(
method = "repeatedcv", number = 5, repeats = 10,
summaryFunction = defaultSummary)
# train a pls model
set.seed(1002)
pls.fit = train(
x, y, method = "pls", tuneLength = 10, trControl = ctrl,
metric = "RMSE", preProc = c("center", "scale"))
# print cross-validation result
print(pls.fit)
The cross-validation result is:
## Partial Least Squares
##
## 297 samples
## 10 predictors
##
## Pre-processing: centered, scaled
## Resampling: Cross-Validated (5 fold, repeated 10 times)
##
## Summary of sample sizes: 237, 237, 237, 238, 239, 238, ...
##
## Resampling results across tuning parameters:
##
## ncomp RMSE Rsquared RMSE SD Rsquared SD
## 1 104 0.884 9.44 0.0285
## 2 86.4 0.92 6.99 0.0194
## 3 83.8 0.924 6.56 0.0185
## 4 79.6 0.931 6.98 0.0194
## 5 76.3 0.937 7.45 0.0187
## 6 74.7 0.94 6.85 0.0162
## 7 73.7 0.941 6.75 0.0159
## 8 73.5 0.942 6.5 0.0142
## 9 72.5 0.944 6.18 0.0137
##
## RMSE was used to select the optimal model using the smallest value.
## The final value used for the model was ncomp = 9.
We see that the RMSE of the PLS regression model was decreasing when the number of principal components (ncomp) was increasing. We can plot the components and RMSE to helps us select the desired number of principal components used in the model.
# # of components vs RMSE
print(plot(pls.fit, asp = 0.5))
Figure 2: Number of principal components vs. RMSE for the PLS regression model
From Figure 2, we consider that selecting six or seven components is acceptable. At last, we plot the experimental RIs and the predicted RIs to see if the model fits well on the training set (Figure 3):
# plot experimental RIs vs predicted RIs
plot(y, predict(pls.fit, x), xlim = range(y), ylim = range(y),
xlab = "Experimental RIs", ylab = "Predicted RIs")
abline(a = 0, b = 1)
Figure 3: Experimental RIs vs. Predicted RIs
In Silico Toxicity Classification for Drug Discovery
In the perspective of quantitative pharmacology, the successful discovery of novel drugs depends on the pharmacokinetics properties, like absorption, distribution, metabolism, and excretion. In addition, the potential toxicity of chemical compounds is taken into account. QSAR or QSPR methods are usually employed to predict the ADME/T qualities of potential drug candidates.
In D-S Cao et al. (2012), quantitative structure-toxicity relationship (QSTR) models were established for classifying five toxicity datasets. Here we use the maximum recommended daily dose dataset (FDAMDD) from FDA Center for Drug Evaluation and Research as the benchmark dataset.
First, load the drug molecules stored in a SMILES file into R:
library("Rcpi")
fdamdd.smi = system.file("vignettedata/FDAMDD.smi", package = "Rcpi")
fdamdd.csv = system.file("vignettedata/FDAMDD.csv", package = "Rcpi")
x.mol = readMolFromSmi(fdamdd.smi, type = "mol")
x.smi = readMolFromSmi(fdamdd.smi, type = "text")
y = as.factor(paste0("class", scan(fdamdd.csv)))
The object x.mol is used for computing the MACCS and E-state fingerprints; the object x.smi is used for computing the FP4 fingerprints. The 0-1 class labels stored in file FDAMDD.csv indicates whether the drug molecule has high toxicity or not.
Then we calculate three different types of molecular fingerprints (E-state, MACCS, and FP4) for the drug molecules:
# calculate molecular fingerprints
x1 = extractDrugEstateComplete(x.mol)
x2 = extractDrugMACCSComplete(x.mol)
x3 = extractDrugOBFP4(x.smi, type = "smile")
As the nature of fingerprint-based structure representation, the calculated 0-1 matrix x1, x2, and x3 will be very sparse. Since there are several columns have nearly exactly the same value for all the molecules, we should remove them with nearZeroVar() in caret before modeling, and split our training set and test set:
library("caret")
# remove near-zero variance variables
x1 = x1[, -nearZeroVar(x1)]
x2 = x2[, -nearZeroVar(x2)]
x3 = x3[, -nearZeroVar(x3)]
# split training and test set
set.seed(1003)
tr.idx = sample(1:nrow(x1), round(nrow(x1) * 0.75))
te.idx = setdiff(1:nrow(x1), tr.idx)
x1.tr = x1[tr.idx, ]
x1.te = x1[te.idx, ]
x2.tr = x2[tr.idx, ]
x2.te = x2[te.idx, ]
x3.tr = x3[tr.idx, ]
x3.te = x3[te.idx, ]
y.tr = y[tr.idx]
y.te = y[te.idx]
On the training sets, we will train three classification models separately using SVM (RBF kernel), using the kernlab package and the caret package. The cross-validation setting is 5-fold repeated CV (repeat for 10 times).
# svm classification on training sets
library("kernlab")
# cross-validation settings
ctrl = trainControl(method = "cv", number = 5, repeats = 10,
classProbs = TRUE,
summaryFunction = twoClassSummary)
# SVM with RBF kernel
svm.fit1 = train(
x1.tr, y.tr, method = "svmRadial", trControl = ctrl,
metric = "ROC", preProc = c("center", "scale"))
svm.fit2 = train(
x2.tr, y.tr, method = "svmRadial", trControl = ctrl,
metric = "ROC", preProc = c("center", "scale"))
svm.fit3 = train(
x3.tr, y.tr, method = "svmRadial", trControl = ctrl,
metric = "ROC", preProc = c("center", "scale"))
# print cross-validation result
print(svm.fit1)
print(svm.fit2)
print(svm.fit3)
The training result when using E-state fingerprints:
## Support Vector Machines with Radial Basis Function Kernel
##
## 597 samples
## 23 predictors
## 2 classes: "class0", "class1"
##
## Pre-processing: centered, scaled
## Resampling: Cross-Validated (5 fold)
##
## Summary of sample sizes: 478, 479, 477, 477, 477
##
## Resampling results across tuning parameters:
##
## C ROC Sens Spec ROC SD Sens SD Spec SD
## 0.25 0.797 0.7 0.765 0.0211 0.0442 0.00666
## 0.5 0.808 0.696 0.79 0.0173 0.059 0.0236
## 1 0.812 0.703 0.781 0.0191 0.0664 0.0228
##
## Tuning parameter "sigma" was held constant at a value of 0.02921559
## ROC was used to select the optimal model using the largest value.
## The final values used for the model were sigma = 0.0292 and C = 1.
We can see that, after removing the near zero variance predictors, there were only 23 predictors left for the original 79 E-state fingerprints.
The training result when using MACCS keys:
## Support Vector Machines with Radial Basis Function Kernel
##
## 597 samples
## 126 predictors
## 2 classes: "class0", "class1"
##
## Pre-processing: centered, scaled
## Resampling: Cross-Validated (5 fold)
##
## Summary of sample sizes: 477, 477, 477, 478, 479
##
## Resampling results across tuning parameters:
##
## C ROC Sens Spec ROC SD Sens SD Spec SD
## 0.25 0.834 0.715 0.775 0.0284 0.0994 0.0589
## 0.5 0.848 0.726 0.79 0.0299 0.065 0.0493
## 1 0.863 0.769 0.793 0.0307 0.0229 0.0561
##
## Tuning parameter "sigma" was held constant at a value of 0.004404305
## ROC was used to select the optimal model using the largest value.
## The final values used for the model were sigma = 0.0044 and C = 1.
There are 126 predictors left for the original 166 MACCS keys after removing the near-zero variance predictors. The model performance by AUC values is slightly better than using the E-state fingerprints.
The training result when using FP4 fingerprints:
## Support Vector Machines with Radial Basis Function Kernel
##
## 597 samples
## 58 predictors
## 2 classes: "class0", "class1"
##
## Pre-processing: centered, scaled
## Resampling: Cross-Validated (5 fold)
##
## Summary of sample sizes: 478, 478, 477, 478, 477
##
## Resampling results across tuning parameters:
##
## C ROC Sens Spec ROC SD Sens SD Spec SD
## 0.25 0.845 0.769 0.746 0.0498 0.0458 0.0877
## 0.5 0.856 0.744 0.777 0.0449 0.0148 0.0728
## 1 0.863 0.751 0.777 0.0428 0.036 0.0651
##
## Tuning parameter "sigma" was held constant at a value of 0.01077024
## ROC was used to select the optimal model using the largest value.
## The final values used for the model were sigma = 0.0108 and C = 1.
There are 58 predictors left for the original 512 FP4 fingerprints after the screening. The model performance by AUC values is almost the same comparing to using MACCS keys, and better than using E-state fingerprints.
We predict on the test sets with the established models, and plot the ROC curves in one figure, as is shown in Figure 4.
# predict on test set
svm.pred1 = predict(svm.fit1, newdata = x1.te, type = "prob")[, 1]
svm.pred2 = predict(svm.fit2, newdata = x2.te, type = "prob")[, 1]
svm.pred3 = predict(svm.fit3, newdata = x3.te, type = "prob")[, 1]
# generate colors
library("RColorBrewer")
pal = brewer.pal(3, "Set1")
# ROC curves of different fingerprints
library("pROC")
plot(smooth(roc(y.te, svm.pred1)), col = pal[1], grid = TRUE)
plot(smooth(roc(y.te, svm.pred2)), col = pal[2], grid = TRUE, add = TRUE)
plot(smooth(roc(y.te, svm.pred3)), col = pal[3], grid = TRUE, add = TRUE)
Figure 4: Smoothed ROC curves for different fingerprint types
Clustering of Molecules Based on Structural Similarities
Apart from supervised methods (classification and regression), unsupervised approaches, like clustering, is also widely applied in the quantitative research of drugs.
In reality, there are usually too many chemical compounds available for identifying drug-like molecules. Thus it would be attractive using clustering methods to aid the selection of a representative subset of all available compounds. For a clustering approach that groups compounds together by their structural similarity, applying the principle similar compounds have similar properties means that we only need to test the representative compounds from each individual cluster, rather than do the time-consuming complete set of experiments, and this should be sufficient to understand the structure-activity relationships of the whole compound set.
The Rcpi package provides easy-to-use functions for computing the similarity between small molecules derived by molecular fingerprints and maximum common substructure search.
As a example, the SDF file tyrphostin.sdf below is a database composed by searching tyrphostin in PubChem and filtered by _Lipinski“s rule of five_. We load this SDF file using readMolFromSDF():
library("Rcpi")
mols = readMolFromSDF(system.file(
"compseq/tyrphostin.sdf", package = "Rcpi"))
Then we compute the E-state fingerprints for all the molecules using extractDrugEstate(), and calculate their pairwise similarity matrix with calcDrugFPSim():
simmat = diag(length(mols))
for (i in 1:length(mols)) {
for (j in i:length(mols)) {
fp1 = extractDrugEstate(mols[[i]])
fp2 = extractDrugEstate(mols[[j]])
tmp = calcDrugFPSim(fp1, fp2, fptype = "compact", metric = "tanimoto")
simmat[i, j] = tmp
simmat[j, i] = tmp
}
}
With the computed similarity matrix simmat, we will try to cluster the molecules using hierarchical clustering, then visualize the clustering result:
mol.hc = hclust(as.dist(1 - simmat), method = "ward.D")
library("ape") # tree visualization of clusters
clus5 = cutree(mol.hc, 5) # cut dendrogram into 5 clusters
# generate colors
library("RColorBrewer")
pal5 = brewer.pal(5, "Set1")
plot(as.phylo(mol.hc), type = "fan",
tip.color = pal5[clus5],
label.offset = 0.1, cex = 0.7)
The hierarchical clustering result for these molecules is shown in Figure 5.
Figure 5: Tree visualization of the molecular clustering result
Structure-Based Chemical Similarity Searching
Structure-based chemical similarity searching ranks molecules in a database by their similarity degree to one query molecule structure. The numerical similarity value is usually computed based on the molecular fingerprints with selected metrics or by maximum common structure search. It is one of the core techniques for ligand-based virtual screening in drug discovery.
The SDF file DB00530.sdf below is retrieved from DrugBank, the drug ID DB00530 is Erlotinib, which is a reversible tyrosine kinase inhibitor. Given this molecule as the query molecule, we will do a similarity searching in the database tyrphostin.sdf presented in the last subsection.
library("Rcpi")
mol = system.file("compseq/DB00530.sdf", package = "Rcpi")
moldb = system.file("compseq/tyrphostin.sdf", package = "Rcpi")
We can do parallelized drug molecular similarity search with the searchDrug() function in Rcpi. Here we choose the search criterion to be MACCS keys with cosine similarity, FP2 fingerprints with tanimoto similarity, and maximum common substructure search with tanimoto similarity.
rank1 = searchDrug(
mol, moldb, cores = 4, method = "fp",
fptype = "maccs", fpsim = "tanimoto")
rank2 = searchDrug(
mol, moldb, cores = 4, method = "fp",
fptype = "fp2", fpsim = "cosine")
rank3 = searchDrug(
mol, moldb, cores = 4, method = "mcs",
mcssim = "tanimoto")
The returned search result is stored as a numerical vector, the name of each element is the molecule number in the database, and the value is the similarity value between the query molecule and this molecule. We print the top search results here:
head(rank1)
## 92 100 83 101 1 36
## 0.6491228 0.6491228 0.5882353 0.5660377 0.5000000 0.4861111
head(rank2)
## 100 92 83 101 94 16
## 0.8310005 0.8208663 0.5405856 0.5033150 0.4390790 0.4274081
head(rank3)
## 92 100 23 39 94 64
## 0.7000000 0.7000000 0.4000000 0.4000000 0.4000000 0.3783784
The Rcpi package also integrated the functionality of converting molecular file formats. For example, we can convert the SDF files to SMILES files using convMolFormat(). Since the #92 molecule ranks the highest in the three searches performed, we will calculate the similarity derived by maximum common substructure search between the query molecule and the #92 molecule using calcDrugMCSSim():
# convert SDF format to SMILES format
convMolFormat(
infile = mol, outfile = "DB00530.smi", from = "sdf", to = "smiles")
convMolFormat(
infile = moldb, outfile = "tyrphostin.smi", from = "sdf", to = "smiles")
smi1 = readLines("DB00530.smi")
smi2 = readLines("tyrphostin.smi")[92] # select the #92 molecule
calcDrugMCSSim(smi1, smi2, type = "smile", plot = TRUE)
The MCS search result is stored in a list, which contains the original MCS result provided by the fmcsR package Wang et al. (2013), the Tanimoto coefficient and the overlap coefficient.
## [[1]]
## An instance of "MCS"
## Number of MCSs: 1
## 530: 29 atoms
## 4705: 22 atoms
## MCS: 18 atoms
## Tanimoto Coefficient: 0.54545
## Overlap Coefficient: 0.81818
##
## [[2]]
## Tanimoto_Coefficient
## 0.5454545
##
## [[3]]
## Overlap_Coefficient
## 0.8181818
By using calcDrugMCSSim(..., plot = TRUE), the maximum common substructure of the two molecules is presented in Figure 6.
Figure 6: Maximum common structure of the query molecule and the #92 molecule in the chemical database
