{"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":30684,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"![BELKA](https://i.imgur.com/NXE4Cgo.jpeg)\n\nIf you are wondering what's going on, just open Google Translate, translate 'squirrel' from English to Russian and listen how it sounds. Also, 'белки' is the plural for both 'proteins' and 'squirrels', the only difference is what syllable is stressed. This is my favourite source of puns and it's impossible to resist and not to add squirrel picture here. Especially this one - I've found the original without text in OpenImages dataset under Alpaca category several years ago. That's what typically happens to alpaca when they meet squirrels. Or not.","metadata":{}},{"cell_type":"markdown","source":"# Table of contents\n\n* [Train data 'compression'](#split_data)\n * [Code](#split_data_code)\n * [Results](#split_data_results)\n* [Unique molecules in building blocks 1-3](#unique_molecules)\n * [Code](#unique_molecules_code)\n * [Intersections between unique SMILES appearing in blocks 1-3](#unique_molecules_intersections)\n * [Collect unique SMILES strings appearing in building blocks](#unique_molecules_collect_smiles)\n* [Assessment of building blocks composition](#composition_plots)\n * [Fixed building block, many binding molecules](#composition_plots_1block)\n * [Pairs of building blocks: composition in binding molecules](#composition_plots_pairs)\n * [Building blocks composition in non-binding molecules](#composition_plots_pairs_not_binding)\n* [Filter negative data once more](#filter_negative_data)\n* [Compute fingerprints for all unique SMILES strings appearing in building blocks (unfinished)](#compute_fingerprints)\n* [Fragment analysis](#fragment_analysis)\n * [Fragment extraction: example 1](#fragment_extraction_example1)\n * [Fragment extraction: example 2](#fragment_extraction_example2)\n * [Fragment extraction and match: code](#fragment_extraction_code)\n * [Fragment extraction and match: partial processing](#fragment_extraction_processing)\n* [Conclusions](#conclusions)","metadata":{}},{"cell_type":"markdown","source":"If you'll check out public EDA notebooks, i.e., these two:\n\nhttps://www.kaggle.com/code/seshurajup/eda-smiles\n\nhttps://www.kaggle.com/code/taichiuemura/belka-eda\n\n-you will see that there are only 3 protein names in both train and test datasets. \n\nIf you'll look at the data closely, you might notice that typically each small molecule appear more than once in our train data: there are 3 consecutive rows with with the same molecule and one of the given proteins.\n\n# Train data 'compression'\n\nLet's 'compress' the provided train.parquet file by splitting it to several parts and removing duplicates!\n\n## Code","metadata":{}},{"cell_type":"code","source":"!pip install supervenn\n!pip install upsetplot\n!pip install mapchiral\n!pip install rdkit","metadata":{"_kg_hide-output":true,"_kg_hide-input":true,"execution":{"iopub.status.busy":"2024-04-28T14:04:43.998967Z","iopub.execute_input":"2024-04-28T14:04:43.999406Z","iopub.status.idle":"2024-04-28T14:06:03.193664Z","shell.execute_reply.started":"2024-04-28T14:04:43.999373Z","shell.execute_reply":"2024-04-28T14:06:03.192206Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# This Python 3 environment comes with many helpful analytics libraries installed\n# It is defined by the kaggle/python Docker image: https://github.com/kaggle/docker-python\n# For example, here's several helpful packages to load\n\nimport numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\n\nfrom tqdm import tqdm\nfrom pathlib import Path\nimport itertools\nfrom collections import defaultdict\n\nimport pyarrow as pa\nimport pyarrow.parquet as pq\nfrom supervenn import supervenn\n\nfrom rdkit.Chem import rdFMCS\nfrom rdkit import Chem\nfrom mapchiral.mapchiral import encode, jaccard_similarity\n\nimport matplotlib.pyplot as plt\nimport seaborn as sns\nimport os\nimport warnings\nwarnings.filterwarnings('ignore')\n\nfor dirname, _, filenames in os.walk('../input'):\n for filename in filenames[:3]:\n print(os.path.join(dirname, filename))\n\n# You can write up to 20GB to the current directory (/kaggle/working/) that gets preserved as output when you create a version using \"Save & Run All\" \n# You can also write temporary files to /kaggle/temp/, but they won't be saved outside of the current session","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","_kg_hide-output":true,"_kg_hide-input":true,"execution":{"iopub.status.busy":"2024-04-28T14:06:03.196036Z","iopub.execute_input":"2024-04-28T14:06:03.196492Z","iopub.status.idle":"2024-04-28T14:06:06.133033Z","shell.execute_reply.started":"2024-04-28T14:06:03.196445Z","shell.execute_reply":"2024-04-28T14:06:06.131682Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"SAVEDIR = Path(\"train\")\nSAVEDIR.mkdir(exist_ok=True)\n\nDRAFT_MODE = False # this should be set to False to enable long computations\n ","metadata":{"execution":{"iopub.status.busy":"2024-04-28T14:06:06.134688Z","iopub.execute_input":"2024-04-28T14:06:06.135376Z","iopub.status.idle":"2024-04-28T14:06:06.141919Z","shell.execute_reply.started":"2024-04-28T14:06:06.135333Z","shell.execute_reply":"2024-04-28T14:06:06.140676Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"protein_names = [\"BRD4\", \"HSA\", \"sEH\"]\nsmiles_col_names = [\n \"buildingblock1_smiles\", \n \"buildingblock2_smiles\", \n \"buildingblock3_smiles\", \n \"molecule_smiles\"\n]","metadata":{"execution":{"iopub.status.busy":"2024-04-28T14:06:06.144791Z","iopub.execute_input":"2024-04-28T14:06:06.145252Z","iopub.status.idle":"2024-04-28T14:06:06.160061Z","shell.execute_reply.started":"2024-04-28T14:06:06.145215Z","shell.execute_reply":"2024-04-28T14:06:06.158887Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We'll need to define the schemas for each of the parquet files we want to create:","metadata":{}},{"cell_type":"code","source":"schema_single_target = pa.schema([\n (col_name, pa.string())\n for col_name in smiles_col_names\n])\n\nschema_mixed_target = pa.schema([\n (col_name, pa.string())\n for col_name in smiles_col_names\n] + [(prot_name, pa.int8()) for prot_name in protein_names])\n\nschema_unprocessed = pa.schema([\n (col_name, pa.string())\n for col_name in smiles_col_names\n] + [\n ('protein_name', pa.string()),\n ('binds', pa.int8())\n])\n\nKEYS_INFO = [\n (\"all0\", schema_single_target), \n (\"all1\", schema_single_target), \n (\"all_mixed\", schema_mixed_target),\n ('unprocessed', schema_unprocessed)\n]","metadata":{"execution":{"iopub.status.busy":"2024-04-28T14:06:06.161557Z","iopub.execute_input":"2024-04-28T14:06:06.162017Z","iopub.status.idle":"2024-04-28T14:06:06.173717Z","shell.execute_reply.started":"2024-04-28T14:06:06.161972Z","shell.execute_reply":"2024-04-28T14:06:06.172510Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Next I will define helper class to read data from `train.parquet` in batches, split it to parts, append them to several different files","metadata":{}},{"cell_type":"code","source":"class PqDataSaverWrapper:\n def __init__(self, savedir=\".\", keys_info=KEYS_INFO, prefix=\"train\"):\n self.file_handlers = None\n if isinstance(savedir, str):\n savedir = Path(savedir)\n self.savedir = savedir\n self.keys_info = keys_info\n self.prefix = prefix\n\n def add_data(self, df_dict):\n if self.file_handlers is None:\n raise Exception(\"Attempt to call 'add_data' outside of with-block\")\n for key in df_dict:\n if not key in self.file_handlers:\n raise Exception(f\"Key should be equal to one of the provided in __init__ method, got '{key}' instead\")\n if df_dict[key].shape[0] > 0:\n pq_batch = pa.RecordBatch.from_pandas(df_dict[key])\n self.file_handlers[key].write(pq_batch)\n\n def __enter__(self):\n self.file_handlers = {}\n for key, schema in self.keys_info:\n filename = self.savedir / f\"{self.prefix}_{key}.parquet\"\n if filename.exists():\n raise Exception(f\"File '{filename}' for key '{key}' already exists, remove it if you want to regenerate everything\")\n self.file_handlers[key] = pq.ParquetWriter(filename.as_posix(), schema)\n return self\n\n def __exit__(self, exc_type, exc_value, traceback):\n for pq_writer in self.file_handlers.values():\n pq_writer.close()\n self.file_handlers = None\n","metadata":{"execution":{"iopub.status.busy":"2024-04-28T14:06:06.174920Z","iopub.execute_input":"2024-04-28T14:06:06.175369Z","iopub.status.idle":"2024-04-28T14:06:06.194435Z","shell.execute_reply.started":"2024-04-28T14:06:06.175329Z","shell.execute_reply":"2024-04-28T14:06:06.192986Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Now let's start processing `train.parquet`","metadata":{}},{"cell_type":"code","source":"train_pq = pq.ParquetFile(\"../input/leash-BELKA/train.parquet\")\n\npartial_df = None\nadditional_records = []\nwith PqDataSaverWrapper(savedir=SAVEDIR, keys_info=KEYS_INFO) as pq_wrapper:\n for i in tqdm(np.arange(train_pq.num_row_groups), total=train_pq.num_row_groups):\n if i > 10 and DRAFT_MODE:\n break\n batch = train_pq.read_row_group(i)\n batch_df = batch.to_pandas()\n if partial_df is not None and partial_df.shape[0] > 0:\n batch_df = pd.concat([partial_df, batch_df], ignore_index=True)\n partial_df = None\n batch_df = batch_df.pivot_table(\n columns=\"protein_name\", \n values=\"binds\", \n index=smiles_col_names).reset_index(drop=False)\n batch_df.columns.name = None\n missing_protein_data_ids = batch_df[protein_names].isnull().any(axis=1)\n\n # filename = SAVEDIR / \"all_present\"/ f\"part_{i}.csv\"\n new_data = {}\n if (~missing_protein_data_ids).any():\n all0_ids = batch_df[protein_names].max(1) == 0\n all1_ids = batch_df[protein_names].min(1) == 1\n new_data['all0'] = batch_df[(~missing_protein_data_ids) & all0_ids][smiles_col_names].reset_index(drop=True)\n new_data['all1'] = batch_df[(~missing_protein_data_ids) & all1_ids][smiles_col_names].reset_index(drop=True)\n all_mixed_df = batch_df[(~missing_protein_data_ids) & (~all1_ids) & (~all0_ids)].reset_index(drop=True)\n all_mixed_df[protein_names] = all_mixed_df[protein_names].astype(np.int8)\n new_data['all_mixed'] = all_mixed_df\n\n if missing_protein_data_ids.any():\n partial_df = batch_df[missing_protein_data_ids].melt(\n id_vars=smiles_col_names,\n value_vars=protein_names,\n value_name=\"binds\",\n var_name=\"protein_name\"\n )\n partial_df = partial_df[~partial_df.binds.isnull()].reset_index(drop=True)\n partial_df['binds'] = partial_df['binds'].astype(np.int8)\n pq_wrapper.add_data(new_data)\n\n if partial_df is not None and partial_df.shape[0] > 0:\n pq_wrapper.add_data({\n 'unprocessed': partial_df\n })\n \n","metadata":{"execution":{"iopub.status.busy":"2024-04-28T14:06:06.195834Z","iopub.execute_input":"2024-04-28T14:06:06.196831Z","iopub.status.idle":"2024-04-28T14:06:46.333593Z","shell.execute_reply.started":"2024-04-28T14:06:06.196798Z","shell.execute_reply":"2024-04-28T14:06:46.332270Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Results\n\nWe've processed `train.parquet` file by splitting it to several parts and removing duplicates. As a result, we have the separate files for the following:\n\n- small molecules which don't bind to any of the given proteins (I'll call this file `train_all0.parquet`);\n- small molecules which bind to all the given proteins (`train_all1.parquet`);\n- small molecules which bind to 1 or 2 of the given proteins (`train_all_mixed.parquet`);\n- unprocessed data (`train_unprocessed.parquet`): if there are small molecules which don't have information for any of our 3 proteins, they are in this file","metadata":{}},{"cell_type":"code","source":"train_all1_pq = pq.ParquetFile((SAVEDIR/\"train_all1.parquet\").as_posix())\nprint(\"Number of small molecules which bind to all 3 proteins:\",\n train_all1_pq.scan_contents())","metadata":{"execution":{"iopub.status.busy":"2024-04-28T14:08:23.402250Z","iopub.execute_input":"2024-04-28T14:08:23.402819Z","iopub.status.idle":"2024-04-28T14:08:23.411463Z","shell.execute_reply.started":"2024-04-28T14:08:23.402772Z","shell.execute_reply":"2024-04-28T14:08:23.410281Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_all0_pq = pq.ParquetFile((SAVEDIR/\"train_all0.parquet\").as_posix())\nprint(\"Number of small molecules which don't bind to any of 3 proteins:\",\n train_all0_pq.scan_contents())\n","metadata":{"execution":{"iopub.status.busy":"2024-04-28T14:08:25.811087Z","iopub.execute_input":"2024-04-28T14:08:25.811553Z","iopub.status.idle":"2024-04-28T14:08:26.088543Z","shell.execute_reply.started":"2024-04-28T14:08:25.811517Z","shell.execute_reply":"2024-04-28T14:08:26.087392Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_all_mixed_pq = pq.ParquetFile((SAVEDIR/\"train_all_mixed.parquet\").as_posix())\nprint(\"Number of small molecules which bind to 1 or 2 out of 3 proteins:\",\n train_all_mixed_pq.scan_contents())","metadata":{"execution":{"iopub.status.busy":"2024-04-28T14:08:26.761903Z","iopub.execute_input":"2024-04-28T14:08:26.764291Z","iopub.status.idle":"2024-04-28T14:08:26.777528Z","shell.execute_reply.started":"2024-04-28T14:08:26.764245Z","shell.execute_reply":"2024-04-28T14:08:26.776471Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df = pd.read_parquet(SAVEDIR / \"train_unprocessed.parquet\")\nprint(\"Number of unprocessed molecules:\", df.shape[0])\n# no rows here means that all the molecules in train dataset have information \n# on binding with all 3 proteins","metadata":{"execution":{"iopub.status.busy":"2024-04-28T14:08:27.099877Z","iopub.execute_input":"2024-04-28T14:08:27.100724Z","iopub.status.idle":"2024-04-28T14:08:27.172349Z","shell.execute_reply.started":"2024-04-28T14:08:27.100689Z","shell.execute_reply":"2024-04-28T14:08:27.171385Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The resulting files with positive samples (molecules binding to 1+ proteins - files `train_all1.parquet` and `train_all_mixed.parquet`) fit into memory and can be opened with pandas. However, the negative samples (`train_all0.parquet`) are the majority of all the available data and better represent 'infinite mini-universe'. To view it or filter the data later we might want to collect unique SMILES string from each of the building blocks.","metadata":{}},{"cell_type":"markdown","source":"# Unique molecules in building blocks 1-3\n## Code\n\nLet's extract unique SMILES strings from the `train_all0.parquet` file, building block fields.","metadata":{}},{"cell_type":"code","source":"UNIQUE_DIR = Path(\"unique_smiles\")\nUNIQUE_DIR.mkdir(exist_ok=True)","metadata":{"execution":{"iopub.status.busy":"2024-04-28T14:08:29.990119Z","iopub.execute_input":"2024-04-28T14:08:29.990545Z","iopub.status.idle":"2024-04-28T14:08:29.996145Z","shell.execute_reply.started":"2024-04-28T14:08:29.990513Z","shell.execute_reply":"2024-04-28T14:08:29.995011Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"smiles_collections = {} ","metadata":{"execution":{"iopub.status.busy":"2024-04-28T14:08:30.493795Z","iopub.execute_input":"2024-04-28T14:08:30.494201Z","iopub.status.idle":"2024-04-28T14:08:30.499243Z","shell.execute_reply.started":"2024-04-28T14:08:30.494171Z","shell.execute_reply":"2024-04-28T14:08:30.498012Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_all0_pq = pq.ParquetFile((SAVEDIR/\"train_all0.parquet\").as_posix())\n\nfor k in np.arange(3):\n smiles_list = []\n col_name = f\"buildingblock{k+1}_smiles\"\n print(\"Processing\", col_name, \"...\")\n for i in tqdm(np.arange(train_all0_pq.num_row_groups), total=train_all0_pq.num_row_groups):\n batch = train_all0_pq.read_row_group(i, columns=[col_name])\n batch_df = batch.to_pandas()\n batch_df.drop_duplicates(inplace=True)\n smiles_list.extend(batch_df.to_dict('records'))\n df = pd.DataFrame(smiles_list).drop_duplicates()[col_name]\n print(f\"Result: {df.shape[0]} unique SMILES strings\\n\")\n smiles_collections[(\"all0\", col_name)] = df\n df.to_csv(\n UNIQUE_DIR / f\"train_all0_{col_name}.csv\", index=None, header=None\n )","metadata":{"execution":{"iopub.status.busy":"2024-04-28T14:08:31.009362Z","iopub.execute_input":"2024-04-28T14:08:31.009748Z","iopub.status.idle":"2024-04-28T14:08:32.333169Z","shell.execute_reply.started":"2024-04-28T14:08:31.009721Z","shell.execute_reply":"2024-04-28T14:08:32.332310Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_all1_df = pd.read_parquet(SAVEDIR/\"train_all1.parquet\")\nprint(\"Small molecules binding to all 3 targets\")\nprint(\"Total:\", train_all1_df.shape[0], \"small molecules\")\nfor k in np.arange(3):\n smiles_list = []\n col_name = f\"buildingblock{k+1}_smiles\"\n df = train_all1_df[col_name].drop_duplicates()\n print(f\"{col_name}: {df.shape[0]} unique SMILES strings\")\n smiles_collections[(\"all1\", col_name)] = df\n df.to_csv(\n UNIQUE_DIR / f\"train_all1_{col_name}.csv\", index=None, header=None\n )","metadata":{"execution":{"iopub.status.busy":"2024-04-28T14:08:32.334806Z","iopub.execute_input":"2024-04-28T14:08:32.335966Z","iopub.status.idle":"2024-04-28T14:08:32.351451Z","shell.execute_reply.started":"2024-04-28T14:08:32.335908Z","shell.execute_reply":"2024-04-28T14:08:32.350361Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_all_mixed_df = pd.read_parquet(SAVEDIR / \"train_all_mixed.parquet\")\nprint(\"Small molecules binding to 1 or 2 of 3 targets\")\nprint(\"Total:\", train_all_mixed_df.shape[0], \"small molecules\")\nfor k in np.arange(3):\n smiles_list = []\n col_name = f\"buildingblock{k+1}_smiles\"\n df = train_all_mixed_df[col_name].drop_duplicates()\n print(f\"{col_name}: {df.shape[0]} unique SMILES strings\")\n smiles_collections[(\"all_mixed\", col_name)] = df\n df.to_csv(\n UNIQUE_DIR / f\"train_all_mixed_{col_name}.csv\", index=None, header=None\n )\n \n for protein_name in protein_names:\n df = train_all_mixed_df.loc[train_all_mixed_df[protein_name] > 0, col_name].drop_duplicates()\n print(f\" - including {df.shape[0]} SMILES blocks in molecules which bind to {protein_name}\")\n smiles_collections[(protein_name, col_name)] = df","metadata":{"execution":{"iopub.status.busy":"2024-04-28T14:08:32.352751Z","iopub.execute_input":"2024-04-28T14:08:32.353113Z","iopub.status.idle":"2024-04-28T14:08:32.401882Z","shell.execute_reply.started":"2024-04-28T14:08:32.353082Z","shell.execute_reply":"2024-04-28T14:08:32.400773Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# test_df = pd.read_parquet(\"../input/leash-BELKA/test.parquet\")\n# test_df_collected = test_df.pivot_table(\n# columns=\"protein_name\", \n# values=\"id\", \n# index=smiles_col_names).reset_index(drop=False)\n# test_df_collected = test_df_collected[protein_names].isnull().reset_index().groupby(protein_names).count()\n# import upsetplot\n# # test_df_collected.values\n# upsetplot.plot(\n# test_df_collected, \n# sum_over='index', #show_counts=True,\n# #orientation=\"vertical\"\n# ); #, subset_size='count')","metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"execution":{"iopub.status.busy":"2024-04-28T14:08:32.756766Z","iopub.execute_input":"2024-04-28T14:08:32.757179Z","iopub.status.idle":"2024-04-28T14:08:32.762330Z","shell.execute_reply.started":"2024-04-28T14:08:32.757145Z","shell.execute_reply":"2024-04-28T14:08:32.761127Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Intersections between unique SMILES appearing in blocks 1-3\n\nYes, every molecule might have multiple SMILES representations (it's isomeric SMILES in our case, which means more fun). Let's ignore this for now and check if the provided representations of our building blocks duplicate each other.","metadata":{}},{"cell_type":"code","source":"# {k: v.values for k, v in smiles_collections.items()}\nblock_column_names = [\"buildingblock1_smiles\", \"buildingblock2_smiles\", \"buildingblock3_smiles\"]","metadata":{"execution":{"iopub.status.busy":"2024-04-28T14:08:34.296122Z","iopub.execute_input":"2024-04-28T14:08:34.296897Z","iopub.status.idle":"2024-04-28T14:08:34.303117Z","shell.execute_reply.started":"2024-04-28T14:08:34.296855Z","shell.execute_reply":"2024-04-28T14:08:34.301580Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"all0_smiles_subsets = [set(list(smiles_collections[(\"all0\", col)].values)) for col in block_column_names]\nsupervenn(all0_smiles_subsets, block_column_names, side_plots=True);\n","metadata":{"execution":{"iopub.status.busy":"2024-04-28T14:08:34.994927Z","iopub.execute_input":"2024-04-28T14:08:34.996095Z","iopub.status.idle":"2024-04-28T14:08:35.556206Z","shell.execute_reply.started":"2024-04-28T14:08:34.996046Z","shell.execute_reply":"2024-04-28T14:08:35.555004Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"\nall1_smiles_subsets = [\n set(list(smiles_collections[(\"all1\", col)].values)) for col in block_column_names]\nif np.max([len(s) for s in all1_smiles_subsets]) > 0:\n supervenn(all1_smiles_subsets, block_column_names, side_plots=True);\n","metadata":{"execution":{"iopub.status.busy":"2024-04-28T14:20:13.961795Z","iopub.execute_input":"2024-04-28T14:20:13.962748Z","iopub.status.idle":"2024-04-28T14:20:13.970407Z","shell.execute_reply.started":"2024-04-28T14:20:13.962710Z","shell.execute_reply":"2024-04-28T14:20:13.969230Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"all_mixed_smiles_subsets = [set(list(smiles_collections[(\"all_mixed\", col)].values)) for col in block_column_names]\nsupervenn(all_mixed_smiles_subsets, block_column_names, side_plots=True);","metadata":{"execution":{"iopub.status.busy":"2024-04-28T14:20:16.166195Z","iopub.execute_input":"2024-04-28T14:20:16.166647Z","iopub.status.idle":"2024-04-28T14:20:16.716160Z","shell.execute_reply.started":"2024-04-28T14:20:16.166612Z","shell.execute_reply":"2024-04-28T14:20:16.714936Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"all_smiles_subsets = [\n set([\n x\n for subset_name in [\"all0\", \"all1\", \"all_mixed\"]\n for x in smiles_collections[(subset_name, col)].values\n ])\n for col in block_column_names \n]\nsupervenn(all_smiles_subsets, block_column_names, side_plots=True);","metadata":{"execution":{"iopub.status.busy":"2024-04-28T14:08:36.597390Z","iopub.status.idle":"2024-04-28T14:08:36.598232Z","shell.execute_reply.started":"2024-04-28T14:08:36.597929Z","shell.execute_reply":"2024-04-28T14:08:36.597971Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Collect unique SMILES strings appearing in building blocks","metadata":{}},{"cell_type":"code","source":"unique_smiles_train = np.concatenate([\n smiles_collections[(subset_name, col)].values\n for col in block_column_names for subset_name in [\"all0\", \"all1\", \"all_mixed\"]\n])\n\nunique_smiles_train = np.unique(unique_smiles_train)\nprint(\n \"Number of unique SMILES appearing in building blocks in train dataset:\", \n unique_smiles_train.shape[0]\n)","metadata":{"execution":{"iopub.status.busy":"2024-04-28T14:08:37.586824Z","iopub.execute_input":"2024-04-28T14:08:37.587579Z","iopub.status.idle":"2024-04-28T14:08:37.599880Z","shell.execute_reply.started":"2024-04-28T14:08:37.587532Z","shell.execute_reply":"2024-04-28T14:08:37.598700Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"for protein_name in protein_names:\n smiles_list = []\n for k in np.arange(3):\n col_name = f\"buildingblock{k+1}_smiles\"\n smiles_list.append(smiles_collections[(protein_name, col_name)].values)\n smiles_list = np.unique(np.concatenate(smiles_list))\n print(protein_name, smiles_list.shape[0])","metadata":{"execution":{"iopub.status.busy":"2024-04-28T14:08:38.068468Z","iopub.execute_input":"2024-04-28T14:08:38.069310Z","iopub.status.idle":"2024-04-28T14:08:38.082486Z","shell.execute_reply.started":"2024-04-28T14:08:38.069264Z","shell.execute_reply":"2024-04-28T14:08:38.081223Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"test_df = pd.read_parquet(\"../input/leash-BELKA/test.parquet\")\nunique_smiles_test = np.concatenate([test_df[f\"buildingblock{i+1}_smiles\"].unique() for i in np.arange(3)])\nunique_smiles_test = np.unique(unique_smiles_test)\nprint(\"Number of unique SMILES strings in blocks 1-3, test dataset:\", \n unique_smiles_test.shape[0])","metadata":{"execution":{"iopub.status.busy":"2024-04-28T14:08:38.675881Z","iopub.execute_input":"2024-04-28T14:08:38.676400Z","iopub.status.idle":"2024-04-28T14:08:40.712392Z","shell.execute_reply.started":"2024-04-28T14:08:38.676363Z","shell.execute_reply":"2024-04-28T14:08:40.711139Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"unique_smiles = np.union1d(unique_smiles_test, unique_smiles_train)\nprint(\"Number of unique SMILES strings in blocks 1-3, train and test datasets:\", unique_smiles.shape[0])","metadata":{"execution":{"iopub.status.busy":"2024-04-28T14:08:40.714628Z","iopub.execute_input":"2024-04-28T14:08:40.715091Z","iopub.status.idle":"2024-04-28T14:08:40.730089Z","shell.execute_reply.started":"2024-04-28T14:08:40.715051Z","shell.execute_reply":"2024-04-28T14:08:40.728686Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"pd.DataFrame(unique_smiles).to_csv(UNIQUE_DIR/\"all_smiles.csv\", index=None, header=False)","metadata":{"execution":{"iopub.status.busy":"2024-04-28T14:08:40.731866Z","iopub.execute_input":"2024-04-28T14:08:40.733093Z","iopub.status.idle":"2024-04-28T14:08:40.746748Z","shell.execute_reply.started":"2024-04-28T14:08:40.733049Z","shell.execute_reply":"2024-04-28T14:08:40.745481Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Assessment of building blocks composition\n\nHaving list of all unique smiles appearing in the dataset, we can use it to build index and to see how often certain building block appear together. ","metadata":{}},{"cell_type":"code","source":"smiles2index = {smiles: i for i, smiles in enumerate(unique_smiles)}\n\nprotein_names = [\"BRD4\", \"HSA\", \"sEH\"]\n\nbuilding_blocks_columns = [\n \"buildingblock1_smiles\", \n \"buildingblock2_smiles\", \n \"buildingblock3_smiles\"\n]\nindex_columns = [\"index1\", \"index2\", \"index3\"]\npair_bb_columns = [list(pair) for pair in itertools.combinations(index_columns, 2)]\npair_bb_columns","metadata":{"execution":{"iopub.status.busy":"2024-04-28T14:08:40.798610Z","iopub.execute_input":"2024-04-28T14:08:40.799046Z","iopub.status.idle":"2024-04-28T14:08:40.811174Z","shell.execute_reply.started":"2024-04-28T14:08:40.799010Z","shell.execute_reply":"2024-04-28T14:08:40.809441Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_mixed_df = pd.read_parquet(SAVEDIR / \"train_all_mixed.parquet\")\ntrain_all1_df = pd.read_parquet(SAVEDIR / \"train_all1.parquet\")\n# add column with corresponding index for smiles strings from building blocks 1-3\nfor i, column in enumerate(building_blocks_columns):\n train_all1_df[f\"index{i+1}\"] = train_all1_df[column].apply(lambda x: smiles2index[x])\n train_mixed_df[f\"index{i+1}\"] = train_mixed_df[column].apply(lambda x: smiles2index[x])","metadata":{"execution":{"iopub.status.busy":"2024-04-28T14:08:41.252709Z","iopub.execute_input":"2024-04-28T14:08:41.253102Z","iopub.status.idle":"2024-04-28T14:08:41.342788Z","shell.execute_reply.started":"2024-04-28T14:08:41.253073Z","shell.execute_reply":"2024-04-28T14:08:41.341579Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Fixed building block, many binding molecules\n\nLet's check how many variants of binding molecules are there for each (fixed) building block.","metadata":{}},{"cell_type":"code","source":"all1_blocks, all1_counts = np.unique(\n np.concatenate([\n train_all1_df[index_column].values \n for index_column in index_columns\n ]), \n return_counts=True\n)\n# for each of the proteins we will count the number of molecules separately\nsingle_blocks_positive_counter = {\n name: defaultdict(int) for name in protein_names\n}\n\nfor protein_name in tqdm(protein_names):\n ids = train_mixed_df[protein_name] == 1\n for index_column in index_columns:\n new_blocks, new_counts = np.unique(\n train_mixed_df.loc[ids, index_column].values, \n return_counts=True\n )\n for block, count in zip(new_blocks, new_counts):\n single_blocks_positive_counter[protein_name][block] += count\n\n for block, count in zip(all1_blocks, all1_counts):\n single_blocks_positive_counter[protein_name][block] += count\n","metadata":{"execution":{"iopub.status.busy":"2024-04-28T14:08:42.321828Z","iopub.execute_input":"2024-04-28T14:08:42.322279Z","iopub.status.idle":"2024-04-28T14:08:42.351148Z","shell.execute_reply.started":"2024-04-28T14:08:42.322245Z","shell.execute_reply":"2024-04-28T14:08:42.350040Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"for protein_name, data in single_blocks_positive_counter.items():\n fig = plt.figure(figsize=(10, 2))\n ax = plt.gca()\n sns.histplot(\n x=list(data.values()), \n element=\"step\", \n fill=False, binwidth=1,\n label=protein_name,\n ax=ax,\n common_norm=False\n )\n \n plt.xlabel(\"Number of small molecules sharing the same (single) building block\")\n plt.title(f\"Histogram for binding molecules in train dataset, {protein_name}\")\n # beautification: remove tick 0 and add tick 1\n ticks, labels = plt.xticks()\n ticks = [1] + [int(t) for t in ticks if t > 0]\n labels = [str(t) for t in ticks]\n plt.xticks(ticks=ticks, labels=labels)\n plt.show()\n","metadata":{"execution":{"iopub.status.busy":"2024-04-28T14:08:42.811684Z","iopub.execute_input":"2024-04-28T14:08:42.812564Z","iopub.status.idle":"2024-04-28T14:08:43.463583Z","shell.execute_reply.started":"2024-04-28T14:08:42.812521Z","shell.execute_reply":"2024-04-28T14:08:43.462333Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"block_indices = np.arange(len(unique_smiles))\nprint(\"Number of unique SMILES in building blocks which appear in binding molecules only once\")\nfor protein_name, data in single_blocks_positive_counter.items():\n values = np.asarray(list(data.values()))\n print(\n f\"- for protein {protein_name}:\", np.sum([values == 1])\n )","metadata":{"execution":{"iopub.status.busy":"2024-04-28T14:08:43.465789Z","iopub.execute_input":"2024-04-28T14:08:43.466237Z","iopub.status.idle":"2024-04-28T14:08:43.474480Z","shell.execute_reply.started":"2024-04-28T14:08:43.466196Z","shell.execute_reply":"2024-04-28T14:08:43.473326Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Pairs of building blocks: composition in binding molecules \n\nLet's fix each unique pair of molecules appearing in building blocks and see how many variants of the 3rd molecule are there. We might also want to take into account that order of building blocks in molecules is arbitrary, i.e. molecules with building blocks (A, B, C) and (D, B, A) might both appear in dataset and they both have (A, B) as their parts. Also, each molecule built with blocks (A, B, C) will be taken into account 3 times: as it consists of pairs (A, B), (B, C), (A, C). I assume that each pair of blocks (A, B) appear in sorted order, as it is the same as (B, A) pair.","metadata":{}},{"cell_type":"code","source":"all1_pairs, all1_counts = np.unique(\n np.concatenate([\n train_all1_df[list(pair)].apply(lambda row: tuple(sorted(row)), axis=1).values \n for pair in pair_bb_columns\n ]), \n return_counts=True\n)","metadata":{"execution":{"iopub.status.busy":"2024-04-28T14:08:44.756535Z","iopub.execute_input":"2024-04-28T14:08:44.756910Z","iopub.status.idle":"2024-04-28T14:08:44.766394Z","shell.execute_reply.started":"2024-04-28T14:08:44.756874Z","shell.execute_reply":"2024-04-28T14:08:44.765045Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# for each of the proteins we will count the number of molecules separately\npair_blocks_positive_counter = {\n name: defaultdict(int) for name in protein_names\n}\n\nfor protein_name in tqdm(protein_names):\n ids = train_mixed_df[protein_name] == 1\n for column_pair in pair_bb_columns:\n new_pairs, new_counts = np.unique(\n train_mixed_df.loc[ids, column_pair].apply(lambda row: tuple(sorted(row)), axis=1).values, \n return_counts=True\n )\n for pair, count in zip(new_pairs, new_counts):\n pair_blocks_positive_counter[protein_name][pair] += count\n\n for pair, count in zip(all1_pairs, all1_counts):\n pair_blocks_positive_counter[protein_name][pair] += count\n","metadata":{"execution":{"iopub.status.busy":"2024-04-28T14:08:45.255799Z","iopub.execute_input":"2024-04-28T14:08:45.256667Z","iopub.status.idle":"2024-04-28T14:08:46.382519Z","shell.execute_reply.started":"2024-04-28T14:08:45.256626Z","shell.execute_reply":"2024-04-28T14:08:46.381382Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig = plt.figure(figsize=(10, 10))\nax = plt.gca()\n# yticks = []\nfor protein_name, data in pair_blocks_positive_counter.items():\n # yticks.append(np.sum([xs == 1 for xs in data.values()]))\n # ones = np.sum([xs == 1 for xs in data.values()])\n # twos = np.sum([xs == 2 for xs in data.values()])\n # total = np.sum(list(data.values()))\n # print(f\"{protein_name}: {ones}, {twos}, {total} total\")\n sns.histplot(\n x=list(data.values()), element=\"step\", \n log_scale=(False, True), \n fill=False, binwidth=10,\n binrange=(1, 900),\n label=protein_name,\n ax=ax,\n common_norm=False\n # cumulative=True,\n # stat=\"percent\"\n )\n \n #plt.title(protein_name)\nplt.xlabel(\"Number of small molecules sharing the same pair of building blocks\")\nplt.title(\"Histogram for binding molecules in train dataset\")\nplt.axvline(x=1, ls='--')\nplt.xticks([1, 200, 400, 600, 800])\n#plt.yticks(yticks)\nplt.legend()\n\naxins = ax.inset_axes(\n [0.4, 0.5, 0.45, 0.45],\n xlim=(1, 10),\n xticks=np.arange(1, 10, 2))\n\nfor protein_name, data in pair_blocks_positive_counter.items():\n # yticks.append(np.sum([xs == 1 for xs in data.values()]))\n # ones = np.sum([xs == 1 for xs in data.values()])\n # twos = np.sum([xs == 2 for xs in data.values()])\n # total = len(list(data.values()))\n # print(f\"{protein_name}: {ones}, {twos}, {total} total\")\n sns.histplot(\n x=list(data.values()), element=\"step\", \n # log_scale=(False, True), \n fill=False, binwidth=1,\n binrange=(1, 10),\n label=protein_name,\n ax=axins,\n \n #cumulative=True\n )\n axins.set_title(\"closer look with non-logarithmic y-scale:\")\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2024-04-28T14:08:46.384893Z","iopub.execute_input":"2024-04-28T14:08:46.385664Z","iopub.status.idle":"2024-04-28T14:08:47.801430Z","shell.execute_reply.started":"2024-04-28T14:08:46.385621Z","shell.execute_reply":"2024-04-28T14:08:47.800301Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig = plt.figure(figsize=(5, 5))\nax = plt.gca()\n# yticks = []\nfor protein_name, data in pair_blocks_positive_counter.items():\n # yticks.append(np.sum([xs == 1 for xs in data.values()]))\n # ones = np.sum([xs == 1 for xs in data.values()])\n # twos = np.sum([xs < 5 for xs in data.values()])\n # total = len(list(data.values()))\n # print(f\"{protein_name}: {ones}, {twos}, {total} total\", ones/total, (twos)/total)\n sns.histplot(\n x=list(data.values()), element=\"step\", \n # log_scale=(False, True), \n fill=False, binwidth=1,\n binrange=(1, 900),\n label=protein_name,\n ax=ax,\n common_norm=False,\n cumulative=True,\n stat=\"percent\"\n )\nplt.xlim((1, 10))\nplt.legend()\nplt.axhline(y=80, ls='--')\nplt.xlabel(\"Number of small molecules sharing the same pair of building blocks\")\nplt.title(\"Histogram for binding molecules in train dataset: cumulative with percentage\")\n\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2024-04-28T14:08:47.803171Z","iopub.execute_input":"2024-04-28T14:08:47.803620Z","iopub.status.idle":"2024-04-28T14:08:48.394192Z","shell.execute_reply.started":"2024-04-28T14:08:47.803579Z","shell.execute_reply":"2024-04-28T14:08:48.393180Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The histograms for binding molecules shown above demonstrate the following:\n1. For protein `sEH` more than 80% of the binding small molecules sharing the same pair of building blocks have 2 or less alternatives (variants of the 3rd building block)\n2. For other 2 proteins more than 80% of the binding small molecules sharing the same pair of building blocks have 5 or less alternatives (variants of the 3rd building block)\n3. The right part of the plot indicates that there are molecules where binding occurs because of the pair of blocks that we fixed and doesn't really depend on the choice of the third one.","metadata":{}},{"cell_type":"markdown","source":"## Building blocks composition in non-binding molecules\n\n(This part is slow to compute, I might want to remove or skip it in the future)","metadata":{}},{"cell_type":"code","source":"train_all0_pq = pq.ParquetFile((SAVEDIR / \"train_all0.parquet\").as_posix())","metadata":{"execution":{"iopub.status.busy":"2024-04-28T14:08:48.498409Z","iopub.execute_input":"2024-04-28T14:08:48.498802Z","iopub.status.idle":"2024-04-28T14:08:48.505131Z","shell.execute_reply.started":"2024-04-28T14:08:48.498772Z","shell.execute_reply":"2024-04-28T14:08:48.503865Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"pairs_blocks_negative_counter = {\n name: defaultdict(int) for name in protein_names\n}\nfor protein_name in tqdm(protein_names):\n ids = train_mixed_df[protein_name] == 0\n if DRAFT_MODE:\n ids = ids & (train_mixed_df.index < 10000)\n for column_pair in pair_bb_columns:\n new_pairs, new_counts = np.unique(\n train_mixed_df.loc[ids, column_pair].apply(lambda row: tuple(sorted(row)), axis=1).values, \n return_counts=True\n )\n for pair, count in zip(new_pairs, new_counts):\n pairs_blocks_negative_counter[protein_name][pair] += count","metadata":{"execution":{"iopub.status.busy":"2024-04-28T14:08:48.961402Z","iopub.execute_input":"2024-04-28T14:08:48.961774Z","iopub.status.idle":"2024-04-28T14:08:49.801004Z","shell.execute_reply.started":"2024-04-28T14:08:48.961746Z","shell.execute_reply":"2024-04-28T14:08:49.799497Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# # the following takes about an hour to compute, uncomment if you want to see the whole picture on all the negative data\n# for rg in tqdm(np.arange(train_all0_pq.num_row_groups), total=train_all0_pq.num_row_groups):\n# batch = train_all0_pq.read_row_group(rg, columns=building_blocks_columns)\n# batch_df = batch.to_pandas()\n# for i, column in enumerate(building_blocks_columns):\n# batch_df[f\"index{i+1}\"] = batch_df[column].apply(lambda x: smiles2index[x])\n# for column_pair in pair_bb_columns:\n# new_pairs, new_counts = np.unique(\n# batch_df[list(column_pair)].apply(lambda row: tuple(sorted(row)), axis=1).values, \n# return_counts=True\n# )\n# for pair, count in zip(new_pairs, new_counts):\n# pairs_blocks_negative_counter[protein_name][pair] += count","metadata":{"execution":{"iopub.status.busy":"2024-04-28T14:08:49.803273Z","iopub.execute_input":"2024-04-28T14:08:49.804209Z","iopub.status.idle":"2024-04-28T14:08:49.810312Z","shell.execute_reply.started":"2024-04-28T14:08:49.804159Z","shell.execute_reply":"2024-04-28T14:08:49.808994Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"for protein_name, data in pairs_blocks_negative_counter.items():\n sns.histplot(\n x=list(data.values()), element=\"step\", \n log_scale=(False, True), \n fill=False, \n binwidth=10,\n label=protein_name\n )\n \n #plt.title(protein_name)\nplt.axvline(x=1, ls='--')\nplt.legend()\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2024-04-28T14:08:50.215938Z","iopub.execute_input":"2024-04-28T14:08:50.216348Z","iopub.status.idle":"2024-04-28T14:08:51.010931Z","shell.execute_reply.started":"2024-04-28T14:08:50.216319Z","shell.execute_reply":"2024-04-28T14:08:51.009778Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Filter negative data once more\n\nWe might want to filter negative samples which are 'too far' from positive samples. But what it really means? One of the possible approaches to define 'too far' by counting the comparing blocks in molecules: if we have non-binding molecule (A, B, C) and there are binding molecules which are composed from pairs of blocks (A, B), or (A, B), or (B, C) - it might be closer to positive data comparing to molecule (D, E, F) if there is no molecules with some of the pairs of blocks (D, E), (E, F), (D, E) in the binding molecules dataset.\n\nLet's split the negative samples to the following parts:\n1. Molecules which contain 3 pairs of building blocks present in positive samples;\n2. Molecules which contain 2 pairs of building blocks present in positive samples;\n3. Molecules which contain exactly 1 pair of building blocks\n4. Molecules which doesn't contain any pair of building blocks\n\nAll these datasets will be protein-specific, as a result we will have 12 files.\n\n","metadata":{}},{"cell_type":"code","source":"schema_single_target = pa.schema([\n (col_name, pa.string())\n for col_name in smiles_col_names\n])\n\nKEYS_INFO2 = [\n (f\"{protein_name}_{i}pairs_building_blocks\", schema_single_target)\n for protein_name in protein_names\n for i in np.arange(4)\n]","metadata":{"execution":{"iopub.status.busy":"2024-04-28T14:08:51.373935Z","iopub.execute_input":"2024-04-28T14:08:51.374660Z","iopub.status.idle":"2024-04-28T14:08:51.380653Z","shell.execute_reply.started":"2024-04-28T14:08:51.374624Z","shell.execute_reply":"2024-04-28T14:08:51.379752Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"protein_info = {\n protein_name: defaultdict(int)\n for protein_name in protein_names\n}\nwith PqDataSaverWrapper(savedir=SAVEDIR, keys_info=KEYS_INFO2, prefix=\"negative\") as pq_wrapper:\n for i in tqdm(np.arange(train_all0_pq.num_row_groups), total=train_all0_pq.num_row_groups):\n batch = train_all0_pq.read_row_group(i)\n batch_df = batch.to_pandas()\n for i, column in enumerate(building_blocks_columns):\n batch_df[f\"index{i+1}\"] = batch_df[column].apply(lambda x: smiles2index[x])\n batch_df[\"pair12\"] = [tuple(sorted(pair)) for pair in batch_df[['index1', 'index2']].values]\n batch_df[\"pair23\"] = [tuple(sorted(pair)) for pair in batch_df[['index2', 'index3']].values]\n batch_df[\"pair13\"] = [tuple(sorted(pair)) for pair in batch_df[['index1', 'index3']].values]\n new_data = {}\n for protein_name in protein_names:\n batch_df[f\"{protein_name}_pair12\"] = batch_df.pair12.apply(lambda x: x in pair_blocks_positive_counter[protein_name])\n batch_df[f\"{protein_name}_pair23\"] = batch_df.pair23.apply(lambda x: x in pair_blocks_positive_counter[protein_name])\n batch_df[f\"{protein_name}_pair13\"] = batch_df.pair13.apply(lambda x: x in pair_blocks_positive_counter[protein_name])\n # batch_df[f\"{protein_name}_match\"] = batch_df[f\"{protein_name}_pair12\"].values | batch_df[f\"{protein_name}_pair23\"].values | batch_df[f\"{protein_name}_pair13\"].values\n batch_df[f\"{protein_name}_match\"] = batch_df[[f\"{protein_name}_pair12\", f\"{protein_name}_pair23\", f\"{protein_name}_pair13\"]].sum(1)\n matches, counts = np.unique(batch_df[f\"{protein_name}_match\"], return_counts=True)\n # the following aggregates info on the whole train0.parquet\n for m, c in zip(matches, counts):\n protein_info[protein_name][m] += c\n ids = batch_df[f\"{protein_name}_match\"] == m\n if ids.sum() > 0:\n # break\n key_name = f\"{protein_name}_{m}pairs_building_blocks\"\n new_data[key_name] = batch_df.loc[ids, smiles_col_names].reset_index(drop=True)\n if len(new_data) > 0:\n pq_wrapper.add_data(new_data)\n\n # break","metadata":{"execution":{"iopub.status.busy":"2024-04-28T14:08:51.807403Z","iopub.execute_input":"2024-04-28T14:08:51.809893Z","iopub.status.idle":"2024-04-28T14:10:11.001033Z","shell.execute_reply.started":"2024-04-28T14:08:51.809854Z","shell.execute_reply":"2024-04-28T14:10:10.999964Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig = plt.figure(figsize=(5, 5))\nax = plt.gca()\nwidth = 1\n\nfor offset, (protein, data) in enumerate(protein_info.items()):\n xs = np.asarray(sorted(data))\n ys = [data[x] for x in xs]\n p = ax.bar(xs*width*4+offset, ys, label=protein, alpha=0.4, width=width)#, edgecolor='purple', color='None')\nax.set_xticks(np.arange(1,14, 4))\nax.set_xticklabels([0, 1, 2, 3])\nax.set_title(\"Number of molecules with matched pairs of blocks in train_all0.parquet\")\nplt.legend()\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2024-04-28T14:10:11.003142Z","iopub.execute_input":"2024-04-28T14:10:11.003482Z","iopub.status.idle":"2024-04-28T14:10:11.428284Z","shell.execute_reply.started":"2024-04-28T14:10:11.003453Z","shell.execute_reply":"2024-04-28T14:10:11.427070Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We might decide to omit the molecules from negative samples as being 'too far' if they have 0 number of matched pair of blocks (leftmost group of bars in the barplot above).","metadata":{}},{"cell_type":"markdown","source":"# Compute fingerprints for all unique SMILES strings appearing in building blocks (unfinished)","metadata":{}},{"cell_type":"markdown","source":"Our dataset includes isomeric SMILES, if we want to take this into account we could use fingerprints which takes chirality into account. I will use [map4c fingerprint](https://github.com/reymond-group/mapchiral) here.\n\nOne of the possible goals of this is train data filtering: it's easier to compute the fingerprints for 2k block smiles than to entire set of molecules.\n\nLater we might use these to define an approximate distance between pairs of molecules which share the same pair of building blocks, and differ only in one building block, i.e. (A, B, C) and (A, B, D) - as a distance between fingerprints of C and D.","metadata":{}},{"cell_type":"markdown","source":"https://www.daylight.com/dayhtml/doc/theory/theory.smiles.html\n> Disconnected compounds are written as individual structures separated by a \".\" (period). The order in which ions or ligands are listed is arbitrary. \n\nWe might want to clean smiles from building blocks before computing fingerprints.\n","metadata":{}},{"cell_type":"code","source":"smiles_with_disconnected = [smiles for smiles in unique_smiles if smiles.find(\".\") >= 0]\nprint(len(smiles_with_disconnected), \"unique smiles strings with disconnected structures\")","metadata":{"execution":{"iopub.status.busy":"2024-04-28T14:10:11.430999Z","iopub.execute_input":"2024-04-28T14:10:11.431931Z","iopub.status.idle":"2024-04-28T14:10:11.438695Z","shell.execute_reply.started":"2024-04-28T14:10:11.431897Z","shell.execute_reply":"2024-04-28T14:10:11.437493Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"smiles, counts = np.unique([x for smiles in smiles_with_disconnected for x in smiles.split(\".\") if len(x) < 5], return_counts=True)\nidx = np.argsort(counts)[::-1]\nsmiles[idx], counts[idx] # most common ones","metadata":{"execution":{"iopub.status.busy":"2024-04-28T14:10:11.441361Z","iopub.execute_input":"2024-04-28T14:10:11.441820Z","iopub.status.idle":"2024-04-28T14:10:11.460850Z","shell.execute_reply.started":"2024-04-28T14:10:11.441788Z","shell.execute_reply":"2024-04-28T14:10:11.459465Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"[smiles for smiles in smiles_with_disconnected if 'O' in smiles.split(\".\")]","metadata":{"execution":{"iopub.status.busy":"2024-04-28T14:10:11.462220Z","iopub.execute_input":"2024-04-28T14:10:11.462562Z","iopub.status.idle":"2024-04-28T14:10:11.475386Z","shell.execute_reply.started":"2024-04-28T14:10:11.462534Z","shell.execute_reply":"2024-04-28T14:10:11.473526Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Chem.MolFromSmiles(\"Cc1nc2[nH]cc(C(=O)O)c(=O)n2n1.O\") # just water - might be removed\n# Chem.MolFromSmiles(\"Cl.O.O=C(O)Cc1cn2ccsc2n1\") # water and HCl - might be removed","metadata":{"execution":{"iopub.status.busy":"2024-04-28T14:10:11.476477Z","iopub.execute_input":"2024-04-28T14:10:11.477092Z","iopub.status.idle":"2024-04-28T14:10:11.488700Z","shell.execute_reply.started":"2024-04-28T14:10:11.477048Z","shell.execute_reply":"2024-04-28T14:10:11.487432Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print([smiles for smiles in smiles_with_disconnected if 'Br' in smiles.split(\".\")])\n# Chem.MolFromSmiles(\"Br.NCc1cccc(Br)n1\")\n# Chem.MolFromSmiles(\"Br.Br.NCC1CCCN1c1cccnn1\")","metadata":{"execution":{"iopub.status.busy":"2024-04-28T14:10:11.489847Z","iopub.execute_input":"2024-04-28T14:10:11.490243Z","iopub.status.idle":"2024-04-28T14:10:11.503132Z","shell.execute_reply.started":"2024-04-28T14:10:11.490211Z","shell.execute_reply":"2024-04-28T14:10:11.502063Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"len([smiles for smiles in smiles_with_disconnected if 'Cl' in smiles.split(\".\")])","metadata":{"execution":{"iopub.status.busy":"2024-04-28T14:10:11.504228Z","iopub.execute_input":"2024-04-28T14:10:11.505078Z","iopub.status.idle":"2024-04-28T14:10:11.524590Z","shell.execute_reply.started":"2024-04-28T14:10:11.505026Z","shell.execute_reply":"2024-04-28T14:10:11.523407Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"unique_smiles_cleaned = [\".\".join([x for x in smiles.split(\".\") if x not in [\"Br\", \"Cl\", \"O\"]]) for smiles in unique_smiles]","metadata":{"execution":{"iopub.status.busy":"2024-04-28T14:10:11.527009Z","iopub.execute_input":"2024-04-28T14:10:11.528115Z","iopub.status.idle":"2024-04-28T14:10:11.538533Z","shell.execute_reply.started":"2024-04-28T14:10:11.528068Z","shell.execute_reply":"2024-04-28T14:10:11.537205Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print([smiles for smiles in unique_smiles_cleaned if smiles.find(\".\") >=0])\nChem.MolFromSmiles(\"NCC(O)c1ccc(Cl)s1.O=C(O)C(F)(F)F\")","metadata":{"execution":{"iopub.status.busy":"2024-04-28T14:10:11.542492Z","iopub.execute_input":"2024-04-28T14:10:11.542870Z","iopub.status.idle":"2024-04-28T14:10:11.565978Z","shell.execute_reply.started":"2024-04-28T14:10:11.542838Z","shell.execute_reply":"2024-04-28T14:10:11.564810Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"Chem.MolFromSmiles(\"Nc1ccsc1.O=C(O)C(=O)O\")","metadata":{"execution":{"iopub.status.busy":"2024-04-28T14:10:11.567452Z","iopub.execute_input":"2024-04-28T14:10:11.568551Z","iopub.status.idle":"2024-04-28T14:10:11.587529Z","shell.execute_reply.started":"2024-04-28T14:10:11.568506Z","shell.execute_reply":"2024-04-28T14:10:11.586275Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Not sure what to do with these two disconnected molecules, for now let's keep them","metadata":{}},{"cell_type":"code","source":"np.unique(unique_smiles_cleaned, return_counts=True)[1].max() \n# 1 means that no duplicates appeared after cleaning!","metadata":{"execution":{"iopub.status.busy":"2024-04-28T14:10:11.589175Z","iopub.execute_input":"2024-04-28T14:10:11.589825Z","iopub.status.idle":"2024-04-28T14:10:11.602231Z","shell.execute_reply.started":"2024-04-28T14:10:11.589779Z","shell.execute_reply":"2024-04-28T14:10:11.600886Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from rdkit import Chem\nfrom mapchiral.mapchiral import encode, jaccard_similarity\n\nunique_molecules = [Chem.MolFromSmiles(smiles) for smiles in unique_smiles_cleaned]\n\nfingerprints = []\nfor mol in tqdm(unique_molecules):\n fp = encode(mol, max_radius=2, n_permutations=2048, mapping=False)\n fingerprints.append(fp)\nfingerprints = np.stack(fingerprints)\n\n# let's save them in case we will want to reuse them in the future\nnp.savez_compressed(\n \"unique_smiles/map4c_fingerprints.npz\", \n fp=fingerprints, \n blocks=unique_molecules,\n blocks_cleaned=unique_smiles_cleaned\n)","metadata":{"execution":{"iopub.status.busy":"2024-04-28T14:10:42.987618Z","iopub.execute_input":"2024-04-28T14:10:42.989250Z","iopub.status.idle":"2024-04-28T14:11:06.420919Z","shell.execute_reply.started":"2024-04-28T14:10:42.989196Z","shell.execute_reply":"2024-04-28T14:11:06.419724Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Fragment analysis\n\nSum of building blocks is not equal to the molecule for which the binding status is measured. However, we can use the blocks to split the molecule to the triazine core and fragments and understand how it was constructed: \n- At what positions the blocks connected to triazine core? \n- Are these positions always the same or not ?\n- Is the linker always connected the same way or not?","metadata":{}},{"cell_type":"markdown","source":"For all our future experiments we'll be using triazine core:","metadata":{}},{"cell_type":"code","source":"def make_query(mol):\n # the following lines are needed to make * work with 'GetSubstructMatches' method - \n p = Chem.AdjustQueryParameters.NoAdjustments()\n p.makeDummiesQueries = True\n mol = Chem.AdjustQueryProperties(mol, p)\n mol.UpdatePropertyCache()\n return mol\n\ntriazine_core = Chem.MolFromSmiles(\"c(*)1nc(*)nc(*)n1\")\ntriazine_core = make_query(triazine_core)\ntriazine_core","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2024-04-28T14:11:12.643377Z","iopub.execute_input":"2024-04-28T14:11:12.643776Z","iopub.status.idle":"2024-04-28T14:11:12.660241Z","shell.execute_reply.started":"2024-04-28T14:11:12.643746Z","shell.execute_reply":"2024-04-28T14:11:12.659044Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Fragment extraction: example 1\nLet's take a closer look how can we split the molecule to fragments and match it with building blocks","metadata":{}},{"cell_type":"code","source":"example_index = 0\nexample_smiles = train_mixed_df.loc[example_index, \"molecule_smiles\"]\n# example_building_blocks_smiles = train_mixed_df.loc[example_index, building_blocks_columns].values\nprint(\"Molecule:\", example_smiles)\nmolecule = Chem.MolFromSmiles(example_smiles)\nChem.Draw.MolsToGridImage((molecule,), subImgSize=(500,500))","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2024-04-28T14:11:19.324349Z","iopub.execute_input":"2024-04-28T14:11:19.324739Z","iopub.status.idle":"2024-04-28T14:11:19.380915Z","shell.execute_reply.started":"2024-04-28T14:11:19.324712Z","shell.execute_reply":"2024-04-28T14:11:19.379606Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"all_matches = []\nfor match in molecule.GetSubstructMatches(triazine_core):\n all_matches.append(match)\n \n\n# Chem.Draw.MolsToImage(all_splits, subImgSize=(400, 600), legends=('core','another one'))\nChem.Draw.MolsToGridImage(\n [molecule]*len(all_matches), \n subImgSize=(300, 300),\n highlightAtomLists=all_matches\n)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2024-04-28T14:11:21.987007Z","iopub.execute_input":"2024-04-28T14:11:21.987382Z","iopub.status.idle":"2024-04-28T14:11:22.019049Z","shell.execute_reply.started":"2024-04-28T14:11:21.987353Z","shell.execute_reply":"2024-04-28T14:11:22.017997Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Fragment extraction: example 2\nFragment which looks like triazine core might appear in molecule more than once!","metadata":{}},{"cell_type":"code","source":"# below is the molecule which corresponds to the following code:\n# example_index = 1083\n# example_smiles = train_mixed_df.loc[example_index, \"molecule_smiles\"]\n# example_building_blocks_smiles = train_mixed_df.loc[example_index, building_blocks_columns].values\n\nexample_smiles = \"C#CC[C@@H](CC(=O)N[Dy])Nc1nc(NCc2cncc(F)c2)nc(Nc2nc(C)nc(OC)n2)n1\"\nexample_building_blocks_smiles = [\n 'C#CC[C@@H](CC(=O)O)NC(=O)OCC1c2ccccc2-c2ccccc21', \n 'Cl.Cl.NCc1cncc(F)c1',\n 'COc1nc(C)nc(N)n1'\n]\n\n# example_smiles = \"C#CCOc1ccc(CNc2nc(NCc3ccc4[nH]c(C)cc4c3)nc(N[C@@H](CC#C)CC(=O)N[Dy])n2)cc1\"\n# example_building_blocks_smiles = [\n# 'C#CC[C@@H](CC(=O)O)NC(=O)OCC1c2ccccc2-c2ccccc21', \n# 'C#CCOc1ccc(CN)cc1.Cl',\n# 'Cc1cc2cc(CN)ccc2[nH]1'\n# ]\nprint(\"Molecule:\", example_smiles)\n\nmolecule = Chem.MolFromSmiles(example_smiles)\n\nChem.Draw.MolsToGridImage((molecule,), subImgSize=(500,500))","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2024-04-28T14:11:22.578384Z","iopub.execute_input":"2024-04-28T14:11:22.579042Z","iopub.status.idle":"2024-04-28T14:11:22.633143Z","shell.execute_reply.started":"2024-04-28T14:11:22.579010Z","shell.execute_reply":"2024-04-28T14:11:22.631966Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(\"Building blocks\")\n# print(example_building_blocks_smiles)\nbb_molecules = [Chem.MolFromSmiles(sm)for sm in example_building_blocks_smiles]\nChem.Draw.MolsToGridImage(\n bb_molecules, subImgSize=(500,500), \n legends=example_building_blocks_smiles)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2024-04-28T14:11:22.963768Z","iopub.execute_input":"2024-04-28T14:11:22.964478Z","iopub.status.idle":"2024-04-28T14:11:23.029743Z","shell.execute_reply.started":"2024-04-28T14:11:22.964439Z","shell.execute_reply":"2024-04-28T14:11:23.028375Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Possible candidates for 'the triazine core' which was used to build the molecule from building blocks are highlighted with red color:","metadata":{}},{"cell_type":"code","source":"all_matches = []\nfor match in molecule.GetSubstructMatches(triazine_core):\n all_matches.append(match)\n \n\n# Chem.Draw.MolsToImage(all_splits, subImgSize=(400, 600), legends=('core','another one'))\nChem.Draw.MolsToGridImage(\n [molecule]*len(all_matches), \n subImgSize=(400, 400),\n highlightAtomLists=all_matches\n)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2024-04-28T14:11:23.896403Z","iopub.execute_input":"2024-04-28T14:11:23.896839Z","iopub.status.idle":"2024-04-28T14:11:23.949235Z","shell.execute_reply.started":"2024-04-28T14:11:23.896804Z","shell.execute_reply":"2024-04-28T14:11:23.948244Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"How can we identify which triazine core was used for molecule construction? One possible way is to split the molecule to to fragments, score them against building blocks and find the match between fragments and blocks with greatest possible score.","metadata":{}},{"cell_type":"code","source":"print(\"Fragments with triazine core removed:\")\nall_splits = []\nfor match in molecule.GetSubstructMatches(triazine_core):\n fragments = Chem.ReplaceCore(molecule, triazine_core, match, replaceDummies=True, requireDummyMatch=True)\n all_splits.append(fragments)\nChem.Draw.MolsToGridImage(all_splits, subImgSize=(500,500), legends=(\"fragments set 1\", \"fragments set 2\"))","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2024-04-28T14:11:25.073638Z","iopub.execute_input":"2024-04-28T14:11:25.074037Z","iopub.status.idle":"2024-04-28T14:11:25.137842Z","shell.execute_reply.started":"2024-04-28T14:11:25.074005Z","shell.execute_reply":"2024-04-28T14:11:25.136699Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from rdkit.Chem import rdFMCS\n\nfragments_set = list(Chem.GetMolFrags(all_splits[0], asMols=True, sanitizeFrags=True))\n\nmol_pairs = []\nall_hit_atoms = []\nall_hit_bonds = []\nlegends = []\nfor i in range(len(fragments_set)):\n fragment = fragments_set[i]\n for j in range(3):\n bb_mol = bb_molecules[j]\n \n match_result = rdFMCS.FindMCS(\n [fragment, bb_mol], \n matchValences=True,\n ringMatchesRingOnly=True,\n completeRingsOnly=True\n )\n pattern = Chem.MolFromSmarts(match_result.smartsString)\n \n hit_atoms = list(fragment.GetSubstructMatch(pattern))\n hit_bonds = []\n for bond in pattern.GetBonds():\n aid1 = hit_atoms[bond.GetBeginAtomIdx()]\n aid2 = hit_atoms[bond.GetEndAtomIdx()]\n hit_bonds.append(fragment.GetBondBetweenAtoms(aid1,aid2).GetIdx())\n all_hit_atoms.append(hit_atoms)\n all_hit_bonds.append(hit_bonds)\n \n hit_atoms = list(bb_mol.GetSubstructMatch(pattern))\n hit_bonds = []\n for bond in pattern.GetBonds():\n aid1 = hit_atoms[bond.GetBeginAtomIdx()]\n aid2 = hit_atoms[bond.GetEndAtomIdx()]\n hit_bonds.append(bb_mol.GetBondBetweenAtoms(aid1,aid2).GetIdx())\n all_hit_atoms.append(hit_atoms)\n all_hit_bonds.append(hit_bonds)\n \n d = max(fragment.GetNumAtoms(), bb_mol.GetNumAtoms())\n has_linker1 = np.any([a.GetSymbol() == 'Dy' for a in fragment.GetAtoms()])\n if has_linker1:\n d = fragment.GetNumAtoms()\n has_linker2 = np.any([a.GetSymbol() == 'Dy' for a in bb_mol.GetAtoms()])\n if has_linker2:\n d = bb_mol.GetNumAtoms()\n score = match_result.numAtoms / d\n mol_pairs.extend([fragment, bb_mol])\n legends.extend([\"\", f\"fragment {i} and building block {j}: score {score:3.3f}\"])\n \n\nChem.Draw.MolsToGridImage(\n mol_pairs, \n subImgSize=(300,300), \n legends=legends,\n molsPerRow=6,\n highlightAtomLists=all_hit_atoms\n)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2024-04-28T14:11:25.633440Z","iopub.execute_input":"2024-04-28T14:11:25.633822Z","iopub.status.idle":"2024-04-28T14:11:25.823758Z","shell.execute_reply.started":"2024-04-28T14:11:25.633793Z","shell.execute_reply":"2024-04-28T14:11:25.822647Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Let's check out scores computed for the second set of fragments:","metadata":{}},{"cell_type":"code","source":"fragments_set = list(Chem.GetMolFrags(all_splits[1], asMols=True, sanitizeFrags=True))\n\nmol_pairs = []\nall_hit_atoms = []\nall_hit_bonds = []\nlegends = []\nfor i in range(len(fragments_set)):\n fragment = fragments_set[i]\n for j in range(3):\n bb_mol = bb_molecules[j]\n \n match_result = rdFMCS.FindMCS(\n [fragment, bb_mol], \n matchValences=True,\n ringMatchesRingOnly=True,\n completeRingsOnly=True\n )\n pattern = Chem.MolFromSmarts(match_result.smartsString)\n \n hit_atoms = list(fragment.GetSubstructMatch(pattern))\n hit_bonds = []\n for bond in pattern.GetBonds():\n aid1 = hit_atoms[bond.GetBeginAtomIdx()]\n aid2 = hit_atoms[bond.GetEndAtomIdx()]\n hit_bonds.append(fragment.GetBondBetweenAtoms(aid1,aid2).GetIdx())\n all_hit_atoms.append(hit_atoms)\n all_hit_bonds.append(hit_bonds)\n \n hit_atoms = list(bb_mol.GetSubstructMatch(pattern))\n hit_bonds = []\n for bond in pattern.GetBonds():\n aid1 = hit_atoms[bond.GetBeginAtomIdx()]\n aid2 = hit_atoms[bond.GetEndAtomIdx()]\n hit_bonds.append(bb_mol.GetBondBetweenAtoms(aid1,aid2).GetIdx())\n all_hit_atoms.append(hit_atoms)\n all_hit_bonds.append(hit_bonds)\n \n d = max(fragment.GetNumAtoms(), bb_mol.GetNumAtoms())\n has_linker1 = np.any([a.GetSymbol() == 'Dy' for a in fragment.GetAtoms()])\n if has_linker1:\n d = fragment.GetNumAtoms()\n has_linker2 = np.any([a.GetSymbol() == 'Dy' for a in bb_mol.GetAtoms()])\n if has_linker2:\n d = bb_mol.GetNumAtoms()\n score = match_result.numAtoms / d\n mol_pairs.extend([fragment, bb_mol])\n legends.extend([\"\", f\"fragment {i} and building block {j}: score {score:3.3f}\"])\n \n\nChem.Draw.MolsToGridImage(\n mol_pairs, \n subImgSize=(300,300), \n legends=legends,\n molsPerRow=6,\n highlightAtomLists=all_hit_atoms\n)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2024-04-28T14:11:26.901986Z","iopub.execute_input":"2024-04-28T14:11:26.902412Z","iopub.status.idle":"2024-04-28T14:11:27.045114Z","shell.execute_reply.started":"2024-04-28T14:11:26.902359Z","shell.execute_reply":"2024-04-28T14:11:27.043901Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"It is clear that after we locate triazine core we should get 3 fragments; each of the fragments should be scored against building blocks 1-3. After that we'll solve this as a linear sum assignment problem as implemented in [scipy](https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.linear_sum_assignment.html) and consider the returned assignment to be the best possible match.","metadata":{}},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Fragment extraction and match: code\nThe finalized code for molecule splitting relies on rdkit's `FindMCS` method to find matches, it might not work in some cases, despite that, it's the 'current best' one.","metadata":{}},{"cell_type":"code","source":"from scipy.optimize import linear_sum_assignment\n\n\ndef make_query(mol):\n p = Chem.AdjustQueryParameters.NoAdjustments()\n p.makeDummiesQueries = True\n mol = Chem.AdjustQueryProperties(mol, p)\n mol.UpdatePropertyCache()\n return mol\n\n\ndef remove_dummy_number(m):\n star = Chem.MolFromSmiles(\"*\")\n new_mol, = Chem.ReplaceSubstructs(m, star, star, replaceAll=True)\n return new_mol\n\n\ndef split_by_triazine_core(smiles, debug=False):\n molecule = Chem.MolFromSmiles(smiles)\n triazine_core = make_query(Chem.MolFromSmiles(\"c1ncncn1\"))\n for match in molecule.GetSubstructMatches(triazine_core):\n fragments = Chem.ReplaceCore(\n molecule, \n triazine_core,\n match, \n replaceDummies=True, \n requireDummyMatch=False\n )\n fragments = [\n remove_dummy_number(x) \n for x in Chem.GetMolFrags(fragments, asMols=True, sanitizeFrags=True)\n ]\n fragments = [Chem.MolToSmiles(x) for x in fragments]\n if len(fragments) != 3:\n if debug:\n print(fragments)\n continue\n yield fragments\n\n\ndef score_pair(mol1, mol2, use_mcs=True, debug=False):\n if use_mcs:\n from rdkit.Chem import rdFMCS\n match_result = rdFMCS.FindMCS(\n [mol1, mol2],\n matchValences=True, \n ringMatchesRingOnly=True,\n completeRingsOnly=True,\n bondCompare=rdFMCS.BondCompare.CompareOrderExact\n )\n d = max(mol1.GetNumAtoms(), mol2.GetNumAtoms())\n has_linker1 = np.any([a.GetSymbol() == 'Dy' for a in mol1.GetAtoms()])\n if has_linker1:\n d = mol1.GetNumAtoms()\n has_linker2 = np.any([a.GetSymbol() == 'Dy' for a in mol2.GetAtoms()])\n if has_linker2:\n d = mol2.GetNumAtoms()\n return match_result.numAtoms/d\n\n matches = mol2.GetSubstructMatches(mol1, useChirality=True)\n if len(matches) == 0:\n return 0.\n r = max([len(m) for m in matches])\n return r / d\n\n\ndef score_sets(frag_mols, bb_mols, debug=False):\n scores = [[\n score_pair(a, b, debug=debug) for b in bb_mols\n ]\n for a in frag_mols\n ]\n frag_ind, bb_ind = linear_sum_assignment(scores, maximize=True)\n if debug:\n print(\"\\nScores:\")\n print(scores)\n print(frag_ind, bb_ind)\n score = np.sum([scores[i][j] for i, j in zip(frag_ind, bb_ind)])\n if debug:\n print(\"\\nScore:\")\n print([scores[i][j] for i, j in zip(frag_ind, bb_ind)])\n print(score)\n # print(scores)\n # print(row_ind, col_ind)\n idx = np.argsort(bb_ind)\n bb_ind = bb_ind[idx]\n frag_ind = frag_ind[idx]\n # ordered molecular fragments corresponding to building blocks 0, 1, 2\n return (frag_ind, bb_ind), score \n\n\ndef match_fragments(smiles, smiles_bb, debug=False):\n mol_bb = [Chem.MolFromSmiles(bb, sanitize=True) for bb in smiles_bb]\n mol_bb = [Chem.AddHs(mol) for mol in mol_bb]\n max_score = 0\n indices = None\n assignment = None\n for fragments in split_by_triazine_core(smiles):\n if debug:\n print(\"split\", fragments)\n frag_mols = [Chem.MolFromSmiles(smiles) for smiles in fragments]\n frag_mols = [Chem.AddHs(mol) for mol in frag_mols]\n (frag_ind, bb_ind), score = score_sets(frag_mols, mol_bb, debug=debug)\n if debug:\n print(\"Score\", score)\n if score > max_score:\n max_score = score\n indices = (bb_ind, frag_ind)\n assignment = [fragments[i] for i in frag_ind]\n return assignment\n\n\n# usage example: the fragments in the output correspond to matching molecules in building blocks\nmatch_fragments(example_smiles, example_building_blocks_smiles)","metadata":{"execution":{"iopub.status.busy":"2024-04-28T14:11:31.640917Z","iopub.execute_input":"2024-04-28T14:11:31.641365Z","iopub.status.idle":"2024-04-28T14:11:31.685669Z","shell.execute_reply.started":"2024-04-28T14:11:31.641334Z","shell.execute_reply":"2024-04-28T14:11:31.684242Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Fragment extraction and match: partial processing\n\nWe are mostly interested in binding molecules from the train dataset. It is slow to compute for all the molecules, so I take a small subset here (with one fixed `buildingblock1_smiles`). If my fragmentation code is more or less valid, it will allow to compute fragments which were extracted from building blocks.\n\nLater I will run this on all binding molecules, all binding molecules for one of the proteins, and probably on the test data - if it is not too time consuming. If it takes a lot of time, I'll add a link to the separate notebook here.","metadata":{}},{"cell_type":"code","source":"buildingblock1_smiles = train_mixed_df.buildingblock1_smiles.unique()[2]\nbuildingblock1_smiles = train_mixed_df.buildingblock1_smiles.value_counts()[5:].index[0]\ntrain_mixed_df.buildingblock1_smiles.value_counts()[5:]","metadata":{"execution":{"iopub.status.busy":"2024-04-28T14:11:32.931536Z","iopub.execute_input":"2024-04-28T14:11:32.931993Z","iopub.status.idle":"2024-04-28T14:11:32.953503Z","shell.execute_reply.started":"2024-04-28T14:11:32.931957Z","shell.execute_reply":"2024-04-28T14:11:32.952246Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"ids = train_mixed_df.buildingblock1_smiles == buildingblock1_smiles\nbuildingblock1_smiles, ids.sum()","metadata":{"execution":{"iopub.status.busy":"2024-04-28T14:11:33.504074Z","iopub.execute_input":"2024-04-28T14:11:33.504523Z","iopub.status.idle":"2024-04-28T14:11:33.517453Z","shell.execute_reply.started":"2024-04-28T14:11:33.504489Z","shell.execute_reply":"2024-04-28T14:11:33.516282Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# all_smiles = set()\nfor index in tqdm(train_mixed_df[ids].index[:10000]):\n smiles = train_mixed_df.loc[index, \"molecule_smiles\"]\n smiles_bb = [train_mixed_df.loc[index, f\"buildingblock{i+1}_smiles\"] for i in range(3)]\n match_result = match_fragments(smiles, smiles_bb)\n if match_result is None:\n print(f\"Couldn't split molecule '{smiles}' to fragments\")\n continue\n b1, b2, b3 = match_result\n train_mixed_df.loc[index, \"frag1\"] = b1\n train_mixed_df.loc[index, \"frag2\"] = b2\n train_mixed_df.loc[index, \"frag3\"] = b3\n # all_smiles.update([b1, b2, b3])","metadata":{"execution":{"iopub.status.busy":"2024-04-28T14:11:34.972895Z","iopub.execute_input":"2024-04-28T14:11:34.973291Z","iopub.status.idle":"2024-04-28T14:12:18.590447Z","shell.execute_reply.started":"2024-04-28T14:11:34.973262Z","shell.execute_reply":"2024-04-28T14:12:18.589276Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"ids = ids & (train_mixed_df[[\"frag1\", \"frag2\", \"frag3\"]].isnull().sum(1) < 1)\ndf3 = train_mixed_df[ids].groupby(\"buildingblock3_smiles\", as_index=False).agg({\n 'frag3': lambda x: np.unique(x)\n})\ndf3.frag3.apply(len).value_counts()","metadata":{"execution":{"iopub.status.busy":"2024-04-28T14:12:22.642197Z","iopub.execute_input":"2024-04-28T14:12:22.642633Z","iopub.status.idle":"2024-04-28T14:12:22.712701Z","shell.execute_reply.started":"2024-04-28T14:12:22.642599Z","shell.execute_reply":"2024-04-28T14:12:22.711459Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df2 = train_mixed_df[ids].groupby(\"buildingblock2_smiles\", as_index=False).agg({\n 'frag2': lambda x: np.unique(x)\n})\ndf2.frag2.apply(len).value_counts()","metadata":{"execution":{"iopub.status.busy":"2024-04-28T14:12:23.537910Z","iopub.execute_input":"2024-04-28T14:12:23.539069Z","iopub.status.idle":"2024-04-28T14:12:23.593896Z","shell.execute_reply.started":"2024-04-28T14:12:23.539028Z","shell.execute_reply":"2024-04-28T14:12:23.593013Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df1 = train_mixed_df[ids].groupby(\"buildingblock1_smiles\", as_index=False).agg({\n 'frag1': lambda x: np.unique(x)\n})\ndf1.frag1.apply(len).value_counts()","metadata":{"execution":{"iopub.status.busy":"2024-04-28T14:12:25.797379Z","iopub.execute_input":"2024-04-28T14:12:25.797783Z","iopub.status.idle":"2024-04-28T14:12:25.816199Z","shell.execute_reply.started":"2024-04-28T14:12:25.797752Z","shell.execute_reply":"2024-04-28T14:12:25.815306Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"If any of the above groups contain more than one variant of fragment for building block, it means that something is not matched correctly OR there are more than one way to use some specific block as a part of the molecule.\n\nThe following (hidden) fragments of code were used for debug purposes.","metadata":{}},{"cell_type":"code","source":"max_index = df3.frag3.apply(len).idxmax()\nbuildingblock3_smiles = df3.loc[max_index, \"buildingblock3_smiles\"]\n\nids2 = ids & (train_mixed_df[\"buildingblock3_smiles\"] == buildingblock3_smiles)\nfrag_columns = [f\"frag{i+1}\" for i in range(3)]\ndf3 = train_mixed_df.loc[ids2, building_blocks_columns+[\"molecule_smiles\"]+frag_columns].groupby(\"frag3\", as_index=False).agg(\"first\")\ndf3.head()","metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"execution":{"iopub.status.busy":"2024-04-28T14:12:36.249277Z","iopub.execute_input":"2024-04-28T14:12:36.249720Z","iopub.status.idle":"2024-04-28T14:12:36.285689Z","shell.execute_reply.started":"2024-04-28T14:12:36.249688Z","shell.execute_reply":"2024-04-28T14:12:36.284598Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def molecular_pair_with_matches(smiles1, smiles2):\n mol1 = Chem.MolFromSmiles(smiles1)\n mol1 = make_query(mol1)\n mol2 = Chem.MolFromSmiles(smiles2)\n mol2 = make_query(mol2)\n match_result = rdFMCS.FindMCS(\n [mol1, mol2], \n matchValences=True,\n ringMatchesRingOnly=True,\n completeRingsOnly=True\n )\n pattern = Chem.MolFromSmarts(match_result.smartsString)\n \n hit_atoms1 = list(mol1.GetSubstructMatch(pattern))\n hit_bonds1 = []\n for bond in pattern.GetBonds():\n aid1 = hit_atoms1[bond.GetBeginAtomIdx()]\n aid2 = hit_atoms1[bond.GetEndAtomIdx()]\n hit_bonds1.append(mol1.GetBondBetweenAtoms(aid1, aid2).GetIdx())\n\n hit_atoms2 = list(mol2.GetSubstructMatch(pattern))\n hit_bonds2 = []\n for bond in pattern.GetBonds():\n aid1 = hit_atoms2[bond.GetBeginAtomIdx()]\n aid2 = hit_atoms2[bond.GetEndAtomIdx()]\n hit_bonds2.append(mol2.GetBondBetweenAtoms(aid1, aid2).GetIdx())\n \n return mol1, hit_atoms1, hit_bonds1, mol2, hit_atoms2, hit_bonds2\n \n\nall_molecules = []\nall_hit_atoms = []\nall_legends = []\nfor i, row in df3.iterrows():\n for k in range(3):\n frag_smiles = row[f\"frag{k+1}\"]\n bb_smiles = row[f\"buildingblock{k+1}_smiles\"]\n mol1, hit_atoms1, hit_bonds1, mol2, hit_atoms2, hit_bonds2 = molecular_pair_with_matches(frag_smiles, bb_smiles)\n all_molecules.extend([mol1, mol2])\n all_hit_atoms.extend([hit_atoms1, hit_atoms2])\n all_legends.extend([f\"fragment {k+1}\", f\"building block {k+1}\"])\n # break\n \nChem.Draw.MolsToGridImage(\n all_molecules, \n subImgSize=(300,300), \n legends=all_legends,\n molsPerRow=6,\n highlightAtomLists=all_hit_atoms\n)","metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"execution":{"iopub.status.busy":"2024-04-28T14:12:38.322548Z","iopub.execute_input":"2024-04-28T14:12:38.323031Z","iopub.status.idle":"2024-04-28T14:12:38.420373Z","shell.execute_reply.started":"2024-04-28T14:12:38.322984Z","shell.execute_reply":"2024-04-28T14:12:38.419230Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"smiles = df3['molecule_smiles'][0]\nsmiles_bb = df3[building_blocks_columns].values[0]\nmatch_fragments(smiles, smiles_bb)","metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"execution":{"iopub.status.busy":"2024-04-28T14:12:39.693018Z","iopub.execute_input":"2024-04-28T14:12:39.693552Z","iopub.status.idle":"2024-04-28T14:12:39.715677Z","shell.execute_reply.started":"2024-04-28T14:12:39.693518Z","shell.execute_reply":"2024-04-28T14:12:39.714429Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"mol_bb = [Chem.MolFromSmiles(bb) for bb in smiles_bb]\nall_splits = []\nfor fragments in split_by_triazine_core(smiles):\n frag_mols = [Chem.MolFromSmiles(smiles) for smiles in fragments]\n frag_mols = [(mol) for mol in frag_mols]\n (frag_ind, bb_ind), score = score_sets(frag_mols, mol_bb)\n all_splits.append(frag_mols)\n print(bb_ind, frag_ind, score)","metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"execution":{"iopub.status.busy":"2024-04-28T14:12:40.865673Z","iopub.execute_input":"2024-04-28T14:12:40.866106Z","iopub.status.idle":"2024-04-28T14:12:40.883494Z","shell.execute_reply.started":"2024-04-28T14:12:40.866074Z","shell.execute_reply":"2024-04-28T14:12:40.882528Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from rdkit.Chem import rdFMCS\n\n# fragments_set = list(Chem.GetMolFrags(all_splits[0], asMols=True, sanitizeFrags=True))\nfragments_set = all_splits[0]\n\nmol_pairs = []\nall_hit_atoms = []\nall_hit_bonds = []\nlegends = []\nfor i in range(len(fragments_set)):\n fragment = Chem.AddHs(fragments_set[i])\n for j in range(3):\n bb_mol = Chem.AddHs(mol_bb[j])\n \n match_result = rdFMCS.FindMCS(\n [fragment, bb_mol], \n matchValences=True,\n ringMatchesRingOnly=True,\n completeRingsOnly=True,\n bondCompare=rdFMCS.BondCompare.CompareOrderExact\n )\n pattern = Chem.MolFromSmarts(match_result.smartsString)\n \n hit_atoms = list(fragment.GetSubstructMatch(pattern))\n hit_bonds = []\n for bond in pattern.GetBonds():\n aid1 = hit_atoms[bond.GetBeginAtomIdx()]\n aid2 = hit_atoms[bond.GetEndAtomIdx()]\n hit_bonds.append(fragment.GetBondBetweenAtoms(aid1,aid2).GetIdx())\n all_hit_atoms.append(hit_atoms)\n all_hit_bonds.append(hit_bonds)\n \n hit_atoms = list(bb_mol.GetSubstructMatch(pattern))\n hit_bonds = []\n for bond in pattern.GetBonds():\n aid1 = hit_atoms[bond.GetBeginAtomIdx()]\n aid2 = hit_atoms[bond.GetEndAtomIdx()]\n hit_bonds.append(bb_mol.GetBondBetweenAtoms(aid1,aid2).GetIdx())\n all_hit_atoms.append(hit_atoms)\n all_hit_bonds.append(hit_bonds)\n \n d = max(fragment.GetNumAtoms(), bb_mol.GetNumAtoms())\n has_linker1 = np.any([a.GetSymbol() == 'Dy' for a in fragment.GetAtoms()])\n if has_linker1:\n d = fragment.GetNumAtoms()\n has_linker2 = np.any([a.GetSymbol() == 'Dy' for a in bb_mol.GetAtoms()])\n if has_linker2:\n d = bb_mol.GetNumAtoms()\n score = match_result.numAtoms / d\n mol_pairs.extend([fragment, bb_mol])\n legends.extend([\"\", f\"fragment {i} and building block {j}: score {score:3.3f}\"])\n \n\nChem.Draw.MolsToGridImage(\n mol_pairs, \n subImgSize=(300,300), \n legends=legends,\n molsPerRow=6,\n highlightAtomLists=all_hit_atoms\n)","metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"execution":{"iopub.status.busy":"2024-04-28T14:12:41.735096Z","iopub.execute_input":"2024-04-28T14:12:41.735490Z","iopub.status.idle":"2024-04-28T14:12:42.008484Z","shell.execute_reply.started":"2024-04-28T14:12:41.735453Z","shell.execute_reply":"2024-04-28T14:12:42.007188Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"I've used that code to process all the binding molecules in train dataset and all the molecules in test dataset. The result of this (time-consuming) computation is collected in the following notebook: https://www.kaggle.com/code/latticetower/belka-is-nuts-from-blocks-to-fragments \n","metadata":{}},{"cell_type":"markdown","source":"# Conclusions\n\n## Building blocks\n1. For protein `sEH` more than 80% of the binding small molecules sharing the same pair of building blocks have 2 or less alternatives (variants of the 3rd building block)\n2. For other 2 proteins more than 80% of the binding small molecules sharing the same pair of building blocks have 5 or less alternatives (variants of the 3rd building block)\n3. The right part of the histograms for binding molecules indicates that there are molecules where binding occurs because of the pair of blocks and doesn't really depend on the choice of the third one.\n\n## Fragment composition\n1. Each building block (in the positive part of the train dataset) corresponds to 1 particular fragment\n2. (In the positive part of the train dataset) all DNA linkers are attached to the fragment which corresponds to the molecule listed in building block 1","metadata":{}},{"cell_type":"code","source":"","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Next steps/TODO\n- view the train data mini-universe\n\nI guess that I might add this to the notebook with precomputed fragments\n\n- check if there are molecules with 2 or 3 equal building blocks\n- check if there are stereoisomers. Does this affect binding?\n","metadata":{"_kg_hide-input":true}},{"cell_type":"code","source":"","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}