Rcpi: R/Bioconductor Package as an Integrated Informatics Platform for Drug Discovery

Nan Xiao <https://nanx.me>

2018-07-16

1 Introduction

The Rcpi package Cao et al. (2015) offers an R/Bioconductor package emphasizing the comprehensive integration of bioinformatics and chemoinformatics into a molecular informatics platform for drug discovery.

Rcpi implemented and integrated the state-of-the-art protein sequence descriptors and molecular descriptors/fingerprints with R. For protein sequences, the Rcpi package can

For small molecules, the Rcpi package can

By combining various types of descriptors for drugs and proteins in different methods, interaction descriptors representing protein-protein or compound-protein interactions can be conveniently generated with Rcpi, including:

Several useful auxiliary utilities are also shipped with Rcpi:

The computed protein sequence descriptors, molecular descriptors/fingerprints, interaction descriptors and pairwise similarities are widely used in various research fields relevant to drug disvery, such as, bioinformatics, chemoinformatics, proteochemometrics and chemogenomics.

1.1 Installation

To install the Rcpi package, use:

To make the Rcpi package fully functional (especially the Open Babel related functionalities), we recommend the users also install the Enhances packages by:

Several dependencies of the Rcpi package may require some system-level libraries, check the corresponding manuals of these packages for detailed installation guides.

1.2 How to Cite Rcpi

If you feel Rcpi is useful in your research, please cite our paper:

Dong-Sheng Cao, Nan Xiao, Qing-Song Xu, and Alex F. Chen. (2015). Rcpi: R/Bioconductor package to generate various descriptors of proteins, compounds and their interactions. Bioinformatics 31 (2), 279-281.

BibTeX entry:

2 Applications in Bioinformatics

For bioinformatics research, Rcpi calculates commonly used descriptors and proteochemometric (PCM) modeling descriptors for protein sequences. Rcpi also computes pairwise similarities derived by GO semantic similarity and sequence alignment.

2.1 Predicting Protein Subcellular Localization

Protein subcellular localization prediction involves the computational prediction of where a protein resides in a cell. It is an important component of bioinformatics-based prediction of protein function and genome annotation, and can also aid us to identify novel drug targets.

Here we use the subcellular localization dataset of human proteins presented in the study of Chou and Shen (2008) for a demonstration. The complete dataset includes 3,134 protein sequences (2,750 different proteins), classified into 14 human subcellular locations. We selected two classes of proteins as our benchmark dataset. Class 1 contains 325 extracell proteins, and class 2 includes 307 mitochondrion proteins.

First, we load the Rcpi package, then read the protein sequences stored in two separated FASTA files with readFASTA():

The loaded sequences are stored as two lists in R, and each component in the list is a character string representing one protein sequence. In this case, there are 325 extracell protein sequences and 306 mitonchon protein sequences:

To assure that the protein sequences only have the twenty standard amino acid types which is required for the descriptor computation, we use the checkProt() function to do the amino acid type sanity checking and remove the non-standard protein sequences:

Two protein sequences were removed from each class. For the remaining sequences, we calculate the amphiphilic pseudo amino acid composition (APAAC) descriptor Chou (2005) and create class labels for the random forest classification modeling.

In Rcpi, the functions of commonly used descriptors for protein sequences and proteochemometric (PCM) modeling descriptors are named after extractProt...() and extractPCM...().

Next, we will split the data into a 75% training set and a 25% test set.

We will train a random forest classification model on the training set with 5-fold cross-validation, using the package.

The training result is:

With the model trained on the training set, we predict on the test set and plot the ROC curve with the package, as is shown in Figure 1.

ROC curve for the test set of protein subcellular localization data

Figure 1: ROC curve for the test set of protein subcellular localization data

The area under the ROC curve (AUC) is:

3 Applications in Chemoinformatics

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.

3.1 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:

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:

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).

The cross-validation result is:

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.

Number of principal components vs. RMSE for the PLS regression model

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):

Experimental RIs vs. Predicted RIs

Figure 3: Experimental RIs vs. Predicted RIs

3.2 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:

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:

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:

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).

The training result when using E-state fingerprints:

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:

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:

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.

Smoothed ROC curves for different fingerprint types

Figure 4: Smoothed ROC curves for different fingerprint types

3.3 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():

Then we compute the E-state fingerprints for all the molecules using extractDrugEstate(), and calculate their pairwise similarity matrix with calcDrugFPSim():

With the computed similarity matrix simmat, we will try to cluster the molecules using hierarchical clustering, then visualize the clustering result:

The hierarchical clustering result for these molecules is shown in Figure 5.

Tree visualization of the molecular clustering result

Figure 5: Tree visualization of the molecular clustering result

3.4 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.

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.

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:

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():

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.

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

4 Applications in Chemogenomics

For chemogenomics modeling, Rcpi calculates compound-protein interaction (CPI) descriptors and protein-protein interaction (PPI) descriptors.

4.1 Predicting Drug-Target Interaction by Integrating Chemical and Genomic Spaces

The prediction of novel interactions between drugs and target proteins is a key area in genomic drug discovery. In this example, we use the G Protein-Coupled Receptor (GPCR) dataset provided by Yamanishi et al. (2008) as our benchmark dataset.

