{"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":30699,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"

Putting on my chemist robe and glasses :)

Let us first figure out what they want from us and what we want from our code.

\n\n\nThe dataset comprises binary classification data representing whether small molecules bind to three different protein targets. Each entry includes SMILES representations of molecule structures and binary labels for binding to each protein target. The competition data, provided by Leash Biosciences, consists of approximately **98M** training examples per protein, **200K** validation examples per protein, and **360K** test molecules per protein. *Keep in mind that dataset is very imbalanced: roughly 0.5% of examples are classified as binders*

\n\n

So, what is SMILES?


\n\nSMILES (Simplified Molecular Input Line Entry System) is a notation system for representing chemical structures in a computer-readable format. Developed with funding from the U.S. Environmental Protection Agency, it offers a flexible and easily learned approach to representing molecules. SMILES notation follows five basic syntax rules:\n\n* Atoms and Bonds: Atoms are represented by their atomic symbols, with lowercase letters indicating aromatic atoms. Bonds are denoted by symbols (- for single, = for double, # for triple, * for aromatic, and . for disconnected structures).\n* Simple Chains: Chains of atoms are represented by combining atomic symbols and bond symbols. Hydrogen atoms are suppressed unless explicitly stated.\n* Branches: Branches from chains are enclosed in parentheses and placed directly after the atom to which they are connected.\n* Rings: Ring structures are identified by using numbers to denote the opening and closing ring atoms. Different numbers are used for each ring, and bond symbols may precede the ring closure number.\n* Charged Atoms: Charges on atoms are indicated by placing the atom symbol within brackets, enclosing the charge.
\n\n\n

What are the targets?


\n\nThe dataset includes three protein targets:\n

EPHX2 (sEH):


\n* This target refers to epoxide hydrolase 2, encoded by the EPHX2 genetic locus. Its protein product, commonly known as soluble epoxide hydrolase (sEH), is an enzyme that catalyzes certain chemical reactions and hydrolyzes phosphate groups. It is a potential drug target for conditions like high blood pressure and diabetes. The dataset includes screening data obtained from Leash Biosciences, along with structural information for model evaluation.\n

BRD4:


\n* Bromodomain 4 is encoded by the BRD4 locus. Its protein product, also named BRD4, plays a role in gene transcription regulation by binding to histones in the nucleus. This protein is implicated in cancer progression, and inhibiting its activity has been explored as a therapeutic strategy. The dataset includes screening data from Leash Biosciences and structural information for model evaluation.\n

ALB (HSA):


\n* Serum albumin, encoded by the ALB locus, is the most abundant protein in blood. It regulates osmotic pressure and transports various molecules, including drugs, hormones, and fatty acids. Predicting the binding of small molecules to albumin is crucial for drug development, as it impacts drug distribution and effectiveness. The dataset includes screening data from Leash Biosciences and structural information for model evaluation.\n","metadata":{}},{"cell_type":"markdown","source":"

Now let's get the things done! :D


\n\nWe shall start with the installation of necessary tools.

\n**DuckDB is an open-source, in-memory database system optimized for fast analytical querying and low memory usage. It's lightweight, embeddable, and seamlessly integrates with popular programming languages like Python and R. DuckDB excels in executing rapid SQL queries on large datasets, making it ideal for data exploration and interactive analysis tasks.**","metadata":{}},{"cell_type":"code","source":"!pip install duckdb","metadata":{"_kg_hide-output":true,"execution":{"iopub.status.busy":"2024-04-22T16:34:04.028473Z","iopub.execute_input":"2024-04-22T16:34:04.028838Z","iopub.status.idle":"2024-04-22T16:34:18.706684Z","shell.execute_reply.started":"2024-04-22T16:34:04.028806Z","shell.execute_reply":"2024-04-22T16:34:18.705231Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**RDKit is a widely-used open-source toolkit for cheminformatics. It offers various functions for working with chemical structures and data, including molecular representation, substructure searching, and compound library handling. RDKit is popular in pharmaceutical research and drug discovery due to its comprehensive features and easy integration with Python.**","metadata":{}},{"cell_type":"code","source":"!pip install rdkit","metadata":{"_kg_hide-output":true,"execution":{"iopub.status.busy":"2024-04-22T16:34:40.837391Z","iopub.execute_input":"2024-04-22T16:34:40.837753Z","iopub.status.idle":"2024-04-22T16:34:58.022609Z","shell.execute_reply.started":"2024-04-22T16:34:40.837726Z","shell.execute_reply":"2024-04-22T16:34:58.021203Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\nimport numpy as np\nfrom tqdm import tqdm\nimport duckdb\nfrom rdkit import Chem\nfrom rdkit.Chem import AllChem\nfrom sklearn.model_selection import train_test_split\nfrom sklearn.preprocessing import OneHotEncoder\nfrom xgboost import XGBClassifier\nfrom catboost import CatBoostClassifier\nfrom lightgbm import LGBMClassifier\nfrom sklearn.metrics import average_precision_score\n\nfrom sklearn.svm import SVC\nfrom sklearn.ensemble import StackingClassifier, VotingClassifier, BaggingClassifier, RandomForestClassifier\nfrom sklearn.neural_network import MLPClassifier\nfrom sklearn.linear_model import LogisticRegression\nfrom sklearn.calibration import CalibratedClassifierCV","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2024-04-22T16:55:25.772449Z","iopub.execute_input":"2024-04-22T16:55:25.772849Z","iopub.status.idle":"2024-04-22T16:55:28.807907Z","shell.execute_reply.started":"2024-04-22T16:55:25.772818Z","shell.execute_reply":"2024-04-22T16:55:28.806913Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\ntrain_path = '/kaggle/input/leash-BELKA/train.parquet'\ntest_path = '/kaggle/input/leash-BELKA/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 30000)\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-04-22T16:55:30.767538Z","iopub.execute_input":"2024-04-22T16:55:30.768352Z","iopub.status.idle":"2024-04-22T16:56:17.568853Z","shell.execute_reply.started":"2024-04-22T16:55:30.768315Z","shell.execute_reply":"2024-04-22T16:56:17.566745Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df.head()","metadata":{"execution":{"iopub.status.busy":"2024-04-22T16:56:28.892186Z","iopub.execute_input":"2024-04-22T16:56:28.893127Z","iopub.status.idle":"2024-04-22T16:56:28.908540Z","shell.execute_reply.started":"2024-04-22T16:56:28.893070Z","shell.execute_reply":"2024-04-22T16:56:28.907305Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Further, we convert SMILES representations of molecules into RDKit molecule objects and then generates ECFP fingerprints for each molecule, storing the results in the DataFrame. \n\n**ECFP, or Extended Connectivity Fingerprints**, are molecular fingerprints used in cheminformatics to encode structural features of molecules. They represent connectivity patterns of atoms by hashing local environments within a specified radius. ECFP fingerprints are widely used for similarity searching and activity prediction in drug discovery and cheminformatics due to their efficiency and robustness.\n\n

