TP komet : practical work

In this TP, we will present the LCIdb database and the Komet algorithm.

  1. Downloading utils and the dataset (train/val/test) from Github

  2. Calculation of molecule features using a subsample of train molecules (MorganFP kernel approximated via Nystrom approximation)

  3. Calculation of protein features (SVD on LAkernel)

  4. Train/Testing with a chosen lambda

  5. Make some predictions

You will find in the Komet github page, in the Test folder, the correction of this TP, executed in Colab.

[1]:

%load_ext autoreload
%autoreload 2

# Importing the current libraries

import matplotlib.pyplot as plt
import pandas as pd
import numpy as np
import zipfile
import pickle
import time

from sklearn.metrics import  average_precision_score,  roc_curve, confusion_matrix, precision_score, recall_score, auc

[2]:

# Importing the pytorch libraries and define the device

import torch
import torch.optim as optim

device = torch.device("cuda:0" if torch.cuda.is_available() else "cpu")
device_cpu = torch.device("cpu")
device_cpu = device
print( device )

mytype = torch.float16 # to save memory (only on GPU)
mytype = torch.float32

cpu

0. Presentation of the rdkit library

[3]:

!pip install rdkit

Requirement already satisfied: rdkit in /home/docs/checkouts/readthedocs.org/user_builds/komet/envs/latest/lib/python3.9/site-packages (2025.3.3)
Requirement already satisfied: numpy in /home/docs/checkouts/readthedocs.org/user_builds/komet/envs/latest/lib/python3.9/site-packages (from rdkit) (2.0.2)
Requirement already satisfied: Pillow in /home/docs/checkouts/readthedocs.org/user_builds/komet/envs/latest/lib/python3.9/site-packages (from rdkit) (11.2.1)

Question 1 Consult the following link and draw the aspirin molecule in the code box below.

[ ]:

1. Download the data from a GitHub repo.

Download utils

[4]:

!wget -q https://github.com/Guichaoua/komet/raw/main/komet/komet.py

[5]:

import komet

cpu

Question 2 In the file icon on the left, download the komet.py file and observe its structure.

Download train/val/test

[6]:

!mkdir data/
!wget -q https://github.com/Guichaoua/komet/raw/main/data/LCIdb/Orphan/dict_ind2fasta_LCIdb.data
!mv dict_ind2fasta_LCIdb.data data/
!wget -q https://github.com/Guichaoua/komet/raw/main/data/LCIdb/Orphan/dict_ind2fasta_LCIdb_K_prot.data
!mv dict_ind2fasta_LCIdb_K_prot.data data/
!wget -q https://github.com/Guichaoua/komet/raw/main/data/LCIdb/Orphan/test.csv
!mv test.csv data/
!wget -q https://github.com/Guichaoua/komet/raw/main/data/LCIdb/Orphan/train.csv.zip
!mv train.csv.zip data/
!wget -q https://github.com/Guichaoua/komet/raw/main/data/LCIdb/Orphan/val.csv
!mv val.csv data/

Question 3 Check that the files have been uploaded successfully in data file.

[7]:

dataset_dir = "data/"

# load data
train = komet.load_df("train.csv.zip",dataset_dir)
val = komet.load_df("val.csv",dataset_dir)
test = komet.load_df("test.csv",dataset_dir)

# dataframe full has all smiles and fasta sequences
full = pd.concat([train, val, test])
print("full shape",full.shape)

