{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.10.13","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"none","dataSources":[{"sourceId":67356,"databundleVersionId":8006601,"sourceType":"competition"},{"sourceId":8196925,"sourceType":"datasetVersion","datasetId":4855347}],"dockerImageVersionId":30698,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"!pip install rdkit","metadata":{"execution":{"iopub.status.busy":"2024-06-26T19:04:18.157793Z","iopub.execute_input":"2024-06-26T19:04:18.158987Z","iopub.status.idle":"2024-06-26T19:04:33.127804Z","shell.execute_reply.started":"2024-06-26T19:04:18.158942Z","shell.execute_reply":"2024-06-26T19:04:33.126414Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"!pip install seaborn==0.13.2","metadata":{"execution":{"iopub.status.busy":"2024-06-26T19:04:33.130353Z","iopub.execute_input":"2024-06-26T19:04:33.130731Z","iopub.status.idle":"2024-06-26T19:04:45.891432Z","shell.execute_reply.started":"2024-06-26T19:04:33.130697Z","shell.execute_reply":"2024-06-26T19:04:45.889909Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"!pip install duckdb","metadata":{"execution":{"iopub.status.busy":"2024-06-26T19:04:45.893401Z","iopub.execute_input":"2024-06-26T19:04:45.893834Z","iopub.status.idle":"2024-06-26T19:04:58.623941Z","shell.execute_reply.started":"2024-06-26T19:04:45.893782Z","shell.execute_reply":"2024-06-26T19:04:58.622422Z"},"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 os\nimport gc\nimport numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\nimport matplotlib.pyplot as plt\nimport seaborn as sns\nimport lightgbm as lgb\nfrom sklearn.model_selection import train_test_split\nfrom sklearn.metrics import average_precision_score\n\nimport pandas as pd\nimport rdkit\nfrom rdkit import Chem \nfrom rdkit.Chem import (rdMolDescriptors, Descriptors, \n                        rdmolfiles, rdFingerprintGenerator)\nfrom sklearn.metrics.pairwise import cosine_similarity\nfrom sklearn.preprocessing import StandardScaler\n\nimport duckdb\n\noriginal_data_dir = \"/kaggle/input/leash-BELKA\"\ntrain_file = \"train.parquet\"\ntest_file = \"test.csv\"\nHSA_pdb_filename  = \"/kaggle/input/protein-pdb-data/ALB-HSA.pdb\"\nBRD4_pdb_filename = \"/kaggle/input/protein-pdb-data/BRD4.pdb\"\nsEH_pdb_filename  = \"/kaggle/input/protein-pdb-data/sEH.pdb\"\n\nchemical_radious = 2\nchemical_nbits = 1024\nchunk_size = 100_000      # Incremental Chunks of data\ndata_cap = 200_000        # Total data for each binding option (i.e *2)\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\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","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2024-06-26T19:04:58.627593Z","iopub.execute_input":"2024-06-26T19:04:58.628153Z","iopub.status.idle":"2024-06-26T19:05:02.454530Z","shell.execute_reply.started":"2024-06-26T19:04:58.628103Z","shell.execute_reply":"2024-06-26T19:05:02.452674Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"##### Thanks to Mehra Kazeminia, et al. and @andrewdblevins for the following code segment which elegently reads a balanced dataset from the original data:\n\nhttps://www.kaggle.com/code/mehrankazeminia/p-1-6-belka-eda-data-separation","metadata":{}},{"cell_type":"code","source":"train_path = os.path.join(original_data_dir, train_file)\ncon = duckdb.connect()","metadata":{"execution":{"iopub.status.busy":"2024-06-26T19:05:02.456503Z","iopub.execute_input":"2024-06-26T19:05:02.457360Z","iopub.status.idle":"2024-06-26T19:05:02.474241Z","shell.execute_reply.started":"2024-06-26T19:05:02.457310Z","shell.execute_reply":"2024-06-26T19:05:02.472953Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"raw","source":"train_df = con.query(f\"\"\"(SELECT *\n                            FROM parquet_scan('{train_path}')\n                            WHERE binds = 0 AND protein_name = 'HSA'\n                            ORDER BY random()\n                            LIMIT {rows})\n                            UNION ALL\n                            (SELECT *\n                            FROM parquet_scan('{train_path}')\n                            WHERE binds = 1 AND protein_name = 'HSA'\n                            ORDER BY random()\n                            LIMIT {rows})\n                            \n                            UNION ALL\n                            (SELECT *\n                            FROM parquet_scan('{train_path}')\n                            WHERE binds = 0 AND protein_name = 'BRD4'\n                            ORDER BY random()\n                            LIMIT {rows})\n                            UNION ALL\n                            (SELECT *\n                            FROM parquet_scan('{train_path}')\n                            WHERE binds = 1 AND protein_name = 'BRD4'\n                            ORDER BY random()\n                            LIMIT {rows})\n                            \n                            UNION ALL\n                            (SELECT *\n                            FROM parquet_scan('{train_path}')\n                            WHERE binds = 0 AND protein_name = 'sEH'\n                            ORDER BY random()\n                            LIMIT {rows})\n                            UNION ALL\n                            (SELECT *\n                            FROM parquet_scan('{train_path}')\n                            WHERE binds = 1 AND protein_name = 'sEH'\n                            ORDER BY random()\n                            LIMIT {rows})\"\"\").df()\n\ncon.close()","metadata":{"execution":{"iopub.status.busy":"2024-06-26T16:20:00.050494Z","iopub.execute_input":"2024-06-26T16:20:00.051208Z","iopub.status.idle":"2024-06-26T16:22:35.493323Z","shell.execute_reply.started":"2024-06-26T16:20:00.051133Z","shell.execute_reply":"2024-06-26T16:22:35.489218Z"}}},{"cell_type":"raw","source":"train_df[['protein_name', 'binds', 'molecule_smiles']].groupby(by=['protein_name', 'binds']).count()","metadata":{"execution":{"iopub.status.busy":"2024-06-26T16:22:35.502696Z","iopub.execute_input":"2024-06-26T16:22:35.503917Z","iopub.status.idle":"2024-06-26T16:22:37.027283Z","shell.execute_reply.started":"2024-06-26T16:22:35.503740Z","shell.execute_reply":"2024-06-26T16:22:37.023788Z"}}},{"cell_type":"raw","source":"train_df.to_csv(os.path.join('train_BELKA_balanced.csv'), index=False)","metadata":{"execution":{"iopub.status.busy":"2024-06-26T16:22:37.030828Z","iopub.execute_input":"2024-06-26T16:22:37.031910Z","iopub.status.idle":"2024-06-26T16:23:08.728846Z","shell.execute_reply.started":"2024-06-26T16:22:37.031814Z","shell.execute_reply":"2024-06-26T16:23:08.725321Z"}}},{"cell_type":"raw","source":"train_df.shape, train_df.columns","metadata":{"execution":{"iopub.status.busy":"2024-06-26T16:23:08.734545Z","iopub.execute_input":"2024-06-26T16:23:08.735234Z","iopub.status.idle":"2024-06-26T16:23:08.748663Z","shell.execute_reply.started":"2024-06-26T16:23:08.735169Z","shell.execute_reply":"2024-06-26T16:23:08.747090Z"}}},{"cell_type":"raw","source":"plt.gca().set_facecolor('lightgray')\ntrain_df['binds'].value_counts(normalize=True).plot(kind='barh', figsize=(12,1), color=['darkcyan','red'])\n\npd.DataFrame(data= {'Number': train_df['binds'].value_counts(), \n                    'Percent': train_df['binds'].value_counts(normalize=True)})","metadata":{"execution":{"iopub.status.busy":"2024-06-26T16:23:08.750510Z","iopub.execute_input":"2024-06-26T16:23:08.751031Z","iopub.status.idle":"2024-06-26T16:23:09.142640Z","shell.execute_reply.started":"2024-06-26T16:23:08.750988Z","shell.execute_reply":"2024-06-26T16:23:09.141347Z"}}},{"cell_type":"raw","source":"plt.gca().set_facecolor('lightgray')\ntrain_df['protein_name'].value_counts(normalize=True).plot(kind='barh', figsize=(12,1.5), color=['darkcyan','red','violet'])\n\npd.DataFrame(data= {'Number': train_df['protein_name'].value_counts(), \n                    'Percent': train_df['protein_name'].value_counts(normalize=True)})","metadata":{"execution":{"iopub.status.busy":"2024-06-26T16:23:09.144151Z","iopub.execute_input":"2024-06-26T16:23:09.144504Z","iopub.status.idle":"2024-06-26T16:23:10.787808Z","shell.execute_reply.started":"2024-06-26T16:23:09.144474Z","shell.execute_reply":"2024-06-26T16:23:10.786402Z"}}},{"cell_type":"raw","source":"plt.gca().set_facecolor('lightgray')\ntrain_df.groupby('protein_name')['binds'].value_counts(normalize=True).plot(kind='barh', figsize=(12,3), color=['darkcyan','red'])\n\npd.DataFrame(data= {'Number': train_df.groupby('protein_name')['binds'].value_counts(), \n                    'Percent': train_df.groupby('protein_name')['binds'].value_counts(normalize=True)})","metadata":{"execution":{"iopub.status.busy":"2024-06-26T16:23:10.789496Z","iopub.execute_input":"2024-06-26T16:23:10.789917Z","iopub.status.idle":"2024-06-26T16:23:12.341675Z","shell.execute_reply.started":"2024-06-26T16:23:10.789868Z","shell.execute_reply":"2024-06-26T16:23:12.340259Z"}}},{"cell_type":"code","source":"#train_df = pd.read_csv(\"/kaggle/working/train_BELKA_balanced.csv\")","metadata":{"execution":{"iopub.status.busy":"2024-06-26T19:05:02.476050Z","iopub.execute_input":"2024-06-26T19:05:02.477045Z","iopub.status.idle":"2024-06-26T19:05:02.486483Z","shell.execute_reply.started":"2024-06-26T19:05:02.476987Z","shell.execute_reply":"2024-06-26T19:05:02.485137Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### UPdate to our effort so far\nAfter a number of Iterartions, it appears the most important factor in binding affinity \nis the MorganFingerPrint bit vector. After adding this feature, the data got too large for memory and\nwe had to read the data in chunks and proces them one chunk at a time.\n\nOnce completed all three train data files (15GB), we realized the MFP array had increased each datafile to the point where we had no further room to process and split the test data or even be able to zip the three train data !!!\n\nThis brings to today, when we decided to reduce, by half, the data we read and see of we can get a better result by adding the MFP to the model.","metadata":{}},{"cell_type":"code","source":"# Define the Morgan Finger Print Generator\nmfpgen = rdFingerprintGenerator.GetMorganGenerator(radius=chemical_radious,\n                                                   fpSize=chemical_nbits)","metadata":{"execution":{"iopub.status.busy":"2024-06-26T19:05:02.488381Z","iopub.execute_input":"2024-06-26T19:05:02.488792Z","iopub.status.idle":"2024-06-26T19:05:02.497200Z","shell.execute_reply.started":"2024-06-26T19:05:02.488752Z","shell.execute_reply":"2024-06-26T19:05:02.495813Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def gen_mol(smiles):\n    return Chem.MolFromSmiles(smiles)","metadata":{"execution":{"iopub.status.busy":"2024-06-26T19:05:02.498838Z","iopub.execute_input":"2024-06-26T19:05:02.499728Z","iopub.status.idle":"2024-06-26T19:05:02.514641Z","shell.execute_reply.started":"2024-06-26T19:05:02.499681Z","shell.execute_reply":"2024-06-26T19:05:02.512961Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def gen_prot_mol(pdb_file):\n    return Chem.MolFromPDBFile(pdb_file)","metadata":{"execution":{"iopub.status.busy":"2024-06-26T19:05:02.516257Z","iopub.execute_input":"2024-06-26T19:05:02.517160Z","iopub.status.idle":"2024-06-26T19:05:02.527375Z","shell.execute_reply.started":"2024-06-26T19:05:02.517118Z","shell.execute_reply":"2024-06-26T19:05:02.526091Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"raw","source":"smiles = train_df.loc[0, 'molecule_smiles']\nmol = gen_mol(smiles)\nnp.array(rdMolDescriptors.GetMorganFingerprintAsBitVect(mol, chemical_radious, chemical_nbits))","metadata":{"execution":{"iopub.status.busy":"2024-06-21T01:32:11.099360Z","iopub.execute_input":"2024-06-21T01:32:11.099776Z","iopub.status.idle":"2024-06-21T01:32:11.130468Z","shell.execute_reply.started":"2024-06-21T01:32:11.099745Z","shell.execute_reply":"2024-06-21T01:32:11.129191Z"}}},{"cell_type":"code","source":"# Function to compute molecular descriptors\n# mol : rdkit.Chem.rdchem.Mol - The Molecule Class\n# TPSA: Topological polar surface area\ndef compute_descriptors_from_mol(mol):\n    descriptors = {}\n    if mol is not None:\n        descriptors['NumRotatableBonds'] = np.array(Descriptors.NumRotatableBonds(mol)).astype(np.int8)\n        descriptors['NumHDonors'] = np.array(Descriptors.NumHDonors(mol)).astype(np.int8)\n        descriptors['NumHAcceptors'] = np.array(Descriptors.NumHAcceptors(mol)).astype(np.int8)\n        descriptors['ExactMolWt'] = float(Descriptors.ExactMolWt(mol))\n        descriptors['TPSA'] = float(Descriptors.TPSA(mol))\n        descriptors['NumAromaticRings'] = np.array(Descriptors.NumAromaticRings(mol)).astype(np.int8)\n        descriptors['NumSaturatedRings'] = np.array(Descriptors.NumSaturatedRings(mol)).astype(np.int8)\n        descriptors['NumAliphaticRings'] = np.array(Descriptors.NumAliphaticRings(mol)).astype(np.int8)\n        descriptors['NumHeteroatoms'] = float(Descriptors.NumHeteroatoms(mol))\n        descriptors['NumValenceElectrons'] = np.array(Descriptors.NumValenceElectrons(mol)).astype(np.int8)\n        descriptors['NumRadicalElectrons'] = np.array(Descriptors.NumRadicalElectrons(mol)).astype(np.int8)\n        descriptors['mfp_bv'] = list(mfpgen.GetFingerprint(mol))\n        descriptors['mfp'] = mfpgen.GetFingerprint(mol)\n   \n    return descriptors","metadata":{"execution":{"iopub.status.busy":"2024-06-26T19:05:02.533115Z","iopub.execute_input":"2024-06-26T19:05:02.533594Z","iopub.status.idle":"2024-06-26T19:05:02.547276Z","shell.execute_reply.started":"2024-06-26T19:05:02.533557Z","shell.execute_reply":"2024-06-26T19:05:02.545955Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Function to compute chemical similarity\n# sbv : Sparce Bit Vector\ndef compute_similarity_spv(sbv1, sbv2):\n    similarity = rdkit.DataStructs.cDataStructs.CosineSimilarity(sbv1, sbv2)\n    return similarity","metadata":{"execution":{"iopub.status.busy":"2024-06-26T19:05:02.548691Z","iopub.execute_input":"2024-06-26T19:05:02.549081Z","iopub.status.idle":"2024-06-26T19:05:02.557732Z","shell.execute_reply.started":"2024-06-26T19:05:02.549051Z","shell.execute_reply":"2024-06-26T19:05:02.556288Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## ------  Testing storage of MorganFingerPrint as a list object -------------------->>","metadata":{}},{"cell_type":"raw","source":" pdf = train_df[train_df['protein_name'] == 'sEH'].reset_index(drop=True)","metadata":{"execution":{"iopub.status.busy":"2024-06-26T13:19:56.388020Z","iopub.execute_input":"2024-06-26T13:19:56.388613Z","iopub.status.idle":"2024-06-26T13:19:56.402432Z","shell.execute_reply.started":"2024-06-26T13:19:56.388554Z","shell.execute_reply":"2024-06-26T13:19:56.400730Z"}}},{"cell_type":"raw","source":"mol_1 = gen_mol(pdf.loc[0, 'molecule_smiles'])\nmfp_bv_1 = mfpgen.GetFingerprint(mol_1)\nmfp_1 = list(mfp_bv_1)\n#mfp_1 = np.array(mfp_bv_1, dtype=np.int8).tolist()\n\nmol_2 = gen_mol(pdf.loc[1, 'molecule_smiles'])\nmfp_bv_2 = mfpgen.GetFingerprint(mol_2)\nmfp_2 = np.array(mfp_bv_2, dtype=np.int8).tolist()\n\ncs = compute_similarity_spv(mfp_bv_1, mfp_bv_2)\nmol_1, mfp_1, mfp_bv_1, mol_2, mfp_bv_2, mfp_2, cs","metadata":{}},{"cell_type":"raw","source":"df = pd.read_csv('/kaggle/working/train_sEH.csv')","metadata":{"execution":{"iopub.status.busy":"2024-06-26T15:34:02.841448Z","iopub.execute_input":"2024-06-26T15:34:02.842012Z","iopub.status.idle":"2024-06-26T15:34:02.957999Z","shell.execute_reply.started":"2024-06-26T15:34:02.841972Z","shell.execute_reply":"2024-06-26T15:34:02.956327Z"}}},{"cell_type":"raw","source":"df['molecule_mfp_bv'].values","metadata":{}},{"cell_type":"markdown","source":"## -------------------------->>","metadata":{}},{"cell_type":"code","source":"def process_molecule(df, protein_filename) -> pd.DataFrame:\n    # Compute molecular descriptors and extract molecule objects\n    molecule_descriptors = df['molecule_smiles'].apply(gen_mol).apply(compute_descriptors_from_mol)\n    \n    # Convert dictionaries of descriptors to separate columns\n    df = pd.concat([df, pd.DataFrame(list(molecule_descriptors)).add_prefix('molecule_')], axis=1)\n    \n    #print(\"after 1st concatinations, columns names are : \\n\", df.info())\n    \n    prot_mol = Chem.MolFromPDBFile(protein_filename)\n    protein_descriptors= compute_descriptors_from_mol(prot_mol)\n    \n    # Create a dataframe with a single row for the protine decriptors\n    prot_df= pd.DataFrame([protein_descriptors], index=[0])\n    # Repeat the protein descriptors for each row in the original dataframe\n    prot_df = pd.concat([prot_df]*len(df), ignore_index=True)\n    # Now add this DataFrame to the df so there are as many copies of \n    # protein decriptors as there are moleciules in the data\n    df = pd.concat([df, prot_df.add_prefix('protein_')], axis=1)\n    \n    #print(\"after 2nd concatinations, columns names are : \\n\", df.info())\n\n    # Calculating similarity for each molecule ExplicitBitVector in the subset\n    #df['cosine_sim'] = [compute_similarity_spv(mfp, protein_mfp) for mfp in molecule_mfp]\n    df['cosine_sim'] = df.apply(lambda row: compute_similarity_spv(row['molecule_mfp'], row['protein_mfp']), axis=1)\n    \n    # define datatypes\n\n    df['binds'] = df['binds'].astype(np.int8)\n    df['molecule_NumRotatableBonds']   = df['molecule_NumRotatableBonds'].astype(np.int8)\n    df['molecule_NumHDonors']          = df['molecule_NumHDonors'].astype(np.int8)\n    df['molecule_NumHAcceptors']       = df['molecule_NumHAcceptors'].astype(np.int8)\n    df['molecule_NumAromaticRings']    = df['molecule_NumAromaticRings'].astype(np.int8)\n    df['molecule_NumSaturatedRings']   = df['molecule_NumSaturatedRings'].astype(np.int8)\n    df['molecule_NumAliphaticRings']   = df['molecule_NumAliphaticRings'].astype(np.int8)\n    df['molecule_NumValenceElectrons'] = df['molecule_NumValenceElectrons'].astype(np.int8)\n    df['molecule_NumRadicalElectrons'] = df['molecule_NumRadicalElectrons'].astype(np.int8)\n    #df['molecule_mfp_bv']              = df['molecule_mfp_bv'].astype(np.int8)\n\n    df['molecule_ExactMolWt']     = df['molecule_ExactMolWt'].astype(float)\n    df['molecule_TPSA']           = df['molecule_TPSA'].astype(float)\n    df['molecule_NumHeteroatoms'] = df['molecule_NumHeteroatoms'].astype(float)\n    #df['molecule_mfp']            = df['molecule_mfp'].astype(np.float32)\n    \n    df['protein_NumRotatableBonds']   = df['protein_NumRotatableBonds'].astype(np.int8)\n    df['protein_NumHDonors']          = df['protein_NumHDonors'].astype(np.int8)\n    df['protein_NumHAcceptors']       = df['protein_NumHAcceptors'].astype(np.int8)\n    df['protein_NumAromaticRings']    = df['protein_NumAromaticRings'].astype(np.int8)\n    df['protein_NumSaturatedRings']   = df['protein_NumSaturatedRings'].astype(np.int8)\n    df['protein_NumAliphaticRings']   = df['protein_NumAliphaticRings'].astype(np.int8)\n    df['protein_NumValenceElectrons'] = df['protein_NumValenceElectrons'].astype(np.int8)\n    df['protein_NumRadicalElectrons'] = df['protein_NumRadicalElectrons'].astype(np.int8)\n    #df['protein_mfp_bv']              = df['protein_mfp_bv'].astype(np.int8)\n\n    df['protein_ExactMolWt']     = df['protein_ExactMolWt'].astype(float)\n    df['protein_TPSA']           = df['protein_TPSA'].astype(float)\n    df['protein_NumHeteroatoms'] = df['protein_NumHeteroatoms'].astype(float)\n    #df['protein_mfp']            = df['protein_mfp'].astype(np.float32)\n    \n    df['cosine_sim']              = df['cosine_sim'].astype(float)\n    \n    return df\n","metadata":{"execution":{"iopub.status.busy":"2024-06-26T19:05:02.559652Z","iopub.execute_input":"2024-06-26T19:05:02.560170Z","iopub.status.idle":"2024-06-26T19:05:02.578268Z","shell.execute_reply.started":"2024-06-26T19:05:02.560119Z","shell.execute_reply":"2024-06-26T19:05:02.576966Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Suppressin gwarnings when plotting correlation data\nimport warnings\n\nwarnings.filterwarnings(\"ignore\", category=UserWarning)\nwarnings.filterwarnings(\"ignore\", \"is_categorical_dtype\")\nwarnings.filterwarnings(\"ignore\", \"use_inf_as_na\")","metadata":{"execution":{"iopub.status.busy":"2024-06-26T19:05:02.579988Z","iopub.execute_input":"2024-06-26T19:05:02.580999Z","iopub.status.idle":"2024-06-26T19:05:02.595622Z","shell.execute_reply.started":"2024-06-26T19:05:02.580955Z","shell.execute_reply":"2024-06-26T19:05:02.594391Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"proteins = ['sEH', 'BRD4', 'HSA']\ncolumns_to_drop = [\"buildingblock1_smiles\", \n                       \"buildingblock2_smiles\", \n                       \"buildingblock3_smiles\", \n                       \"protein_name\"]","metadata":{"execution":{"iopub.status.busy":"2024-06-26T19:05:02.596954Z","iopub.execute_input":"2024-06-26T19:05:02.597380Z","iopub.status.idle":"2024-06-26T19:05:02.607138Z","shell.execute_reply.started":"2024-06-26T19:05:02.597343Z","shell.execute_reply":"2024-06-26T19:05:02.606085Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"NaN_Columns = ['binds',\n       'molecule_NumRotatableBonds', 'molecule_NumHDonors',\n       'molecule_NumHAcceptors', 'molecule_ExactMolWt', 'molecule_TPSA',\n       'molecule_NumAromaticRings', 'molecule_NumSaturatedRings',\n       'molecule_NumAliphaticRings',\n       'molecule_NumValenceElectrons', 'molecule_NumRadicalElectrons',\n       'cosine_sim']\nInf_columns = ['molecule_ExactMolWt', 'molecule_TPSA', 'molecule_NumHeteroatoms']","metadata":{"execution":{"iopub.status.busy":"2024-06-26T19:05:02.608975Z","iopub.execute_input":"2024-06-26T19:05:02.609499Z","iopub.status.idle":"2024-06-26T19:05:02.622187Z","shell.execute_reply.started":"2024-06-26T19:05:02.609455Z","shell.execute_reply":"2024-06-26T19:05:02.621148Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Pair plots\npair_plot_columns = [\n       'molecule_NumRotatableBonds', 'molecule_NumHDonors',\n       'molecule_NumHAcceptors', 'molecule_ExactMolWt', 'molecule_TPSA',\n       'molecule_NumAromaticRings', 'molecule_NumSaturatedRings',\n       'molecule_NumAliphaticRings', 'molecule_NumHeteroatoms',\n       'molecule_NumValenceElectrons', 'cosine_sim']","metadata":{"execution":{"iopub.status.busy":"2024-06-26T19:05:02.623557Z","iopub.execute_input":"2024-06-26T19:05:02.623959Z","iopub.status.idle":"2024-06-26T19:05:02.634982Z","shell.execute_reply.started":"2024-06-26T19:05:02.623926Z","shell.execute_reply":"2024-06-26T19:05:02.633781Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def corr_plot(df, prot_name) :\n    correlation_with_target = df[pair_plot_columns + ['binds']].corr()['binds'].drop('binds')\n\n    # Plot the correlations with the target\n    plt.figure(figsize=(10, 6))\n    correlation_with_target.sort_values().plot(kind='barh', color='skyblue')\n    plt.title(f'Correlation of Features with Binding Afinity ({prot_name})')\n    plt.xlabel('Correlation coefficient')\n    plt.ylabel('Features')\n    plt.show()","metadata":{"execution":{"iopub.status.busy":"2024-06-26T19:05:02.636483Z","iopub.execute_input":"2024-06-26T19:05:02.636958Z","iopub.status.idle":"2024-06-26T19:05:02.648241Z","shell.execute_reply.started":"2024-06-26T19:05:02.636924Z","shell.execute_reply":"2024-06-26T19:05:02.646983Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def read_data(protein_name, chunk_size, offset):\n    train_df = con.query(f\"\"\"(SELECT *\n                            FROM parquet_scan('{train_path}')\n                            WHERE binds = 0 AND protein_name = ('{protein_name}')\n                            ORDER BY id\n                            LIMIT {chunk_size} OFFSET {offset})\n                            UNION ALL\n                            (SELECT *\n                            FROM parquet_scan('{train_path}')\n                            WHERE binds = 1 AND protein_name = ('{protein_name}')\n                            ORDER BY id\n                            LIMIT {chunk_size} OFFSET {offset})\"\"\").df()\n    return train_df","metadata":{"execution":{"iopub.status.busy":"2024-06-26T19:09:27.126276Z","iopub.execute_input":"2024-06-26T19:09:27.126821Z","iopub.status.idle":"2024-06-26T19:09:27.135171Z","shell.execute_reply.started":"2024-06-26T19:09:27.126785Z","shell.execute_reply":"2024-06-26T19:09:27.133222Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def process_chunk(pdf, p):\n    \n    print(f\" Processing chunk of data for protein {p}\")\n    pdf.dropna(axis=0, how='any', inplace=True)\n\n    print(\" N/A values : \\n\", pdf.isna().any())\n\n    # Shuffle the data\n    pdf = pdf.sample(frac=1, random_state=42).reset_index(drop=True)\n    \n    print(\"adding fearures ...\")\n    if p == 'HSA':\n        pdb_filename = f\"/kaggle/input/protein-pdb-data/ALB-{p}.pdb\"\n    else:\n        pdb_filename = f\"/kaggle/input/protein-pdb-data/{p}.pdb\"\n    pdf = process_molecule(pdf, pdb_filename)\n    \n    # Drop the BBloch columns\n    print(\"Droping unwanted columns\")\n    pdf.drop(columns=columns_to_drop, inplace=True)\n\n    print(f\"train data file with featires for protein {p}:\\n\\n {pdf.info()}\")\n    \n    nan_values_exist = pdf[NaN_Columns].isna().any()    # Any NAN?\n    inf_values_exist = np.isinf(pdf[Inf_columns]).any() # Amy Inf?\n\n    print(\"NaN values exist:\\n\", nan_values_exist)\n    print(\"Infinite values exist:\\n\", inf_values_exist)\n    print(\"lenngth of data : \", len(pdf))\n    \n    \n    # Plot some sample data\n    # As the data is too large for pair plot, take 20,000 samples\n    if(len(pdf) > 10_000):\n        binds_1_sample = pdf[pdf['binds'] == 1].sample(n=10000, random_state=42)\n        binds_0_sample = pdf[pdf['binds'] == 0].sample(n=10000, random_state=42)\n        sampled_df = pd.concat([binds_1_sample, binds_0_sample])\n        sampled_df.reset_index(drop=True, inplace=True)\n        corr_plot(sampled_df, p)\n\n        del sampled_df\n    gc.collect()\n    \n    return pdf","metadata":{"execution":{"iopub.status.busy":"2024-06-26T19:09:29.842490Z","iopub.execute_input":"2024-06-26T19:09:29.842964Z","iopub.status.idle":"2024-06-26T19:09:29.854328Z","shell.execute_reply.started":"2024-06-26T19:09:29.842930Z","shell.execute_reply":"2024-06-26T19:09:29.853174Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"\nfor p in proteins:\n    \n    offset = 0\n    total_rows_processed = 0\n    first_chunk = True\n\n    while total_rows_processed < data_cap:\n        pdf = {}\n        chunk_df = {}\n        \n        remaining_rows = data_cap - total_rows_processed\n        current_chunk_size = min(chunk_size, remaining_rows)\n        \n        # Read a chunk of data\n        chunk_df = read_data(p, current_chunk_size, offset)\n        \n        if chunk_df.empty:\n            break\n        \n        # Process the chunk\n        pdf = process_chunk(chunk_df, p)\n        \n        fname = f\"/kaggle/working/train_{p}.csv\"\n    \n        print(f\"\\n\\n Writing {total_rows_processed + current_chunk_size} training data for protein {p}\")\n    \n        # Write the processed chunk to CSV\n        if first_chunk:\n            pdf.to_csv(fname, index=False, mode='w', header=True)\n            first_chunk = False\n        else:\n            pdf.to_csv(fname, index=False, mode='a', header=False)\n        \n        # Update the offset and total rows processed\n        offset += current_chunk_size\n        total_rows_processed += current_chunk_size\n","metadata":{"execution":{"iopub.status.busy":"2024-06-26T19:09:34.903220Z","iopub.execute_input":"2024-06-26T19:09:34.903643Z","iopub.status.idle":"2024-06-26T21:20:35.637835Z","shell.execute_reply.started":"2024-06-26T19:09:34.903611Z","shell.execute_reply":"2024-06-26T21:20:35.635945Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"con.close()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import zipfile\nimport os\n\npath = \"/kaggle/working\"\nhsa = os.path.join(path, 'train_HSA.csv')\nbrd = os.path.join(path, 'train_BRD4.csv')\nseh = os.path.join(path, 'train_sEH.csv')\nfiles_to_compress = [hsa, brd, seh]\n\n# Name of the output zip file\noutput_zip_file = 'compressed_BELKA_tain_data.zip'\n\n# Create a ZipFile object in write mode\nwith zipfile.ZipFile(output_zip_file, 'w') as zipf:\n    for file in files_to_compress:\n        # Add file to the zip file\n        zipf.write(file, os.path.basename(file))\n\nprint(f'Files compressed successfully into {output_zip_file}')","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#hsa_train_df = pd.read_csv(\"/kaggle/working/train_HSA.csv\")","metadata":{"execution":{"iopub.status.busy":"2024-06-26T19:05:02.826415Z","iopub.status.idle":"2024-06-26T19:05:02.827066Z","shell.execute_reply.started":"2024-06-26T19:05:02.826776Z","shell.execute_reply":"2024-06-26T19:05:02.826798Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"raw","source":"# Now let's look at some graphs on hsa\nimport warnings\n# Suppress the warning\nwarnings.filterwarnings(\"ignore\", category=UserWarning)\nsns.scatterplot(data=hsa_train_df, x=hsa_train_df.index, y='cosine_sim', hue='binds')\nplt.legend(loc='center left', bbox_to_anchor=(1.02, 0.5))  # Place legend outside plot area on the right\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2024-06-13T02:22:53.556984Z","iopub.execute_input":"2024-06-13T02:22:53.557525Z","iopub.status.idle":"2024-06-13T02:23:10.706675Z","shell.execute_reply.started":"2024-06-13T02:22:53.557469Z","shell.execute_reply":"2024-06-13T02:23:10.705246Z"}}},{"cell_type":"markdown","source":"##### Get some pair plots to see if there are any correlations apparent","metadata":{}},{"cell_type":"raw","source":"from IPython.display import Markdown\nimport sys\n\ndef pair_Plot(df):\n\n    hue = 'binds'\n    vars_per_line = 5\n    #all_vars = list(sampled_df.columns.symmetric_difference([hue]))\n    all_vars = pair_plot_columns\n\n    for var in all_vars:\n        rest_vars = list(all_vars)\n        rest_vars.remove(var)\n        display(Markdown(f\"#### {var}\"))\n        while rest_vars:\n            line_vars = rest_vars[:vars_per_line]\n            del rest_vars[:vars_per_line]\n            line_var_names = \", \".join(line_vars)\n            display(Markdown(f\"##### {var} vs {line_var_names}\"))\n            sns.pairplot(df, x_vars=line_vars, y_vars=[var], hue=hue, palette='bright', )\n            plt.show()\n            plt.close()\n","metadata":{"execution":{"iopub.status.busy":"2024-04-26T14:14:23.254309Z","iopub.execute_input":"2024-04-26T14:14:23.254761Z","iopub.status.idle":"2024-04-26T14:17:06.716781Z","shell.execute_reply.started":"2024-04-26T14:14:23.254729Z","shell.execute_reply":"2024-04-26T14:17:06.715449Z"}}},{"cell_type":"raw","source":"# Look at Principal Component Analysis of the data\nfrom sklearn.decomposition import PCA\nfrom sklearn.preprocessing import StandardScaler\nfrom mpl_toolkits.mplot3d import Axes3D\n\nx = hsa_train_df.loc[:, pair_plot_columns].values\nx = StandardScaler().fit_transform(x)\n\n# Applying PCA\npca = PCA(n_components=3)\nprincipalComponents = pca.fit_transform(x)\n# Create a DataFrame with the PCA results\npca_df = pd.DataFrame(data=principalComponents, columns=['PC1', 'PC2', 'PC3'])\n\n# 3D Plot\nfig = plt.figure(figsize=(10, 8))\nax = fig.add_subplot(111, projection='3d')\n#ax = Axes3D(fig)\nax.scatter(pca_df['PC1'], pca_df['PC2'], pca_df['PC3'], c='blue', marker='o')\n\n# Add labels\nax.set_xlabel('Principal Component 1')\nax.set_ylabel('Principal Component 2')\nax.set_zlabel('Principal Component 3')\n\n# Title\nax.set_title('3D PCA Plot')\n\n# Show plot\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2024-06-14T00:22:44.988351Z","iopub.execute_input":"2024-06-14T00:22:44.989241Z","iopub.status.idle":"2024-06-14T00:23:03.099567Z","shell.execute_reply.started":"2024-06-14T00:22:44.989199Z","shell.execute_reply":"2024-06-14T00:23:03.098286Z"}}},{"cell_type":"markdown","source":"","metadata":{}}]}