{"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":"markdown","source":"<a id = \"contents\"></a>\n## Table of Contents\n\n- #### [Introduction](#intro)\n- #### [RNA](#RNA)\n    - [What does RNA do in simple terms?](#RNA_do)\n    - [What is the structure of RNA?](#RNA_structure)\n    - [What are the important fundamental properties of RNA structure?](#RNA_properties)\n    - [Why is a distinction drawn between RNA's primary and secondary structures?](#RNA_distinction)\n    - [Can an RNA sequence have the same primary structure and different secondary structures?](#RNA_diff)\n    - [What is the significance of RNA’s secondary structure?](#RNA_sec)    \n- #### [Initial EDA](#init_EDA)\n    - [How can the training data be used to build a model that predicts the structure of a RNA molecule?](#trng)\n    - [How can you visualize a reactivity profile?](#react_viz)\n    - [Setting up an RNA Science Environment](#RNA_env)\n    - [Starters EDA + Fastai + RNN](#EDA_start)\n        - [An introduction to Codons](#codon)\n            - [History](#cdn_hist)\n            - [Analysis](#cdn_anal)\n            - [Frequency plots](#cdn_freq)\n            - [Heat map](#cdn_heat)\n            - [Nucleotide level heatmap](#nuc_heat)\n            - [Composition of A,U,G,C Nucleotides](#nuc_comp)\n            - [Bias due to Codon Pair Interaction](#nuc_bias)\n            - [3D Scatter Plot of Codon Counts for AUG, GCA, and GCU Codons](#nuc_3d)\n        - [Embeddings (Word2Vec)](#embed)\n        - [An approach using average Reactivity](#react)\n        - [An approach using Seq-to-Seq RNN](#rnn)\n    - [RNA Starter [0.186 LB]](#RNA_LB)\n    - [RNA Starter submission [.0186 LB]](#RNA_sub)","metadata":{}},{"cell_type":"markdown","source":"<a id = \"intro\"></a>\n## Introduction\n\nI like Kaggle competitions that address interesting and challenging areas about which I know very little. This competition fits the bill very nicely\n\nChatGPT has been used to understand the background to the competition better. Questions in bold were posed\n\n[Return to Contents](#contents)","metadata":{}},{"cell_type":"markdown","source":"<a id = \"RNA\"></a>\n## RNA","metadata":{}},{"cell_type":"markdown","source":"<a id = \"RNA_do\"></a>\n## What does RNA do in simple terms?\n\nImagine your body is like a big factory. In this factory, there are instruction manuals that tell the machines how to make various products.\n\nDNA is like the master instruction manual stored in a safe vault (the nucleus of the cell). But you don’t want to risk taking the master manual out every time you need to make something. So, you make a copy of the relevant page and give it to the factory floor.\n\nRNA is like that copied page. It carries the instructions from the master manual (DNA) to the factory floor (the rest of the cell) where the products (proteins and other molecules) are made.\n\nIn short, RNA helps read the instructions from our DNA and ensure that our body's machinery makes the right products.\n\n[Return to Contents](#contents)","metadata":{}},{"cell_type":"markdown","source":"<a id = \"RNA_structure\"></a>\n## What is the structure of RNA?\n\n[Link to LibreTexts Structure and Function of DNA](https://bio.libretexts.org/Bookshelves/Microbiology/Microbiology_(OpenStax)/10%3A_Biochemistry_of_the_Genome/10.03%3A_Structure_and_Function_of_RNA)\n\nRNA is typically single stranded and is made of ribonucleotides that are linked by phosphodiester bonds.\n\n### What is a ribonucleotide?\n\nA ribonucleotide is the basic building block of RNA (ribonucleic acid). It consists of three components:\n\n1. **A ribose sugar**: Unlike the deoxyribose sugar found in DNA (deoxyribonucleic acid), ribose has a hydroxyl (-OH) group attached to the 2' carbon atom.\n\n2. **A phosphate group**: This group is bound to the 5' carbon atom of the ribose sugar. In RNA chains, ribonucleotides are linked together via these phosphate groups, which form the backbone of the RNA strand.\n\n3. **A nitrogenous base**: There are four primary nitrogenous bases in RNA – adenine (A), cytosine (C), guanine (G), and uracil (U). In DNA, uracil is replaced by thymine (T).\n\nWhen ribonucleotides link together via phosphodiester bonds between their phosphate group and the 3' OH group of the ribose sugar of the next ribonucleotide, they form the RNA polymer. The sequence of the nitrogenous bases in this RNA chain carries the genetic information necessary for various cellular functions, especially in the synthesis of proteins.\n\n### What is a phosphodiester bond?\n\nA phosphodiester bond is a type of chemical bond that plays a fundamental role in the structure of nucleic acids, including DNA (deoxyribonucleic acid) and RNA (ribonucleic acid). These bonds are responsible for linking together the individual nucleotide units within the nucleic acid polymer.\n\nHere's a breakdown of what a phosphodiester bond is and how it forms:\n\n1. **Nucleotides**: Nucleic acids, such as DNA and RNA, are composed of long chains (polymers) of smaller subunits called nucleotides. Each nucleotide consists of three main components: a sugar molecule, a phosphate group, and a nitrogenous base.\n\n2. **Sugar Molecule**: In DNA, the sugar is deoxyribose, and in RNA, it is ribose. These sugars have a backbone-like structure and serve as the structural framework for the nucleic acid.\n\n3. **Phosphate Group**: Attached to the sugar molecule is a phosphate group. This phosphate group consists of a central phosphorus atom (P) bonded to oxygen atoms (O). Two of the oxygen atoms are single-bonded to the sugar's carbon atoms, while a third oxygen atom is double-bonded to one of the carbon atoms. The oxygen atom that is double-bonded to the carbon atom is also bonded to the phosphorus atom, forming a phosphodiester linkage.\n\n4. **Phosphodiester Linkage**: A phosphodiester bond is formed when the phosphate group of one nucleotide reacts with the hydroxyl (OH) group on the sugar molecule of another nucleotide. This reaction involves the removal of a water molecule (dehydration or condensation reaction) and the formation of a covalent bond between the phosphate and sugar molecules. The result is a covalent linkage between the 3' carbon of one sugar molecule and the 5' carbon of the next sugar molecule in the chain.\n\n5. **Backbone Formation**: Repeated formation of phosphodiester bonds between adjacent nucleotides results in the formation of the sugar-phosphate backbone of the nucleic acid. The nitrogenous bases project outward from this backbone and provide the sequence information.\n\nIn summary, a phosphodiester bond is a covalent bond that joins the sugar molecules of adjacent nucleotides in a nucleic acid chain (DNA or RNA). These bonds create a continuous and stable backbone structure for the molecule, with the nitrogenous bases extending from this backbone to encode genetic information. The sequence of nucleotides linked by phosphodiester bonds carries the genetic code and is essential for the storage and transmission of genetic information in living organisms.\n\n[Return to Contents](#contents)","metadata":{}},{"cell_type":"markdown","source":"<a id = \"RNA_properties\"></a>\n## What are the important fundamental properties of RNA structure?\n\nRNA (ribonucleic acid) is a molecule that plays a crucial role in various cellular processes, including protein synthesis, gene regulation, and information transfer. The structure of RNA is fundamental to its function, and several important properties characterize RNA structure:\n\n1. **Nucleotides**: RNA is composed of a chain of nucleotides, which are the building blocks of the molecule. Each nucleotide consists of three components: a sugar molecule (ribose in RNA), a phosphate group, and a nitrogenous base. The four nitrogenous bases found in RNA are adenine (A), uracil (U), cytosine (C), and guanine (G).\n\n2. **Single-Stranded**: Unlike DNA (deoxyribonucleic acid), which typically forms a double helix, RNA is usually single-stranded. However, it can form secondary structures through intramolecular base-pairing, leading to complex three-dimensional structures.\n\n3. **Base Pairing**: In RNA, base pairing occurs between complementary bases, just as it does in DNA. Adenine (A) pairs with uracil (U), and cytosine (C) pairs with guanine (G). This base pairing is essential for the formation of secondary structures such as stem-loops and hairpins.\n\n4. **Secondary Structure**: RNA can fold back on itself, creating secondary structures through hydrogen bonding between complementary bases. Common secondary structures in RNA include stem-loops, bulges, internal loops, and pseudoknots. These structures are crucial for functions like RNA folding, stability, and interaction with other molecules.\n\n5. **Tertiary Structure**: In addition to secondary structures, RNA can also adopt complex three-dimensional shapes through interactions between distant parts of the molecule. These tertiary structures are essential for many RNA functions, such as ribozyme catalysis and binding to proteins or other molecules.\n\n6. **Base Stacking**: Stacking interactions between adjacent bases contribute to the stability and compactness of RNA structures. These interactions help RNA molecules maintain their three-dimensional shapes.\n\n7. **Functional Groups**: The ribose sugar in RNA contains a hydroxyl group (OH) at the 2' position, which is not present in DNA. This hydroxyl group can participate in chemical reactions, making RNA more chemically versatile than DNA.\n\n8. **Base Modifications**: RNA can undergo various chemical modifications to its bases. These modifications can affect RNA stability, structure, and function. Examples include methylation and pseudouridylation.\n\n9. **Non-canonical Base Pairing**: RNA can form non-canonical base pairs, which are interactions between bases other than the standard A-U and C-G pairs. These non-canonical base pairs can contribute to the diversity and stability of RNA structures.\n\n10. **Flexibility**: RNA is a flexible molecule that can adopt different conformations depending on its sequence and environmental conditions. This flexibility is essential for its various functions, including recognition and binding to other molecules.\n\nUnderstanding these fundamental properties of RNA structure is crucial for unraveling the diverse roles that RNA plays in cells, including its involvement in gene expression, RNA interference, and catalysis by ribozymes.\n\nSee also:\n\n[Anna Marie Pyle (Yale U./HHMI) Part 1: RNA Structure](https://www.youtube.com/watch?v=WCrlm18KQ48)  \n[Anna Marie Pyle (Yale U./HHMI) Part 2: Inside an RNA Splicing Machine](https://www.youtube.com/watch?v=ESXo3fTThBI)  \n[Anna Marie Pyle (Yale U./HHMI) Part 3: RNA Helicases and RNA-triggered Signaling Proteins](https://www.youtube.com/watch?v=LjCDqL8n5F0)  \n\n[Return to Contents](#contents)","metadata":{}},{"cell_type":"markdown","source":"<a id = \"RNA_distinction\"></a>\n## Why is a distinction drawn between RNA's primary and secondary structures?\n\nThe distinction between primary and secondary structures (as well as tertiary and quaternary structures) in biomolecules, like RNA and proteins, is important because it helps describe and differentiate levels of organization and complexity in these molecules. Each level provides a different perspective and understanding of the molecule's properties and functions. Here's why distinguishing between RNA's primary and secondary structures is particularly important:\n\n1. **Definition and Simplicity**:\n\n    - Primary Structure refers to the linear sequence of nucleotides in the RNA molecule, connected by phosphodiester bonds. It is the simplest level of organization and dictates the higher order structures that the RNA molecule can adopt.\n    \n    - Secondary Structure refers to the local spatial arrangement of the RNA molecule due to base pairing (usually through hydrogen bonding) between complementary nucleotides. \n    \n    Motifs, in a biological context, refer to recurring, recognizable patterns or elements within molecules that usually have a specific structure and function. Common motifs in RNA secondary structure include stems (or helices), loops, hairpins, bulges, and internal loops.  They can be found in both DNA/RNA sequences and protein structures. Motifs are important because their recurrence in different molecules or within the same molecule often indicates a conserved function. Identifying and studying these motifs can provide insights into the molecular mechanisms and evolution of biological processes. Furthermore, in bioinformatics and computational biology, detecting and analyzing motifs can be crucial for predicting the function of newly discovered molecules or understanding regulatory networks within cells.\n    \n2. **Determination and Prediction**:\n\n    While the primary structure is determined directly by sequencing the RNA molecule, the secondary structure can be predicted based on thermodynamic principles and is sometimes confirmed through experimental techniques like X-ray crystallography, nuclear magnetic resonance (NMR) spectroscopy, or cryo-electron microscopy.\n\n3. **Functional Implications**:\n\n    The primary structure determines the identity and order of nucleotides, but it's the secondary structure (and higher-level structures) that often gives the RNA its functional properties. For instance, the ability of tRNA to bring amino acids to the ribosome is deeply intertwined with its cloverleaf secondary structure.\n\n4. **Evolutionary Conservation**:\n\n    While primary structures can vary between related RNA molecules across species, secondary structures (or the motifs they form) might be conserved, indicating their functional importance.\n\n5. **Stability and Dynamics**:\n\n    The primary structure provides the potential for various secondary structures. However, the actual secondary structure adopted by an RNA molecule in a cell can depend on factors like the local environment, interactions with other molecules, and cellular conditions. These secondary structures can be dynamic, with the RNA molecule shifting between different conformations based on cellular needs.\n\n6. **Framework for Further Complexity**:\n\n    Secondary structures serve as the foundation for tertiary structures, where the overall three-dimensional shape of the RNA molecule is formed through interactions between different parts of the secondary structure.\n\nIn essence, the distinction between primary and secondary structures provides a hierarchical framework to understand the architecture and function of RNA molecules. It allows scientists to discuss, analyze, and predict how changes in the nucleotide sequence (primary structure) might affect the molecule's shape and function by altering its secondary and higher-order structures.\n\n[Return to Contents](#contents)","metadata":{}},{"cell_type":"markdown","source":"<a id =  \"RNA_diff\"></a>\n## Can an RNA sequence have the same primary structure and different secondary structures?\n\nYes, an RNA sequence can have the same primary structure (the linear sequence of nucleotides) and adopt different secondary structures under various conditions. The secondary structure of RNA refers to the local folding of the molecule into helices, loops, bulges, and junctions, largely driven by base-pairing interactions, particularly between adenine-uracil (A-U) and guanine-cytosine (G-C) pairs.\n\nSeveral factors can influence the secondary structure of a given RNA sequence:\n\n1. **Ionic Conditions**: Changes in the ionic composition of the solution, particularly the concentration of divalent cations like magnesium, can significantly affect RNA folding.\n\n2. **Temperature**: Heating and cooling can shift the equilibrium between different folded states.\n\n3. **Molecular Crowding**: The presence of other macromolecules can affect RNA folding through excluded volume effects.\n\n4. **RNA Concentration**: At higher concentrations, RNA molecules might form dimers or higher-order structures, leading to alternative folding patterns.\n\n5. **Post-transcriptional Modifications**: Modifications such as methylation can change the folding propensity of RNA.\n\n6. **Binding Partners**: The presence of proteins, small molecules, or other RNA molecules can stabilize certain conformations over others.\n\n7. **pH**: Changes in pH can affect protonation states of nucleotides, influencing their hydrogen bonding patterns.\n\n8. **Mutational Changes**: Though the primary sequence remains the same, single-nucleotide polymorphisms (SNPs) or post-transcriptional edits can change a nucleotide's identity, potentially altering the secondary structure.\n\nDue to these variables, RNA molecules can exhibit structural plasticity and exist in an ensemble of conformations that are in dynamic equilibrium. The ability of RNA to adopt multiple conformations is also fundamental to its role in the cell, where its function can be regulated by changes in structure, allowing it to respond to cellular signals and environmental changes.\n\nMoreover, RNA molecules, especially ribozymes and regulatory RNAs, can undergo conformational changes as part of their functional cycle. For example, riboswitches change their secondary structure in response to the binding of small molecules, which in turn affects gene expression.\n\nRNA folding is often studied using techniques such as nuclear magnetic resonance (NMR) spectroscopy, X-ray crystallography, cryo-electron microscopy, and various biochemical assays to understand the conditions and dynamics of these structural transitions.\n\n[Return to Contents](#contents)","metadata":{}},{"cell_type":"markdown","source":"<a id =  \"RNA_sec\"></a>\n## What is the significance of RNA’s secondary structure?\n\nRNA's secondary structure holds immense significance in biology for several reasons:\n\n1. **Functionality**: RNA molecules, unlike DNA, often play active roles in cellular processes. The functionality of many RNA molecules is closely tied to their secondary structures. For example, the well-defined secondary structures of tRNA molecules are crucial for their role in protein synthesis.\n\n2. **Molecular Recognition**: RNA's secondary structure provides a platform for specific molecular interactions. The ribosome, which is largely composed of ribosomal RNA (rRNA), utilizes its structured RNA components to recognize and interact with messenger RNA (mRNA) and transfer RNA (tRNA) during protein synthesis.\n\n3. **Catalytic Activity**: Some RNA molecules can act as enzymes, called ribozymes. The catalytic activity of ribozymes is often directly related to their secondary and tertiary structures. A famous example is the self-splicing intron.\n\n4. **Regulation of Gene Expression**: In both prokaryotes and eukaryotes, RNA secondary structures can play roles in regulating gene expression. For instance:\n\n    - Riboswitches are segments of mRNA molecules that can bind small molecules, leading to structural changes in the RNA. These structural changes can regulate various processes, such as transcription or translation.\n    - RNA interference (RNAi) mechanisms involve small RNA molecules, like siRNA or miRNA, which can base-pair with target mRNA molecules, leading to their degradation or blocking their translation. The secondary structure of these RNA molecules can influence their stability, loading into protein complexes, and specificity.\n    \n    \n5. **Viruses**: In many viruses, RNA serves as the genetic material. The secondary structures in viral RNAs can be critical for their replication, packaging, and infectivity.\n\n6. **Evolutionary Conservation**: Secondary structures are often conserved in evolution, even when the primary sequence can vary. This conservation indicates the importance of the structure in the function of the RNA molecule.\n\n7. **Alternative Structures & Plasticity**: Some RNA molecules can adopt multiple stable secondary structures. This structural plasticity can be harnessed for regulatory purposes. For instance, an RNA molecule might adopt one structure under certain conditions and a different structure under others, leading to different functional outcomes.\n\n8. **Experimental Relevance**: Understanding RNA secondary structure is crucial in many experimental procedures. For instance, **when predicting the potential function or interactions of a novel RNA molecule, researchers will often first predict its secondary structure**. Similarly, targeting RNA with drugs or other molecules often requires a detailed understanding of its structure.\n\n9. **Technological Applications**: RNA's ability to fold into specific structures has been harnessed in the field of nanotechnology. Researchers can design RNA molecules to fold into desired shapes, which can then serve as scaffolds for building nano-scale devices or structures.\n\nIn summary, RNA's secondary structure is integral to its function, and understanding these structures is crucial for a comprehensive grasp of molecular biology, cellular regulation, evolution, and biotechnology.\n\n[Return to Contents](#contents)\n","metadata":{}},{"cell_type":"markdown","source":"<a id =  \"init_EDA\"></a>\n## Initial EDA\n[Return to Contents](#contents)","metadata":{}},{"cell_type":"markdown","source":"<a id =  \"trng\"></a>\n## How can the training data be used to build a model that predicts the structure of a RNA molecule?\n(There is one experimental profile per row)\n\n**sequence_id** - (string) An arbitrary identifier like 8cdfeef009ea for each sequence.\n\n**sequence** - (string) Describes the RNA sequence, a string of A, C, G, and U.\n\n**experiment_type** - (string) Either DMS_MaP or 2A3_MaP to describe the type of chemical mapping experiment that was used to generate each profile. References: DMS, 2A3.\n\n**dataset_name** - (string) arbitrary name of high throughput sequencing dataset from which the reactivity profile was extracted.\n\n**reads** - (integer) Number of reads in the high throughput sequencing experiment that were assigned to the RNA sequence, and whose mutations were tabulated to compile the reactivity profile. (These values do not need to be predicted in this competition.) \n\n**reactivity profile**. (These values do not need to be predicted in this competition.) \n\n**signal_to_noise** - (float) Signal/noise value for the profile, defined as mean( measurement value over probed nts )/mean( statistical error in measurement value over probed nts). (These values do not need to be predicted in this competition.)\n\n**SN_filter** - (Boolean) 0 or 1 depending on whether the profile has signal_to_noise≥ 1.00 and reads≥ 100. For evaluation, only sequences whose DMS_MaP and 2A3_MaP profiles both pass this filter will be used to score submissions. (These values do not need to be predicted in this competition.)\n\n**reactivity_0001, reactivity_0002,…** - (float) An array of floating point numbers of the train data, should have the same length as the RNA sequence, which defines the reactivity profile for the RNA. These are the type of data that need to be predicted in this competition. For sequences shorter than the maximum RNA length, positions that go beyond the sequence length have null. Several positions near the beginning and end of the sequence also cannot be probed due to technical reasons, and their reactivity values are null.\n\n**reactivity_error_0001, reactivity_error_0002,…** - (float) An array of floating point numbers, should have the same length as the corresponding reactivity_ columns, calculated errors in experimental values obtained in reactivity derived from counting statistics in the high-throughput sequencing experiment. (These values do not need to be predicted in this competition.)\n\n\n**Given the dataset you've described, predicting the reactivity profile of RNA sequences is essentially a sequence regression problem**. The reactivity profile essentially reflects how each nucleotide in an RNA molecule interacts with its environment. Predicting this can give insights into the structure and function of the RNA.\nHere's a suggested approach to predict the reactivity profile of RNA sequences using this data:\n\n1. **Data Preprocessing**:\n    - One-Hot Encoding: Encode the RNA sequences (A, C, G, U) as one-hot vectors. For a sequence like \"AC\", the one-hot encoding would be [1, 0, 0, 0], [0, 1, 0, 0].\n    - Padding Sequences: Since the sequences can have different lengths, pad the sequences and their corresponding reactivity profiles to a consistent length (using the maximum sequence length from the dataset).\n    - Normalize Data: Normalize continuous data columns (like `signal_to_noise`) for better model performance.\n    - Filtering: Use only rows where `SN_filter` is 1 for training, as these rows pass the signal-to-noise and reads threshold.\n    \n    \n2. **Model Architecture**:\n    - Input Layer: Start with an input layer for the one-hot encoded sequences.\n    - Embedding Layer (optional): Though we have one-hot encoded sequences, embedding layers can still capture intricate relationships.\n    - Recurrent Layers: LSTMs or GRUs are suitable due to their capability to handle sequential data. Multiple stacked layers might improve performance.\n    - Dense Layers: Add one or more dense layers to capture non-linear relationships.\n    - Output Layer: This should be a layer with the same length as the maximum sequence length, predicting the reactivity for each nucleotide in the sequence.\n    \n    \n3. **Incorporate Experimental Data**:\n    - Consider using `experiment_type` and `dataset_name` as additional inputs. One can one-hot encode these categorical variables and merge them with the main sequence input, possibly after a dense layer. \n    \n    \n4. **Training**:\n    - Loss Function: Mean squared error or mean absolute error are appropriate loss functions for regression tasks.\n    - Validation: Use a subset of the data to validate the model's performance during training. Monitor the validation loss to detect overfitting.\n    - Regularization: To prevent overfitting, consider adding dropout layers in your network.\n    \n    \n5. **Post-Processing**:\n    - For sequences shorter than the maximum RNA length, positions that go beyond the sequence length should be ignored or set to null.\n    - Reactivity values near the beginning and end of sequences, which are null due to technical reasons, should be handled accordingly in predictions.\n    \n    \n6. **Evaluation**:\n    - Compare the predicted reactivity profiles with the actual profiles using metrics like RMSE (Root Mean Squared Error) or MAE (Mean Absolute Error) to gauge the performance of your model.\n    \nRemember, model architecture and hyperparameters would need tuning. Using techniques like cross-validation, Bayesian optimization, or grid search can be useful.\n\n[Return to Contents](#contents)","metadata":{}},{"cell_type":"markdown","source":"<a id =  \"react_viz\"></a>\n## How can you visualize a reactivity profile?\n\nVisualizing a reactivity profile typically involves generating a plot where the x-axis represents the nucleotide position along the RNA sequence, and the y-axis represents the reactivity value. There are several ways and tools that can be used to visualize a reactivity profile:\n\n1. **Basic Line Plot**: A simple line plot can effectively represent reactivity values across an RNA sequence. Each nucleotide position is plotted along the x-axis, with its corresponding reactivity value plotted on the y-axis.\n\n2. **Heatmap**: For longer sequences or when comparing multiple profiles, a heatmap can be useful. In this representation, the RNA sequence is usually on the x-axis, different conditions or experiments on the y-axis, and reactivity values are represented by color intensity.\n\n3. **Overlay on Secondary Structure**: If the secondary structure of the RNA is known or predicted, reactivity values can be overlaid on this structure. This provides a visual representation of which regions (e.g., loops, stems, bulges) are more or less reactive. \n\n4. **Software and Tools**:\n   - **RNAstructure**: This is a suite of tools for RNA sequence analysis and structure prediction. It can incorporate SHAPE reactivity data and visualize it alongside predicted secondary structures.\n   - **RNAStructViz**: It's a tool designed for visualization and comparison of RNA secondary structures.\n   - **R and Python libraries**: Both R and Python have extensive plotting libraries (like ggplot2 for R and Matplotlib or Seaborn for Python) that can be used to generate plots for reactivity profiles.\n\n5. **Interactive Web Tools**: Some online platforms allow for interactive visualization of RNA reactivity data. These can be particularly useful for exploration and sharing data with collaborators.\n\n6. **Annotations**: When visualizing a reactivity profile, it's often helpful to annotate regions of interest, such as known or predicted structural motifs, functional domains, or binding sites. This provides context and can help in the interpretation of the data.\n\nWhen visualizing reactivity profiles, the goal is to present the data in a way that facilitates understanding of RNA structure and dynamics. By correlating reactivity data with other information (like sequence motifs or known structures), one can gain deeper insights into RNA function and behavior.\n\n### Interactive web tools\n\nThere are several interactive web tools that can be used to visualize RNA structures, reactivity profiles, and related data. Here are some notable ones:\n\n1. **RNAfold Web Server**: This is part of the ViennaRNA Web Services. It allows users to predict RNA secondary structures and can incorporate reactivity data. The visualizations are interactive, enabling users to explore predicted structures in detail.\n   - [Link to RNAfold](http://rna.tbi.univie.ac.at/cgi-bin/RNAWebSuite/RNAfold.cgi)\n\n2. **RNAstructure Web Server**: This tool predicts RNA secondary structures with the option to include SHAPE reactivity data. While the main software suite is downloadable, they also provide a web interface for some of the tools.\n   - [Link to RNAstructure Web Server](https://rna.urmc.rochester.edu/RNAstructureWeb/)\n\n3. **VARNA**: While primarily a standalone application, VARNA also offers a web interface. It's a tool for the automated drawing, visualization, and annotation of RNA secondary structures.\n   - [Link to VARNA](http://varna.lri.fr/)\n\n4. **RNALive**: This is an interactive web server for modeling RNA 3D structures. It allows users to generate 3D models based on sequence and pairing constraints.\n   - [Link to RNALive](https://rnajournal.cshlp.org/content/22/11/1629.full)\n\n5. **RiboVision**: RiboVision is a visualization and analysis tool for the simultaneous display of multiple layers of diverse information on primary (1D), secondary (2D), and three-dimensional (3D) structures of ribosomes. While it's focused on ribosomes, it's an excellent example of interactive visualization for complex RNA structures.\n   - [Link to RiboVision](http://apollo.chemistry.gatech.edu/RiboVision/)\n\n6. **RNA-Puzzles**: This is a platform for the evaluation of RNA three-dimensional structure prediction. It includes tools for visualization and comparison of predicted structures with experimental ones.\n   - [Link to RNA-Puzzles](http://rnasociety.org/rna-puzzles/)\n\nWhen using these tools, it's essential to note the source of the reactivity data and any potential limitations or assumptions inherent in the software. Always consult the original publications or documentation associated with each tool for a full understanding of its capabilities and best practices.\n\n[Return to Contents](#contents)","metadata":{}},{"cell_type":"markdown","source":"## Data Files","metadata":{}},{"cell_type":"code","source":"import gc\nimport numpy as np\nimport pandas as pd\nimport seaborn as sns\nimport matplotlib.pyplot as plt\nfrom collections import Counter, defaultdict\nimport plotly.express as px\nfrom gensim.models import Word2Vec\nfrom fastai.tabular.all import *\nfrom tensorflow.keras.preprocessing.sequence import pad_sequences\nimport torch.nn as nn\nimport torch.optim as optim\nfrom torch.utils.data import DataLoader, TensorDataset, random_split\nimport matplotlib.patches as mpatches\n\nsns.set()","metadata":{"execution":{"iopub.status.busy":"2023-10-29T12:44:39.185389Z","iopub.execute_input":"2023-10-29T12:44:39.185744Z","iopub.status.idle":"2023-10-29T12:44:58.624496Z","shell.execute_reply.started":"2023-10-29T12:44:39.185717Z","shell.execute_reply":"2023-10-29T12:44:58.623526Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### train_data.csv","metadata":{}},{"cell_type":"code","source":"train_data = pd.read_csv('/kaggle/input/stanford-ribonanza-rna-folding/train_data.csv', nrows = 5)\ntrain_data.head()","metadata":{"execution":{"iopub.status.busy":"2023-10-29T12:46:16.403749Z","iopub.execute_input":"2023-10-29T12:46:16.404564Z","iopub.status.idle":"2023-10-29T12:46:16.490932Z","shell.execute_reply.started":"2023-10-29T12:46:16.404529Z","shell.execute_reply":"2023-10-29T12:46:16.489577Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df = pd.DataFrame(train_data)\n\n# Extract the column names and their values for the first row\nfirst_row_values = df.iloc[0].to_dict()\n\n# Print the results\nfor column, value in first_row_values.items():\n    print(f\"{column}: {value}\")\n   \n#for col in df.columns:\n#    print(col, \":\", df[col].head(5).tolist())","metadata":{"execution":{"iopub.status.busy":"2023-10-29T13:10:38.602322Z","iopub.execute_input":"2023-10-29T13:10:38.602745Z","iopub.status.idle":"2023-10-29T13:10:38.614237Z","shell.execute_reply.started":"2023-10-29T13:10:38.602716Z","shell.execute_reply":"2023-10-29T13:10:38.612997Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"<a id =  \"RNA_env\"></a>\n## Setting up an RNA Science Environment\n### (This is based on [RNA Science Computational Environment](https://www.kaggle.com/code/brainbowrna/rna-science-computational-environment) by [Thomas](https://www.kaggle.com/brainbowrna))\n\nThe `Arnie` Python library is a tool for RNA sequence design and free energy calculation. It provides a unified interface to several RNA structure prediction packages. Typically, with RNA design, the goal is to predict how an RNA sequence will fold or to design a sequence that will fold into a desired structure. Arnie facilitates this by giving users a common way to interact with multiple backend prediction algorithms.\n\nSome functionalities you might find in such a library include:\n1. **Free Energy Calculations**: Predicting the stability of a given RNA structure.\n2. **Secondary Structure Prediction**: Determining how an RNA sequence will fold.\n3. **Sequence Design**: Finding an RNA sequence that will fold into a specified structure.\n\nThe `arnie` library in Python is used for predicting RNA secondary structures and evaluating free energies of RNA molecules. It provides a convenient interface to various RNA folding tools and algorithms. By integrating these tools, `arnie` facilitates the task of predicting RNA secondary structures for researchers and developers working in the RNA field.\n\nSome key features include:\n\n1. **Unified Interface**: It allows users to seamlessly switch between different folding algorithms without having to handle the nuances of each tool separately.  \n2. **Easy Installation**: `arnie` is designed to simplify the installation process for many of the integrated RNA folding tools, which sometimes can be challenging to install on their own.\n3. **Support for Multiple Algorithms**: The library provides support for several popular RNA folding algorithms, such as ViennaRNA, NUPACK, and CONTRAfold.\n4. **Utilities**: It also offers some utility functions for handling RNA sequences, such as reverse complementation, format conversion, etc.\n\nResearchers and developers might use `arnie` to predict RNA structures, analyze RNA sequences, or even develop new algorithms in the domain of RNA bioinformatics.\n\nTo get started or for detailed information, users would usually consult the official documentation, Github repository, or other resources associated with the library. If you're planning to work with RNA structures, `Arnie` and similar tools can be invaluable.","metadata":{}},{"cell_type":"code","source":"!pip install arnie","metadata":{"execution":{"iopub.status.busy":"2023-10-29T12:11:21.186716Z","iopub.execute_input":"2023-10-29T12:11:21.187345Z","iopub.status.idle":"2023-10-29T12:11:33.722166Z","shell.execute_reply.started":"2023-10-29T12:11:21.187314Z","shell.execute_reply":"2023-10-29T12:11:33.720750Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Arnie needs at least one secondary structure predictor, so let's install EternaFold. [Eternafold](https://www.nature.com/articles/s41592-022-01605-0) is a leading prediction package that was trained using sequences collected via the citizen science game [Eterna](http://eternagame.org).\n\n`EternaFold` is an RNA secondary structure prediction algorithm that originated from the Eterna project. The Eterna project is an online platform where players can design RNA sequences to fold into specific structures. The goal is to solve real-world biochemistry challenges through crowdsourcing. As players on the platform create and test RNA designs, the accumulated data is then used to refine RNA structure prediction algorithms.\n\nHere's a brief overview of EternaFold:\n\n1. **Data-Driven**: Unlike some other RNA folding algorithms that are mainly based on physical models or thermodynamics, EternaFold is driven by the data gathered from the Eterna game. This means that it incorporates the collective intelligence and insights from thousands of players.\n\n2. **Refinement Over Traditional Models**: EternaFold tweaks and refines the energy model parameters based on the results of the RNA designs made by Eterna players. This allows it to potentially make better predictions than some other algorithms in certain scenarios.\n\n3. **Collaboration with ViennaRNA**: At its core, EternaFold leverages the ViennaRNA package but modifies its energy parameters based on Eterna game data.\n\n    The ViennaRNA package is a collection of software tools for the prediction and comparison of RNA secondary structures. Developed and maintained by the Institute for Theoretical Chemistry at the University of Vienna, it is one of the most widely used packages for RNA structure prediction and analysis. The tools within the package are based on the principles of thermodynamics, providing energy minimization methods to predict the most stable structures of RNA sequences.\n\n    Key features and functionalities of the ViennaRNA package include:\n\n    1. **RNA Secondary Structure Prediction**: The primary use of the package is to predict the secondary structure of a given RNA sequence based on minimum free energy (MFE) principles.\n\n    2. **Suboptimal Folding**: Beyond just the MFE structure, the package can predict suboptimal foldings that are close in energy to the most stable structure.\n\n    3. **Partition Function and Base Pairing Probabilities**: The package can calculate the partition function and the base pairing probabilities for RNA sequences, offering insights into the ensemble of possible structures.\n\n    4. **Comparative Structure Prediction**: By leveraging multiple sequence alignments, the package can predict a consensus structure for a set of related RNA sequences.\n\n    5. **Structure Comparisons**: Tools are provided to compare different predicted structures and assess their similarity.\n\n    6. **Local Structure Prediction**: The package can predict local secondary structures within a longer RNA sequence.\n\n    7. **Design Tools**: The package includes utilities for inverse folding, which means designing RNA sequences that will fold into a desired structure.\n\n    8. **Other Utilities**: The package comes with various additional tools for tasks such as sequence alignment, structure visualization, and more.\n\n    The ViennaRNA package is highly influential in the field of RNA bioinformatics and is a staple tool for researchers working on RNA structure and function. The algorithms it uses, rooted in thermodynamic principles, provide a robust and reliable way to predict and analyze RNA structures.\n\n4. **Research and Applications**: The insights from EternaFold and the Eterna game have been used in research to improve our understanding of RNA folding and function. Additionally, the results have applications in the design of RNA-based therapeutics and diagnostics.\n\nIn essence, EternaFold represents an exciting melding of community-driven science, gamification, and bioinformatics, leading to advancements in the field of RNA research.\n\nEterna players provided many of the sequences in the data for this competition.","metadata":{}},{"cell_type":"code","source":"# Install Eternafold\n!conda config --set auto_update_conda false\n!conda install -c bioconda eternafold --yes","metadata":{"execution":{"iopub.status.busy":"2023-10-29T12:11:43.370398Z","iopub.execute_input":"2023-10-29T12:11:43.370766Z","iopub.status.idle":"2023-10-29T12:13:49.206980Z","shell.execute_reply.started":"2023-10-29T12:11:43.370737Z","shell.execute_reply":"2023-10-29T12:13:49.205598Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Ordinarily, the EternaFold conda package will automatically set necessary environment variables, but Kaggle's conda install works a little differently. Let's set them manually here using `%env`.","metadata":{}},{"cell_type":"code","source":"%env ETERNAFOLD_PATH=/opt/conda/bin/eternafold-bin\n%env ETERNAFOLD_PARAMETERS=/opt/conda/lib/eternafold-lib/parameters/EternaFoldParams.v1","metadata":{"execution":{"iopub.status.busy":"2023-10-29T12:15:14.393442Z","iopub.execute_input":"2023-10-29T12:15:14.393984Z","iopub.status.idle":"2023-10-29T12:15:14.403615Z","shell.execute_reply.started":"2023-10-29T12:15:14.393938Z","shell.execute_reply":"2023-10-29T12:15:14.402342Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Now that we have a predictor, we can make structure predictions about a given sequence. For example, let's look at an example Hammerhead ribozyme sequence. We can use arnie's `mfe`, or Minimum Free Energy, function to predict a secondary structure for this RNA sequence. The structure will be represented in \"dot-bracket\" notation, where `.` is an unpaired base and `()` represent two paired bases.","metadata":{}},{"cell_type":"code","source":"from arnie.mfe import mfe\nsequence = \"CGCUGUCUGUACUUGUAUCAGUACACUGACGAGUCCCUAAAGGACGAAACAGCG\"\n#sequence=\"GGGAACGACUCGAGUAGAGUCGAAAAACGUUGAUAUGGAUUUACUCCGAGGAGACGAACUACCACGAACAGGGGAAACUCUACCCGUGGCGUCUCCGUUUGACGAGUAAGUCCUAAGUCAACAUGCCACGCGGGUCCUUCGGGACCCGCAAAAGAAACAACAACAACAAC\"\nstructure = mfe(sequence,package=\"eternafold\")\nprint(structure)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Dot bracket notation can be a little hard to read if you're new to RNA structures. Let's visualize the structure in another way. We're going to install `draw_rna`, a Das Lab tool that let's us plot RNA structures in 2D.","metadata":{}},{"cell_type":"code","source":"!pip install draw_rna","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from draw_rna.ipynb_draw import draw_struct\ndraw_struct(sequence, structure)","metadata":{"execution":{"iopub.status.busy":"2023-10-29T12:27:32.638243Z","iopub.execute_input":"2023-10-29T12:27:32.638613Z","iopub.status.idle":"2023-10-29T12:27:33.864261Z","shell.execute_reply.started":"2023-10-29T12:27:32.638586Z","shell.execute_reply":"2023-10-29T12:27:33.863100Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from arnie.bpps import bpps\nbpps(sequence,package=\"eternafold\")","metadata":{"execution":{"iopub.status.busy":"2023-10-29T12:27:45.138349Z","iopub.execute_input":"2023-10-29T12:27:45.138781Z","iopub.status.idle":"2023-10-29T12:27:45.162749Z","shell.execute_reply.started":"2023-10-29T12:27:45.138749Z","shell.execute_reply":"2023-10-29T12:27:45.161431Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from arnie.mfe import mfe\nsequence=\"GGGAACGACUCGAGUAGAGUCGAAAAACGUUGAUAUGGAUUUACUCCGAGGAGACGAACUACCACGAACAGGGGAAACUCUACCCGUGGCGUCUCCGUUUGACGAGUAAGUCCUAAGUCAACAUGCCACGCGGGUCCUUCGGGACCCGCAAAAGAAACAACAACAACAAC\"\nstructure = mfe(sequence,package=\"eternafold\")\nprint(structure)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from draw_rna.ipynb_draw import draw_struct\ndraw_struct(sequence, structure)","metadata":{"execution":{"iopub.status.busy":"2023-10-29T12:21:00.639690Z","iopub.execute_input":"2023-10-29T12:21:00.640037Z","iopub.status.idle":"2023-10-29T12:21:03.901113Z","shell.execute_reply.started":"2023-10-29T12:21:00.640014Z","shell.execute_reply":"2023-10-29T12:21:03.899703Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from arnie.mfe import mfe\nsequence=\"GGGAACGACUCGAGUAGAGUCGAAAAUUGAUCUUUCAGUAGCUUCUACCUAUUUUUUAGUCCGCAUCUUGCAAGAUAAGACUGGCGACUUUAUGUCUACAAUUAUUACUUCCUGCCAAACUGCUGUACGCUUUGCCGUUCGCGGCAAAGAAAAGAAACAACAACAACAAC\"\nstructure = mfe(sequence,package=\"eternafold\")\nprint(structure)","metadata":{"execution":{"iopub.status.busy":"2023-10-29T12:23:21.758318Z","iopub.execute_input":"2023-10-29T12:23:21.758679Z","iopub.status.idle":"2023-10-29T12:23:21.848330Z","shell.execute_reply.started":"2023-10-29T12:23:21.758649Z","shell.execute_reply":"2023-10-29T12:23:21.847629Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from draw_rna.ipynb_draw import draw_struct\ndraw_struct(sequence, structure)","metadata":{"execution":{"iopub.status.busy":"2023-10-29T12:23:46.438758Z","iopub.execute_input":"2023-10-29T12:23:46.439267Z","iopub.status.idle":"2023-10-29T12:23:48.408138Z","shell.execute_reply.started":"2023-10-29T12:23:46.439230Z","shell.execute_reply":"2023-10-29T12:23:48.407364Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Arnie provides other functions for structure prediction. We can generate a 'Base Pair Probablility' matrix that predicts the probability of every possible base pairing (e.g, how likely is base 1 to pair with base 2, base 3, base 4...). ","metadata":{}},{"cell_type":"code","source":"from arnie.bpps import bpps\nbpps(sequence,package=\"eternafold\")","metadata":{"execution":{"iopub.status.busy":"2023-10-29T12:24:57.827418Z","iopub.execute_input":"2023-10-29T12:24:57.827743Z","iopub.status.idle":"2023-10-29T12:24:57.920140Z","shell.execute_reply.started":"2023-10-29T12:24:57.827714Z","shell.execute_reply":"2023-10-29T12:24:57.918936Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"[Return to Contents](#contents)","metadata":{}},{"cell_type":"markdown","source":"<a id =  \"EDA_start\"></a>\n## Starters EDA + Fastai + RNN\n### (This is based on [Starters EDA + Fastai + RNN](https://www.kaggle.com/code/wrecked22/starters-eda-fastai-rnn) by [Shashank Gupta](https://www.kaggle.com/wrecked22))","metadata":{}},{"cell_type":"code","source":"import gc\nimport numpy as np\nimport pandas as pd\nimport seaborn as sns\nimport matplotlib.pyplot as plt\nfrom collections import Counter, defaultdict\nimport plotly.express as px\nfrom gensim.models import Word2Vec\nfrom fastai.tabular.all import *\nfrom tensorflow.keras.preprocessing.sequence import pad_sequences\nimport torch.nn as nn\nimport torch.optim as optim\nfrom torch.utils.data import DataLoader, TensorDataset, random_split\nimport matplotlib.patches as mpatches\n\nsns.set()","metadata":{"execution":{"iopub.status.busy":"2023-10-29T10:09:30.012625Z","iopub.execute_input":"2023-10-29T10:09:30.012967Z","iopub.status.idle":"2023-10-29T10:09:52.339021Z","shell.execute_reply.started":"2023-10-29T10:09:30.012940Z","shell.execute_reply":"2023-10-29T10:09:52.338135Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"rna_sequence_data = pd.read_csv('/kaggle/input/stanford-ribonanza-rna-folding/train_data.csv', nrows = 25000)\nrna_sequence_data.head()","metadata":{"execution":{"iopub.status.busy":"2023-10-29T10:09:52.341257Z","iopub.execute_input":"2023-10-29T10:09:52.342395Z","iopub.status.idle":"2023-10-29T10:09:53.859895Z","shell.execute_reply.started":"2023-10-29T10:09:52.342351Z","shell.execute_reply":"2023-10-29T10:09:53.858870Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"rna_sequence_data = rna_sequence_data.drop(columns=rna_sequence_data.filter(like='error').columns)\nrna_sequence_data = rna_sequence_data.drop(['sequence_id', 'dataset_name', 'reads', 'SN_filter'], axis=1)\nreactivity_cols = rna_sequence_data.filter(like='reactivity').columns\nrna_sequence_data[reactivity_cols] = rna_sequence_data[reactivity_cols].clip(lower=0)\nrna_sequence_data['reactivity'] = rna_sequence_data[reactivity_cols].mean(axis=1)\nrna_sequence_data = rna_sequence_data.drop(columns=reactivity_cols)\nrna_sequence_data","metadata":{"execution":{"iopub.status.busy":"2023-10-29T10:09:53.861058Z","iopub.execute_input":"2023-10-29T10:09:53.861349Z","iopub.status.idle":"2023-10-29T10:09:54.172707Z","shell.execute_reply.started":"2023-10-29T10:09:53.861323Z","shell.execute_reply":"2023-10-29T10:09:54.171629Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Optionally, to de-fragment the DataFrame:\nrna_sequence_data = rna_sequence_data.copy()","metadata":{"execution":{"iopub.status.busy":"2023-10-29T10:09:54.175165Z","iopub.execute_input":"2023-10-29T10:09:54.175499Z","iopub.status.idle":"2023-10-29T10:09:54.181402Z","shell.execute_reply.started":"2023-10-29T10:09:54.175469Z","shell.execute_reply":"2023-10-29T10:09:54.180429Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"<a id =  \"codon\"></a>\n## An introduction to Codons","metadata":{}},{"cell_type":"markdown","source":"<a id =  \"cdn_hist\"></a>\n### History\n\nThe concept of the codon is fundamental to our understanding of molecular biology and genetics. Codons are the \"words\" of the genetic code, and they instruct cells on which amino acids to use when building proteins. Here's a brief overview of the history and discovery of codons:\n\n1. **The Central Dogma of Molecular Biology**: Proposed by Francis Crick in 1958, the central dogma describes the flow of genetic information in cells from DNA to RNA to protein. Crick's idea was that DNA encodes information, RNA decodes or transcribes that information, and proteins are the end products.\n\n2. **Triplets and the Genetic Code**: By the early 1960s, researchers had deduced that the genetic code was composed of nucleotide triplets, which means every three nucleotides in DNA or RNA specify a particular amino acid. This was supported by the fact that there are 4^3 (or 64) possible combinations of the four nucleotides (A, T, G, and C in DNA; A, U, G, and C in RNA), which is more than enough to code for the 20 common amino acids.\n\n3. **Deciphering the Code**: In the 1960s, Marshall Nirenberg, Heinrich Matthaei, and their colleagues made significant strides in breaking the genetic code. By using synthetic RNA molecules in bacterial systems, they determined which amino acid was specified by each RNA triplet or codon. Their pioneering experiments began with poly-U RNA (a sequence of uracil nucleotides) directing the synthesis of a polypeptide made solely of phenylalanine. This confirmed that the RNA codon UUU specified the amino acid phenylalanine.\n\n4. **64 Codons for 20 Amino Acids**: Of the 64 possible codons, 61 specify amino acids and three are \"stop\" codons, signaling the end of protein synthesis. Because there are more codons than there are amino acids, the genetic code is described as \"degenerate\" or \"redundant.\" This means that most amino acids are specified by more than one codon. For example, both UCU and UCC codons code for the amino acid serine.\n\n5. **tRNA and Anticodons**: The link between codons and amino acids is provided by a set of molecules called transfer RNAs (tRNAs). Each tRNA molecule carries a specific amino acid and has a three-nucleotide sequence called an anticodon that can base pair with a complementary codon in mRNA. This mechanism ensures the correct insertion of amino acids into the growing polypeptide chain during protein synthesis.\n\n6. **Evolution and Universality of the Genetic Code**: The genetic code is remarkably consistent across all known living organisms, from bacteria to humans, suggesting its ancient origins and the advantage of its particular design. However, some minor variations do exist, especially in certain organelles and some organisms, but these are exceptions rather than the rule.\n\n7. **Applications and Implications**: Understanding the genetic code has been central to the biotechnology revolution. It has allowed scientists to produce recombinant proteins, develop gene therapies, and even artificially synthesize genes.\n\nIn summary, the discovery and understanding of codons have been pivotal in the field of molecular biology. They have provided insight into how genetic information is stored, read, and translated into functional proteins in living organisms.\n\n[Return to Contents](#contents)","metadata":{}},{"cell_type":"markdown","source":"<a id =  \"cdn_anal\"></a>\n### Analysis\n\nAnalyzing codons involves determining the sequence of nucleotides in DNA or RNA and then interpreting that sequence in terms of the genetic code. Several methodologies and steps are involved in the process of analyzing codons, depending on the context and purpose of the analysis:\n\n1. **DNA Sequencing**:\n    * The first step to analyzing codons is often determining the nucleotide sequence of a segment of DNA. Modern methods like next-generation sequencing (NGS) allow researchers to sequence large fragments of DNA rapidly.\n2. **Transcription and Translation**:\n    - Once you have the DNA sequence, it can be transcribed to its corresponding RNA sequence (replacing thymine with uracil).\n    - This RNA sequence can then be \"translated\" into amino acids using the genetic code.\n3. **Bioinformatics Tools**:\n    - A variety of software tools and databases exist to help researchers analyze codons. These tools can quickly translate nucleotide sequences into amino acid sequences, identify open reading frames (ORFs), predict protein domains, and more.\n    - Tools like BLAST (Basic Local Alignment Search Tool) can be used to compare a given nucleotide sequence with databases to find similar sequences and annotate gene functions.\n4. **Codon Usage Analysis**:\n    - Different organisms have preferences for certain synonymous codons (different codons that code for the same amino acid) over others. This is called codon bias.\n    - Analyzing codon usage can provide insights into the evolutionary pressures on a gene, the likely level of protein expression, and can be important when expressing a gene in a heterologous system (like producing a human protein in bacteria).\n5. **Mutation Analysis**:\n    - Codons can be analyzed to identify mutations, either single nucleotide polymorphisms (SNPs) or larger mutations.\n    - Depending on the location and type of mutation, it might lead to synonymous changes (no change in the amino acid), nonsynonymous changes (change in the amino acid), or can introduce a premature stop codon.\n6. **tRNA Adaptation Index (tAI)**:\n    - tAI is a measure used to predict the speed of translation based on the availability of tRNAs for each codon. By analyzing codon usage in the context of available tRNAs in an organism, it's possible to make predictions about protein translation efficiency.\n7. **Experimental Approaches**:\n    - In addition to computational analyses, experimental methods like ribosome profiling can give insights into translation efficiency and how often certain codons are being read by ribosomes.\n8. **Structural and Functional Analysis**:\n    - Once codons are translated into amino acid sequences, these sequences can be further analyzed for potential functional domains, post-translational modification sites, or modeled for 3D structure using various bioinformatics tools.  \n\nIn essence, the analysis of codons is an interdisciplinary endeavor, merging molecular biology techniques with computational and statistical analyses to interpret genetic information and understand its implications for protein function, gene regulation, evolution, and more.\n\n[Return to Contents](#contents)","metadata":{}},{"cell_type":"markdown","source":"<a id =  \"cdn_freq\"></a>\n### Frequency plots\n\nA codon frequency plot (often referred to in the context of \"codon usage\" or \"codon bias\") provides a graphical representation of how often each codon is used to encode its corresponding amino acid in a particular organism, tissue, or gene. The purpose and interpretation of such a plot are as follows:\n\n**Purpose**\n\n1. **Understanding Evolutionary Pressures**:\n   - Codon usage patterns can provide insights into evolutionary pressures on an organism. For instance, genes under strong selective pressures might display codon usage patterns optimized for efficient and accurate translation.  \n2. **Optimizing Gene Expression**:\n   - Understanding codon usage is crucial when expressing a gene from one organism in another (heterologous expression). For efficient protein expression, the gene might need to be \"codon-optimized\" for the host organism to match its codon usage preferences.  \n3. **Inference about tRNA Availability**:\n   - Codons that are frequently used in an organism likely have corresponding tRNAs that are also abundant. Conversely, rare codons might slow down translation due to less available tRNAs.  \n4. **Studying Gene and Genome Evolution**:\n   - Differences in codon usage patterns among genes can hint at different evolutionary origins, such as horizontal gene transfer events.\n5. **Identifying Highly Expressed Genes**:\n   - In some organisms, highly expressed genes tend to use a specific subset of codons more frequently. By analyzing codon frequencies, it's possible to make inferences about gene expression levels.\n\n**Interpretation**\n\n1. **Axis Values**:\n   - Typically, on the x-axis, you'll find the 61 sense codons (or a subset if looking at specific amino acids). The y-axis typically displays the frequency of each codon.  \n2. **Peak Codons**:\n   - Codons that appear more frequently (higher y-values) are preferred codons for that organism or dataset. They might correlate with higher concentrations of the corresponding tRNA in the cell.  \n3. **Valley Codons**:\n   - Codons that appear less frequently (lower y-values) are less preferred and might be used more slowly in translation because of the limited availability of their corresponding tRNAs.  \n4. **Comparison Across Organisms or Tissues**:\n   - By overlaying codon frequency plots of different organisms or tissues, one can identify differences in codon preferences. Such differences might be due to varying evolutionary pressures, different tRNA pools, or specific adaptive requirements.  \n5. **Multiple Codons for an Amino Acid**:\n   - Amino acids encoded by multiple codons might show a bias towards one or a few of those codons. Such biases can reflect evolutionary optimizations for translation efficiency and accuracy.\n6. **Consider the Context**:\n   - Codon usage can be influenced by various factors like GC content of the genome, mutational biases, and selection for translational efficiency. Interpretation should consider the broader genomic context.\n\nIn summary, a codon frequency plot is a valuable tool for understanding the molecular biology and evolutionary pressures acting on organisms. Proper interpretation requires an integration of knowledge from genetics, evolutionary biology, and molecular biology.","metadata":{}},{"cell_type":"code","source":"overall_frequency = Counter()\n\nfor index, row in rna_sequence_data.iterrows():\n    sequence = row['sequence']\n    codons = [sequence[i:i+3] for i in range(0, len(sequence), 3)]\n    frequency = Counter(codons)\n    \n    overall_frequency.update(frequency)\n\ncustom_palette = sns.color_palette(\"husl\", len(overall_frequency))\n\nplt.figure(figsize=(16, 6))\nsns.barplot(x=list(overall_frequency.keys()), y=list(overall_frequency.values()), palette=custom_palette)\nplt.xticks(rotation=90)\nplt.xlabel('Codon')\nplt.ylabel('Frequency')\nplt.title('Overall Codon Frequency')\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-10-29T10:09:54.182607Z","iopub.execute_input":"2023-10-29T10:09:54.182933Z","iopub.status.idle":"2023-10-29T10:09:57.647320Z","shell.execute_reply.started":"2023-10-29T10:09:54.182906Z","shell.execute_reply":"2023-10-29T10:09:57.646284Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"[Return to Contents](#contents)","metadata":{}},{"cell_type":"markdown","source":"<a id =  \"cdn_heat\"></a>\n### Heat map\n\nA codon heat map is a graphical representation that shows the frequency or bias of codon usage across different genes, samples, or organisms. Instead of using a line or bar (as in a codon frequency plot), it uses colors to represent the values. Here's how to interpret it and its purposes:\n\n**Purpose**\n\n1. **Comparison Across Multiple Genes or Organisms**: Heat maps allow for an easy visualization of codon usage differences across multiple genes or organisms simultaneously.\n2. **Highlighting Patterns**: By visualizing codon usage data in a color-coded matrix format, researchers can quickly identify patterns, trends, or outliers that might not be evident in tabular data.\n3. **Assessing Codon Optimization**: When expressing a gene from one organism in another, understanding codon preferences is essential. A heat map can provide a quick snapshot of which codons might need optimization for efficient expression.\n4. **Studying Evolutionary and Adaptive Changes**: By comparing codon usage across different organisms or conditions, researchers can make inferences about evolutionary changes and possible adaptive strategies.\n\n**Interpretation**\n\n1. **Color Gradient**: Heat maps use a color gradient to represent values. In the context of codon usage, darker or more intense colors might represent higher frequencies (or other metrics like relative adaptiveness), while lighter colors denote lower values. Always refer to the accompanying color legend or scale to interpret the values correctly.\n2. **Rows and Columns**: Typically, each row in the heat map represents a specific gene or organism, and each column represents a specific codon. The intersection of a row and column will be color-coded based on the frequency of that codon for the corresponding gene or organism.\n3. **Patterns Across Rows or Columns**:\n   - If a particular codon column has consistently dark colors across all rows, that codon might be universally preferred across the analyzed genes or organisms.\n   - Conversely, if a specific row displays a unique color pattern compared to others, the corresponding gene or organism might have a distinct codon usage profile.\n4. **Clustering**: Often, heat maps are accompanied by dendrograms or clustering patterns, which group similar entities together. In the context of codon heat maps:\n   - Rows (genes or organisms) with similar codon usage profiles will cluster together.\n   - Columns (codons) that have similar usage patterns across all genes or organisms will cluster together.\n5. **Annotations**: Sometimes, additional information might be provided alongside the heat map, such as GC content, gene expression levels, or known evolutionary relationships. These annotations can provide context and help in drawing more nuanced conclusions from the map.\nIn essence, a codon heat map allows researchers to visually compare codon usage profiles across multiple entities in a compact and intuitive format. The patterns and differences observed in the heat map can provide insights into evolutionary processes, translational efficiency, and other molecular biology phenomena.","metadata":{}},{"cell_type":"code","source":"sequences = rna_sequence_data['sequence']\n\ncodon_freqs = []\nfor sequence in sequences:\n    codons = [sequence[i:i+3] for i in range(0, len(sequence), 3)]\n    frequency = Counter(codons)\n    codon_freqs.append(frequency)\n\ndf = pd.DataFrame(codon_freqs).fillna(0)\n\nplt.figure(figsize=(12, 8))\nsns.heatmap(df, cmap=\"YlGnBu\")\nplt.title('Heatmap of Codon Frequency')\nplt.xlabel('Codon')\n\ninterval = 10000\nplt.yticks(range(0, len(sequences), interval), range(0, len(sequences), interval))\n\nplt.ylabel('Sequence Index')\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-10-29T10:09:57.648820Z","iopub.execute_input":"2023-10-29T10:09:57.649217Z","iopub.status.idle":"2023-10-29T10:10:01.621141Z","shell.execute_reply.started":"2023-10-29T10:09:57.649178Z","shell.execute_reply":"2023-10-29T10:10:01.620129Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The RNA codon \"ACA\" serves as an instruction within messenger RNA (mRNA) during the process of protein synthesis. Specifically, \"ACA\" is a codon that specifies the incorporation of the amino acid threonine into a growing polypeptide chain during translation. \n\nHere's a breakdown of its purpose:\n\n1. **Protein Synthesis**: The main role of mRNA is to act as a template for the synthesis of proteins in the ribosomes. Codons in the mRNA, like \"ACA\", provide the sequential information to build a protein by specifying which amino acids should be added and in which order.\n\n2. **tRNA Recognition**: During translation, a specific transfer RNA (tRNA) molecule that carries the amino acid threonine will have an anticodon that's complementary to the \"ACA\" codon. This tRNA will recognize and bind to the \"ACA\" codon in the mRNA, ensuring that threonine gets added to the growing protein chain at the correct position.\n\n3. **Ensuring Genetic Fidelity**: The correct reading and translation of codons like \"ACA\" ensure that proteins are synthesized accurately according to the genetic information in the DNA. Mistakes in reading or translating codons can lead to non-functional or even harmful proteins.\n\nIn summary, the RNA codon \"ACA\" is an essential part of the molecular machinery that ensures the correct translation of genetic information from DNA to proteins.","metadata":{}},{"cell_type":"markdown","source":"<a id =  \"nuc_heat\"></a>\n### Nucleotide level heatmap\nA nucleotide-level heatmap is a graphical representation that visualizes nucleotide sequence data using colors. By representing nucleotide frequencies or other metrics using a color-coded system, the heatmap provides insights into patterns and variations within DNA or RNA sequences. Here's the purpose and interpretation of such a heatmap:\n\n**Purpose**\n\n1. **Visualizing Sequence Variation**:\n   - Heatmaps can reveal regions of high or low variability within a set of sequences. For instance, in a set of aligned sequences, positions that are conserved (i.e., the same nucleotide across all sequences) might be colored differently from positions with a lot of variability.\n2. **Comparing Multiple Sequences**:\n   - If you're comparing multiple sequences, a nucleotide-level heatmap can help you quickly identify regions of similarity and difference.\n3. **Highlighting Motifs and Features**:\n    - Certain nucleotide patterns or motifs might be of interest. A heatmap can help visualize the prevalence and location of these motifs across different sequences.\n4. **Analyzing Sequence Evolution**:\n   - By comparing nucleotide-level heatmaps of sequences from different organisms or populations, you can draw inferences about their evolutionary relationships and mutational patterns.\n5. **Studying DNA Methylation or Other Modifications**:\n   - In some cases, the heatmap might represent not just nucleotide identity but other features, like the frequency of DNA methylation at each cytosine position.\n\n**Interpretation:**\n\n1. **Color Legend**: Always begin by understanding the color legend associated with the heatmap. It will tell you what each color represents, be it nucleotide identity, frequency, level of conservation, or other metrics.\n2. **Rows and Columns**:\n   - Typically, each row of the heatmap might represent a different sequence or a different sample.\n   - Each column might represent a specific position within the sequence.\n3. **Conserved vs. Variable Regions**:\n   - If certain columns (positions) have the same color across all rows (sequences), this suggests conservation at those positions.\n   - Conversely, columns with a mix of colors suggest variability or mutations at those nucleotide positions.\n4. **Patterns and Clusters**:\n   - Look for overarching patterns or clusters in the heatmap. Clusters of similar colors can indicate related sequences or shared evolutionary history.\n   - Patterns might also suggest the presence of repeated motifs, regulatory elements, or structurally important regions.\n5. **External Annotations**:\n   - Sometimes, heatmaps are complemented by additional annotations like gene locations, functional domains, or known regulatory elements. This additional information can help interpret patterns observed in the heatmap.\n\nIn summary, a nucleotide-level heatmap provides a visual summary of nucleotide sequence data, making it easier to discern patterns, variations, and features within and across sequences. They are especially valuable in genomics and bioinformatics for visualizing large-scale sequence data.","metadata":{}},{"cell_type":"code","source":"# nucleotides mapping\nmapping = {\n    'A': 1,\n    'C': 2,\n    'G': 3,\n    'U': 4,\n    '.': 0\n}\n\nSEQUENCE_LENGTH = 170\nnumerical_sequences = rna_sequence_data['sequence'].apply(lambda x: [mapping[n] for n in x[:SEQUENCE_LENGTH]])\nrna_map = np.array(numerical_sequences.tolist())\n\nCMAP = 'Set3'\n\nplt.figure(figsize=(16,4))\nim = plt.imshow(rna_map, aspect='auto', cmap=CMAP, vmin=0, vmax=4, interpolation='nearest')\n\nrna_values = {\n    1: 'A',\n    2: 'C',\n    3: 'G',\n    4: 'U',\n    0: 'Pad'\n}\n\ncolor_values = [1, 2, 3, 4, 0]\ncolors = [im.cmap(im.norm(color_value)) for color_value in color_values]\npatches = [mpatches.Patch(color=colors[i], label=rna_values[color_value]) for i, color_value in enumerate(color_values)]\nplt.legend(handles=patches, bbox_to_anchor=(0.65, -0.1), ncol=5, borderaxespad=0.)\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-10-29T10:10:01.622400Z","iopub.execute_input":"2023-10-29T10:10:01.622701Z","iopub.status.idle":"2023-10-29T10:10:02.986313Z","shell.execute_reply.started":"2023-10-29T10:10:01.622674Z","shell.execute_reply":"2023-10-29T10:10:02.985247Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"[Return to Contents](#contents)","metadata":{}},{"cell_type":"markdown","source":"<a id =  \"nuc_comp\"></a>\n### Composition of A,U,G,C Nucleotides","metadata":{}},{"cell_type":"code","source":"nucleotide_counts = Counter(sequence)\n\nplt.pie(nucleotide_counts.values(), labels=nucleotide_counts.keys(), autopct='%1.1f%%')\nplt.title('Nucleotide Composition Pie Chart')\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-10-29T10:10:02.987588Z","iopub.execute_input":"2023-10-29T10:10:02.987903Z","iopub.status.idle":"2023-10-29T10:10:03.142675Z","shell.execute_reply.started":"2023-10-29T10:10:02.987874Z","shell.execute_reply":"2023-10-29T10:10:03.141364Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**GC Bond Counts on Sliding window = 30**\n\nThis refers to an analysis technique used to assess the GC content across a nucleotide sequence.\n\n1. **GC Content**: GC content refers to the percentage of guanine (G) and cytosine (C) bases in a DNA sequence relative to the total number of bases. The GC content can influence DNA's physical and biological properties, like its melting temperature and potentially its gene expression level.\n2. **Sliding Window**: A sliding window analysis involves examining a sequence in chunks, moving a fixed-sized window (in this case, 30 nucleotides) along the sequence one nucleotide at a time (or sometimes more than one, depending on the step size). For each position of the window, you calculate the metric of interest — in this case, the count of GC bonds.\n3. **Interpretation**:\n   - By analyzing the GC bond counts using a sliding window of 30 nucleotides, you can visualize how the GC content varies along the length of a sequence. \n   - Areas with higher GC bond counts in their windows have a higher GC content and might form more stable DNA regions due to the three hydrogen bonds formed between G and C (compared to the two bonds between adenine (A) and thymine (T)).\n   - Conversely, areas with lower GC bond counts are AT-rich. AT-rich regions might be more prone to denaturing at lower temperatures and might be associated with specific biological functions or regulatory regions.\n4. **Applications**:\n   - **Genomic Signatures**: Different organisms might have varying GC content. Some bacteria, for example, have high GC content, while others are AT-rich. Sliding window analysis can help identify regions that deviate from the norm, suggesting potential horizontal gene transfer events.\n   - **Isochore Identification**: Eukaryotic genomes often contain regions, called isochores, that have relatively homogeneous GC content. Identifying these can provide insights into genome structure and evolution.\n   - **Gene Prediction and Annotation**: Some coding regions might have different GC content compared to non-coding regions. Thus, analyzing GC content can provide hints about potential coding sequences.\n   - **Regulatory Elements**: Some regulatory elements, like promoters, might be associated with specific GC content profiles.\n\nIn summary, analyzing \"GC Bond Counts on Sliding window = 30\" provides a way to assess and visualize the variability of GC content across a DNA sequence, offering insights into its structure, function, and evolutionary history.\n\n- GC bonds are stronger than AT and AU bonds.\n- Strong GC Bonds are stable in high-temperature regions. (Organisms having more GC Bonds in DNA or RNA are more likely to survive in high temperature regions.)\n- GC bond abundant regions are more challenging to replicate leading to errors in replication.","metadata":{}},{"cell_type":"code","source":"window_size = 30 \ngc_content = [sequence[i:i+window_size].count('G') + sequence[i:i+window_size].count('C') for i in range(len(sequence) - window_size + 1)]\n\nplt.figure(figsize=(16, 6))\nplt.plot(gc_content)\nplt.xlabel('Window Start Position')\nplt.ylabel('GC Content')\nplt.title('Sliding Window GC Content')\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-10-29T10:10:03.144265Z","iopub.execute_input":"2023-10-29T10:10:03.144746Z","iopub.status.idle":"2023-10-29T10:10:03.507421Z","shell.execute_reply.started":"2023-10-29T10:10:03.144708Z","shell.execute_reply":"2023-10-29T10:10:03.506336Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"[Return to Contents](#contents)","metadata":{}},{"cell_type":"markdown","source":"<a id =  \"nuc_bias\"></a>\n### Bias due to Codon Pair Interaction\nBias due to codon pair interaction refers to the non-random association between pairs of codons in a given genome, which goes beyond the bias expected from individual codon frequencies. This phenomenon suggests that some pairs of codons are preferred or avoided in a sequence, not just because of the inherent bias of individual codons, but due to their interaction when they are adjacent.\n\nThere are several potential reasons for the existence of codon pair bias:\n\n1. **Translational Efficiency**: It's hypothesized that certain pairs of codons can be translated more efficiently than others. This could be due to the spatial arrangement and availability of tRNAs in the ribosome. Some codon pairs might allow for faster tRNA accommodation and peptide bond formation, leading to more efficient translation.\n\n2. **mRNA Structure**: The local structure of mRNA, such as stem-loops or other secondary structures, can influence the pairing of adjacent codons. Some codon pairs might stabilize or destabilize specific mRNA structures, which can subsequently impact translation.\n\n3. **Avoidance of Frame-shifting**: Certain codon pairs might be avoided because they promote ribosomal frame-shifting, leading to erroneous protein synthesis. \n\n4. **Genome Evolution and Horizontal Gene Transfer**: Genes acquired by horizontal gene transfer might have a codon and codon pair usage different from the native genes of the organism. Over time, the foreign genes might undergo amelioration, adopting the codon pair preferences of the host genome.\n\n5. **Protein Folding**: The speed of translation can influence protein folding. If certain codon pairs lead to pauses or changes in the translation rate, they can influence how the nascent polypeptide exits the ribosome and folds into its functional form. This can, in turn, create a selection for or against specific codon pairs based on the protein's functional requirements.\n\nRecognizing and understanding codon pair bias has practical applications:\n\n1. **Synthetic Biology and Genetic Engineering**: When designing genes for heterologous expression, considering codon pair bias can optimize protein production and function. This is especially important in the production of therapeutic proteins and industrial enzymes.\n\n2. **Vaccine Development**: Some attenuated viruses are designed using unfavorable codon pairs, ensuring they replicate less efficiently, leading to a weaker, controlled infection that can serve as a vaccine.\n\n3. **Evolutionary Studies**: Analyzing codon pair biases can provide insights into evolutionary pressures on genomes and hint at the functional constraints and requirements of proteins.\n\nIn summary, while individual codon usage bias is a recognized phenomenon, codon pair interaction adds another layer of complexity, suggesting that not only individual codons but also the context in which they appear can have functional and evolutionary implications.","metadata":{}},{"cell_type":"code","source":"overall_pair_freqs = Counter()\n\nfor seq in rna_sequence_data['sequence']:\n    pairs = [(seq[i:i+3], seq[i+3:i+6]) for i in range(0, len(seq) - 3, 3)]\n    overall_pair_freqs.update(pairs)\n\nindex = list(set([pair[0] for pair in overall_pair_freqs.keys()]))\ncolumns = list(set([pair[1] for pair in overall_pair_freqs.keys()]))\npair_df = pd.DataFrame(index=index, columns=columns).fillna(0)\n\nfor pair, freq in overall_pair_freqs.items():\n    pair_df.at[pair[0], pair[1]] = freq\n    \nplt.figure(figsize=(10, 6))\nsns.heatmap(pair_df, cmap='viridis')\nplt.title('Overall Codon Pair Bias')\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-10-29T10:10:03.511346Z","iopub.execute_input":"2023-10-29T10:10:03.511696Z","iopub.status.idle":"2023-10-29T10:10:05.538581Z","shell.execute_reply.started":"2023-10-29T10:10:03.511665Z","shell.execute_reply":"2023-10-29T10:10:05.537626Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"[Return to Contents](#contents)","metadata":{}},{"cell_type":"markdown","source":"<a id =  \"nuc_3d\"></a>\n### 3D Scatter Plot of Codon Counts for AUG, GCA, and GCU Codons","metadata":{}},{"cell_type":"code","source":"codon_counts = defaultdict(list)\n\nfor seq in rna_sequence_data['sequence']:\n    counter = Counter([seq[i:i+3] for i in range(0, len(seq), 3)])\n    \n    for codon in counter:\n        codon_counts[codon].append(counter[codon])\n\nfor codon, counts in codon_counts.items():\n    while len(counts) < len(rna_sequence_data['sequence']):\n        counts.append(0)\n\ncodon_count = pd.DataFrame(dict(codon_counts))\n\nfig = px.scatter_3d(codon_count, x='AUG', y='GCA', z='GCU')\nfig.show()","metadata":{"execution":{"iopub.status.busy":"2023-10-29T10:10:05.539851Z","iopub.execute_input":"2023-10-29T10:10:05.540219Z","iopub.status.idle":"2023-10-29T10:10:13.188616Z","shell.execute_reply.started":"2023-10-29T10:10:05.540180Z","shell.execute_reply":"2023-10-29T10:10:13.187584Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"[Return to Contents](#contents)","metadata":{}},{"cell_type":"markdown","source":"<a id =  \"embed\"></a>\n### Embeddings (Word2Vec)\nWord2Vec is a popular neural network-based technique used to produce vector representations of words called embeddings. These embeddings capture semantic meanings based on the context in which words appear in the text. The primary idea behind Word2Vec is that words that occur in similar contexts tend to have similar meanings.\n\nHere's a more detailed breakdown:\n\n1. **Embeddings**:\nEmbeddings are dense vector representations of words in a continuous vector space. Instead of representing words as discrete symbols (e.g., \"apple\" or \"orange\"), embeddings represent them as high-dimensional vectors (e.g., [0.23, -0.56, ...]). This representation allows for capturing semantic relationships between words.\n2. **Two Main Architectures**:\nWord2Vec can be implemented using one of two primary neural network architectures:\n    - **Continuous Bag-of-Words (CBOW)**: Given a context (a set of surrounding words), CBOW predicts the target word. For example, for the sentence \"The cat sat on the ___\", CBOW would use \"The\", \"cat\", \"sat\", and \"on\" to predict the word \"mat\".\n    - **Skip-Gram**: It's the opposite of CBOW. Given a target word, it predicts the surrounding context words. Using the same example, Skip-Gram would use \"mat\" to predict \"The\", \"cat\", \"sat\", and \"on\".\n3. **Semantic Relationships**:\nOne of the remarkable properties of Word2Vec embeddings is that they can capture various semantic and syntactic relationships. A classic example is the vector arithmetic that demonstrates relationships like \\[ \\text{vector}(\\text{\"king\"}) - \\text{vector}(\\text{\"man\"}) + \\text{vector}(\\text{\"woman\"}) \\approx \\text{vector}(\\text{\"queen\"}) \\]\n4. **Training**:\nWord2Vec models are trained on large corpora to learn word embeddings that capture contextual similarities. The underlying neural network adjusts the embeddings so that semantically similar words have vectors that are close in the vector space.\n5. **Applications**:\nWord2Vec embeddings have various applications in natural language processing tasks, including:\n    - **Semantic Similarity**: Determining how similar two pieces of text are.\n    - **Text Classification**: Assigning predefined categories to text.\n    - **Named Entity Recognition**: Identifying entities like names, places, and dates in the text.\n    - **Word Analogies**: Finding words that have a similar relationship as a provided word pair.\n6. **Extensions and Variants**:\nSince the introduction of Word2Vec, various other word embedding models have been proposed, such as FastText (which considers subword information) and GloVe (which is based on factorizing the word co-occurrence matrix).\n\nIn summary, Word2Vec embeddings are dense vector representations of words that capture semantic and syntactic information based on word co-occurrences in large text corpora. These embeddings have revolutionized many areas of natural language processing by providing a rich, dense representation of text that captures contextual meanings.","metadata":{}},{"cell_type":"code","source":"def sequence_to_kmers(sequence, k=3):\n    \"\"\"Convert sequence to overlapping k-mers.\"\"\"\n    return [sequence[i:i+k] for i in range(len(sequence) - k + 1)]\n\nkmers = rna_sequence_data['sequence'].apply(sequence_to_kmers)\n\nmodel = Word2Vec(sentences=kmers, vector_size=100, window=5, min_count=1, workers=4)\n\ndef sequence_to_vector(sequence, model):\n    kmers = sequence_to_kmers(sequence)\n    vectors = [model.wv[kmer] for kmer in kmers if kmer in model.wv]\n    return sum(vectors) / len(vectors)","metadata":{"execution":{"iopub.status.busy":"2023-10-29T10:10:13.190497Z","iopub.execute_input":"2023-10-29T10:10:13.191364Z","iopub.status.idle":"2023-10-29T10:10:23.283490Z","shell.execute_reply.started":"2023-10-29T10:10:13.191325Z","shell.execute_reply":"2023-10-29T10:10:23.282633Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Explanation**\n\nThe code is using the Word2Vec model, originally designed for natural language processing, to convert RNA sequence data into a dense vector representation.\n\n1. **sequence_to_kmers function**:\n    - This function takes a sequence (presumably a string of RNA nucleotides) and a value `k` (default is 3).\n    - It returns overlapping k-mers (subsequences of length `k`) from the given sequence.\n    - Example: For the sequence \"AUCGAU\" and k=3, the resulting k-mers would be ['AUC', 'UCG', 'CGA', 'GAU'].\n2. **kmers extraction**:\n    - `rna_sequence_data['sequence'].apply(sequence_to_kmers)` implies there's a dataframe `rna_sequence_data` with a column 'sequence' containing RNA sequences.\n    - For each sequence in the dataframe, the code extracts overlapping k-mers and stores them in the `kmers` variable.\n3. **Word2Vec model training**:\n    - The code then trains a Word2Vec model using the k-mers as sentences.\n    - `vector_size=100` specifies that each k-mer will be represented as a 100-dimensional vector.\n    - `window=5` implies that, during training, the model considers 5 neighboring k-mers as context for each k-mer.\n    - `min_count=1` ensures that all k-mers, even those that appear only once, are considered during training.\n    - `workers=4` means that the training will be parallelized across 4 CPU cores for faster computation.\n4. **sequence_to_vector function**:\n    - This function converts a given RNA sequence into a vector representation.\n    - For each k-mer in the sequence, if the k-mer exists in the trained Word2Vec model, its vector representation is fetched.\n    - The function then returns the average of all k-mer vectors for that sequence, resulting in a single 100-dimensional vector for the entire sequence.\n\nThe approach aims to capture the context in which different k-mers (subsequences) appear in the RNA sequences. By converting sequences into dense vectors, one can potentially utilize these vectors for downstream machine learning tasks, such as clustering, classification, or similarity searches. Using Word2Vec for sequence data is a creative way of borrowing techniques from NLP to understand and represent biological sequences.","metadata":{}},{"cell_type":"code","source":"rna_sequence_data['vectors'] = rna_sequence_data['sequence'].apply(lambda x: sequence_to_vector(x, model))","metadata":{"execution":{"iopub.status.busy":"2023-10-29T10:10:23.284653Z","iopub.execute_input":"2023-10-29T10:10:23.284958Z","iopub.status.idle":"2023-10-29T10:10:39.925793Z","shell.execute_reply.started":"2023-10-29T10:10:23.284931Z","shell.execute_reply":"2023-10-29T10:10:39.924751Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Explanation**\n\nAdds a new column to the `rna_sequence_data` dataframe named 'vectors'. This new column will contain the vector representations of the RNA sequences found in the 'sequence' column of the dataframe. \n\n1. **rna_sequence_data['sequence'].apply(...)**: Applies a function to each entry in the 'sequence' column of the `rna_sequence_data` dataframe.\n\n2. **lambda x: sequence_to_vector(x, model)**: This is an anonymous function (lambda function) that takes an RNA sequence `x` and returns its vector representation using the `sequence_to_vector` function and the trained `model`.\n\n3. **rna_sequence_data['vectors'] = ...**: The resulting vector representations are stored in a new 'vectors' column in the `rna_sequence_data` dataframe.\n\nAfter running this line of code, for each RNA sequence in the 'sequence' column, you will have a corresponding vector representation in the 'vectors' column. These vectors can be used for various downstream tasks, such as clustering similar sequences, visualizing the sequence space, or serving as input features for machine learning models.\n","metadata":{}},{"cell_type":"code","source":"rna_sequence_data = rna_sequence_data.dropna(subset=['reactivity'])\n\ncols_to_drop = ['sequence', 'experiment_type', 'signal_to_noise']\nexperiment_2A3_MaP = rna_sequence_data[rna_sequence_data['experiment_type'] == '2A3_MaP'].drop(columns=cols_to_drop).reset_index(drop=True)\nexperiment_DMS_MaP = rna_sequence_data[rna_sequence_data['experiment_type'] == 'DMS_MaP'].drop(columns=cols_to_drop).reset_index(drop=True)","metadata":{"execution":{"iopub.status.busy":"2023-10-29T10:10:39.927199Z","iopub.execute_input":"2023-10-29T10:10:39.927630Z","iopub.status.idle":"2023-10-29T10:10:39.954077Z","shell.execute_reply.started":"2023-10-29T10:10:39.927591Z","shell.execute_reply":"2023-10-29T10:10:39.953012Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Explanation**\n\n1. **rna_sequence_data.dropna(subset=['reactivity'])**:\n    - This line drops (removes) any rows from the dataframe where the 'reactivity' column has missing (NaN) values.\n    - The `subset` parameter specifies that only NaN values in the 'reactivity' column should be considered.\n    - The dataframe is then updated to reflect this change.\n2. **cols_to_drop**:\n    - This is a list containing column names that you want to drop later on. The columns specified are 'sequence', 'experiment_type', and 'signal_to_noise'.\n3. **experiment_2A3_MaP**:\n    - This line first filters the dataframe to only include rows where the 'experiment_type' column has the value '2A3_MaP'.\n    - Then, it drops the columns specified in `cols_to_drop`.\n    - The `.reset_index(drop=True)` ensures that the row indices are reset (0, 1, 2, ...), and the old indices are dropped rather than being added as a new column.\n    - The resulting dataframe is stored in `experiment_2A3_MaP`.\n4. **experiment_DMS_MaP**:\n    - This line performs the same operations as the previous one, but filters for rows where the 'experiment_type' column has the value 'DMS_MaP'.\n    - The resulting dataframe is stored in `experiment_DMS_MaP`.\nIn essence, you're cleaning and splitting the `rna_sequence_data` dataframe based on experiment types, resulting in two separate dataframes: `experiment_2A3_MaP` and `experiment_DMS_MaP`. Both of these dataframes will not have the 'sequence', 'experiment_type', and 'signal_to_noise' columns, and will only contain rows corresponding to their respective experiment types.","metadata":{}},{"cell_type":"code","source":"experiment_DMS_MaP.head()","metadata":{"execution":{"iopub.status.busy":"2023-10-29T10:10:39.955335Z","iopub.execute_input":"2023-10-29T10:10:39.955650Z","iopub.status.idle":"2023-10-29T10:10:39.973123Z","shell.execute_reply.started":"2023-10-29T10:10:39.955622Z","shell.execute_reply":"2023-10-29T10:10:39.972067Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"[Return to Contents](#contents)","metadata":{}},{"cell_type":"markdown","source":"<a id =  \"react\"></a>\n### An approach using average Reactivity","metadata":{}},{"cell_type":"code","source":"vector_data = pd.DataFrame(experiment_2A3_MaP['vectors'].to_list(), columns=[f'feature_{i}' for i in range(100)])\nexp_2a3 = pd.concat([experiment_2A3_MaP, vector_data], axis=1).drop(columns='vectors')\n\ncont_names = [f'feature_{i}' for i in range(100)]\ndep_var = 'reactivity'\n\nsplits = RandomSplitter(valid_pct=0.3)(range_of(exp_2a3))  # 70% train, 30% validation\nto = TabularPandas(exp_2a3, procs=[Normalize], cont_names=cont_names, y_names=dep_var, splits=splits)\n\ndls = to.dataloaders(bs=64)\nlearn = tabular_learner(dls, layers=[300,200], metrics=mae)\n\nlearn.lr_find(start_lr=1e-7, end_lr=10, num_it=100)\n\nlearn.fit_one_cycle(1, 1e-2)\n\nlearn.recorder.plot_lr_find()\n\nlearn.recorder.plot_loss()","metadata":{"execution":{"iopub.status.busy":"2023-10-29T10:10:39.974473Z","iopub.execute_input":"2023-10-29T10:10:39.974887Z","iopub.status.idle":"2023-10-29T10:10:53.099424Z","shell.execute_reply.started":"2023-10-29T10:10:39.974858Z","shell.execute_reply":"2023-10-29T10:10:53.098470Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Explanation**\n\nThis code is using the fastai library to prepare and execute a machine learning workflow to predict the 'reactivity' column in the `experiment_2A3_MaP` dataframe.\n\n1. **Vector data extraction**:\n    - The vectors in the 'vectors' column of `experiment_2A3_MaP` are converted into a DataFrame named `vector_data`, where each dimension of the vector becomes a separate column.\n    - Columns are named as `feature_0`, `feature_1`, ..., `feature_99`.\n2. **Concatenate and clean-up**:\n    - The `vector_data` dataframe is concatenated to `experiment_2A3_MaP` to form a new dataframe, and the original 'vectors' column is dropped.\n3. **Defining continuous features and dependent variable**:\n    - Continuous feature names are stored in `cont_names` (all the vector feature columns).\n    - The target variable (what you're trying to predict) is set to 'reactivity'.\n4. **Creating Train-Validation split**:\n    - Data is split into training and validation sets, with 70% of the data for training and 30% for validation using `RandomSplitter`.\n5. **Pre-processing and DataLoader creation**:\n    - The `TabularPandas` class from fastai is used to preprocess the data. Continuous features are normalized (mean-centered and scaled by standard deviation).\n    - `dataloaders` method creates DataLoaders for training and validation datasets with batch size set to 64.\n6. **Model initialization**:\n    - A tabular learner with two hidden layers (300 and 200 neurons) is initialized. The model will use Mean Absolute Error (MAE) as its metric.\n7. **Learning rate finder**:\n    - `learn.lr_find()` is a utility from fastai that helps in finding an optimal learning rate by increasing the learning rate gradually for each batch and recording how the loss changes.    \n8. **Train the model**:\n    - The model is trained for one epoch using the `fit_one_cycle` method with a learning rate of 1e-2.\n9. **Plott Learning Rate vs Loss**:\n    - The `plot_lr_find()` method plots the learning rates against the losses recorded during the learning rate finder process, helping in visually selecting an optimal learning rate.\n10. **Plot Training Loss**:\n    - The `plot_loss()` method plots the training and validation losses during training. This helps in diagnosing if the model is underfitting, overfitting, or performing well.\n    \nThis workflow uses the vectors derived from RNA sequences to train a model to predict 'reactivity'. The utility of fastai is demonstrated through its ease in preprocessing, model creation, training, and evaluation, all with just a few lines of code.","metadata":{}},{"cell_type":"code","source":"del vector_data\ndel exp_2a3\ndel experiment_DMS_MaP\ndel experiment_2A3_MaP\ndel df\ngc.collect()","metadata":{"execution":{"iopub.status.busy":"2023-10-29T10:10:53.100768Z","iopub.execute_input":"2023-10-29T10:10:53.101135Z","iopub.status.idle":"2023-10-29T10:10:53.554872Z","shell.execute_reply.started":"2023-10-29T10:10:53.101101Z","shell.execute_reply":"2023-10-29T10:10:53.553892Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Explanation**\n\n1. **`del vector_data`, `del exp_2a3`, etc.**:\n    - These lines use the `del` statement to delete the named variables or objects (`vector_data`, `exp_2a3`, `experiment_DMS_MaP`, `experiment_2A3_MaP`, and `df`) from the Python environment. This means that these objects are no longer accessible, and the memory they were using is now available to be reclaimed by the garbage collector.    \n2. **`gc.collect()`**:\n    - The `gc` module provides an interface to the garbage collection facility for reference cycles (cycles of references where objects refer to each other, preventing them from being garbage collected under normal conditions). The `gc.collect()` method forces a garbage collection, ensuring that all unreferenced memory is freed up.\n    - By explicitly calling `gc.collect()`, you're making sure that any memory that can be reclaimed is indeed reclaimed. This can be especially useful in environments like Jupyter Notebooks where memory usage can grow over time as different computations are performed.\n\nThe motivation for such code is often to free up memory, which can be crucial when working with large datasets or performing memory-intensive operations. If you're running into memory issues or just want to ensure you're being efficient with memory usage, periodically deleting unused objects and forcing garbage collection can be a good practice.\n\n[Return to Contents](#contents)","metadata":{}},{"cell_type":"markdown","source":"<a id =  \"rnn\"></a>\n### An approach using Seq-to-Seq RNN\n\nA Seq-to-Seq (Sequence-to-Sequence) RNN (Recurrent Neural Network) is a type of deep learning model designed to convert sequences from one domain (e.g., sentences in English) into sequences in another domain (e.g., sentences in French). This model architecture has been particularly popular for tasks like machine translation, speech recognition, and text summarization.\n\nHere's a basic outline of how a Seq-to-Seq RNN works:\n\n1. **Encoder**:\n    - The encoder processes the input sequence and compresses the information into a context (or hidden state).\n    - This is typically achieved using layers of RNNs (like LSTM or GRU cells) which process the sequence step-by-step, maintaining an evolving state that reflects the information seen so far.\n\n2. **Context Vector**:\n    - At the end of the encoder's processing, the final state of the encoder represents the context for the entire input sequence. This context vector aims to capture all the necessary information from the input sequence to produce the correct output sequence.\n\n3. **Decoder**:\n    - The decoder takes the context vector and produces the output sequence step-by-step.\n    - Like the encoder, it also typically uses layers of RNNs. It starts generating the output sequence from the context vector, producing one element of the output sequence at a time.\n\n4. **Teacher Forcing**:\n    - During training, rather than using the decoder's previous predictions as input for subsequent steps, the true previous output in the target sequence is fed into the decoder. This approach, called teacher forcing, can make training faster and more stable.\n\n5. **Attention Mechanism** (optional but often used):\n    - The basic Seq-to-Seq model assumes that the encoder produces a fixed-size context vector that captures all the information of the input sequence, which can be a limitation, especially for long sequences.\n    - An attention mechanism allows the decoder to focus on different parts of the input sequence at each step of output generation, providing a dynamic context instead of a fixed one. This has been shown to significantly improve performance for many sequence-to-sequence tasks.\n\n**Example - Machine Translation**:\nImagine you're trying to translate the English sentence \"How are you?\" to French, which is \"Comment ça va?\".\n\n1. The encoder processes the English words one by one and creates a context vector.\n2. The decoder starts with this context vector and begins producing the French translation, word by word.\n3. With attention, at each step of producing the French words, the decoder can 'focus' on different words in the English sentence. For instance, when translating \"ça\", the attention mechanism might focus on \"are\" in the English sentence.\n\nSeq-to-Seq models, especially with attention mechanisms, have been the backbone of many state-of-the-art systems in various NLP tasks. However, with the advent of Transformer-based architectures like BERT and GPT, the attention mechanism's concept has been further generalized, leading to even more powerful models.","metadata":{}},{"cell_type":"markdown","source":"### Data pre-processing","metadata":{}},{"cell_type":"code","source":"rna_sequence_data = pd.read_csv('/kaggle/input/stanford-ribonanza-rna-folding/train_data.csv', nrows = 150000)\nrna_sequence_data.head()","metadata":{"execution":{"iopub.status.busy":"2023-10-29T10:10:53.556119Z","iopub.execute_input":"2023-10-29T10:10:53.556429Z","iopub.status.idle":"2023-10-29T10:11:01.519891Z","shell.execute_reply.started":"2023-10-29T10:10:53.556403Z","shell.execute_reply":"2023-10-29T10:11:01.518816Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"rna_sequence_data = rna_sequence_data.drop(columns=rna_sequence_data.filter(like='error').columns)\nrna_sequence_data = rna_sequence_data.drop(['sequence_id', 'dataset_name', 'reads', 'SN_filter'], axis=1)\nreactivity_cols = rna_sequence_data.filter(like='reactivity').columns\nrna_sequence_data[reactivity_cols] = rna_sequence_data[reactivity_cols].clip(lower=0)\nrow_means = rna_sequence_data[reactivity_cols].mean(axis=1)\nrna_sequence_data[reactivity_cols] = rna_sequence_data[reactivity_cols].apply(lambda col: col.fillna(row_means))\nrna_sequence_data[reactivity_cols] = rna_sequence_data[reactivity_cols].fillna(0)","metadata":{"execution":{"iopub.status.busy":"2023-10-29T10:11:01.521204Z","iopub.execute_input":"2023-10-29T10:11:01.521552Z","iopub.status.idle":"2023-10-29T10:11:03.950671Z","shell.execute_reply.started":"2023-10-29T10:11:01.521512Z","shell.execute_reply":"2023-10-29T10:11:03.949609Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Explanation**\n\n1. **Remove Error Columns**:\n    ```python\n    rna_sequence_data = rna_sequence_data.drop(columns=rna_sequence_data.filter(like='error').columns)\n    ```\n    This line drops all columns with names containing the word \"error\". The `filter(like='error')` method returns columns with names that contain the substring 'error'.\n\n2. **Drop Specific Columns**:\n    ```python\n    rna_sequence_data = rna_sequence_data.drop(['sequence_id', 'dataset_name', 'reads', 'SN_filter'], axis=1)\n    ```\n    Here, specific columns ('sequence_id', 'dataset_name', 'reads', and 'SN_filter') are dropped from the DataFrame.\n\n3. **Identify Reactivity Columns**:\n    ```python\n    reactivity_cols = rna_sequence_data.filter(like='reactivity').columns\n    ```\n    This line identifies and stores all columns with names containing the word \"reactivity\".\n\n4. **Clip Reactivity Values**:\n    ```python\n    rna_sequence_data[reactivity_cols] = rna_sequence_data[reactivity_cols].clip(lower=0)\n    ```\n    For all the columns identified as \"reactivity\" columns, any values below 0 are clipped to 0. This ensures that there are no negative values in these columns.\n\n5. **Replace NaN with Row-wise Mean for Reactivity Columns**:\n    ```python\n    row_means = rna_sequence_data[reactivity_cols].mean(axis=1)\n    rna_sequence_data[reactivity_cols] = rna_sequence_data[reactivity_cols].apply(lambda col: col.fillna(row_means))\n    ```\n    If any cells in the \"reactivity\" columns have missing values (`NaN`), those missing values are replaced with the mean value of the respective row for those columns.\n\n6. **Fill Remaining NaN Values with 0**:\n    ```python\n    rna_sequence_data[reactivity_cols] = rna_sequence_data[reactivity_cols].fillna(0)\n    ```\n    Any remaining missing values (`NaN`) in the \"reactivity\" columns are replaced with 0.\n\nIn summary, the code is cleaning and transforming the `rna_sequence_data` DataFrame by removing unwanted columns, clipping values to ensure non-negativity, and dealing with missing values in a specific manner for columns related to \"reactivity\". These preprocessing steps are often crucial before performing further analyses or modeling tasks.\n","metadata":{}},{"cell_type":"code","source":"def encode_nucleotide(nucleotide):\n    mapping = {'A': [1, 0, 0, 0],\n               'U': [0, 1, 0, 0],\n               'G': [0, 0, 1, 0],\n               'C': [0, 0, 0, 1]}\n    return mapping[nucleotide]\n\ndef encode_sequence(sequence):\n    return [encode_nucleotide(n) for n in sequence]\n\nrna_sequence_data['encoded_sequence'] = rna_sequence_data['sequence'].apply(encode_sequence)\n\nrna_sequence_data['padded_sequence'] = rna_sequence_data['encoded_sequence'].apply(lambda x: pad_sequences([x], maxlen=457, padding='post')[0])\n\ndef pad_reactivity_sequence(sequence, max_length=457, padding_value=-0.1):\n    # Calculate how many padding values are needed\n    padding_length = max_length - len(sequence)\n    \n    # Create the padded sequence\n    padded_sequence = list(sequence) + [padding_value] * padding_length\n    \n    return padded_sequence","metadata":{"execution":{"iopub.status.busy":"2023-10-29T10:11:03.952272Z","iopub.execute_input":"2023-10-29T10:11:03.952708Z","iopub.status.idle":"2023-10-29T10:12:26.189204Z","shell.execute_reply.started":"2023-10-29T10:11:03.952669Z","shell.execute_reply":"2023-10-29T10:12:26.188116Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Optionally, to de-fragment the DataFrame:\nrna_sequence_data = rna_sequence_data.copy()","metadata":{"execution":{"iopub.status.busy":"2023-10-29T10:12:26.190485Z","iopub.execute_input":"2023-10-29T10:12:26.190765Z","iopub.status.idle":"2023-10-29T10:12:26.589644Z","shell.execute_reply.started":"2023-10-29T10:12:26.190740Z","shell.execute_reply":"2023-10-29T10:12:26.588565Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Explanation**\n\nThis code provides a series of functions and operations to pre-process RNA sequence data and prepare it for modeling or further analyses.\n\n1. **Encoding Nucleotides**:\n    - The `encode_nucleotide` function is a straightforward one-hot encoding for nucleotides. Each nucleotide (`A`, `U`, `G`, `C`) is represented as a list of length 4, where only one position is `1` and the rest are `0`s. This is a common method in bioinformatics for converting nucleotide sequences to a numerical format suitable for machine learning models.\n    \n2. **Encoding Entire Sequences**:\n    - The `encode_sequence` function takes an entire RNA sequence (a string) and applies the `encode_nucleotide` function to each nucleotide in the sequence. The result is a list of lists, where each inner list is a one-hot encoded representation of a nucleotide.\n    - The `rna_sequence_data['encoded_sequence']` line applies this encoding to the entire 'sequence' column of the DataFrame.\n\n3. **Padding Sequences**:\n    - Sequence data often comes in variable lengths, but many machine learning models (especially deep learning models) require fixed-length input. To handle this, sequences are often padded to a consistent length.\n    - The `rna_sequence_data['padded_sequence']` line pads each encoded sequence to a maximum length of 457. If a sequence is shorter than this, it is extended with zeros at the end (i.e., 'post' padding). The `pad_sequences` function (presumably from TensorFlow's Keras API) is used to perform this padding.\n    \n4. **Padding Reactivity Sequences**:\n    - The `pad_reactivity_sequence` function is designed to pad sequences (specifically sequences of reactivity values, it seems) to a consistent length. If a sequence is shorter than the specified maximum length (`max_length`), it is extended by adding a specific padding value (`padding_value`, which defaults to `-0.1`) to the end. This function differs from the previous padding function in that it's designed for sequences of reactivity values, not encoded nucleotides.\n\nThis series of operations prepares the RNA sequence data in a format that's often suitable for feeding into machine learning models, especially deep learning models like RNNs, CNNs, or Transformers.","metadata":{}},{"cell_type":"code","source":"# 1. Consolidate Reactivity Columns\nreactivity_columns = [f'reactivity_{i:04}' for i in range(1, 207)]  # Adjust this based on the number of columns you have\n\nrna_sequence_data['reactivity'] = rna_sequence_data[reactivity_columns].values.tolist()\n\n# Drop the individual reactivity columns as they're now redundant (optional)\nrna_sequence_data.drop(columns=reactivity_columns, inplace=True)\n\n# Apply the padding function\nrna_sequence_data['padded_reactivity'] = rna_sequence_data['reactivity'].apply(pad_reactivity_sequence)\n\n# 1. Store original sequence lengths\nrna_sequence_data['original_length'] = rna_sequence_data['reactivity'].apply(len)\n\n# 2. Inverse transform function\ndef inverse_transform(padded_sequence, original_length):\n    return padded_sequence[:original_length]\n\n# Apply the inverse transform\nrna_sequence_data['retrieved_reactivity'] = rna_sequence_data.apply(lambda row: inverse_transform(row['padded_reactivity'], row['original_length']), axis=1)","metadata":{"execution":{"iopub.status.busy":"2023-10-29T10:12:26.591425Z","iopub.execute_input":"2023-10-29T10:12:26.591750Z","iopub.status.idle":"2023-10-29T10:12:33.574128Z","shell.execute_reply.started":"2023-10-29T10:12:26.591721Z","shell.execute_reply":"2023-10-29T10:12:33.573341Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Explanation**\n\nThe provided code executes a series of transformations on the `rna_sequence_data` DataFrame to handle RNA sequence reactivity values. The operations can be understood as follows:\n\n1. **Consolidate Reactivity Columns**:\n   - This portion of the code attempts to consolidate multiple columns of reactivity data into a single column where each entry is a list of reactivity values.\n   - The line `reactivity_columns = [f'reactivity_{i:04}' for i in range(1, 207)]` generates a list of expected column names that store individual reactivity values. It assumes that these columns are named 'reactivity_0001', 'reactivity_0002', ..., 'reactivity_0206'.\n   - `rna_sequence_data['reactivity']` is created as a new column in the DataFrame. Each entry of this column is a list that consolidates the reactivity values from the individual columns.\n   - The individual reactivity columns are then dropped to reduce redundancy and save memory using the `drop` method.\n\n2. **Padding the Reactivity Sequences**:\n   - The code applies the `pad_reactivity_sequence` function (from the previous code block) to the consolidated reactivity sequences. The resultant padded sequences are stored in the 'padded_reactivity' column.\n\n3. **Store Original Sequence Lengths**:\n   - Before transforming (or padding) sequences, it's often beneficial to keep track of their original lengths. This allows for the potential to reverse any transformations applied to the sequences.\n   - The `rna_sequence_data['original_length']` line creates a new column that stores the length of each reactivity sequence before padding.\n\n4. **Inverse Transformation**:\n   - After processing sequences (like padding), there may come a time when you wish to retrieve the original sequence. An inverse transformation function serves this purpose.\n   - The `inverse_transform` function trims a padded sequence back to its original length. It takes the padded sequence and its original length as arguments.\n   - The `rna_sequence_data['retrieved_reactivity']` line applies this inverse transformation to all rows in the DataFrame, effectively retrieving the original sequences before they were padded.\n\nThe operations in this code are useful when preparing sequence data for machine learning while also ensuring the possibility to revert to the original data when required.","metadata":{}},{"cell_type":"code","source":"rna_2A3 = rna_sequence_data[rna_sequence_data['experiment_type'] == '2A3_MaP']\nrna_dms = rna_sequence_data[rna_sequence_data['experiment_type'] == 'DMS_MaP']\n\ndel rna_sequence_data\ngc.collect()","metadata":{"execution":{"iopub.status.busy":"2023-10-29T10:12:33.575274Z","iopub.execute_input":"2023-10-29T10:12:33.575593Z","iopub.status.idle":"2023-10-29T10:12:37.234885Z","shell.execute_reply.started":"2023-10-29T10:12:33.575566Z","shell.execute_reply":"2023-10-29T10:12:37.233706Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"device = torch.device(\"cuda:0\" if torch.cuda.is_available() else \"cpu\")\n\nclass RNASeq2SeqModel(nn.Module):\n    def __init__(self, embedding_dim, hidden_dim, num_layers):\n        super(RNASeq2SeqModel, self).__init__()\n        \n        # Input is 4 due to one-hot encoding of A, U, G, C\n        self.embedding = nn.Embedding(4, embedding_dim)  # Embedding layer\n        self.rnn = nn.LSTM(embedding_dim, hidden_dim, num_layers, batch_first=True)\n        self.fc = nn.Linear(hidden_dim, 1)  # Predicting reactivity for each nucleotide\n        \n    def forward(self, x):\n        x = torch.argmax(x, dim=-1)  # Convert one-hot to indices for embedding\n        x = self.embedding(x)\n        rnn_out, _ = self.rnn(x)\n        output = self.fc(rnn_out)\n        return output\n\n    \nsequence_matrix = rna_2A3['padded_sequence']\nsequence_tensor = torch.FloatTensor(sequence_matrix)\nreactivities = torch.FloatTensor(rna_2A3['padded_reactivity'])","metadata":{"execution":{"iopub.status.busy":"2023-10-29T10:12:37.236283Z","iopub.execute_input":"2023-10-29T10:12:37.236659Z","iopub.status.idle":"2023-10-29T10:14:09.625389Z","shell.execute_reply.started":"2023-10-29T10:12:37.236627Z","shell.execute_reply":"2023-10-29T10:14:09.624263Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"This warning indicates that creating a PyTorch tensor directly from a list of numpy arrays can be inefficient. This inefficiency stems from the need to individually convert each numpy array in the list to a tensor and then combine them. The recommended approach is to first convert the list of numpy arrays into a single numpy array, and then convert this combined array into a tensor.\n\nHere's a step-by-step example to illustrate the recommendation:\n\n1. **List of numpy arrays**:\n   Let's assume you have the following list of numpy arrays:\n   ```python\n   import numpy as np\n\n   list_of_arrays = [np.array([1, 2]), np.array([3, 4]), np.array([5, 6])]\n   ```\n\n2. **Convert list to a single numpy array**:\n   Use `numpy.array()` to combine the list into one array:\n   ```python\n   combined_array = np.array(list_of_arrays)\n   ```\n\n3. **Convert numpy array to tensor**:\n   Now, use PyTorch to convert the combined numpy array to a tensor:\n   ```python\n   import torch\n\n   tensor = torch.tensor(combined_array)\n   ```\n\nBy following these steps, you can efficiently convert a list of numpy arrays to a PyTorch tensor. This method reduces overhead and speeds up the conversion process.","metadata":{}},{"cell_type":"markdown","source":"**Explanation**\n\nThe provided code is the beginning of an implementation of a sequence-to-sequence (Seq2Seq) model using PyTorch. The model's primary objective is to predict RNA reactivity values from RNA sequences. Here's a breakdown of the code:\n\n1. **Device Selection**:\n   - The line `device = torch.device(\"cuda:0\" if torch.cuda.is_available() else \"cpu\")` checks whether a GPU with CUDA support is available. If it is, PyTorch computations will be directed to use the GPU, otherwise, the CPU will be used. This allows for faster computations when a GPU is available.\n2. **RNASeq2SeqModel Class**:\n   - This class defines the RNA Seq2Seq model. Here are its key components:   \n     - **Embedding Layer**: This layer transforms the one-hot encoded nucleotides (A, U, G, C) into dense vectors of size `embedding_dim`. The `nn.Embedding` layer in PyTorch requires that the input be indices (not one-hot vectors), which is why there's a `torch.argmax` in the forward method to convert from one-hot vectors to indices.\n     - **RNN Layer**: The LSTM (a type of RNN) is used to process the sequence. The LSTM will read the sequence and maintain hidden states that summarize the information seen so far.\n     - **Fully Connected Layer**: This layer maps from the LSTM's hidden state to a predicted reactivity value for each nucleotide.   \n3. **Data Preparation**:\n   - The `sequence_matrix` extracts the 'padded_sequence' column from the dataset (which seems to be named `rna_2A3` here). This matrix contains the one-hot encoded RNA sequences.\n   - The `sequence_tensor` and `reactivities` lines convert the sequence matrix and padded reactivity values to PyTorch tensors. These tensors can be fed into the model during training.\n\nTo complete the implementation, you would still need to:\n- Define a loss function (e.g., Mean Squared Error for regression tasks like this).\n- Choose an optimizer (e.g., Adam).\n- Implement the training loop where you feed sequences into the model, compute the loss by comparing the model's predictions to the true reactivity values, and update the model's weights.\n\nAdditionally, ensuring that the data and model are moved to the selected device (`device`) is also important for efficient computations, especially when using a GPU.","metadata":{}},{"cell_type":"markdown","source":"### Training RNN Seq-to-Seq on 2A3_MaP data","metadata":{}},{"cell_type":"code","source":"# Create Dataset\ndataset = TensorDataset(sequence_tensor, reactivities)\n\n# Split data into train and validation sets (80% train, 20% validation for this example)\ntrain_size = int(0.8 * len(dataset))\nval_size = len(dataset) - train_size\ntrain_dataset, val_dataset = random_split(dataset, [train_size, val_size])\n\ntrain_loader = DataLoader(train_dataset, batch_size=64, shuffle=True)\nval_loader = DataLoader(val_dataset, batch_size=64, shuffle=False)\n\n# Hyperparameters\nembedding_dim = 16\nhidden_dim = 32\nnum_layers = 2\nlearning_rate = 0.001\nnum_epochs = 100\n\n# Instantiate model\nmodel_2a3 = RNASeq2SeqModel(embedding_dim=16, hidden_dim=32, num_layers=2).to(device)\ncriterion = nn.MSELoss()\noptimizer = torch.optim.Adam(model_2a3.parameters(), lr=0.001)\n\nfor epoch in range(num_epochs):\n    \n    model_2a3.train()\n    total_loss = 0.0\n    \n    for batch_seq, batch_react in train_loader:\n        \n        batch_seq, batch_react = batch_seq.to(device), batch_react.to(device) \n        \n        # Zero the parameter gradients\n        optimizer.zero_grad()\n        \n        # Forward pass\n        outputs = model_2a3(batch_seq)\n        \n        # Compute loss\n        loss = criterion(outputs.squeeze(), batch_react)\n        \n        # Backward pass and optimize\n        loss.backward()\n        optimizer.step()\n        \n        total_loss += loss.item()\n    \n    # Validate and compute MSE for the validation set\n    model_2a3.eval()\n    val_loss = 0.0\n    with torch.no_grad():\n        for batch_seq, batch_react in val_loader:\n            batch_seq, batch_react = batch_seq.to(device), batch_react.to(device)  # Transfer data to GPU\n            outputs = model_2a3(batch_seq)\n            loss = criterion(outputs.squeeze(), batch_react)\n            val_loss += loss.item()\n        \n    print(f\"Epoch {epoch+1}/{num_epochs}, Training Loss: {total_loss/len(train_loader)}, Validation MSE: {val_loss/len(val_loader)}\")\n\nprint(\"Training finished.\")","metadata":{"execution":{"iopub.status.busy":"2023-10-29T10:19:45.813357Z","iopub.execute_input":"2023-10-29T10:19:45.814529Z","iopub.status.idle":"2023-10-29T10:31:22.049672Z","shell.execute_reply.started":"2023-10-29T10:19:45.814489Z","shell.execute_reply":"2023-10-29T10:31:22.048670Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"del dataset\ndel train_dataset\ndel val_dataset\ndel train_loader\ndel val_loader\ndel sequence_tensor\ndel reactivities\n\ngc.collect()","metadata":{"execution":{"iopub.status.busy":"2023-10-29T10:37:34.995431Z","iopub.execute_input":"2023-10-29T10:37:34.995894Z","iopub.status.idle":"2023-10-29T10:37:38.584599Z","shell.execute_reply.started":"2023-10-29T10:37:34.995861Z","shell.execute_reply":"2023-10-29T10:37:38.583636Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Training RNN Seq-to-Seq on DMS_MaP data","metadata":{}},{"cell_type":"code","source":"sequence_matrix_dms = np.stack(rna_dms['padded_sequence'].to_numpy())  # Convert Series of lists to a numpy array\nsequence_tensor_dms = torch.FloatTensor(sequence_matrix_dms)\nreactivities_dms = torch.FloatTensor(np.stack(rna_dms['padded_reactivity'].to_numpy()))\n\n# Create Dataset for rna_dms\ndataset_dms = TensorDataset(sequence_tensor_dms, reactivities_dms)\n\n# Split data into train and validation sets (80% train, 20% validation for this example)\ntrain_size_dms = int(0.8 * len(dataset_dms))\nval_size_dms = len(dataset_dms) - train_size_dms\ntrain_dataset_dms, val_dataset_dms = random_split(dataset_dms, [train_size_dms, val_size_dms])\n\ntrain_loader_dms = DataLoader(train_dataset_dms, batch_size=64, shuffle=True)\nval_loader_dms = DataLoader(val_dataset_dms, batch_size=64, shuffle=False)","metadata":{"execution":{"iopub.status.busy":"2023-10-29T10:43:53.207481Z","iopub.execute_input":"2023-10-29T10:43:53.208453Z","iopub.status.idle":"2023-10-29T10:43:57.565551Z","shell.execute_reply.started":"2023-10-29T10:43:53.208420Z","shell.execute_reply":"2023-10-29T10:43:57.564722Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Explanation**\n\nThe code converts RNA sequence and reactivity data into PyTorch tensors, creates a dataset using those tensors and then splits the dataset into training and validation sets. It then loads these datasets into data loaders that will be used for training and validation.\n\n1. **Converting RNA Data to Numpy Array and then to PyTorch Tensors**:\n   ```python\n   sequence_matrix_dms = np.stack(rna_dms['padded_sequence'].to_numpy())  # Convert Series of lists to a numpy array\n   sequence_tensor_dms = torch.FloatTensor(sequence_matrix_dms)\n   reactivities_dms = torch.FloatTensor(np.stack(rna_dms['padded_reactivity'].to_numpy()))\n   ```\n   In this segment, RNA sequences are first stacked into a numpy array, and then this numpy array is converted into a PyTorch tensor of type `FloatTensor`. The same is done for the reactivities.\n\n2. **Creating a Dataset for RNA Data**:\n   ```python\n   dataset_dms = TensorDataset(sequence_tensor_dms, reactivities_dms)\n   ```\n   Here, the RNA sequence tensor and reactivity tensor are combined into a `TensorDataset`. This will allow easy batching and indexing of the dataset when training a machine learning model.\n\n3. **Splitting the Data into Training and Validation Sets**:\n   ```python\n   train_size_dms = int(0.8 * len(dataset_dms))\n   val_size_dms = len(dataset_dms) - train_size_dms\n   train_dataset_dms, val_dataset_dms = random_split(dataset_dms, [train_size_dms, val_size_dms])\n   ```\n   In this part, the data is split into training and validation sets. 80% of the data is used for training and 20% for validation. The `random_split` function ensures that the split is random.\n\n4. **Creating DataLoaders**:\n   ```python\n   train_loader_dms = DataLoader(train_dataset_dms, batch_size=64, shuffle=True)\n   val_loader_dms = DataLoader(val_dataset_dms, batch_size=64, shuffle=False)\n   ```\n   DataLoaders are created for both the training and validation sets. These are essential for efficiently iterating over the datasets in batches during training and validation. The training set is shuffled (`shuffle=True`) to ensure the model receives data in a different order each epoch, helping the model generalize better. The validation set doesn't need to be shuffled (`shuffle=False`), as the order of data doesn't affect validation metrics.\n\nThis is a standard workflow for preprocessing and setting up data for training deep learning models using PyTorch. The next steps (not shown in the code) would typically involve defining a model architecture, setting up a loss function and optimizer, and then training the model using the provided dataloaders.","metadata":{}},{"cell_type":"code","source":"# Instantiate another model for rna_dms\nmodel_dms = RNASeq2SeqModel(embedding_dim=embedding_dim, hidden_dim=hidden_dim, num_layers=num_layers).to(device)\noptimizer_dms = torch.optim.Adam(model_dms.parameters(), lr=learning_rate)\n\n# Training Loop for rna_dms\nfor epoch in range(num_epochs):\n    model_dms.train()\n    total_loss = 0.0\n    for batch_seq, batch_react in train_loader_dms:\n        batch_seq, batch_react = batch_seq.to(device), batch_react.to(device) \n        optimizer_dms.zero_grad()\n        \n        # Forward pass\n        outputs = model_dms(batch_seq)\n        \n        # Compute loss\n        loss = criterion(outputs.squeeze(), batch_react)\n        \n        # Backward pass and optimize\n        loss.backward()\n        optimizer_dms.step()\n        \n        total_loss += loss.item()\n    \n    # Validate and compute MSE for the validation set\n    model_dms.eval()\n    val_loss = 0.0\n    with torch.no_grad():\n        for batch_seq, batch_react in val_loader_dms:\n            batch_seq, batch_react = batch_seq.to(device), batch_react.to(device)  # Transfer data to GPU\n            outputs = model_dms(batch_seq)\n            loss = criterion(outputs.squeeze(), batch_react)\n            val_loss += loss.item()\n        \n    \n    print(f\"[DMS] Epoch {epoch+1}/{num_epochs}, Training Loss: {total_loss/len(train_loader_dms)}, Validation MSE: {val_loss/len(val_loader_dms)}\")\n\nprint(\"Training for DMS finished.\")","metadata":{"execution":{"iopub.status.busy":"2023-10-29T10:49:09.287852Z","iopub.execute_input":"2023-10-29T10:49:09.288614Z","iopub.status.idle":"2023-10-29T11:00:25.813369Z","shell.execute_reply.started":"2023-10-29T10:49:09.288573Z","shell.execute_reply":"2023-10-29T11:00:25.812326Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Explantion**\n\nThis is a simple training loop for a deep learning model in PyTorch.\n\n1. **Model and Optimizer Instantiation**:\n   ```python\n   model_dms = RNASeq2SeqModel(embedding_dim=embedding_dim, hidden_dim=hidden_dim, num_layers=num_layers).to(device)\n   optimizer_dms = torch.optim.Adam(model_dms.parameters(), lr=learning_rate)\n   ```\n   An instance of an `RNASeq2SeqModel` is created with specified parameters. This model is then moved to a computing device (`device`), which could be a CPU or GPU. The Adam optimizer is then initialized with the model's parameters and a given learning rate.\n\n2. **Training Loop**:\n   The model is trained for `num_epochs`. For each epoch:\n   - The model is set to training mode using `model_dms.train()`.\n   - A total loss variable is initialized to zero.\n   - Each batch of data is loaded from the training data loader (`train_loader_dms`).\n   - The optimizer's gradients are zeroed out using `optimizer_dms.zero_grad()`.\n   - The model makes predictions (forward pass) using the input sequences.\n   - The loss is computed using a loss function (`criterion`) which compares the model's predictions to the true reactivities.\n   - Gradients are calculated via backpropagation using `loss.backward()`.\n   - The optimizer updates the model's parameters using `optimizer_dms.step()`.\n   - The training loss for the batch is added to the total loss.\n\n3. **Validation Loop**:\n   After processing all batches for an epoch:\n   - The model is set to evaluation mode using `model_dms.eval()`.\n   - A total validation loss variable is initialized to zero.\n   - Each batch of data is loaded from the validation data loader (`val_loader_dms`).\n   - The model makes predictions for the validation set.\n   - The loss for the validation batch is computed and added to the total validation loss.\n   - Gradients are not calculated during the validation phase, so the code is wrapped in a `torch.no_grad()` context to save memory and speed up computations.\n\n4. **Printing Epoch Details**:\n   The average training and validation losses for the epoch are printed.\n\n5. A message is printed indicating the completion of training.\n\nFor this training loop to work, certain elements need to be predefined:\n- `RNASeq2SeqModel`: The model architecture.\n- `device`: A variable specifying the compute device, e.g., `torch.device(\"cuda\" if torch.cuda.is_available() else \"cpu\")`.\n- `criterion`: The loss function. Given the usage of \"MSE\" in the print statement, it might be the mean squared error loss.\n\nThis is a standard training loop structure for deep learning in PyTorch. The code seems correct, given that the necessary contextual elements (like the model definition, loss function, and device) are defined elsewhere.","metadata":{}},{"cell_type":"code","source":"torch.save(model_dms, \"/kaggle/working/models/model_dms.pth\")\ntorch.save(model_2a3, \"/kaggle/working/models/model_2a3.pth\")","metadata":{"execution":{"iopub.status.busy":"2023-10-29T11:12:01.942182Z","iopub.execute_input":"2023-10-29T11:12:01.943055Z","iopub.status.idle":"2023-10-29T11:12:01.951278Z","shell.execute_reply.started":"2023-10-29T11:12:01.943021Z","shell.execute_reply":"2023-10-29T11:12:01.950142Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"[Return to Contents](#contents)","metadata":{}},{"cell_type":"markdown","source":"### Viewing the contents of a .pth file\n\nA `.pth` file is typically a PyTorch model checkpoint or saved tensors. To view the contents of a `.pth` file, you'll first load it using PyTorch and then inspect its contents. Here's how you can do it:\n\n1. **Install and Import Necessary Libraries**:\n   First, ensure you have PyTorch installed. If not, install it using pip:\n   ```\n   pip install torch torchvision\n   ```\n   Then, import the necessary libraries in your Python script or Jupyter Notebook:\n   ```python\n   import torch\n   ```\n\n2. **Load the `.pth` File**:\n   Use the `torch.load()` function to load the contents of the `.pth` file:\n   ```python\n   checkpoint = torch.load('path_to_file.pth', map_location='cpu')\n   ```\n\n   Here, `map_location='cpu'` ensures that the tensors are loaded onto the CPU. This is particularly useful if the checkpoint was saved on a GPU and you're trying to load it on a machine without a GPU.\n\n3. **Inspect the Contents**:\n   The content and structure of the checkpoint file depend on how it was saved.\n\n   - If it's a model checkpoint, it might contain the model's state_dict or even optimizer state, and other metadata:\n     ```python\n     print(checkpoint.keys())\n     ```\n\n   - If the `.pth` file just contains a model's state_dict, you can view it like this:\n     ```python\n     for param_tensor in checkpoint:\n         print(param_tensor, checkpoint[param_tensor].size())\n     ```\n\n   - If it's a simple tensor saved in the `.pth` file, you can print it directly:\n     ```python\n     print(checkpoint)\n     ```\n\n4. **(Optional) Load into Model**:\n   If the `.pth` file is a saved model checkpoint, and you have the corresponding model architecture defined, you can load the weights into the model like this:\n   ```python\n   model = YourModelArchitecture()  # Replace with your model's architecture\n   model.load_state_dict(checkpoint['state_dict'])  # Assuming state_dict is a key in your checkpoint\n   ```\n\nRemember, the exact structure of the checkpoint will depend on how it was saved. It's important to know the structure (e.g., whether it contains just the model weights or other information) when trying to load and inspect it.","metadata":{}},{"cell_type":"markdown","source":"<a id =  \"RNA_LB\"></a>\n## RNA Starter [0.186 LB]\n### (This is based on [RNA starter 0.186 LB](https://www.kaggle.com/code/iafoss/rna-starter-0-186-lb) by [Iafoss](https://www.kaggle.com/iafoss))\n\nProvides a simple baseline that may be used as a starting point for further experiments. Improvment of the baseline may include:\n\n- Use of proper loss function to incorporate SN_filter = 0 samples as well as reactivity errors into training\n- Model improvement and use of additional data, e.g. Ribonanza_bpp_files\n\nFinally, working on this competition keep in mind that train/public LB have different sequence length distribution from private LB, i.e. 115-206 vs. 207-457. Therefore, to avoid a strong shakeup one may need to look into performance vs. sequence length end ensure generalizability.","metadata":{}},{"cell_type":"code","source":"import pandas as pd\nimport os, gc\nimport numpy as np\nfrom sklearn.model_selection import KFold\n\nimport torch\nimport torch.nn as nn\nimport torch.nn.functional as F\nfrom torch.utils.data import Dataset, DataLoader","metadata":{"execution":{"iopub.status.busy":"2023-10-30T08:35:34.323091Z","iopub.execute_input":"2023-10-30T08:35:34.323359Z","iopub.status.idle":"2023-10-30T08:35:38.935889Z","shell.execute_reply.started":"2023-10-30T08:35:34.323333Z","shell.execute_reply":"2023-10-30T08:35:38.935156Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Explanation**\n\n1. `import os, gc`: These imports are for the standard Python libraries `os` (for interacting with the file system) and `gc` (for garbage collection), respectively.\n2. `from sklearn.model_selection import KFold`: This imports the `KFold` class from the scikit-learn library. `KFold` is a method for cross-validation, which is used to split data into multiple subsets (folds) for training and testing machine learning models.\n3. `import torch`, `import torch.nn as nn`, `import torch.nn.functional as F`: \n    - `torch` is the core PyTorch library.\n    - `torch.nn` is the module containing neural network-related classes and functions.\n    - `torch.nn.functional` is used for various activation functions, loss functions, and other operations in PyTorch.\n4. `from torch.utils.data import Dataset, DataLoader`: These imports are for PyTorch's data handling utilities.\n   - `Dataset` is an abstract class that you can subclass to create custom datasets for training and evaluation.\n   - `DataLoader` is a class used to load data from a dataset into batches, making it suitable for training deep learning models in mini-batches.\n\nOverall, these imports set up the environment for working with data and deep learning using libraries like Pandas, NumPy, scikit-learn, and PyTorch. The specific functionalities provided by these libraries will be used throughout the code, but the code implementation itself is not shown in this snippet.","metadata":{}},{"cell_type":"code","source":"# Fix fastai bug to enable fp16 training with dictionaries\n\nimport torch\nfrom fastai.vision.all import *\ndef flatten(o):\n    \"Concatenate all collections and items as a generator\"\n    for item in o:\n        if isinstance(o, dict): yield o[item]; continue\n        elif isinstance(item, str): yield item; continue\n        try: yield from flatten(item)\n        except TypeError: yield item\n\nfrom torch.cuda.amp import GradScaler, autocast\n@delegates(GradScaler)\nclass MixedPrecision(Callback):\n    \"Mixed precision training using Pytorch's `autocast` and `GradScaler`\"\n    order = 10\n    def __init__(self, **kwargs): self.kwargs = kwargs\n    def before_fit(self): \n        self.autocast,self.learn.scaler,self.scales = autocast(),GradScaler(**self.kwargs),L()\n    def before_batch(self): self.autocast.__enter__()\n    def after_pred(self):\n        if next(flatten(self.pred)).dtype==torch.float16: self.learn.pred = to_float(self.pred)\n    def after_loss(self): self.autocast.__exit__(None, None, None)\n    def before_backward(self): self.learn.loss_grad = self.scaler.scale(self.loss_grad)\n    def before_step(self):\n        \"Use `self` as a fake optimizer. `self.skipped` will be set to True `after_step` if gradients overflow. \"\n        self.skipped=True\n        self.scaler.step(self)\n        if self.skipped: raise CancelStepException()\n        self.scales.append(self.scaler.get_scale())\n    def after_step(self): self.learn.scaler.update()\n\n    @property \n    def param_groups(self): \n        \"Pretend to be an optimizer for `GradScaler`\"\n        return self.opt.param_groups\n    def step(self, *args, **kwargs): \n        \"Fake optimizer step to detect whether this batch was skipped from `GradScaler`\"\n        self.skipped=False\n    def after_fit(self): self.autocast,self.learn.scaler,self.scales = None,None,None\n        \nimport fastai\nfastai.callback.fp16.MixedPrecision = MixedPrecision","metadata":{"execution":{"iopub.status.busy":"2023-10-30T08:40:49.543595Z","iopub.execute_input":"2023-10-30T08:40:49.544112Z","iopub.status.idle":"2023-10-30T08:40:51.456647Z","shell.execute_reply.started":"2023-10-30T08:40:49.544078Z","shell.execute_reply":"2023-10-30T08:40:51.455673Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Explanation**\n\nThis is a patch to the `MixedPrecision` callback class in the `fastai` library to enable mixed precision (fp16) training with dictionaries. Modifications have been made to handle cases where predictions (`pred`) might be dictionaries and implemented the necessary changes to make it work with the `autocast` and `GradScaler` functions provided by PyTorch for mixed precision training.\n\nIf fastai` is already installed, this code can be run to patch the existing `MixedPrecision` class. After doing so, any calls to `MixedPrecision` within `fastai` will use the modified version, enabling support for dictionaries.","metadata":{}},{"cell_type":"code","source":"def seed_everything(seed):\n    random.seed(seed)\n    os.environ['PYTHONHASHSEED'] = str(seed)\n    np.random.seed(seed)\n    torch.manual_seed(seed)\n    torch.cuda.manual_seed(seed)\n    torch.backends.cudnn.deterministic = True\n    torch.backends.cudnn.benchmark = True","metadata":{"execution":{"iopub.status.busy":"2023-10-30T08:41:54.893687Z","iopub.execute_input":"2023-10-30T08:41:54.894072Z","iopub.status.idle":"2023-10-30T08:41:54.899817Z","shell.execute_reply.started":"2023-10-30T08:41:54.894042Z","shell.execute_reply":"2023-10-30T08:41:54.898902Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fname = 'example0'\nPATH = '/kaggle/input/stanford-ribonanza-rna-folding-converted/'\nOUT = './'\nbs = 256\nnum_workers = 2\nSEED = 2023\nnfolds = 4\ndevice = 'cuda' if torch.cuda.is_available() else 'cpu'","metadata":{"execution":{"iopub.status.busy":"2023-10-30T08:42:16.202983Z","iopub.execute_input":"2023-10-30T08:42:16.203331Z","iopub.status.idle":"2023-10-30T08:42:16.208582Z","shell.execute_reply.started":"2023-10-30T08:42:16.203303Z","shell.execute_reply":"2023-10-30T08:42:16.207568Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Explanation**\n\nThe function `seed_everything` aims to set seeds for various random number generators to ensure that the code produces the same results every time it's run, given the same seed. This is especially useful when you want reproducibility in your experiments.\n\n1. `random.seed(seed)`: Sets the seed for Python's built-in random module.\n2. `os.environ['PYTHONHASHSEED'] = str(seed)`: This is used to ensure the reproducibility of hash functions between different runs and Python processes. It's especially important when dealing with Python dictionaries or sets.\n3. `np.random.seed(seed)`: Sets the seed for NumPy's random number generator.\n4. `torch.manual_seed(seed)`: Sets the seed for generating random numbers for both CPU and CUDA with PyTorch.\n5. `torch.cuda.manual_seed(seed)`: Specifically sets the seed for generating random numbers in CUDA. Note that if you have multiple GPUs, `torch.cuda.manual_seed_all(seed)` can be used to set the seed on all devices.\n6. `torch.backends.cudnn.deterministic = True`: Ensures that every time you run your code on a GPU, it uses the same algorithms, providing deterministic results.\n7. `torch.backends.cudnn.benchmark = True`: Enables the inbuilt cuDNN auto-tuner to find the best algorithm to use for the hardware. This can optimize runtime, but it might lead to slight variations in results due to non-deterministic algorithms. If you prioritize reproducibility over speed, you might consider setting this to `False`.\n\nWhen using this function:\n\n- Always call `seed_everything` at the beginning of your script or notebook before any random operations to ensure reproducibility.\n- Keep in mind that there are other factors outside of the code (like GPU type and driver version) that can affect the reproducibility of results in deep learning experiments. \n- If using other libraries (e.g., TensorFlow), you might need to extend the function to seed their random generators as well.","metadata":{}},{"cell_type":"code","source":"fname = 'example0'\nPATH = '/kaggle/input/stanford-ribonanza-rna-folding-converted/'\nOUT = './'\nbs = 256\nnum_workers = 2\nSEED = 2023\nnfolds = 4\ndevice = 'cuda' if torch.cuda.is_available() else 'cpu'","metadata":{"execution":{"iopub.status.busy":"2023-10-30T08:42:29.139151Z","iopub.execute_input":"2023-10-30T08:42:29.139886Z","iopub.status.idle":"2023-10-30T08:42:29.144567Z","shell.execute_reply.started":"2023-10-30T08:42:29.139855Z","shell.execute_reply":"2023-10-30T08:42:29.143728Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Explanation**\n\nSeveral variables that are often used as configuration or setup parameters in a typical machine learning or data processing script have been set-up.\n\n1. `fname`: This appears to be the filename (or part of it) that you might be working with. In this case, the name is 'example0'.\n\n2. `PATH`: This is the path to your dataset or input files.\n\n3. `OUT`: This specifies the directory where you'll save any output files or models. It's set to the current directory ('./').\n\n4. `bs`: This stands for \"batch size\". \n\n5. `num_workers`: This is often used with PyTorch's DataLoader. It specifies the number of subprocesses to use for data loading. More workers can speed up data loading at the cost of using more memory. \n\n6. `SEED`: This is a seed value, which can be used to ensure reproducibility in random operations.\n\n7. `nfolds`: This refers to the number of folds in k-fold cross-validation, a technique used to evaluate machine learning models by dividing the data into 'k' subsets (or \"folds\") and training/testing on them in a rotating fashion.\n\n8. `device`: This determines where your computations will take place in PyTorch - on a CUDA-enabled GPU ('cuda') if one is available, or on the CPU ('cpu') otherwise.","metadata":{}},{"cell_type":"markdown","source":"### Data\n\nThe primary training data is provided in train_data.csv, which contains 821,840 RNA sequences and the corresponding reactivity measurements with 2A3_MaP and DMS_MaP methods. The reactivity is reported in columns reactivity_0001 - reactivity_0206 and is set to NaN for the first 26 and the last 21 nucleotides as well as padding for sequences shorter than 206. For faster loading and effective RAM use, the data was convertedinto a float32 parquet file.\n\nEvaluation in the competition is performed only on samples with SN_filter = 1 for both measurement methods. In this example, training is only performed on samples with SN_filter = 1. This gives a noticeable CV boost but uses only 1/4 of the data (i.e. training on noisy SN_filter = 0 data degrades performance). A proper consideration of all data as well as reactivity errors may boost the performance.\n\nThis example uses a simple CV Kfold split. However, given a mismatch in the RNA length between train/public LB vs. private LB data, it may be important to verify the effect of the sequence length to avoid a significant shakeup at the private LB.\n\nOne of the tricks in the NLP community, which is used here, is length matching batch sampling: composing batches of samples of approximately the same length to minimize the overhead caused by padding tokens.","metadata":{}},{"cell_type":"code","source":"class RNA_Dataset(Dataset):\n    def __init__(self, df, mode='train', seed=2023, fold=0, nfolds=4, \n                 mask_only=False, **kwargs):\n        self.seq_map = {'A':0,'C':1,'G':2,'U':3}\n        self.Lmax = 206\n        df['L'] = df.sequence.apply(len)\n        df_2A3 = df.loc[df.experiment_type=='2A3_MaP']\n        df_DMS = df.loc[df.experiment_type=='DMS_MaP']\n        \n        split = list(KFold(n_splits=nfolds, random_state=seed, \n                shuffle=True).split(df_2A3))[fold][0 if mode=='train' else 1]\n        df_2A3 = df_2A3.iloc[split].reset_index(drop=True)\n        df_DMS = df_DMS.iloc[split].reset_index(drop=True)\n        \n        m = (df_2A3['SN_filter'].values > 0) & (df_DMS['SN_filter'].values > 0)\n        df_2A3 = df_2A3.loc[m].reset_index(drop=True)\n        df_DMS = df_DMS.loc[m].reset_index(drop=True)\n        \n        self.seq = df_2A3['sequence'].values\n        self.L = df_2A3['L'].values\n        \n        self.react_2A3 = df_2A3[[c for c in df_2A3.columns if \\\n                                 'reactivity_0' in c]].values\n        self.react_DMS = df_DMS[[c for c in df_DMS.columns if \\\n                                 'reactivity_0' in c]].values\n        self.react_err_2A3 = df_2A3[[c for c in df_2A3.columns if \\\n                                 'reactivity_error_0' in c]].values\n        self.react_err_DMS = df_DMS[[c for c in df_DMS.columns if \\\n                                'reactivity_error_0' in c]].values\n        self.sn_2A3 = df_2A3['signal_to_noise'].values\n        self.sn_DMS = df_DMS['signal_to_noise'].values\n        self.mask_only = mask_only\n        \n    def __len__(self):\n        return len(self.seq)  \n    \n    def __getitem__(self, idx):\n        seq = self.seq[idx]\n        if self.mask_only:\n            mask = torch.zeros(self.Lmax, dtype=torch.bool)\n            mask[:len(seq)] = True\n            return {'mask':mask},{'mask':mask}\n        seq = [self.seq_map[s] for s in seq]\n        seq = np.array(seq)\n        mask = torch.zeros(self.Lmax, dtype=torch.bool)\n        mask[:len(seq)] = True\n        seq = np.pad(seq,(0,self.Lmax-len(seq)))\n        \n        react = torch.from_numpy(np.stack([self.react_2A3[idx],\n                                           self.react_DMS[idx]],-1))\n        react_err = torch.from_numpy(np.stack([self.react_err_2A3[idx],\n                                               self.react_err_DMS[idx]],-1))\n        sn = torch.FloatTensor([self.sn_2A3[idx],self.sn_DMS[idx]])\n        \n        return {'seq':torch.from_numpy(seq), 'mask':mask}, \\\n               {'react':react, 'react_err':react_err,\n                'sn':sn, 'mask':mask}","metadata":{"execution":{"iopub.status.busy":"2023-10-30T08:42:43.646288Z","iopub.execute_input":"2023-10-30T08:42:43.647026Z","iopub.status.idle":"2023-10-30T08:42:43.663633Z","shell.execute_reply.started":"2023-10-30T08:42:43.646994Z","shell.execute_reply":"2023-10-30T08:42:43.662631Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Explanation**\n\nThis is a `Dataset` class definition named `RNA_Dataset`. It is customized for handling RNA-related data.\n\n1. **Initialization (`__init__`)**:\n    - The `seq_map` dictionary maps RNA nucleotides to numerical values.\n    - Some pre-processing steps are applied to the input dataframe `df`. The dataframe contains details about RNA sequences and their reactivity profiles. Sequences are grouped by `experiment_type`.\n    - The dataset is further split into training and validation sets using `KFold` cross-validation, based on the `mode` argument.\n    - Only sequences that pass certain filtering criteria are retained.\n    - Data arrays for sequence, reactivity, error, and signal-to-noise ratio are prepared.\n\n2. **Length (`__len__`)**:\n    - This function returns the number of sequences in the dataset.\n\n3. **Get Item (`__getitem__`)**:\n    - Given an index `idx`, this function returns the corresponding RNA sequence and its associated reactivity information.\n    - The function outputs a dictionary with sequence data and another with reactivity-related data.\n    - Sequences are converted to numerical format using the `seq_map`.\n    - If `mask_only` is `True`, the function only returns a mask representing the valid length of the sequence.\n    - Otherwise, it returns the processed sequence, mask, reactivity profiles, associated errors, and signal-to-noise ratios. \n\nTo use this `Dataset` class:\n\n- You'd initialize it by passing a DataFrame with RNA sequence data and relevant configurations.\n- With a DataLoader, you can then use it to retrieve batches of data for training or evaluation in PyTorch. \n\nThis is a great representation of how PyTorch allows for custom datasets to be created and tailored to specific needs, ensuring flexibility in handling various types of data.","metadata":{}},{"cell_type":"code","source":"class LenMatchBatchSampler(torch.utils.data.BatchSampler):\n    def __iter__(self):\n        buckets = [[]] * 100\n        yielded = 0\n\n        for idx in self.sampler:\n            s = self.sampler.data_source[idx]\n            if isinstance(s,tuple): L = s[0][\"mask\"].sum()\n            else: L = s[\"mask\"].sum()\n            L = max(1,L // 16) \n            if len(buckets[L]) == 0:  buckets[L] = []\n            buckets[L].append(idx)\n            \n            if len(buckets[L]) == self.batch_size:\n                batch = list(buckets[L])\n                yield batch\n                yielded += 1\n                buckets[L] = []\n                \n        batch = []\n        leftover = [idx for bucket in buckets for idx in bucket]\n\n        for idx in leftover:\n            batch.append(idx)\n            if len(batch) == self.batch_size:\n                yielded += 1\n                yield batch\n                batch = []\n\n        if len(batch) > 0 and not self.drop_last:\n            yielded += 1\n            yield batch\n            \ndef dict_to(x, device='cuda'):\n    return {k:x[k].to(device) for k in x}\n\ndef to_device(x, device='cuda'):\n    return tuple(dict_to(e,device) for e in x)","metadata":{"execution":{"iopub.status.busy":"2023-10-30T08:43:13.277372Z","iopub.execute_input":"2023-10-30T08:43:13.277713Z","iopub.status.idle":"2023-10-30T08:43:13.288113Z","shell.execute_reply.started":"2023-10-30T08:43:13.277685Z","shell.execute_reply":"2023-10-30T08:43:13.287205Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Explanation**\n\nThis is a custom batch sampler and some utility functions for PyTorch.\n\n1. **`LenMatchBatchSampler`**:\n   This is a custom batch sampler. In many NLP and sequence tasks, sequences of similar lengths are batched together to minimize the amount of padding, which can accelerate training.\n   \n    - The sampler organizes data points (sequences) into buckets based on their length (number of non-masked positions).\n    - When a bucket has enough data points to form a batch (as specified by `self.batch_size`), it yields a batch of indices.\n    - After processing all data points, the sampler yields any remaining batches from the leftover items.\n    - This sampler prioritizes batching sequences of similar length together, which is efficient especially when working with RNNs as it minimizes computations on padded tokens.\n\n2. **`dict_to`**:\n   A utility function to send all tensors within a dictionary to a specified device (e.g., a GPU).\n\n3. **`to_device`**:\n   This function works similarly to `dict_to` but handles a tuple of dictionaries. It sends all tensors within each dictionary in the tuple to the specified device.\n\nWhen training a model, it's common to use utility functions like these to ensure that the data and the model are on the same device (usually a GPU for faster training). The custom batch sampler is particularly useful when dealing with variable-length sequences, ensuring computational efficiency by minimizing the amount of unnecessary padding.","metadata":{}},{"cell_type":"code","source":"class DeviceDataLoader:\n    def __init__(self, dataloader, device='cuda'):\n        self.dataloader = dataloader\n        self.device = device\n    \n    def __len__(self):\n        return len(self.dataloader)\n    \n    def __iter__(self):\n        for batch in self.dataloader:\n            yield tuple(dict_to(x, self.device) for x in batch)","metadata":{"execution":{"iopub.status.busy":"2023-10-30T08:43:26.230807Z","iopub.execute_input":"2023-10-30T08:43:26.231142Z","iopub.status.idle":"2023-10-30T08:43:26.237117Z","shell.execute_reply.started":"2023-10-30T08:43:26.231107Z","shell.execute_reply":"2023-10-30T08:43:26.236097Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Explanation**\n\nThe `DeviceDataLoader` class is a wrapper around a standard PyTorch `DataLoader`. Its main purpose is to ensure that each batch of data fetched from the `DataLoader` is immediately sent to a specified device (usually a GPU, denoted by `'cuda'`). This is particularly useful when working with deep learning models in PyTorch, as both the model and its input data need to be on the same device for computation.\n\n1. **Initialization (`__init__`)**:\n    - Takes in a PyTorch `DataLoader` object and an optional device specification (default is `'cuda'`).\n    - Stores these as instance variables.\n\n2. **Length (`__len__`)**:\n    - Simply returns the length of the original `DataLoader`. This would correspond to the number of batches the `DataLoader` can provide.\n\n3. **Iterator (`__iter__`)**:\n    - Defines how the `DeviceDataLoader` should iterate over batches.\n    - For each batch fetched from the original `DataLoader`, it uses the previously discussed `dict_to` function to send each part of the batch to the specified device.\n    - Yields the device-transferred batch.\n\nBy wrapping your data loader with `DeviceDataLoader`, the data preparation step is streamlined you make sure each batch is on the correct device right when it is fetched, before it is passed to the model. This minimizes the chances of \"device mismatch\" errors and centralizes the device handling logic, making the training loop cleaner and less error-prone.","metadata":{}},{"cell_type":"markdown","source":"### Model\nTransformers are ideal for the considered task because they naturally capture long dependences in RNN, which define the secondary structure of the molecule and the corresponding chemical reactivity. For illustration purposes, a simple S-size transformer model is provided below.","metadata":{}},{"cell_type":"code","source":"class SinusoidalPosEmb(nn.Module):\n    def __init__(self, dim=16, M=10000):\n        super().__init__()\n        self.dim = dim\n        self.M = M\n\n    def forward(self, x):\n        device = x.device\n        half_dim = self.dim // 2\n        emb = math.log(self.M) / half_dim\n        emb = torch.exp(torch.arange(half_dim, device=device) * (-emb))\n        emb = x[...,None] * emb[None,...]\n        emb = torch.cat((emb.sin(), emb.cos()), dim=-1)\n        return emb","metadata":{"execution":{"iopub.status.busy":"2023-10-30T08:43:40.583090Z","iopub.execute_input":"2023-10-30T08:43:40.583979Z","iopub.status.idle":"2023-10-30T08:43:40.590307Z","shell.execute_reply.started":"2023-10-30T08:43:40.583945Z","shell.execute_reply":"2023-10-30T08:43:40.589337Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Explanation**\n\nThe `SinusoidalPosEmb` class is an implementation of sinusoidal positional embeddings, which is commonly used in Transformer architectures to provide information about the position of tokens in a sequence. Let's dissect the module:\n\n**Purpose**:\nPositional embeddings are crucial for architectures like the Transformer because they operate on set-like data (the order doesn't inherently matter), so positional information isn't preserved. By adding positional embeddings to the input data, the network can understand and use the order of the tokens.\n\n**Implementation**:\n\n1. **Initialization (`__init__`)**:\n    - `dim`: The dimensionality of the positional embeddings.\n    - `M`: A constant used to control the frequency of the sinusoids. Default is 10,000, as typically used in Transformer literature.    \n2. **Forward Method (`forward`)**:\n    - `x`: The input tensor whose last dimension will represent positions for which positional embeddings need to be calculated.\n    - The method calculates sinusoidal positional embeddings for the input tensor `x`.\n    - `emb` values are calculated such that they represent a series of frequencies in a geometric progression. This is used to produce a unique combination of sine and cosine values for each position.\n    - For each position in `x`, sine values are computed for the first half of dimensions and cosine values for the second half, effectively encoding the positional information into the embedding space.\n\n**Usage**:\nThis can be used as part of the input pipeline for a Transformer or any other architecture where the inherent order of the data isn't captured. You'd typically add these positional embeddings to the token embeddings before feeding them into the model.\n\n**Example**:\nIf you have a sequence of tokens represented by embeddings, you can use this module to add positional information to them. For instance, if your token embeddings are of dimension 16 and you have a sequence length of 10, the input tensor `x` to this module would be of shape `(10,)`, and the output would be of shape `(10, 16)`. This output can then be element-wise added to the token embeddings of the sequence to include positional information.","metadata":{}},{"cell_type":"code","source":"class RNA_Model(nn.Module):\n    def __init__(self, dim=192, depth=12, head_size=32, **kwargs):\n        super().__init__()\n        self.emb = nn.Embedding(4,dim)\n        self.pos_enc = SinusoidalPosEmb(dim)\n        self.transformer = nn.TransformerEncoder(\n            nn.TransformerEncoderLayer(d_model=dim, nhead=dim//head_size, dim_feedforward=4*dim,\n                dropout=0.1, activation=nn.GELU(), batch_first=True, norm_first=True), depth)\n        self.proj_out = nn.Linear(dim,2)\n    \n    def forward(self, x0):\n        mask = x0['mask']\n        Lmax = mask.sum(-1).max()\n        mask = mask[:,:Lmax]\n        x = x0['seq'][:,:Lmax]\n        \n        pos = torch.arange(Lmax, device=x.device).unsqueeze(0)\n        pos = self.pos_enc(pos)\n        x = self.emb(x)\n        x = x + pos\n        \n        x = self.transformer(x, src_key_padding_mask=~mask)\n        x = self.proj_out(x)\n        \n        return x","metadata":{"execution":{"iopub.status.busy":"2023-10-30T08:43:50.391769Z","iopub.execute_input":"2023-10-30T08:43:50.392558Z","iopub.status.idle":"2023-10-30T08:43:50.401074Z","shell.execute_reply.started":"2023-10-30T08:43:50.392527Z","shell.execute_reply":"2023-10-30T08:43:50.400133Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Explanation**\n\nA PyTorch neural network model called `RNA_Model` is defined. \n\n1. **Initialization**: In the constructor (`__init__` method), the architecture of the model is defined.\n\n    - `dim`: This parameter sets the dimension of the embeddings and hidden layers in the model. It's set to 192 by default.    \n    - `depth`: This parameter determines the number of layers in the Transformer encoder. It's set to 12 by default.    \n    - `head_size`: This parameter controls the size of the attention heads in the Transformer encoder. It's set to 32 by default.\n2. **Embedding Layer**: The model starts with an embedding layer (`self.emb`) that embeds input sequences. It's designed for sequences with four different tokens (perhaps representing different nucleotides in RNA sequences) and embeds them into a `dim`-dimensional space.\n3. **Positional Encoding**: The model uses a sinusoidal positional encoding (`SinusoidalPosEmb`) to add positional information to the embedded sequences. This helps the model understand the order of tokens in the input sequence.\n4. **Transformer Encoder**: The core of the model is a stack of Transformer encoder layers (`nn.TransformerEncoder`). It uses `depth` layers of `nn.TransformerEncoderLayer`, each with the following characteristics:\n    \n    - `d_model`: The dimension of model embeddings, which is set to `dim`.    \n    - `nhead`: The number of attention heads. It's calculated as `dim // head_size`.    \n    - `dim_feedforward`: The dimension of the feedforward network inside the Transformer layers, set to `4 * dim`.    \n    - `dropout`: Dropout rate set to 0.1.    \n    - `activation`: The activation function used within the model, which is the GELU (Gaussian Error Linear Unit) activation function.    \n    - `batch_first` and `norm_first`: These parameters specify whether batch normalization should be applied before or after the attention and feedforward layers.\n5. **Output Projection Layer**: After the Transformer encoder, there's a linear projection layer (`self.proj_out`) that maps the hidden representations back to a 2-dimensional output. This suggests that the model is designed for a binary classification task where it needs to output two classes or values.\n6. **Forward Method**: In the `forward` method, the input `x0` is expected to be a dictionary containing two keys:\n   \n   - `'mask'`: A mask for sequence padding, which appears to be used to handle variable-length sequences.   \n   - `'seq'`: The input sequences.\n   \n   The code processes the input as follows:   \n   - Trims the mask and sequences to the maximum length found in the batch.   \n   - Computes positional encodings and adds them to the embedded sequences.   \n   - Passes the sequences through the Transformer encoder, using the mask to mask out padded elements.   \n   - Applies the projection layer to obtain the final output.\n\nOverall, this model is a variant of the Transformer architecture customized for a specific sequence classification task involving RNA sequences. This model can be used for tasks like RNA sequence classification, and you can adjust the hyperparameters such as `dim`, `depth`, and `head_size` as needed.","metadata":{}},{"cell_type":"markdown","source":"### Loss & Metric\nThe metric accumulates all predictions and then performs average to be consistent with the competition metric. However, the difference with a simple batch-based average is negligible.","metadata":{}},{"cell_type":"code","source":"def loss(pred,target):\n    p = pred[target['mask'][:,:pred.shape[1]]]\n    y = target['react'][target['mask']].clip(0,1)\n    loss = F.l1_loss(p, y, reduction='none')\n    loss = loss[~torch.isnan(loss)].mean()\n    \n    return loss","metadata":{"execution":{"iopub.status.busy":"2023-10-30T08:44:06.746774Z","iopub.execute_input":"2023-10-30T08:44:06.747658Z","iopub.status.idle":"2023-10-30T08:44:06.752997Z","shell.execute_reply.started":"2023-10-30T08:44:06.747624Z","shell.execute_reply":"2023-10-30T08:44:06.752032Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Explanation**\n\nDefines a custom loss function called `loss` involving predictions (`pred`) and target values (`target`) with associated masks.\n\n1. **Input Parameters**:\n   - `pred`: The predicted values produced by the model.\n   - `target`: A dictionary containing various elements, including masks and target values.\n2. **Extracting Relevant Predictions and Targets**:\n   - `p = pred[target['mask'][:,:pred.shape[1]]]`: This line extracts predicted values (`p`) based on the mask from the `target` dictionary. The mask appears to be used to select specific elements from the predictions.   \n   - `y = target['react'][target['mask']].clip(0,1)`: This line extracts target values (`y`) from the `'react'` key in the `target` dictionary, again based on the mask. Additionally, it clips the values between 0 and 1 using `.clip(0, 1)` to ensure they fall within that range.\n3. **Loss Calculation**:\n   - `loss = F.l1_loss(p, y, reduction='none')`: This line calculates the element-wise L1 loss (mean absolute error) between the predicted values `p` and the clipped target values `y`. The `reduction='none'` argument means that it doesn't compute the mean over all elements yet, leaving the loss as a tensor of the same shape as `p` and `y`.\n4. **Handling NaN Values**:\n   - `loss = loss[~torch.isnan(loss)].mean()`: This line removes elements with NaN (Not-a-Number) values from the calculated loss tensor and then computes the mean of the remaining elements. This step is necessary to handle cases where NaN values might occur in the loss due to specific conditions in the data.\n5. **Return Value**:\n   - The computed loss value is returned as the final result.\n\nIn summary, this custom loss function calculates the mean absolute error (L1 loss) between selected elements of the predicted values and the clipped target values while handling NaN values that might arise during the computation. This loss function is useful when working with tasks where only specific elements of the predictions and targets are relevant and when you want to penalize differences between predicted and target values in an absolute manner.","metadata":{}},{"cell_type":"code","source":"class MAE(Metric):\n    def __init__(self): \n        self.reset()\n        \n    def reset(self): \n        self.x,self.y = [],[]\n        \n    def accumulate(self, learn):\n        x = learn.pred[learn.y['mask'][:,:learn.pred.shape[1]]]\n        y = learn.y['react'][learn.y['mask']].clip(0,1)\n        self.x.append(x)\n        self.y.append(y)\n\n    @property\n    def value(self):\n        x,y = torch.cat(self.x,0),torch.cat(self.y,0)\n        loss = F.l1_loss(x, y, reduction='none')\n        loss = loss[~torch.isnan(loss)].mean()\n        return loss","metadata":{"execution":{"iopub.status.busy":"2023-10-30T08:44:16.244631Z","iopub.execute_input":"2023-10-30T08:44:16.244977Z","iopub.status.idle":"2023-10-30T08:44:16.252810Z","shell.execute_reply.started":"2023-10-30T08:44:16.244947Z","shell.execute_reply":"2023-10-30T08:44:16.251863Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Explanation**\n\nDefines a custom metric class called `MAE` (Mean Absolute Error) which is used to calculate the mean absolute error between predicted values and target values while handling masks and NaN values.\n\n1. **Initialization**:\n   - In the `__init__` method, the class initializes its state by calling the `reset` method.\n2. **`reset` Method**:\n   - The `reset` method initializes two empty lists, `self.x` and `self.y`. These lists are used to accumulate predicted (`x`) and target (`y`) values during training or evaluation.\n3. **`accumulate` Method**:\n   - The `accumulate` method is used to accumulate predicted and target values for later calculation of the mean absolute error. It takes a `learn` argument, which presumably contains predicted (`learn.pred`) and target (`learn.y`) values.   \n   - Inside this method:\n     - `x` is calculated by extracting predicted values based on the mask from `learn.pred`.\n     - `y` is calculated by extracting target values based on the mask from `learn.y`, and it's clipped between 0 and 1.\n     - `x` and `y` are appended to the respective lists `self.x` and `self.y`.\n4. **`value` Property**:\n   - The `value` property calculates the mean absolute error (MAE) between all accumulated predicted (`self.x`) and target (`self.y`) values.   \n   - It first concatenates the accumulated `self.x` and `self.y` lists along the 0th dimension to create tensors `x` and `y`.   \n   - It then calculates the element-wise L1 loss (mean absolute error) between `x` and `y` using `F.l1_loss`. The `reduction='none'` argument means that it doesn't compute the mean over all elements yet, leaving the loss as a tensor of the same shape as `x` and `y`.   \n   - It removes elements with NaN values from the loss tensor and computes the mean of the remaining elements, handling cases where NaN values might occur during the computation.","metadata":{}},{"cell_type":"markdown","source":"### Training","metadata":{}},{"cell_type":"code","source":"seed_everything(SEED)\nos.makedirs(OUT, exist_ok=True)\ndf = pd.read_parquet(os.path.join(PATH,'train_data.parquet'))\n\nfor fold in [0]: # running multiple folds at kaggle may cause OOM\n    ds_train = RNA_Dataset(df, mode='train', fold=fold, nfolds=nfolds)\n    ds_train_len = RNA_Dataset(df, mode='train', fold=fold, \n                nfolds=nfolds, mask_only=True)\n    sampler_train = torch.utils.data.RandomSampler(ds_train_len)\n    len_sampler_train = LenMatchBatchSampler(sampler_train, batch_size=bs,\n                drop_last=True)\n    dl_train = DeviceDataLoader(torch.utils.data.DataLoader(ds_train, \n                batch_sampler=len_sampler_train, num_workers=num_workers,\n                persistent_workers=True), device)\n\n    ds_val = RNA_Dataset(df, mode='eval', fold=fold, nfolds=nfolds)\n    ds_val_len = RNA_Dataset(df, mode='eval', fold=fold, nfolds=nfolds, \n               mask_only=True)\n    sampler_val = torch.utils.data.SequentialSampler(ds_val_len)\n    len_sampler_val = LenMatchBatchSampler(sampler_val, batch_size=bs, \n               drop_last=False)\n    dl_val= DeviceDataLoader(torch.utils.data.DataLoader(ds_val, \n               batch_sampler=len_sampler_val, num_workers=num_workers), device)\n    gc.collect()\n\n    data = DataLoaders(dl_train,dl_val)\n    model = RNA_Model()   \n    model = model.to(device)\n    learn = Learner(data, model, loss_func=loss,cbs=[GradientClip(3.0)],\n                metrics=[MAE()]).to_fp16() \n    #fp16 doesn't help at P100 but gives x1.6-1.8 speedup at modern hardware\n\n    learn.fit_one_cycle(32, lr_max=5e-4, wd=0.05, pct_start=0.02)\n    torch.save(learn.model.state_dict(),os.path.join(OUT,f'{fname}_{fold}.pth'))\n    gc.collect()","metadata":{"execution":{"iopub.status.busy":"2023-10-30T08:54:39.489068Z","iopub.execute_input":"2023-10-30T08:54:39.489748Z","iopub.status.idle":"2023-10-30T11:43:01.743182Z","shell.execute_reply.started":"2023-10-30T08:54:39.489711Z","shell.execute_reply":"2023-10-30T11:43:01.742042Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Explanation**\n\n1. **Setting Seed and Creating Output Directory**:\n   - `seed_everything(SEED)`: This function is used to set a random seed for reproducibility.\n   - `os.makedirs(OUT, exist_ok=True)`: It creates the output directory specified by the `OUT` variable if it doesn't already exist. This directory is likely used to save model checkpoints.\n2. **Loading Data**:\n   - `df = pd.read_parquet(os.path.join(PATH,'train_data.parquet'))`: This line reads a parquet file containing training data, presumably RNA sequence data, into a Pandas DataFrame.\n3. **Training Loop (for a single fold)**:\n   - The training loop appears to run for a single fold, as indicated by the loop `[0]`.\n4. **Creating Training Dataset and DataLoader**:\n   - `ds_train` and `ds_train_len` are created as instances of the `RNA_Dataset` class, presumably designed for handling the training data. `ds_train` is for the full training dataset, while `ds_train_len` appears to be used for calculating the length of sequences.\n   - `sampler_train` and `len_sampler_train` are used to create a custom data sampler and batch sampler for training data. These samplers seem to be related to handling sequences of variable lengths.\n   - `dl_train` is a `DeviceDataLoader` that wraps a PyTorch `DataLoader` for the training dataset. It is designed to load data batches onto the specified `device` (e.g., GPU).\n5. **Creating Validation Dataset and DataLoader**:\n   - Similar to the training dataset, `ds_val`, `ds_val_len`, `sampler_val`, `len_sampler_val`, and `dl_val` are created for the validation dataset. These components follow a similar pattern as the training dataset.\n6. **Creating DataLoaders**:\n   - `data` is an instance of `DataLoaders` that combines both the training and validation data loaders (`dl_train` and `dl_val`). It's used for training and evaluation.\n7. **Model Initialization**:\n   - `model = RNA_Model()`: An instance of the `RNA_Model` class is created which is a neural network model designed for RNA sequence analysis.\n8. **Moving Model to Device**:\n   - `model = model.to(device)`: The model is moved to the specified `device` (e.g., GPU) for training.\n9. **Creating a Learner**:\n   - `learn` is an instance of the `Learner` class, which is part of a deep learning framework, possibly fastai. It is configured with various settings, including loss function, metrics, and training callbacks.\n10. **Training**:\n    - `learn.fit_one_cycle(32, lr_max=5e-4, wd=0.05, pct_start=0.02)`: The training loop runs for 32 epochs (`fit_one_cycle` is a training method in fastai). It specifies the maximum learning rate, weight decay, and percentage of training iterations to increase the learning rate.\n11. **Saving Model Checkpoint**:\n    - `torch.save(learn.model.state_dict(), os.path.join(OUT, f'{fname}_{fold}.pth'))`: The trained model's state dictionary is saved to a file in the output directory. The file name appears to include the fold number.\n12. **Garbage Collection**:\n    - `gc.collect()`: This is used to perform garbage collection, potentially freeing up memory resources.\n    \nIt's important to note that this code snippet represents a single fold of a training loop. In practice, when dealing with cross-validation or multiple folds, you would typically iterate over different folds, train the model on each fold, and possibly aggregate the results from multiple folds for evaluation and model selection. Additionally, specific details about the `RNA_Model`, `RNA_Dataset`, and other custom classes used in this code would depend on their implementations, which are not provided here.\n\n[Return to Contents](#contents)","metadata":{}},{"cell_type":"markdown","source":"<a id =  \"RNA_sub\"></a>\n## RNA starter submission\n### (This is based on [RNA starter submission - 0.186 LB](https://www.kaggle.com/code/iafoss/rna-starter-submission-0-186-lb) by [Iafoss](https://www.kaggle.com/iafoss))","metadata":{}},{"cell_type":"code","source":"import pandas as pd\nimport os, gc\nimport numpy as np\nfrom tqdm.notebook import tqdm\nimport math\nfrom sklearn.model_selection import KFold\n\nimport torch\nimport torch.nn as nn\nimport torch.nn.functional as F\nfrom torch.utils.data import Dataset, DataLoader","metadata":{"execution":{"iopub.status.busy":"2023-10-30T15:20:25.748132Z","iopub.execute_input":"2023-10-30T15:20:25.748396Z","iopub.status.idle":"2023-10-30T15:20:30.206273Z","shell.execute_reply.started":"2023-10-30T15:20:25.748371Z","shell.execute_reply":"2023-10-30T15:20:30.205490Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Explanation**\n\n1. `from tqdm.notebook import tqdm`: Imports the `tqdm` module from the `tqdm.notebook` package. `tqdm` is a library for displaying progress bars in loops, making it easier to track the progress of time-consuming operations.\n2. `import math`: Imports the Python math module, which provides mathematical functions and constants.\n3. `from sklearn.model_selection import KFold`: Imports the `KFold` class from scikit-learn. `KFold` is a method for cross-validation, commonly used for splitting data into multiple subsets for training and testing machine learning models.\n4. `import torch`, `import torch.nn as nn`, `import torch.nn.functional as F`: Imports PyTorch and related modules.\n   - `torch` is the core PyTorch library.\n   - `torch.nn` is the module containing classes and functions for defining neural networks.\n   - `torch.nn.functional` provides various activation functions, loss functions, and other operations used in deep learning with PyTorch.\n5. `from torch.utils.data import Dataset, DataLoader`: Imports PyTorch's data handling utilities.\n   - `Dataset` is an abstract class used for creating custom datasets for training and evaluation.\n   - `DataLoader` is a class that loads data from a dataset into mini-batches, making it suitable for training deep learning models.","metadata":{}},{"cell_type":"code","source":"MODELS = ['/kaggle/input/rna-starter-jdm/example0_0.pth']\nPATH = '/kaggle/input/stanford-ribonanza-rna-folding-converted/'\nbs = 256\nnum_workers = 2\ndevice = 'cuda' if torch.cuda.is_available() else 'cpu'","metadata":{"execution":{"iopub.status.busy":"2023-10-30T15:20:42.068118Z","iopub.execute_input":"2023-10-30T15:20:42.068626Z","iopub.status.idle":"2023-10-30T15:20:42.101814Z","shell.execute_reply.started":"2023-10-30T15:20:42.068592Z","shell.execute_reply":"2023-10-30T15:20:42.100515Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Data","metadata":{}},{"cell_type":"code","source":"class RNA_Dataset_Test(Dataset):\n    def __init__(self, df, mask_only=False, **kwargs):\n        self.seq_map = {'A':0,'C':1,'G':2,'U':3}\n        df['L'] = df.sequence.apply(len)\n        self.Lmax = df['L'].max()\n        self.df = df\n        self.mask_only = mask_only\n        \n    def __len__(self):\n        return len(self.df)  \n    \n    def __getitem__(self, idx):\n        id_min, id_max, seq = self.df.loc[idx, ['id_min','id_max','sequence']]\n        mask = torch.zeros(self.Lmax, dtype=torch.bool)\n        L = len(seq)\n        mask[:L] = True\n        if self.mask_only: return {'mask':mask},{}\n        ids = np.arange(id_min,id_max+1)\n        \n        seq = [self.seq_map[s] for s in seq]\n        seq = np.array(seq)\n        seq = np.pad(seq,(0,self.Lmax-L))\n        ids = np.pad(ids,(0,self.Lmax-L), constant_values=-1)\n        \n        return {'seq':torch.from_numpy(seq), 'mask':mask}, \\\n               {'ids':ids}\n            \ndef dict_to(x, device='cuda'):\n    return {k:x[k].to(device) for k in x}\n\ndef to_device(x, device='cuda'):\n    return tuple(dict_to(e,device) for e in x)","metadata":{"execution":{"iopub.status.busy":"2023-10-30T15:20:50.288406Z","iopub.execute_input":"2023-10-30T15:20:50.289256Z","iopub.status.idle":"2023-10-30T15:20:50.303711Z","shell.execute_reply.started":"2023-10-30T15:20:50.289212Z","shell.execute_reply":"2023-10-30T15:20:50.302612Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Explanation**\n\nDefines a custom PyTorch dataset class called `RNA_Dataset_Test`.\n\n1. **Initialization (`__init__` method)**:\n   - `df`: This is a Pandas DataFrame containing the RNA sequence data. It's expected that the DataFrame has columns named `'id_min'`, `'id_max'`, and `'sequence'`.\n   - `mask_only`: A boolean parameter that, if set to `True`, indicates that the dataset should only return masks, without any sequence information.\n2. **`__len__` method**:\n   - Returns the length of the dataset, which is the number of rows in the DataFrame `df`.\n3. **`__getitem__` method**:\n   - This method is called to retrieve an item from the dataset at a given index `idx`.\n4. **Data Processing in `__getitem__`**:\n   - Extracts the `'id_min'`, `'id_max'`, and `'sequence'` values for the sample at index `idx` from the DataFrame `df`.\n   - Initializes a boolean mask tensor (`mask`) with a length of `self.Lmax`, which is the maximum sequence length in the dataset.\n   - Sets the first `L` elements of the mask to `True`, where `L` is the length of the sequence.\n   - If `self.mask_only` is `True`, it returns only the mask as input data.\n5. **Data Preparation if not Mask-Only**:\n   - If `self.mask_only` is `False`, it proceeds to prepare sequence data and additional information.\n   - It creates a NumPy array `ids` containing a range of values from `id_min` to `id_max`.\n   - It maps the characters in the sequence to integers using `self.seq_map` (which maps 'A' to 0, 'C' to 1, 'G' to 2, and 'U' to 3).\n   - Pads both the `seq` and `ids` arrays with zeros to match the length of `self.Lmax` if the sequence is shorter than `self.Lmax`.\n6. **Return Value**:\n   - Returns a dictionary containing input data and a dictionary containing additional information.\n   - For input data, it includes a `'seq'` tensor with the mapped sequence values and a `'mask'` tensor representing the mask.\n   - For additional information, it includes an `'ids'` tensor containing the IDs.\n\nAdditionally, there are two utility functions defined below the `RNA_Dataset_Test` class:\n\n- `dict_to(x, device='cuda')`: Converts a dictionary of tensors to the specified device (default is 'cuda' for GPU).\n- `to_device(x, device='cuda')`: Converts a tuple of dictionaries (typically input and target data) to the specified device.\n\nThis dataset class is useful for preparing RNA sequence data for testing or evaluation, allowing you to extract sequence information, apply masks, and map characters to numeric values. It's important to note that the implementation of `self.seq_map` is crucial for correctly mapping RNA sequence characters to integers.","metadata":{}},{"cell_type":"code","source":"class DeviceDataLoader:\n    def __init__(self, dataloader, device='cuda'):\n        self.dataloader = dataloader\n        self.device = device\n    \n    def __len__(self):\n        return len(self.dataloader)\n    \n    def __iter__(self):\n        for batch in self.dataloader:\n            yield tuple(dict_to(x, self.device) for x in batch)","metadata":{"execution":{"iopub.status.busy":"2023-10-30T15:21:18.072994Z","iopub.execute_input":"2023-10-30T15:21:18.073330Z","iopub.status.idle":"2023-10-30T15:21:18.079660Z","shell.execute_reply.started":"2023-10-30T15:21:18.073306Z","shell.execute_reply":"2023-10-30T15:21:18.078660Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The `DeviceDataLoader` class is a wrapper around a standard PyTorch `DataLoader` designed to load data batches onto a specified device (e.g., GPU). It facilitates moving data to the desired device during the iteration process.\n\n- **Initialization (`__init__` method)**:\n  - `dataloader`: This is the PyTorch `DataLoader` that you want to wrap.\n  - `device`: Specifies the target device (e.g., 'cuda' for GPU or 'cpu' for CPU). The default value is 'cuda', indicating GPU usage if available.\n- **`__len__` method**:\n  - Returns the length of the `DataLoader` object that is being wrapped. This allows you to use `len(DeviceDataLoader)` to get the number of batches in the original `DataLoader`.\n- **`__iter__` method**:\n  - This method allows the `DeviceDataLoader` to be used as an iterator, making it compatible with `for` loops.\n  - It iterates over batches from the original `DataLoader`.\n  - For each batch, it yields a tuple where each element of the batch (which is typically a dictionary of tensors) is converted to the specified device using the `dict_to` function. This ensures that all tensors in the batch are moved to the target device.\n\nThis `DeviceDataLoader` class provides a convenient way to iterate through data batches while automatically moving the data to the specified device, making it easier to train deep learning models on the GPU when needed. The `dict_to` function mentioned in this class is likely a utility function defined elsewhere in your code for converting a dictionary of tensors to the specified device.","metadata":{}},{"cell_type":"markdown","source":"### Model","metadata":{}},{"cell_type":"code","source":"class SinusoidalPosEmb(nn.Module):\n    def __init__(self, dim=16, M=10000):\n        super().__init__()\n        self.dim = dim\n        self.M = M\n\n    def forward(self, x):\n        device = x.device\n        half_dim = self.dim // 2\n        emb = math.log(self.M) / half_dim\n        emb = torch.exp(torch.arange(half_dim, device=device) * (-emb))\n        emb = x[...,None] * emb[None,...]\n        emb = torch.cat((emb.sin(), emb.cos()), dim=-1)\n        return emb","metadata":{"execution":{"iopub.status.busy":"2023-10-30T15:21:25.637554Z","iopub.execute_input":"2023-10-30T15:21:25.638300Z","iopub.status.idle":"2023-10-30T15:21:25.645099Z","shell.execute_reply.started":"2023-10-30T15:21:25.638267Z","shell.execute_reply":"2023-10-30T15:21:25.644010Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Explanation**\n\nThis class is a PyTorch module for generating sinusoidal positional embeddings. Sinusoidal positional embeddings are often used in sequence-related tasks, including natural language processing and time series analysis, to provide the model with information about the position or order of elements in a sequence.\n\n- **Initialization (`__init__` method)**:\n  - `dim`: The dimension of the positional embeddings. It determines how many different frequencies of sinusoids are used to represent position information. The default value is 16.\n  - `M`: A constant used to scale the positional embeddings. It defaults to 10,000, which is a common value used in practice.\n- **`forward` method**:\n  - This method computes the positional embeddings based on the input `x`.\n  - `device = x.device`: It retrieves the device (CPU or GPU) where the input tensor `x` resides.\n  - `half_dim = self.dim // 2`: It calculates half of the specified dimension, which will be used to create both sine and cosine components.\n  - `emb = math.log(self.M) / half_dim`: It computes a scaling factor for the sinusoidal embeddings.\n  - `emb = torch.exp(torch.arange(half_dim, device=device) * (-emb))`: It computes the exponential values for the sine and cosine components. These values are used to generate the periodic patterns.\n  - `emb = x[..., None] * emb[None, ...]`: It creates positional embeddings by multiplying the input `x` (which typically represents position indices) by the computed exponential values.\n  - `emb = torch.cat((emb.sin(), emb.cos()), dim=-1)`: It concatenates the sine and cosine components along the last dimension to form the final positional embeddings. This is a common practice in positional encoding.\n\nThe `forward` method takes an input tensor `x` representing positions and returns the corresponding sinusoidal positional embeddings. These embeddings are used to augment the input data in sequence-related tasks, providing the model with information about the positions of elements in the sequence. The use of sine and cosine components ensures that the embeddings have distinct and meaningful patterns for different positions.","metadata":{}},{"cell_type":"code","source":"class RNA_Model(nn.Module):\n    def __init__(self, dim=192, depth=12, head_size=32, **kwargs):\n        super().__init__()\n        self.emb = nn.Embedding(4,dim)\n        self.pos_enc = SinusoidalPosEmb(dim)\n        self.transformer = nn.TransformerEncoder(\n            nn.TransformerEncoderLayer(d_model=dim, nhead=dim//head_size, dim_feedforward=4*dim,\n                dropout=0.1, activation=nn.GELU(), batch_first=True, norm_first=True), depth)\n        self.proj_out = nn.Linear(dim,2)\n    \n    def forward(self, x0):\n        mask = x0['mask']\n        Lmax = mask.sum(-1).max()\n        mask = mask[:,:Lmax]\n        x = x0['seq'][:,:Lmax]\n        \n        pos = torch.arange(Lmax, device=x.device).unsqueeze(0)\n        pos = self.pos_enc(pos)\n        x = self.emb(x)\n        x = x + pos\n        \n        x = self.transformer(x, src_key_padding_mask=~mask)\n        x = self.proj_out(x)\n        \n        return x","metadata":{"execution":{"iopub.status.busy":"2023-10-30T15:21:37.050197Z","iopub.execute_input":"2023-10-30T15:21:37.051138Z","iopub.status.idle":"2023-10-30T15:21:37.062046Z","shell.execute_reply.started":"2023-10-30T15:21:37.051094Z","shell.execute_reply":"2023-10-30T15:21:37.060978Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Explanation**\n\nThe code defines a PyTorch neural network model called `RNA_Model`.\n\n1. **Initialization (`__init__` method)**:\n   - `dim`: This parameter sets the dimension of the embeddings and hidden layers in the model. It's set to 192 by default.\n   - `depth`: This parameter determines the number of layers in the Transformer encoder. It's set to 12 by default.\n   - `head_size`: This parameter controls the size of the attention heads in the Transformer encoder. It's calculated as `dim // head_size` by default.\n2. **Embedding Layer (`self.emb`)**:\n   - The model starts with an embedding layer for the input sequences. It embeds input tokens into a `dim`-dimensional space. In this case, it seems to be designed for sequences with four different tokens (possibly representing different nucleotides in RNA sequences).\n3. **Positional Encoding (`self.pos_enc`)**:\n   - This layer is used to add positional information to the embedded sequences. It uses the `SinusoidalPosEmb` module, which generates sinusoidal positional embeddings. These embeddings help the model understand the order of tokens in the input sequence.\n4. **Transformer Encoder (`self.transformer`)**:\n   - The core of the model is a stack of Transformer encoder layers (`nn.TransformerEncoder`). It uses `depth` layers of `nn.TransformerEncoderLayer`, each with the following characteristics:\n     - `d_model`: The dimension of model embeddings, which is set to `dim`.\n     - `nhead`: The number of attention heads. It's calculated as `dim // head_size`.\n     - `dim_feedforward`: The dimension of the feedforward network inside the Transformer layers, set to `4 * dim`.\n     - `dropout`: Dropout rate set to 0.1.\n     - `activation`: The activation function used within the model, which is the GELU (Gaussian Error Linear Unit) activation function.\n     - `batch_first` and `norm_first`: These parameters specify whether batch normalization should be applied before or after the attention and feedforward layers.\n5. **Output Projection Layer (`self.proj_out`)**:\n   - After the Transformer encoder, there's a linear projection layer that maps the hidden representations back to a 2-dimensional output. This suggests that the model is designed for a binary classification task where it needs to output two classes or values.\n6. **Forward Method (`forward`)**:\n   - The `forward` method defines how data flows through the model during inference.\n   - It takes an input dictionary `x0` which is expected to contain:\n     - `'mask'`: A mask for sequence padding, used to handle variable-length sequences.\n     - `'seq'`: The input sequences.\n   - It trims the mask and sequences to the maximum length found in the batch.\n   - Computes positional encodings and adds them to the embedded sequences.\n   - Passes the sequences through the Transformer encoder, using the mask to mask out padded elements.\n   - Applies the projection layer to obtain the final output, which is returned.\n\nOverall, this model is a variant of the Transformer architecture customized for a specific sequence classification task, likely involving RNA sequences. You can use this model for tasks like RNA sequence classification, and you can adjust the hyperparameters such as `dim`, `depth`, and `head_size` as needed.","metadata":{}},{"cell_type":"markdown","source":"### Inference","metadata":{}},{"cell_type":"code","source":"df_test = pd.read_parquet(os.path.join(PATH,'test_sequences.parquet'))\nds = RNA_Dataset_Test(df_test)\ndl = DeviceDataLoader(torch.utils.data.DataLoader(ds, batch_size=bs, \n               shuffle=False, drop_last=False, num_workers=num_workers), device)\ndel df_test\ngc.collect()\n\nmodels = []\nfor m in MODELS:\n    model = RNA_Model()   \n    model = model.to(device)\n    model.load_state_dict(torch.load(m,map_location=torch.device('cpu')))\n    model.eval()\n    models.append(model)","metadata":{"execution":{"iopub.status.busy":"2023-10-30T15:21:50.493446Z","iopub.execute_input":"2023-10-30T15:21:50.493823Z","iopub.status.idle":"2023-10-30T15:21:58.547059Z","shell.execute_reply.started":"2023-10-30T15:21:50.493794Z","shell.execute_reply":"2023-10-30T15:21:58.546042Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Explanation**\n\nThis is part of the inference pipeline for a machine learning model on RNA sequence data.\n\n1. **Loading Test Data**:\n   - `df_test = pd.read_parquet(os.path.join(PATH,'test_sequences.parquet'))`: Reads test RNA sequence data from a Parquet file into a Pandas DataFrame `df_test`.\n2. **Creating a Test Dataset and DataLoader**:\n   - `ds = RNA_Dataset_Test(df_test)`: Initializes a test dataset using the `RNA_Dataset_Test` class, likely designed for processing RNA sequence data.\n   - `dl = DeviceDataLoader(torch.utils.data.DataLoader(ds, batch_size=bs, shuffle=False, drop_last=False, num_workers=num_workers), device)`: Creates a test data loader using `DeviceDataLoader`. It wraps a PyTorch `DataLoader` for the test dataset, specifying the batch size, no shuffling (`shuffle=False`), not dropping the last incomplete batch (`drop_last=False`), and the number of workers for data loading (`num_workers`).\n3. **Memory Management and Cleanup**:\n   - `del df_test` and `gc.collect()`: Deletes the `df_test` DataFrame and triggers garbage collection to free up memory. This is a common practice to conserve memory after loading large datasets.\n4. **Loading Trained Models**:\n   - `models = []`: Initializes an empty list to store trained models.\n   - Iterates over a list of model paths specified by `MODELS`.\n   - For each model path (`m`), it:\n     - Initializes a new instance of the `RNA_Model`.\n     - Moves the model to the specified device (`device`) using `.to(device)`.\n     - Loads the model's trained weights from the file at path `m` using `torch.load`. The `map_location=torch.device('cpu')` argument ensures that the model is loaded onto the CPU, regardless of the original device it was trained on.\n     - Sets the model to evaluation mode using `.eval()`.\n     - Appends the loaded and prepared model to the `models` list.\n\nOverall, this code prepares the test dataset, sets up a data loader for batch processing, loads multiple trained models from specified paths, and ensures that the models are ready for inference on the specified device. It appears to be part of a larger pipeline for testing the models on RNA sequence data, likely for classification or prediction tasks.","metadata":{}},{"cell_type":"code","source":"ids,preds = [],[]\nfor x,y in tqdm(dl):\n    with torch.no_grad(),torch.cuda.amp.autocast():\n        p = torch.stack([torch.nan_to_num(model(x)) for model in models]\n                        ,0).mean(0).clip(0,1)\n        \n    for idx, mask, pi in zip(y['ids'].cpu(), x['mask'].cpu(), p.cpu()):\n        ids.append(idx[mask])\n        preds.append(pi[mask[:pi.shape[0]]])\n\nids = torch.concat(ids)\npreds = torch.concat(preds)\n\ndf = pd.DataFrame({'id':ids.numpy(), 'reactivity_DMS_MaP':preds[:,1].numpy(), \n                   'reactivity_2A3_MaP':preds[:,0].numpy()})\ndf.to_csv('submission.csv', index=False, float_format='%.4f') # 6.5GB\ndf.head()","metadata":{"execution":{"iopub.status.busy":"2023-10-30T15:22:11.045297Z","iopub.execute_input":"2023-10-30T15:22:11.045623Z","iopub.status.idle":"2023-10-30T15:47:35.509454Z","shell.execute_reply.started":"2023-10-30T15:22:11.045597Z","shell.execute_reply":"2023-10-30T15:47:35.507708Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Explanation**\n\nThis is the final step in the inference pipeline and **generates an output file of 19.5GB**\n\n1. **Initialization**:\n   - `ids` and `preds` are initialized as empty lists. These lists will be used to store IDs and corresponding predictions.\n2. **Iterating Through Data Loader (`dl`)**:\n   - The code iterates over batches of data using a `for` loop with `tqdm(dl)` for progress tracking.\n   - Inside the loop, it processes each batch of data:\n     - `x` contains the input data (likely sequences and masks).\n     - `y` contains additional information (IDs).\n   - The `torch.no_grad()` context manager is used to disable gradient computation during inference, which is more memory-efficient.\n3. **Inference with Models**:\n   - For each batch of data, it performs inference with multiple models:\n     - `torch.stack([torch.nan_to_num(model(x)) for model in models], 0)`: It iterates over the `models` list and applies each model to the input data `x`. The `torch.nan_to_num` function is used to replace NaN values in model predictions with zeros (for numerical stability).\n     - `.mean(0)`: Calculates the mean prediction across all models. This is often used in ensemble methods to combine predictions from multiple models.\n     - `.clip(0, 1)`: Clips the predictions to the range [0, 1] to ensure they are valid probabilities.\n4. **Collecting IDs and Predictions**:\n   - For each batch, it collects the IDs, masks, and predictions:\n     - `ids.append(idx[mask])`: Appends the IDs corresponding to the non-masked positions to the `ids` list.\n     - `preds.append(pi[mask[:pi.shape[0]]])`: Appends the predictions corresponding to the non-masked positions to the `preds` list. It uses the mask to exclude padded positions.\n5. **Concatenating and Processing Results**:\n   - After processing all batches, it concatenates the collected IDs and predictions:\n     - `ids = torch.concat(ids)`: Concatenates the lists of IDs into a single tensor.\n     - `preds = torch.concat(preds)`: Concatenates the lists of predictions into a single tensor.\n6. **Creating a DataFrame and Saving Results**:\n   - It creates a Pandas DataFrame `df` with columns for IDs and predictions.\n   - The DataFrame is then saved to a CSV file named `'submission.csv'`. The `float_format='%.4f'` argument specifies the format for saving floating-point numbers with four decimal places.\n7. **Displaying the First Rows of the DataFrame**:\n   - Finally, it displays the first few rows of the DataFrame to provide an overview of the saved results.\n\nThis code is responsible for running inference on the test data, collecting the predictions, and saving them in a submission file in CSV format. The saved CSV file can be used for evaluation and submission in a machine learning competition or task related to RNA sequence analysis.","metadata":{}},{"cell_type":"markdown","source":"[Return to Contents](#contents)","metadata":{}},{"cell_type":"markdown","source":"## Reference material\n\n### Competition code\n[RNA Science Computational Environment](https://www.kaggle.com/code/julianmacnamara/rna-science-computational-environment/edit)  \n[Starters EDA + Fastai + RNN](https://www.kaggle.com/code/julianmacnamara/starters-eda-fastai-rnn/edit)  \n[SR RNA reactivity](https://www.kaggle.com/code/julianmacnamara/sr-rna-reactivity-learn-eda-baseline/edit)  \n[Standford Ribonanza | RNA theory | short EDA](https://www.kaggle.com/code/julianmacnamara/standford-ribonanza-rna-theory-short-eda/edit)\n\n### Recurrent neural networks and sequence regression\n\n[RNN-from-scratch](https://www.kaggle.com/code/julianmacnamara/rnn-from-scratch/edit)  \n[Intro to Recurrent Neural Networks LSTM | GRU](https://www.kaggle.com/code/julianmacnamara/intro-to-recurrent-neural-networks-lstm-gru/edit)  \n[Recurrent Neural Network: Regression](https://www.kaggle.com/code/julianmacnamara/recurrent-neural-network-regression/edit)  \n[Recurrent Neural Network](https://www.kaggle.com/code/julianmacnamara/recurrent-neural-network/edit)  \n[Time Series Analysis using LSTM Keras](https://www.kaggle.com/code/julianmacnamara/time-series-analysis-using-lstm-keras/edit)  \n[Time Series Prediction with LSTM Recurrent Neural Networks in Python with Keras](https://machinelearningmastery.com/time-series-prediction-lstm-recurrent-neural-networks-python-keras/)","metadata":{}}]}