Predict drugs’ protein interaction profiles or proteins’ molecular interaction profiles with Komet and LCIdb
Downloading utils from Github
Downloading dataset from Zenodo
Predict proteins’ molecular interaction profiles with Komet and LCIdb
Predict drugs’ protein interaction profiles with Komet and LCIdb
[1]:
%load_ext autoreload
%autoreload 2
import torch
import torch.optim as optim
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
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
[2]:
!pip install rdkit
Collecting rdkit
Downloading rdkit-2023.9.5-cp310-cp310-manylinux_2_17_x86_64.manylinux2014_x86_64.whl (34.4 MB)
━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 34.4/34.4 MB 25.9 MB/s eta 0:00:00
Requirement already satisfied: numpy in /usr/local/lib/python3.10/dist-packages (from rdkit) (1.25.2)
Requirement already satisfied: Pillow in /usr/local/lib/python3.10/dist-packages (from rdkit) (9.4.0)
Installing collected packages: rdkit
Successfully installed rdkit-2023.9.5
1. Download the data from a GitHub repo.
[2]:
!wget -q https://github.com/Guichaoua/komet/raw/main/komet/komet.py
zsh:1: command not found: wget
[3]:
import komet
cpu
2. Download the data from Zenodo.
You can directly download train_arr.pkl or make your training dataset.
[5]:
LCIdb_path = 'data/LCIdb_v2.csv'
LCIdb = pd.read_csv(LCIdb_path,low_memory=False)
LCIdb.head()
[5]:
| smiles | fasta | ChEMBL ID | PubChem ID | IUPHAR ID | Ligand names | Target | uniprot | mean | mean pIC50 | ... | mean pKd | min | min pIC50 | min pKi | min pKd | max | max pIC50 | max pKi | max pKd | score | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | BrC(=C\c1ccccc1)/C=N/n1cnnc1 | MNPTLILAAFCLGIASATLTFDHSLEAQWTKWKAMHNRLYGMNEEG... | CHEMBL3190095 | 6861939.0 | NaN | sid24812872 | ctsl | P07711 | 4.2 | 4.2 | ... | NaN | 4.2 | 4.2 | NaN | NaN | 4.2 | 4.2 | NaN | NaN | 0.5 |
| 1 | BrC(Cn1ncc2c(N3CCCC3)ncnc21)c1ccccc1 | MGSNKSKPKDASQRRRSLEPAENVHGAGGGAFPASQTPSKPASADG... | CHEMBL382012 | 11617657.0 | NaN | 1-(2-bromo-2-phenylethyl)-4-(pyrrolidin-1-yl)-... | src | P12931 | 5.0 | NaN | ... | NaN | 5.0 | NaN | 5.0 | NaN | 5.0 | NaN | 5.0 | NaN | 0.5 |
| 2 | BrC(Cn1ncc2c(N3CCCCC3)ncnc21)c1ccccc1 | MGSNKSKPKDASQRRRSLEPAENVHGAGGGAFPASQTPSKPASADG... | CHEMBL203922 | 11696703.0 | NaN | 1-(2-bromo-2-phenylethyl)-4-(piperidin-1-yl)-1... | src | P12931 | 5.8 | NaN | ... | NaN | 5.8 | NaN | 5.8 | NaN | 5.8 | NaN | 5.8 | NaN | 0.5 |
| 3 | BrC(Cn1ncc2c(NCc3ccccc3)ncnc21)c1ccccc1 | MGSNKSKPKDASQRRRSLEPAENVHGAGGGAFPASQTPSKPASADG... | CHEMBL382153 | 11545984.0 | NaN | n-benzyl-1-(2-bromo-2-phenylethyl)-1h-pyrazolo... | src | P12931 | 5.5 | NaN | ... | NaN | 5.5 | NaN | 5.5 | NaN | 5.5 | NaN | 5.5 | NaN | 0.5 |
| 4 | BrC1=CC2CNCCC2S1 | MSGADRSPNAGAAPDSAPGQAAVASAYQRFEPRAYLRNNYAPPRGD... | CHEMBL398475 | 44447166.0 | NaN | 2-bromo-3a,4,5,6,7,7a-hexahydrothieno[3,2-c]py... | pnmt | P11086 | 5.9 | NaN | ... | NaN | 5.9 | NaN | 5.9 | NaN | 5.9 | NaN | 5.9 | NaN | 0.5 |
5 rows × 21 columns
[6]:
train = komet.make_trains_full('data/LCIdb_v2.csv','/Users/gguichaoua/Dropbox/gwenn/these/2024_01_Komet_JCIM/Zenodo_v2/',nb_trains=5)
nombre de smiles: 271180
nombre de fasta: 2060
/Users/gguichaoua/Dropbox/gwenn/these/chemogenomique/komet/docs/source/vignettes/komet.py:2619: FutureWarning: In a future version of pandas all arguments of DataFrame.drop except for the argument 'labels' will be keyword-only.
df = df.drop(indsmiles_index_with_nan,0)
/Users/gguichaoua/Dropbox/gwenn/these/chemogenomique/komet/docs/source/vignettes/komet.py:2622: FutureWarning: In a future version of pandas all arguments of DataFrame.drop except for the argument 'labels' will be keyword-only.
df = df.drop(indfasta_index_with_nan,0)
matrice d'interactions: (2060, 271180)
train (793596, 3)
train (793596, 3)
train (793596, 3)
train (793596, 3)
train (793596, 3)
[7]:
train_0 = train[0]
train_0.head()
[7]:
| SMILES | Target Sequence | Label | |
|---|---|---|---|
| 0 | Brc1cccc(-n2cnc(-c3ccccc3Nc3ccncc3)n2)c1 | MQSKVLLAVALWLCVETRAASVGLPSVSLDLPRLSIQKDILTIKAN... | 1 |
| 1 | Brc1cccc(-n2cnc(-c3ccccc3Nc3ccncc3)n2)c1 | MGMACLTMTEMEGTSTSSIYQNGDISGNANSMKQIDPVLQVYLYHS... | 0 |
| 2 | Brc1cccc(-n2cnc(-c3ccccc3Nc3ccncc3)n2)c1 | MVSYWDTGVLLCALLSCLLLTGSSSGSKLKDPELSLKGTQHIMQAG... | 1 |
| 3 | Brc1cccc(Nc2ncc(-c3ccccc3)s2)n1 | MQSKVLLAVALWLCVETRAASVGLPSVSLDLPRLSIQKDILTIKAN... | 1 |
| 4 | Brc1cccc(Oc2ccc3c(/C=C/c4ccccc4)n[nH]c3c2)c1 | MQSKVLLAVALWLCVETRAASVGLPSVSLDLPRLSIQKDILTIKAN... | 1 |
3. Predict proteins’ molecular interaction profiles with Komet and LCIdb
You need to download to Zenodo and save in ./data, the 3 files dict_ind2fasta_all.data, dict_ind2fasta_all_K_prot.data and fasta_uniprot_family.csv
[72]:
name_drug = 'aspirin'
smiles_drug = 'CC1=CC(=O)C2=C(C1=O)C=CC=C2'
df_pred_drug,pred = komet.predict_drug_profile(train,smiles_drug,mM = 3000,dM = 1000,lamb = 1e-6)
df_pred_drug.head()
The drug is not in the dataset
kernel prot shape torch.Size([2060, 2060])
protein features shape torch.Size([2060, 2060])
mol kernel shape torch.Size([3000, 271181])
L-BFGS time: 622.5025 seconds
mol kernel shape torch.Size([3000, 271181])
L-BFGS time: 624.5909 seconds
mol kernel shape torch.Size([3000, 271181])
L-BFGS time: 620.6984 seconds
mol kernel shape torch.Size([3000, 271181])
L-BFGS time: 631.1737 seconds
mol kernel shape torch.Size([3000, 271181])
L-BFGS time: 637.6766 seconds
[72]:
| Target | uniprot | family | Label_known | Proba_predicted_mean | Proba_predicted_std | Proba_predicted_min | Proba_predicted_max | |
|---|---|---|---|---|---|---|---|---|
| 0 | ces1 | P23141 | type-B carboxylesterase/lipase family | -1 | 0.999757 | 0.000130 | 0.999517 | 0.999869 |
| 1 | ces2 | O00748 | type-B carboxylesterase/lipase family | -1 | 0.996522 | 0.002068 | 0.992532 | 0.998258 |
| 2 | maoa | P21397 | flavin monoamine oxidase family | -1 | 0.994963 | 0.003355 | 0.988908 | 0.998240 |
| 3 | maob | P27338 | flavin monoamine oxidase family | -1 | 0.991253 | 0.004065 | 0.983757 | 0.996147 |
| 4 | prkch | C0HM02 | Not known | -1 | 0.983708 | 0.010395 | 0.966951 | 0.995497 |
Remarks :
In Drugbank, for aspirin, we can see 16 targets but only 2 bind to aspirin (ptgs1 and ptgs2). For others, the pharmacological action is unknown. Our prediction finds these 2 connections and rejects the others.
[83]:
df_pred_drug[df_pred_drug['uniprot'].isin(['P23219','P35354'])]
[83]:
| Target | uniprot | family | Label_known | Proba_predicted_mean | Proba_predicted_std | Proba_predicted_min | Proba_predicted_max | |
|---|---|---|---|---|---|---|---|---|
| 12 | ptgs2 | P35354 | prostaglandin G/H synthase family | -1 | 0.902107 | 0.020121 | 0.884362 | 0.936080 |
| 24 | ptgs1 | P23219 | prostaglandin G/H synthase family | -1 | 0.710457 | 0.095295 | 0.590268 | 0.854953 |
[84]:
df_pred_drug[df_pred_drug['uniprot'].isin(['P23219','P35354','Q04828','Q13131','P54646','Q9Y478','O43741','P54619','Q9UGJ0','Q9UGI9','P25101','P04637','P11021','P51812', 'P25963', 'P98066' ,'P29466', 'P42574', 'O14920','P28482','Q8TD08','P27361','P31152','Q16659','Q13164','P24385','P01106','P12004','Q99519'])]
[84]:
| Target | uniprot | family | Label_known | Proba_predicted_mean | Proba_predicted_std | Proba_predicted_min | Proba_predicted_max | |
|---|---|---|---|---|---|---|---|---|
| 12 | ptgs2 | P35354 | prostaglandin G/H synthase family | -1 | 0.902107 | 0.020121 | 0.884362 | 0.936080 |
| 24 | ptgs1 | P23219 | prostaglandin G/H synthase family | -1 | 0.710457 | 0.095295 | 0.590268 | 0.854953 |
| 163 | akr1c1 | Q04828 | aldo/keto reductase family | -1 | 0.120188 | 0.101206 | 0.014647 | 0.294680 |
| 219 | casp3 | P42574 | peptidase family | -1 | 0.066869 | 0.045929 | 0.018455 | 0.142134 |
| 779 | ednra | P25101 | G-protein coupled receptor 1 family | -1 | 0.005633 | 0.006959 | 0.000550 | 0.019016 |
| 838 | mapk1 | P28482 | protein kinase superfamily | -1 | 0.004710 | 0.008011 | 0.000240 | 0.020719 |
| 849 | neu1 | Q99519 | glycosyl hydrolase 33 family | -1 | 0.004558 | 0.000592 | 0.003846 | 0.005532 |
| 884 | casp1 | P29466 | peptidase family | -1 | 0.004312 | 0.001560 | 0.002841 | 0.006690 |
| 1006 | myc | P01106 | Not known | -1 | 0.003082 | 0.000511 | 0.002452 | 0.003981 |
| 1033 | ccnd1 | P24385 | cyclin family | -1 | 0.002879 | 0.004249 | 0.000073 | 0.011333 |
| 1092 | prkab1 | Q9Y478 | 5'-AMP-activated protein kinase beta subunit f... | -1 | 0.002471 | 0.001259 | 0.001092 | 0.004723 |
| 1194 | mapk3 | P27361 | protein kinase superfamily | -1 | 0.001931 | 0.003030 | 0.000172 | 0.007980 |
| 1228 | prkab2 | O43741 | 5'-AMP-activated protein kinase beta subunit f... | -1 | 0.001795 | 0.000655 | 0.001019 | 0.002916 |
| 1361 | mapk15 | Q8TD08 | protein kinase superfamily | -1 | 0.001129 | 0.001246 | 0.000206 | 0.003584 |
| 1364 | prkag3 | Q9UGI9 | 5'-AMP-activated protein kinase gamma subunit ... | -1 | 0.001119 | 0.000392 | 0.000756 | 0.001847 |
| 1412 | prkag1 | P54619 | 5'-AMP-activated protein kinase gamma subunit ... | -1 | 0.000964 | 0.000485 | 0.000586 | 0.001857 |
| 1421 | prkag2 | Q9UGJ0 | 5'-AMP-activated protein kinase gamma subunit ... | -1 | 0.000927 | 0.000288 | 0.000653 | 0.001455 |
| 1437 | prkaa2 | P54646 | protein kinase superfamily | -1 | 0.000860 | 0.001134 | 0.000149 | 0.003118 |
| 1502 | tp53 | P04637 | p53 family | -1 | 0.000642 | 0.000262 | 0.000431 | 0.001142 |
| 1578 | mapk7 | Q13164 | protein kinase superfamily | -1 | 0.000457 | 0.000400 | 0.000093 | 0.001194 |
| 1720 | hspa5 | P11021 | heat shock protein 70 family | -1 | 0.000219 | 0.000040 | 0.000171 | 0.000289 |
| 1731 | prkaa1 | Q13131 | protein kinase superfamily | -1 | 0.000207 | 0.000280 | 0.000044 | 0.000766 |
| 1862 | ikbkb | O14920 | protein kinase superfamily | -1 | 0.000085 | 0.000113 | 0.000012 | 0.000310 |
| 1938 | rps6ka3 | P51812 | protein kinase superfamily | -1 | 0.000044 | 0.000067 | 0.000005 | 0.000177 |
4. Predict drugs’ protein interaction profiles with Komet and LCIdb
You need to download to Zenodo and save in ./data, the 3 files dict_ind2fasta_all.data, dict_ind2fasta_all_K_prot.data and fasta_uniprot_family.csv
[132]:
df_family = pd.read_csv("data/fasta_uniprot_family.csv")
[133]:
name_protein = 'ptgs1'
fasta_protein = df_family[df_family["Target"] == "ptgs1"]["fasta"].values[0]
df_pred_protein = komet.predict_protein_profile(train,fasta_protein,mM = 3000,dM = 1000,lamb = 1e-6)
df_pred_protein.head()
kernel prot shape torch.Size([2060, 2060])
protein features shape torch.Size([2060, 2060])
mol kernel shape torch.Size([3000, 271180])
L-BFGS time: 598.0121 seconds
mol kernel shape torch.Size([3000, 271180])
L-BFGS time: 589.7209 seconds
mol kernel shape torch.Size([3000, 271180])
L-BFGS time: 607.9472 seconds
mol kernel shape torch.Size([3000, 271180])
L-BFGS time: 591.4509 seconds
mol kernel shape torch.Size([3000, 271180])
L-BFGS time: 618.3177 seconds
[133]:
| SMILES | Label | Proba_predicted_mean | Proba_predicted_std | Proba_predicted_min | Proba_predicted_max | |
|---|---|---|---|---|---|---|
| 250550 | COc1ccc(-n2nc(CC(C(=O)O)c3ccccc3)cc2-c2ccc(Cl)... | -1.0 | 0.999925 | 0.000031 | 0.999883 | 0.999976 |
| 6791 | CC(C)(C)c1cc(C(=O)c2ccccc2)cc2c1OCC2(C)C | -1.0 | 0.999882 | 0.000064 | 0.999802 | 0.999956 |
| 250552 | COc1ccc(-n2nc(CC(C(=O)O)c3ccccc3OC)cc2-c2ccc(C... | -1.0 | 0.999855 | 0.000060 | 0.999773 | 0.999951 |
| 250542 | COc1ccc(-n2nc(CC(C(=O)O)c3ccc(Cl)cc3)cc2-c2ccc... | -1.0 | 0.999781 | 0.000079 | 0.999669 | 0.999913 |
| 250555 | COc1ccc(C(Cc2cc(-c3ccc(Cl)cc3)n(-c3ccc(OC)cc3)... | -1.0 | 0.999781 | 0.000079 | 0.999669 | 0.999913 |
[140]:
df_pred_protein.head(15)
[140]:
| index | SMILES | Label | Proba_predicted_mean | Proba_predicted_std | Proba_predicted_min | Proba_predicted_max | |
|---|---|---|---|---|---|---|---|
| 0 | 250550 | COc1ccc(-n2nc(CC(C(=O)O)c3ccccc3)cc2-c2ccc(Cl)... | -1.0 | 0.999925 | 0.000031 | 0.999883 | 0.999976 |
| 1 | 6791 | CC(C)(C)c1cc(C(=O)c2ccccc2)cc2c1OCC2(C)C | -1.0 | 0.999882 | 0.000064 | 0.999802 | 0.999956 |
| 2 | 250552 | COc1ccc(-n2nc(CC(C(=O)O)c3ccccc3OC)cc2-c2ccc(C... | -1.0 | 0.999855 | 0.000060 | 0.999773 | 0.999951 |
| 3 | 250542 | COc1ccc(-n2nc(CC(C(=O)O)c3ccc(Cl)cc3)cc2-c2ccc... | -1.0 | 0.999781 | 0.000079 | 0.999669 | 0.999913 |
| 4 | 250555 | COc1ccc(C(Cc2cc(-c3ccc(Cl)cc3)n(-c3ccc(OC)cc3)... | -1.0 | 0.999781 | 0.000079 | 0.999669 | 0.999913 |
| 5 | 8619 | COc1ccc2c(c1)c(CC(=O)O)c(C)n2C(=O)c1ccc(Cl)cc1 | -1.0 | 0.999779 | 0.000118 | 0.999595 | 0.999910 |
| 6 | 250549 | COc1ccc(-n2nc(CC(C(=O)O)c3ccccc3)cc2-c2ccc(C)c... | -1.0 | 0.999701 | 0.000266 | 0.999175 | 0.999877 |
| 7 | 250540 | COc1ccc(-n2nc(CC(C(=O)O)c3ccc(C)cc3)cc2-c2ccc(... | -1.0 | 0.999654 | 0.000119 | 0.999457 | 0.999832 |
| 8 | 250541 | COc1ccc(-n2nc(CC(C(=O)O)c3ccc(Cl)cc3)cc2-c2ccc... | -1.0 | 0.999654 | 0.000119 | 0.999457 | 0.999832 |
| 9 | 250544 | COc1ccc(-n2nc(CC(C(=O)O)c3cccc(C)c3)cc2-c2ccc(... | -1.0 | 0.999624 | 0.000180 | 0.999328 | 0.999848 |
| 10 | 250548 | COc1ccc(-n2nc(CC(C(=O)O)c3cccc(OC)c3)cc2-c2ccc... | -1.0 | 0.999605 | 0.000160 | 0.999359 | 0.999774 |
| 11 | 140471 | Cc1c(CC(=O)O)c2ccccc2n1C(=O)c1ccc(Cl)cc1 | 1.0 | 0.999450 | 0.000264 | 0.999154 | 0.999857 |
| 12 | 210947 | COc1ccc2c(c1)c(CC(=O)Oc1cccnc1)c(C)n2C(=O)c1cc... | -1.0 | 0.999449 | 0.000210 | 0.999166 | 0.999669 |
| 13 | 250551 | COc1ccc(-n2nc(CC(C(=O)O)c3ccccc3OC)cc2-c2ccc(C... | -1.0 | 0.999325 | 0.000578 | 0.998312 | 0.999904 |
| 14 | 210925 | COc1ccc(OC(=O)Cc2c(C)n(C(=O)c3ccc(Cl)cc3)c3ccc... | -1.0 | 0.999261 | 0.000204 | 0.999020 | 0.999576 |
[141]:
# the protein is in the dataset : 1 (DTI+), 0 (DTI-), -1 (unknown DTI)
df_pred_protein.value_counts('Label')
[141]:
Label
-1.0 271032
0.0 74
1.0 74
dtype: int64
[143]:
for i in range(15):
print("rank",i+1,df_pred_protein['SMILES'].values[i],df_pred_protein['Proba_predicted_mean'].values[i])
rank 1 COc1ccc(-n2nc(CC(C(=O)O)c3ccccc3)cc2-c2ccc(Cl)cc2)cc1 0.9999254941940308
rank 2 CC(C)(C)c1cc(C(=O)c2ccccc2)cc2c1OCC2(C)C 0.999882173538208
rank 3 COc1ccc(-n2nc(CC(C(=O)O)c3ccccc3OC)cc2-c2ccc(Cl)cc2)cc1 0.9998553276062012
rank 4 COc1ccc(-n2nc(CC(C(=O)O)c3ccc(Cl)cc3)cc2-c2ccc(Cl)cc2)cc1 0.9997808814048768
rank 5 COc1ccc(C(Cc2cc(-c3ccc(Cl)cc3)n(-c3ccc(OC)cc3)n2)C(=O)O)cc1 0.9997808814048768
rank 6 COc1ccc2c(c1)c(CC(=O)O)c(C)n2C(=O)c1ccc(Cl)cc1 0.9997786998748779
rank 7 COc1ccc(-n2nc(CC(C(=O)O)c3ccccc3)cc2-c2ccc(C)cc2)cc1 0.9997011065483093
rank 8 COc1ccc(-n2nc(CC(C(=O)O)c3ccc(C)cc3)cc2-c2ccc(Cl)cc2)cc1 0.999653947353363
rank 9 COc1ccc(-n2nc(CC(C(=O)O)c3ccc(Cl)cc3)cc2-c2ccc(C)cc2)cc1 0.999653947353363
rank 10 COc1ccc(-n2nc(CC(C(=O)O)c3cccc(C)c3)cc2-c2ccc(Cl)cc2)cc1 0.9996240496635437
rank 11 COc1ccc(-n2nc(CC(C(=O)O)c3cccc(OC)c3)cc2-c2ccc(Cl)cc2)cc1 0.9996052384376526
rank 12 Cc1c(CC(=O)O)c2ccccc2n1C(=O)c1ccc(Cl)cc1 0.9994504809379577
rank 13 COc1ccc2c(c1)c(CC(=O)Oc1cccnc1)c(C)n2C(=O)c1ccc(Cl)cc1 0.9994487404823303
rank 14 COc1ccc(-n2nc(CC(C(=O)O)c3ccccc3OC)cc2-c2ccc(C)cc2)cc1 0.9993251919746399
rank 15 COc1ccc(OC(=O)Cc2c(C)n(C(=O)c3ccc(Cl)cc3)c3ccc(OC)cc23)cc1 0.9992609858512879
Remarks :
For the predictions :
For the second drug CHEMBL275120, the binding was found: https://pubchem.ncbi.nlm.nih.gov/compound/10852345#section=Biological-Test-Results
For the 6th drug (Indomethacin), the binding was found: https://pubchem.ncbi.nlm.nih.gov/compound/3715#section=BioAssay-Results
For the 12th durg (CHEMBL1241992), the binding was known