{"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":"nvidiaTeslaT4","dataSources":[{"sourceId":67356,"databundleVersionId":8006601,"sourceType":"competition"},{"sourceId":8859975,"sourceType":"datasetVersion","datasetId":4881163}],"dockerImageVersionId":30733,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":true}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"!pip install duckdb rdkit mordred","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import pandas as pd\nimport joblib\nimport duckdb\n\nfrom rdkit import Chem\nfrom rdkit.Chem import AllChem, Descriptors, Crippen, rdMolDescriptors\nfrom sklearn.preprocessing import MinMaxScaler, MaxAbsScaler\nimport numpy as np\nfrom mordred import Calculator, descriptors\n\nimport multiprocessing\nimport concurrent.futures\n","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_path = '/kaggle/input/leash-BELKA/train.parquet'\n# test_path = '/kaggle/input/leash-BELKA/test.parquet'\nrange = 10000 #20000000\ncon = duckdb.connect()\n\ndf = con.query(f\"\"\"(SELECT *\n                        FROM parquet_scan('{train_path}')\n                        WHERE binds = 0\n                        ORDER BY random()\n                        LIMIT {range})\n                        UNION ALL\n                        (SELECT *\n                        FROM parquet_scan('{train_path}')\n                        WHERE binds = 1)\"\"\").df()\n\ncon.close()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# df = pd.read_csv('/kaggle/input/leash-bio-belka-chunk/somepart_normalized/somepart_normalized.csv')","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from rdkit import Chem\nfrom rdkit.Chem import AllChem, Descriptors, Crippen, rdMolDescriptors\nfrom sklearn.preprocessing import MinMaxScaler, MaxAbsScaler\nfrom mordred import Calculator, descriptors\n","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Initialize Mordred calculator with all descriptors\ncalc = Calculator(descriptors, ignore_3D=True)\n\ndef mol_descriptors(smiles_list, pre=''):\n    moldescriptors = []\n    for smile in smiles_list:\n        mol = Chem.MolFromSmiles(smile)\n        if mol:\n            rdkit_desc = {\n                pre+'MolecularWeight': Descriptors.MolWt(mol),\n                pre+'LogP': Descriptors.MolLogP(mol),\n                pre+'TPSA': Descriptors.TPSA(mol),\n                pre+'NumRotatableBonds': Descriptors.NumRotatableBonds(mol),\n                pre+'NumHDonors': Descriptors.NumHDonors(mol),\n                pre+'NumHAcceptors': Descriptors.NumHAcceptors(mol),\n                pre+'NumRings': Descriptors.RingCount(mol),\n                pre+'NumAromaticRings': Descriptors.NumAromaticRings(mol),\n                pre+'ExactMass': Descriptors.ExactMolWt(mol),\n                pre+'HeavyAtomCount': Descriptors.HeavyAtomCount(mol),\n                pre+'NumValenceElectrons': Descriptors.NumValenceElectrons(mol),\n                pre+'FractionCSP3': Descriptors.FractionCSP3(mol),\n                pre+'MolMR': Descriptors.MolMR(mol),\n                pre+'FormalCharge': Chem.GetFormalCharge(mol),\n                pre+'NumAliphaticRings': Descriptors.NumAliphaticRings(mol),\n                pre+'NumSaturatedRings': Descriptors.NumSaturatedRings(mol),\n                pre+'NumHeteroatoms': Descriptors.NumHeteroatoms(mol),\n                pre+'NumSaturatedCarbocycles': Descriptors.NumSaturatedCarbocycles(mol),\n                pre+'NumAliphaticHeterocycles': Descriptors.NumAliphaticHeterocycles(mol),\n                pre+'NumAromaticHeterocycles': Descriptors.NumAromaticHeterocycles(mol),\n                pre+'AtomCount': mol.GetNumAtoms(),\n                pre+'NumSingleBonds': len([bond for bond in mol.GetBonds() if bond.GetBondType() == Chem.rdchem.BondType.SINGLE]),\n                pre+'NumDoubleBonds': len([bond for bond in mol.GetBonds() if bond.GetBondType() == Chem.rdchem.BondType.DOUBLE]),\n                pre+'NumTripleBonds': len([bond for bond in mol.GetBonds() if bond.GetBondType() == Chem.rdchem.BondType.TRIPLE]),\n                pre+'NumAromaticBonds': len([bond for bond in mol.GetBonds() if bond.GetIsAromatic()]),\n                pre+'MolecularConnectivityIndex': Descriptors.MolLogP(mol),  # Example\n                pre+'Kier_Hall_Alpha': Descriptors.Kappa3(mol),  # Example\n                pre+'HOMO': Descriptors.MaxAbsEStateIndex(mol),  # Example (approximation)\n                pre+'LUMO': Descriptors.MinAbsEStateIndex(mol),  # Example (approximation)\n            }\n            \n            # Calculate Mordred descriptors\n            mordred_desc = calc(mol)\n            mordred_desc = {f\"{pre}{str(key)}\": value for key, value in mordred_desc.items()}\n            \n            # Combine RDKit and Mordred descriptors\n            rdkit_desc.update(mordred_desc)\n            moldescriptors.append(rdkit_desc)\n            \n        else:\n            # If molecule is invalid, append None for all descriptors\n            rdkit_desc = {key: None for key in [\n                pre+'MolecularWeight',\n                pre+'LogP',\n                pre+'TPSA',\n                pre+'NumRotatableBonds',\n                pre+'NumHDonors',\n                pre+'NumHAcceptors',\n                pre+'NumRings',\n                pre+'NumAromaticRings',\n                pre+'ExactMass',\n                pre+'HeavyAtomCount',\n                pre+'NumValenceElectrons',\n                pre+'FractionCSP3',\n                pre+'MolMR',\n                pre+'FormalCharge',\n                pre+'NumAliphaticRings',\n                pre+'NumSaturatedRings',\n                pre+'NumHeteroatoms',\n                pre+'NumSaturatedCarbocycles',\n                pre+'NumAliphaticHeterocycles',\n                pre+'NumAromaticHeterocycles',\n                pre+'AtomCount',\n                pre+'NumSingleBonds',\n                pre+'NumDoubleBonds',\n                pre+'NumTripleBonds',\n                pre+'NumAromaticBonds',\n                pre+'MolecularConnectivityIndex',\n                pre+'Kier_Hall_Alpha',\n                pre+'HOMO',\n                pre+'LUMO',\n            ]}\n            \n            # Add None for Mordred descriptors\n            mordred_desc = {f\"{pre}{str(desc)}\": None for desc in calc.descriptors}\n            rdkit_desc.update(mordred_desc)\n            \n            moldescriptors.append(rdkit_desc)\n            \n    return moldescriptors\n\ndef parallelize_dataframe(df, molecule_smiles_cat, pre='', n_cores=80):\n    df_split = np.array_split(df, n_cores)\n    smiles_lists = [chunk[molecule_smiles_cat].tolist() for chunk in df_split]\n    \n    with concurrent.futures.ProcessPoolExecutor(max_workers=n_cores) as executor:\n        results = executor.map(mol_descriptors, smiles_lists, [pre]*n_cores)\n    \n    results = [item for sublist in results for item in sublist]\n    \n    descriptors_df = pd.DataFrame(results)\n    return pd.concat([df.reset_index(drop=True), descriptors_df.reset_index(drop=True)], axis=1)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def normalize_df(df, columns_to_scale, scalers):\n    for column in columns_to_scale:\n        if column not in scalers:\n            scaler = MaxAbsScaler()\n            scalers[column] = scaler  # Save the scaler in case you need it later\n            df[column] = scaler.fit_transform(df[[column]])\n        else:\n            df[column] = scalers[column].transform(df[[column]])\n    return df","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# multiprocessing.cpu_count()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"p1 = parallelize_dataframe(df, 'buildingblock1_smiles', 'b1')","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"p2 = parallelize_dataframe(b1, 'buildingblock2_smiles', 'b2')","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"b3 = parallelize_dataframe(b2, 'buildingblock3_smiles', 'b3')","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"scalers = {}\nnormalized = normalize_df(b3, [\n  'b1MolecularWeight',\n  'b1LogP',\n  'b1TPSA',\n  'b1NumRotatableBonds',\n  'b1NumHDonors',\n  'b1NumHAcceptors',\n  'b1NumRings',\n  'b1NumAromaticRings',\n  'b1ExactMass',\n  'b1HeavyAtomCount',\n  'b1NumValenceElectrons',\n  'b1FractionCSP3',\n  'b1MolMR',\n  'b1FormalCharge',\n  'b1NumAliphaticRings',\n  'b1NumSaturatedRings',\n  'b1NumHeteroatoms',\n  'b1NumHeterocycles',\n  'b1NumSaturatedCarbocycles',\n  'b1NumAliphaticHeterocycles',\n  'b1NumAromaticHeterocycles',\n  'b1BondCount',\n  'b1AtomCount',\n  'b1NumSingleBonds',\n  'b1NumDoubleBonds',\n  'b1NumTripleBonds',\n  'b1NumAromaticBonds',\n\n  'b2MolecularWeight',\n  'b2LogP',\n  'b2TPSA',\n  'b2NumRotatableBonds',\n  'b2NumHDonors',\n  'b2NumHAcceptors',\n  'b2NumRings',\n  'b2NumAromaticRings',\n  'b2ExactMass',\n  'b2HeavyAtomCount',\n  'b2NumValenceElectrons',\n  'b2FractionCSP3',\n  'b2MolMR',\n  'b2FormalCharge',\n  'b2NumAliphaticRings',\n  'b2NumSaturatedRings',\n  'b2NumHeteroatoms',\n  'b2NumHeterocycles',\n  'b2NumSaturatedCarbocycles',\n  'b2NumAliphaticHeterocycles',\n  'b2NumAromaticHeterocycles',\n  'b2BondCount',\n  'b2AtomCount',\n  'b2NumSingleBonds',\n  'b2NumDoubleBonds',\n  'b2NumTripleBonds',\n  'b2NumAromaticBonds',\n\n  'b3MolecularWeight',\n  'b3LogP',\n  'b3TPSA',\n  'b3NumRotatableBonds',\n  'b3NumHDonors',\n  'b3NumHAcceptors',\n  'b3NumRings',\n  'b3NumAromaticRings',\n  'b3ExactMass',\n  'b3HeavyAtomCount',\n  'b3NumValenceElectrons',\n  'b3FractionCSP3',\n  'b3MolMR',\n  'b3FormalCharge',\n  'b3NumAliphaticRings',\n  'b3NumSaturatedRings',\n  'b3NumHeteroatoms',\n  'b3NumHeterocycles',\n  'b3NumSaturatedCarbocycles',\n  'b3NumAliphaticHeterocycles',\n  'b3NumAromaticHeterocycles',\n  'b3BondCount',\n  'b3AtomCount',\n  'b3NumSingleBonds',\n  'b3NumDoubleBonds',\n  'b3NumTripleBonds',\n  'b3NumAromaticBonds',\n], scalers)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import pandas as pd\nimport tensorflow as tf\nfrom tensorflow.keras.models import Sequential, Model\nfrom tensorflow.keras.layers import Dense, Dropout, BatchNormalization, Input, Concatenate, Embedding, GlobalAveragePooling1D, Conv1D, GlobalMaxPooling1D, Attention\n# from tensorflow.keras.regularizers import l2, l1\nfrom tensorflow.keras.callbacks import LearningRateScheduler, EarlyStopping\nimport matplotlib.pyplot as plt\nfrom sklearn.metrics import confusion_matrix\nimport numpy as np\nimport sklearn.metrics\n# from sklearn.preprocessing import RobustScaler, MinMaxScaler, MaxAbsScaler\nimport json\nfrom tensorflow.keras.optimizers import Adam\nfrom tensorflow.keras.callbacks import EarlyStopping, ReduceLROnPlateau\nfrom sklearn.utils.class_weight import compute_class_weight\n","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### SMILES Vectorization","metadata":{}},{"cell_type":"code","source":"smile_vector = tf.keras.layers.TextVectorization(max_tokens=None, split='character', standardize=None, output_sequence_length=70)\nsmile_vector.set_vocabulary([\"\", \"[UNK]\", \"c\", \"C\", \"1\", \")\", \"(\", \"O\", \"2\", \"N\", \"=\", \"n\", \"-\", \"l\", \"]\", \"[\", \"@\", \"H\", \"F\", \".\", \"3\", \"s\", \"B\", \"r\", \"S\", \"#\", \"+\", \"o\", \"I\", \"4\", \"/\", \"5\", \"i\"])","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# df['binds'].values.flatten()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### encode the proteins","metadata":{}},{"cell_type":"code","source":"insert_loc = df.columns.get_loc('protein_name')\ndata = pd.concat([df.iloc[:, :insert_loc], pd.get_dummies(df.loc[:, ['protein_name']], dtype=int), df.iloc[:, insert_loc+1:]], axis=1)\nlen(data)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# data = data[:10]","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# target_protein_name = 'HSA'\n# data = normalized[normalized['protein_name'] == target_protein_name].sample(frac=1)\n# len(data)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"data[['protein_name_BRD4', 'protein_name_sEH', 'protein_name_HSA']]","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"buildingblock1_smiles_vec = smile_vector(data['buildingblock1_smiles'].values)\nbuildingblock2_smiles_vec = smile_vector(data['buildingblock2_smiles'].values)\nbuildingblock3_smiles_vec = smile_vector(data['buildingblock3_smiles'].values)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"p1 = int(6*len(data)/10)\np2 = int(9*len(data)/10)\nprint(p1, p2, len(data))","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Balance Data w/ Weights\nclass_weights = compute_class_weight('balanced', classes=np.unique(data['binds'].values.flatten()), y=data['binds'].values.flatten())\nclass_weight_dict = dict(enumerate(class_weights))","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"vector_length = buildingblock3_smiles_vec.shape[1]\nvector_length","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"data_cols = [\n    'protein_name_BRD4', 'protein_name_sEH', 'protein_name_HSA',\n 'b1LogP',\n 'b1MolMR',\n 'b1TPSA',\n 'b1FractionCSP3',\n 'b1NumHeteroatoms',\n 'b1MolecularWeight',\n 'b1ExactMass',\n 'b1NumRotatableBonds',\n 'b1NumValenceElectrons',\n 'b1BondCount',\n 'b1NumHAcceptors',\n\n 'b2LogP',\n 'b2MolMR',\n 'b2TPSA',\n 'b2FractionCSP3',\n 'b2NumValenceElectrons',\n 'b2MolecularWeight',\n 'b2ExactMass',\n 'b2NumRotatableBonds',\n 'b2BondCount',\n 'b2NumHeteroatoms',\n 'b2NumHAcceptors',\n \n 'b3LogP',\n 'b3MolMR',\n 'b3TPSA',\n 'b3FractionCSP3',\n 'b3NumHAcceptors',\n 'b3MolecularWeight',\n 'b3ExactMass',\n 'b3BondCount',\n 'b3NumValenceElectrons',\n 'b3NumRotatableBonds',\n 'b3NumHeteroatoms',\n 'b3NumAromaticBonds',\n]","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"data[data_cols][:p1].values","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Define inputs\ninput1 = Input(shape=(vector_length,))\ninput2 = Input(shape=(vector_length,))\ninput3 = Input(shape=(vector_length,))\ninput4 = Input(shape=(len(data_cols),))\n\n# Embedding layer\nembedding_dim = 128  # Increase embedding dimension\nembedded1 = Embedding(input_dim=len(smile_vector.get_vocabulary()), output_dim=embedding_dim)(input1)\nembedded2 = Embedding(input_dim=len(smile_vector.get_vocabulary()), output_dim=embedding_dim)(input2)\nembedded3 = Embedding(input_dim=len(smile_vector.get_vocabulary()), output_dim=embedding_dim)(input3)\n\n# Convolutional and pooling layers\nconv_filters = 128\nconv1 = Conv1D(filters=conv_filters, kernel_size=5, activation='relu', padding='same')(embedded1)\nconv2 = Conv1D(filters=conv_filters, kernel_size=5, activation='relu', padding='same')(embedded2)\nconv3 = Conv1D(filters=conv_filters, kernel_size=5, activation='relu', padding='same')(embedded3)\n\n# Attention layers\nattention1 = Attention()([conv1, conv1])\nattention2 = Attention()([conv2, conv2])\nattention3 = Attention()([conv3, conv3])\n\n# Pooling layers\npooled1 = GlobalMaxPooling1D()(attention1)\npooled2 = GlobalMaxPooling1D()(attention2)\npooled3 = GlobalMaxPooling1D()(attention3)\n\n# Concatenate inputs\nconcat = Concatenate()([pooled1, pooled2, pooled3])\n# concat = Concatenate()([embedded1, embedded2, embedded3])\n\n# Dense layers\ndense1 = Dense(256, activation='relu')(concat)\ndropout1 = Dropout(0.3)(dense1)\n\n# Concatenate with numerical input\nconcat_with_numerical = Concatenate()([dropout1, input4])\n\ndense2 = Dense(128, activation='relu')(concat_with_numerical)\ndropout2 = Dropout(0.1)(dense2)\nbatch_norm = BatchNormalization()(dropout2)\ndense3 = Dense(64, activation='relu')(batch_norm)\noutput = Dense(1, activation='sigmoid')(dense3)  # For binary classification\n\n# Create the model\nmodel = Model(inputs=[input1, input2, input3, input4], outputs=output)\n\n# Compile the model\nmodel.compile(optimizer='adam', loss='binary_crossentropy', metrics=['accuracy'])\n\n# Callbacks for early stopping and learning rate reduction\nearly_stopping = EarlyStopping(monitor='val_loss', patience=5, restore_best_weights=True)\nreduce_lr = ReduceLROnPlateau(monitor='val_loss', factor=0.5, patience=2)\n\n# Train the model with validation data\nhistory = model.fit(\n  [buildingblock1_smiles_vec[:p1], buildingblock2_smiles_vec[:p1], buildingblock3_smiles_vec[:p1], data[data_cols][:p1].values],\n  data['binds'][:p1].values.flatten(),\n  validation_data=(\n      [buildingblock1_smiles_vec[p1:p2], buildingblock2_smiles_vec[p1:p2], buildingblock3_smiles_vec[p1:p2], data[data_cols][p1:p2].values],\n      data['binds'][p1:p2].values.flatten()),\n  batch_size=30,  # Decrease batch size for more updates\n  epochs=100,\n  class_weight=class_weight_dict,\n  callbacks=[early_stopping, reduce_lr])","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Evaluate the model\nloss, accuracy = model.evaluate([buildingblock1_smiles_vec[p2:], buildingblock2_smiles_vec[p2:], buildingblock3_smiles_vec[p2:], data[data_cols][p2:]], data['binds'][p2:].values.flatten())\nprint(f'Test Accuracy: {accuracy:.4f}')\nprint(f'Test Loss: {loss:.4f}')\n","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# model.save('/kaggle/working/'+target_protein_name+'.keras')","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"model.save('/kaggle/working/all.keras')","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"testds = pd.read_csv('/kaggle/input/leash-bio-belka-chunk/test_normalized/test_normalized.csv')\n","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"testbuildingblock1_smiles_vec = smile_vector(testds['buildingblock1_smiles'].values)\ntestbuildingblock2_smiles_vec = smile_vector(testds['buildingblock2_smiles'].values)\ntestbuildingblock3_smiles_vec = smile_vector(testds['buildingblock3_smiles'].values)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"predictions = model.predict([testbuildingblock1_smiles_vec, testbuildingblock2_smiles_vec, testbuildingblock3_smiles_vec, testds[data_cols]])\ntestds['binds'] = predictions.flatten()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"testds[['id', 'binds']].to_csv('submition.csv', index=False)","metadata":{"trusted":true},"execution_count":null,"outputs":[]}]}