Recently I’m learning conformal prediction and today I used the ready to use model for target prediction which is trained by chembl_24 (little bit old) with lightGBM.
To use the image in your local environment it is very easy just type following code.
docker run -p 8080:8080 chembl/mcp
After typing the code, docker image will run and server will be launched. Now target prediction server is ready to use. ;)
Original github repo describes how to use the server from terminal with curl command.
BTW, I would like to use it from python. It is easy to do it by using requests package.
My example is below.
example code
The model predicts that the compound has activity against Xanthine dehydrogenase. (Maybe the compound was included training data…)
It is interesting point that conformal prediction returns not only active/non active label but also both/empty label. It means the model predicts posi and nega with uncertainly.
To curate data and build model requires several steps and time. So ready to use system is user friendly. If one is allowed to wish so much I would like to build model with new version of chembl so I will be happy if source trainig code will be available. ;)
Fortunately ChEMBL provides common DB format dump file and RDKit has postgres DB cartridge. It means that you can search compound with rdkit functionality in postgres.
BTW, to handle the database, sqlalchemy which is ORMapper is very useful. So is it nice if we can communicate ChEMBL DB with sqlalchemy and rdkit function? I think yes! Fortunately razi provides us the methods to communicate postgres DB with chemoinformatics methods.
However, I couldn’t find how to communicate DB which has rdk schema already.
So I tried to do it today. Luckly, @ymasaKit_ provides very useful tutorial in his repo! The document is written in Japanese but it’s useful because code is written in python. ;)
OK let’s start today’s code! At first, install postgresql, rdkit-cartridge and ChEMBL27. I used conda for convenience.
Now I installed chembl27 to my postgresql DB. Then make rdk schema. Details are described in rdkit original document. URL is below. https://www.rdkit.org/docs/Cartridge.html
chembl_27=# create extension if not exists rdkit;
chembl_27=# create schema rdk;
chembl_27=# select * into rdk.mols from (select molregno,mol_from_ctab(molfile::cstring) m from compound_structures) tmp where m is not null;
chembl_27=# create index molidx on rdk.mols using gist(m);
CREATE INDEX
chembl_27=# alter table rdk.mols add primary key (molregno);
ALTER TABLE
chembl_27=# select molregno,torsionbv_fp(m) as torsionbv,morganbv_fp(m) as mfp2,featmorganbv_fp(m) as ffp2 into rdk.fps from rdk.mols;
chembl_27=# create index fps_ttbv_idx on rdk.fps using gist(torsionbv);
CREATE INDEX
chembl_27=# create index fps_mfp2_idx on rdk.fps using gist(mfp2);
CREATE INDEX
chembl_27=# create index fps_ffp2_idx on rdk.fps using gist(ffp2);
CREATE INDEX
chembl_27=# alter table rdk.fps add primary key (molregno);
ALTER TABLE
Now I could make rdk.mols and rdk.fps table in my chembl_27 DB.
Let’s move to python!. To communicate postgresql with sqlalchemy, psycopg2 should be installed at first. Map mols and fps table at first. rdk schema has custom column such as rdkit mol object and fingerprints so autoload option should not use. And extend_existing=True is used for mapping mol column.
from rdkit import Chem
from rdkit.Chem import Draw
from rdkit.Chem.Draw import IPythonConsole
from sqlalchemy import create_engine, MetaData
from sqlalchemy.orm import sessionmaker
from sqlalchemy import Table
from sqlalchemy import Column, BIGINT, Index
from sqlalchemy.orm import mapper
from sqlalchemy.ext.declarative import declarative_base
from razi import rdkit_postgresql
from razi.rdkit_postgresql.types import Mol
from razi.rdkit_postgresql.types import Bfp
db = create_engine('postgres://iwatobipen@localhost:5432/chembl_27',
#echo=True
)
metadata = MetaData(schema='rdk')
Base = declarative_base()
class Mols(Base):
__table__ = Table('mols',
metadata,
Column('molregno', BIGINT, primary_key=True),
Column('m', Mol),
extend_existing=True,
)
class Fps(Base):
__table__ = Table('fps',
metadata,
Column('molregno', BIGINT, primary_key=True),
Column('torsionbv', Bfp),
Column('mfp2', Bfp),
Column('ffp2', Bfp),
extend_existing=True)
Now I could map rdk.mol, rdk.fps table to my defined class. Next I tested some chmoinformatics related query.
First example is retrieve data from rdk.mols table.
ession = sessionmaker(bind=db)
session = Session()
cpds = session.query(Mols)
from IPython.display import display
for row in cpds[:4]:
display(row.m)
print(row.molregno)
Next, retrieve compounds with substructure query.
cpds = session.query(Mols)
cpds = cpds.filter(Mols.m.hassubstruct('n1cccc2ccccc12'))
for row in cpds[:4]:
display(row.m)
print(row.molregno)
Seems work fine. ;) Next example is query combination substructure and molecular weight. Razi supports queries which are defined in rdkit-postgres cartridge. Substructure query is same as above.
from razi.rdkit_postgresql.functions import mol_amw
cpds = session.query(Mols)
cpds = cpds.filter(Mols.m.hassubstruct('n1cccc2ccccc12')).filter(mol_amw(Mols.m) <= 200)
for row in cpds[:4]:
display(row.m)
print(row.molregno)
And last example is similarity search with rdk.fps table. At first I changed threshold with SQL statement. And then search fps table and display matched compounds.
res = fps.filter(Fps.ffp2.tanimoto_sml(fp))
regno = []
for row in res:
print(row.molregno)
regno.append(row.molregno)
cpds = session.query(Mols)
cpds = cpds.filter(Mols.molregno.in_(regno))
for cpd in cpds[:10]:
display(cpd.m)
All examples seems work fine. Razi and posgresql combination is useful because user can handle database without SQL statement and more pythonic.
Another work, @fmkz__ developed pychembldb. It is useful package to communicate ChEMBLDB with python. So I think it is interesting that pychembldb enhancement with razi. ;)
I posted about conformal prediction with python and rdkit some days ago. After that I could get very informative advice from @kjelljorner. Thanks a lot! His advice was below.
Kjell Jorner @kjelljorner3dReplying to @iwatobipen I can recommend the cross conformal prediction or bootstrapped conformal prediction (also in nonconformist) to avoid having to put aside data for calibration. There is also a related method called the jackknife+ which can be applied to any method:
Fortunately nonconformist has some implementations of ACP. So I tried to use acp module with same dataset.
For convenience, shared my gist link below instead of pasting code on the blog post.
In the post, I tested AggregatedCp and BootstrapConfomralClassifier. These methods internally make train and calibration set with defined sampling method. In the default setting, aggregatedCp uses bootstrapsampler. So both method is almost same I think. And results are almost same.
By using acp method, used don’t need to separate train data to train and calibration data before training. It seems natural and useful way to use conformal prediction against chemoinformatics problems.
BTW, in the these area, there are many skillful people in the world and they have an open mind. I feel writing the post is useful for me because it keeps my idea, tips code(this is main reason for keeping the blog) but also gives opportunity for having a discussion about chemoinformatics and drug discovery. ;)
He introduced basic concept of conformal prediction and demonstration with excellent work flow of KNIME. I recommend to check the site. ;)
Conformal prediction is not new method. It can estimate confidence of predicted values. Traditional predictive model can predict probability, or class of given task (classification task) but can’t estimate how confidence the predicted values are.
I like KIME but I love python too so I searched package to conduct conformal prediction with python. The package name is ‘nonconformist‘. The packge depends scikit-learn but it don’t support recent version of scikit-learn. So some modification is required (PR is original repo, hope it will marge soon). So, I modified original code and uploaded my repository. https://github.com/iwatobipen/nonconformist
Now ready to use nonconformist. Following example data comes from rdkit provided solubility data set and use ICP (inductive conformal prediction). At first import packages.
import os
import pandas as pd
import numpy as np
from rdkit import Chem
from rdkit.Chem import RDConfig
from rdkit.Chem import DataStructs
from rdkit.Chem import AllChem
from rdkit.Chem.Draw import IPythonConsole
from rdkit.Chem import Draw
from sklearn.ensemble import RandomForestClassifier
from nonconformist.nc import ClassifierNc
from nonconformist.nc import ClassifierAdapter
from nonconformist.icp import IcpClassifier
from nonconformist.evaluation import ClassIcpCvHelper
from nonconformist.evaluation import cross_val_score
Then load train and test data set. For comformal prediction, calibration data is required so I divided train and calibration data. And calculate molecular fingerprint. In the knime tutorial I mentioned above, lots of train calibration data is generated for train model. Following code, I made one train – calibration table(data).
train = os.path.join(RDConfig.RDDocsDir, 'Book/data/solubility.train.sdf')
test = os.path.join(RDConfig.RDDocsDir, 'Book/data/solubility.test.sdf')
trainmol = [m for m in Chem.SDMolSupplier(train)]
testmol = [m for m in Chem.SDMolSupplier(test)]
label2cls = {'(A) low':0, '(B) medium':1, '(C) high':2}
def fp2arr(fp):
arr = np.zeros((1,))
DataStructs.ConvertToNumpyArray(fp, arr)
return arr
trainfps = [AllChem.GetMorganFingerprintAsBitVect(m, 2, 1024) for m in trainmol]
trainfps = np.array([fp2arr(fp) for fp in trainfps])
testfps = [AllChem.GetMorganFingerprintAsBitVect(m, 2, 1024) for m in testmol]
testfps = np.array([fp2arr(fp) for fp in testfps])
train_cls = [label2cls[m.GetProp('SOL_classification')] for m in trainmol]
train_cls = np.array(train_cls)
test_cls = [label2cls[m.GetProp('SOL_classification')] for m in testmol]
test_cls = np.array(test_cls)
print(trainfps.shape, train_cls.shape, testfps.shape, test_cls.shape)
#train data is divided to train and calibration data
ids = np.random.permutation(train_cls.size)
# Use first 700 data for train and second set is used for calibration
trainX, calibX = trainfps[ids[:700],:],trainfps[ids[700:],:]
trainY, calibY = train_cls[ids[:700]],train_cls[ids[700:]]
testX = testfps
testY = test_cls
Now ready to train with nonconformist. I used RandomForestClassifier for base classifier. And made IcpClassifier object with random forest.
Interestingly, classification with conformal prediction result has many combination of classes. My data set has 3 classes, so traditional classifier will returns 3 classes [1,0,0], [0,1,0], [0,0,1]. However conformal predictor compares significanse and each predicted p-value of class. And if the p-value > significance, it labeled true. So the predictor has possibility of return following combinations. [0,0,0], [0,0,1], [0,1,1], [1,1,1], [1,0,0], [0,1,0],[1,1,0],[1,0,1]..
Fortunately nonconformist has cross validation method. To do that, ClassIcpCvHelper is used.
icpmodel = ClassIcpCvHelper(icp)
res = cross_val_score(icpmodel,
trainfps,
train_cls,
iterations=10,
scoring_funcs=[class_avg_c],
significance_levels=[0.05, 0.1, 0.2],
)
#it will take few sec-min
Evaluation method is difficult to under stand for me. For example, class_n_correct method returns ratio of predicted true / ground truth. But, conformal prediction return multiple true labels, so max values will be over number of true labels.
Conformal prediction is good tool for QSAR, I would like to learn knime workflow of CP in this week end.
The author shared their code! Thanks. I have interest the code. So I installed the code and used the package. Following code is almost same as original code but reactions are different.
The model can annotate complex multi component ‘Ugi’ reaction. But failed to map simple Diels Alder reaction and Claisen Rearrangement.
I would like to test more reactions for evaluate the model. It seems very useful for me.
And now I’m reading the article. Thanks for sharing nice code and article!
Recently there are many publications about de-novo molecular generator which mainly use Deep Learning. One problem of the approach is that generated molecules are not systematic so it’s difficult to synthesis them with parallel chemistry. So sometime chemists dislike the proposal from generated form the method I think.
Rule or Rxn or MMP based molecule generation is another approach to do that. It’s based on more chemist friend rules. They are not new but useful method and also related approaches are still reported in these days.
The author proposed new workflow for molecular framework mutation it seems like MMP approach, it degrade molecule to fragment with local context (radi1-5) for making interchangable fragment database, like MMP key-value structures are stored. And the data is used for ‘MUTATE’, ‘GROW’ and ‘LINK’ for new structure generation.
They proposed similar approach for molecular generation with fragment database.
Compared these approach, I think main difference is that CreM can set context radius. The setting affects feature of generated molecules.
In the fig4, and fig5 of the J. Cheminform article, the author shows properties of generated molecules with different radius. For example novelty, diversity score is decreased when large radius(5) is used. It means that more context similar compounds are generated with the setting.
As you know, CreM author disclosed the implementation so let’s use it. It is easy to install crem has very few dependency, just rdkit and gaucamol(optional). At first I installed CReM with pip and get ready to use DB.
By using the dataset, it took few minutes for structure generation. After generating the molecule rdkit can render mols with Drawing function.
It seems that radius=1 generates more diverse compounds set. It is easy to use for molecular generation.
Ok now we can use deep learning based and rule based structure generator. Each methods has pros and cons. As author said that CReM can generate chemistry reasonable structure but can’t generate new rings which isn’t fragment db.
Which is good proposal for medchem new structure constructed from know fragments or new structure with novel fragments?
It’s depends on situation but novel fragments requires new chemistry or many wet experiments. AI driven drug discovery can’t replace all wet experiments to dry experiments. Which molecules do you make at first and next, experimental design is key for the many projects.
Here is an example for draw molecules with FG highlighting. Improved version of the code can be found in last part of the post.
RDKit has fdef file for common functional groups. So we can use it for highlighting atoms as same as pharmacohore. Today’s example is just for my memo. Because I often forget many useful options.
MolDraw2DSVG object can draw molecules with DrawMolecules method. So it is easy to draw molecules as high quality SVG image.
Greg suggested me following improved version of the code! Thanks a lot!!
