{"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":"gpu","dataSources":[{"sourceId":41875,"databundleVersionId":5521661,"sourceType":"competition"},{"sourceId":9075045,"sourceType":"datasetVersion","datasetId":5474447}],"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":true}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"!pip install biopython\n!pip install obonet","metadata":{"execution":{"iopub.status.busy":"2024-07-31T15:22:19.654823Z","iopub.execute_input":"2024-07-31T15:22:19.655097Z","iopub.status.idle":"2024-07-31T15:22:46.518925Z","shell.execute_reply.started":"2024-07-31T15:22:19.655071Z","shell.execute_reply":"2024-07-31T15:22:46.517995Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# This Python 3 environment comes with many helpful analytics libraries installed\n# It is defined by the kaggle/python Docker image: https://github.com/kaggle/docker-python\n# For example, here's several helpful packages to load\n\nimport numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\n\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\nfrom Bio import SeqIO\nimport networkx\nimport obonet","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2024-07-31T15:22:46.520935Z","iopub.execute_input":"2024-07-31T15:22:46.521238Z","iopub.status.idle":"2024-07-31T15:22:47.272301Z","shell.execute_reply.started":"2024-07-31T15:22:46.52121Z","shell.execute_reply":"2024-07-31T15:22:47.271558Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def import_sequences(fasta_path):\n    records = list(SeqIO.parse(fasta_path, \"fasta\"))\n    data = []\n    for record in records:\n        data.append({\"id\": record.id, \"name\": record.name, \"description\": record.description, \"sequence\": str(record.seq)})\n    return pd.DataFrame(data)","metadata":{"execution":{"iopub.status.busy":"2024-07-31T15:22:47.27335Z","iopub.execute_input":"2024-07-31T15:22:47.273758Z","iopub.status.idle":"2024-07-31T15:22:47.27919Z","shell.execute_reply.started":"2024-07-31T15:22:47.273733Z","shell.execute_reply":"2024-07-31T15:22:47.278284Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sequences = import_sequences(\"/kaggle/input/cafa-5-protein-function-prediction/Train/train_sequences.fasta\")\nsequences","metadata":{"execution":{"iopub.status.busy":"2024-07-31T15:22:47.281077Z","iopub.execute_input":"2024-07-31T15:22:47.28133Z","iopub.status.idle":"2024-07-31T15:22:51.605845Z","shell.execute_reply.started":"2024-07-31T15:22:47.281308Z","shell.execute_reply":"2024-07-31T15:22:51.604996Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"targets = pd.read_csv(\"/kaggle/input/cafa-5-protein-function-prediction/Train/train_terms.tsv\", sep=\"\\t\")\ntargets","metadata":{"execution":{"iopub.status.busy":"2024-07-31T15:22:51.607175Z","iopub.execute_input":"2024-07-31T15:22:51.60753Z","iopub.status.idle":"2024-07-31T15:22:55.672193Z","shell.execute_reply.started":"2024-07-31T15:22:51.607498Z","shell.execute_reply":"2024-07-31T15:22:55.671242Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"go_graph = obonet.read_obo(\"/kaggle/input/cafa-5-protein-function-prediction/Train/go-basic.obo\")\nnames = networkx.get_node_attributes(go_graph, \"name\")\ndf_terms = pd.json_normalize(names).T\ndf_terms[\"id\"] = df_terms.index\ndf_terms = df_terms.rename(columns={0:\"term\"})\ndf_terms = df_terms.reset_index(drop=True)[[\"id\", \"term\"]]\ndf_terms","metadata":{"execution":{"iopub.status.busy":"2024-07-31T15:22:55.673445Z","iopub.execute_input":"2024-07-31T15:22:55.673737Z","iopub.status.idle":"2024-07-31T15:23:14.23475Z","shell.execute_reply.started":"2024-07-31T15:22:55.673712Z","shell.execute_reply":"2024-07-31T15:23:14.233873Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"id_dna_binding = df_terms[df_terms.term == \"DNA binding\"].id.values[0]\nentry_dna_binding = targets[targets.term == id_dna_binding]\ndna_binding_sequences = sequences.merge(entry_dna_binding, left_on=\"id\", right_on=\"EntryID\")\ndna_binding_sequences = dna_binding_sequences[[\"sequence\"]].drop_duplicates().reset_index(drop=True)\ndna_binding_sequences","metadata":{"execution":{"iopub.status.busy":"2024-07-31T15:23:14.235931Z","iopub.execute_input":"2024-07-31T15:23:14.23621Z","iopub.status.idle":"2024-07-31T15:23:15.228405Z","shell.execute_reply.started":"2024-07-31T15:23:14.236185Z","shell.execute_reply":"2024-07-31T15:23:15.227536Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Preprocessing","metadata":{}},{"cell_type":"code","source":"alphabet = \"ACDEFGHIKLMNPQRSTVWY\"","metadata":{"execution":{"iopub.status.busy":"2024-07-31T15:23:15.229512Z","iopub.execute_input":"2024-07-31T15:23:15.229773Z","iopub.status.idle":"2024-07-31T15:23:15.234002Z","shell.execute_reply.started":"2024-07-31T15:23:15.229752Z","shell.execute_reply":"2024-07-31T15:23:15.233146Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def eval_sequences(sequence):\n    if set(alphabet) == set(sequence):\n        return True\n    return False\ndna_binding_sequences[\"is_canon\"] = dna_binding_sequences.map(eval_sequences)\ndna_binding_sequences = dna_binding_sequences[dna_binding_sequences.is_canon].drop(columns=[\"is_canon\"]).reset_index(drop=True)\ndna_binding_sequences","metadata":{"execution":{"iopub.status.busy":"2024-07-31T15:23:15.235057Z","iopub.execute_input":"2024-07-31T15:23:15.235331Z","iopub.status.idle":"2024-07-31T15:23:15.310815Z","shell.execute_reply.started":"2024-07-31T15:23:15.235308Z","shell.execute_reply":"2024-07-31T15:23:15.309723Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import math\ndna_binding_sequences[\"length\"] = dna_binding_sequences.sequence.str.len()\nq3 = dna_binding_sequences.length.quantile(0.75)\nq1 = dna_binding_sequences.length.quantile(0.25)\nmax_length = q3 + 1.5*(q3 - q1)\na = math.log2(max_length) - 1\nmax_length = math.pow(2, math.floor(a))\ndna_binding_sequences = dna_binding_sequences[dna_binding_sequences.length <= max_length]\ndna_binding_sequences = dna_binding_sequences.drop(columns=[\"length\"])\nprint(max_length)\ndna_binding_sequences = dna_binding_sequences.dropna()\ndna_binding_sequences","metadata":{"execution":{"iopub.status.busy":"2024-07-31T15:23:15.313856Z","iopub.execute_input":"2024-07-31T15:23:15.314163Z","iopub.status.idle":"2024-07-31T15:23:15.337088Z","shell.execute_reply.started":"2024-07-31T15:23:15.314141Z","shell.execute_reply":"2024-07-31T15:23:15.336279Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"encoding = pd.read_csv(\"/kaggle/input/clustering-encoding-residues/clustering_encoders.csv\", index_col=0).T\nencoding_length = 8\nencoding","metadata":{"execution":{"iopub.status.busy":"2024-07-31T15:23:15.338157Z","iopub.execute_input":"2024-07-31T15:23:15.3384Z","iopub.status.idle":"2024-07-31T15:23:15.368438Z","shell.execute_reply.started":"2024-07-31T15:23:15.338378Z","shell.execute_reply":"2024-07-31T15:23:15.367556Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from sklearn.preprocessing import MinMaxScaler\nscaler = MinMaxScaler()\nencoding_data = scaler.fit_transform(encoding)\nencoding = pd.DataFrame(encoding_data, columns=encoding.columns)\nencoding[\" \"] = [0, 0, 0, 0, 0, 0, 0, 0]\nencoding","metadata":{"execution":{"iopub.status.busy":"2024-07-31T15:23:15.36949Z","iopub.execute_input":"2024-07-31T15:23:15.369744Z","iopub.status.idle":"2024-07-31T15:23:15.828301Z","shell.execute_reply.started":"2024-07-31T15:23:15.369722Z","shell.execute_reply":"2024-07-31T15:23:15.827398Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def encode_sequence(sequence):\n    data = []\n    for res in sequence:\n        encoded_res = encoding[res].values\n        data.append(encoded_res)\n    actual_length = len(sequence)\n    if max_length > actual_length:\n        for _ in range(int(max_length) - actual_length):\n            data.append(np.zeros(8))\n    return np.array(data)\n\ndef encode_column(column):\n    data = []\n    for seq in column:\n        encoded_seq = encode_sequence(seq)\n        data.append(encoded_seq)\n    return np.array(data)\nencoded_data = encode_column(dna_binding_sequences.sequence)\nencoded_data.shape","metadata":{"execution":{"iopub.status.busy":"2024-07-31T15:23:15.829594Z","iopub.execute_input":"2024-07-31T15:23:15.829959Z","iopub.status.idle":"2024-07-31T15:23:24.229185Z","shell.execute_reply.started":"2024-07-31T15:23:15.829929Z","shell.execute_reply":"2024-07-31T15:23:24.228241Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## VAE","metadata":{}},{"cell_type":"code","source":"import keras\nimport tensorflow as tf\nimport tensorflow.keras.backend as K \nimport tensorflow.keras.backend as K \nfrom tensorflow.keras import datasets\nfrom tensorflow.keras import layers, models, callbacks","metadata":{"execution":{"iopub.status.busy":"2024-07-31T15:23:24.230725Z","iopub.execute_input":"2024-07-31T15:23:24.231108Z","iopub.status.idle":"2024-07-31T15:23:35.701221Z","shell.execute_reply.started":"2024-07-31T15:23:24.231077Z","shell.execute_reply":"2024-07-31T15:23:35.700463Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"class Sampling(layers.Layer):\n    def call(self, inputs):\n        z_mean, z_log_var = inputs\n        batch = tf.shape(z_mean)[0]\n        dim = tf.shape(z_mean)[1]\n        epsilon = K.random_normal(shape = (batch, dim))\n        return z_mean + tf.exp(0.5*z_log_var) * epsilon","metadata":{"execution":{"iopub.status.busy":"2024-07-31T15:23:35.702333Z","iopub.execute_input":"2024-07-31T15:23:35.702942Z","iopub.status.idle":"2024-07-31T15:23:35.708553Z","shell.execute_reply.started":"2024-07-31T15:23:35.702916Z","shell.execute_reply":"2024-07-31T15:23:35.707591Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"vae_encoder_input = layers.Input(shape=(int(max_length), encoding_length), name=\"vae_encoder_input\")\n\nx = layers.Conv1D(64, 3, strides=2, padding=\"same\")(vae_encoder_input)\nx = layers.BatchNormalization()(x)\nx = layers.LeakyReLU()(x)\n\nx = layers.Conv1D(64, 3, strides=2, padding=\"same\")(x)\nx = layers.BatchNormalization()(x)\nx = layers.LeakyReLU()(x)\n\nx = layers.Conv1D(64, 3, strides=2, padding=\"same\")(x)\nx = layers.BatchNormalization()(x)\nx = layers.LeakyReLU()(x)\n\nx = layers.Conv1D(64, 3, strides=2, padding=\"same\")(x)\nx = layers.BatchNormalization()(x)\nx = layers.LeakyReLU()(x)\n\nx = layers.Conv1D(64, 3, strides=2, padding=\"same\")(x)\nx = layers.BatchNormalization()(x)\nx = layers.LeakyReLU()(x)\n\nshape_before_flattening = K.int_shape(x)[1:]\nx = layers.Flatten()(x)\nz_mean = layers.Dense(200, name=\"z_mean\")(x)\nz_log_var = layers.Dense(200, name=\"z_log_var\")(x)\nz = Sampling()([z_mean, z_log_var])\nencoder = models.Model(vae_encoder_input, [z_mean, z_log_var, z], name=\"vae_encoder\")\nencoder.summary()","metadata":{"execution":{"iopub.status.busy":"2024-07-31T15:23:35.709922Z","iopub.execute_input":"2024-07-31T15:23:35.710425Z","iopub.status.idle":"2024-07-31T15:23:36.537355Z","shell.execute_reply.started":"2024-07-31T15:23:35.710393Z","shell.execute_reply":"2024-07-31T15:23:36.536526Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"decoder_input = layers.Input(shape=(200,), name=\"decoder_input\")\nx = layers.Dense(np.prod(shape_before_flattening))(decoder_input)\nx = layers.BatchNormalization()(x)\nx = layers.LeakyReLU()(x)\n\nx = layers.Reshape(shape_before_flattening)(x)\nx = layers.Conv1DTranspose(64, 3, strides=2, padding=\"same\")(x)\nx = layers.BatchNormalization()(x)\nx = layers.LeakyReLU()(x)\n\nx = layers.Conv1DTranspose(64, 3, strides=2, padding=\"same\")(x)\nx = layers.BatchNormalization()(x)\nx = layers.LeakyReLU()(x)\n\nx = layers.Conv1DTranspose(64, 3, strides=2, padding=\"same\")(x)\nx = layers.BatchNormalization()(x)\nx = layers.LeakyReLU()(x)\n\nx = layers.Conv1DTranspose(64, 3, strides=2, padding=\"same\")(x)\nx = layers.BatchNormalization()(x)\nx = layers.LeakyReLU()(x)\n\nx = layers.Conv1DTranspose(64, 3, strides=2, padding=\"same\")(x)\nx = layers.BatchNormalization()(x)\nx = layers.LeakyReLU()(x)\n\ndecoder_output = layers.Conv1D(encoding_length, 3, strides=1, activation=\"sigmoid\", padding=\"same\", name=\"decoder_output\")(x)\ndecoder = models.Model(decoder_input, decoder_output)\ndecoder.summary()\n\n","metadata":{"execution":{"iopub.status.busy":"2024-07-31T15:23:36.538555Z","iopub.execute_input":"2024-07-31T15:23:36.538866Z","iopub.status.idle":"2024-07-31T15:23:36.683234Z","shell.execute_reply.started":"2024-07-31T15:23:36.538839Z","shell.execute_reply":"2024-07-31T15:23:36.682336Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"\n\nfrom keras import metrics, losses, optimizers\ntf.config.run_functions_eagerly(True)\nclass VAE(models.Model):\n    def __init__(self, encoder, decoder, **kwargs):\n        super(VAE, self).__init__(**kwargs)\n        self.encoder = encoder\n        self.decoder = decoder\n        self.total_loss_tracker = metrics.Mean(name=\"total_loss\")\n        self.reconstruction_loss_tracker = metrics.Mean(\n            name=\"reconstruction_loss\"\n        )\n        self.kl_loss_tracker = metrics.Mean(name=\"kl_loss\")\n\n    @property\n    def metrics(self):\n        return [\n            self.total_loss_tracker,\n            self.reconstruction_loss_tracker,\n            self.kl_loss_tracker,\n        ]\n\n    def call(self, inputs):\n        \"\"\"Call the model on a particular input.\"\"\"\n        z_mean, z_log_var, z = encoder(inputs)\n        reconstruction = decoder(z)\n        return z_mean, z_log_var, reconstruction\n\n    def train_step(self, data):\n        \"\"\"Step run during training.\"\"\"\n        with tf.GradientTape() as tape:\n            z_mean, z_log_var, reconstruction = self(data)\n            #Binary cross entropy loss\n            reconstruction_loss = tf.reduce_mean(\n                2000\n                * losses.binary_crossentropy(\n                    data, reconstruction, axis=(1, 2)\n                )\n            )\n            ## Kullback-Leibler divergence loss\n            kl_loss = tf.reduce_mean(\n                tf.reduce_sum(\n                    -0.5\n                    * (1 + z_log_var - tf.square(z_mean) - tf.exp(z_log_var)),\n                    axis=1,\n                )\n            )\n            #Just sum boths losses to get the total loss\n            total_loss = reconstruction_loss + kl_loss\n\n        grads = tape.gradient(total_loss, self.trainable_weights)\n        self.optimizer.apply_gradients(zip(grads, self.trainable_weights))\n\n        self.total_loss_tracker.update_state(total_loss)\n        self.reconstruction_loss_tracker.update_state(reconstruction_loss)\n        self.kl_loss_tracker.update_state(kl_loss)\n\n        return {m.name: m.result() for m in self.metrics}\n\n","metadata":{"execution":{"iopub.status.busy":"2024-07-31T15:23:36.684425Z","iopub.execute_input":"2024-07-31T15:23:36.684778Z","iopub.status.idle":"2024-07-31T15:23:36.696687Z","shell.execute_reply.started":"2024-07-31T15:23:36.684747Z","shell.execute_reply":"2024-07-31T15:23:36.695827Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"vae = VAE(encoder, decoder)","metadata":{"execution":{"iopub.status.busy":"2024-07-31T15:23:36.69799Z","iopub.execute_input":"2024-07-31T15:23:36.6984Z","iopub.status.idle":"2024-07-31T15:23:36.720457Z","shell.execute_reply.started":"2024-07-31T15:23:36.698369Z","shell.execute_reply":"2024-07-31T15:23:36.719701Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"optimizer = optimizers.Adam()\nvae.compile(optimizer=optimizer)","metadata":{"execution":{"iopub.status.busy":"2024-07-31T15:23:36.7214Z","iopub.execute_input":"2024-07-31T15:23:36.721628Z","iopub.status.idle":"2024-07-31T15:23:36.733543Z","shell.execute_reply.started":"2024-07-31T15:23:36.721608Z","shell.execute_reply":"2024-07-31T15:23:36.732722Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from keras.callbacks import *\nfilename = \"/kaggle/working/model.weights.h5\"\ncallbacks = [\n    EarlyStopping(monitor=\"total_loss\", patience=4, mode=\"min\"),\n    ReduceLROnPlateau(monitor=\"total_loss\", patience=2, mode=\"min\"),\n    ModelCheckpoint(filename, monitor=\"total_loss\", mode=\"min\",save_weights_only=True, save_best_only=True)\n]","metadata":{"execution":{"iopub.status.busy":"2024-07-31T15:23:36.73472Z","iopub.execute_input":"2024-07-31T15:23:36.735532Z","iopub.status.idle":"2024-07-31T15:23:36.740744Z","shell.execute_reply.started":"2024-07-31T15:23:36.735508Z","shell.execute_reply":"2024-07-31T15:23:36.739942Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"vae.fit(\n    encoded_data,\n    epochs=20,\n    batch_size=100,\n    shuffle=True,\n    callbacks=callbacks\n)","metadata":{"execution":{"iopub.status.busy":"2024-07-31T15:23:36.74178Z","iopub.execute_input":"2024-07-31T15:23:36.742082Z","iopub.status.idle":"2024-07-31T15:26:46.408678Z","shell.execute_reply.started":"2024-07-31T15:23:36.742053Z","shell.execute_reply":"2024-07-31T15:26:46.407533Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"vae.load_weights(filename)","metadata":{"execution":{"iopub.status.busy":"2024-07-31T15:27:07.796415Z","iopub.execute_input":"2024-07-31T15:27:07.797066Z","iopub.status.idle":"2024-07-31T15:27:07.977619Z","shell.execute_reply.started":"2024-07-31T15:27:07.797035Z","shell.execute_reply":"2024-07-31T15:27:07.976851Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from scipy.spatial.distance import cityblock\ndef decode_sequence(values):\n    sequence = \"\"\n    for row in values:\n        distances = np.array([cityblock(row, encoding[res]) for res in encoding.columns])\n        res = encoding.columns[np.argmin(distances)]\n        sequence+=res\n    return sequence\n\nfor index in range(0, 10):\n    result = vae(np.array([encoded_data[index]]))\n    encoded_sequence = result[-1][0]\n    print(\"Original sequence:\", dna_binding_sequences.sequence.values[index], sep=\"\\n\")\n    print(\"Reconstructed sequence:\",decode_sequence(encoded_sequence), sep=\"\\n\")\n    print()","metadata":{"execution":{"iopub.status.busy":"2024-07-31T15:27:09.852757Z","iopub.execute_input":"2024-07-31T15:27:09.853447Z","iopub.status.idle":"2024-07-31T15:27:17.569142Z","shell.execute_reply.started":"2024-07-31T15:27:09.853414Z","shell.execute_reply":"2024-07-31T15:27:17.568198Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}