[09:42:41] Explicit valence for atom # 5 P, 7, is greater than permitted
[09:42:41] Explicit valence for atom # 23 P, 7, is greater than permitted
[09:42:41] Explicit valence for atom # 19 P, 7, is greater than permitted
[09:42:42] Explicit valence for atom # 10 P, 7, is greater than permitted
[09:42:42] Explicit valence for atom # 5 P, 7, is greater than permitted
[09:42:42] Explicit valence for atom # 16 P, 7, is greater than permitted
[09:42:42] Explicit valence for atom # 20 P, 7, is greater than permitted
[09:42:43] Explicit valence for atom # 5 P, 7, is greater than permitted
[09:42:43] Explicit valence for atom # 9 P, 7, is greater than permitted
[09:42:43] Explicit valence for atom # 5 P, 7, is greater than permitted
[09:42:43] Explicit valence for atom # 10 P, 7, is greater than permitted
[09:42:44] Explicit valence for atom # 17 P, 7, is greater than permitted
[09:42:44] Explicit valence for atom # 10 P, 7, is greater than permitted
[09:42:45] Explicit valence for atom # 10 P, 7, is greater than permitted
[09:42:49] Explicit valence for atom # 13 P, 7, is greater than permitted
[09:42:49] Explicit valence for atom # 7 P, 7, is greater than permitted
[09:42:50] Explicit valence for atom # 13 P, 7, is greater than permitted
[09:42:50] Explicit valence for atom # 14 P, 7, is greater than permitted
[09:42:50] Explicit valence for atom # 3 P, 7, is greater than permitted
[09:42:51] Explicit valence for atom # 1 P, 7, is greater than permitted
[09:42:52] Explicit valence for atom # 17 P, 7, is greater than permitted
[09:42:52] Explicit valence for atom # 3 P, 7, is greater than permitted
[09:42:52] Explicit valence for atom # 1 P, 7, is greater than permitted
[09:42:54] Explicit valence for atom # 1 P, 7, is greater than permitted
[09:42:55] Explicit valence for atom # 11 P, 7, is greater than permitted
[09:42:55] Explicit valence for atom # 14 P, 7, is greater than permitted
[09:42:55] Explicit valence for atom # 1 P, 7, is greater than permitted
[09:42:55] Explicit valence for atom # 18 P, 7, is greater than permitted
[09:42:55] Explicit valence for atom # 46 P, 7, is greater than permitted
[09:42:55] Explicit valence for atom # 37 P, 7, is greater than permitted
[09:42:55] Explicit valence for atom # 1 P, 7, is greater than permitted
[09:42:55] Explicit valence for atom # 1 P, 7, is greater than permitted
[09:42:55] Explicit valence for atom # 13 P, 7, is greater than permitted
[09:42:55] Explicit valence for atom # 3 P, 7, is greater than permitted
[09:42:56] Explicit valence for atom # 16 P, 7, is greater than permitted

number of smiles to clean: 35
train.csv shape (236396, 3)

[09:42:56] Explicit valence for atom # 43 P, 7, is greater than permitted
[09:42:58] Explicit valence for atom # 8 P, 7, is greater than permitted
[09:42:58] Explicit valence for atom # 20 P, 7, is greater than permitted
[09:42:59] Explicit valence for atom # 11 P, 7, is greater than permitted

number of smiles to clean: 4
val.csv shape (21839, 3)

[09:43:01] Explicit valence for atom # 19 P, 7, is greater than permitted
[09:43:03] Explicit valence for atom # 18 P, 7, is greater than permitted
[09:43:03] Explicit valence for atom # 16 P, 7, is greater than permitted
[09:43:03] Explicit valence for atom # 16 P, 7, is greater than permitted
[09:43:03] Explicit valence for atom # 1 P, 7, is greater than permitted
[09:43:03] Explicit valence for atom # 16 P, 7, is greater than permitted
[09:43:03] Explicit valence for atom # 18 P, 7, is greater than permitted
[09:43:03] Explicit valence for atom # 9 P, 7, is greater than permitted
[09:43:04] Explicit valence for atom # 17 P, 7, is greater than permitted
[09:43:04] Explicit valence for atom # 8 P, 7, is greater than permitted
[09:43:05] Explicit valence for atom # 4 P, 7, is greater than permitted

number of smiles to clean: 11
test.csv shape (44985, 3)
full shape (303220, 3)

Question 4

  1. Go to the komet.py file to see what the load_df function does.

  2. What is the type of train? Show first lines. How many values ​​for the Label?

[ ]:

2. Calculation of molecule features using a subsample of train molecules (molecule kernel approximated via Nystrom approximation)

[8]:

#### MOLECULE####

list_smiles = full[['SMILES']].drop_duplicates().values.flatten()
nM = len(list_smiles)
print("number of different smiles (mol):",nM)

