{"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":"# Scikit-fingerprints\nHi! I've joined this competition about a week ago. So far I've met a few challenges (I guess similarly to all other competitors) with the task we've been given and one of them is the sheer amount of data. Many people here try to use molecular fingerprints for obtaining features from the molecules. Of course, RDKit is currently the industry standard and offers some algorithms to do that. The problem is, that the notebooks I've seen so far, all use the sequential processing, which takes a lot of time. Also, the RDKit documentation and interface is not the most convenient to use :)\n\nLuckily, not so long ago me and a few colleagues of mine developed a library for this exact purpose: **scikit-fingerprints** (https://github.com/scikit-fingerprints/scikit-fingerprints). It has already been released and is being actively extended. \n\n### Main benefits of this library are:\n- parallel computation\n- 20+ implemented fingerprint algorithms\n- scikit-learn-like interface, which allows you to plug it directly into your pipeline\n- clear documentation (https://scikit-fingerprints.github.io/scikit-fingerprints/)\n- library is pip-installable\n\nThis short notebook is created to demonstrate how you can utilize scikit-fingerprints library to process your data using molecular fingerprints in this competition. As an extra, we'll plug the results into a *RandomForestClassifier* for validation, so you'll get the complete workflow.\n\n### Conclusions from this notebook (TL;DR):\n- Scikit-fingerprints runtime was **3min 35s**, whereas sequential calculations with RDKit took **30min 24s** - that's about **x10 times** faster!\n- Calculation of fingerprints can be easily done, just by passing list of SMILES strings into the transformer (2 lines of code)\n- I needed to lower the *DATA_SIZE*, because sequential computing appeared to be even more memory consuming than parallel. However, 2 mil rows are enough to show the speed difference, and that was the point.\n\nThere are lots of different fingerprint algorithms for you to try in this library, but keep in mind the memory limitations in this challenge! Some of them require extensive calculations, so it may not be possible for you to use them in this particular competition.","metadata":{}},{"cell_type":"code","source":"!pip install -qq duckdb rdkit scikit-fingerprints","metadata":{"execution":{"iopub.status.busy":"2024-05-01T17:50:51.005005Z","iopub.execute_input":"2024-05-01T17:50:51.005397Z","iopub.status.idle":"2024-05-01T17:51:13.406705Z","shell.execute_reply.started":"2024-05-01T17:50:51.005366Z","shell.execute_reply":"2024-05-01T17:51:13.404937Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import gc\nimport duckdb\nimport numpy as np\nimport pandas as pd\n\nfrom rdkit import Chem\nfrom rdkit.Chem import AllChem\nfrom skfp.fingerprints import AtomPairFingerprint, ECFPFingerprint, MACCSFingerprint\nfrom sklearn.ensemble import RandomForestClassifier\nfrom sklearn.metrics import average_precision_score\nfrom sklearn.model_selection import train_test_split\nfrom tqdm import tqdm","metadata":{"execution":{"iopub.status.busy":"2024-05-01T17:51:13.411334Z","iopub.execute_input":"2024-05-01T17:51:13.411859Z","iopub.status.idle":"2024-05-01T17:51:29.862561Z","shell.execute_reply.started":"2024-05-01T17:51:13.411808Z","shell.execute_reply":"2024-05-01T17:51:29.861216Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"DATASET_DIR = '/kaggle/input/leash-BELKA/'\ntrain_path = DATASET_DIR + 'train.parquet'\nDATA_SIZE = 2_000_000","metadata":{"execution":{"iopub.status.busy":"2024-05-01T17:51:29.864091Z","iopub.execute_input":"2024-05-01T17:51:29.864744Z","iopub.status.idle":"2024-05-01T17:51:29.870689Z","shell.execute_reply.started":"2024-05-01T17:51:29.864702Z","shell.execute_reply":"2024-05-01T17:51:29.869470Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\n\n# As you've probably already seen this strategy in many other notebooks, we'll use duckdb.\n# Lets fetch some random data directly from the provided parquet file:\n\ncon = duckdb.connect()\n\ndf_train = con.query(f\"\"\"(SELECT *\n                        FROM parquet_scan('{train_path}')\n                        ORDER BY random()\n                        LIMIT {DATA_SIZE})\n                        \"\"\").df()\ncon.close()","metadata":{"execution":{"iopub.status.busy":"2024-05-01T17:52:29.326176Z","iopub.execute_input":"2024-05-01T17:52:29.326698Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_train.head()","metadata":{"execution":{"iopub.status.busy":"2024-05-01T17:51:29.893258Z","iopub.status.idle":"2024-05-01T17:51:29.893731Z","shell.execute_reply.started":"2024-05-01T17:51:29.893489Z","shell.execute_reply":"2024-05-01T17:51:29.893506Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Secondly, let's get the molecules in SMILES format to our X variable, and protein bindings into y. \n# Also, we should divide this data into 3 targets, based on the *protein_name*.\n\n# ------------------------------------------------------------------------\n#                                   Split\n# ------------------------------------------------------------------------\n\ntrain_data, valid_data = train_test_split(df_train, test_size=0.2, random_state=0, shuffle=True)\n\n# ------------------------------------------------------------------------\n#                                 Training\n# ------------------------------------------------------------------------\n\ndata_seh = train_data.loc[train_data[\"protein_name\"] == \"sEH\"]\ndata_brd4 = train_data.loc[train_data[\"protein_name\"] == \"BRD4\"]\ndata_hsa = train_data.loc[train_data[\"protein_name\"] == \"HSA\"]\n\nX_seh = data_seh[\"molecule_smiles\"]\nX_brd4 = data_brd4[\"molecule_smiles\"]\nX_hsa = data_hsa[\"molecule_smiles\"]\n\ny_seh = data_seh[\"binds\"]\ny_brd4 = data_brd4[\"binds\"]\ny_hsa = data_hsa[\"binds\"]\n\n# ------------------------------------------------------------------------\n#                                 Validation\n# ------------------------------------------------------------------------\n\ndata_valid_seh = valid_data.loc[valid_data[\"protein_name\"] == \"sEH\"]\ndata_valid_brd4 = valid_data.loc[valid_data[\"protein_name\"] == \"BRD4\"]\ndata_valid_hsa = valid_data.loc[valid_data[\"protein_name\"] == \"HSA\"]\n\nX_valid_seh = data_valid_seh[\"molecule_smiles\"]\nX_valid_brd4 = data_valid_brd4[\"molecule_smiles\"]\nX_valid_hsa = data_valid_hsa[\"molecule_smiles\"]\n\ny_valid_seh = data_valid_seh[\"binds\"]\ny_valid_brd4 = data_valid_brd4[\"binds\"]\ny_valid_hsa = data_valid_hsa[\"binds\"]","metadata":{"execution":{"iopub.status.busy":"2024-05-01T17:51:29.895425Z","iopub.status.idle":"2024-05-01T17:51:29.895960Z","shell.execute_reply.started":"2024-05-01T17:51:29.895719Z","shell.execute_reply":"2024-05-01T17:51:29.895739Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Now let's create instances of a few fingerprint classes. In our example we'll only use ECFP, but I encourage you to try also other ones.\n# Note, that you should pass the n_jobs = -1 argument, so that the computation will run in parallel.\n\nfp_ap_transformer = AtomPairFingerprint(fp_size=2048, n_jobs=-1)\nfp_ecfp_transformer = ECFPFingerprint(fp_size=2048, radius=3, n_jobs=-1)\nfp_maccs_transformer = MACCSFingerprint(n_jobs=-1)","metadata":{"execution":{"iopub.status.busy":"2024-05-01T17:51:29.897933Z","iopub.status.idle":"2024-05-01T17:51:29.898761Z","shell.execute_reply.started":"2024-05-01T17:51:29.898484Z","shell.execute_reply":"2024-05-01T17:51:29.898505Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The big problem in this competition is the limited RAM memory. Parallel computation is faster, but at the moment of making calculations it requires more memory than sequential solutions. Because of this fact, here is a small helper function: it divides the data into a smaller batches, so that the memory won't go over the limit. Feel free to modify it by adding other fingerprints.","metadata":{}},{"cell_type":"code","source":"ELEMENTS_PER_WORKER = 100_000\n\ndef calculate_fp(x, fp_name: str):\n    if fp_name == 'ecfp':\n        fp_transformer = fp_ecfp_transformer\n    elif fp_name == 'ap':\n        fp_transformer = fp_ap_transformer\n    elif fp_name == 'maccs':\n        fp_transformer = fp_maccs_transformer\n    else:\n        raise Exception('Wrong fingerprint name!')\n    \n    middle_parts = []\n    k_splits = x.shape[0] // ELEMENTS_PER_WORKER\n        \n    for i in tqdm(range(k_splits)):\n        middle_parts.append(fp_transformer.transform(x[i * ELEMENTS_PER_WORKER: (i + 1) * ELEMENTS_PER_WORKER]))\n        \n    if x.shape[0] % ELEMENTS_PER_WORKER > 0:   \n        middle_parts.append(fp_transformer.transform(x[k_splits * ELEMENTS_PER_WORKER:]))\n    \n    return np.concatenate(middle_parts)","metadata":{"execution":{"iopub.status.busy":"2024-05-01T17:51:29.900794Z","iopub.status.idle":"2024-05-01T17:51:29.901191Z","shell.execute_reply.started":"2024-05-01T17:51:29.901004Z","shell.execute_reply":"2024-05-01T17:51:29.901021Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Fingerprint computation times:\nNow the main part - let's firstly calculate ECFP using scikit-fingerprints and later using RDKit.\n\nIf you don't have the patience to wait for the cells below to complete, that's okay :)\n\nScikit-fingerprints completed in 3min 35s and sequential computation with RDKit in 30min 24s - that's a **x10 times** speed-up!","metadata":{}},{"cell_type":"code","source":"%%time\n# ------------------------------------------------------------------------\n#                               Scikit-fingerprints\n# ------------------------------------------------------------------------\n\nX_seh_transformed = calculate_fp(X_seh, fp_name='ecfp')","metadata":{"execution":{"iopub.status.busy":"2024-05-01T17:51:29.903391Z","iopub.status.idle":"2024-05-01T17:51:29.904003Z","shell.execute_reply.started":"2024-05-01T17:51:29.903716Z","shell.execute_reply":"2024-05-01T17:51:29.903742Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Delete the result, to free up memory for the next step\ndel X_seh_transformed\ngc.collect()","metadata":{"execution":{"iopub.status.busy":"2024-05-01T17:51:29.906653Z","iopub.status.idle":"2024-05-01T17:51:29.907123Z","shell.execute_reply.started":"2024-05-01T17:51:29.906905Z","shell.execute_reply":"2024-05-01T17:51:29.906923Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\n# ------------------------------------------------------------------------\n#                                 Sequential RDKit\n# ------------------------------------------------------------------------\n# (Remember to transform SMILES into Mol first)\n\nX_seh_transformed_sequentially = X_seh.apply(\n    lambda molecule: list(AllChem.GetMorganFingerprintAsBitVect(Chem.MolFromSmiles(molecule), radius=3, nBits=2048))\n)","metadata":{"execution":{"iopub.status.busy":"2024-05-01T17:51:29.908842Z","iopub.status.idle":"2024-05-01T17:51:29.909682Z","shell.execute_reply.started":"2024-05-01T17:51:29.909434Z","shell.execute_reply":"2024-05-01T17:51:29.909454Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"del X_seh_transformed_sequentially\ngc.collect()","metadata":{"execution":{"iopub.status.busy":"2024-05-01T17:51:29.911958Z","iopub.status.idle":"2024-05-01T17:51:29.912368Z","shell.execute_reply.started":"2024-05-01T17:51:29.912181Z","shell.execute_reply":"2024-05-01T17:51:29.912198Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Full preprocessing and training using RandomForestClassifier:\nTo finish this notebook, let's process the data once more - this time for each target - and at last fit it into a *RandomForestClassifier*.","metadata":{}},{"cell_type":"code","source":"X_seh_transformed = calculate_fp(X_seh, fp_name='ecfp')\nX_brd4_transformed = calculate_fp(X_brd4, fp_name='ecfp')\nX_hsa_transformed = calculate_fp(X_hsa, fp_name='ecfp')\n\nX_valid_seh_transformed = calculate_fp(X_valid_seh, fp_name='ecfp')\nX_valid_brd4_transformed = calculate_fp(X_valid_brd4, fp_name='ecfp')\nX_valid_hsa_transformed = calculate_fp(X_valid_hsa, fp_name='ecfp')","metadata":{"execution":{"iopub.status.busy":"2024-05-01T17:51:29.913863Z","iopub.status.idle":"2024-05-01T17:51:29.914737Z","shell.execute_reply.started":"2024-05-01T17:51:29.914488Z","shell.execute_reply":"2024-05-01T17:51:29.914506Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"clf_seh = RandomForestClassifier(\n    class_weight='balanced',\n    n_jobs=-1,\n    random_state=0\n)\nclf_brd4 = RandomForestClassifier(\n    class_weight='balanced',\n    n_jobs=-1,\n    random_state=0\n)\nclf_hsa = RandomForestClassifier(\n    class_weight='balanced',\n    n_jobs=-1,\n    random_state=0\n)","metadata":{"execution":{"iopub.status.busy":"2024-05-01T17:51:29.916251Z","iopub.status.idle":"2024-05-01T17:51:29.917359Z","shell.execute_reply.started":"2024-05-01T17:51:29.917136Z","shell.execute_reply":"2024-05-01T17:51:29.917156Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nclf_seh.fit(X_seh_transformed, y_seh)\nclf_brd4.fit(X_brd4_transformed, y_brd4)\nclf_hsa.fit(X_hsa_transformed, y_hsa)","metadata":{"execution":{"iopub.status.busy":"2024-05-01T17:51:29.918999Z","iopub.status.idle":"2024-05-01T17:51:29.919398Z","shell.execute_reply.started":"2024-05-01T17:51:29.919207Z","shell.execute_reply":"2024-05-01T17:51:29.919224Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"preds_seh = clf_seh.predict_proba(X_valid_seh_transformed)[:, 1]\npreds_brd4 = clf_brd4.predict_proba(X_valid_brd4_transformed)[:, 1]\npreds_hsa = clf_hsa.predict_proba(X_valid_hsa_transformed)[:, 1]","metadata":{"execution":{"iopub.status.busy":"2024-05-01T17:51:29.920882Z","iopub.status.idle":"2024-05-01T17:51:29.921285Z","shell.execute_reply.started":"2024-05-01T17:51:29.921087Z","shell.execute_reply":"2024-05-01T17:51:29.921103Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"preds_concatenated = np.concatenate([\n    preds_seh,\n    preds_brd4,\n    preds_hsa\n])\n\ny_valid_concatenated = np.concatenate([\n    y_valid_seh,\n    y_valid_brd4,\n    y_valid_hsa\n])\n\nprint(f'Average precision score: {average_precision_score(y_valid_concatenated, preds_concatenated, average=\"micro\") * 100 :.2f}%')","metadata":{"execution":{"iopub.status.busy":"2024-05-01T17:51:29.923381Z","iopub.status.idle":"2024-05-01T17:51:29.923794Z","shell.execute_reply.started":"2024-05-01T17:51:29.923588Z","shell.execute_reply":"2024-05-01T17:51:29.923603Z"},"trusted":true},"execution_count":null,"outputs":[]}]}