{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"import pandas as pd\nimport numpy as np","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2023-10-31T17:17:25.515908Z","iopub.execute_input":"2023-10-31T17:17:25.517165Z","iopub.status.idle":"2023-10-31T17:17:25.522154Z","shell.execute_reply.started":"2023-10-31T17:17:25.517123Z","shell.execute_reply":"2023-10-31T17:17:25.520891Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Exploring all unique SMILES","metadata":{}},{"cell_type":"code","source":"def find_unique_vals(values):\n    unique_vals = []\n    \n    for i in values:\n        if i in unique_vals:\n            pass\n        else:\n            unique_vals.append(i)\n    \n    return unique_vals","metadata":{"execution":{"iopub.status.busy":"2023-10-31T17:17:25.526292Z","iopub.execute_input":"2023-10-31T17:17:25.526692Z","iopub.status.idle":"2023-10-31T17:17:25.536758Z","shell.execute_reply.started":"2023-10-31T17:17:25.526661Z","shell.execute_reply":"2023-10-31T17:17:25.535799Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df = pd.read_parquet('/kaggle/input/open-problems-single-cell-perturbations/de_train.parquet')\nsmiles_value = df['SMILES'].values\n\nfor i in find_unique_vals(smiles_value):\n    print(i)","metadata":{"execution":{"iopub.status.busy":"2023-10-31T17:17:25.538832Z","iopub.execute_input":"2023-10-31T17:17:25.539936Z","iopub.status.idle":"2023-10-31T17:17:26.914246Z","shell.execute_reply.started":"2023-10-31T17:17:25.539891Z","shell.execute_reply":"2023-10-31T17:17:26.912804Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Converting into 2D array","metadata":{}},{"cell_type":"code","source":"map_array = np.zeros((256,256)) \n# Making the 256,256 array, 256 Size was chosen because after converting all the x,y points from\n# Negative to positive, The minimum value was 0 and maximum was 200, So 256 fits comfortably","metadata":{"execution":{"iopub.status.busy":"2023-10-31T17:17:26.915935Z","iopub.execute_input":"2023-10-31T17:17:26.916390Z","iopub.status.idle":"2023-10-31T17:17:26.921886Z","shell.execute_reply.started":"2023-10-31T17:17:26.916348Z","shell.execute_reply":"2023-10-31T17:17:26.920728Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Installing Pysmiles","metadata":{}},{"cell_type":"code","source":"!pip install pysmiles","metadata":{"execution":{"iopub.status.busy":"2023-10-31T17:17:26.923613Z","iopub.execute_input":"2023-10-31T17:17:26.924484Z","iopub.status.idle":"2023-10-31T17:17:40.663037Z","shell.execute_reply.started":"2023-10-31T17:17:26.924452Z","shell.execute_reply":"2023-10-31T17:17:40.661702Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from pysmiles import read_smiles\nimport matplotlib.pyplot as plt\nimport networkx as nx","metadata":{"execution":{"iopub.status.busy":"2023-10-31T17:17:40.666747Z","iopub.execute_input":"2023-10-31T17:17:40.667174Z","iopub.status.idle":"2023-10-31T17:17:40.673028Z","shell.execute_reply.started":"2023-10-31T17:17:40.667140Z","shell.execute_reply":"2023-10-31T17:17:40.671819Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"mol = read_smiles(\"COC(=O)N(C)c1c(N)nc(-c2nn(Cc3ccccc3F)c3ncccc23)nc1N\") \n# Reading SMILES code, this returns a NetworkX graph which we can use to calculate the position of atoms\n# in a 2D plane, Atoms are represented as Nodes and Bonds are represented as Edges\n\n\npos = nx.spring_layout(mol)\n# spring_layout is an algorithm which will return x, y points of all the nodes in the node graph\n# pos is a dictionary where keys are nodes and values are x, y coordinates","metadata":{"execution":{"iopub.status.busy":"2023-10-31T17:17:40.674653Z","iopub.execute_input":"2023-10-31T17:17:40.675132Z","iopub.status.idle":"2023-10-31T17:17:40.707854Z","shell.execute_reply.started":"2023-10-31T17:17:40.675091Z","shell.execute_reply":"2023-10-31T17:17:40.706972Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"What mol looks like: ","metadata":{}},{"cell_type":"code","source":"nx.draw(mol, with_labels=True, node_color='lightblue', font_weight='bold')\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-10-31T17:17:40.709323Z","iopub.execute_input":"2023-10-31T17:17:40.709860Z","iopub.status.idle":"2023-10-31T17:17:40.993954Z","shell.execute_reply.started":"2023-10-31T17:17:40.709828Z","shell.execute_reply":"2023-10-31T17:17:40.992704Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Adding Atom number to 2D array","metadata":{}},{"cell_type":"code","source":"# Storing all Node locations in a dictionary\nnode_location = {}\n\n\nfor node in mol.nodes():\n    x,y = pos[node]\n    x,y = (x*100) + 100, (y*100) + 100 # Multiplying by 100 because the coordinates are too small \n    # Adding 100 because the lowest coordinate of both axes are -100 (After multiplying by 100)\n    x,y = int(x), int(y)\n    # Converting into int so that we can easily put it into the array\n    \n    node_location[node] = [[x,y], mol.nodes[node][\"element\"]]\n    # We are storing x, y coordinates as well as the element notation, Ex: C, N, O etc","metadata":{"execution":{"iopub.status.busy":"2023-10-31T17:17:40.995689Z","iopub.execute_input":"2023-10-31T17:17:40.996081Z","iopub.status.idle":"2023-10-31T17:17:41.003197Z","shell.execute_reply.started":"2023-10-31T17:17:40.996050Z","shell.execute_reply":"2023-10-31T17:17:41.002021Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Now we can add this to our map array","metadata":{}},{"cell_type":"code","source":"taken_points = []\n# This array is to there to aid in the future algorithm where we have to map out bonds ourseleves\n# Since NetworkX graph edges (Bonds) are represented as Ideal straight lines so they don't have\n# x, y points","metadata":{"execution":{"iopub.status.busy":"2023-10-31T17:17:41.004555Z","iopub.execute_input":"2023-10-31T17:17:41.004898Z","iopub.status.idle":"2023-10-31T17:17:41.015561Z","shell.execute_reply.started":"2023-10-31T17:17:41.004868Z","shell.execute_reply":"2023-10-31T17:17:41.014535Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from PyAstronomy import pyasl\n\nan = pyasl.AtomicNo()\n\n# This helps us to convert Element notation to Atomic Number","metadata":{"execution":{"iopub.status.busy":"2023-10-31T17:17:41.016816Z","iopub.execute_input":"2023-10-31T17:17:41.017122Z","iopub.status.idle":"2023-10-31T17:17:41.029976Z","shell.execute_reply.started":"2023-10-31T17:17:41.017095Z","shell.execute_reply":"2023-10-31T17:17:41.029167Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"for i in node_location.keys():\n    x = node_location[i][0][0]\n    y = node_location[i][0][1]\n    \n    \n    atomic_number = an.getAtomicNo(node_location[i][1])\n    \n    map_array[x][y] = atomic_number\n    \n    taken_points.append([x,y])","metadata":{"execution":{"iopub.status.busy":"2023-10-31T17:17:41.030861Z","iopub.execute_input":"2023-10-31T17:17:41.031179Z","iopub.status.idle":"2023-10-31T17:17:41.040601Z","shell.execute_reply.started":"2023-10-31T17:17:41.031153Z","shell.execute_reply":"2023-10-31T17:17:41.039740Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Plotting the graph","metadata":{}},{"cell_type":"code","source":"plt.imshow(map_array, cmap='viridis')\nplt.colorbar()\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-10-31T17:17:41.041861Z","iopub.execute_input":"2023-10-31T17:17:41.042175Z","iopub.status.idle":"2023-10-31T17:17:41.420843Z","shell.execute_reply.started":"2023-10-31T17:17:41.042148Z","shell.execute_reply":"2023-10-31T17:17:41.419697Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### We have added Atomic numbers of the atoms involved in the SMILES notation to our array","metadata":{}},{"cell_type":"markdown","source":"## Adding Bonds to 2D array","metadata":{}},{"cell_type":"markdown","source":"#### Since edges in NetworkX graph are ideal lines, they don't have a set of x, y coordinates so we will have to make our own algorithm for it","metadata":{}},{"cell_type":"markdown","source":"##### Chosen Algorithm\n\n1) Loop through each edge\n\n2) Target x,y are the coordinates of the atom which the bond goes to (second place in array)\n\n3) We start with coordinates of the atom itself, We get all adjacent points (Or nearby points) which can only be 8 Points in a 2D plane\n\n4) Calculate the point out of the 8 points which has the least distance to the target x, y points\n\n5) Choose that point and add the value of the bond in it\n\n6) Repeat until there's no available nearby points or the target x, y points are in nearby points\n\n\nNote: Available points are those which aren't taken by an existing bond or atom, Nearby points are all nearby points regardless of their availability","metadata":{}},{"cell_type":"code","source":"def get_available_nearby_points(x, y, limit):\n    return_array = []\n    \n    point_1 = [x + 1, y]\n    point_2 = [x - 1, y]\n    \n    point_3 = [x, y + 1]\n    point_4 = [x, y - 1]\n    \n    point_5 = [x - 1, y + 1]\n    point_6 = [x + 1, y + 1]\n    \n    point_7 = [x - 1, y - 1]\n    point_8 = [x + 1, y - 1]\n    \n    if point_1[0] >= 0 and point_1[0] <= limit and point_1[1] >= 0 and point_1[1] <= limit:\n        if point_1 not in taken_points:\n            return_array.append(point_1)\n    \n\n    if point_2[0] >= 0 and point_2[0] <= limit and point_2[1] >= 0 and point_2[1] <= limit:\n        if point_2 not in taken_points:\n            return_array.append(point_2)\n\n\n    if point_3[0] >= 0 and point_3[0] <= limit and point_3[1] >= 0 and point_3[1] <= limit:\n        if point_3 not in taken_points:\n            return_array.append(point_3)\n    \n\n    if point_4[0] >= 0 and point_4[0] <= limit and point_4[1] >= 0 and point_4[1] <= limit:\n        if point_4 not in taken_points:\n            return_array.append(point_4)\n\n    \n    if point_5[0] >= 0 and point_5[0] <= limit and point_5[1] >= 0 and point_5[1] <= limit:\n        if point_5 not in taken_points:\n            return_array.append(point_5)\n    \n\n    if point_6[0] >= 0 and point_6[0] <= limit and point_6[1] >= 0 and point_6[1] <= limit:\n        if point_6 not in taken_points:\n            return_array.append(point_6)\n    \n\n    if point_7[0] >= 0 and point_7[0] <= limit and point_7[1] >= 0 and point_7[1] <= limit:\n        if point_7 not in taken_points:\n            return_array.append(point_7)\n    \n\n    if point_8[0] >= 0 and point_8[0] <= limit and point_8[1] >= 0 and point_8[1] <= limit:\n        if point_8 not in taken_points:\n            return_array.append(point_8)\n    \n\n    return return_array\n\n\ndef get_nearby_points(x, y, limit):\n    return_array = []\n    \n    point_1 = [x + 1, y]\n    point_2 = [x - 1, y]\n    \n    point_3 = [x, y + 1]\n    point_4 = [x, y - 1]\n    \n    point_5 = [x - 1, y + 1]\n    point_6 = [x + 1, y + 1]\n    \n    point_7 = [x - 1, y - 1]\n    point_8 = [x + 1, y - 1]\n    \n    if point_1[0] >= 0 and point_1[0] <= limit and point_1[1] >= 0 and point_1[1] <= limit:\n        return_array.append(point_1)\n    \n\n    if point_2[0] >= 0 and point_2[0] <= limit and point_2[1] >= 0 and point_2[1] <= limit:\n        return_array.append(point_2)\n\n\n    if point_3[0] >= 0 and point_3[0] <= limit and point_3[1] >= 0 and point_3[1] <= limit:\n        return_array.append(point_3)\n    \n\n    if point_4[0] >= 0 and point_4[0] <= limit and point_4[1] >= 0 and point_4[1] <= limit:\n        return_array.append(point_4)\n\n    \n    if point_5[0] >= 0 and point_5[0] <= limit and point_5[1] >= 0 and point_5[1] <= limit:\n        return_array.append(point_5)\n    \n\n    if point_6[0] >= 0 and point_6[0] <= limit and point_6[1] >= 0 and point_6[1] <= limit:\n        return_array.append(point_6)\n    \n\n    if point_7[0] >= 0 and point_7[0] <= limit and point_7[1] >= 0 and point_7[1] <= limit:\n        return_array.append(point_7)\n    \n\n    if point_8[0] >= 0 and point_8[0] <= limit and point_8[1] >= 0 and point_8[1] <= limit:\n        return_array.append(point_8)\n    \n\n    return return_array\n    ","metadata":{"execution":{"iopub.status.busy":"2023-10-31T17:17:41.423452Z","iopub.execute_input":"2023-10-31T17:17:41.423832Z","iopub.status.idle":"2023-10-31T17:17:41.593422Z","shell.execute_reply.started":"2023-10-31T17:17:41.423762Z","shell.execute_reply":"2023-10-31T17:17:41.592085Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import math\n\n\nfor edge in mol.edges():\n    current_atom = edge[0]\n    target_atom = edge[1]\n    \n    target_x,target_y = node_location[target_atom][0]\n    \n    current_atom_x, current_atom_y = node_location[current_atom][0]\n    \n    available_nearby_points = get_available_nearby_points(current_atom_x, current_atom_y, 256)\n    nearby_points = get_nearby_points(current_atom_x, current_atom_y, 256)\n    \n    \n    while [target_x,target_y] not in nearby_points:\n        least_distance = None\n        selected_point = None\n        \n        if len(available_nearby_points) == 0:\n            break\n        \n        for i in available_nearby_points:\n            distance = math.sqrt( ((target_x - i[0]) ** 2) + ((target_y - i[1]) ** 2) )\n            if least_distance == None:\n                least_distance = distance\n                selected_point = i\n            \n            if least_distance > distance:\n                least_distance = distance\n                selected_point = i\n        \n        \n        \n        taken_points.append(selected_point)\n        map_array[selected_point[0]][selected_point[1]] = int(mol.edges[edge][\"order\"]) / 3\n        # We are dividing all the order values by 3 because the highest value is 3 for a carbon atom\n        # There might be able to exist a bond order of 4 for carbon but that's not a normal organic\n        # Substance which you can use\n        # If we just use bond values as 1, 1.5 or 2 etc then it will be equivalent to the neural network\n        # Thinking the bonds are just hydrogen atoms which we don't want\n  \n        \n        available_nearby_points = get_available_nearby_points(selected_point[0], selected_point[1], 256)\n        nearby_points = get_nearby_points(selected_point[0], selected_point[1], 256)\n        \n    ","metadata":{"execution":{"iopub.status.busy":"2023-10-31T17:17:41.597451Z","iopub.execute_input":"2023-10-31T17:17:41.597930Z","iopub.status.idle":"2023-10-31T17:17:41.648901Z","shell.execute_reply.started":"2023-10-31T17:17:41.597897Z","shell.execute_reply":"2023-10-31T17:17:41.647689Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Plotting graph to show the bonds","metadata":{}},{"cell_type":"code","source":"plt.imshow(map_array, cmap='viridis')\nplt.colorbar()\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-10-31T17:17:41.650349Z","iopub.execute_input":"2023-10-31T17:17:41.651008Z","iopub.status.idle":"2023-10-31T17:17:42.031150Z","shell.execute_reply.started":"2023-10-31T17:17:41.650967Z","shell.execute_reply":"2023-10-31T17:17:42.029909Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Applying max pooling to reduce the array size to 64,64","metadata":{}},{"cell_type":"code","source":"def max_pooling_2d(arr, size):\n    rows, cols = arr.shape\n    pooled = np.zeros((rows // size, cols // size))\n    for i in range(0, rows, size):\n        for j in range(0, cols, size):\n            pooled[i // size, j // size] = np.max(arr[i:i+size, j:j+size])\n    return pooled","metadata":{"execution":{"iopub.status.busy":"2023-10-31T17:17:42.032550Z","iopub.execute_input":"2023-10-31T17:17:42.032902Z","iopub.status.idle":"2023-10-31T17:17:42.042280Z","shell.execute_reply.started":"2023-10-31T17:17:42.032872Z","shell.execute_reply":"2023-10-31T17:17:42.040785Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"map_array = max_pooling_2d(map_array, 4)","metadata":{"execution":{"iopub.status.busy":"2023-10-31T17:17:42.043811Z","iopub.execute_input":"2023-10-31T17:17:42.044459Z","iopub.status.idle":"2023-10-31T17:17:42.088278Z","shell.execute_reply.started":"2023-10-31T17:17:42.044416Z","shell.execute_reply":"2023-10-31T17:17:42.087126Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.imshow(map_array, cmap='viridis')\nplt.colorbar()\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-10-31T17:17:42.089692Z","iopub.execute_input":"2023-10-31T17:17:42.090579Z","iopub.status.idle":"2023-10-31T17:17:42.459946Z","shell.execute_reply.started":"2023-10-31T17:17:42.090546Z","shell.execute_reply":"2023-10-31T17:17:42.458845Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(map_array.shape)","metadata":{"execution":{"iopub.status.busy":"2023-10-31T17:18:53.149012Z","iopub.execute_input":"2023-10-31T17:18:53.149830Z","iopub.status.idle":"2023-10-31T17:18:53.154812Z","shell.execute_reply.started":"2023-10-31T17:18:53.149771Z","shell.execute_reply":"2023-10-31T17:18:53.153689Z"},"trusted":true},"execution_count":null,"outputs":[]}]}