{"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":30673,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# RDKit Descriptors\n## Introduction\n\nRDKit includes many descriptors that may be useful for predicting binding affinity, such as the number of a sepecific fragment and TPSA, just to name a few. \n\nLet's try to evaluate using the RDKit Descriptors and compare with ECFP4.\n\nThe code below is partly based on Andrew D. Blevins's notebook [Leash Tutorial - ECFPs and Random Forest](https://www.kaggle.com/code/andrewdblevins/leash-tutorial-ecfps-and-random-forest) \n","metadata":{}},{"cell_type":"code","source":"!pip install rdkit --quiet\n!pip install duckdb --quiet","metadata":{"execution":{"iopub.status.busy":"2024-04-06T22:20:44.299507Z","iopub.execute_input":"2024-04-06T22:20:44.299993Z","iopub.status.idle":"2024-04-06T22:21:11.738358Z","shell.execute_reply.started":"2024-04-06T22:20:44.299950Z","shell.execute_reply":"2024-04-06T22:21:11.735797Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Preparation","metadata":{}},{"cell_type":"code","source":"import duckdb\nimport pandas as pd\nfrom 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\nimport copy\nfrom rdkit.Chem import Descriptors\nfrom rdkit.ML.Descriptors.MoleculeDescriptors import MolecularDescriptorCalculator\n\nn_train = 300 # the number of train data to sample\ntest_size = 0.2\nrandom_state = 42\nclf = RandomForestClassifier(n_estimators=100, random_state=random_state)","metadata":{"execution":{"iopub.status.busy":"2024-04-06T22:21:11.741795Z","iopub.execute_input":"2024-04-06T22:21:11.742335Z","iopub.status.idle":"2024-04-06T22:21:15.504030Z","shell.execute_reply.started":"2024-04-06T22:21:11.742289Z","shell.execute_reply":"2024-04-06T22:21:15.502616Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_path = '/kaggle/input/leash-BELKA/train.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 {n_train})\n                        UNION ALL\n                        (SELECT *\n                        FROM parquet_scan('{train_path}')\n                        WHERE binds = 1\n                        ORDER BY random()\n                        LIMIT {n_train})\"\"\").df()\n\ncon.close()","metadata":{"execution":{"iopub.status.busy":"2024-04-06T22:21:15.505610Z","iopub.execute_input":"2024-04-06T22:21:15.506344Z","iopub.status.idle":"2024-04-06T22:21:54.318411Z","shell.execute_reply.started":"2024-04-06T22:21:15.506308Z","shell.execute_reply":"2024-04-06T22:21:54.316728Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Generate ECFPs\ndef generate_ecfp(molecule, radius=2, bits=1024):\n    if molecule is None:\n        return None\n    return list(AllChem.GetMorganFingerprintAsBitVect(molecule, radius, nBits=bits))\n\ndescriptor_list = [dname[0] for dname in Descriptors._descList]\ncalculator = MolecularDescriptorCalculator(descriptor_list)\n\n# Generate rdkit descriptor\ndef generate_rdkitdescriptor(molecule):\n    if molecule is None:\n        return None\n    return list(calculator.CalcDescriptors(molecule))","metadata":{"execution":{"iopub.status.busy":"2024-04-06T22:21:54.321240Z","iopub.execute_input":"2024-04-06T22:21:54.321667Z","iopub.status.idle":"2024-04-06T22:21:54.331338Z","shell.execute_reply.started":"2024-04-06T22:21:54.321627Z","shell.execute_reply":"2024-04-06T22:21:54.329684Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Let's create DataFrame with the RDKit Descriptors!","metadata":{}},{"cell_type":"code","source":"# Convert SMILES to RDKit molecules\ndf['molecule'] = df['molecule_smiles'].apply(Chem.MolFromSmiles)\n\ndef create_descriptor_df_for_test():\n    _df = pd.DataFrame(columns=['smiles'] + descriptor_list)\n\n    for smiles in df['molecule_smiles']:\n        mol = Chem.MolFromSmiles(smiles)\n        descriptors = generate_rdkitdescriptor(mol)\n        data = [smiles] + descriptors\n        _df.loc[len(_df)] = data\n    return _df\n\ndf_rdkitdescriptor = create_descriptor_df_for_test()\nprint(df_rdkitdescriptor.head(2))","metadata":{"execution":{"iopub.status.busy":"2024-04-06T22:21:54.333475Z","iopub.execute_input":"2024-04-06T22:21:54.333850Z","iopub.status.idle":"2024-04-06T22:22:19.449805Z","shell.execute_reply.started":"2024-04-06T22:21:54.333818Z","shell.execute_reply":"2024-04-06T22:22:19.448753Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Some columns include inf or nan. Let's remove at the moment.","metadata":{}},{"cell_type":"code","source":"import numpy as np\ncols_to_remove = df_rdkitdescriptor.columns[df_rdkitdescriptor.isin([np.nan, np.inf, -np.inf]).any(axis=0)].tolist()\nprint(f\"columns to remove: {cols_to_remove}\")\n# redefine descriptor calculator\ndescriptor_list = [dname[0] for dname in Descriptors._descList if dname[0] not in cols_to_remove]\ncalculator = MolecularDescriptorCalculator(descriptor_list)\n# Generate rdkit descriptor\ndef generate_rdkitdescriptor(molecule):\n    if molecule is None:\n        return None\n    return list(calculator.CalcDescriptors(molecule))","metadata":{"execution":{"iopub.status.busy":"2024-04-06T22:23:47.626317Z","iopub.execute_input":"2024-04-06T22:23:47.627438Z","iopub.status.idle":"2024-04-06T22:23:47.661474Z","shell.execute_reply.started":"2024-04-06T22:23:47.627382Z","shell.execute_reply":"2024-04-06T22:23:47.660432Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Train and Evaluate 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\ndef prepare_dataset(generate_feature):\n    # Combine ECFPs and one-hot encoded protein_name\n    _df = df.copy()\n    _df['feature'] = _df['molecule'].apply(generate_feature)\n    X = [feature + protein for feature, protein in zip(_df['feature'].tolist(), protein_onehot.tolist())]\n    y = _df['binds'].tolist()\n\n    # Split the data into train and test sets\n    X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=test_size, random_state=random_state)\n    return X_train, X_test, y_train, y_test\n\ndef fit_clf(clf, X_train, y_train):\n    # Create and train the model\n    clf_copy = copy.deepcopy(clf)\n    clf_copy.fit(X_train, y_train)\n    return clf_copy\n\ndef evaluate_map(generate_feature, clf):\n    X_train, X_test, y_train, y_test = prepare_dataset(generate_feature)\n    clf_fitted = fit_clf(clf, X_train, y_train)\n    \n    # Make predictions on the test set\n    y_pred_proba = clf_fitted.predict_proba(X_test)[:, 1]  # Probability of the positive class\n\n    # Calculate the mean average precision\n    map_score = average_precision_score(y_test, y_pred_proba)\n    return map_score, clf_fitted","metadata":{"execution":{"iopub.status.busy":"2024-04-06T22:24:01.260369Z","iopub.execute_input":"2024-04-06T22:24:01.261717Z","iopub.status.idle":"2024-04-06T22:24:01.271406Z","shell.execute_reply.started":"2024-04-06T22:24:01.261654Z","shell.execute_reply":"2024-04-06T22:24:01.270477Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nprint(\"Calculating Mean Average Precision (mAP) for each feature...\")\nmap_score_ecfp, _ = evaluate_map(generate_feature=generate_ecfp, clf=clf)\nprint(f\"ECFP: {map_score_ecfp:.2f}\")\nmap_score_rdkitdescriptor, clf_rdkitdescriptor = evaluate_map(generate_feature=generate_rdkitdescriptor, clf=clf)\nprint(f\"rdkit descriptor: {map_score_rdkitdescriptor:.2f}\")","metadata":{"execution":{"iopub.status.busy":"2024-04-06T22:24:02.065310Z","iopub.execute_input":"2024-04-06T22:24:02.066429Z","iopub.status.idle":"2024-04-06T22:24:14.132225Z","shell.execute_reply.started":"2024-04-06T22:24:02.066387Z","shell.execute_reply":"2024-04-06T22:24:14.131279Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"For this dataset, ECFPs performs better than the RDKit Descriptors.  \nFeature selection or further treatments may be necessary for further improvement.","metadata":{}}]}