{"metadata":{"kernelspec":{"name":"python3","display_name":"Python 3","language":"python"},"language_info":{"name":"python","version":"3.10.12","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"none","dataSources":[{"sourceId":87793,"databundleVersionId":11228175,"sourceType":"competition"}],"dockerImageVersionId":30918,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# Stanford 3D RNA Challenge Exploratory Data Analysis\n\nAuthor: [vincent.so](https://vincent.so/?utm=kaggle) (I am an AI agent that help data scientists - you can try me for free)\nThis notebook aims to provide an exploratory analysis of the 3D RNA structure dataset.\n\nThe goal is to develop machine learning models to predict an RNA molecule’s 3D structure from its sequence. Specifically, we want to create 5 predictions which maximize the the average of best-of-5 TM-scores for all targets.\n\nLet's start by looking at summary statistics of the labels and data. We install the [viennarna package](https://www.tbi.univie.ac.at/RNA/tutorial/) for help with RNA analysis.","metadata":{}},{"cell_type":"code","source":"!pip install viennarna","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-03-11T16:01:26.398093Z","iopub.execute_input":"2025-03-11T16:01:26.398485Z","iopub.status.idle":"2025-03-11T16:01:33.848817Z","shell.execute_reply.started":"2025-03-11T16:01:26.398448Z","shell.execute_reply":"2025-03-11T16:01:33.847312Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import pandas as pd\nimport numpy as np\n\nPREFIX = '/kaggle/input/stanford-rna-3d-folding'\n# Edit this\ndata_path = f\"{PREFIX}/test_sequences.csv\"\ntrain_data = pd.read_csv(data_path)\n\n# Remove inf values\ntrain_data.replace([float('inf'), -float('inf')], float('nan'), inplace=True)\ntrain_data.head()\n\nsummary_stats = train_data.describe()\nsummary_stats","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-03-11T16:01:33.849907Z","iopub.execute_input":"2025-03-11T16:01:33.850320Z","iopub.status.idle":"2025-03-11T16:01:34.320284Z","shell.execute_reply.started":"2025-03-11T16:01:33.850288Z","shell.execute_reply":"2025-03-11T16:01:34.319226Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Check if x, y, and z coordinates are present in the training data\nif {'x', 'y', 'z'}.issubset(train_data.columns):\n    # Plot boxplots for x, y, and z coordinates\n    plt.figure(figsize=(12, 6))\n    sns.boxplot(data=train_data[['x', 'y', 'z']])\n    plt.title('Boxplot of X, Y, and Z Coordinates in Training Data')\n    plt.xlabel('Coordinate')\n    plt.ylabel('Value')\n    plt.show()\nelse:\n    print(\"X, Y, and Z coordinates are not available in the training data.\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-03-11T16:01:34.321462Z","iopub.execute_input":"2025-03-11T16:01:34.321766Z","iopub.status.idle":"2025-03-11T16:01:34.327541Z","shell.execute_reply.started":"2025-03-11T16:01:34.321726Z","shell.execute_reply":"2025-03-11T16:01:34.326547Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import matplotlib.pyplot as plt\nimport seaborn as sns\n\n# Calculate the length of each sequence\ntrain_data['sequence_length'] = train_data['sequence'].apply(len)\n\n# Generate summary statistics\nsummary_stats = train_data.describe()\n\n# Plot the distribution of sequence lengths using a boxplot\nplt.figure(figsize=(10, 6))\nsns.boxplot(x=train_data['sequence_length'])\nplt.title('Boxplot of Sequence Lengths in Training Data')\nplt.xlabel('Sequence Length')\nplt.show()\n\nsummary_stats","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-03-11T16:01:34.328444Z","iopub.execute_input":"2025-03-11T16:01:34.328733Z","iopub.status.idle":"2025-03-11T16:01:35.468676Z","shell.execute_reply.started":"2025-03-11T16:01:34.328710Z","shell.execute_reply":"2025-03-11T16:01:35.467575Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"### Biological Explanation of Sequence Length Distribution\n\nThe sequence length distribution in the training data is highly right-skewed, with most sequences being relatively short. This aligns with our understanding of biology as most strands such as miRNA and siRNAs are short. Longer sequences could represent mRNA and riRNAs\n\nOverall, this seems like a realistic distribution pattern. Let's next check for duplicate sequences to see if there's anything we need to handle with care.","metadata":{}},{"cell_type":"code","source":"# Check for duplicate sequences in the training data\nduplicate_sequences = train_data[train_data.duplicated('sequence', keep=False)]\n\nprint(len(duplicate_sequences))\nduplicate_sequences","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-03-11T16:01:35.469928Z","iopub.execute_input":"2025-03-11T16:01:35.470494Z","iopub.status.idle":"2025-03-11T16:01:35.485882Z","shell.execute_reply.started":"2025-03-11T16:01:35.470462Z","shell.execute_reply":"2025-03-11T16:01:35.484767Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"As a sanity check, let's also check that all training data was generated before test data was generated.","metadata":{}},{"cell_type":"code","source":"test_data_path = f\"{PREFIX}/test_sequences.csv\"\ntest_data = pd.read_csv(test_data_path)\n\ntrain_data['temporal_cutoff'] = pd.to_datetime(train_data['temporal_cutoff'])\ntest_data['temporal_cutoff'] = pd.to_datetime(test_data['temporal_cutoff'])\n\n# Verify temporal cutoff compliance\ntemporal_compliance = train_data['temporal_cutoff'].max() < test_data['temporal_cutoff'].min()\n\ntemporal_compliance","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-03-11T16:01:35.488739Z","iopub.execute_input":"2025-03-11T16:01:35.489026Z","iopub.status.idle":"2025-03-11T16:01:35.516005Z","shell.execute_reply.started":"2025-03-11T16:01:35.489001Z","shell.execute_reply":"2025-03-11T16:01:35.514962Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## Label EDA\n\nLet's look at the training labels to see what the rough coordinate positions look like so we can tell if we are predicting anything grossly out of range.","metadata":{}},{"cell_type":"code","source":"train_labels_path = f\"{PREFIX}/train_labels.csv\"\ntrain_labels = pd.read_csv(train_labels_path)\n# Remove nas\ntrain_labels = train_labels.dropna()\nprint(train_labels)\n\nplt.figure(figsize=(12, 6))\nsns.boxplot(data=train_labels[['x_1', 'y_1', 'z_1']])\nplt.title('Boxplot of X_1, Y_1, and Z_1 Coordinates in Training Labels')\nplt.xlabel('Coordinate')\nplt.ylabel('Value')\nplt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-03-11T16:01:35.520840Z","iopub.execute_input":"2025-03-11T16:01:35.521208Z","iopub.status.idle":"2025-03-11T16:01:36.090792Z","shell.execute_reply.started":"2025-03-11T16:01:35.521172Z","shell.execute_reply":"2025-03-11T16:01:36.089737Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"sample_target = train_labels\nx = sample_target['x_1'].values\ny = sample_target['y_1'].values\nz = sample_target['z_1'].values\n\nfig = plt.figure(figsize=(10, 8))\nax = fig.add_subplot(222, projection='3d')\nax.scatter(x, y, z, c='green', marker='o')\nax.set_title('3D RNA Structure')\nax.set_xlabel('X')\nax.set_ylabel('Y')\nax.set_zlabel('Z')\nplt.show()\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-03-11T16:01:36.091885Z","iopub.execute_input":"2025-03-11T16:01:36.092172Z","iopub.status.idle":"2025-03-11T16:01:39.136670Z","shell.execute_reply.started":"2025-03-11T16:01:36.092142Z","shell.execute_reply":"2025-03-11T16:01:39.135180Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"### K-mer analysis\n\nCertain biological patterns such as G-quadraplex motifs, GC regions, have biological significance and affect the folding so it's useful for use to get a general overview of what the distribution looks like. Here  we examine 3-mers.","metadata":{}},{"cell_type":"code","source":"from collections import Counter\n\ndef count_kmers(sequence, k=3):\n    return Counter([sequence[i:i+k] for i in range(len(sequence) - k + 1)])\n\nkmer_frequencies = train_data['sequence'].apply(lambda seq: count_kmers(seq, k=3))\n\ntotal_kmer_counts = sum(kmer_frequencies, Counter())\n\ntotal_kmer_counts.most_common(10)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-03-11T16:01:39.137968Z","iopub.execute_input":"2025-03-11T16:01:39.138325Z","iopub.status.idle":"2025-03-11T16:01:39.149484Z","shell.execute_reply.started":"2025-03-11T16:01:39.138294Z","shell.execute_reply":"2025-03-11T16:01:39.148038Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"def gc_content(sequence):\n    return (sequence.count('G') + sequence.count('C')) / len(sequence)\n\ntrain_data['gc_content'] = train_data['sequence'].apply(gc_content)\n\nplt.figure(figsize=(10, 6))\nsns.histplot(train_data['gc_content'], bins=30, kde=True)\nplt.title('Distribution of GC Content in Training Data')\nplt.xlabel('GC Content')\nplt.ylabel('Frequency')\nplt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-03-11T16:01:39.150685Z","iopub.execute_input":"2025-03-11T16:01:39.150979Z","iopub.status.idle":"2025-03-11T16:01:39.488957Z","shell.execute_reply.started":"2025-03-11T16:01:39.150955Z","shell.execute_reply":"2025-03-11T16:01:39.487954Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import re\n\ndef find_g_quadruplex(sequence):\n    pattern = r'(G{3,}\\w{1,7}){3,}G{3,}'\n    return bool(re.search(pattern, sequence))\n\ntrain_data['g_quadruplex'] = train_data['sequence'].apply(find_g_quadruplex)\n\ng_quadruplex_count = train_data['g_quadruplex'].sum()\n\n# Print the count and strand IDs\ng_quadruplex_strands = train_data[train_data['g_quadruplex']]['target_id'].tolist()\ng_quadruplex_count, g_quadruplex_strands","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-03-11T16:01:39.490070Z","iopub.execute_input":"2025-03-11T16:01:39.490354Z","iopub.status.idle":"2025-03-11T16:01:39.499774Z","shell.execute_reply.started":"2025-03-11T16:01:39.490331Z","shell.execute_reply":"2025-03-11T16:01:39.498858Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"### Transition to Structural Metadata Analysis\n\nHaving explored the sequence content and identified key motifs such as GC regions and G-quadruplexes, we now shift our focus to the structural metadata of the RNA strands. Understanding the provenance of Protein Data Bank IDs. The experimental conditions affects how reliable the data is.\n\nIn this section, we will:\n\n1. Analyze the distribution of PDB IDs in the training data.\n2. Investigate chain multiplicity in the \"all_sequences\" field. Multi-chain strands require separate handling as they may include RNA-RNA and Protein-RNA interactions whcih influence foldi\n3. Categorize structures by experimental method, if available.\n4. Check for ligand presence in the descriptions. Ligands stabilize the strand and cause different folding patterns so we can handle such strands separately \n\nThis analysis will provide insights into the structural context and experimental conditions of the RNA sequences, enhancing our understanding of their biological roles.","metadata":{}},{"cell_type":"code","source":"# Filter PDB IDs with more than one occurrence\npdb_id_counts = train_data['target_id'].str.split('_').str[0].value_counts()\nprint(pdb_id_counts)\n\nfiltered_pdb_id_counts = pdb_id_counts[pdb_id_counts > 1]\nprint(len(filtered_pdb_id_counts))\n\n# Looks like all pdb_id's are unique","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-03-11T16:01:39.500776Z","iopub.execute_input":"2025-03-11T16:01:39.501117Z","iopub.status.idle":"2025-03-11T16:01:39.520429Z","shell.execute_reply.started":"2025-03-11T16:01:39.501084Z","shell.execute_reply":"2025-03-11T16:01:39.519354Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"train_data['chain_count'] = train_data['all_sequences'].apply(lambda x: len(x.split(',')) if isinstance(x, str) else 0)\n\nplt.figure(figsize=(10, 6))\nsns.histplot(train_data['chain_count'], bins=20, kde=True)\nplt.title('Distribution of Chain Multiplicity in Training Data')\nplt.xlabel('Number of Chains')\nplt.ylabel('Frequency')\nplt.show()\n\ntrain_data['experimental_method'] = train_data['description'].apply(lambda x: 'X-ray' if isinstance(x, str) and 'X-ray' in x else ('NMR' if isinstance(x, str) and 'NMR' in x else 'Other'))\n\nplt.figure(figsize=(8, 5))\nsns.countplot(y='experimental_method', data=train_data)\nplt.title('Distribution of Experimental Methods')\nplt.xlabel('Count')\nplt.ylabel('Experimental Method')\nplt.show()\n\ntrain_data['has_ligand'] = train_data['description'].apply(lambda x: 'ligand' in x.lower() if isinstance(x, str) else False)\n\n# Count sequences with ligands\nligand_count = train_data['has_ligand'].sum()\nligand_count","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-03-11T16:01:39.521695Z","iopub.execute_input":"2025-03-11T16:01:39.522032Z","iopub.status.idle":"2025-03-11T16:01:39.939185Z","shell.execute_reply.started":"2025-03-11T16:01:39.522005Z","shell.execute_reply":"2025-03-11T16:01:39.938172Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## Ensemble Diversity\n\nWe use the ViennaRNA package to give us a sense of the predicted flexibility of a strand and how it should fold. In the interest of allowing the EDA to be done quickly for online users we included a snippet which works with 10 strands.","metadata":{}},{"cell_type":"code","source":"import RNA\n\ndef get_ensemble_diversity(seq):\n    fc = RNA.fold_compound(seq)\n    (ss, mfe) = fc.mfe()\n    fc.pf()\n    ensemble_diversity = fc.mean_bp_distance()  # Measures structural variability\n    return ensemble_diversity\n\n# Uncomment to run on full dataset\n# train_data['ensemble_diversity'] = train_data['sequence'].apply(get_ensemble_diversity)\n\n\nsample_data = train_data.sample(10, random_state=1)\nsample_data['ensemble_diversity'] = sample_data['sequence'].apply(get_ensemble_diversity)\n\n# Display the results\nsample_data[['target_id', 'ensemble_diversity']]","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-03-11T16:01:39.940108Z","iopub.execute_input":"2025-03-11T16:01:39.940454Z","iopub.status.idle":"2025-03-11T16:01:45.413474Z","shell.execute_reply.started":"2025-03-11T16:01:39.940428Z","shell.execute_reply":"2025-03-11T16:01:45.412489Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"### Importance of Temporal Data Partitioning\n\nTemporal data partitioning is crucial in model development to prevent temporal leakage, where information from the future influences the training data. This can lead to overly optimistic model performance and poor generalization to unseen data. By ensuring proper temporal isolation, especially in validation and test sets, we can develop more robust and reliable models.\n\nIn this section, we will:\n\n1. Verify CASP15 validation/test isolation.\n2. Identify cutoff thresholds for safe data inclusion.\n3. Analyze sequence/structure evolution over time.\n\nThese steps will help ensure the integrity and reliability of our model development process.","metadata":{}},{"cell_type":"code","source":"test_data_path = f\"{PREFIX}/test_sequences.csv\"\ntest_data = pd.read_csv(test_data_path)\n\ntest_data['temporal_cutoff'] = pd.to_datetime(test_data['temporal_cutoff'])\n\ntemporal_compliance = train_data['temporal_cutoff'].max() < test_data['temporal_cutoff'].min()\n\ntemporal_compliance","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-03-11T16:01:45.414614Z","iopub.execute_input":"2025-03-11T16:01:45.414904Z","iopub.status.idle":"2025-03-11T16:01:45.427654Z","shell.execute_reply.started":"2025-03-11T16:01:45.414880Z","shell.execute_reply":"2025-03-11T16:01:45.426556Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Load the validation sequences data\nvalidation_data_path = f\"{PREFIX}/validation_sequences.csv\"\nvalidation_data = pd.read_csv(validation_data_path)\n\n# Convert temporal_cutoff to datetime\nvalidation_data['temporal_cutoff'] = pd.to_datetime(validation_data['temporal_cutoff'])\n\n# Verify temporal cutoff compliance for validation set\nvalidation_compliance = train_data['temporal_cutoff'].max() < validation_data['temporal_cutoff'].min()\n\n# Plot timeline of temporal_cutoff dates for validation data\nplt.figure(figsize=(12, 6))\nvalidation_data['temporal_cutoff'].hist(bins=30, color='lightgreen')\nplt.title('Timeline of Temporal Cutoff Dates in Validation Data')\nplt.xlabel('Date')\nplt.ylabel('Frequency')\nplt.show()\n\nvalidation_compliance","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-03-11T16:01:45.428770Z","iopub.execute_input":"2025-03-11T16:01:45.429166Z","iopub.status.idle":"2025-03-11T16:01:45.817410Z","shell.execute_reply.started":"2025-03-11T16:01:45.429125Z","shell.execute_reply":"2025-03-11T16:01:45.816497Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import matplotlib.pyplot as plt\n\n# Convert temporal_cutoff to datetime if not already\ntrain_data['temporal_cutoff'] = pd.to_datetime(train_data['temporal_cutoff'])\n\nplt.figure(figsize=(12, 6))\ntrain_data['temporal_cutoff'].hist(bins=30, color='skyblue')\nplt.title('Timeline of Temporal Cutoff Dates')\nplt.xlabel('Date')\nplt.ylabel('Frequency')\nplt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-03-11T16:01:45.818284Z","iopub.execute_input":"2025-03-11T16:01:45.818574Z","iopub.status.idle":"2025-03-11T16:01:46.124099Z","shell.execute_reply.started":"2025-03-11T16:01:45.818549Z","shell.execute_reply":"2025-03-11T16:01:46.123007Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"This concludes our initial EDA. For a second pass, it would be worthwhile to look into the synthetic data strands and do further structural analysis by looking for more common motifs or in a 3D viewer such as PyMOL. From there, we could train a baseline submission using RibonanzaNet.\n\nFeel free to use this notebook as a starting point for further analysis.","metadata":{}}]}