A drug-target interaction network can be naturally modeled as a bipartite graph, where the nodes are target proteins or drug molecules and edges (only drugs and proteins could be connected by edges) represent drug-target interactions. Initially, the graph only contains edges describing the real drug-target interactions determined by experiments or other ways. In this example, all real drug-target interaction pairs (i.e., 635 drug-target interactions) are used as the positive samples. For negative samples we select random, non-interacting pairs from these drugs and proteins. They are constructed as follows:

  1. Separate the pairs in the above positive samples into single drugs and proteins;
  2. Re-couple these singles into pairs in a way that none of them occurs in the corresponding positive dataset.

Ten generated negative sets were used in Dong-Sheng Cao et al. (2012), here we only used one of them for a demonstration. The drug ID and target ID is stored in GPCR.csv. The first column is KEGG protein ID, and the second column is KEGG drug ID. The first 635 rows form the positive set, and the last 635 rows form the negative set.

Take a look at the data:

We will visualize the network first. Figure 7 shows the connection pattern for the GPCR drug-target interaction network in the form of an arc diagram.

Arc diagram visualization of the GPCR drug-target interaction network

Figure 7: Arc diagram visualization of the GPCR drug-target interaction network

An arc diagram visualize the nodes in the network in a one-dimensional layout, while using circular arcs to represent edges. With a good ordering of nodes, it is easy to identify cliques and bridges.

Next, we will download the target protein sequences (in FASTA format) and drug molecule (in SMILES format) from the KEGG database, in parallel:

If the connection is slow or accidentally interrupts, just try a few more times until success.

The functions in Rcpi named after getMolFrom...() and getSmiFrom...() supports the parallelized retrieval of (drug) molecules from PubChem, ChEMBL, CAS, KEGG, and DrugBank. The functions named after getSeqFrom...(), getFASTAFrom...(), and getPDBFrom...() supports the parallelized retrieval of proteins from UniProt, KEGG and RCSB PDB. The functions getDrug() and getProt() are two integrated wrapper functions for downloading the molecules and protein sequences from these online databases.

After the sequences were downloaded, we can calculate the protein sequence descriptors and molecular descriptors for the targets and drugs:

Since the descriptors is only for the deduplicated drug and target list, we need to generate the full descriptor matrix for the training data:

Generate drug-target interaction descriptors using getCPI().

The pairwise interaction is another useful type of representation in drug-target prediction, protein-protein interaction prediction and related research. Rcpi also provides getPPI() to generate protein-protein interaction descriptors. getPPI() provides three types of interactions while getCPI() provides two types. The argument type is used to control this.

4.1.1 Compound-protein interaction (CPI) descriptors

For compound descriptor vector \(d_1^{1 \times p_1}\) and the protein descriptor vector \(d_2^{1 \times p_2}\), there are two methods for construction of descriptor vector \(d\) for compound-protein interaction:

  1. type = "combine" - combine the two feature matrix, \(d\) has \(p_1 + p_2\) columns;
  2. type = "tensorprod" - column-by-column (pseudo)-tensor product type interactions, \(d\) has \(p_1 \times p_2\) columns.

4.1.2 Protein-protein interaction (PPI) descriptors

For interaction protein A and protein B, let \(d_1^{1 \times p}\) and \(d_2^{1 \times p}\) be the descriptor vectors. There are three methods to construct the protein-protein interaction descriptor \(d\):

  1. type = "combine" - combine the two descriptor matrix, \(d\) has \(p + p\) columns;
  2. type = "tensorprod" - column-by-column (pseudo)-tensor product type interactions, \(d\) has \(p \times p\) columns;
  3. type = "entrywise" - entrywise product and entrywise sum of the two matrices, then combine them, \(d\) has \(p + p\) columns.

Train a random forest classification model with 5-fold repeated CV:

Predict on the training set (for demonstration purpose only) and plot ROC curve.

The ROC curve is shown in Figure 8.

ROC curve for predicting on the training set of the GPCR drug-target interaction dataset using random forest

Figure 8: ROC curve for predicting on the training set of the GPCR drug-target interaction dataset using random forest

References

Cao,D.-S. et al. (2012) Large-scale prediction of drug-target interactions using protein sequences and drug topological structures. Analytica chimica acta, 752, 1–10.

Cao,D.-S. et al. (2015) Rcpi: R/Bioconductor package to generate various descriptors of proteins, compounds and their interactions. Bioinformatics, 31, 279–281.

Cao,D.-S. et al. (2012) In silico toxicity prediction by support vector machine and smiles representation-based string kernel. SAR and QSAR in Environmental Research, 23, 141–153.

Chou,K.-C. (2005) Using amphiphilic pseudo amino acid composition to predict enzyme subfamily classes. Bioinformatics, 21, 10–19.

Chou,K.-C. and Shen,H.-B. (2008) Cell-ploc: A package of web servers for predicting subcellular localization of proteins in various organisms. Nature Protocols, 3, 153–162.

Horan,K. and Girke,T. (2013) ChemmineOB: R interface to a subset of openbabel functionalities.

Steinbeck,C. et al. (2003) The chemistry development kit (cdk): An open-source java library for chemo-and bioinformatics. Journal of chemical information and computer sciences, 43, 493–500.

Wang,Y. et al. (2013) FmcsR: Mismatch tolerant maximum common substructure searching in r. Bioinformatics, 29, 2792–2794.

Yamanishi,Y. et al. (2008) Prediction of drug–target interaction networks from the integration of chemical and genomic spaces. Bioinformatics, 24, i232–i240.

Yan,J. et al. (2012) Comparison of quantitative structure–retention relationship models on four stationary phases with different polarity for a diverse set of flavor compounds. Journal of Chromatography A, 1223, 118–125.