{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.10.13","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"none","dataSources":[{"sourceId":67356,"databundleVersionId":8006601,"sourceType":"competition"},{"sourceId":8042988,"sourceType":"datasetVersion","datasetId":4740586}],"dockerImageVersionId":30684,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"This notebook is intended to explore the necessary domain of applicability of our models.\n\nOften in medicinal chemistry, a core will stay constant (or at least similar) and the groups appended to the core will be different. We often call compounds having the same core (or pharmacophore once you understand the chemical matter better) as a \"chemical series\".\n\nSee below for the competition organizers description of the data.\n\n\"[train/test].[csv/parquet] - The train or test data, available in both the csv and parquet formats.\n\n...\n\nmolecule_smiles - The structure of the fully assembled molecule, in SMILES. This includes the **three building blocks and the triazine core**. Note we use a [Dy] as the stand-in for the DNA linker.\n\n...\"\n\nNote the fact that this description mentioned a triazine core for molecule SMILES. This implies that the train and the test set will be molecules where each building block reacts with a triazine (this is classic DEL chemistry) where the triazine acts as the central core of the molecule. \n\nIf the entire train and test set are from a \"triazine chemical series\", then the competition organizers would be asking us to predict the binary binding of molecules to one of three proteins when new building blocks are used, but the core is the same. This in itself could be a very useful model (particularly as it would eventually becomes a binding affinity regressor when sufficient binding affinity data was obtained) and it wouldn't be uncommon to see a model like this used in practice on a medicinal chemistry program where the hits were identified from a DEL screen.\n\nHowever, it is often the case that when optimizing a chemical series that some issues associated with the core arises. Perhaps a tox liability is identified, or there's a problem with the physical properties that are inherent to the core. In this instance, medicinal chemists will undergo an excercise we call \"scaffold hopping\". In \"scaffold hopping\" the core of the molecule is swapped out for alternatives in the hopes of identifying a new \"chemical series\" that avoids the limitations of it's predecessor. This task can be rather challenging, but it is the \"holy grail\" (in my personal opinion) of any quantitative structure activity relationship (commonly called QSAR) model. To be clear: switching out the core of the molecule is often much more difficult than switching out the R-groups on a molecule.\n\nTo summarize: binding affinity is often a function of the core + R groups, where R groups are the parts attached to the core. We know that we have many R groups (building blocks) in this dataset, but how many cores do we have? The answer to this question will inform us whether or not we're being asked to do scaffold hopping, or if we're building a model for a triazine chemical series.","metadata":{}},{"cell_type":"code","source":"import numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\nimport os","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2024-04-12T02:20:28.634041Z","iopub.execute_input":"2024-04-12T02:20:28.634600Z","iopub.status.idle":"2024-04-12T02:20:28.641450Z","shell.execute_reply.started":"2024-04-12T02:20:28.634551Z","shell.execute_reply":"2024-04-12T02:20:28.639508Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"!pip install mapply","metadata":{"execution":{"iopub.status.busy":"2024-04-12T02:20:31.325255Z","iopub.execute_input":"2024-04-12T02:20:31.325629Z","iopub.status.idle":"2024-04-12T02:20:47.287998Z","shell.execute_reply.started":"2024-04-12T02:20:31.325600Z","shell.execute_reply":"2024-04-12T02:20:47.286849Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import mapply","metadata":{"execution":{"iopub.status.busy":"2024-04-12T02:20:51.482358Z","iopub.execute_input":"2024-04-12T02:20:51.483210Z","iopub.status.idle":"2024-04-12T02:20:51.693984Z","shell.execute_reply.started":"2024-04-12T02:20:51.483167Z","shell.execute_reply":"2024-04-12T02:20:51.692999Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"!pip install rdkit","metadata":{"execution":{"iopub.status.busy":"2024-04-12T02:20:52.813172Z","iopub.execute_input":"2024-04-12T02:20:52.814296Z","iopub.status.idle":"2024-04-12T02:21:09.966306Z","shell.execute_reply.started":"2024-04-12T02:20:52.814259Z","shell.execute_reply":"2024-04-12T02:21:09.964792Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#RDKit is the cheminformatics package that will help us with substructure matching\nfrom rdkit import Chem\nfrom rdkit.Chem import Draw","metadata":{"execution":{"iopub.status.busy":"2024-04-12T02:25:53.367043Z","iopub.execute_input":"2024-04-12T02:25:53.367471Z","iopub.status.idle":"2024-04-12T02:25:53.372548Z","shell.execute_reply.started":"2024-04-12T02:25:53.367442Z","shell.execute_reply":"2024-04-12T02:25:53.371448Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Let's read in the test dataset. Special thanks to GreySnow for making this lightweight test.csv!\nSee: https://www.kaggle.com/code/shlomoron/belka-shrinking-the-dataset & https://www.kaggle.com/code/shlomoron/belka-shrunken-train-set-loading","metadata":{}},{"cell_type":"code","source":"dtypes_test = {'buildingblock1_smiles': np.int16, 'buildingblock2_smiles': np.int16, 'buildingblock3_smiles': np.int16,\n          'is_BRD4':np.byte, 'is_HSA':np.byte, 'is_sEH':np.byte}\n\ntest = pd.read_csv('/kaggle/input/belka-shrunken-train-set/test.csv', dtype = dtypes_test)","metadata":{"execution":{"iopub.status.busy":"2024-04-12T02:21:10.139801Z","iopub.execute_input":"2024-04-12T02:21:10.140250Z","iopub.status.idle":"2024-04-12T02:21:13.432358Z","shell.execute_reply.started":"2024-04-12T02:21:10.140208Z","shell.execute_reply":"2024-04-12T02:21:13.431415Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#Here is a free website where you can generate SMILES or SMARTS strings: https://pubchem.ncbi.nlm.nih.gov//edit3/index.html\n#Let's take a look at what a triazine looks like.\ntriazine = Chem.MolFromSmiles('C1=NC=NC=N1')\ntriazine","metadata":{"execution":{"iopub.status.busy":"2024-04-12T02:21:13.434554Z","iopub.execute_input":"2024-04-12T02:21:13.434899Z","iopub.status.idle":"2024-04-12T02:21:13.451458Z","shell.execute_reply.started":"2024-04-12T02:21:13.434864Z","shell.execute_reply":"2024-04-12T02:21:13.450486Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#Checking to make sure our matching works and highlighting the triazine in the core of the molecule.\nmol = Chem.MolFromSmiles(test['molecule_smiles'].iloc[0])\nprint(mol.HasSubstructMatch(triazine))\nmol.GetSubstructMatch(triazine)\nmol","metadata":{"execution":{"iopub.status.busy":"2024-04-12T02:29:22.158273Z","iopub.execute_input":"2024-04-12T02:29:22.159162Z","iopub.status.idle":"2024-04-12T02:29:22.179979Z","shell.execute_reply.started":"2024-04-12T02:29:22.159120Z","shell.execute_reply":"2024-04-12T02:29:22.178833Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#Thanks to Hengck23 for the tip about mapply that allows this analysis to run in a Kaggle kernel!\n#See: https://www.kaggle.com/competitions/leash-BELKA/discussion/492846\nmapply.init(\n    n_workers=-1,\n    progressbar=True,\n)","metadata":{"execution":{"iopub.status.busy":"2024-04-12T02:22:10.358766Z","iopub.execute_input":"2024-04-12T02:22:10.359843Z","iopub.status.idle":"2024-04-12T02:22:10.365604Z","shell.execute_reply.started":"2024-04-12T02:22:10.359797Z","shell.execute_reply":"2024-04-12T02:22:10.364276Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#A simple function that converts a SMILES string to a RDKit molecule object, then checks if it has a triazine substructure match.\ndef check_for_triazine(x):\n    check = Chem.MolFromSmiles(x).HasSubstructMatch(triazine)\n    return check","metadata":{"execution":{"iopub.status.busy":"2024-04-12T02:22:12.042012Z","iopub.execute_input":"2024-04-12T02:22:12.042459Z","iopub.status.idle":"2024-04-12T02:22:12.049169Z","shell.execute_reply.started":"2024-04-12T02:22:12.042417Z","shell.execute_reply":"2024-04-12T02:22:12.047383Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#Let's create a column that encodes whether the molecule for that row has a triazine core or not.\n#This is going to take some time...\ntest['triazine'] = test['molecule_smiles'].mapply(check_for_triazine)","metadata":{"execution":{"iopub.status.busy":"2024-04-12T02:22:13.494526Z","iopub.execute_input":"2024-04-12T02:22:13.494964Z","iopub.status.idle":"2024-04-12T02:24:29.335003Z","shell.execute_reply.started":"2024-04-12T02:22:13.494930Z","shell.execute_reply":"2024-04-12T02:24:29.333806Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#Check to see if we have any rows that are not triazines in the dataset.\ntest['triazine'].unique()","metadata":{"execution":{"iopub.status.busy":"2024-04-12T02:24:29.356166Z","iopub.execute_input":"2024-04-12T02:24:29.356540Z","iopub.status.idle":"2024-04-12T02:24:29.368871Z","shell.execute_reply.started":"2024-04-12T02:24:29.356510Z","shell.execute_reply":"2024-04-12T02:24:29.367623Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#Get a dataframe where there are no triazines in the data.\nnot_triazines = test[test['triazine']==False]","metadata":{"execution":{"iopub.status.busy":"2024-04-12T02:26:54.868809Z","iopub.execute_input":"2024-04-12T02:26:54.869486Z","iopub.status.idle":"2024-04-12T02:26:54.916960Z","shell.execute_reply.started":"2024-04-12T02:26:54.869453Z","shell.execute_reply":"2024-04-12T02:26:54.915720Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"percent_new_scaffolds = 100*len(not_triazines)/len(test)\nprint(f\"Approximately {round(percent_new_scaffolds, 2)} percent of the test set is asking us to predict whether or not non-triazine molecules bind our proteins of interest\")","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#Let's see what some of the non-triazine containing molecules in the test set look like.\nDraw.MolsToGridImage([Chem.MolFromSmiles(x) for x in not_triazines['molecule_smiles'].sample(12)], molsPerRow=4, subImgSize=(400,300))","metadata":{"execution":{"iopub.status.busy":"2024-04-12T02:27:28.696244Z","iopub.execute_input":"2024-04-12T02:27:28.697176Z","iopub.status.idle":"2024-04-12T02:27:28.856245Z","shell.execute_reply.started":"2024-04-12T02:27:28.697132Z","shell.execute_reply":"2024-04-12T02:27:28.855057Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_dtypes = {'buildingblock1_smiles': np.int16, 'buildingblock2_smiles': np.int16, 'buildingblock3_smiles': np.int16,\n          'binds_BRD4':np.byte, 'binds_HSA':np.byte, 'binds_sEH':np.byte}\n\ntrain = pd.read_csv('/kaggle/input/belka-shrunken-train-set/train.csv', dtype = train_dtypes)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#Let's create a column that encodes whether the molecule for that row has a triazine core or not.\n#This is going to take some time...\ntrain['triazine'] = train['molecule_smiles'].mapply(check_for_triazine)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#Let's see if there are any non-triazine containing molecules in the train set (i.e. if there are any rows marked False)\ntrain['triazine'].unique()","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"As we can see the train data only has triazine cores, but the test data has other unique cores. Building validation sets where we identify how well we are doing at scaffold hopping may be a bit tricky, since the train set only has triazine cores for us to work with.\n\nEarly in the competition, Hengck23 found some data from another paper. See: https://www.kaggle.com/competitions/leash-BELKA/discussion/491908\n\nI was able to identify the DNA attachment points and convert the SMILES from this dataset into SMILES that are consistent with our Kaggle data (i.e. the DNA attachment point is signified with [Dy]). See: https://www.kaggle.com/code/chemdatafarmer/additional-seh-data\n\nBest of luck in this competition!","metadata":{}}]}