{"cells": [{"cell_type": "markdown", "metadata": {}, "source": "# Introduction\nNormally, I start with an exploratory data analysis (EDA) script and dive straight into the data. However, this time the topic is more complex, requiring a deeper understanding before analysis is possible. To build a solid foundation, let\u2019s start with the fundamental question:\n\n**What is RNA? Why is it important? And why does its 3D structure matter?**\n\n---\n\n## What is RNA?\nRNA (Ribonucleic Acid) is a **nucleic acid**, a type of biological macromolecule, similar to DNA but with key differences in structure and function. It plays a crucial role in genetics and cellular biology.\n\n### Differences between RNA and DNA\n| Feature          | DNA                      | RNA                      |\n|-----------------|--------------------------|--------------------------|\n| Sugar           | Deoxyribose               | Ribose                   |\n| Strands         | Double-stranded (mostly, forming a **double helix**)  | Single-stranded (mostly) |\n| Bases          | A, T, C, G                 | A, **U**, C, G (Uracil replaces Thymine) |\n| Stability       | More stable               | Less stable (prone to degradation) |\n| Function        | **Genetic storage, replication** (stores hereditary information in cells) | **Protein synthesis, regulation, catalysis** (actively involved in gene expression) |\n\n### **DNA as the Blueprint of Life**\nDNA is **the genetic material in all known living organisms**. It carries the instructions for the growth, development, functioning, and reproduction of cells. The **double-helix structure** ensures that genetic information is stored safely and can be accurately copied during **cell division** (mitosis and meiosis). Every cell in an organism contains the same DNA, which acts as the \"blueprint\" for making proteins that define cellular function.\n\n---\n\n## **Types of RNA**\nRNA comes in different forms, each with distinct functions that work together to **convert genetic information into proteins**. The three main types of RNA directly involved in **protein synthesis** are:\n\n### **1. Messenger RNA (mRNA) \u2013 The Genetic Courier**\n- **Function:** Carries genetic instructions from DNA to ribosomes to synthesize proteins. It is a temporary copy of a gene used to direct protein synthesis.\n- **Example:** **mRNA vaccines** (e.g., Pfizer-BioNTech COVID-19 vaccine) work by introducing a synthetic mRNA sequence that cells use to produce a viral protein, training the immune system to recognize and fight the virus.\n\n### **2. Transfer RNA (tRNA) \u2013 The Adapter Molecule**\n- **Function:** Brings amino acids to the ribosome, matching them with the corresponding mRNA codon.\n- **Structure:** It has a **cloverleaf** shape and an **anticodon region** that pairs with mRNA codons.\n- **Example:** If mRNA has the codon `AUG` (Methionine), tRNA carries the corresponding amino acid (Methionine) and adds it to the growing protein chain.\n\n### **3. Ribosomal RNA (rRNA) \u2013 The Protein Factory**\n- **Function:** Forms the core of ribosomes and catalyzes peptide bond formation between amino acids.\n- **Highly conserved:** Used in evolutionary studies to compare species.\n- **Example:** The ribosome consists of both rRNA and proteins, acting as the **assembly site for proteins**.\n\n### **4. Small Nuclear RNA (snRNA) \u2013 The Gene Editor**\n- **Function:** Involved in **RNA processing**, especially **splicing**, which removes non-coding sequences (**introns**) from pre-mRNA before it becomes mature mRNA.\n- **Example:** Helps form the **spliceosome**, a molecular machine that cuts out introns and joins exons together, allowing different proteins to be produced from the same gene.\n\n### **5. MicroRNA (miRNA) & Small Interfering RNA (siRNA) \u2013 The Gene Silencers**\n- **Function:** Regulate gene expression by binding to mRNA and preventing translation or promoting degradation.\n- **Example:** miRNAs are involved in **cancer suppression** by regulating oncogenes (cancer-causing genes). siRNA is used in **gene therapy** to silence disease-causing genes.\n\n---\n\n## **Why is RNA important?**\nRNA is fundamental in biology and medicine because it acts as a bridge between genetic information (DNA) and functional molecules (proteins).\n\n### **1. Role in Protein Synthesis**\nRNA is central to the **central dogma of molecular biology**:\n\n**DNA \u2192 RNA \u2192 Protein**\n\n- **Transcription:** DNA is transcribed into mRNA in the nucleus.\n- **Translation:** mRNA is translated into proteins by ribosomes (rRNA + tRNA) in the cytoplasm.\n- **Example of the process:**\n  1. DNA: `TAC GGC TAA`\n  2. mRNA: `AUG CCG AUU`\n  3. tRNA brings the corresponding amino acids: Methionine - Proline - Isoleucine\n  4. The ribosome links these amino acids together to form a protein.\n\n### **2. Gene Regulation & Cell Function**\nRNA is not just a passive messenger; it actively regulates gene expression.\n- **miRNA and siRNA** silence genes by binding to mRNA and blocking translation.\n- **Long non-coding RNA (lncRNA)** influences chromatin structure and gene expression.\n- **Riboswitches** change their conformation to regulate gene activity based on environmental conditions.\n\n### **3. Importance in Evolution & Medicine**\n- RNA was likely the **first genetic material** in early life forms (**RNA world hypothesis**).\n- Many viruses, including **SARS-CoV-2**, use RNA as their genetic material.\n- **RNA-based drugs** and **RNA vaccines** are revolutionizing medicine.\n\n---\n## **What is 3D RNA folding?**\nUnlike DNA, RNA does not stay linear. It **folds into complex structures** that determine its function.\n\n### **Levels of RNA Structure**\n1. **Primary Structure**: The sequence of nucleotides (A, U, C, G).\n2. **Secondary Structure**: Local folding patterns (stems, loops, bulges).\n3. **Tertiary Structure**: The full **3D shape** stabilized by interactions.\n\n### **Key Folding Elements**\n- **Hairpin loops**: Short double-stranded stems with a loop at the end.\n- **Pseudoknots**: Complex folds where loops interact with other RNA regions.\n- **Coaxial stacking**: Stacking of helices for stability.\n\nRNA folding is **not random**\u2014it follows **thermodynamic principles** to find the most stable conformation.\n\n---\n\n## **Why do we need to study 3D RNA folding?**\nThe **3D structure determines RNA function**, affecting its ability to interact with other molecules.\n\n### **1. RNA-Protein Interactions**\n- Ribosomes, CRISPR, and viral replication depend on RNA-protein binding.\n\n### **2. Drug Design & Synthetic Biology**\n- Many drugs target RNA structures (e.g., antibiotics targeting rRNA).\n- Artificially designed RNA molecules can be used in therapy.\n\n---\n\n## **How does RNA 3D folding work?**\nRNA folds due to various **molecular forces** and stabilizing factors.\n\n### **1. Key Forces in RNA Folding**\n- **Hydrogen bonding**: Base pairing (A-U, G-C).\n- **Stacking interactions**: \u03c0-\u03c0 interactions between bases.\n- **Metal ion stabilization**: Mg\u00b2\u207a ions stabilize tertiary structures.\n\n### **2. Experimental Methods to Study RNA Folding**\n- **X-ray crystallography**: High-resolution but difficult for RNA.\n- **NMR spectroscopy**: Good for small RNAs.\n- **Cryo-EM**: Ideal for large RNA complexes like ribosomes.\n\n### **3. Computational Approaches**\n- **RNAfold**: Predicts 2D structures based on thermodynamics.\n- **Rosetta & AlphaFold RNA**: AI-driven 3D structure prediction.\n- **Molecular dynamics simulations**: Simulates RNA folding over time.\n\n---\n\n## **Challenges & Future Directions**\nDespite progress, RNA structure prediction remains a challenge.\n\n### **1. Predicting RNA Structures is Difficult**\n- Experimental methods are slow and expensive.\n- Computational models still struggle with large, complex RNAs.\n\n### **2. Combining Experimental & Computational Approaches**\n- Integrating lab data (e.g., Cryo-EM) with AI models improves accuracy.\n\n### **3. AI & Deep Learning in RNA Folding**\n- **AlphaFold RNA** is emerging as a powerful tool.\n- Machine learning could revolutionize RNA-based drug discovery.\n\n---\n### **Modern Approaches in RNA Structure Prediction**\nPredicting RNA 3D structures is a complex challenge that combines experimental and computational techniques. In recent years, **AI-driven models** have significantly improved RNA structure prediction.\n\n- **RNAfold**: A widely used thermodynamic model that predicts RNA secondary structures based on minimum free energy.\n- **Rosetta**: A computational tool that models RNA tertiary structures using fragment-based assembly and molecular dynamics.\n- **AlphaFold RNA**: A deep learning-based approach, inspired by AlphaFold for proteins, aiming to predict RNA 3D structures with high accuracy.\n\nBy integrating **machine learning, molecular simulations, and experimental data**, these tools are paving the way for more precise RNA structural insights, benefiting **drug discovery, synthetic biology, and gene regulation studies**.\n\n### **Example: RNAfold in Drug Discovery**\nOne practical application of RNA structure prediction is in **drug discovery**. Many **antiviral drugs** target RNA structures to disrupt viral replication. For example, **RNAfold** has been used to predict the secondary structures of **SARS-CoV-2 RNA elements**, helping researchers identify potential drug-binding sites. This approach aids in the design of **RNA-targeted therapeutics**, which could be crucial for treating RNA viruses and genetic disorders.\n\n\n\n## **Final Thoughts**\nUnderstanding RNA and its 3D folding is crucial for biology, medicine, and technology. By combining experimental and computational approaches, researchers can unlock new possibilities in **drug design, genetic engineering, and synthetic biology**.\n\n"}, {"cell_type": "markdown", "metadata": {}, "source": "# EDA"}, {"cell_type": "code", "execution_count": null, "metadata": {"_cell_guid": "b1076dfc-b9ad-4769-8c92-a6c4dae69d19", "_uuid": "8f2839f25d086af736a60e9eeb907d3b93b6e0e5", "trusted": true}, "outputs": [], "source": "import numpy as np \nimport pandas as pd \nimport seaborn as sns\nimport matplotlib.pyplot as plt\nfrom IPython.display import display, Markdown\nimport scipy.stats as stats\nimport plotly.graph_objects as go\nfrom mpl_toolkits.mplot3d import Axes3D\nimport plotly.express as px\nfrom collections import Counter\nimport subprocess\nimport re\nimport networkx as nx\nimport torch\n\n\n# Input data files are available in the read-only \"../input/\" directory\n# For example, running this (by clicking run or pressing Shift+Enter) will list all files under the input directory\n\nimport os\nfor dirname, _, filenames in os.walk('/kaggle/input'):\n    for filename in filenames:\n        print(os.path.join(dirname, filename))\n\n# Decide between local or kaggle cloud storage         \nKAGGLE_ENV = 'kaggle' in os.listdir('/')\ndata_path = '/kaggle/input' if KAGGLE_ENV else '../kaggle/input'\n\nif KAGGLE_ENV:\n    !pip install torch_geometric torch_scatter torch_sparse torch_cluster torch_spline_conv -f https://data.pyg.org/whl/torch-2.0.0+cpu.html\n    !apt-get update && apt-get install -y vienna-rna\n\nfrom torch_geometric.data import Data\n    \nfor dirname, _, filenames in os.walk(data_path):\n    for filename in filenames:\n        print(os.path.join(dirname, filename)) "}, {"cell_type": "markdown", "metadata": {}, "source": "## Load Data"}, {"cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": "file_paths = {\n    \"df_sample_submission_sunny\": data_path + \"/standford-rna-3d-folding-sunny/sample_submission.csv\",\n    \"df_test_sequences_sunny\": data_path + \"/standford-rna-3d-folding-sunny/test_sequences.csv\",\n    \"df_train_labels_sunny\": data_path + \"/standford-rna-3d-folding-sunny/train_labels.csv\",\n    \"df_train_sequences_sunny\": data_path + \"/standford-rna-3d-folding-sunny/train_sequences.csv\",\n    \"df_validation_labels\": data_path + \"/stanford-rna-3d-folding/validation_labels.csv\",\n    \"df_sample_submission\": data_path + \"/stanford-rna-3d-folding/sample_submission.csv\",\n    \"df_test_sequences\": data_path + \"/stanford-rna-3d-folding/test_sequences.csv\",\n    \"df_validation_sequences\": data_path + \"/stanford-rna-3d-folding/validation_sequences.csv\",\n    \"df_train_labels\": data_path + \"/stanford-rna-3d-folding/train_labels.csv\",\n    \"df_train_sequences\": data_path + \"/stanford-rna-3d-folding/train_sequences.csv\"\n}\n\nfor var_name, path in file_paths.items():\n    try:\n        globals()[var_name] = pd.read_csv(path)\n        print(f\"{var_name} load, {globals()[var_name].shape[0]} rows, {globals()[var_name].shape[1]} columns.\")\n    except FileNotFoundError:\n        print(f\"file not found: {path}\")\n    except Exception as e:\n        print(f\"error .. {path}: {e}\")"}, {"cell_type": "markdown", "metadata": {}, "source": "## df_train (RNA sequences)\nStart with the RNA sequences (df_train) and have a look about the features.\nLet's start with the first row:\n- target_id , 1SCL_A, the name of the RNA is 1SCL and the letter A is about the chain (RNA can have multiple chains). This RNA just have one Chain, https://www.rcsb.org/structure/1SCL.\n- sequence, An RNA sequence is the linear order of nucleotides in an RNA strand. It consists of the bases AUCG. This come first and the 3D folding happens later!\n- temporal cutoff - publish/ discovered date\n- description: This is the scientific description from the Protein Data Bank (PDB) (https://www.rcsb.org/). Each description includes the RNA type, biological function, specific organism, and experimental method used to determine the structure."}, {"cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": "display('df_train_sequences')\ndisplay(df_train_sequences.head())"}, {"cell_type": "markdown", "metadata": {}, "source": "## df_train_labels (3D RNA Structure)\n- **ID**: Unique identifier for each nucleotide in the RNA sequence.  \n  - **Focus on 1SCL**: The RNA structure **\"1SCL\"** consists of a **single chain (Chain A) with 29 nucleotides**. The sequence length defines the numbering **1\u201329**.  \n- **Resname**: The **Residue Name**, representing the type of nucleotide (**A, G, U, C**).  \n- **x1, y1, z1**: **3D Cartesian coordinates** of the **C1' (C1-prime) carbon atom**. These coordinates define the **spatial position** of the nucleotide in the 3D RNA structure.  \n  - **Focus on C1'**: The **C1' carbon atom serves as the central reference point** for each nucleotide in structural studies.  \n\n\n"}, {"cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": "display('df_train_labels')\ndisplay(df_train_labels.head(50))"}, {"cell_type": "markdown", "metadata": {}, "source": "## The Rest of the Data\nThe **validation datasets** (`df_validation_sequence` and `df_validation_labels`) have the **same structure** as the training dataframes.  \n\nThe **test dataset** (`df_test_sequence`) is identical to `df_sequence`, but the **3D structure labels** need to be predicted.\n"}, {"cell_type": "markdown", "metadata": {}, "source": "## Distribution (Boxplot) of the lenght of the RNA sequences"}, {"cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": "df_train_sequences[\"length\"] = df_train_sequences[\"sequence\"].str.len()\n\n\nmedian_value = df_train_sequences[\"length\"].median()\nq1 = df_train_sequences[\"length\"].quantile(0.25)  \nq3 = df_train_sequences[\"length\"].quantile(0.75)  \niqr = q3 - q1 \n\nlower_whisker = max(df_train_sequences[\"length\"].min(), q1 - 1.5 * iqr)\nupper_whisker = min(df_train_sequences[\"length\"].max(), q3 + 1.5 * iqr)\n\nmin_value = df_train_sequences[\"length\"].min()\nmax_value = df_train_sequences[\"length\"].max()\n\nstats_df = pd.DataFrame({\n    \"Statistic\": [\"Min\", \"Lower Whisker\", \"Q1 (25%)\", \"Median (50%)\", \"Q3 (75%)\", \"Upper Whisker\", \"Max\"],\n    \"Value\": [min_value, lower_whisker, q1, median_value, q3, upper_whisker, max_value]\n})\n\nprint(stats_df)\nplt.figure(figsize=(10, 5))\nsns.boxplot(x=df_train_sequences[\"length\"])\nplt.title(\"Distribution of the Length of RNA Sequences\")\nplt.xlabel(\"Sequence Length\")\nplt.show()"}, {"cell_type": "markdown", "metadata": {}, "source": "## Visualize the 3D RNA Structure"}, {"cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [], "source": "# If this not work at kaggle, please pull it to your local machine and run it there\ndf_1SCL = df_train_labels[df_train_labels[\"ID\"].str.startswith(\"1SCL_A_\")].iloc[:29]\n\nfig = go.Figure()\n\nfig.add_trace(go.Scatter3d(\n    x=df_1SCL[\"x_1\"], \n    y=df_1SCL[\"y_1\"], \n    z=df_1SCL[\"z_1\"],\n    mode='markers',\n    marker=dict(size=6, color='blue', opacity=0.8),\n    text=df_1SCL[\"resname\"],  # show by hovering over point the niklotide type\n    name=\"C1' Carbon Atoms\"\n))\n\nfig.update_layout(\n    title=\"Interactive 3D RNA Structure - 1SCL\",\n    scene=dict(\n        xaxis_title=\"X Coordinate\",\n        yaxis_title=\"Y Coordinate\",\n        zaxis_title=\"Z Coordinate\"\n    )\n)\nfig.show()"}, {"cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [], "source": "# Alternative 3D-Scatterplot for 1SCL RNA-Structure\ndf_1SCL = df_train_labels[df_train_labels[\"ID\"].str.startswith(\"1SCL_A_\")].iloc[:29]\n\nfig = plt.figure(figsize=(8, 6))\nax = fig.add_subplot(111, projection='3d')\n\nax.scatter(df_1SCL[\"x_1\"], df_1SCL[\"y_1\"], df_1SCL[\"z_1\"], c=\"blue\", marker=\"o\", label=\"C1' Carbon Atoms\")\n\nax.set_xlabel(\"X Coordinate\")\nax.set_ylabel(\"Y Coordinate\")\nax.set_zlabel(\"Z Coordinate\")\nax.set_title(\"3D RNA Structure - 1SCL (C1' Carbon Atoms)\")\n\nplt.legend()\nplt.show()"}, {"cell_type": "markdown", "metadata": {}, "source": "## Visualize all C1' Atoms (get a picture and look for outliers)\nYou can see there is one big cloud in the middle and a lot of many small clouds around it. But if you look carefully you can see there is one bigger cloud an the edge.. I got curious and began to analyze.. This is (I hope I'm right - the RNA 4V4G).https://www.rcsb.org/structure/4V4G\n<p align=\"center\">\n  <img src=\"https://cdn.rcsb.org/images/structures/4v4g_assembly-1.jpeg\" alt=\"Proteinstruct\" width=\"400\">\n</p>\n\n\n"}, {"cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [], "source": "# only in X-Y plane\nplt.figure(figsize=(8, 6))\nsns.scatterplot(x=df_train_labels[\"x_1\"], y=df_train_labels[\"y_1\"])\nplt.title(\"Scatterplot of RNA C1' Atoms (X-Y Plane)\")\nplt.xlabel(\"X Coordinate\")\nplt.ylabel(\"Y Coordinate\")\nplt.show()\n"}, {"cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [], "source": "# Visualization of RNA C1' Atoms in 3D\nif {\"x_1\", \"y_1\", \"z_1\"}.issubset(df_train_labels.columns):\n    fig = plt.figure(figsize=(10, 8))\n    ax = fig.add_subplot(111, projection='3d')\n\n    ax.scatter(df_train_labels[\"x_1\"], df_train_labels[\"y_1\"], df_train_labels[\"z_1\"], alpha=0.6)\n\n    ax.set_title(\"3D Scatterplot of RNA C1' Atoms\")\n    ax.set_xlabel(\"X Coordinate\")\n    ax.set_ylabel(\"Y Coordinate\")\n    ax.set_zlabel(\"Z Coordinate\")\n\n    plt.show()\n\nelse:\n    print(\"Missing columns: The columns 'x_1', 'y_1', and 'z_1' were not found.\")"}, {"cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [], "source": "# Visualize RNA C1' Atoms in 3D with Plotly (interactive), if it is not working, please run it on your local machine\nfig = px.scatter_3d(df_train_labels, x=\"x_1\", y=\"y_1\", z=\"z_1\", opacity=0.7)\nfig.show()"}, {"cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [], "source": "def find_atoms_in_region(df, x_target, y_target, z_target, tolerance=100):\n    \"\"\"\n    Filter all RNA atoms that are within a given tolerance range around a target coordinate.\n    \n    Parameters:\n    df (DataFrame): DataFrame containing atomic coordinates (must have 'x', 'y', 'z' columns).\n    x_target (float): X-coordinate of the region of interest.\n    y_target (float): Y-coordinate of the region of interest.\n    z_target (float): Z-coordinate of the region of interest.\n    tolerance (float): Range around the target coordinates to search for atoms.\n    \n    Returns:\n    DataFrame: Filtered DataFrame with atoms in the specified region.\n    \"\"\"\n    filtered_atoms = df[\n        (df[\"x_1\"].between(x_target - tolerance, x_target + tolerance)) &\n        (df[\"y_1\"].between(y_target - tolerance, y_target + tolerance)) &\n        (df[\"z_1\"].between(z_target - tolerance, z_target + tolerance))\n    ]\n    return filtered_atoms\n\nregion_atoms = find_atoms_in_region(df_train_labels, -730, -100, 400)\nprint(region_atoms) # there you can see how i found the 4V4G. Am I right?"}, {"cell_type": "markdown", "metadata": {}, "source": ""}, {"cell_type": "markdown", "metadata": {}, "source": "## K-mer Frequency Analysis in RNA Structure Prediction\nWhy look at K-mers?\n\nRNA folding is driven by sequence patterns. K-mers\u2014short sequence fragments\u2014help us spot motifs that influence structure.\nHow does this help our competition?\n- Better features for ML models\n- Identifies key motifs affecting folding\n- Makes predictions more interpretable\nBy analyzing K-mer frequencies, we capture useful patterns that could improve RNA 3D structure predictions. "}, {"cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [], "source": "# K-mer- Function, k=4\ndef get_kmers(sequence, k=3):\n    return [sequence[i:i+k] for i in range(len(sequence) - k + 1)]\n\nkmer_counts = Counter()\nfor seq in df_train_sequences[\"sequence\"]:\n    kmer_counts.update(get_kmers(seq, k=4))\nprint(kmer_counts.most_common(10))"}, {"cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [], "source": "# k=8\nkmer_counts = Counter()\nfor seq in df_train_sequences[\"sequence\"]:\n    kmer_counts.update(get_kmers(seq, k=8))\n\nprint(kmer_counts.most_common(10))"}, {"cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [], "source": "# k=12\nkmer_counts = Counter()\nfor seq in df_train_sequences[\"sequence\"]:\n    kmer_counts.update(get_kmers(seq, k=12))  \n\nprint(kmer_counts.most_common(10))"}, {"cell_type": "markdown", "metadata": {}, "source": "**How Do We Use K-mer Results in 3D RNA Folding?**\n\nK-mer analysis helps us identify patterns in RNA sequences that influence 3D structure formation. But what do we do with the results?\n\nFirst, we can use K-mer frequencies as features for machine learning models. By converting sequences into numerical representations, we can check if specific K-mer patterns correlate with stable or unstable 3D structures. This can improve the predictive power of our models.\n\nNext, K-mers can reveal structural motifs. Many RNA folding elements, such as hairpins or pseudoknots, are linked to specific sequence patterns. By analyzing frequent K-mers, we might detect biologically relevant motifs that contribute to RNA stability. Comparing them with predictions from RNAfold or 3D modeling tools can help validate their importance.\n\nAnother approach is clustering sequences based on K-mer distributions. If sequences with similar K-mer profiles tend to form similar 3D structures, this could provide valuable insights into RNA folding principles.\n\nFinally, if we have ground truth 3D structures, we can compare our K-mer results to see if specific patterns consistently appear in certain structural elements. This validation step ensures that our findings are meaningful rather than random.\nKey Takeaways:\n\n- K-mers as ML features: Improve structure prediction models.\n- Identify structural motifs: Detect sequence patterns linked to RNA folding.\n- Cluster similar sequences: Explore relationships between K-mers and 3D structures.\n- Validate with real data: Check if observed patterns align with known structures.\n\nBy leveraging K-mer analysis, we can uncover sequence-structure relationships that enhance RNA 3D folding predictions"}, {"cell_type": "markdown", "metadata": {}, "source": "## GC Content Analysis  \n\nGC content represents the proportion of **guanine (G) and cytosine (C) nucleotides** in an RNA sequence. It is an important factor in RNA stability because **G\u2261C base pairs have three hydrogen bonds**, making them stronger than **A=U base pairs (two hydrogen bonds)**.\n\n**Why is GC content important?**  \n- **RNA stability:** Higher GC content generally leads to more stable secondary structures.  \n- **Thermodynamics:** GC-rich sequences often have lower **free energy** (more stable folds).  \n- **RNA function:** Certain RNA types have characteristic GC content distributions.  \n\nBelow, we analyze the **distribution of GC content** in our dataset."}, {"cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [], "source": "# Calculate GC content\ndf_train_sequences[\"GC_content\"] = df_train_sequences[\"sequence\"].apply(\n    lambda seq: (seq.count(\"G\") + seq.count(\"C\")) / len(seq)\n)\n\nplt.figure(figsize=(10, 5))\nsns.histplot(df_train_sequences[\"GC_content\"], bins=30, kde=True)\nplt.title(\"GC Content Distribution\")\nplt.xlabel(\"GC Content\")\nplt.ylabel(\"Frequency\")\nplt.show()\n\nprint(df_train_sequences[\"GC_content\"].describe())"}, {"cell_type": "markdown", "metadata": {}, "source": "-Most sequences are GC-rich, which suggests increased RNA stability.\n- A few sequences are entirely A/U or entirely G/C, which might have special properties.\n- GC content varies but mostly stays within a normal biological range (40-70%)."}, {"cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [], "source": "# Merge sequence-based GC content with structural labels\ndf_combined = df_train_sequences.copy()\ndf_combined[\"ID\"] = df_train_labels[\"ID\"]  # Assuming IDs match\ndf_combined[\"resid\"] = df_train_labels[\"resid\"]\ndf_combined[\"x_1\"] = df_train_labels[\"x_1\"]\ndf_combined[\"y_1\"] = df_train_labels[\"y_1\"]\ndf_combined[\"z_1\"] = df_train_labels[\"z_1\"]\n\n# Compute distance of each nucleotide from the origin (as a rough spatial feature)\ndf_combined[\"distance_from_origin\"] = np.sqrt(\n    df_combined[\"x_1\"]**2 + df_combined[\"y_1\"]**2 + df_combined[\"z_1\"]**2\n)\n\n# Scatter plot of GC content vs. 3D structure distance\nplt.figure(figsize=(8, 5))\nsns.scatterplot(x=df_combined[\"GC_content\"], y=df_combined[\"distance_from_origin\"])\nplt.xlabel(\"GC Content\")\nplt.ylabel(\"Distance from Origin (3D Space)\")\nplt.title(\"GC Content vs. 3D RNA Spatial Structure\")\nplt.show()\n\n# Compute correlation\ncorrelation = df_combined[[\"GC_content\", \"distance_from_origin\"]].corr()\nprint(\"Correlation between GC content and 3D structure distance:\\n\", correlation)"}, {"cell_type": "markdown", "metadata": {}, "source": "- The correlation coefficient is -0.09868, which is very close to 0.\n- A weak negative correlation suggests that GC content has little to no direct impact on the 3D spatial positioning of nucleotides.\n\nRNA folding is driven by more than just GC content\n- While GC base pairs increase stability, tertiary interactions (e.g., pseudoknots, loops) determine final 3D shape.\n- Other forces (like protein interactions, metal ion binding) influence structure.\n\nNext step for that would be to check about GC content affects base-pairing density..."}, {"cell_type": "markdown", "metadata": {}, "source": "## Using RNA Secondary Structure for Machine Learning & Graph Neural Networks (GNNs)\nRNA is not just a linear sequence; it folds into complex secondary and tertiary structures that determine its function. The secondary structure, which can be predicted using RNAfold, serves as an intermediate step before the full 3D folding and provides critical information about base pairings, loops, and structural stability.\n\nThis section explores how RNA secondary structure can be utilized for Machine Learning and Graph Neural Networks (GNNs):\n- RNA as a Graph: Nucleotides as nodes, base pairings as edges.\n- Folding energy as a feature: Lower values indicate more stable structures!\n- RNA classification: Using structural features in ML models (XGBoost, Random Forest).\n- GNN Training: Predicting RNA structures using Graph Convolutional Networks (GCNs).\n\nThese approaches enable a deeper analysis of RNA sequences and help identify structural patterns, which is essential for applications such as RNA-based drug development and functional RNA classification."}, {"cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [], "source": "def predict_rna_structure(sequence):\n    \"\"\"Predicts RNA secondary structure and free energy using RNAfold.\"\"\"\n    process = subprocess.run(\n        [\"RNAfold\"], input=sequence, capture_output=True, text=True\n    )\n    output = process.stdout.strip().split(\"\\n\")\n    if len(output) < 2:\n        return None, None\n\n    structure = output[1].split(\" \")[0]  # Extract secondary structure\n    energy_match = re.search(r\"-?\\d+\\.\\d+\", output[1])  # Extract energy\n    energy = float(energy_match.group()) if energy_match else None\n\n    return structure, energy\n\ndef rna_to_graph(sequence, structure):\n    \"\"\"Converts RNA sequence and secondary structure into a graph representation.\"\"\"\n    G = nx.Graph()\n    for i, nucleotide in enumerate(sequence):\n        G.add_node(i, nucleotide=nucleotide)\n\n    # Sequential edges\n    for i in range(len(sequence) - 1):\n        G.add_edge(i, i + 1)\n\n    # Base-pairing edges\n    stack = []\n    for i, char in enumerate(structure):\n        if char == \"(\":\n            stack.append(i)\n        elif char == \")\":\n            if stack:\n                j = stack.pop()\n                G.add_edge(i, j)\n\n    return G\n\n# Get RNA sequence (1SCL example)\nsequence = df_train_sequences[\"sequence\"].iloc[0]\n\n# Predict secondary structure & energy\nstructure, energy = predict_rna_structure(sequence)\n\nif structure:\n    print(f\"RNA Structure: {structure} | Energy: {energy} kcal/mol\")\n\n    # Convert to graph and PyTorch Geometric format\n    rna_graph = rna_to_graph(sequence, structure)\n    edge_index = torch.tensor(list(rna_graph.edges), dtype=torch.long).t().contiguous()\n    x = torch.tensor([[ord(n)] for n in sequence], dtype=torch.float)  \n\n    data = Data(x=x, edge_index=edge_index, energy=torch.tensor([energy], dtype=torch.float))\n    print(data)\nelse:\n    print(\"Error: Could not generate RNA secondary structure.\")\n"}, {"cell_type": "markdown", "metadata": {}, "source": "The secondary structure is represented in dot-bracket notation, where:\n- ( and ) indicate base-paired nucleotides.\n- . represents unpaired nucleotides (loops or bulges\n\nThe energy value -9.7 kcal/mol represents the stability of the RNA fold:\n- More negative values indicate a more stable RNA structure.\n- Near-zero or positive values suggest an unstable or unfolded RNA.\n\nNodes and Edges\n- x=[29, 1]: There are 29 nodes, each represented by a single feature (e.g., ASCII encoding of A, U, G, C).\n- edge_index=[2, 37]: The graph contains 37 edges, indicating the number of nucleotide connections (both sequential and base-pairing).\n- energy=[1]: The energy feature is stored as a graph-level attribute, which can be used for stability prediction in ML models."}, {"cell_type": "code", "execution_count": 19, "metadata": {}, "outputs": [], "source": "print(f\"Total Nodes: {data.num_nodes}\")\nprint(f\"Total Edges: {data.num_edges}\")\nprint(\"Edge List:\")\nprint(data.edge_index.t().tolist())  # Print all edges to inspect them\n"}, {"cell_type": "markdown", "metadata": {}, "source": "- Total Nodes (29): Each nucleotide in the RNA sequence is a node in the graph.\n- Total Edges (37): These include both sequential connections (linking adjacent nucleotides) and base-pairing connections (from the dot-bracket notation)."}, {"cell_type": "markdown", "metadata": {}, "source": "- Each pair of numbers [i, j] represents an edge between nucleotide i and nucleotide j.\n- The graph contains:\n    - Sequential edges (e.g., [0,1], [1,2], [2,3]), which connect adjacent nucleotides.\n    - Base-pairing edges (e.g., [0,28], [1,27], [2,26]), derived from the dot-bracket structure.\n\nThis edge list forms the connectivity structure required for GNN training, allowing the model to learn from both local (sequential) and long-range (base-pairing) interactions."}, {"cell_type": "markdown", "metadata": {}, "source": "**Summary: Why Is This Useful?**\n- Captures RNA folding structure in a format suitable for Machine Learning.\n- Encodes long-range interactions, which are critical for RNA function.\n- Prepares RNA data for GNN-based classification and prediction models."}, {"cell_type": "markdown", "metadata": {}, "source": "# Classic EDA Overview"}, {"cell_type": "code", "execution_count": 20, "metadata": {}, "outputs": [], "source": "def data_overview(data, target):\n    # Overview\n    display(Markdown(\"## Data Overview\"))\n    \n    display(Markdown(\"### General Information\"))\n    display(Markdown(f\"- Number of rows and columns: {data.shape[0]} x {data.shape[1]}\"))\n    display(Markdown(\"- Column names:\"))\n    display(list(data.columns))\n\n    display(Markdown(\"### Data Types & Missing Values\"))\n    missing = data.isnull().sum()\n    dtypes = pd.DataFrame(data.dtypes, columns=[\"Data Type\"])\n    missing_df = pd.DataFrame(missing, columns=[\"Missing Values\"])\n    overview_df = dtypes.join(missing_df)\n    display(overview_df.style.background_gradient(cmap=\"coolwarm\"))\n\n    display(Markdown(\"### Classic head of Data\"))\n    display(data.head().style.set_properties(**{\"background-color\": \"#f5f5f5\"}))\n\n    display(Markdown(\"### Statistical Summary (describe)\"))\n    display(data.describe().T.style.background_gradient(cmap=\"viridis\"))\n\n    # Target variable analysis\n    if target is not None:\n        display(Markdown(f\"## Target Variable: `{target}`\"))\n        sns.set_style(\"whitegrid\")  \n        sns.set_palette(\"viridis\")   \n\n        fig, ax = plt.subplots(1, 2, figsize=(14, 5))\n\n        #   Absolute frequency barplot\n        sns.barplot(x=data[target].value_counts().index, \n                y=data[target].value_counts(), \n                ax=ax[0])  \n\n        ax[0].set_title(\"Absolute Frequency\", fontsize=12, fontweight=\"bold\")\n        ax[0].set_ylabel(\"Count\")\n        ax[0].set_xlabel(target)\n        ax[0].grid(axis=\"y\", linestyle=\"--\", alpha=0.5)  \n\n        # Percentage distribution barplot\n        sns.barplot(x=data[target].value_counts().index, \n                    y=data[target].value_counts(normalize=True), \n                    ax=ax[1])  \n\n        ax[1].set_title(\"Percentage Distribution\", fontsize=12, fontweight=\"bold\")\n        ax[1].set_ylabel(\"Percentage\")\n        ax[1].set_xlabel(target)\n        ax[1].grid(axis=\"y\", linestyle=\"--\", alpha=0.5)\n\n    \n\n        for spine in [\"top\", \"right\"]:\n            ax[0].spines[spine].set_visible(False)\n            ax[1].spines[spine].set_visible(False)\n\n    plt.tight_layout()\n    plt.show()"}, {"cell_type": "code", "execution_count": 21, "metadata": {}, "outputs": [], "source": "data_overview(df_train_sequences, target=None)"}, {"cell_type": "code", "execution_count": 22, "metadata": {}, "outputs": [], "source": "data_overview(df_train_sequences_sunny, target=None)"}, {"cell_type": "code", "execution_count": 23, "metadata": {}, "outputs": [], "source": "data_overview(df_train_labels, target=None)"}, {"cell_type": "code", "execution_count": 24, "metadata": {}, "outputs": [], "source": "data_overview(df_train_labels_sunny, target=None)"}, {"cell_type": "markdown", "metadata": {}, "source": "# Open Points\n- Check about the PCA to reduce the dimensionality of the RNA dataset.\n- ... Any Suggestions?"}], "metadata": {"kaggle": {"accelerator": "none", "dataSources": [{"databundleVersionId": 10008389, "sourceId": 84895, "sourceType": "competition"}], "isGpuEnabled": false, "isInternetEnabled": true, "language": "python", "sourceType": "notebook"}, "kernelspec": {"display_name": "base", "language": "python", "name": "python3"}, "language_info": {"codemirror_mode": {"name": "ipython", "version": 3}, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.11.5"}}, "nbformat": 4, "nbformat_minor": 4}