","metadata":{}},{"cell_type":"code","source":"%%time\n# Convert SMILES to RDKit molecules\ndf['molecule'] = df['molecule_smiles'].apply(Chem.MolFromSmiles)\n\n# 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\ndf['ecfp'] = df['molecule'].apply(generate_ecfp)","metadata":{"execution":{"iopub.status.busy":"2024-04-22T16:56:32.228178Z","iopub.execute_input":"2024-04-22T16:56:32.230812Z","iopub.status.idle":"2024-04-22T16:58:07.739894Z","shell.execute_reply.started":"2024-04-22T16:56:32.230772Z","shell.execute_reply":"2024-04-22T16:58:07.737888Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df['ecfp'].head()","metadata":{"execution":{"iopub.status.busy":"2024-04-22T16:59:11.549391Z","iopub.execute_input":"2024-04-22T16:59:11.549805Z","iopub.status.idle":"2024-04-22T16:59:11.562350Z","shell.execute_reply.started":"2024-04-22T16:59:11.549775Z","shell.execute_reply":"2024-04-22T16:59:11.561273Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"RANDOM_STATE = 42","metadata":{"execution":{"iopub.status.busy":"2024-04-22T16:59:13.627464Z","iopub.execute_input":"2024-04-22T16:59:13.628176Z","iopub.status.idle":"2024-04-22T16:59:13.632868Z","shell.execute_reply.started":"2024-04-22T16:59:13.628133Z","shell.execute_reply":"2024-04-22T16:59:13.631811Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\n# 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, validation and test sets\nX_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=RANDOM_STATE)\n\nX_train, X_val, y_train, y_val = train_test_split(X_train, y_train, test_size=0.25, random_state=RANDOM_STATE)","metadata":{"execution":{"iopub.status.busy":"2024-04-22T16:59:15.836886Z","iopub.execute_input":"2024-04-22T16:59:15.837290Z","iopub.status.idle":"2024-04-22T16:59:17.596719Z","shell.execute_reply.started":"2024-04-22T16:59:15.837260Z","shell.execute_reply":"2024-04-22T16:59:17.595616Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# # Define base models\n# rf_model = RandomForestClassifier()\n# xgb_model = XGBClassifier()\n# lgbm_model = LGBMClassifier()\n# catboost_model = CatBoostClassifier()\n# svm_model = SVC(kernel='rbf') # RBF kernel SVM\n# mlp_model = MLPClassifier() # DNN\n\n# # Define ensemble methods\n# estimators = [\n# ('rf', rf_model),\n# ('xgb', xgb_model),\n# ('lgbm', lgbm_model),\n# ('catboost', catboost_model),\n# ('svm', svm_model),\n# ('mlp', mlp_model)\n# ]\n\n# # Stacking\n# stacking_model = StackingClassifier(estimators=estimators, final_estimator=LogisticRegression())\n\n# # Blending\n# blending_model = VotingClassifier(estimators=estimators, voting='soft')\n\n# # Bagging\n# bagging_model = BaggingClassifier(base_estimator=rf_model, n_estimators=10)\n\n# # Train and evaluate each model\n# models = {\n# 'Random Forest': rf_model,\n# 'XGBoost': xgb_model,\n# 'LightGBM': lgbm_model,\n# 'CatBoost': catboost_model,\n# #'SVM (RBF Kernel)': svm_model,\n# 'DNN': mlp_model,\n# 'Stacking': stacking_model,\n# 'Blending': blending_model,\n# 'Bagging': bagging_model\n# }\n\n# for name, model in models.items():\n# model.fit(X_train, y_train)\n# y_pred = model.predict_proba(X_test)[:, 1]\n# accuracy = average_precision_score(y_test, y_pred)\n# print(f'{name} Accuracy: {accuracy}')","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"NUMBER_OF_MODELS = 3\n# Create and train the ensemble models\ncatboost_model = CatBoostClassifier(iterations=500, random_state=RANDOM_STATE)\nlgbm_model = LGBMClassifier(random_state=RANDOM_STATE)\nxgb_model = XGBClassifier(random_state=RANDOM_STATE, n_jobs=-1,)\n\nmodels = [catboost_model, lgbm_model, xgb_model]\nfor model in models:\n model.fit(X_train, y_train)","metadata":{"_kg_hide-output":true,"execution":{"iopub.status.busy":"2024-04-22T16:59:32.006406Z","iopub.execute_input":"2024-04-22T16:59:32.007152Z","iopub.status.idle":"2024-04-22T17:01:07.815723Z","shell.execute_reply.started":"2024-04-22T16:59:32.007104Z","shell.execute_reply":"2024-04-22T17:01:07.814741Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"predictions_test = []\n# Calibrate probabilities on validation set\nfor model in models:\n calibrated_model = CalibratedClassifierCV(model, method='sigmoid', cv='prefit')\n calibrated_model.fit(X_val, y_val)\n # Evaluate calibrated model on test set\n calibrated_probabilities = calibrated_model.predict_proba(X_test)[:, 1]\n predictions_test.append(calibrated_probabilities)\n \n# Ensemble predictions for the test set\nensemble_predictions_test = np.mean(predictions_test, axis=0)","metadata":{"execution":{"iopub.status.busy":"2024-04-22T17:01:16.879293Z","iopub.execute_input":"2024-04-22T17:01:16.879685Z","iopub.status.idle":"2024-04-22T17:01:53.407142Z","shell.execute_reply.started":"2024-04-22T17:01:16.879656Z","shell.execute_reply":"2024-04-22T17:01:53.404703Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Calculate the mean average precision\nmap_score = average_precision_score(y_test, ensemble_predictions_test)\nprint(f\"Mean Average Precision (mAP): {map_score:.8f}\")","metadata":{"execution":{"iopub.status.busy":"2024-04-22T17:01:58.990177Z","iopub.execute_input":"2024-04-22T17:01:58.990870Z","iopub.status.idle":"2024-04-22T17:01:59.017413Z","shell.execute_reply.started":"2024-04-22T17:01:58.990793Z","shell.execute_reply":"2024-04-22T17:01:59.015976Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import os\n# Process the test.parquet file chunk by chunk sinc the file is huge\ntest_file = '/kaggle/input/leash-BELKA/test.csv'\noutput_file = 'submission.csv' # Specify the path and filename for the output file\n\n\n\n# Read the test.csv file into a pandas DataFrame\n# Wrap the loop with tqdm to display progress\nfor df_test in tqdm(pd.read_csv(test_file, chunksize=10_000)):\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 predictions_test = []\n # Calibrate probabilities on validation set\n for model in models:\n calibrated_model = CalibratedClassifierCV(model, method='sigmoid', cv='prefit')\n calibrated_model.fit(X_val, y_val)\n # Evaluate calibrated model on test set\n calibrated_probabilities = calibrated_model.predict_proba(X_test)[:, 1]\n predictions_test.append(calibrated_probabilities)\n\n # Ensemble predictions for the test set\n ensemble_predictions_test = np.mean(predictions_test, axis=0)\n\n# # Predict the probabilities using ensemble models\n# probabilities_catboost = catboost_model.predict_proba(X_test)[:, 1]\n# probabilities_lgbm = lgbm_model.predict_proba(X_test)[:, 1]\n# probabilities_xgb = xgb_model.predict_proba(X_test)[:, 1]\n\n# # Average the predictions from the ensemble\n# probabilities_ensemble = (probabilities_catboost + probabilities_lgbm + probabilities_xgb) / NUMBER_OF_MODELS\n \n # Append the probabilities to the list\n\n# Create a DataFrame with 'id' and 'probability' columns\n output_df = pd.DataFrame({'id': df_test['id'], 'binds': ensemble_predictions_test})\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))","metadata":{"execution":{"iopub.status.busy":"2024-04-22T17:02:02.674495Z","iopub.execute_input":"2024-04-22T17:02:02.675215Z","iopub.status.idle":"2024-04-22T19:00:14.292189Z","shell.execute_reply.started":"2024-04-22T17:02:02.675169Z","shell.execute_reply":"2024-04-22T19:00:14.290891Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}