# add indsmiles in train, val, test
#dict_ind2smiles = {i:list_smiles[i] for i in range(nM)}
dict_smiles2ind = {list_smiles[i]:i for i in range(nM)}
train['indsmiles'] = train['SMILES'].apply(lambda x:dict_smiles2ind[x] )
val['indsmiles'] = val['SMILES'].apply(lambda x: dict_smiles2ind[x])
test['indsmiles'] = test['SMILES'].apply(lambda x: dict_smiles2ind[x])

number of different smiles (mol): 143206

[9]:

# molecule kernel_first step : compute Morgan FP for each smiles of all the dataset
MorganFP = komet.Morgan_FP(list_smiles)

Question 5 What are the dimensions of MorganFP ? Plot Morgan FP for the first molecule.

For more details on Morgan FP.

[ ]:



[10]:

# Choice of the parameters for the Nystrom approximation and the reduction dimension of the features
mM = 3000 #all mol to compute the mol kernel for medium-scale database
dM = 1000 #all dim for the mol features for medium-scale database

# In case there are less molecules than the number of molecules to compute the Nystrom approximation
mM = min(mM,nM) # number of molecule to compute nystrom
dM = min(dM,nM) # final dimension of features for molecules
print("mM",mM,"dM",dM)

mM 3000 dM 1000

[11]:

# compute the Nystrom approximation of the mol kernel and the features of the Kronecker kernel (features normalized and calculated on all mol contained in the dataset (train/val/test))
X_cn = komet.Nystrom_X_cn(mM,dM,nM,MorganFP,n_max = max(train['indsmiles']))
print("mol features shape",X_cn.shape)

mol kernel shape torch.Size([3000, 143206])
mol features shape torch.Size([143206, 1000])

Question 6

  1. Understand the different steps of the function Nystrom_X_cn, in particular explain the formula for Tanimoto kernel.

  2. What are the dimensions of X_cn ? Plot features for the first molecule.

[ ]:

3. Calculation of protein features

[12]:

# Index of the protein in the dataset
fasta = full[['Target Sequence']].drop_duplicates().values.flatten() # fasta sequence on the dataset, in the same order as the dataset
print("number of different Fasta (protein):",len(fasta))
# add ind_fasta dans train, val et test
train['indfasta'] = train['Target Sequence'].apply(lambda x: np.where(fasta==x)[0][0])
val['indfasta'] = val['Target Sequence'].apply(lambda x: np.where(fasta==x)[0][0])
test['indfasta'] = test['Target Sequence'].apply(lambda x:  np.where(fasta==x)[0][0])

number of different Fasta (protein): 2069

[13]:

# Load Protein kernel and dictionary of index
dict_ind2fasta_LCIdb = pickle.load(open(dataset_dir + "dict_ind2fasta_LCIdb.data", 'rb'))
dict_fasta2ind_LCIdb = {fasta:ind for ind,fasta in dict_ind2fasta_LCIdb.items()}
with open(dataset_dir + "dict_ind2fasta_LCIdb_K_prot.data", 'rb') as f:
    KP_LCIdb = pickle.load(f)

[14]:

# Protein kernel for the dataset
I_fasta = [int(dict_fasta2ind_LCIdb[fasta[i]]) for i in range(len(fasta))] # index of fasta in the precomputed dict and protein kernel, in the same order as the dataset
KP = KP_LCIdb[I_fasta,:][:,I_fasta]
KP = torch.tensor(KP, dtype=mytype).to(device)

[15]:

# computation of feature for protein (no nystrom, just SVD)
rP = KP.shape[0]#min(KP.shape[0],500)
U, Lambda, VT = torch.svd(KP)
Y = U[:,:rP] @ torch.diag(torch.sqrt(Lambda[:rP]))

# nomramlisation of the features
n_max = max(train['indfasta'])
Y_c = Y - Y[:n_max,:].mean(axis = 0)
Y_cn = Y_c / torch.norm(Y_c,dim = 1)[:,None]

Question 7

  1. What are the dimensions of KP? and Y_cn?

  2. Can we easily modify the dimensions of protein features ?

[ ]:

Load interactions index for train and test

[16]:

# Train
I, J, y = komet.load_datas(train)
n = len(I)
print("len(train)",n)

len(train) 236396

Question 7 What are I, J, y ?

[ ]:



