{"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":"gpu","dataSources":[{"sourceId":67356,"databundleVersionId":8006601,"sourceType":"competition"}],"dockerImageVersionId":30673,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":true}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# Leash Tutorial - ECFPs and Random Forest\n## Introduction\n\nThere are many ways to represent molecules for machine learning. \n\nIn this tutorial we will go through one of the simplest: ECFPs [[1]](https://pubs.acs.org/doi/10.1021/ci100050t) and Random Forest. This technique is surprisingly powerful, and on previous benchmarks often gets uncomfortably close to the state of the art.\n\nFirst molecule graphs are broken into bags of subgraphs of varying sizes.\n\n![ecfp featurizing process (chemaxon)](https://docs.chemaxon.com/display/docs/images/download/attachments/1806333/ecfp_generation.png)\n\nThen the bag of subgraphs is hashed into a bit vector\n\n![hashing process (chemaxon)](https://docs.chemaxon.com/display/docs/images/download/attachments/1806333/ecfp_folding.png)\n\nThis can be thought of as analogous to the [hashing trick](https://en.wikipedia.org/wiki/Feature_hashing) [[2]](https://alex.smola.org/papers/2009/Weinbergeretal09.pdf) on bag of words for NLP problems, from the days before transformers. \n\nRDKit, an open-source cheminformatics tool, is used for generating ECFP features. It facilitates the creation of hashed bit vectors, streamlining the process. We can install it as follows:","metadata":{}},{"cell_type":"code","source":"!pip install rdkit","metadata":{"execution":{"iopub.status.busy":"2024-11-06T16:42:00.657744Z","iopub.execute_input":"2024-11-06T16:42:00.658118Z","iopub.status.idle":"2024-11-06T16:42:18.171578Z","shell.execute_reply.started":"2024-11-06T16:42:00.658086Z","shell.execute_reply":"2024-11-06T16:42:18.170437Z"},"trusted":true},"execution_count":1,"outputs":[{"name":"stdout","text":"Collecting rdkit\n Downloading rdkit-2024.3.5-cp310-cp310-manylinux_2_28_x86_64.whl.metadata (3.9 kB)\nRequirement already satisfied: numpy in /opt/conda/lib/python3.10/site-packages (from rdkit) (1.26.4)\nRequirement already satisfied: Pillow in /opt/conda/lib/python3.10/site-packages (from rdkit) (9.5.0)\nDownloading rdkit-2024.3.5-cp310-cp310-manylinux_2_28_x86_64.whl (33.1 MB)\n\u001b[2K \u001b[90m━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━\u001b[0m \u001b[32m33.1/33.1 MB\u001b[0m \u001b[31m46.6 MB/s\u001b[0m eta \u001b[36m0:00:00\u001b[0m:00:01\u001b[0m00:01\u001b[0m\n\u001b[?25hInstalling collected packages: rdkit\nSuccessfully installed rdkit-2024.3.5\n","output_type":"stream"}]},{"cell_type":"markdown","source":"The training set is pretty big, but we can treat the parquet files as databases using duckdb. We will use this to sample down to a smaller dataset for demonstration purposes. Lets install duckdb as well.","metadata":{}},{"cell_type":"code","source":"!pip install duckdb","metadata":{"execution":{"iopub.status.busy":"2024-11-06T16:42:18.173819Z","iopub.execute_input":"2024-11-06T16:42:18.174246Z","iopub.status.idle":"2024-11-06T16:42:33.370803Z","shell.execute_reply.started":"2024-11-06T16:42:18.174204Z","shell.execute_reply":"2024-11-06T16:42:33.369693Z"},"trusted":true},"execution_count":2,"outputs":[{"name":"stdout","text":"Collecting duckdb\n Downloading duckdb-1.1.3-cp310-cp310-manylinux_2_17_x86_64.manylinux2014_x86_64.whl.metadata (762 bytes)\nDownloading duckdb-1.1.3-cp310-cp310-manylinux_2_17_x86_64.manylinux2014_x86_64.whl (20.1 MB)\n\u001b[2K \u001b[90m━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━\u001b[0m \u001b[32m20.1/20.1 MB\u001b[0m \u001b[31m68.7 MB/s\u001b[0m eta \u001b[36m0:00:00\u001b[0m:00:01\u001b[0m00:01\u001b[0m\n\u001b[?25hInstalling collected packages: duckdb\nSuccessfully installed duckdb-1.1.3\n","output_type":"stream"}]},{"cell_type":"markdown","source":"## Data Preparation\n\nThe training and testing data paths are defined for the .parquet files. We use duckdb to scan search through the large training sets. Just to get started lets sample out an equal number of positive and negatives. \n\nThis query selects an equal number of samples where binds equals 0 (non-binding) and 1 (binding), limited to 30,000 each, to avoid model bias towards a particular class.","metadata":{}},{"cell_type":"code","source":"import duckdb\nimport pandas as pd\n\ntrain_path = '/kaggle/input/leash-predict-chemical-bindings/train.parquet'\ntest_path = '/kaggle/input/leash-predict-chemical-bindings/test.parquet'\n\ncon = duckdb.connect()\n\ndf = con.query(f\"\"\"(SELECT *\n FROM parquet_scan('{train_path}')\n WHERE binds = 0\n ORDER BY random()\n LIMIT 100000)\n UNION ALL\n (SELECT *\n FROM parquet_scan('{train_path}')\n WHERE binds = 1\n ORDER BY random()\n LIMIT 30000)\"\"\").df()\n\ncon.close()","metadata":{"execution":{"iopub.status.busy":"2024-11-06T16:42:33.372218Z","iopub.execute_input":"2024-11-06T16:42:33.372567Z","iopub.status.idle":"2024-11-06T16:43:28.664289Z","shell.execute_reply.started":"2024-11-06T16:42:33.372534Z","shell.execute_reply":"2024-11-06T16:43:28.663264Z"},"trusted":true},"execution_count":3,"outputs":[{"output_type":"display_data","data":{"text/plain":"FloatProgress(value=0.0, layout=Layout(width='auto'), style=ProgressStyle(bar_color='black'))","application/vnd.jupyter.widget-view+json":{"version_major":2,"version_minor":0,"model_id":"f1d59ac23298471eb271137bd617a6f9"}},"metadata":{}}]},{"cell_type":"code","source":"df.head()","metadata":{"execution":{"iopub.status.busy":"2024-11-06T16:43:28.666718Z","iopub.execute_input":"2024-11-06T16:43:28.667211Z","iopub.status.idle":"2024-11-06T16:43:28.68723Z","shell.execute_reply.started":"2024-11-06T16:43:28.667179Z","shell.execute_reply":"2024-11-06T16:43:28.686023Z"},"trusted":true},"execution_count":4,"outputs":[{"execution_count":4,"output_type":"execute_result","data":{"text/plain":" id buildingblock1_smiles \\\n0 48709474 Cc1cc(Br)c(NC(=O)OCC2c3ccccc3-c3ccccc32)c(C(=O... \n1 191685433 O=C(Nc1ccc(C(=O)O)cc1F)OCC1c2ccccc2-c2ccccc21 \n2 172433919 O=C(Nc1cc(F)c(Br)cc1C(=O)O)OCC1c2ccccc2-c2ccccc21 \n3 62222940 Cn1cc(C[C@@H](NC(=O)OCC2c3ccccc3-c3ccccc32)C(=... \n4 7355871 C=CCC(NC(=O)OCC1c2ccccc2-c2ccccc21)C(=O)O \n\n buildingblock2_smiles buildingblock3_smiles \\\n0 NCCCCN1CCCC1 Nc1ccncc1[N+](=O)[O-] \n1 Nc1cnccn1 COc1cc2nc(Cl)nc(N)c2cc1OC \n2 COc1nc(Br)ccc1N Nc1nc[nH]n1 \n3 CCON(CC)C(=O)CN.Cl Cc1nc(CN)oc1C.Cl.Cl \n4 NCc1ccccc1CS(=O)(=O)N1CCOCC1 NCc1ccc(N2CCOCC2)cc1C(F)(F)F \n\n molecule_smiles protein_name binds \n0 Cc1cc(Br)c(Nc2nc(NCCCCN3CCCC3)nc(Nc3ccncc3[N+]... HSA 0 \n1 COc1cc2nc(Cl)nc(Nc3nc(Nc4cnccn4)nc(Nc4ccc(C(=O... HSA 0 \n2 COc1nc(Br)ccc1Nc1nc(Nc2nc[nH]n2)nc(Nc2cc(F)c(B... BRD4 0 \n3 CCON(CC)C(=O)CNc1nc(NCc2nc(C)c(C)o2)nc(N[C@H](... BRD4 0 \n4 C=CCC(Nc1nc(NCc2ccccc2CS(=O)(=O)N2CCOCC2)nc(NC... BRD4 0 ","text/html":"
\n\n\n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n
idbuildingblock1_smilesbuildingblock2_smilesbuildingblock3_smilesmolecule_smilesprotein_namebinds
048709474Cc1cc(Br)c(NC(=O)OCC2c3ccccc3-c3ccccc32)c(C(=O...NCCCCN1CCCC1Nc1ccncc1[N+](=O)[O-]Cc1cc(Br)c(Nc2nc(NCCCCN3CCCC3)nc(Nc3ccncc3[N+]...HSA0
1191685433O=C(Nc1ccc(C(=O)O)cc1F)OCC1c2ccccc2-c2ccccc21Nc1cnccn1COc1cc2nc(Cl)nc(N)c2cc1OCCOc1cc2nc(Cl)nc(Nc3nc(Nc4cnccn4)nc(Nc4ccc(C(=O...HSA0
2172433919O=C(Nc1cc(F)c(Br)cc1C(=O)O)OCC1c2ccccc2-c2ccccc21COc1nc(Br)ccc1NNc1nc[nH]n1COc1nc(Br)ccc1Nc1nc(Nc2nc[nH]n2)nc(Nc2cc(F)c(B...BRD40
362222940Cn1cc(C[C@@H](NC(=O)OCC2c3ccccc3-c3ccccc32)C(=...CCON(CC)C(=O)CN.ClCc1nc(CN)oc1C.Cl.ClCCON(CC)C(=O)CNc1nc(NCc2nc(C)c(C)o2)nc(N[C@H](...BRD40
47355871C=CCC(NC(=O)OCC1c2ccccc2-c2ccccc21)C(=O)ONCc1ccccc1CS(=O)(=O)N1CCOCC1NCc1ccc(N2CCOCC2)cc1C(F)(F)FC=CCC(Nc1nc(NCc2ccccc2CS(=O)(=O)N2CCOCC2)nc(NC...BRD40
\n
"},"metadata":{}}]},{"cell_type":"markdown","source":"## Feature Preprocessing\n\nLets grab the smiles for the fully assembled molecule `molecule_smiles` and generate ecfps for it. We could choose different radiuses or bits, but 2 and 1024 is pretty standard.","metadata":{}},{"cell_type":"code","source":"from rdkit import Chem\nfrom rdkit.Chem import AllChem\nfrom imblearn.ensemble import BalancedRandomForestClassifier\nfrom sklearn.model_selection import train_test_split\nfrom sklearn.metrics import average_precision_score\nfrom sklearn.preprocessing import OneHotEncoder\n\n# Convert SMILES to RDKit molecules\ndf['molecule'] = df['molecule_smiles'].apply(Chem.MolFromSmiles)\n\n# Initialize the Morgan fingerprint generator\nmorgan_generator = AllChem.GetMorganGenerator(radius=2, fpSize=1024)\n\n# Generate ECFPs using MorganGenerator\ndef generate_ecfp(molecule):\n if molecule is None:\n return None\n return list(morgan_generator.GetFingerprint(molecule))\n\ndf['ecfp'] = df['molecule'].apply(generate_ecfp)","metadata":{"execution":{"iopub.status.busy":"2024-11-06T16:45:27.298052Z","iopub.execute_input":"2024-11-06T16:45:27.299107Z","iopub.status.idle":"2024-11-06T16:49:52.859631Z","shell.execute_reply.started":"2024-11-06T16:45:27.299064Z","shell.execute_reply":"2024-11-06T16:49:52.858492Z"},"trusted":true},"execution_count":5,"outputs":[]},{"cell_type":"markdown","source":"## Train Model","metadata":{}},{"cell_type":"code","source":"# One-hot encode the protein_name\nonehot_encoder = OneHotEncoder(sparse_output=False)\nprotein_onehot = onehot_encoder.fit_transform(df['protein_name'].values.reshape(-1, 1))\n\n# Combine ECFPs and one-hot encoded protein_name\nX = [ecfp + protein for ecfp, protein in zip(df['ecfp'].tolist(), protein_onehot.tolist())]\ny = df['binds'].tolist()\n\n# Split the data into train and test sets\nX_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)\n\n# Create and train the random forest model\nrf_model = BalancedRandomForestClassifier(n_estimators=100, random_state=42)\nrf_model.fit(X_train, y_train)\n\n# Make predictions on the test set\ny_pred_proba = rf_model.predict_proba(X_test)[:, 1] # Probability of the positive class\n\n# Calculate the mean average precision\nmap_score = average_precision_score(y_test, y_pred_proba)\nprint(f\"Mean Average Precision (mAP): {map_score:.2f}\")","metadata":{"execution":{"iopub.status.busy":"2024-11-06T16:50:20.144155Z","iopub.execute_input":"2024-11-06T16:50:20.14484Z","iopub.status.idle":"2024-11-06T16:51:18.372881Z","shell.execute_reply.started":"2024-11-06T16:50:20.144798Z","shell.execute_reply":"2024-11-06T16:51:18.371734Z"},"trusted":true},"execution_count":6,"outputs":[{"name":"stderr","text":"/opt/conda/lib/python3.10/site-packages/imblearn/ensemble/_forest.py:576: FutureWarning: The default of `sampling_strategy` will change from `'auto'` to `'all'` in version 0.13. This change will follow the implementation proposed in the original paper. Set to `'all'` to silence this warning and adopt the future behaviour.\n warn(\n/opt/conda/lib/python3.10/site-packages/imblearn/ensemble/_forest.py:588: FutureWarning: The default of `replacement` will change from `False` to `True` in version 0.13. This change will follow the implementation proposed in the original paper. Set to `True` to silence this warning and adopt the future behaviour.\n warn(\n/opt/conda/lib/python3.10/site-packages/imblearn/ensemble/_forest.py:600: FutureWarning: The default of `bootstrap` will change from `True` to `False` in version 0.13. This change will follow the implementation proposed in the original paper. Set to `False` to silence this warning and adopt the future behaviour.\n warn(\n","output_type":"stream"},{"name":"stdout","text":"Mean Average Precision (mAP): 0.90\n","output_type":"stream"}]},{"cell_type":"markdown","source":"Look at that Average Precision score. We did amazing! \n\nActually no, we just overfit. This is likely recurring theme for this competition. It is easy to predict molecules that come from the same corner of chemical space, but generalizing to new molecules is extremely difficult.","metadata":{}},{"cell_type":"markdown","source":"## Test Prediction\n\n The trained Random Forest model is then used to predict the binding probabilities. These predictions are saved to a CSV file, which serves as the submission file for the Kaggle competition.","metadata":{}},{"cell_type":"code","source":"import os\n\n# Process the test.parquet file chunk by chunk\ntest_file = '/kaggle/input/leash-predict-chemical-bindings/test.csv'\noutput_file = 'submission.csv' # Specify the path and filename for the output file\n\n# Read the test.parquet file into a pandas DataFrame\nfor df_test in pd.read_csv(test_file, chunksize=100000):\n\n # Generate ECFPs for the molecule_smiles\n df_test['molecule'] = df_test['molecule_smiles'].apply(Chem.MolFromSmiles)\n df_test['ecfp'] = df_test['molecule'].apply(generate_ecfp)\n\n # One-hot encode the protein_name\n protein_onehot = onehot_encoder.transform(df_test['protein_name'].values.reshape(-1, 1))\n\n # Combine ECFPs and one-hot encoded protein_name\n X_test = [ecfp + protein for ecfp, protein in zip(df_test['ecfp'].tolist(), protein_onehot.tolist())]\n\n # Predict the probabilities\n probabilities = rf_model.predict_proba(X_test)[:, 1]\n\n # Create a DataFrame with 'id' and 'probability' columns\n output_df = pd.DataFrame({'id': df_test['id'], 'binds': probabilities})\n\n # Save the output DataFrame to a CSV file\n output_df.to_csv(output_file, index=False, mode='a', header=not os.path.exists(output_file))\n","metadata":{"execution":{"iopub.status.busy":"2024-11-06T16:52:05.400279Z","iopub.execute_input":"2024-11-06T16:52:05.400716Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}