{"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":"nvidiaTeslaT4","dataSources":[{"sourceId":67356,"databundleVersionId":8006601,"sourceType":"competition"},{"sourceId":170362385,"sourceType":"kernelVersion"}],"dockerImageVersionId":30674,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":true}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"!pip install rdkit","metadata":{"execution":{"iopub.status.busy":"2024-11-07T01:43:09.488652Z","iopub.execute_input":"2024-11-07T01:43:09.489079Z","iopub.status.idle":"2024-11-07T01:43:21.874596Z","shell.execute_reply.started":"2024-11-07T01:43:09.489048Z","shell.execute_reply":"2024-11-07T01:43:21.873313Z"},"trusted":true},"execution_count":null,"outputs":[]},{"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-07T01:43:21.877037Z","iopub.execute_input":"2024-11-07T01:43:21.877619Z","iopub.status.idle":"2024-11-07T01:43:34.273446Z","shell.execute_reply.started":"2024-11-07T01:43:21.877576Z","shell.execute_reply":"2024-11-07T01:43:34.272197Z"},"trusted":true},"execution_count":null,"outputs":[]},{"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-chemica|l-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-07T01:43:34.275244Z","iopub.execute_input":"2024-11-07T01:43:34.275613Z","iopub.status.idle":"2024-11-07T01:44:16.783712Z","shell.execute_reply.started":"2024-11-07T01:43:34.275579Z","shell.execute_reply":"2024-11-07T01:44:16.782797Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df.head()","metadata":{"execution":{"iopub.status.busy":"2024-11-07T01:44:16.788217Z","iopub.execute_input":"2024-11-07T01:44:16.788525Z","iopub.status.idle":"2024-11-07T01:44:16.803262Z","shell.execute_reply.started":"2024-11-07T01:44:16.788499Z","shell.execute_reply":"2024-11-07T01:44:16.802255Z"},"trusted":true},"execution_count":null,"outputs":[]},{"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 sklearn.ensemble import RandomForestClassifier\nfrom sklearn.model_selection import train_test_split\nfrom sklearn.metrics import average_precision_score\nfrom sklearn.preprocessing import OneHotEncoder\nfrom rdkit.Chem import Descriptors\nfrom imblearn.ensemble import BalancedRandomForestClassifier","metadata":{"execution":{"iopub.status.busy":"2024-11-07T01:44:16.804547Z","iopub.execute_input":"2024-11-07T01:44:16.804848Z","iopub.status.idle":"2024-11-07T01:44:16.955839Z","shell.execute_reply.started":"2024-11-07T01:44:16.804821Z","shell.execute_reply":"2024-11-07T01:44:16.955071Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Convert SMILES to RDKit molecules\ndf['molecule'] = df['molecule_smiles'].apply(Chem.MolFromSmiles)\n\n# Define a function to calculate molecular descriptors\ndef calculate_descriptors(molecule):\n    if molecule is None:\n        return [None] * 5  # Adjust the number of descriptors if you add/remove\n    return [\n        Descriptors.MolWt(molecule),                # Molecular weight\n        Descriptors.MolLogP(molecule),              # Octanol-water partition coefficient (logP)\n        Descriptors.NumHAcceptors(molecule),        # Number of hydrogen bond acceptors\n        Descriptors.NumHDonors(molecule),           # Number of hydrogen bond donors\n        Descriptors.TPSA(molecule)                  # Topological polar surface area\n    ]\n\n# Apply descriptor calculation\ndescriptor_columns = ['MolWt', 'MolLogP', 'NumHAcceptors', 'NumHDonors', 'TPSA']\ndf[descriptor_columns] = df['molecule'].apply(calculate_descriptors).apply(pd.Series)","metadata":{"execution":{"iopub.status.busy":"2024-11-07T01:44:16.956965Z","iopub.execute_input":"2024-11-07T01:44:16.957745Z","iopub.status.idle":"2024-11-07T01:47:02.412634Z","shell.execute_reply.started":"2024-11-07T01:44:16.957716Z","shell.execute_reply":"2024-11-07T01:47:02.411793Z"},"trusted":true},"execution_count":null,"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 descriptors and one-hot encoded protein_name\nX = [list(descriptors) + list(protein) for descriptors, protein in zip(df[descriptor_columns].values, protein_onehot)]\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 = RandomForestClassifier(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}\")\n","metadata":{"execution":{"iopub.status.busy":"2024-11-07T01:47:02.414873Z","iopub.execute_input":"2024-11-07T01:47:02.415225Z","iopub.status.idle":"2024-11-07T01:47:22.766245Z","shell.execute_reply.started":"2024-11-07T01:47:02.415187Z","shell.execute_reply":"2024-11-07T01:47:22.765310Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Encoder with ECFPs","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-chemica|l-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-07T01:47:22.767577Z","iopub.execute_input":"2024-11-07T01:47:22.767882Z","iopub.status.idle":"2024-11-07T01:48:10.741400Z","shell.execute_reply.started":"2024-11-07T01:47:22.767853Z","shell.execute_reply":"2024-11-07T01:48:10.740496Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# 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)\n","metadata":{"execution":{"iopub.status.busy":"2024-11-07T01:48:10.743105Z","iopub.execute_input":"2024-11-07T01:48:10.743866Z","iopub.status.idle":"2024-11-07T01:52:00.734960Z","shell.execute_reply.started":"2024-11-07T01:48:10.743828Z","shell.execute_reply":"2024-11-07T01:52:00.734097Z"},"trusted":true},"execution_count":null,"outputs":[]},{"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 = RandomForestClassifier(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-07T01:52:00.737427Z","iopub.execute_input":"2024-11-07T01:52:00.737730Z","iopub.status.idle":"2024-11-07T01:53:45.011646Z","shell.execute_reply.started":"2024-11-07T01:52:00.737703Z","shell.execute_reply":"2024-11-07T01:53:45.010432Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Using BalancedRandomForestClassifier","metadata":{}},{"cell_type":"code","source":"# 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-07T01:53:45.013166Z","iopub.execute_input":"2024-11-07T01:53:45.013588Z","iopub.status.idle":"2024-11-07T01:54:40.612607Z","shell.execute_reply.started":"2024-11-07T01:53:45.013549Z","shell.execute_reply":"2024-11-07T01:54:40.611540Z"},"trusted":true},"execution_count":null,"outputs":[]}]}