[17]:

# Test
I_test, J_test, y_test = komet.load_datas(test)
n_test = len(I_test)
print("len(test)",n_test)

len(test) 44985

Question 8 This dataset is named Orphan because none of the molecules present in a pair of the train, are present in the test, same for proteins. Can you check it ?

[ ]:

Training/Testing with a choosen lambda

[18]:

#### TRAINING ####
lamb = 1e-6
# train the model
w_bfgs,b_bfgs,history_lbfgs_SVM = komet.SVM_bfgs(X_cn,Y_cn,y,I,J,lamb,niter=50)
# compute a probability using weights (Platt scaling)
s,t,history_lbfgs_Platt = komet.compute_proba_Platt_Scalling(w_bfgs,X_cn,Y_cn,y,I,J,niter=20)

L-BFGS time: 463.6205 seconds

Question 9

  1. Explains the different steps of the function SVM_bfgs.

  2. Plot history_lbfgs_SVM and explain the graphic.

[ ]:



[19]:

#### TRAIN ####
# we compute a probability using weights (Platt scaling)
m_train,y_pred_train, proba_pred_train = komet.compute_proba(w_bfgs,b_bfgs,s,t,X_cn,Y_cn,I,J)
# we compute the results
acc1_train,au_Roc_train,au_PR_train,thred_optim_train,acc_best_train,cm_train,FP_train = komet.results(y.cpu(),y_pred_train.cpu(),proba_pred_train.cpu())
print(f"roc AUC = {au_Roc_train:.4f}")
print(f"AUPR = {au_PR_train:.4f}")
print(f"accuracy (threshold 0.5)= {acc1_train:.4f}")
print(f"best threshold = {thred_optim_train:.4f}")
print(f"accuracy (best threshold)= {acc_best_train:.4f}")
print(f"false positive (best threshold)= {FP_train:.4f}")

roc AUC = 0.9997
AUPR = 0.9995
accuracy (threshold 0.5)= 0.9977
best threshold = 0.4499
accuracy (best threshold)= 0.9977
false positive (best threshold)= 0.0030

[20]:

#### TEST ####
# we compute a probability using weights (Platt scaling)
m,y_pred, proba_pred = komet.compute_proba(w_bfgs,b_bfgs,s,t,X_cn,Y_cn,I_test,J_test)
# we compute the results
acc1,au_Roc,au_PR,thred_optim,acc_best,cm,FP = komet.results(y_test.cpu(),y_pred.cpu(),proba_pred.cpu())
print(f"roc AUC = {au_Roc:.4f}")
print(f"AUPR = {au_PR:.4f}")
print(f"accuracy (threshold 0.5)= {acc1:.4f}")
print(f"best threshold = {thred_optim:.4f}")
print(f"accuracy (best threshold)= {acc_best:.4f}")
print(f"false positive (best threshold)= {FP:.4f}")

roc AUC = 0.8818
AUPR = 0.8971
accuracy (threshold 0.5)= 0.6867
best threshold = 0.0009
accuracy (best threshold)= 0.8145
false positive (best threshold)= 0.1447

[21]:

# plot distribution (density) of p when y_test=1
plt.hist(proba_pred.cpu().numpy()[y_test.cpu().numpy()==1],bins=100,alpha=0.8,color='green',label='y_test=1');
plt.hist(proba_pred.cpu().numpy()[y_test.cpu()==-1],bins=100,alpha=0.5,color='red',label='y_test=-1');
plt.legend()
plt.title('Distribution of predicted probability')
plt.show()
../_images/vignettes_komet_TP_49_0.png 
[22]:

# plot confusion matrix
from sklearn.metrics import ConfusionMatrixDisplay
labels = [-1., 1.]
disp = ConfusionMatrixDisplay(confusion_matrix=cm,
                              display_labels=labels)
disp.plot()
plt.show()
../_images/vignettes_komet_TP_50_0.png

Question 10 Comments the code and different results.

[ ]:

Question 11

  1. Is the aspirin molecule in train/val/test ?

  2. Train the model with the full dataset.

  3. Predict for the aspirin molecule, protein(s) which bind with.

  4. Compare with Targets of the Drugbank database.

[ ]: