{"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"}],"dockerImageVersionId":30698,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# This notebook is intended to demonstrate some useful cheminformatic transformations to help you with handling the Dy in the molecules.","metadata":{}},{"cell_type":"code","source":"import pandas as pd\nimport numpy as np\nimport os","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2024-04-26T00:43:53.201083Z","iopub.execute_input":"2024-04-26T00:43:53.201604Z","iopub.status.idle":"2024-04-26T00:43:54.618964Z","shell.execute_reply.started":"2024-04-26T00:43:53.201550Z","shell.execute_reply":"2024-04-26T00:43:54.615991Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"!pip install mapply","metadata":{"_kg_hide-output":true,"execution":{"iopub.status.busy":"2024-04-26T00:43:54.621974Z","iopub.execute_input":"2024-04-26T00:43:54.622527Z","iopub.status.idle":"2024-04-26T00:44:12.342451Z","shell.execute_reply.started":"2024-04-26T00:43:54.622494Z","shell.execute_reply":"2024-04-26T00:44:12.340723Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import mapply","metadata":{"execution":{"iopub.status.busy":"2024-04-26T00:44:12.345407Z","iopub.execute_input":"2024-04-26T00:44:12.345853Z","iopub.status.idle":"2024-04-26T00:44:12.612166Z","shell.execute_reply.started":"2024-04-26T00:44:12.345817Z","shell.execute_reply":"2024-04-26T00:44:12.611060Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"!pip install rdkit","metadata":{"_kg_hide-output":true,"execution":{"iopub.status.busy":"2024-04-26T00:44:12.615114Z","iopub.execute_input":"2024-04-26T00:44:12.616951Z","iopub.status.idle":"2024-04-26T00:44:30.958761Z","shell.execute_reply.started":"2024-04-26T00:44:12.616907Z","shell.execute_reply":"2024-04-26T00:44:30.957551Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from rdkit import Chem\nfrom rdkit.Chem import Draw, AllChem","metadata":{"execution":{"iopub.status.busy":"2024-04-26T00:44:30.960941Z","iopub.execute_input":"2024-04-26T00:44:30.961634Z","iopub.status.idle":"2024-04-26T00:44:31.208859Z","shell.execute_reply.started":"2024-04-26T00:44:30.961581Z","shell.execute_reply":"2024-04-26T00:44:31.207503Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#Let's load the molecule_smiles column of the test set, since it's smaller than train and this is just a demo.\ntest = pd.read_csv('/kaggle/input/leash-BELKA/test.csv',usecols=['molecule_smiles'],)","metadata":{"execution":{"iopub.status.busy":"2024-04-26T00:44:31.210632Z","iopub.execute_input":"2024-04-26T00:44:31.211142Z","iopub.status.idle":"2024-04-26T00:44:37.735740Z","shell.execute_reply.started":"2024-04-26T00:44:31.211100Z","shell.execute_reply":"2024-04-26T00:44:37.734383Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#Let's take a look at the first molecule in the dataframe\ntest_molecule = Chem.MolFromSmiles(test['molecule_smiles'].iloc[0])\ntest_molecule","metadata":{"execution":{"iopub.status.busy":"2024-04-26T00:44:37.737296Z","iopub.execute_input":"2024-04-26T00:44:37.737730Z","iopub.status.idle":"2024-04-26T00:44:37.769333Z","shell.execute_reply.started":"2024-04-26T00:44:37.737696Z","shell.execute_reply":"2024-04-26T00:44:37.767851Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Looks good! It's not exactly a real molecule though. That Dy is useful for helping us keep track of where the DNA is attached to this molecule, but it's likely going to get in our way if we want to convert this to a 3D structure that has some physical meaning.\n\nLet's replace the Dy with a methyl group (-CH3).","metadata":{}},{"cell_type":"code","source":"#We need to generate a molecule object that we can use to create our methyl group.\nMe = Chem.MolFromSmiles('C')\n\n#Now let's create a molecule object that represents the Dy atom. We can use this as a pattern we'd like to replace (below).\nDy = Chem.MolFromSmiles('[Dy]')\n\n#Let's use these groups to replace the Dy on our test molecule with CH3.\nnew_mol = AllChem.ReplaceSubstructs(test_molecule, Dy, Me)\nnew_mol","metadata":{"execution":{"iopub.status.busy":"2024-04-26T00:44:37.771260Z","iopub.execute_input":"2024-04-26T00:44:37.771754Z","iopub.status.idle":"2024-04-26T00:44:37.780218Z","shell.execute_reply.started":"2024-04-26T00:44:37.771718Z","shell.execute_reply":"2024-04-26T00:44:37.779214Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Note that the output of the above cell is a tuple with one entry. This is a helpful hint that this function could have returned more than one possibility. In this case, the replacement was pretty simple so we can just access the first element and call it a day.","metadata":{}},{"cell_type":"code","source":"#Access the first element of the tuple\nnew_mol = new_mol[0]\n#Sanitization is just a good habit after replacements/reactions as sometimes things get a little wonky.\nChem.SanitizeMol(new_mol)\n#Let's check out our new molecule\nnew_mol","metadata":{"execution":{"iopub.status.busy":"2024-04-26T00:44:37.781775Z","iopub.execute_input":"2024-04-26T00:44:37.782409Z","iopub.status.idle":"2024-04-26T00:44:37.801328Z","shell.execute_reply.started":"2024-04-26T00:44:37.782374Z","shell.execute_reply":"2024-04-26T00:44:37.799719Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Great, the Dy was replaced with a methyl group. A methyl group is a nice replacement if you want to do some 3D structure generations because it doesn't have many conformers (keeps things simple), but it might not do a good job of representing the linker used to attach DNA to our molecules. I haven't tested this hypothesis, but it's on my list of things to do. Regardless, knowing how to handle more complex transforms might be valuable to you in this competition.\n\nLet's take a look at something a bit more complicated and replace the Dy group with an ethylene glycol derived linker (a little closer to what was likely used as a linker on DNA).\n\nHere's a useful link for generating SMILES or SMARTS strings from things you can draw: https://pubchem.ncbi.nlm.nih.gov//edit3/index.html","metadata":{}},{"cell_type":"code","source":"#Let's create a linker molecule, just like we did above for the methyl replacement.\nlinker = Chem.MolFromSmiles('C(COC)OCC')\nlinker","metadata":{"execution":{"iopub.status.busy":"2024-04-26T00:44:37.805813Z","iopub.execute_input":"2024-04-26T00:44:37.806230Z","iopub.status.idle":"2024-04-26T00:44:37.824108Z","shell.execute_reply.started":"2024-04-26T00:44:37.806197Z","shell.execute_reply":"2024-04-26T00:44:37.822736Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#Let's try replacing our Dy with this group\nnew_mol = AllChem.ReplaceSubstructs(test_molecule, Dy, linker)\nnew_mol","metadata":{"execution":{"iopub.status.busy":"2024-04-26T00:44:37.825756Z","iopub.execute_input":"2024-04-26T00:44:37.826684Z","iopub.status.idle":"2024-04-26T00:44:37.839170Z","shell.execute_reply.started":"2024-04-26T00:44:37.826638Z","shell.execute_reply":"2024-04-26T00:44:37.837584Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"That's neat. It only made one new molecule, even though we didn't have a good way of telling it where to attach. Something doesn't seem right.","metadata":{}},{"cell_type":"code","source":"#Let's visualize our new molecule\nChem.SanitizeMol(new_mol[0])\nnew_mol[0]","metadata":{"execution":{"iopub.status.busy":"2024-04-26T00:44:37.840701Z","iopub.execute_input":"2024-04-26T00:44:37.841137Z","iopub.status.idle":"2024-04-26T00:44:37.862416Z","shell.execute_reply.started":"2024-04-26T00:44:37.841098Z","shell.execute_reply":"2024-04-26T00:44:37.860884Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Hmm...That's certainly a perfectly viable molecule someone could make, but it wasn't really what we were going for.\n\nTo fix this, we'll need to be a bit more specific. While the ReplaceSubstructs function is a very useful RDKit functionality, it works best when you're doing simple transformations. For more complex transformations, we'll need to learn a little about how to run virtual reactions.","metadata":{}},{"cell_type":"code","source":"#Let's define a reaction object using a SMARTS string.\nrxn = AllChem.ReactionFromSmarts('[*:1][Dy].[At][*:2]>>[*:1][*:2]')\nrxn","metadata":{"execution":{"iopub.status.busy":"2024-04-26T00:44:37.864358Z","iopub.execute_input":"2024-04-26T00:44:37.864880Z","iopub.status.idle":"2024-04-26T00:44:37.884109Z","shell.execute_reply.started":"2024-04-26T00:44:37.864837Z","shell.execute_reply":"2024-04-26T00:44:37.882629Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Let's break down that SMARTS string a little bit. \n\nThe >> delineates what's on the left and the right sides of the reaction.\n\nthe * means \"any atom\" and the ':number' maps the reactants to the products so the algorithm knows what is going where.\n\nSo, what we're saying here is essentially this: any atom bound to Dy will form a new bond to any atom bound to At.\n\nOf course, this isn't a real world reaction; it's just virtual. But that's what's great about it! Since it's not a reaction that exists in the real world, we can use it as a virtual reaction that allows us to do virtual transformations that are very specific!\n\nIf you want to learn more about SMARTS patterns, this is a great place to start: https://www.daylight.com/dayhtml/doc/theory/theory.smarts.html\n\nLet's use this reaction to attach our linker.","metadata":{}},{"cell_type":"code","source":"#First, we need to get a linker that has an At on it, signifying where we'd like it to react under our reaction scheme.\nnew_linker = Chem.MolFromSmiles('COCCOCC[At]')\nnew_linker","metadata":{"execution":{"iopub.status.busy":"2024-04-26T00:44:37.886251Z","iopub.execute_input":"2024-04-26T00:44:37.886734Z","iopub.status.idle":"2024-04-26T00:44:37.902069Z","shell.execute_reply.started":"2024-04-26T00:44:37.886699Z","shell.execute_reply":"2024-04-26T00:44:37.900982Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#Now let's use our reaction object to transform our test molecule with our linker.\n#Note: the way we defined the reaction, the molecule with Dy should be on the left and the molecule with At should be on the right.\nnew_mol = rxn.RunReactants((test_molecule, new_linker))\nnew_mol","metadata":{"execution":{"iopub.status.busy":"2024-04-26T00:44:37.904112Z","iopub.execute_input":"2024-04-26T00:44:37.904967Z","iopub.status.idle":"2024-04-26T00:44:37.916057Z","shell.execute_reply.started":"2024-04-26T00:44:37.904922Z","shell.execute_reply":"2024-04-26T00:44:37.914824Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We see now that our new_mol is a tuple of tuples. Having a tuple of tuples as an output is pretty useful for more complex reactions as there are often many places your reaction could theoretically have taken place, but we were really specific with our virtual imaginary reaction so we know there will only be one product here.","metadata":{}},{"cell_type":"code","source":"#Let's take a look at the molecule\nnew_mol = new_mol[0][0]\n#Good habit to sanitize after reactions\nChem.SanitizeMol(new_mol)\nnew_mol","metadata":{"execution":{"iopub.status.busy":"2024-04-26T00:44:37.918082Z","iopub.execute_input":"2024-04-26T00:44:37.918674Z","iopub.status.idle":"2024-04-26T00:44:37.939322Z","shell.execute_reply.started":"2024-04-26T00:44:37.918632Z","shell.execute_reply":"2024-04-26T00:44:37.937712Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Excellent! We can now replace our Dy atom with whatever we want!\n\nLet's write a little function that will apply this to a dataframe of our choice.","metadata":{}},{"cell_type":"code","source":"def run_reaction(SMILES, linker, reaction):\n    \"\"\"\n    This function takes a SMILES string, a linker, and a reaction object and runs a virtual reaction between the SMILES input molecule and the linker.\n    \n    SMILES(str): The SMILES string of the molecule you'd like to transform.\n    \n    linker(obj): The linker you'd like to attach to the molecule, as a RDKit mol object.\n    \n    reaction(obj): The reaction object that represents the virtual reaction.\n    ___________\n    returns(str): The SMILES string for the transformed molecule.\n    \n    \"\"\"\n    #Convert the SMILES to a molecule object\n    mol = Chem.MolFromSmiles(SMILES)\n    \n    #Run the reaction between the molecule object and the linker\n    products = reaction.RunReactants((mol,linker))\n    \n    #Get the first product of the reaction. If this was defined similarly to above, it should be the only product\n    product = products[0][0]\n    \n    #Good hygene\n    Chem.SanitizeMol(product)\n    \n    #Return the molecule object as a SMILES string\n    return Chem.MolToSmiles(product)","metadata":{"execution":{"iopub.status.busy":"2024-04-26T00:44:37.941777Z","iopub.execute_input":"2024-04-26T00:44:37.942294Z","iopub.status.idle":"2024-04-26T00:44:37.951770Z","shell.execute_reply.started":"2024-04-26T00:44:37.942251Z","shell.execute_reply":"2024-04-26T00:44:37.950008Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#Let's use mapply to transform the test set in parallel\n#Thanks to Hengck23 for the tip about mapply!\n#See: https://www.kaggle.com/competitions/leash-BELKA/discussion/492846\n\n#Initialize mapply to use all of our CPU cores and provide us a progress bar.\nmapply.init(\n    n_workers=-1,\n    progressbar=True,\n)","metadata":{"execution":{"iopub.status.busy":"2024-04-26T00:44:37.953445Z","iopub.execute_input":"2024-04-26T00:44:37.954342Z","iopub.status.idle":"2024-04-26T00:44:37.965161Z","shell.execute_reply.started":"2024-04-26T00:44:37.954307Z","shell.execute_reply":"2024-04-26T00:44:37.963763Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#Let's run our virtual reaction using the new_linker and rxn we previously wrote. This takes a little while...\ntest['molecule_smiles'] = test['molecule_smiles'].mapply(lambda x: run_reaction(x, new_linker, rxn))","metadata":{"execution":{"iopub.status.busy":"2024-04-26T00:44:37.966741Z","iopub.execute_input":"2024-04-26T00:44:37.967116Z","iopub.status.idle":"2024-04-26T00:56:47.257664Z","shell.execute_reply.started":"2024-04-26T00:44:37.967086Z","shell.execute_reply":"2024-04-26T00:56:47.255922Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#Let's take a look at some of the molecules we've created\nDraw.MolsToGridImage([Chem.MolFromSmiles(x) for x in test['molecule_smiles'].sample(12)], molsPerRow=4, subImgSize=(400,300))","metadata":{"execution":{"iopub.status.busy":"2024-04-26T00:56:47.260245Z","iopub.execute_input":"2024-04-26T00:56:47.260813Z","iopub.status.idle":"2024-04-26T00:56:47.546965Z","shell.execute_reply.started":"2024-04-26T00:56:47.260742Z","shell.execute_reply":"2024-04-26T00:56:47.545268Z"},"trusted":true},"execution_count":null,"outputs":[]}]}