{"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"},{"sourceId":8196925,"sourceType":"datasetVersion","datasetId":4855347},{"sourceId":8752384,"sourceType":"datasetVersion","datasetId":5205859}],"dockerImageVersionId":30732,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"Three binary targets (protein binding affinity):\n- sEH\n- BRD4\n- HSA (ALB)\n\nA number of features that need to be extracted from a molecule (given in SMILES format) through some calculations.\n","metadata":{}},{"cell_type":"markdown","source":"## Environment setup, package installs, library imports\n---","metadata":{}},{"cell_type":"code","source":"#!pip install rdkit","metadata":{"execution":{"iopub.status.busy":"2024-06-22T15:39:47.108720Z","iopub.execute_input":"2024-06-22T15:39:47.109233Z","iopub.status.idle":"2024-06-22T15:39:47.115136Z","shell.execute_reply.started":"2024-06-22T15:39:47.109194Z","shell.execute_reply":"2024-06-22T15:39:47.113888Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#!nvidia-smi  # this should display information about available GPUs","metadata":{"execution":{"iopub.status.busy":"2024-06-22T15:39:47.119256Z","iopub.execute_input":"2024-06-22T15:39:47.119722Z","iopub.status.idle":"2024-06-22T15:39:47.128283Z","shell.execute_reply.started":"2024-06-22T15:39:47.119682Z","shell.execute_reply":"2024-06-22T15:39:47.126939Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#%load_ext cudf.pandas\n\nimport os\nimport gc\nimport numpy as np \nimport pandas as pd\nimport joblib\n\n# import rdkit\n# from rdkit.Chem import (AllChem, MolFromSmiles, rdMolDescriptors, \n#                         Descriptors, rdmolfiles, MolFromPDBFile)\n \nimport matplotlib.pyplot as plt\nimport seaborn as sns\n# import lightgbm as lgb\nfrom sklearn.model_selection import cross_val_score, GridSearchCV\nfrom sklearn.metrics import average_precision_score\nfrom sklearn.metrics.pairwise import cosine_similarity\nfrom sklearn.metrics import (accuracy_score, precision_score, \n                             recall_score, f1_score, classification_report)\nfrom xgboost import XGBClassifier\nimport shap\n\nfor dirname, _, filenames in os.walk('/kaggle/input'):\n    for filename in filenames:\n        print(os.path.join(dirname, filename))","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2024-06-22T15:39:47.130384Z","iopub.execute_input":"2024-06-22T15:39:47.130740Z","iopub.status.idle":"2024-06-22T15:39:47.149196Z","shell.execute_reply.started":"2024-06-22T15:39:47.130711Z","shell.execute_reply":"2024-06-22T15:39:47.147743Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# import cudf\n# print(\"cuDF version: \", cudf.__version__)","metadata":{"execution":{"iopub.status.busy":"2024-06-22T15:39:47.151496Z","iopub.execute_input":"2024-06-22T15:39:47.151956Z","iopub.status.idle":"2024-06-22T15:39:47.158423Z","shell.execute_reply.started":"2024-06-22T15:39:47.151917Z","shell.execute_reply":"2024-06-22T15:39:47.157236Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"base_path = \"/kaggle/input/belka-training-data\"\nsub_file = \"/kaggle/input/leash-BELKA/sample_submission.csv\"\nproteins = ['sEH', 'BRD4', 'HSA']","metadata":{"execution":{"iopub.status.busy":"2024-06-22T15:39:47.160801Z","iopub.execute_input":"2024-06-22T15:39:47.161205Z","iopub.status.idle":"2024-06-22T15:39:47.171329Z","shell.execute_reply.started":"2024-06-22T15:39:47.161174Z","shell.execute_reply":"2024-06-22T15:39:47.170067Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Data loading and sampling\n---\n\nThe train dataset is way too big so it cannot be read in full neither as csv nor parquet. Let´s try the DuckDB trick.... \n\nCurrently using a balanced (-ish) sample of 20K records.","metadata":{}},{"cell_type":"markdown","source":"## Generate features\n---","metadata":{}},{"cell_type":"code","source":"def reduce_dtypes(df, train=False):\n    \n    if train :\n        df['binds'] = df['binds'].astype(np.int8)\n        \n    df['molecule_NumRotatableBonds']   = df['molecule_NumRotatableBonds'].astype(np.int8)\n    df['molecule_NumHDonors']          = df['molecule_NumHDonors'].astype(np.int8)\n    df['molecule_NumHAcceptors']       = df['molecule_NumHAcceptors'].astype(np.int8)\n    df['molecule_NumAromaticRings']    = df['molecule_NumAromaticRings'].astype(np.int8)\n    df['molecule_NumSaturatedRings']   = df['molecule_NumSaturatedRings'].astype(np.int8)\n    df['molecule_NumAliphaticRings']   = df['molecule_NumAliphaticRings'].astype(np.int8)\n    df['molecule_NumValenceElectrons'] = df['molecule_NumValenceElectrons'].astype(np.int8)\n    df['molecule_NumRadicalElectrons'] = df['molecule_NumRadicalElectrons'].astype(np.int8)\n\n    df['molecule_ExactMolWt']     = df['molecule_ExactMolWt'].astype(np.float32)\n    df['molecule_TPSA']           = df['molecule_TPSA'].astype(np.float32)\n    df['molecule_NumHeteroatoms'] = df['molecule_NumHeteroatoms'].astype(np.float32)\n    #df['molecule_mfp']            = df['molecule_mfp'].astype(np.float32)\n    \n    df['protein_NumRotatableBonds']   = df['protein_NumRotatableBonds'].astype(np.int8)\n    df['protein_NumHDonors']          = df['protein_NumHDonors'].astype(np.int8)\n    df['protein_NumHAcceptors']       = df['protein_NumHAcceptors'].astype(np.int8)\n    df['protein_NumAromaticRings']    = df['protein_NumAromaticRings'].astype(np.int8)\n    df['protein_NumSaturatedRings']   = df['protein_NumSaturatedRings'].astype(np.int8)\n    df['protein_NumAliphaticRings']   = df['protein_NumAliphaticRings'].astype(np.int8)\n    df['protein_NumValenceElectrons'] = df['protein_NumValenceElectrons'].astype(np.int8)\n    df['protein_NumRadicalElectrons'] = df['protein_NumRadicalElectrons'].astype(np.int8)\n\n    df['protein_ExactMolWt']     = df['protein_ExactMolWt'].astype(np.float32)\n    df['protein_TPSA']           = df['protein_TPSA'].astype(np.float32)\n    df['protein_NumHeteroatoms'] = df['protein_NumHeteroatoms'].astype(np.float32)\n    #df['protein_mfp']            = df['protein_mfp'].astype(np.float32)\n    \n    df['cosine_sim']              = df['cosine_sim'].astype(np.float32)\n    \n    return df","metadata":{"execution":{"iopub.status.busy":"2024-06-22T15:39:47.173722Z","iopub.execute_input":"2024-06-22T15:39:47.174135Z","iopub.status.idle":"2024-06-22T15:39:47.192029Z","shell.execute_reply.started":"2024-06-22T15:39:47.174101Z","shell.execute_reply":"2024-06-22T15:39:47.190749Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Create dictionary of train datasets and separate data\ntrain_df = {}\n\nfor p in proteins:\n    print(f\"reading train data for {p}\")\n    train_df[p] = pd.read_csv(os.path.join(base_path, f\"train_{p}.csv\"))\n    train_df[p] = reduce_dtypes(train_df[p], train=True)\n    print(train_df[p].info())","metadata":{"execution":{"iopub.status.busy":"2024-06-22T15:39:47.193498Z","iopub.execute_input":"2024-06-22T15:39:47.193951Z","iopub.status.idle":"2024-06-22T15:40:16.553374Z","shell.execute_reply.started":"2024-06-22T15:39:47.193919Z","shell.execute_reply":"2024-06-22T15:40:16.552146Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Run grid-search cross-validation and obtain optimal parameters for each of the models\n---","metadata":{}},{"cell_type":"code","source":"features = ['molecule_NumRotatableBonds', 'molecule_NumHDonors', 'molecule_NumHAcceptors', 'molecule_NumAromaticRings',\n            'molecule_NumSaturatedRings', 'molecule_NumAliphaticRings', 'molecule_NumHeteroatoms', 'molecule_NumValenceElectrons',\n            'molecule_NumRadicalElectrons', 'molecule_ExactMolWt', 'molecule_TPSA', \n            \n            'protein_NumRotatableBonds', 'protein_NumHDonors', 'protein_NumHAcceptors', 'protein_NumAromaticRings',\n            'protein_NumSaturatedRings', 'protein_NumAliphaticRings', 'protein_NumHeteroatoms', 'protein_NumValenceElectrons',\n            'protein_NumRadicalElectrons', 'protein_ExactMolWt', 'protein_TPSA',\n            \n            'cosine_sim']\n\ntarget = 'binds'","metadata":{"execution":{"iopub.status.busy":"2024-06-22T15:40:16.555673Z","iopub.execute_input":"2024-06-22T15:40:16.556065Z","iopub.status.idle":"2024-06-22T15:40:16.562723Z","shell.execute_reply.started":"2024-06-22T15:40:16.556011Z","shell.execute_reply":"2024-06-22T15:40:16.561456Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Define dictionaries (keys are the proteins)\nX = {}\ny = {}\nspw = {}\nmodel = {}\noptimal_params = {}","metadata":{"execution":{"iopub.status.busy":"2024-06-22T15:40:16.564405Z","iopub.execute_input":"2024-06-22T15:40:16.564783Z","iopub.status.idle":"2024-06-22T15:40:16.582102Z","shell.execute_reply.started":"2024-06-22T15:40:16.564752Z","shell.execute_reply":"2024-06-22T15:40:16.580797Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import xgboost as xgb\n\noptimal_params = {\n    'HSA' : {\n        'n_estimators': 1500,\n        'learning_rate': 0.01,\n        'max_depth': 7\n    },\n    'BRD4' : {\n        'n_estimators': 2500,\n        'learning_rate': 0.05,\n        'max_depth': 9\n    },\n    'sEH' : {\n        'n_estimators': 2000,\n        'learning_rate': 0.05,\n        'max_depth': 9\n    }\n}\n\noptimal_params","metadata":{"execution":{"iopub.status.busy":"2024-06-22T15:40:16.583658Z","iopub.execute_input":"2024-06-22T15:40:16.584089Z","iopub.status.idle":"2024-06-22T15:40:16.600183Z","shell.execute_reply.started":"2024-06-22T15:40:16.584028Z","shell.execute_reply":"2024-06-22T15:40:16.598951Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"raw","source":"import xgboost as xgb\n\nsX = {}\nsy = {}\nspw = {}\nfor p in proteins:\n    # Take a sample of the train data for locatiung optimal hyper parameters\n    sample_df = train_df[p].sample(75_000)\n    sX[p] = sample_df[features]\n    sy[p] = sample_df[target]\n    spw[p] = sample_df[target].value_counts()[0] / sample_df[target].value_counts()[1]","metadata":{"execution":{"iopub.status.busy":"2024-06-21T18:21:43.095524Z","iopub.execute_input":"2024-06-21T18:21:43.095883Z","iopub.status.idle":"2024-06-21T18:21:43.374003Z","shell.execute_reply.started":"2024-06-21T18:21:43.095849Z","shell.execute_reply":"2024-06-21T18:21:43.372919Z"}}},{"cell_type":"raw","source":"# Parameters grid for the cross-validation exercise\nhyperparam_grid = {\n    'n_estimators': [1500, 2000, 2500],\n    'learning_rate': [0.01, 0.03, 0.05, 0.07, 0.1],\n    'max_depth': [3, 5, 7, 9],\n#     'n_estimators': [400, 600, 800, 1000, 1200, 1500],\n#     'learning_rate': [0.05, 0.07, 0.1, 0.13, 0.15],\n#     'max_depth': [3, 5, 7, 9, 10, 11],\n#     'n_estimators': [100, 150, 200, 250, 300, 350, 400, 450, 500],\n#     'learning_rate': [0.05, 0.07, 0.1, 0.13, 0.15],\n#     'max_depth': [3, 5, 7, 9],\n#     'subsample': [0, 0.2, 0.4, 0.6, 0.8, 1],\n#     'colsample_bytree': [0, 0.2, 0.4, 0.6, 0.8, 1],\n#     'min_child_weight': [0, 0.2, 0.4, 0.6, 0.8, 1],\n    }","metadata":{"execution":{"iopub.status.busy":"2024-06-21T18:21:43.375251Z","iopub.execute_input":"2024-06-21T18:21:43.375576Z","iopub.status.idle":"2024-06-21T18:21:43.381350Z","shell.execute_reply.started":"2024-06-21T18:21:43.375549Z","shell.execute_reply":"2024-06-21T18:21:43.380292Z"}}},{"cell_type":"raw","source":"%%time\n# Create models and run grid search cross-validation\nfor p in proteins:\n    # Instantiate XGBoost model\n    model[p] = XGBClassifier(scale_pos_weight=spw[p],\n                             random_state=13,\n                             tree_method='hist',\n                             #device='cuda'\n                            )\n    print('Model', p)\n    print('Running grid search cross validation....')\n    # Set up the gscv object with 4-fold\n    gs_cv = GridSearchCV(estimator=model[p],\n                         param_grid=hyperparam_grid,\n                         scoring='average_precision',\n                         cv=4,\n                         return_train_score=True,\n                         n_jobs=-1,\n                         verbose=1)\n    \n    gs_cv.fit(sX[p],sy[p])\n    optimal_params[p] = gs_cv.best_params_\n    print('Optimal params for model', p)\n    print(optimal_params[p])\n    print('-'*66)","metadata":{"execution":{"iopub.status.busy":"2024-06-21T18:21:43.385880Z","iopub.execute_input":"2024-06-21T18:21:43.386629Z","iopub.status.idle":"2024-06-21T20:25:02.162730Z","shell.execute_reply.started":"2024-06-21T18:21:43.386589Z","shell.execute_reply":"2024-06-21T20:25:02.160512Z"}}},{"cell_type":"raw","source":"del sX, sy\ngc.collect()","metadata":{"execution":{"iopub.status.busy":"2024-06-21T20:25:02.165553Z","iopub.execute_input":"2024-06-21T20:25:02.166066Z","iopub.status.idle":"2024-06-21T20:25:02.408788Z","shell.execute_reply.started":"2024-06-21T20:25:02.166016Z","shell.execute_reply":"2024-06-21T20:25:02.407637Z"}}},{"cell_type":"raw","source":"# Is XGBoost GPU-Enabled?\nimport xgboost as xgb\nprint(xgb.rabit.get_rank() == 0)","metadata":{"execution":{"iopub.status.busy":"2024-06-21T20:25:02.410499Z","iopub.execute_input":"2024-06-21T20:25:02.410839Z","iopub.status.idle":"2024-06-21T20:25:02.426573Z","shell.execute_reply.started":"2024-06-21T20:25:02.410811Z","shell.execute_reply":"2024-06-21T20:25:02.425203Z"}}},{"cell_type":"markdown","source":"## Fit optimal models with full sample and calculate SHAP values\n---\n\nNote that the column *molecule_NumRadicalElectrons* has been removed from the training set in this run, following an earlier run which showed it was not relevant.","metadata":{}},{"cell_type":"code","source":"# Fit optimal models\nbest_model = {}\nfor p in proteins:\n    X = train_df[p][features]\n    y = train_df[p][target]\n    spw = train_df[p][target].value_counts()[0] / train_df[p][target].value_counts()[1]\n       \n    #booster.set_param({\"device\": \"cuda:0\"})\n    best_model[p] = XGBClassifier(**optimal_params[p],\n                                  scale_pos_weight=spw,\n                                  random_state=13,\n                                  tree_method='hist',\n                                  #device = \"cuda\"\n                                 )    \n    best_model[p].fit(X, y)\n    cv_scores = cross_val_score(best_model[p], X, y, cv=4, scoring='average_precision')\n    print('Cross validation score for model', p)\n    print('Average precision:', str(round(cv_scores.mean(), 3)))\n    print('-'*40)\n    \n    # Save the model\n    joblib.dump(best_model[p], f'{p}_model.pkl')","metadata":{"execution":{"iopub.status.busy":"2024-06-22T15:40:16.601914Z","iopub.execute_input":"2024-06-22T15:40:16.603146Z","iopub.status.idle":"2024-06-22T16:10:08.515538Z","shell.execute_reply.started":"2024-06-22T15:40:16.603099Z","shell.execute_reply":"2024-06-22T16:10:08.514238Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"raw","source":"%%time\n# Calculate and show Shap Values for the full sample\nshap.initjs()\nfor p in proteins:\n    explainer = shap.TreeExplainer(best_model[p])\n    shap_values = explainer.shap_values(X[p])\n    plt.suptitle('Model ' + p, fontsize=15) \n    shap.summary_plot(shap_values, X[p], plot_size=[8,8], plot_type='dot', show=False)\n    plt.show()","metadata":{"execution":{"iopub.status.busy":"2024-06-18T15:19:56.073735Z","iopub.execute_input":"2024-06-18T15:19:56.074248Z"}}},{"cell_type":"markdown","source":"## Processing of test dataset and preparation for  (TO BE CONTINUED...)\n---","metadata":{}},{"cell_type":"raw","source":"# Read the model files\nbest_model = {}\nfor p in proteins:  # Assuming 'proteins' is a list of model names or identifiers\n    print(f\"Reading model {p}\")\n    best_model[p] = joblib.load(f'{p}_model.pkl')","metadata":{"execution":{"iopub.status.busy":"2024-06-15T19:43:24.219578Z","iopub.execute_input":"2024-06-15T19:43:24.220261Z","iopub.status.idle":"2024-06-15T19:43:24.306099Z","shell.execute_reply.started":"2024-06-15T19:43:24.220229Z","shell.execute_reply":"2024-06-15T19:43:24.304900Z"}}},{"cell_type":"code","source":"test_df = {}\nfor p in proteins:\n    print(f\"reading test data for {p}\")\n    test_df[p] = pd.read_csv(os.path.join(base_path, f\"test_{p}.csv\"))\n    test_df[p] = reduce_dtypes(test_df[p], train=False)\n    print(test_df[p].info())","metadata":{"execution":{"iopub.status.busy":"2024-06-22T16:10:08.518625Z","iopub.execute_input":"2024-06-22T16:10:08.518991Z","iopub.status.idle":"2024-06-22T16:10:33.543644Z","shell.execute_reply.started":"2024-06-22T16:10:08.518959Z","shell.execute_reply":"2024-06-22T16:10:33.542398Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Prepare test data\npredictions = {}\nfor p in proteins:\n    X_test = test_df[p][features]\n    predictions[p] = best_model[p].predict_proba(X_test)    ","metadata":{"execution":{"iopub.status.busy":"2024-06-22T16:10:33.545296Z","iopub.execute_input":"2024-06-22T16:10:33.545639Z","iopub.status.idle":"2024-06-22T16:11:38.464094Z","shell.execute_reply.started":"2024-06-22T16:10:33.545610Z","shell.execute_reply":"2024-06-22T16:11:38.462984Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Create submission DataFrame\nsubmission_df = pd.DataFrame()\nfor p in proteins:\n    sub_df = {\n        'id': test_df[p]['id'],\n        'binds': np.max(predictions[p], axis=1)\n    }\n    sub_df = pd.DataFrame(sub_df)\n    submission_df = pd.concat([submission_df, sub_df], ignore_index=True)","metadata":{"execution":{"iopub.status.busy":"2024-06-22T16:11:38.465398Z","iopub.execute_input":"2024-06-22T16:11:38.468273Z","iopub.status.idle":"2024-06-22T16:11:38.496649Z","shell.execute_reply.started":"2024-06-22T16:11:38.468233Z","shell.execute_reply":"2024-06-22T16:11:38.495656Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Restore the original order in test file","metadata":{}},{"cell_type":"code","source":"orig_ids = pd.read_csv('/kaggle/input/leash-BELKA/test.csv', usecols=['id'])\nsubmission_df = pd.merge(orig_ids, submission_df, on='id')","metadata":{"execution":{"iopub.status.busy":"2024-06-22T16:11:38.497945Z","iopub.execute_input":"2024-06-22T16:11:38.498302Z","iopub.status.idle":"2024-06-22T16:11:41.303925Z","shell.execute_reply.started":"2024-06-22T16:11:38.498273Z","shell.execute_reply":"2024-06-22T16:11:41.302470Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(\"length of submisison fiule : \", submission_df.shape)\nprint(submission_df.head())","metadata":{"execution":{"iopub.status.busy":"2024-06-22T16:11:41.305355Z","iopub.execute_input":"2024-06-22T16:11:41.305684Z","iopub.status.idle":"2024-06-22T16:11:41.313513Z","shell.execute_reply.started":"2024-06-22T16:11:41.305657Z","shell.execute_reply":"2024-06-22T16:11:41.312473Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"submission_df.to_csv(\"submission.csv\", index=False)","metadata":{"execution":{"iopub.status.busy":"2024-06-22T16:11:41.314772Z","iopub.execute_input":"2024-06-22T16:11:41.315154Z","iopub.status.idle":"2024-06-22T16:11:45.557963Z","shell.execute_reply.started":"2024-06-22T16:11:41.315112Z","shell.execute_reply":"2024-06-22T16:11:45.556769Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"raw","source":"bind = sub['binds'].value_counts()\nprint(\"distribution\\n\")\nprint(bind)\nprint(f\"percent binding :  {int((bind[1]/bind[0])*100)}%\")","metadata":{"execution":{"iopub.status.busy":"2024-06-21T22:56:10.408023Z","iopub.execute_input":"2024-06-21T22:56:10.408836Z","iopub.status.idle":"2024-06-21T22:56:10.433116Z","shell.execute_reply.started":"2024-06-21T22:56:10.408797Z","shell.execute_reply":"2024-06-21T22:56:10.432023Z"}}},{"cell_type":"markdown","source":"The line below takes ages as the test dataset is huge too - I aborted after 40 minutes without sign of finishing!","metadata":{}},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}