{"metadata":{"kernelspec":{"display_name":"Python 3 (ipykernel)","language":"python","name":"python3"},"language_info":{"codemirror_mode":{"name":"ipython","version":3},"file_extension":".py","mimetype":"text/x-python","name":"python","nbconvert_exporter":"python","pygments_lexer":"ipython3","version":"3.10.12"},"kaggle":{"accelerator":"none","dataSources":[{"sourceId":67356,"databundleVersionId":8006601,"sourceType":"competition"},{"sourceId":8042988,"sourceType":"datasetVersion","datasetId":4740586}],"dockerImageVersionId":30698,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# BELKA EDA and data split\n\nMy main focus in this notebook is to show my creation of a dataset with CV split that intends to mirror the test set split created by the host. I ran this notebook locally so you may not reproduce it on Kaggle, but I still want to show my logic. If you feel I made any mistakes, please share in the comments!\n\nSpecifically, I create CV split based on:\n- Scaffold split (edit - removed in V2 in favor of bb split and external data, credit to @roberthatch)\n- Building block split\n- Random split\n- External data from @chemdatafarmer\n\nThis is based on discussion shared by host [here](https://www.kaggle.com/competitions/leash-BELKA/discussion/491362#2737103).\n\n## DATASET WITH CV SPLIT: [LINK](https://www.kaggle.com/datasets/thedrcat/belka-cv-split) 👈👍🙏\n\nCredits: \n- https://www.kaggle.com/datasets/shlomoron/belka-shrunken-train-set\n- https://www.kaggle.com/code/shlomoron/belka-shrinking-the-dataset\n- https://www.kaggle.com/competitions/leash-BELKA/discussion/493294\n- https://www.kaggle.com/code/chemdatafarmer/additional-seh-data","metadata":{}},{"cell_type":"markdown","source":"## 1. Create building block dict across all building blocks columns\n\nThis is necessary as building blocks from groups 2 & 3 overlap. ","metadata":{}},{"cell_type":"code","source":"import dask.dataframe as dd\n\ntrain = dd.read_parquet('data/train.parquet')","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train.head(1)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"bblocks = set(train.buildingblock1_smiles)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"len(bblocks)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"bblocks2 = set(train.buildingblock2_smiles)\nlen(bblocks2)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"bblocks3 = set(train.buildingblock3_smiles)\nlen(bblocks3)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"shared_blocks = bblocks & (bblocks3 | bblocks2)\nlen(shared_blocks)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"len(bblocks3 & bblocks2)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"bblocksall = set(bblocks | bblocks2 | bblocks3)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"len(bblocksall), 271+693+872","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import numpy as np\nbbdict = {x:i for i,x in enumerate(bblocksall)}","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import pickle\npickle.dump(bbdict, open('bbdict.pickle', 'bw'))","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Let's do the same for test\ntest = dd.read_parquet('data/test.parquet')\nbblocks = set(test.buildingblock1_smiles)\nbblocks2 = set(test.buildingblock2_smiles)\nbblocks3 = set(test.buildingblock3_smiles)\nbblocksall = set(bblocks | bblocks2 | bblocks3)\nbbdicttest = {x:i for i,x in enumerate(bblocksall)}\npickle.dump(bbdicttest, open('bbdict-test.pickle', 'bw'))","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## 2. Loading data\n\nCredit for most of the dataloading functions to Greysnow's [notebook](https://www.kaggle.com/code/shlomoron/belka-shrinking-the-dataset). ","metadata":{}},{"cell_type":"code","source":"dataset_path = 'data/train.parquet'\nimport pandas as pd","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def get_encoded(dataset_path, col, BBs_dict):\n    BBs = pd.read_parquet(dataset_path, engine = 'pyarrow', columns=[col])\n    BBs = BBs[col].to_numpy()\n    BBs_reshaped = np.reshape(BBs, [-1, 3])\n    BBs = BBs_reshaped[:, 0]\n    encoded_BBs = [BBs_dict[x] for x in BBs]\n    encoded_BBs = np.asarray(encoded_BBs, dtype = np.int16)\n    return encoded_BBs\n\nencoded_BBs_1 = get_encoded(dataset_path, 'buildingblock1_smiles', bbdict)\nencoded_BBs_2 = get_encoded(dataset_path, 'buildingblock2_smiles', bbdict)\nencoded_BBs_3 = get_encoded(dataset_path, 'buildingblock3_smiles', bbdict)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def get_molecule_smiles(dataset_path):\n    molecule_smiles = pd.read_parquet(dataset_path, engine = 'pyarrow', columns=['molecule_smiles'])\n    molecule_smiles = molecule_smiles.molecule_smiles.to_numpy()\n    molecule_smiles = np.reshape(molecule_smiles, [-1, 3])\n    if np.mean(molecule_smiles[:, 0] == molecule_smiles[:, 1]) != 1:\n        print('ERROR')\n    if np.mean(molecule_smiles[:, 0] == molecule_smiles[:, 2]) != 1:\n        print('ERROR')\n    return molecule_smiles[:, 0]\n\nmolecule_smiles = get_molecule_smiles(dataset_path)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def get_binds(dataset_path):\n    binds =  pd.read_parquet(dataset_path, engine = 'pyarrow', columns=['binds']).binds.to_numpy()\n    return np.reshape(binds.astype('byte'), [-1, 3])\n\nbinds = get_binds(dataset_path)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"data = {'buildingblock1_smiles':encoded_BBs_1, 'buildingblock2_smiles':encoded_BBs_2, 'buildingblock3_smiles':encoded_BBs_3,\n        'molecule_smiles':molecule_smiles, 'binds_BRD4':binds[:, 0], 'binds_HSA':binds[:, 1], 'binds_sEH':binds[:, 2]}\ndf = pd.DataFrame(data=data)\ndf.head(2)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"pd.set_option('display.float_format', lambda x: '%.3f' % x)\ndf[['binds_BRD4','binds_HSA','binds_sEH']].describe()","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Now do the same for test\ntest_path = 'data/test.parquet'\n\nmolecule_smiles = pd.read_parquet(test_path, engine = 'pyarrow', columns=['molecule_smiles']).molecule_smiles.to_numpy()\nprotein_name = pd.read_parquet(test_path, engine = 'pyarrow', columns=['protein_name']).protein_name.to_numpy()\nfirst_unique_molecule_smiles_indices = []\nmolecule_smiles_unique = {}\nis_BRD4 = {}\nis_HSA = {}\nis_sEH = {}\nfor i,x in enumerate(molecule_smiles):\n    if x not in molecule_smiles_unique:\n        molecule_smiles_unique[x] = [i]\n        first_unique_molecule_smiles_indices.append(i)\n        is_BRD4[x] = False\n        is_HSA[x] = False\n        is_sEH[x] = False\n        if protein_name[i] == 'BRD4':\n            is_BRD4[x] = True\n        if protein_name[i] == 'HSA':\n            is_HSA[x] = True\n        if protein_name[i] == 'sEH':\n            is_sEH[x] = True\n    else:\n        molecule_smiles_unique[x].append(i)\n        if protein_name[i] == 'BRD4':\n            is_BRD4[x] = True\n        if protein_name[i] == 'HSA':\n            is_HSA[x] = True\n        if protein_name[i] == 'sEH':\n            is_sEH[x] = True\nfirst_unique_molecule_smiles_indices = np.asarray(first_unique_molecule_smiles_indices)\nprint(len(is_BRD4))\nprint(np.sum([is_BRD4[x] for x in is_BRD4]))\nprint(np.sum([is_HSA[x] for x in is_HSA]))\nprint(np.sum([is_sEH[x] for x in is_sEH]))","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"molecule_smiles_unique_arr = molecule_smiles[first_unique_molecule_smiles_indices]\nprint(len(np.unique(molecule_smiles_unique_arr)) == len(molecule_smiles_unique_arr))\n\nis_BRD4_arr = np.asarray([is_BRD4[x] for x in molecule_smiles_unique])\nis_HSA_arr = np.asarray([is_HSA[x] for x in molecule_smiles_unique])\nis_sEH_arr = np.asarray([is_sEH[x] for x in molecule_smiles_unique])\n\nprint(np.sum(is_BRD4_arr))\nprint(np.sum(is_HSA_arr))\nprint(np.sum(is_sEH_arr))","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def get_encoded_test(dataset_path, col, BBs_dict):\n    BBs = pd.read_parquet(dataset_path, engine = 'pyarrow', columns=[col])\n    BBs = BBs[col].to_numpy()\n    BBs = BBs[first_unique_molecule_smiles_indices]\n    encoded_BBs = [BBs_dict[x] for x in BBs]\n    encoded_BBs = np.asarray(encoded_BBs, dtype = np.int16)\n    return encoded_BBs\n\nencoded_BBs_1_test = get_encoded_test(test_path, 'buildingblock1_smiles', bbdicttest)\nencoded_BBs_2_test = get_encoded_test(test_path, 'buildingblock2_smiles', bbdicttest)\nencoded_BBs_3_test = get_encoded_test(test_path, 'buildingblock3_smiles', bbdicttest)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"testdata = {'buildingblock1_smiles':encoded_BBs_1_test, 'buildingblock2_smiles':encoded_BBs_2_test,\n        'buildingblock3_smiles':encoded_BBs_3_test,'molecule_smiles':molecule_smiles_unique_arr,\n        'is_BRD4':is_BRD4_arr, 'is_HSA':is_HSA_arr, 'is_sEH':is_sEH_arr}\ntestdf = pd.DataFrame(data=testdata)\ntestdf.to_parquet('test.parquet', index=False)\ntestdf.head(1)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## 3. Building block split\n\nLet's take 17 bbs from block1 and 34 from (block2 & block3) and 2 from (block 3 & !block2). \n\nCredit: https://www.kaggle.com/competitions/leash-BELKA/discussion/496576","metadata":{}},{"cell_type":"code","source":"import random\nrandom.seed(42)\nbb1ids = df.buildingblock1_smiles.unique().tolist()\nbb2ids = df.buildingblock2_smiles.unique().tolist()\nbb3ids = df.buildingblock3_smiles.unique().tolist()\n\ngroup2 = list(set(bb2ids) & set(bb3ids))\ngroup3 = list(set(bb3ids) - set(bb2ids))\n\nbbs1 = random.sample(bb1ids, 17)\nbbs2 = random.sample(group2, 34)\nbbs3 = random.sample(group3, 2)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df['bb_test_noshare'] = False","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df['bb_test_noshare'].iloc[df.buildingblock1_smiles.isin(bbs1) & df.buildingblock2_smiles.isin(bbs2) & df.buildingblock3_smiles.isin(bbs2+bbs3)] = True","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df['bb_test_noshare'].value_counts()","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df['bb_test_mixshare'] = False","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df['bb_test_mixshare'].iloc[df.buildingblock1_smiles.isin(bbs1) | \\\ndf.buildingblock2_smiles.isin(bbs2) | df.buildingblock3_smiles.isin(bbs2+bbs3)] \\\n= True","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df['bb_test_mixshare'].value_counts()","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"test1 = df[df['bb_test_noshare']==True]","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"test1[['binds_BRD4','binds_HSA','binds_sEH']].describe()","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"test1[['binds_BRD4','binds_HSA','binds_sEH']].sum()","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## 4. Random split\n\nThis one should be easy!","metadata":{}},{"cell_type":"code","source":"df['random_test'] = False","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import random\nrandom.seed(42)\nrandom_indices = random.sample(range(len(df)), int(len(df)*0.003))","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df['random_test'].iloc[random_indices] = True","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"test3 = df[df['random_test']==True]","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"test3[['binds_BRD4','binds_HSA','binds_sEH']].describe()","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"len(test3)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## 6. Let's put all tests together\n\nI'll create another `full_test` column to highlight all data points in test. I'd recommend using both combined and granular metrics for each split though!","metadata":{}},{"cell_type":"code","source":"df[['bb_test_noshare', 'bb_test_mixshare', 'random_test']].mean()","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df['full_test'] = df['bb_test_noshare'] + df['bb_test_mixshare'] + df['random_test']\ndf['full_test'] = df['full_test'] > 0","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df['full_test'].value_counts()","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df['full_test'].mean()","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df[df.full_test == True][['binds_BRD4','binds_HSA','binds_sEH']].describe()","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df[df.full_test == True][['binds_BRD4','binds_HSA','binds_sEH']].sum()","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df.head(1).T","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df.to_parquet('train_folds.parquet', index=False)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# End\n\nDataset: https://www.kaggle.com/datasets/thedrcat/belka-cv-split\n\nIf you found this useful, I'll appreciate feedback! ❤️🙏❤️","metadata":{}},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}