AI + Chemistry: Building Drug Discovery Pipelines with Free Tools
An ABL kinase project can be assembled from free tools. The snippets below show the calls. They are not a trained screen: several of them run on two molecules, imatinib and a toy inactive, so they demonstrate the interface and nothing about predictive performance.
The target in the example is ABL kinase, the kinase inhibited by imatinib. Bcr-Abl drives chronic myeloid leukemia. The practical sequence is: pull measured activities, featurize, fit a model, write a 3D structure, look at it.
Data
ChEMBL is the open source of IC50 and Ki values for compounds tested against ABL. The web services or a bulk download feed a table. Pandas cleans it. NumPy does the arithmetic. This sketch starts from a CSV that is already local.
import pandas as pd
df = pd.read_csv("ABL_bioactivity.csv")
RDKit
Morgan fingerprints, a Tanimoto, and a substructure check. The SMILES are imatinib and nilotinib. The SMARTS is the shared aniline–pyrimidine motif used in the example.
from rdkit import Chem, DataStructs
from rdkit.Chem import AllChem
imatinib_smiles = "Cc1ccc(cc1Nc2nccc(n2)c3cccnc3)NC(=O)c4ccc(cc4)CN5CCN(CC5)C"
nilotinib_smiles = "Cc1ccc(cc1Nc2nccc(n2)c3cccnc3)C(=O)Nc4cc(cc(c4)n5cc(nc5)C)C(F)(F)F"
mol1 = Chem.MolFromSmiles(imatinib_smiles)
mol2 = Chem.MolFromSmiles(nilotinib_smiles)
fp1 = AllChem.GetMorganFingerprintAsBitVect(mol1, radius=2, nBits=2048)
fp2 = AllChem.GetMorganFingerprintAsBitVect(mol2, radius=2, nBits=2048)
sim = DataStructs.TanimotoSimilarity(fp1, fp2)
print(f"Tanimoto similarity: {sim:.2f}")
query = Chem.MolFromSmarts("c1ccc(cc1)Nc2nccc(n2)")
match = mol1.HasSubstructMatch(query)
print("Substructure match:", match)
DeepChem
A graph-convolution classifier on two labeled SMILES, fit for 20 epochs. That is enough to see the objects: a featurizer, a NumpyDataset, and GraphConvModel in classification mode with ROC AUC as the metric. It is not an ABL model.
import deepchem as dc
import numpy as np
smiles_list = [
"Cc1ccc(cc1Nc2nccc(n2)c3cccnc3)NC(=O)c4ccc(cc4)CN5CCN(CC5)C", # active
"CCOC(=O)c1ccc(cc1)N", # inactive
]
labels = np.array([1, 0])
featurizer = dc.feat.MolGraphConvFeaturizer()
X = featurizer.featurize(smiles_list)
y = labels
dataset = dc.data.NumpyDataset(X, y)
model = dc.models.GraphConvModel(
n_tasks=1,
mode="classification",
metrics=[dc.metrics.Metric(dc.metrics.roc_auc_score)],
)
model.fit(dataset, nb_epoch=20)
pred_probs = model.predict(dataset)
print(pred_probs)
Open Babel
SMILES to a 3D SDF, with hydrogens.
from openbabel import pybel
smiles = "Cc1ccc(cc1Nc2nccc(n2)c3cccnc3)NC(=O)c4ccc(cc4)CN5CCN(CC5)C"
mol = pybel.readstring("smi", smiles)
mol.addh()
mol.make3D()
mol.write("sdf", "imatinib_3D.sdf", overwrite=True)
PyMOL
Load a kinase structure and the SDF. ABL_kinase.pdb is a filename, not a file this note provides.
import pymol2
with pymol2.PyMOL() as pymol:
cmd = pymol.cmd
cmd.load("ABL_kinase.pdb", "protein")
cmd.load("imatinib_3D.sdf", "ligand")
cmd.hide("everything")
cmd.show("cartoon", "protein")
cmd.show("sticks", "ligand")
cmd.zoom("ligand", 5)
cmd.png("abl_imatinib.png", width=800, height=600, ray=1)
Chemprop
Command line, then the same training call from Python. Both expect abl_activity.csv with a smiles column and an active column. Thirty epochs in the Python dict, classification mode, output directory abl_model.
chemprop_train --data_path abl_activity.csv --smiles_column smiles --target_columns active \
--dataset_type classification --save_dir abl_model
from chemprop.train import run_training
params = {
"data_path": "abl_activity.csv",
"smiles_column": "smiles",
"target_columns": ["active"],
"dataset_type": "classification",
"save_dir": "abl_model",
"epochs": 30,
}
run_training(params)
A random forest on the two fingerprints
This continues from the RDKit fingerprints and the two-row labels array above. A 20% split of two rows is not a validation set. The call is here so the scikit-learn shape is visible.
from sklearn.ensemble import RandomForestClassifier
from sklearn.model_selection import train_test_split
X = [list(map(int, fp.ToBitString())) for fp in [fp1, fp2]]
y = labels
X_train, X_val, y_train, y_val = train_test_split(X, y, test_size=0.2, random_state=0)
model = RandomForestClassifier(n_estimators=100).fit(X_train, y_train)
print("Validation accuracy:", model.score(X_val, y_val))
What each tool is for
ChEMBL and PubChem hold bioactivity. RDKit handles structures. DeepChem and Chemprop are the message-passing and graph-convolution calls. Open Babel converts formats. PyMOL inspects a pose. Pandas, NumPy, and scikit-learn are the table and the baseline. All of these are free. A competitive model still depends on the table you did not sketch here: more than two compounds, a split that respects the chemistry, and a metric you compute on held-out molecules.
