{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.10.13","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"none","dataSources":[{"sourceId":67356,"databundleVersionId":8006601,"sourceType":"competition"},{"sourceId":8213585,"sourceType":"datasetVersion","datasetId":4867995}],"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# Data statistics and vizualisation","metadata":{}},{"cell_type":"markdown","source":"The competition aims to predict the binding between 3 proteins and molecules built from 3 building blocks. Leash Bio measured bindings and labeled the data in-house. Proteins are: \n* cEPHX2 (sEH) commonly named “soluble epoxide hydrolase”. \"**Hydrolases are enzymes that catalyze certain chemical reactions, and EPHX2/sEH also hydrolyzes certain phosphate groups. EPHX2/sEH is a potential drug target for high blood pressure and diabetes progression.**\"\n<div>\n<img src=\"https://cdn.rcsb.org/images/structures/3i28_assembly-1.jpeg\" width=\"300\"/>\n</div>\n* BRD4 named as \"bromodomain 4\". \"**Bromodomains play roles in cancer progression and a number of drugs have been discovered to inhibit their activities**\"\n<div>\n<img src=\"https://cdn.rcsb.org/images/structures/7usk_assembly-1.jpeg\" width=\"300\"/>\n</div>\n* ALB (HSA) named as \"serum albumin\" is the \"most common protein in the blood\", \"is used to drive osmotic pressure (to bring fluid back from tissues into blood vessels) and to transport many ligands, hormones, fatty acids, and more\"\n<div>\n<img src=\"https://cdn.rcsb.org/images/structures/5otb_assembly-1.jpeg\" width=\"300\"/>\n</div>\n\nA tested smile molecule is composed by 3 building blocks, named respecitvely buildingblock{i}_smiles for i in {1,2,3}, and in the training set is composed of the 3 building block and a triazine core (cf [here](https://www.kaggle.com/code/chemdatafarmer/scaffold-exploration)). The target is 'binds' and is 1 if a molecule binds with the protein, 0 otherwise.","metadata":{}},{"cell_type":"markdown","source":"**The aim of this notebook is to devle into data to understand its distribution, and its statistical meaning before starting to implement predictive algorithms. Keep in mind it is still WIP**\nWe will thus try to:\n* Understand distributions and representativity of features and target (i.e building blocks)\n* See the influence of building blocks as priors (if a priori some molecules blocks or patterns will react more with a protein)\n* Try to link building blocks and their corresponding chemical meaning (WIP)\n* Check cleaness of data (for example verify if smiles correspond to the stack of building block) (TODO)\n\n*NB: We use plotly as library to plot, since static notebooks don't display plotly graphs, for display only we used image of plots rather than plots*","metadata":{}},{"cell_type":"markdown","source":"## Setup, import of librairies and data","metadata":{}},{"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\n!pip install --quiet rdkit duckdb plotly kaleido global-chem\n\nimport numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\nimport duckdb\n\nimport math\n\nimport random\nfrom collections import Counter\n\nimport plotly.express as px\nimport plotly.graph_objects as go\nfrom IPython.display import Image\nfrom plotly.subplots import make_subplots\n\n# Input data files are available in the read-only \"../input/\" directory\n# For example, running this (by clicking run or pressing Shift+Enter) will list all files under the input directory\n\nfrom rdkit import Chem\nfrom rdkit.Chem import Draw, AllChem\nfrom rdkit import RDLogger\n\nimport os\nfor dirname, _, filenames in os.walk('/kaggle/input'):\n    for filename in filenames:\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\n\ntrain_path = '/kaggle/input/leash-BELKA/train.parquet'\ntest_path = '/kaggle/input/leash-BELKA/test.parquet'\ncon = duckdb.connect()","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2024-05-03T15:16:17.254148Z","iopub.execute_input":"2024-05-03T15:16:17.255097Z","iopub.status.idle":"2024-05-03T15:16:41.610824Z","shell.execute_reply.started":"2024-05-03T15:16:17.255052Z","shell.execute_reply":"2024-05-03T15:16:41.609633Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"n_limit = 300000\n\nsample = con.query(f\"\"\"SELECT *\n                        FROM parquet_scan('{train_path}')\n                        ORDER BY random()\n                        LIMIT {n_limit}\n                       \"\"\").df()\n\nsample_bind = con.query(f\"\"\"SELECT *\n                        FROM parquet_scan('{train_path}')\n                        WHERE binds = 1\n                        ORDER BY random()\n                        LIMIT {n_limit}\n                        \"\"\").df()\n\nsample_no_bind = con.query(f\"\"\"SELECT *\n                        FROM parquet_scan('{train_path}')\n                        WHERE binds = 0\n                        ORDER BY random()\n                        LIMIT {n_limit}\n                        \"\"\").df()\n\nsample_test = con.query(f\"\"\"SELECT *\n                        FROM parquet_scan('{test_path}')\n                        ORDER BY random()\n                        LIMIT {n_limit}\n                        \"\"\").df()\nsample_balanced = pd.concat([sample_bind,sample_no_bind])\n\n","metadata":{"execution":{"iopub.status.busy":"2024-05-03T15:16:41.612979Z","iopub.execute_input":"2024-05-03T15:16:41.613623Z","iopub.status.idle":"2024-05-03T15:18:21.390293Z","shell.execute_reply.started":"2024-05-03T15:16:41.613583Z","shell.execute_reply":"2024-05-03T15:18:21.389099Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Basic statistics","metadata":{}},{"cell_type":"code","source":"nb_rows = con.query(f\"\"\"SELECT COUNT(*) FROM parquet_scan('{train_path}')\"\"\").df().loc[0,\"count_star()\"]\nprint(f\"There are {'{:e}'.format(nb_rows)} rows in the dataset\")","metadata":{"execution":{"iopub.status.busy":"2024-05-03T15:18:21.392071Z","iopub.execute_input":"2024-05-03T15:18:21.392484Z","iopub.status.idle":"2024-05-03T15:18:21.466984Z","shell.execute_reply.started":"2024-05-03T15:18:21.392444Z","shell.execute_reply":"2024-05-03T15:18:21.465854Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"nb_HSA = con.query(f\"\"\"SELECT COUNT(*) FROM parquet_scan('{train_path}') WHERE protein_name = 'HSA' \"\"\").df().loc[0,\"count_star()\"]\nnb_sEH = con.query(f\"\"\"SELECT COUNT(*) FROM parquet_scan('{train_path}') WHERE protein_name = 'sEH' \"\"\").df().loc[0,\"count_star()\"]\nnb_BRD4 = con.query(f\"\"\"SELECT COUNT(*) FROM parquet_scan('{train_path}') WHERE protein_name = 'BRD4' \"\"\").df().loc[0,\"count_star()\"]\n\n\nprint(f\"There are {'{:e}'.format(nb_HSA)} rows for protein HSA in the dataset\")\nprint(f\"There are {'{:e}'.format(nb_sEH)} rows for protein sEH in the dataset\")\nprint(f\"There are {'{:e}'.format(nb_BRD4)} rows for protein BRD4 in the dataset\")","metadata":{"execution":{"iopub.status.busy":"2024-05-03T15:18:21.470135Z","iopub.execute_input":"2024-05-03T15:18:21.470558Z","iopub.status.idle":"2024-05-03T15:18:28.337668Z","shell.execute_reply.started":"2024-05-03T15:18:21.470527Z","shell.execute_reply":"2024-05-03T15:18:28.336327Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We can thus confirm that we have the same number of examples for each protein and it makes sense! In fact each molecule has been tested onto the three proteins","metadata":{}},{"cell_type":"code","source":"nb_bind = con.query(f\"\"\"SELECT COUNT(*) FROM parquet_scan('{train_path}') WHERE binds = 1 \"\"\").df().loc[0,\"count_star()\"]\nnb_non_bind = con.query(f\"\"\"SELECT COUNT(*) FROM parquet_scan('{train_path}') WHERE binds = 0 \"\"\").df().loc[0,\"count_star()\"]\n\nprint(f\"The ratio of 'non binding examples/binding examples' = {int(nb_non_bind/nb_bind)}\")","metadata":{"execution":{"iopub.status.busy":"2024-05-03T15:18:28.339392Z","iopub.execute_input":"2024-05-03T15:18:28.339832Z","iopub.status.idle":"2024-05-03T15:18:31.080851Z","shell.execute_reply.started":"2024-05-03T15:18:28.339774Z","shell.execute_reply":"2024-05-03T15:18:31.079874Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Distributions","metadata":{}},{"cell_type":"code","source":"blocks_lengths = [sample_no_bind[f'buildingblock{i+1}_smiles'].apply(lambda x: len(x)) for i in range(3)]\n\n\nfig = go.Figure()\n\nfor i in range(3):\n  fig.add_trace(go.Box(y=blocks_lengths[i], name=f'buildingblock{i+1}_smiles'))\nfig.update_layout(title=\"Distribution of length of building blocks\")\n#fig.show()\nImage(fig.to_image(format=\"png\", width=600, height=350, scale=2))\n","metadata":{"execution":{"iopub.status.busy":"2024-05-03T15:18:31.082223Z","iopub.execute_input":"2024-05-03T15:18:31.082623Z","iopub.status.idle":"2024-05-03T15:18:31.760425Z","shell.execute_reply.started":"2024-05-03T15:18:31.082585Z","shell.execute_reply":"2024-05-03T15:18:31.759142Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"One can see that building blocks 1 seem to be a bit longer that the two other ones, maybe it would be appropriate to try to know why","metadata":{}},{"cell_type":"markdown","source":"## Balance of dataset","metadata":{}},{"cell_type":"markdown","source":"### Balance of binding by protein","metadata":{}},{"cell_type":"code","source":"# Get representation of each protein by binds count\nfig = px.histogram(sample, x=\"protein_name\", color=\"binds\", barmode=\"group\",title=\"Protein representation by binds in log scale\", log_y=True)\n#fig.show()\nImage(fig.to_image(format=\"png\", width=600, height=350, scale=2))","metadata":{"execution":{"iopub.status.busy":"2024-05-03T15:18:31.762000Z","iopub.execute_input":"2024-05-03T15:18:31.762691Z","iopub.status.idle":"2024-05-03T15:18:35.691829Z","shell.execute_reply.started":"2024-05-03T15:18:31.762658Z","shell.execute_reply":"2024-05-03T15:18:35.690776Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Balance of building blocks by binding","metadata":{}},{"cell_type":"code","source":"buildingblock_nb = [sample_no_bind[f\"buildingblock{i+1}_smiles\"].value_counts().to_list() for i in range(3)]\nbuildingblock_nb += [sample_bind[f\"buildingblock{i+1}_smiles\"].value_counts().to_list() for i in range(3)]\ncumulative_probs = []\n\nfor i in range(6):\n    counts = buildingblock_nb[i]\n    total_count = sum(counts)\n    cumulative_sum = np.cumsum(counts)\n    cumulative_prob = cumulative_sum / total_count\n    cumulative_probs.append(cumulative_prob)\nbindings_label = [\"no_binding\"]*3 + [\"binding\"]*3\nfig = go.Figure()\n\nfor i in range(6):\n    fig.add_trace(go.Scatter(x=np.arange(len(cumulative_probs[i])) + 1, \n                             y=cumulative_probs[i], mode='lines', \n                             name=f'buildingblock{(i%3)+1}_smiles_{bindings_label[i]}',\n                            legendgroup=i%3))\n\n\n\n    \nfig.update_layout(title='Cumulative Distribution of Building Blocks in the training set',\n                  xaxis_title='Building Block Count',\n                  yaxis_title='Cumulative Probability')\n\n#fig.show()\nImage(fig.to_image(format=\"png\", width=600, height=420, scale=2))","metadata":{"execution":{"iopub.status.busy":"2024-05-03T15:18:35.693197Z","iopub.execute_input":"2024-05-03T15:18:35.693771Z","iopub.status.idle":"2024-05-03T15:18:35.946173Z","shell.execute_reply.started":"2024-05-03T15:18:35.693740Z","shell.execute_reply":"2024-05-03T15:18:35.945148Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**From this curve, for binding 0, one can see that building block 1 seem to be more or less linearly spreaded (thus each building block seems to appear the same number of time than others). For binding 1, it is way more logarithmly spreaded that binding 0. It is a bit different from building blocks 2 & 3 that seem to be a bit more logarithmic and thus some blocks might appear way less frequently than others. It can problematic during training**","metadata":{}},{"cell_type":"code","source":"buildingblock_test = [sample_test[f\"buildingblock{i+1}_smiles\"].value_counts().to_list() for i in range(3)]\ncumulative_probs = []\n\nfor i in range(3):\n    counts = buildingblock_nb[i]\n    total_count = sum(counts)\n    cumulative_sum = np.cumsum(counts)\n    cumulative_prob = cumulative_sum / total_count\n    cumulative_probs.append(cumulative_prob)\nbindings_label = [\"no_binding\"]*3 + [\"binding\"]*3\nfig = go.Figure()\n\nfor i in range(3):\n    fig.add_trace(go.Scatter(x=np.arange(len(cumulative_probs[i])) + 1, \n                             y=cumulative_probs[i], mode='lines', \n                             name=f'buildingblock{(i)+1}_smiles'))\n\n\n\n    \nfig.update_layout(title='Cumulative Distribution of Building Blocks in the test set',\n                  xaxis_title='Building Block Count',\n                  yaxis_title='Cumulative Probability')\n\n#fig.show()\nImage(fig.to_image(format=\"png\", width=600, height=420, scale=2))","metadata":{"execution":{"iopub.status.busy":"2024-05-03T15:18:35.947619Z","iopub.execute_input":"2024-05-03T15:18:35.947925Z","iopub.status.idle":"2024-05-03T15:18:36.081118Z","shell.execute_reply.started":"2024-05-03T15:18:35.947900Z","shell.execute_reply":"2024-05-03T15:18:36.080099Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"In the test set, it seems that representativity of molecules is more or less linear at least for building block 1 (like non-binding, and cumulative distribution are more or less the same as in the non-binding case in the train set). It is surely explained by the fact that in the test set there is the same unbalancy between non-binding and binding examples.","metadata":{}},{"cell_type":"markdown","source":"## Overlapping between train and test set","metadata":{}},{"cell_type":"markdown","source":"In this part, we check which building blocks are within the train and test set or both. We see that all the blocks in the training set are within the test set, but not the contrary...\nModeling can be thus tricky since models will have to generalize quite well!","metadata":{}},{"cell_type":"code","source":"fq_builds_interesect = []\nfq_builds_train = []\nfq_builds_test = []\n\n\nfor i in range(1,4):\n\n    buildi_train_set = con.query(f\"\"\"SELECT DISTINCT buildingblock{i}_smiles\n                            FROM parquet_scan('{train_path}')\n                           \"\"\").df()\n\n    buildi_test_set = con.query(f\"\"\"SELECT DISTINCT buildingblock{i}_smiles\n                            FROM parquet_scan('{test_path}')\n                           \"\"\").df()\n\n    intersect_buildi = set(list(buildi_train_set.values.squeeze())).intersection(list(buildi_test_set.values.squeeze()))\n\n    set_train_build_i = set(list(buildi_train_set.values.squeeze()))\n    set_test_build_i = set(list(buildi_test_set.values.squeeze()))\n\n    intersect_buildi = set_train_build_i.intersection(set_test_build_i)\n\n    train_only_buildi = set_train_build_i.difference(intersect_buildi)\n    test_only_buildi = set_test_build_i.difference(intersect_buildi)\n\n    n_buildi_interesect = len(intersect_buildi)\n    n_buildi_test = len(test_only_buildi)\n    n_buildi_train = len(train_only_buildi)\n    n_buildi_total = n_buildi_test + n_buildi_train + n_buildi_interesect\n\n\n    fq_buildi_interesect = n_buildi_interesect/n_buildi_total\n    fq_buildi_train = n_buildi_train/n_buildi_total\n    fq_buildi_test = n_buildi_test/n_buildi_total\n    \n    fq_builds_interesect.append(fq_buildi_interesect)\n    fq_builds_train.append(fq_buildi_train)\n    fq_builds_test.append(fq_buildi_test)","metadata":{"execution":{"iopub.status.busy":"2024-05-03T15:18:36.084411Z","iopub.execute_input":"2024-05-03T15:18:36.084706Z","iopub.status.idle":"2024-05-03T15:18:52.649863Z","shell.execute_reply.started":"2024-05-03T15:18:36.084681Z","shell.execute_reply":"2024-05-03T15:18:52.648843Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import plotly.graph_objects as go\n\n#fig = go.Figure()\ncircles = []\nsubplot_titles = [f\"Overlapping of building blocks {i+1}\" for i in range(3)]\nfig = make_subplots(rows=1, cols=3, subplot_titles=subplot_titles)\nfor i in range(3):\n\n    text = [f\"Train set <br>{int(100*fq_builds_train[i])}%\", \n            f\"Intersection <br>{int(100*fq_builds_interesect[i])}%\",\n            f\"Test set <br>{int(100*fq_builds_test[i])}%\"]\n    # Create scatter trace of text labels\n    fig.add_trace(go.Scatter(\n        x=[1, 1.75, 2.5],\n        y=[1, 1, 1],\n        text=text,\n        mode=\"text\",\n        textfont=dict(\n            color=\"black\",\n            size=18,\n            family=\"Arail\",\n        )\n    ),row=1, col=i+1)\n\n    # Update axes properties\n    fig.update_xaxes(\n        showticklabels=False,\n        showgrid=False,\n        zeroline=False,\n    )\n\n    fig.update_yaxes(\n        showticklabels=False,\n        showgrid=False,\n        zeroline=False,\n    )\n\n\n    # Add circles\n    circles.append(dict(type=\"circle\",\n        line_color=\"gray\", fillcolor=\"green\",\n        x0=0, y0=0, x1=2, y1=2,opacity=0.3,xref=f\"x{i+1}\",yref=f\"y{i+1}\",\n    ))\n\n    circles.append(dict(type=\"circle\",\n        line_color=\"gray\", fillcolor=\"orange\",\n        x0=1.5, y0=0, x1=3.5, y1=2,opacity=0.3,xref=f\"x{i+1}\",yref=f\"y{i+1}\",\n    ))\n\n\n#fig.update_shapes(opacity=0.3, xref=\"x\", yref=\"y\")\nfig.update_layout(\n    shapes=circles)\n\nfig.update_layout(\n    margin=dict(l=20, r=20, b=100),\n    height=600, width=800,\n    plot_bgcolor=\"white\"\n)\n\nfig.update_layout(title=\"Overlapping of building blocks between train and test sets\")\n\nfig.update_layout(\n    autosize=False,\n    width=2000,\n    height=500)\n\n#fig.show()\nImage(fig.to_image(format=\"png\", width=2000, height=500, scale=1))\n","metadata":{"execution":{"iopub.status.busy":"2024-05-03T15:18:52.653677Z","iopub.execute_input":"2024-05-03T15:18:52.656648Z","iopub.status.idle":"2024-05-03T15:18:52.793979Z","shell.execute_reply.started":"2024-05-03T15:18:52.656605Z","shell.execute_reply":"2024-05-03T15:18:52.793171Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"'''\n# TODO in the future, change the size of circles to really represent overlapping percentage\n\nfig = go.Figure()\n\nfreqs = [int(100*fq_build1_train), int(100*fq_build1_interesect),int(100*fq_build1_test)]\ntexts = [\"Train set only <br>\", \"Intersection <br>\", \"Test set only <br>\"]\ntexts = [text+str(freq)+\"%\" if freq>0 else \"\" for text, freq in zip(texts, freqs) ]\n\n\n\nRad_max = 5\n\n\nRad_interesect = Rad_max*fq_build1_interesect\nXY_interesct = [0]*2+[Rad_interesect]*2\n\nRad_train = Rad_max*(fq_build1_interesect+fq_build1_train)\nXY_train = [0]*2+[Rad_train]*2\n\n\nRad_test = Rad_max*(fq_build1_interesect+fq_build1_test)\nXY_test = [0]*2+[Rad_test/(math.sqrt(2))]*2\n\n\n\n# Create scatter trace of text labels\nfig.add_trace(go.Scatter(\n    x=[(XY_train[0]+XY_train[2])/2, (XY_interesct[0]+XY_interesct[2])/2, (XY_test[0]+XY_test[2])/2],\n    y=[(XY_train[1]+XY_train[3])/2, (XY_interesct[1]+XY_interesct[3])/2, (XY_test[1]+XY_test[3])/2],\n    text=texts,\n    mode=\"text\",\n    textfont=dict(\n        color=\"black\",\n        size=18,\n        family=\"Arail\",\n    )\n))\n\n# Update axes properties\nfig.update_xaxes(\n    showticklabels=False,\n    showgrid=False,\n    zeroline=False,\n)\n\nfig.update_yaxes(\n    showticklabels=False,\n    showgrid=False,\n    zeroline=False,\n)\n\n# Add circles\nfig.add_shape(type=\"circle\",\n    line_color=\"gray\", fillcolor=\"green\",\n    x0=XY_train[0], y0=XY_train[1], x1=XY_train[2], y1=XY_train[3]\n)\nfig.add_shape(type=\"circle\",\n    line_color=\"gray\", fillcolor=\"orange\",\n    x0=XY_interesct[0], y0=XY_interesct[1], x1=XY_interesct[2], y1=XY_interesct[3]\n)\n\nfig.add_shape(type=\"circle\",\n   line_color=\"gray\", fillcolor=\"orange\",\n    x0=XY_test[0], y0=XY_test[1], x1=XY_test[2], y1=XY_test[3]\n)\n\nfig.update_shapes(opacity=0.3, xref=\"x\", yref=\"y\")\n\nfig.update_layout(\n    margin=dict(l=20, r=20, b=100),\n    height=600, width=800,\n    plot_bgcolor=\"white\"\n)\n\nfig.show()\n\n'''","metadata":{"execution":{"iopub.status.busy":"2024-05-03T15:18:52.795562Z","iopub.execute_input":"2024-05-03T15:18:52.796241Z","iopub.status.idle":"2024-05-03T15:18:52.807011Z","shell.execute_reply.started":"2024-05-03T15:18:52.796202Z","shell.execute_reply":"2024-05-03T15:18:52.805905Z"},"jupyter":{"source_hidden":true},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Visualizing molecules","metadata":{}},{"cell_type":"code","source":"###Visualize some of the original data\n#print(sample['molecule_smiles'][0:4].tolist())\nDraw.MolsToGridImage([Chem.MolFromSmiles(x) for x in sample['molecule_smiles'][0:4]], molsPerRow=4, subImgSize=(400,300))","metadata":{"execution":{"iopub.status.busy":"2024-05-03T15:18:52.808736Z","iopub.execute_input":"2024-05-03T15:18:52.810139Z","iopub.status.idle":"2024-05-03T15:18:52.877031Z","shell.execute_reply.started":"2024-05-03T15:18:52.810101Z","shell.execute_reply":"2024-05-03T15:18:52.875827Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"From these examples, we can see the triazine core with the three building blocks around","metadata":{}},{"cell_type":"markdown","source":"# WIP","metadata":{}},{"cell_type":"code","source":"def compute_ngrams(text, n):\n    ngrams = [text[i:i+n] for i in range(len(text) - n + 1)]\n    return ngrams\n\ndef most_common_ngrams(text, n, top_k):\n    ngrams = compute_ngrams(text, n)\n    ngram_counter = Counter(ngrams)\n    return ngram_counter.most_common(top_k)\n\ndef create_fig_n_grams(df, n, top_k):\n    \n    data1 = most_common_ngrams(\"\".join(df[df.binds == 0][\"molecule_smiles\"]),n,10)\n    data2 = most_common_ngrams(\"\".join(df[df.binds == 1][\"molecule_smiles\"]),n,10)\n\n    # Extract x and y values for each list\n    x_values = set([item[0] for item in data1] + [item[0] for item in data2])\n    y_values1 = {item[0]: item[1] for item in data1}\n    y_values2 = {item[0]: item[1] for item in data2}\n\n    # Fill in missing values with zeros\n    for x in x_values:\n        if x not in y_values1:\n            y_values1[x] = 0\n        if x not in y_values2:\n            y_values2[x] = 0\n\n    # Sort x_values alphabetically for better visualization\n    x_values = sorted(x_values)\n\n    # Create traces for each list\n    trace1 = go.Bar(x=x_values, y=[y_values1[x] for x in x_values], name='Binds = 0', marker=dict(color='rgb(55, 83, 109)'))\n    trace2 = go.Bar(x=x_values, y=[y_values2[x] for x in x_values], name='Binds = 1', marker=dict(color='rgb(26, 118, 255)'))\n\n    # Define layout\n    layout = go.Layout(\n        title=f'Most Common {str(n)}-grams',\n        xaxis=dict(title='n-grams'),\n        yaxis=dict(title='Frequency'),\n        barmode='group'\n    )\n\n    # Create figure\n    fig = go.Figure(data=[trace1, trace2], layout=layout)\n    \n    return fig\n\nn = 3\ntop_k = 10\nfig1 = create_fig_n_grams(sample_balanced,n,top_k)\n\n# Show the figure\n#fig1.show()\nImage(fig1.to_image(format=\"png\", width=600, height=420, scale=2))\n","metadata":{"jupyter":{"source_hidden":true,"outputs_hidden":true},"execution":{"iopub.status.busy":"2024-05-03T15:18:52.878566Z","iopub.execute_input":"2024-05-03T15:18:52.879579Z","iopub.status.idle":"2024-05-03T15:19:07.366976Z","shell.execute_reply.started":"2024-05-03T15:18:52.879539Z","shell.execute_reply":"2024-05-03T15:19:07.366210Z"},"collapsed":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"n = 5\ntop_k = 10\nfig2 = create_fig_n_grams(sample_balanced,n,top_k)\n\n# Show the figure\n#fig2.show()\nImage(fig2.to_image(format=\"png\", width=600, height=420, scale=2))\n","metadata":{"jupyter":{"source_hidden":true,"outputs_hidden":true},"execution":{"iopub.status.busy":"2024-05-03T15:19:07.368282Z","iopub.execute_input":"2024-05-03T15:19:07.368831Z","iopub.status.idle":"2024-05-03T15:19:22.332007Z","shell.execute_reply.started":"2024-05-03T15:19:07.368802Z","shell.execute_reply":"2024-05-03T15:19:22.330951Z"},"collapsed":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# 2. WIP: Method using an external db of usual molecular components ","metadata":{}},{"cell_type":"markdown","source":"**This part is still WIP. It aims to identify classical chemical group in the molecules and see if there is some unbalancy and some links with binding or not with the proteins**","metadata":{}},{"cell_type":"code","source":"simple_compounds = pd.read_csv('/kaggle/input/pubchem-simple-compound/PubChem_simple_compound_list.csv').sort_values(by='mw')\nsimple_compounds.head(2)","metadata":{"execution":{"iopub.status.busy":"2024-05-03T15:19:22.332897Z","iopub.execute_input":"2024-05-03T15:19:22.333170Z","iopub.status.idle":"2024-05-03T15:19:22.488697Z","shell.execute_reply.started":"2024-05-03T15:19:22.333148Z","shell.execute_reply":"2024-05-03T15:19:22.487660Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def count_occurrences(strings_list, target_strings):\n    occurrences = {target: 0 for target in target_strings}  # Initialize counts for each target string\n    \n    for string in strings_list:\n        for target in target_strings:\n            occurrences[target] += string.count(target)\n    \n    return occurrences\n\nbasic_compounds = simple_compounds.isosmiles.tolist()[:300]\n\ncompound_occurrences_binds0 = count_occurrences(sample_balanced[sample_balanced.binds == 0].molecule_smiles.tolist(),basic_compounds)\ncompound_occurrences_binds1 = count_occurrences(sample_balanced[sample_balanced.binds == 1].molecule_smiles.tolist(),basic_compounds)","metadata":{"execution":{"iopub.status.busy":"2024-05-03T15:19:22.490020Z","iopub.execute_input":"2024-05-03T15:19:22.490315Z","iopub.status.idle":"2024-05-03T15:20:27.662312Z","shell.execute_reply.started":"2024-05-03T15:19:22.490290Z","shell.execute_reply":"2024-05-03T15:20:27.661173Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"compound_occurrences_binds0_filtered = {key: val for key,val in  sorted(compound_occurrences_binds0.items(), key=lambda item: item[1],reverse=True) if val > 0}\ncompound_occurrences_binds1_filtered = {key: val for key,val in  sorted(compound_occurrences_binds1.items(), key=lambda item: item[1],reverse=True) if val > 0}","metadata":{"execution":{"iopub.status.busy":"2024-05-03T15:20:27.663600Z","iopub.execute_input":"2024-05-03T15:20:27.664055Z","iopub.status.idle":"2024-05-03T15:20:27.670521Z","shell.execute_reply.started":"2024-05-03T15:20:27.664027Z","shell.execute_reply":"2024-05-03T15:20:27.669738Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def plot_occurrences(dict1, dict2):\n    keys = list(dict1.keys())\n    values1 = list(dict1.values())\n    values2 = list(dict2.values())\n\n    trace1 = go.Bar(x=keys, y=values1, name='Binds 0')\n    trace2 = go.Bar(x=keys, y=values2, name='Binds 1')\n\n    layout = go.Layout(\n        title='Occurrences of Strings',\n        xaxis=dict(title='Strings'),\n        yaxis=dict(title='Occurrences', type='log'),  # Setting y-axis to log scale\n        barmode='group'\n    )\n\n    fig = go.Figure(data=[trace1, trace2], layout=layout)\n    fig.show()\n    Image(fig.to_image(format=\"png\", width=600, height=420, scale=2))\n    \nplot_occurrences(compound_occurrences_binds0_filtered,compound_occurrences_binds1_filtered)","metadata":{"execution":{"iopub.status.busy":"2024-05-03T15:20:27.671699Z","iopub.execute_input":"2024-05-03T15:20:27.672381Z","iopub.status.idle":"2024-05-03T15:20:27.883360Z","shell.execute_reply.started":"2024-05-03T15:20:27.672347Z","shell.execute_reply":"2024-05-03T15:20:27.882190Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"con.close()","metadata":{"execution":{"iopub.status.busy":"2024-05-03T15:20:27.885420Z","iopub.execute_input":"2024-05-03T15:20:27.885819Z","iopub.status.idle":"2024-05-03T15:20:27.915211Z","shell.execute_reply.started":"2024-05-03T15:20:27.885769Z","shell.execute_reply":"2024-05-03T15:20:27.914288Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}