{"metadata":{"kernelspec":{"display_name":"Python 3","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.10.12"},"kaggle":{"accelerator":"nvidiaTeslaT4","dataSources":[{"sourceType":"competition","sourceId":87793,"databundleVersionId":11553390,"isSourceIdPinned":false}],"dockerImageVersionId":30919,"isInternetEnabled":false,"language":"python","sourceType":"notebook","isGpuEnabled":true},"papermill":{"default_parameters":{},"duration":144.373007,"end_time":"2025-03-26T03:44:56.636459","environment_variables":{},"exception":null,"input_path":"__notebook__.ipynb","output_path":"__notebook__.ipynb","parameters":{},"start_time":"2025-03-26T03:42:32.263452","version":"2.6.0"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# RNA 3D Structure Prediction Pipeline: Graph-Based Approach 🧬\n\n## Overview 📜\nThis project implements an advanced graph-based pipeline for the Stanford RNA 3D Folding competition, focusing on predicting the three-dimensional structure of RNA molecules from nucleotide sequences. The approach leverages graph representation of RNA molecules to explicitly model both local and long-range interactions critical for accurate structure prediction.\n\n## Key Components 🧩\n\n### 1. Data Processing and Management 📊\n**Memory-Optimized Data Loading**\n- Efficient loading for large RNA datasets\n- Coordinate normalization to handle numerical issues\n- Padding strategies for variable-length sequences\n\n**Feature Engineering**\n- One-hot encoding of RNA sequences (A, C, G, U, N)\n- Sequence composition analysis (GC/AU content)\n- Structure normalization and centralization\n\n### 2. Graph Model and Processing 🌐\n**RNA Graph Representation**\n- Explicit modeling of RNA as a graph structure\n- Nucleotides as nodes with comprehensive feature vectors\n- Multiple edge types (covalent, base-pairing, stacking, tertiary)\n- Weighted connections based on interaction confidence\n\n**Graph Neural Network Architecture**\n- Attention-based message passing between nucleotides\n- Type-specific edge processing for different interactions\n- Ability to model long-range tertiary interactions\n- Multi-head attention for capturing complex relationship patterns\n\n### 3. Advanced Structure Generation Pipeline 🧮\n**Contact Map Prediction**\n- Prediction of base-pairing interactions\n- Confidence scoring for each potential contact\n- Filtering of redundant or mutually exclusive contacts\n\n**Fragment-Based Assembly**\n- Decomposition of sequences into structural elements\n- Library of RNA structural motifs\n- Template-based assembly of structural elements\n- Stem-loop and other motif identification\n\n**3D Coordinate Generation**\n- Reference model initialization with optimized parameters\n- Integration of graph-based contact predictions\n- Fragment-based structure assembly\n- Multi-stage refinement process\n\n**Geometric Refinement**\n- Distance geometry optimization\n- Backbone dihedral adjustment\n- RNA-specific constraint satisfaction\n- Weighted coordinate averaging of multiple models\n\n**Ensemble Generation**\n- Progressive structural variation\n- Size-dependent perturbation scaling\n- Preservation of critical structural features\n- Generation of diverse yet physically plausible models\n\n### 4. Structure Generation and Sampling 🎯\n**Geometric Sampling with Physical Constraints**\n- Correlated noise for natural structural transitions\n- Bond length preservation (3.8 Å typical RNA backbone)\n- Global movement simulation for domain flexibility\n- Sequence-aware geometric constraints\n\n**RNA-Specific Structure Refinement**\n- Fragment-based structural assembly\n- Stem-loop template application\n- GC/AU content-based refinement\n- Backbone angle adjustment to RNA-specific values\n- Natural hinge point detection and rotation\n\n**Distance Geometry Optimization**\n- Explicit modeling of distance constraints\n- Gradient-based refinement of coordinates\n- Prioritization of important structural contacts\n- Weighted optimization for critical interactions\n\n### 5. Submission Creation ✅\n- Processing all test sequences through graph-based pipeline\n- Generation of 5 models per RNA sequence\n- Comprehensive logging and error handling\n- Performance statistics reporting\n- Multiple submission formats (graph-based, standard)\n\n### 6. Execution Mode 🔄\n**Graph-Based Pipeline Mode**\n- Advanced approach using graph representation of RNA\n- Fragment-based assembly with library of structural motifs\n- Distance geometry optimization for constraint satisfaction\n- Explicit modeling of tertiary interactions\n\n## Methodology 🔍\nThe graph-based pipeline employs a specialized approach focusing on RNA tertiary structure:\n\n1. **Data Preparation**: Sequences are converted to one-hot encoding and structures are normalized to ensure numerical stability.\n\n2. **Graph Construction**: RNA sequences are converted into graph representations with different edge types for various interactions:\n - Covalent edges for backbone connections\n - Base-pairing edges for Watson-Crick and wobble pairs\n - Stacking edges for adjacent nucleotides\n - Tertiary interaction edges for complex motifs\n\n3. **Structure Prediction**:\n - Contact map prediction identifies potential base pairs and tertiary interactions\n - Initial coordinates are generated using reference model\n - Fragment-based assembly creates structural templates\n - Coordinates are refined through distance geometry optimization\n - Ensemble of structures is generated with controlled variations\n\n4. **Size-Adaptive Strategy**: Different parameters are applied based on RNA size:\n - Small RNAs (<50 residues): Higher structural diversity with controlled variations\n - Medium RNAs (50-120 residues): Balanced approach with moderate variations\n - Large RNAs (>120 residues): Conservative variations with optimized parameters\n\n5. **Submission Creation**: Final structures are compiled into the required submission format with comprehensive error handling and fallback mechanisms.\n\nThe graph-based approach is particularly effective at capturing the long-range tertiary interactions that are critical for complex RNA folding, resulting in more accurate predictions especially for RNAs with complex 3D architectures.","metadata":{"papermill":{"duration":0.012176,"end_time":"2025-03-26T03:42:35.898136","exception":false,"start_time":"2025-03-26T03:42:35.885960","status":"completed"},"tags":[]}},{"cell_type":"markdown","source":"## Library Imports 📚🔧","metadata":{"papermill":{"duration":0.009894,"end_time":"2025-03-26T03:42:35.918527","exception":false,"start_time":"2025-03-26T03:42:35.908633","status":"completed"},"tags":[]}},{"cell_type":"code","source":"# Standard Library Imports\nimport datetime\nimport gc\nimport hashlib\nimport json\nimport os\nimport random\nimport time\nimport traceback\nimport warnings\nfrom collections import Counter\n\n# Scientific Computing and Numerical Libraries\nimport numpy as np\nimport pandas as pd\n\n# Graph Processing\nimport networkx as nx\n\n# Machine Learning and Deep Learning\nimport torch\nimport torch.nn as nn\n\n# Visualization Libraries\nimport matplotlib.colors as mcolors\nimport matplotlib.pyplot as plt\n\n# Machine Learning Library Import with Error Handling\ntry:\n # TensorFlow and Keras\n import tensorflow as tf\n from tensorflow.keras import layers, models, optimizers\n from tensorflow.keras.callbacks import EarlyStopping\n from tensorflow.keras.layers import (\n BatchNormalization, Bidirectional, Conv1D, \n Dense, Dropout, Flatten, Input, LSTM, Reshape\n )\n from tensorflow.keras.models import Model\n \n # Scikit-learn\n from sklearn.model_selection import train_test_split\n \n # XGBoost\n import xgboost as xgb\n \n ML_AVAILABLE = True\nexcept ImportError:\n print(\"Warning: ML libraries not available. Will use only reference-based methods.\")\n ML_AVAILABLE = False\n\n# Set random seed for reproducibility\nnp.random.seed(0)\n\n# Suppress warnings\nwarnings.filterwarnings('ignore')","metadata":{"execution":{"iopub.status.busy":"2025-04-09T15:13:12.907647Z","iopub.execute_input":"2025-04-09T15:13:12.907961Z","iopub.status.idle":"2025-04-09T15:13:28.413983Z","shell.execute_reply.started":"2025-04-09T15:13:12.907933Z","shell.execute_reply":"2025-04-09T15:13:28.413349Z"},"papermill":{"duration":21.560049,"end_time":"2025-03-26T03:42:57.488491","exception":false,"start_time":"2025-03-26T03:42:35.928442","status":"completed"},"tags":[],"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## 🧬 RNA 3D Structure Prediction and Analysis Pipeline 🔬","metadata":{"papermill":{"duration":0.009708,"end_time":"2025-03-26T03:42:57.508811","exception":false,"start_time":"2025-03-26T03:42:57.499103","status":"completed"},"tags":[]}},{"cell_type":"code","source":"# Directories and files adjusted for the new competition\nDATA_DIR = os.getenv('DATA_DIR', '/kaggle/input/stanford-rna-3d-folding/')\nmain_files = [\n \"train_sequences.csv\", \n \"train_labels.csv\", \n \"validation_sequences.csv\", \n \"validation_labels.csv\", \n \"test_sequences.csv\",\n \"sample_submission.csv\"\n]\n\nDEFAULT_THRESHOLD = 0.4 # Default threshold after analysis\n\ndef optimize_dataframe(df, inplace=False, category_threshold=DEFAULT_THRESHOLD):\n \"\"\"\n Optimizes the DataFrame to save memory.\n \"\"\"\n if category_threshold < 0 or category_threshold > 1:\n raise ValueError(\"category_threshold must be between 0 and 1.\")\n \n if not inplace:\n df = df.copy()\n \n for col in df.columns:\n col_type = df[col].dtype\n if np.issubdtype(col_type, np.integer):\n c_min, c_max = df[col].min(), df[col].max()\n if c_min > np.iinfo(np.int8).min and c_max < np.iinfo(np.int8).max:\n df[col] = df[col].astype(np.int8)\n elif c_min > np.iinfo(np.int16).min and c_max < np.iinfo(np.int16).max:\n df[col] = df[col].astype(np.int16)\n elif c_min > np.iinfo(np.int32).min and c_max < np.iinfo(np.int32).max:\n df[col] = df[col].astype(np.int32)\n elif np.issubdtype(col_type, np.floating):\n if df[col].min() > np.finfo(np.float32).min and df[col].max() < np.finfo(np.float32).max:\n df[col] = df[col].astype(np.float32)\n if col_type == object:\n unique_vals = len(df[col].unique())\n if unique_vals / len(df) < category_threshold:\n df[col] = df[col].astype('category')\n \n return df\n\ndef load_main_data(chunksize=50000):\n \"\"\"\n Loads the main files.\n \"\"\"\n data = {}\n for file_name in main_files:\n file_path = os.path.join(DATA_DIR, file_name)\n if os.path.exists(file_path):\n chunks = pd.read_csv(file_path, on_bad_lines='skip', low_memory=False, chunksize=chunksize)\n dataframes = [optimize_dataframe(chunk, category_threshold=DEFAULT_THRESHOLD) for chunk in chunks]\n data[file_name] = pd.concat(dataframes, ignore_index=True)\n else:\n print(f\"File {file_path} not found!\")\n return data\n\ndef check_data_integrity(original_df, optimized_df):\n \"\"\"\n Checks if the optimization did not alter the data.\n \"\"\"\n try:\n pd.testing.assert_frame_equal(original_df, optimized_df, check_like=True)\n print(\"Integrity check passed: No changes in data after optimization.\")\n except AssertionError as e:\n print(f\"Data integrity check failed: {e}\")\n\ndef check_duplicates(df):\n \"\"\"\n Checks for duplicates in the DataFrame.\n \"\"\"\n duplicates = df[df.duplicated(keep=False)]\n if not duplicates.empty:\n print(f\"Warning: Duplicates found in the dataset. Number of duplicates: {duplicates.shape[0]}\")\n return duplicates\n else:\n print(\"No duplicates found.\")\n return None\n\ndef test_thresholds(df):\n \"\"\"\n Tests different thresholds for DataFrame optimization.\n \"\"\"\n thresholds = np.linspace(0.1, 0.9, 9)\n memory_usages = []\n for threshold in thresholds:\n optimized_df = optimize_dataframe(df.copy(), category_threshold=threshold)\n memory_usages.append(optimized_df.memory_usage(deep=True).sum() / 1024**2)\n return thresholds, memory_usages\n\ndef plot_memory_usage(thresholds, memory_usages):\n \"\"\"\n Plots memory usage versus thresholds.\n \"\"\"\n plt.figure(figsize=(10, 6))\n plt.plot(thresholds, memory_usages, marker='o', linestyle='-')\n plt.title(\"Memory Usage vs. Threshold\")\n plt.xlabel(\"Threshold\")\n plt.ylabel(\"Memory Usage (MB)\")\n plt.grid(True)\n plt.show()\n\ndef analyze_sequence_data(df_sequences):\n \"\"\"\n Analyzes RNA sequence data.\n \"\"\"\n # Basic information\n print(f\"Total sequences: {len(df_sequences)}\")\n print(f\"Available columns: {df_sequences.columns.tolist()}\")\n \n # Sequence analysis\n if 'sequence' in df_sequences.columns:\n # Distribution of sequence lengths\n seq_lengths = df_sequences['sequence'].apply(len)\n print(f\"\\nSequence length statistics:\")\n print(f\"Minimum: {seq_lengths.min()}\")\n print(f\"Maximum: {seq_lengths.max()}\")\n print(f\"Average: {seq_lengths.mean():.2f}\")\n \n # Nucleotide count\n nucleotides = ['A', 'C', 'G', 'U']\n nucleotide_counts = {n: df_sequences['sequence'].str.count(n).sum() for n in nucleotides}\n total_nucleotides = sum(nucleotide_counts.values())\n \n print(\"\\nNucleotide distribution:\")\n for n, count in nucleotide_counts.items():\n print(f\"{n}: {count} ({count/total_nucleotides*100:.2f}%)\")\n \n return df_sequences\n\ndef analyze_label_data(df_labels):\n \"\"\"\n Analyzes 3D coordinate data (labels).\n \"\"\"\n print(f\"Total entries in labels: {len(df_labels)}\")\n print(f\"Available columns: {df_labels.columns.tolist()}\")\n \n # Analysis of 3D coordinates if available\n coord_columns = [col for col in df_labels.columns if col.startswith(('x_', 'y_', 'z_'))]\n if coord_columns:\n print(f\"\\nCoordinate columns found: {len(coord_columns)}\")\n \n # Basic statistics of coordinates\n for i in range(1, 6): # For the 5 possible structures\n x_col = f'x_{i}'\n y_col = f'y_{i}'\n z_col = f'z_{i}'\n \n if x_col in df_labels.columns and y_col in df_labels.columns and z_col in df_labels.columns:\n print(f\"\\nStatistics for structure {i}:\")\n print(f\"X - Mean: {df_labels[x_col].mean():.2f}, Std: {df_labels[x_col].std():.2f}\")\n print(f\"Y - Mean: {df_labels[y_col].mean():.2f}, Std: {df_labels[y_col].std():.2f}\")\n print(f\"Z - Mean: {df_labels[z_col].mean():.2f}, Std: {df_labels[z_col].std():.2f}\")\n \n return df_labels\n\ndef create_submission_template(test_df, sample_submission_df):\n \"\"\"\n Creates a submission template based on test data.\n \"\"\"\n # Check if sample_submission.csv is available\n if sample_submission_df is None:\n print(\"Sample submission file not found. Creating a new template.\")\n \n # Create a new DataFrame for submission\n submission_df = pd.DataFrame()\n \n # Example code to fill the template (adjust as needed)\n ids = []\n resnames = []\n resids = []\n \n for _, row in test_df.iterrows():\n sequence = row['sequence']\n target_id = row['target_id']\n \n for i, nucleotide in enumerate(sequence, 1):\n ids.append(f\"{target_id}_{i}\")\n resnames.append(nucleotide)\n resids.append(i)\n \n submission_df['ID'] = ids\n submission_df['resname'] = resnames\n submission_df['resid'] = resids\n \n # Add coordinate columns (5 structures)\n for i in range(1, 6):\n submission_df[f'x_{i}'] = 0.0\n submission_df[f'y_{i}'] = 0.0\n submission_df[f'z_{i}'] = 0.0\n else:\n submission_df = sample_submission_df.copy()\n print(\"Submission template created based on the provided example.\")\n \n return submission_df\n\ndef main():\n start_time = time.time()\n \n # Load main data\n print(\"Loading main data...\")\n main_data = load_main_data()\n \n # Check which files were loaded\n print(\"\\nLoaded files:\")\n for file_name, df in main_data.items():\n print(f\"- {file_name}: {df.shape if df is not None else 'Not found'}\")\n \n # Analyze training sequence data\n if \"train_sequences.csv\" in main_data:\n print(\"\\n===== Training Sequences Analysis =====\")\n analyze_sequence_data(main_data[\"train_sequences.csv\"])\n \n # Analyze training label data\n if \"train_labels.csv\" in main_data:\n print(\"\\n===== Training Labels Analysis =====\")\n analyze_label_data(main_data[\"train_labels.csv\"])\n \n # Check for duplicates in training data\n if \"train_sequences.csv\" in main_data:\n print(\"\\nChecking for duplicates in training sequences...\")\n check_duplicates(main_data[\"train_sequences.csv\"])\n \n # Create submission template\n if \"test_sequences.csv\" in main_data:\n print(\"\\nCreating submission template...\")\n submission_template = create_submission_template(\n main_data[\"test_sequences.csv\"],\n main_data.get(\"sample_submission.csv\")\n )\n print(f\"Submission template shape: {submission_template.shape}\")\n print(f\"First rows of the submission template:\")\n print(submission_template.head())\n \n # Calculate execution time\n end_time = time.time()\n print(f\"\\nRuntime: {end_time - start_time:.2f} seconds\")\n \n return main_data\n\nif __name__ == '__main__':\n main_data = main()","metadata":{"_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","execution":{"iopub.status.busy":"2025-04-09T15:13:28.415225Z","iopub.execute_input":"2025-04-09T15:13:28.416531Z","iopub.status.idle":"2025-04-09T15:13:28.986664Z","shell.execute_reply.started":"2025-04-09T15:13:28.416496Z","shell.execute_reply":"2025-04-09T15:13:28.985892Z"},"papermill":{"duration":0.867643,"end_time":"2025-03-26T03:42:58.386472","exception":false,"start_time":"2025-03-26T03:42:57.518829","status":"completed"},"tags":[],"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## Directory Explorer & CSV Verification for RNA3D 🗂️🔬","metadata":{"papermill":{"duration":0.010576,"end_time":"2025-03-26T03:42:58.408104","exception":false,"start_time":"2025-03-26T03:42:58.397528","status":"completed"},"tags":[]}},{"cell_type":"code","source":"# Updated main directory\ndir_main = \"/kaggle/input/stanford-rna-3d-folding/\"\n\n# List all files and directories in the main directory\ntry:\n all_files = os.listdir(dir_main)\n print(f\"All files and directories in '{dir_main}':\")\n \n for file in all_files:\n # Check if it's a file or directory\n full_path = os.path.join(dir_main, file)\n type_desc = \"directory\" if os.path.isdir(full_path) else \"file\"\n size = os.path.getsize(full_path) / 1024 # Size in KB\n print(f\" - {file} ({type_desc}, {size:.2f} KB)\")\n \n # If it's a directory, list up to 5 files inside it\n if os.path.isdir(full_path):\n try:\n internal_files = os.listdir(full_path)[:5] # Limit to 5 files\n if internal_files:\n print(f\" First files in '{file}':\")\n for internal_file in internal_files:\n print(f\" * {internal_file}\")\n if len(os.listdir(full_path)) > 5:\n print(f\" * ... and {len(os.listdir(full_path)) - 5} more file(s)\")\n else:\n print(f\" '{file}' is empty\")\n except Exception as e:\n print(f\" Error listing contents of '{file}': {e}\")\nexcept Exception as e:\n print(f\"Error listing directory {dir_main}: {e}\")\n\n# Check the structure of the main CSV files\nmain_files = [\n \"train_sequences.csv\", \n \"train_labels.csv\", \n \"validation_sequences.csv\", \n \"validation_labels.csv\", \n \"test_sequences.csv\",\n \"sample_submission.csv\"\n]\nprint(\"\\nChecking main CSV files:\")\n\nfor file in main_files:\n full_path = os.path.join(dir_main, file)\n if os.path.exists(full_path):\n # Get file size\n size_mb = os.path.getsize(full_path) / (1024 * 1024) # Size in MB\n \n # Read the first lines to check the structure\n try:\n import pandas as pd\n df = pd.read_csv(full_path, nrows=1)\n print(f\"\\n{file} ({size_mb:.2f} MB):\")\n print(f\"Columns: {df.columns.tolist()}\")\n print(f\"Example:\")\n print(df.head())\n except Exception as e:\n print(f\"Error reading {file}: {e}\")\n else:\n print(f\"{file} not found.\")","metadata":{"execution":{"iopub.status.busy":"2025-04-09T15:13:28.988033Z","iopub.execute_input":"2025-04-09T15:13:28.988313Z","iopub.status.idle":"2025-04-09T15:13:29.058857Z","shell.execute_reply.started":"2025-04-09T15:13:28.988291Z","shell.execute_reply":"2025-04-09T15:13:29.058153Z"},"papermill":{"duration":0.094882,"end_time":"2025-03-26T03:42:58.513093","exception":false,"start_time":"2025-03-26T03:42:58.418211","status":"completed"},"tags":[],"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## RNA3D Data Checker 🔍🧬","metadata":{"papermill":{"duration":0.010053,"end_time":"2025-03-26T03:42:58.533702","exception":false,"start_time":"2025-03-26T03:42:58.523649","status":"completed"},"tags":[]}},{"cell_type":"code","source":"# Updated main directory\ndir_main = \"/kaggle/input/stanford-rna-3d-folding/\"\n\ndef load_data():\n \"\"\"\n Loads the main CSV files from the Stanford RNA 3D Folding competition.\n Returns a dictionary with DataFrames.\n \"\"\"\n main_files = [\n \"train_sequences.csv\", \n \"train_labels.csv\", \n \"validation_sequences.csv\", \n \"validation_labels.csv\", \n \"test_sequences.csv\",\n \"sample_submission.csv\"\n ]\n \n data = {}\n for file_name in main_files:\n file_path = os.path.join(dir_main, file_name)\n if os.path.exists(file_path):\n try:\n data[file_name] = pd.read_csv(file_path)\n print(f\"File {file_name} loaded successfully. Shape: {data[file_name].shape}\")\n except Exception as e:\n print(f\"Error loading {file_name}: {e}\")\n else:\n print(f\"File {file_name} not found.\")\n data[file_name] = None\n \n return data\n\ndef compare_columns(main_data):\n \"\"\"\n Compares columns between different DataFrames.\n \"\"\"\n # List all available keys\n print(\"\\nLoaded files:\")\n print(list(main_data.keys()))\n \n # Compare columns between train_sequences.csv and test_sequences.csv\n if \"train_sequences.csv\" in main_data and \"test_sequences.csv\" in main_data:\n train_cols = set(main_data[\"train_sequences.csv\"].columns)\n test_cols = set(main_data[\"test_sequences.csv\"].columns)\n \n print(\"\\nColumns in train_sequences.csv:\")\n print(list(main_data[\"train_sequences.csv\"].columns))\n \n print(\"\\nUnique columns in train_sequences.csv (not present in test_sequences.csv):\")\n print(train_cols - test_cols)\n \n print(\"\\nUnique columns in test_sequences.csv (not present in train_sequences.csv):\")\n print(test_cols - train_cols)\n \n # Compare columns between train_labels.csv and validation_labels.csv\n if \"train_labels.csv\" in main_data and \"validation_labels.csv\" in main_data:\n train_label_cols = set(main_data[\"train_labels.csv\"].columns)\n val_label_cols = set(main_data[\"validation_labels.csv\"].columns)\n \n print(\"\\nColumns in train_labels.csv:\")\n print(list(main_data[\"train_labels.csv\"].columns))\n \n print(\"\\nColumns in validation_labels.csv:\")\n print(list(main_data[\"validation_labels.csv\"].columns))\n \n print(\"\\nUnique columns in validation_labels.csv (not present in train_labels.csv):\")\n print(val_label_cols - train_label_cols)\n \n # Compare columns between validation_labels.csv and sample_submission.csv\n if \"validation_labels.csv\" in main_data and \"sample_submission.csv\" in main_data:\n val_label_cols = set(main_data[\"validation_labels.csv\"].columns)\n sample_cols = set(main_data[\"sample_submission.csv\"].columns)\n \n print(\"\\nColumns in sample_submission.csv:\")\n print(list(main_data[\"sample_submission.csv\"].columns))\n \n print(\"\\nUnique columns in validation_labels.csv (not present in sample_submission.csv):\")\n print(val_label_cols - sample_cols)\n \n print(\"\\nUnique columns in sample_submission.csv (not present in validation_labels.csv):\")\n print(sample_cols - val_label_cols)\n\ndef analyze_structure_format(main_data):\n \"\"\"\n Analyzes the format of 3D structures (coordinates).\n \"\"\"\n if \"validation_labels.csv\" in main_data and main_data[\"validation_labels.csv\"] is not None:\n df = main_data[\"validation_labels.csv\"]\n \n # Find all coordinate columns (x_1, y_1, z_1, etc.)\n coord_cols = [col for col in df.columns if col.startswith(('x_', 'y_', 'z_'))]\n \n # Group by structure\n structures = {}\n for col in coord_cols:\n # Extract structure number (e.g., \"x_1\" -> 1)\n parts = col.split('_')\n if len(parts) == 2:\n struct_num = int(parts[1])\n coord_type = parts[0]\n \n if struct_num not in structures:\n structures[struct_num] = []\n \n structures[struct_num].append(col)\n \n print(\"\\nStructure of the labels file:\")\n print(f\"Total structures found: {len(structures)}\")\n \n # Show details of the first structure\n if structures:\n first_struct = min(structures.keys())\n print(f\"\\nDetails of structure {first_struct}:\")\n print(f\"Columns: {sorted(structures[first_struct])}\")\n \n # Check for missing values\n for col in structures[first_struct]:\n missing = df[col].isna().sum()\n total = len(df)\n print(f\"{col}: {missing} missing values ({missing/total*100:.2f}%)\")\n \n # Check the range of non-missing values for the first structure\n for col in structures[first_struct]:\n non_null = df[col][df[col] != -1.0e+18] # Values that are not -1.0e+18\n if not non_null.empty:\n print(f\"{col} - Range: [{non_null.min():.3f}, {non_null.max():.3f}]\")\n\ndef main():\n # Load the data\n main_data = load_data()\n \n # Compare columns between different files\n compare_columns(main_data)\n \n # Analyze the format of 3D structures\n analyze_structure_format(main_data)\n \n return main_data\n\nif __name__ == '__main__':\n main_data = main()","metadata":{"execution":{"iopub.status.busy":"2025-04-09T15:13:29.060153Z","iopub.execute_input":"2025-04-09T15:13:29.060453Z","iopub.status.idle":"2025-04-09T15:13:29.294091Z","shell.execute_reply.started":"2025-04-09T15:13:29.060431Z","shell.execute_reply":"2025-04-09T15:13:29.293300Z"},"papermill":{"duration":0.270583,"end_time":"2025-03-26T03:42:58.814635","exception":false,"start_time":"2025-03-26T03:42:58.544052","status":"completed"},"tags":[],"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## Integrated RNA3D Sequence and Structure Analyzer 🔬🧬","metadata":{"papermill":{"duration":0.010382,"end_time":"2025-03-26T03:42:58.836058","exception":false,"start_time":"2025-03-26T03:42:58.825676","status":"completed"},"tags":[]}},{"cell_type":"code","source":"# Initialize seed to control randomness\nnp.random.seed(0)\n\n# Directories and files adjusted for the new competition\nDATA_DIR = os.getenv('DATA_DIR', '/kaggle/input/stanford-rna-3d-folding/')\nmain_files = [\n \"train_sequences.csv\", \n \"train_labels.csv\", \n \"validation_sequences.csv\", \n \"validation_labels.csv\", \n \"test_sequences.csv\",\n \"sample_submission.csv\"\n]\n\nDEFAULT_THRESHOLD = 0.4 # Default threshold after analysis\n\ndef optimize_dataframe(df, inplace=False, category_threshold=DEFAULT_THRESHOLD):\n \"\"\"\n Optimizes the DataFrame to save memory.\n \"\"\"\n if category_threshold < 0 or category_threshold > 1:\n raise ValueError(\"category_threshold must be between 0 and 1.\")\n \n if not inplace:\n df = df.copy()\n \n for col in df.columns:\n col_type = df[col].dtype\n if np.issubdtype(col_type, np.integer):\n c_min, c_max = df[col].min(), df[col].max()\n if c_min > np.iinfo(np.int8).min and c_max < np.iinfo(np.int8).max:\n df[col] = df[col].astype(np.int8)\n elif c_min > np.iinfo(np.int16).min and c_max < np.iinfo(np.int16).max:\n df[col] = df[col].astype(np.int16)\n elif c_min > np.iinfo(np.int32).min and c_max < np.iinfo(np.int32).max:\n df[col] = df[col].astype(np.int32)\n elif np.issubdtype(col_type, np.floating):\n # First check if it's not the special value -1.0e+18\n if df[col].min() > np.finfo(np.float32).min and df[col].max() < np.finfo(np.float32).max:\n df[col] = df[col].astype(np.float32)\n if col_type == object:\n unique_vals = len(df[col].unique())\n if unique_vals / len(df) < category_threshold:\n df[col] = df[col].astype('category')\n \n return df\n\ndef load_main_data(chunksize=50000):\n \"\"\"\n Loads the main files.\n \"\"\"\n data = {}\n for file_name in main_files:\n file_path = os.path.join(DATA_DIR, file_name)\n if os.path.exists(file_path):\n chunks = pd.read_csv(file_path, on_bad_lines='skip', low_memory=False, chunksize=chunksize)\n dataframes = [optimize_dataframe(chunk, category_threshold=DEFAULT_THRESHOLD) for chunk in chunks]\n data[file_name] = pd.concat(dataframes, ignore_index=True)\n print(f\"File {file_name} loaded successfully. Shape: {data[file_name].shape}\")\n else:\n print(f\"File {file_path} not found!\")\n data[file_name] = None\n return data\n\ndef filter_columns_by_prefix(df, prefix=\"x_\"):\n \"\"\"\n Filters and counts the number of columns in a DataFrame based on a provided prefix.\n \n :param df: DataFrame where filtering will be applied.\n :param prefix: Prefix to be used for filtering. Ex: \"x_\", \"y_\", \"z_\".\n :return: List of filtered columns.\n \"\"\"\n filtered_columns = [col for col in df.columns if col.startswith(prefix)]\n return filtered_columns\n\ndef count_nucleotides(df, column_name='sequence'):\n \"\"\"\n Counts the frequency of each nucleotide in a specific column of a DataFrame.\n \n :param df: DataFrame containing the sequences.\n :param column_name: Name of the column containing the sequences. Default is 'sequence'.\n :return: Counter object with the nucleotide counts.\n \"\"\"\n from collections import Counter\n\n # Check if the column exists in the DataFrame\n if column_name not in df.columns:\n raise ValueError(f\"Column '{column_name}' not found in DataFrame.\")\n \n # Concatenate all sequences and count nucleotides\n all_sequences = ''.join(df[column_name].tolist())\n nucleotide_counts = Counter(all_sequences)\n \n return nucleotide_counts\n\ndef get_columns_without_missing_values(df):\n \"\"\"\n Returns columns without any missing values in the DataFrame.\n \n :param df: DataFrame to be checked.\n :return: List of columns without missing values.\n \"\"\"\n missing_values = df.isnull().sum()\n return missing_values[missing_values == 0].index.tolist()\n\ndef get_empty_columns(df):\n \"\"\"\n Returns columns that are completely empty in the DataFrame.\n \n :param df: DataFrame to be checked.\n :return: List of empty columns.\n \"\"\"\n missing_values = df.isnull().sum()\n return missing_values[missing_values == df.shape[0]].index.tolist()\n\ndef plot_coord_distributions(df_labels, prefix='x_', max_structures=5):\n \"\"\"\n Plots the distribution of coordinates (x, y, or z) for up to max_structures structures.\n \n :param df_labels: DataFrame containing the coordinates.\n :param prefix: Prefix of columns to be plotted ('x_', 'y_', or 'z_').\n :param max_structures: Maximum number of structures to show.\n \"\"\"\n # Find coordinate columns with the specified prefix\n coord_cols = filter_columns_by_prefix(df_labels, prefix)\n \n # Limit to the maximum number of structures\n coord_cols = sorted(coord_cols)[:max_structures]\n \n if not coord_cols:\n print(f\"No column with prefix '{prefix}' found.\")\n return\n \n # Set up the plot\n fig, axes = plt.subplots(1, len(coord_cols), figsize=(16, 4))\n if len(coord_cols) == 1:\n axes = [axes] # Ensure axes is iterable even with a single subplot\n \n # Plot histograms for each column\n for i, col in enumerate(coord_cols):\n # Filter special values (-1.0e+18) if present\n values = df_labels[col]\n filtered_values = values[values > -1.0e+17] # Cutoff value to filter -1.0e+18\n \n axes[i].hist(filtered_values, bins=30, alpha=0.7)\n axes[i].set_title(f'Distribution of {col}')\n axes[i].set_xlabel('Value')\n axes[i].set_ylabel('Frequency')\n \n plt.tight_layout()\n plt.show()\n\ndef analyze_3d_structure(df_labels):\n \"\"\"\n Analyzes the 3D coordinates of RNA structures.\n \n :param df_labels: DataFrame containing 3D coordinates.\n \"\"\"\n # Find all coordinate columns\n x_cols = filter_columns_by_prefix(df_labels, 'x_')\n y_cols = filter_columns_by_prefix(df_labels, 'y_')\n z_cols = filter_columns_by_prefix(df_labels, 'z_')\n \n print(f\"Number of x columns: {len(x_cols)}\")\n print(f\"Number of y columns: {len(y_cols)}\")\n print(f\"Number of z columns: {len(z_cols)}\")\n \n # Check for missing or special values in coordinates\n special_value = -1.0e+18 # Special value observed in the data\n \n for i, (x_col, y_col, z_col) in enumerate(zip(x_cols, y_cols, z_cols), 1):\n # Count missing or special values\n x_special = (df_labels[x_col] == special_value).sum()\n y_special = (df_labels[y_col] == special_value).sum()\n z_special = (df_labels[z_col] == special_value).sum()\n \n x_null = df_labels[x_col].isnull().sum()\n y_null = df_labels[y_col].isnull().sum()\n z_null = df_labels[z_col].isnull().sum()\n \n # Count how many complete structures exist (all x, y, z are neither special nor null)\n valid_structures = ((df_labels[x_col] != special_value) & \n (df_labels[y_col] != special_value) & \n (df_labels[z_col] != special_value) &\n df_labels[x_col].notnull() & \n df_labels[y_col].notnull() & \n df_labels[z_col].notnull()).sum()\n \n total_rows = len(df_labels)\n \n print(f\"\\nStructure {i}:\")\n print(f\" Special values: x={x_special} ({x_special/total_rows*100:.2f}%), y={y_special} ({y_special/total_rows*100:.2f}%), z={z_special} ({z_special/total_rows*100:.2f}%)\")\n print(f\" Null values: x={x_null} ({x_null/total_rows*100:.2f}%), y={y_null} ({y_null/total_rows*100:.2f}%), z={z_null} ({z_null/total_rows*100:.2f}%)\")\n print(f\" Complete structures: {valid_structures} ({valid_structures/total_rows*100:.2f}%)\")\n \n # Limit analysis to the first 5 structures\n if i >= 5:\n print(\"\\nAnalysis limited to the first 5 structures.\")\n break\n\ndef analyze_sequences(df_sequences):\n \"\"\"\n Analyzes RNA sequences.\n \n :param df_sequences: DataFrame containing the 'sequence' column.\n \"\"\"\n # Basic statistics of the sequence column\n print(\"\\nBasic statistics of the 'sequence' column:\")\n print(df_sequences['sequence'].describe())\n \n # Sequence lengths\n seq_lengths = df_sequences['sequence'].apply(len)\n print(\"\\nSequence length statistics:\")\n print(f\"Minimum: {seq_lengths.min()}\")\n print(f\"Maximum: {seq_lengths.max()}\")\n print(f\"Mean: {seq_lengths.mean():.2f}\")\n print(f\"Median: {seq_lengths.median()}\")\n \n # Nucleotide counts\n nucleotide_counts = count_nucleotides(df_sequences)\n total_nucleotides = sum(nucleotide_counts.values())\n \n print(\"\\nNucleotide distribution:\")\n for nucleotide, count in sorted(nucleotide_counts.items()):\n print(f\"{nucleotide}: {count} ({count/total_nucleotides*100:.2f}%)\")\n \n # Plot length distribution\n plt.figure(figsize=(10, 6))\n plt.hist(seq_lengths, bins=30, alpha=0.7)\n plt.title('Sequence Length Distribution')\n plt.xlabel('Length')\n plt.ylabel('Frequency')\n plt.grid(True, alpha=0.3)\n plt.show()\n\ndef main():\n # Load main data\n main_data = load_main_data()\n\n # Check which files were loaded\n print(\"\\nLoaded files:\")\n for file_name, df in main_data.items():\n if df is not None:\n print(f\"- {file_name}: {df.shape}\")\n \n # Analyze 3D structures in validation_labels.csv\n if \"validation_labels.csv\" in main_data and main_data[\"validation_labels.csv\"] is not None:\n print(\"\\n===== Analysis of 3D Structures (validation_labels.csv) =====\")\n df_labels = main_data[\"validation_labels.csv\"]\n \n # Count coordinate columns\n x_cols = filter_columns_by_prefix(df_labels, 'x_')\n y_cols = filter_columns_by_prefix(df_labels, 'y_')\n z_cols = filter_columns_by_prefix(df_labels, 'z_')\n \n print(f\"There are {len(x_cols)} x_ columns in the DataFrame.\")\n print(f\"There are {len(y_cols)} y_ columns in the DataFrame.\")\n print(f\"There are {len(z_cols)} z_ columns in the DataFrame.\")\n \n # Identify columns without missing values\n columns_without_missing = get_columns_without_missing_values(df_labels)\n print(f\"\\nColumns without missing values: {len(columns_without_missing)}\")\n \n # Identify completely empty columns\n empty_columns = get_empty_columns(df_labels)\n print(f\"Completely empty columns: {len(empty_columns)}\")\n \n # Analyze 3D coordinates in detail\n analyze_3d_structure(df_labels)\n \n # Plot distribution of x, y, z coordinates for the first structures\n print(\"\\nDistribution of X coordinates:\")\n plot_coord_distributions(df_labels, 'x_', max_structures=3)\n print(\"\\nDistribution of Y coordinates:\")\n plot_coord_distributions(df_labels, 'y_', max_structures=3)\n print(\"\\nDistribution of Z coordinates:\")\n plot_coord_distributions(df_labels, 'z_', max_structures=3)\n \n # Analyze sequences in train_sequences.csv\n if \"train_sequences.csv\" in main_data and main_data[\"train_sequences.csv\"] is not None:\n print(\"\\n===== Analysis of Sequences (train_sequences.csv) =====\")\n df_sequences = main_data[\"train_sequences.csv\"]\n \n # First few rows of the sequence column\n print(\"\\nFirst few rows of the 'sequence' column:\")\n print(df_sequences['sequence'].head())\n \n # Data type of the sequence column\n print(\"\\nData type of the 'sequence' column:\")\n print(df_sequences['sequence'].dtype)\n \n # Complete sequence analysis\n analyze_sequences(df_sequences)\n \n return main_data\n\nif __name__ == '__main__':\n main_data = main()","metadata":{"execution":{"iopub.status.busy":"2025-04-09T15:13:29.294890Z","iopub.execute_input":"2025-04-09T15:13:29.295115Z","iopub.status.idle":"2025-04-09T15:13:31.955184Z","shell.execute_reply.started":"2025-04-09T15:13:29.295093Z","shell.execute_reply":"2025-04-09T15:13:31.954290Z"},"papermill":{"duration":2.579715,"end_time":"2025-03-26T03:43:01.426154","exception":false,"start_time":"2025-03-26T03:42:58.846439","status":"completed"},"tags":[],"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## Data Preparation for RNA 3D Structure Prediction 🧬🔍","metadata":{"papermill":{"duration":0.013504,"end_time":"2025-03-26T03:43:01.454026","exception":false,"start_time":"2025-03-26T03:43:01.440522","status":"completed"},"tags":[]}},{"cell_type":"code","source":"# File paths\nDATA_DIR = \"/kaggle/input/stanford-rna-3d-folding/\"\nOUTPUT_DIR = \"/kaggle/working/\"\nos.makedirs(OUTPUT_DIR, exist_ok=True)\n\ndef load_data():\n \"\"\"\n Loads the necessary data for the competition.\n \"\"\"\n data = {}\n \n # Load sequences\n data['train_seq'] = pd.read_csv(os.path.join(DATA_DIR, \"train_sequences.csv\"))\n data['valid_seq'] = pd.read_csv(os.path.join(DATA_DIR, \"validation_sequences.csv\"))\n data['test_seq'] = pd.read_csv(os.path.join(DATA_DIR, \"test_sequences.csv\"))\n \n # Load structures (labels)\n data['train_labels'] = pd.read_csv(os.path.join(DATA_DIR, \"train_labels.csv\"))\n data['valid_labels'] = pd.read_csv(os.path.join(DATA_DIR, \"validation_labels.csv\"))\n \n # Load submission format\n data['sample_submission'] = pd.read_csv(os.path.join(DATA_DIR, \"sample_submission.csv\"))\n \n return data\n\ndef analyze_id_structure(data_dict):\n \"\"\"\n Analyzes the ID structure in different files to understand the correct mapping.\n \"\"\"\n # We'll analyze the specific formats for train and valid\n \n # 1. Analysis of training labels\n train_label_ids = data_dict['train_labels']['ID'].tolist()\n print(f\"Total IDs in training labels: {len(train_label_ids)}\")\n print(f\"Number of unique IDs: {len(set(train_label_ids))}\")\n \n # Try to understand the ID format in the labels file\n train_id_parts = {}\n for id_str in train_label_ids[:100]: # Analyze the first 100\n parts = id_str.split('_')\n num_parts = len(parts)\n if num_parts not in train_id_parts:\n train_id_parts[num_parts] = []\n train_id_parts[num_parts].append(parts)\n \n print(\"\\nID formats found in train_labels:\")\n for num_parts, examples in train_id_parts.items():\n print(f\"\\nFormat with {num_parts} parts:\")\n for i, parts in enumerate(examples[:3]):\n print(f\" Example {i+1}: {parts}\")\n \n # 2. Analysis of training sequences\n train_seq_ids = data_dict['train_seq']['target_id'].tolist()\n print(f\"\\nTotal IDs in training sequences: {len(train_seq_ids)}\")\n print(f\"Number of unique IDs: {len(set(train_seq_ids))}\")\n \n # Try to understand the ID format in the sequences file\n train_seq_id_parts = {}\n for id_str in train_seq_ids[:100]: # Analyze the first 100\n parts = id_str.split('_')\n num_parts = len(parts)\n if num_parts not in train_seq_id_parts:\n train_seq_id_parts[num_parts] = []\n train_seq_id_parts[num_parts].append(parts)\n \n print(\"\\nID formats found in train_sequences:\")\n for num_parts, examples in train_seq_id_parts.items():\n print(f\"\\nFormat with {num_parts} parts:\")\n for i, parts in enumerate(examples[:3]):\n print(f\" Example {i+1}: {parts}\")\n \n # 3. Analysis of validation labels\n valid_label_ids = data_dict['valid_labels']['ID'].tolist()\n print(f\"\\nTotal IDs in validation labels: {len(valid_label_ids)}\")\n print(f\"Number of unique IDs: {len(set(valid_label_ids))}\")\n \n # Count unique sequence IDs in validation labels\n valid_seq_ids_from_labels = set([id_str.split('_')[0] for id_str in valid_label_ids])\n print(f\"Number of unique sequence IDs in validation labels: {len(valid_seq_ids_from_labels)}\")\n print(f\"Examples: {list(valid_seq_ids_from_labels)[:5]}\")\n \n # 4. Analysis of validation sequences\n valid_seq_ids = data_dict['valid_seq']['target_id'].tolist()\n print(f\"\\nTotal IDs in validation sequences: {len(valid_seq_ids)}\")\n print(f\"Number of unique IDs: {len(set(valid_seq_ids))}\")\n print(f\"Examples: {valid_seq_ids[:5]}\")\n \n # 5. Check correspondence between unique IDs\n overlap_valid = set(valid_seq_ids).intersection(valid_seq_ids_from_labels)\n print(f\"\\nCorrespondence between validation sequences and labels: {len(overlap_valid)} of {len(valid_seq_ids)}\")\n \n # 6. Check how sequences and residues relate\n if len(overlap_valid) > 0:\n sample_id = list(overlap_valid)[0]\n sample_seq = data_dict['valid_seq'][data_dict['valid_seq']['target_id'] == sample_id]['sequence'].iloc[0]\n sample_labels = data_dict['valid_labels'][data_dict['valid_labels']['ID'].str.startswith(f\"{sample_id}_\")]\n \n print(f\"\\nAnalysis for sequence ID: {sample_id}\")\n print(f\"Sequence length: {len(sample_seq)}\")\n print(f\"Number of residues in labels: {len(sample_labels)}\")\n \n # Check how residue numbers are related\n residue_numbers = sample_labels['resid'].sort_values().tolist()\n print(f\"First residue numbers: {residue_numbers[:10]}\")\n print(f\"Last residue numbers: {residue_numbers[-10:]}\")\n \n return train_id_parts, train_seq_id_parts, overlap_valid\n\ndef fix_train_mapping(train_seq_df, train_labels_df):\n \"\"\"\n Identifies a correct mapping between train_sequences.csv and train_labels.csv\n using the ID format from the validation file as a reference.\n \n This is necessary because there's no obvious direct correspondence between the IDs.\n \"\"\"\n # First, extract the prefix of the ID from labels (format: XX_Y_Z)\n train_labels_df['seq_id'] = train_labels_df['ID'].apply(lambda x: x.split('_')[0] + '_' + x.split('_')[1])\n \n # Check if this format corresponds to the format of sequence IDs\n seq_ids_set = set(train_seq_df['target_id'])\n label_seq_ids_set = set(train_labels_df['seq_id'])\n \n overlap = seq_ids_set.intersection(label_seq_ids_set)\n print(f\"Overlap after format adjustment: {len(overlap)} of {len(seq_ids_set)}\")\n \n if len(overlap) > 0:\n print(f\"Examples of matching IDs: {list(overlap)[:5]}\")\n return overlap\n \n # If it still doesn't work, we need to analyze the structure in more detail\n print(\"No matches found, checking other formats...\")\n \n # Try other possible formats\n formats_to_try = [\n lambda x: x.split('_')[0], # Only first part\n lambda x: '_'.join(x.split('_')[:2]), # First two parts\n lambda x: x.split('_')[0] + '_' + x.split('_')[1][0], # First part + first letter of second part\n ]\n \n for i, format_func in enumerate(formats_to_try):\n train_labels_df[f'seq_id_{i}'] = train_labels_df['ID'].apply(format_func)\n label_seq_ids_set = set(train_labels_df[f'seq_id_{i}'])\n overlap = seq_ids_set.intersection(label_seq_ids_set)\n print(f\"Format {i}: Overlap = {len(overlap)} of {len(seq_ids_set)}\")\n \n if len(overlap) > 0:\n print(f\"Examples of matching IDs: {list(overlap)[:5]}\")\n return overlap, f'seq_id_{i}'\n \n # If no match is found, create a mapping based on observed patterns\n print(\"No matches found using simple patterns.\")\n print(\"Creating a manual mapping based on data structure...\")\n \n # Group labels by first parts of ID\n train_labels_df['prefix'] = train_labels_df['ID'].apply(lambda x: x.split('_')[0])\n label_groups = train_labels_df.groupby('prefix')\n \n # For each sequence, find the best match based on number of residues\n mapping = {}\n for _, seq_row in train_seq_df.iterrows():\n seq_id = seq_row['target_id']\n seq_length = len(seq_row['sequence'])\n \n best_match = None\n best_diff = float('inf')\n \n for prefix, group in label_groups:\n residue_count = len(group)\n diff = abs(residue_count - seq_length)\n \n if diff < best_diff:\n best_diff = diff\n best_match = prefix\n \n # Consider a match only if the number of residues is close\n if best_diff <= 10: # Tolerance of 10 residues\n mapping[seq_id] = best_match\n \n print(f\"Manual mapping created with {len(mapping)} matches\")\n return mapping\n\ndef create_mapping_valid(valid_seq_df, valid_labels_df):\n \"\"\"\n Creates a mapping between validation sequences and their coordinates.\n \n In this case, the IDs already correspond directly (R1107 -> R1107_1, R1107_2, etc.)\n \"\"\"\n # Check which ID format is used in the validation set\n valid_labels_df['seq_id'] = valid_labels_df['ID'].apply(lambda x: x.split('_')[0])\n \n # Check overlap\n seq_ids = set(valid_seq_df['target_id'])\n label_seq_ids = set(valid_labels_df['seq_id'])\n \n overlap = seq_ids.intersection(label_seq_ids)\n print(f\"Correspondence for validation: {len(overlap)} of {len(seq_ids)}\")\n \n mapping = {}\n for seq_id in overlap:\n # Get sequence\n seq = valid_seq_df[valid_seq_df['target_id'] == seq_id]['sequence'].iloc[0]\n \n # Get all residues for this sequence\n residues = valid_labels_df[valid_labels_df['seq_id'] == seq_id].sort_values('resid')\n \n # Extract coordinates for all structures\n num_structures = 1\n for col in residues.columns:\n if col.startswith('x_'):\n struct_num = int(col.split('_')[1])\n num_structures = max(num_structures, struct_num)\n \n # Initialize structures\n structures = []\n \n for struct_idx in range(1, num_structures + 1):\n coords = []\n has_valid_coords = False\n \n # Check if this structure has coordinates\n if f'x_{struct_idx}' in residues.columns:\n for _, row in residues.iterrows():\n x = row[f'x_{struct_idx}']\n y = row[f'y_{struct_idx}']\n z = row[f'z_{struct_idx}']\n \n # Check if they are valid values\n if abs(x) < 1.0e+17 and abs(y) < 1.0e+17 and abs(z) < 1.0e+17:\n coords.append([x, y, z])\n has_valid_coords = True\n else:\n coords.append([np.nan, np.nan, np.nan])\n \n if has_valid_coords:\n structures.append(coords)\n \n # Add to mapping if there are valid structures\n if structures:\n mapping[seq_id] = {\n 'sequence': seq,\n 'structures': structures\n }\n \n print(f\"Mapping created with {len(mapping)} valid sequences\")\n return mapping\n\ndef create_processed_data(mapping, output_prefix):\n \"\"\"\n Creates and saves processed data from the mapping.\n \n Parameters:\n mapping: Dictionary with the mapping of sequences to structures\n output_prefix: Prefix for output files ('train' or 'valid')\n \n Returns:\n X, y: Arrays for training\n \"\"\"\n if not mapping:\n print(f\"WARNING: No valid mapping for {output_prefix}\")\n return None, None\n \n X_data = []\n y_data = []\n ids = []\n \n for seq_id, data in mapping.items():\n seq = data['sequence']\n structures = data['structures']\n \n # Skip if there are no structures\n if not structures:\n continue\n \n # Use the first valid structure\n structure = structures[0]\n \n # Check if the structure has valid coordinates for all residues\n if len(structure) != len(seq):\n print(f\"WARNING: Difference between sequence length ({len(seq)}) and coordinates ({len(structure)}) for {seq_id}\")\n # If needed, we could consider padding or truncation here\n continue\n \n # Create feature matrix (one-hot encoding)\n features = []\n for nucleotide in seq:\n if nucleotide == 'A':\n features.append([1, 0, 0, 0, 0])\n elif nucleotide == 'C':\n features.append([0, 1, 0, 0, 0])\n elif nucleotide == 'G':\n features.append([0, 0, 1, 0, 0])\n elif nucleotide == 'U':\n features.append([0, 0, 0, 1, 0])\n else:\n features.append([0, 0, 0, 0, 1]) # For unknown nucleotides\n \n X_data.append(np.array(features))\n y_data.append(np.array(structure))\n ids.append(seq_id)\n \n if not X_data:\n print(f\"WARNING: No valid processed data for {output_prefix}\")\n return None, None, []\n \n # Padding to ensure all sequences have the same length\n max_length = max(len(x) for x in X_data)\n X_padded = []\n y_padded = []\n \n for x, y in zip(X_data, y_data):\n if len(x) < max_length:\n x_pad = np.zeros((max_length, 5))\n x_pad[:len(x), :] = x\n \n y_pad = np.zeros((max_length, 3))\n y_pad[:len(y), :] = y\n \n X_padded.append(x_pad)\n y_padded.append(y_pad)\n else:\n X_padded.append(x)\n y_padded.append(y)\n \n X = np.array(X_padded)\n y = np.array(y_padded)\n \n # Save the processed data\n np.save(os.path.join(OUTPUT_DIR, f'X_{output_prefix}.npy'), X)\n np.save(os.path.join(OUTPUT_DIR, f'y_{output_prefix}.npy'), y)\n \n with open(os.path.join(OUTPUT_DIR, f'{output_prefix}_ids.txt'), 'w') as f:\n for id in ids:\n f.write(f\"{id}\\n\")\n \n print(f\"Processed data for {output_prefix}: X.shape = {X.shape}, y.shape = {y.shape}\")\n return X, y, ids\n\ndef explore_sequence_mapping(seq_id, mapping, data_dict):\n \"\"\"\n Explores a mapping example in detail for diagnostics.\n \"\"\"\n if seq_id not in mapping:\n print(f\"WARNING: Sequence ID {seq_id} not found in mapping\")\n return\n \n data = mapping[seq_id]\n seq = data['sequence']\n structures = data['structures']\n \n print(f\"Exploring mapping for sequence: {seq_id}\")\n print(f\"Sequence length: {len(seq)}\")\n print(f\"Number of available structures: {len(structures)}\")\n \n # Detail each structure\n for i, structure in enumerate(structures):\n print(f\"\\nStructure {i+1}:\")\n print(f\" Number of coordinates: {len(structure)}\")\n if len(structure) > 0:\n print(f\" First coordinates: {structure[:3]}\")\n print(f\" Last coordinates: {structure[-3:]}\")\n \n # Check correspondence with the sequence\n if len(structure) != len(seq):\n print(f\" WARNING: Difference between sequence length ({len(seq)}) and coordinates ({len(structure)})\")\n else:\n print(f\" Perfect match between sequence and coordinates\")\n\ndef main():\n # Load the data\n print(\"Loading data...\")\n data_dict = load_data()\n \n # Analyze ID structure to understand the mapping\n print(\"\\nAnalyzing ID structure...\")\n train_id_parts, train_seq_id_parts, overlap_valid = analyze_id_structure(data_dict)\n \n # For validation, the mapping is direct (R1107 -> R1107_1, R1107_2, etc.)\n print(\"\\nCreating mapping for validation data...\")\n valid_mapping = create_mapping_valid(data_dict['valid_seq'], data_dict['valid_labels'])\n \n # Explore a validation mapping example to verify\n if valid_mapping:\n sample_id = list(valid_mapping.keys())[0]\n print(f\"\\nExploring a validation mapping example ({sample_id}):\")\n explore_sequence_mapping(sample_id, valid_mapping, data_dict)\n \n # Create and save processed data for validation\n X_valid, y_valid, valid_ids = create_processed_data(valid_mapping, 'valid')\n \n # Since we couldn't establish a mapping for training,\n # we'll use validation data for training as well (transfer learning)\n print(\"\\nUsing validation data as training (due to lack of direct mapping)...\")\n X_train = X_valid\n y_train = y_valid\n train_ids = valid_ids\n \n if X_train is not None:\n np.save(os.path.join(OUTPUT_DIR, 'X_train.npy'), X_train)\n np.save(os.path.join(OUTPUT_DIR, 'y_train.npy'), y_train)\n \n with open(os.path.join(OUTPUT_DIR, 'train_ids.txt'), 'w') as f:\n for id in train_ids:\n f.write(f\"{id}\\n\")\n \n # Return the processed data\n return {\n 'X_train': X_train,\n 'y_train': y_train,\n 'X_valid': X_valid,\n 'y_valid': y_valid,\n 'valid_mapping': valid_mapping,\n 'valid_ids': valid_ids\n }\n\nif __name__ == \"__main__\":\n processed_data = main()","metadata":{"execution":{"iopub.status.busy":"2025-04-09T15:13:31.956237Z","iopub.execute_input":"2025-04-09T15:13:31.956621Z","iopub.status.idle":"2025-04-09T15:13:37.169843Z","shell.execute_reply.started":"2025-04-09T15:13:31.956587Z","shell.execute_reply":"2025-04-09T15:13:37.169031Z"},"papermill":{"duration":5.431349,"end_time":"2025-03-26T03:43:06.898790","exception":false,"start_time":"2025-03-26T03:43:01.467441","status":"completed"},"tags":[],"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## Heatmap Viewer for RNA Sequences 🔥🧬","metadata":{"papermill":{"duration":0.013025,"end_time":"2025-03-26T03:43:06.925988","exception":false,"start_time":"2025-03-26T03:43:06.912963","status":"completed"},"tags":[]}},{"cell_type":"code","source":"def visualize_rna_heatmap_from_processed_data(processed_data, num_samples=12):\n \"\"\"\n Visualizes a heatmap for RNA sequences using processed data.\n \n Parameters:\n processed_data: Dictionary with processed data returned by the main() function\n num_samples: Number of sequences to visualize\n \"\"\"\n try:\n # Check if we have the necessary data\n if 'X_valid' not in processed_data or processed_data['X_valid'] is None:\n print(\"Validation data not found in processed_data object\")\n return None\n \n # Get the data\n X_valid = processed_data['X_valid']\n print(f\"Data found with format: {X_valid.shape}\")\n \n # Limit to the number of samples\n X_valid_subset = X_valid[:num_samples]\n \n # If we have IDs, use them\n if 'valid_ids' in processed_data and processed_data['valid_ids']:\n valid_ids = processed_data['valid_ids'][:num_samples]\n else:\n valid_ids = [f\"Seq_{i+1}\" for i in range(X_valid_subset.shape[0])]\n \n # Convert one-hot encoding to nucleotide indices\n # Expected format: A=[1,0,0,0,0], C=[0,1,0,0,0], G=[0,0,1,0,0], U=[0,0,0,1,0], N=[0,0,0,0,1]\n sequences_matrix = np.argmax(X_valid_subset, axis=2)\n \n # Replace zeros (padding) with 4 (N/Unknown) when all values are zero\n is_padding = np.all(X_valid_subset == 0, axis=2)\n sequences_matrix[is_padding] = 4\n \n # Define a categorical colormap (distinct colors per nucleotide)\n cmap = mcolors.ListedColormap(['#3498db', '#2ecc71', '#e74c3c', '#9b59b6', '#95a5a6'])\n bounds = [0, 1, 2, 3, 4, 5]\n norm = mcolors.BoundaryNorm(bounds, cmap.N)\n \n # Create figure\n plt.figure(figsize=(20, 10))\n im = plt.imshow(sequences_matrix, cmap=cmap, norm=norm, aspect='auto')\n \n # Add color bar\n cbar = plt.colorbar(im, ticks=[0.5, 1.5, 2.5, 3.5, 4.5])\n cbar.set_label('Nucleotides', fontsize=14)\n cbar.set_ticklabels(['A', 'C', 'G', 'U', 'N/Padding'])\n \n # Add axis labels\n plt.xlabel(\"Position in Sequence\", fontsize=14)\n plt.ylabel(\"RNA Sequences\", fontsize=14)\n \n # Add title\n plt.title(\"RNA Sequences Heatmap\", fontsize=16)\n \n # Add sequence IDs as y-axis labels\n plt.yticks(range(len(valid_ids)), valid_ids, fontsize=10)\n \n # Show only some labels on x-axis to avoid crowding\n sequence_length = sequences_matrix.shape[1]\n step = max(1, sequence_length // 20) # Show at most 20 labels\n plt.xticks(range(0, sequence_length, step), range(1, sequence_length + 1, step))\n \n # Add grid\n plt.grid(False)\n \n # Add information about nucleotide distribution\n all_nucleotides = sequences_matrix.flatten()\n nucleotide_counts = {\n 'A': np.sum(all_nucleotides == 0),\n 'C': np.sum(all_nucleotides == 1),\n 'G': np.sum(all_nucleotides == 2),\n 'U': np.sum(all_nucleotides == 3),\n 'N': np.sum(all_nucleotides == 4)\n }\n \n total_nucleotides = sum(nucleotide_counts.values())\n nucleotide_percentages = {k: (v / total_nucleotides) * 100 for k, v in nucleotide_counts.items()}\n \n # Add text with statistics\n info_text = \"\\n\".join([\n f\"Total sequences visualized: {num_samples}\",\n f\"Maximum length: {sequence_length}\",\n f\"A: {nucleotide_percentages['A']:.1f}%\",\n f\"C: {nucleotide_percentages['C']:.1f}%\",\n f\"G: {nucleotide_percentages['G']:.1f}%\",\n f\"U: {nucleotide_percentages['U']:.1f}%\",\n f\"N/Padding: {nucleotide_percentages['N']:.1f}%\"\n ])\n \n plt.figtext(0.02, 0.02, info_text, fontsize=10, bbox=dict(facecolor='white', alpha=0.8))\n \n # Show the plot\n plt.tight_layout()\n plt.show()\n \n # Optionally, save the plot\n output_dir = '/kaggle/working/'\n plt.savefig(os.path.join(output_dir, 'rna_heatmap.png'), dpi=300)\n print(f\"Heatmap saved to {os.path.join(output_dir, 'rna_heatmap.png')}\")\n \n return sequences_matrix\n except Exception as e:\n print(f\"Error processing data: {e}\")\n return None\n\n# Use the function (assuming processed_data is available)\nvisualize_rna_heatmap_from_processed_data(processed_data)","metadata":{"execution":{"iopub.status.busy":"2025-04-09T15:13:37.170635Z","iopub.execute_input":"2025-04-09T15:13:37.170843Z","iopub.status.idle":"2025-04-09T15:13:37.954372Z","shell.execute_reply.started":"2025-04-09T15:13:37.170824Z","shell.execute_reply":"2025-04-09T15:13:37.953649Z"},"papermill":{"duration":0.854076,"end_time":"2025-03-26T03:43:07.793328","exception":false,"start_time":"2025-03-26T03:43:06.939252","status":"completed"},"tags":[],"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## 🧬 Advanced RNA Structural Prediction: Graph Neural Networks and Conformational Dynamics 🔬","metadata":{"papermill":{"duration":0.018177,"end_time":"2025-03-26T03:43:07.831159","exception":false,"start_time":"2025-03-26T03:43:07.812982","status":"completed"},"tags":[]}},{"cell_type":"markdown","source":"> ## Utility and Preprocessing Functions","metadata":{"papermill":{"duration":0.018692,"end_time":"2025-03-26T03:43:07.869111","exception":false,"start_time":"2025-03-26T03:43:07.850419","status":"completed"},"tags":[]}},{"cell_type":"code","source":"# File paths\nDATA_DIR = \"/kaggle/input/stanford-rna-3d-folding/\"\nOUTPUT_DIR = \"/kaggle/working/\"\nos.makedirs(OUTPUT_DIR, exist_ok=True)\n\nclass StructureWrapper:\n \"\"\"\n Wrapper for RNA structure arrays that provides both quality attributes\n and compatibility with NumPy array operations.\n \"\"\"\n def __init__(self, structure, quality_score=0.5):\n self.structure = structure\n self.quality = {'quality_score': quality_score}\n # Store shape from the underlying structure for numpy compatibility\n self.shape = structure.shape if hasattr(structure, 'shape') else None\n \n def __getitem__(self, idx):\n return self.structure[idx]\n \n def __len__(self):\n return len(self.structure)\n \n # Implement arithmetic operators for numpy compatibility\n def __sub__(self, other):\n \"\"\"Implement subtraction between structures\"\"\"\n if isinstance(other, StructureWrapper):\n # Subtract the underlying structures\n return StructureWrapper(self.structure - other.structure)\n else:\n # Subtract a scalar or numpy array directly\n return StructureWrapper(self.structure - other)\n \n def __add__(self, other):\n \"\"\"Implement addition between structures\"\"\"\n if isinstance(other, StructureWrapper):\n return StructureWrapper(self.structure + other.structure)\n else:\n return StructureWrapper(self.structure + other)\n \n def __mul__(self, other):\n \"\"\"Implement multiplication between structures\"\"\"\n if isinstance(other, StructureWrapper):\n return StructureWrapper(self.structure * other.structure)\n else:\n return StructureWrapper(self.structure * other)\n \n def __truediv__(self, other):\n \"\"\"Implement division between structures\"\"\"\n if isinstance(other, StructureWrapper):\n return StructureWrapper(self.structure / other.structure)\n else:\n return StructureWrapper(self.structure / other)\n \n def __neg__(self):\n \"\"\"Implement negation\"\"\"\n return StructureWrapper(-self.structure)\n \n def __abs__(self):\n \"\"\"Implement absolute value\"\"\"\n return StructureWrapper(abs(self.structure))\n \n # Implement reverse operations (for scalar op structure)\n def __radd__(self, other):\n return StructureWrapper(other + self.structure)\n \n def __rsub__(self, other):\n return StructureWrapper(other - self.structure)\n \n def __rmul__(self, other):\n return StructureWrapper(other * self.structure)\n \n def __rtruediv__(self, other):\n return StructureWrapper(other / self.structure)\n \n # Implement comparison operators\n def __eq__(self, other):\n if isinstance(other, StructureWrapper):\n return self.structure == other.structure\n else:\n return self.structure == other\n \n def __lt__(self, other):\n if isinstance(other, StructureWrapper):\n return self.structure < other.structure\n else:\n return self.structure < other\n \n def __gt__(self, other):\n if isinstance(other, StructureWrapper):\n return self.structure > other.structure\n else:\n return self.structure > other\n \n def __le__(self, other):\n if isinstance(other, StructureWrapper):\n return self.structure <= other.structure\n else:\n return self.structure <= other\n \n def __ge__(self, other):\n if isinstance(other, StructureWrapper):\n return self.structure >= other.structure\n else:\n return self.structure >= other\n \n # Implement numpy compatibility methods\n def __array__(self):\n \"\"\"Allow numpy to automatically convert to array when needed\"\"\"\n import numpy as np\n return np.array(self.structure)\n \n def sum(self, *args, **kwargs):\n \"\"\"Implement sum method for numpy compatibility\"\"\"\n return self.structure.sum(*args, **kwargs)\n \n def mean(self, *args, **kwargs):\n \"\"\"Implement mean method for numpy compatibility\"\"\"\n return self.structure.mean(*args, **kwargs)\n \n def max(self, *args, **kwargs):\n \"\"\"Implement max method for numpy compatibility\"\"\"\n return self.structure.max(*args, **kwargs)\n \n def min(self, *args, **kwargs):\n \"\"\"Implement min method for numpy compatibility\"\"\"\n return self.structure.min(*args, **kwargs)\n \n def reshape(self, *args, **kwargs):\n \"\"\"Implement reshape method for numpy compatibility\"\"\"\n reshaped = self.structure.reshape(*args, **kwargs)\n return StructureWrapper(reshaped, self.quality.get('quality_score', 0.5))\n \n def transpose(self, *args, **kwargs):\n \"\"\"Implement transpose method for numpy compatibility\"\"\"\n transposed = self.structure.transpose(*args, **kwargs)\n return StructureWrapper(transposed, self.quality.get('quality_score', 0.5))\n \n # String representation\n def __repr__(self):\n return f\"StructureWrapper(shape={self.shape}, quality_score={self.quality.get('quality_score', 0.5)})\"\n\nclass ParameterOptimizer:\n \"\"\"\n Meta-learning system for continuous parameter optimization\n based on historical results.\n \"\"\"\n \n def __init__(self, history_file=None):\n \"\"\"\n Initializes the optimizer, optionally loading previous history.\n \n Parameters:\n -----------\n history_file: str, optional\n Path to a file containing the parameter and result history\n \"\"\"\n self.history = []\n if history_file and os.path.exists(history_file):\n self.load_history(history_file)\n \n # Parameter bounds\n self.param_bounds = {\n 'divisor_mean': (3.0, 4.5),\n 'divisor_std': (0.5, 1.5),\n 'noise_base': (0.01, 0.3),\n 'correlation': (0.7, 0.95)\n }\n \n def load_history(self, filename):\n \"\"\"Loads previous parameter and result history\"\"\"\n try:\n with open(filename, 'r') as f:\n self.history = json.load(f)\n except Exception as e:\n print(f\"Error loading history: {str(e)}\")\n \n def save_history(self, filename):\n \"\"\"Saves the current history to a file\"\"\"\n with open(filename, 'w') as f:\n json.dump(self.history, f)\n \n def record_result(self, params, size_category, mode, quality_score):\n \"\"\"\n Records a new result into the history\n \n Parameters:\n -----------\n params: dict\n Parameters used\n size_category: str\n Size category ('small', 'medium', 'large')\n mode: str\n Mode used ('adaptive' or 'fixed')\n quality_score: float\n Quality score obtained\n \"\"\"\n self.history.append({\n 'params': params,\n 'size_category': size_category,\n 'mode': mode,\n 'quality_score': quality_score,\n 'timestamp': datetime.datetime.now().isoformat()\n })\n \n def suggest_parameters(self, size_category, mode):\n \"\"\"\n Suggests optimized parameters based on the history\n for a given size category and mode\n \n Parameters:\n -----------\n size_category: str\n Size category ('small', 'medium', 'large')\n mode: str\n Operation mode ('adaptive' or 'fixed')\n \n Returns:\n --------\n dict: Suggested parameters\n \"\"\"\n # Filter history for the specified category and mode\n relevant_history = [\n entry for entry in self.history \n if entry['size_category'] == size_category and entry['mode'] == mode\n ]\n \n if len(relevant_history) < 5:\n # Not enough history, use default values\n return self._get_default_params(size_category, mode)\n \n # Sort by quality score, from best to worst\n relevant_history.sort(key=lambda x: x['quality_score'], reverse=True)\n \n # Extract parameters from the top N results\n top_n = min(5, len(relevant_history))\n top_params = [entry['params'] for entry in relevant_history[:top_n]]\n \n # Compute weighted average of parameters\n weights = [0.4, 0.25, 0.15, 0.1, 0.1][:top_n] # Weights for top N results\n \n suggested_params = {}\n for param in self.param_bounds.keys():\n if all(param in p for p in top_params):\n weighted_sum = sum(w * p[param] for w, p in zip(weights, top_params))\n suggested_params[param] = weighted_sum / sum(weights[:top_n])\n \n # Ensure suggested parameters are within bounds\n for param, (min_val, max_val) in self.param_bounds.items():\n if param in suggested_params:\n suggested_params[param] = max(min_val, min(suggested_params[param], max_val))\n \n return suggested_params\n \n def _get_default_params(self, size_category, mode):\n \"\"\"Returns default parameters when historical data is insufficient\"\"\"\n # Default values for different size categories and modes\n defaults = {\n 'small': {\n 'adaptive': {'divisor_mean': 3.6, 'divisor_std': 0.9, 'noise_base': 0.15},\n 'fixed': {'noise_base': 0.12, 'correlation': 0.85}\n },\n 'medium': {\n 'adaptive': {'divisor_mean': 3.8, 'divisor_std': 1.0, 'noise_base': 0.1},\n 'fixed': {'noise_base': 0.08, 'correlation': 0.85}\n },\n 'large': {\n 'adaptive': {'divisor_mean': 4.0, 'divisor_std': 1.1, 'noise_base': 0.05},\n 'fixed': {'noise_base': 0.04, 'correlation': 0.9}\n }\n }\n \n return defaults.get(size_category, {}).get(mode, {})\n\ndef normalize_structure(coords):\n \"\"\"\n Centralizes and normalizes the structure.\n \"\"\"\n # Remove padding\n valid_mask = ~np.all(coords == 0, axis=1)\n valid_coords = coords[valid_mask]\n \n # Center at center of mass\n center = np.mean(valid_coords, axis=0)\n centered_coords = coords.copy()\n centered_coords[valid_mask] = valid_coords - center\n \n return centered_coords\n\ndef normalize_coordinates(coords):\n \"\"\"\n Normalizes 3D coordinates of RNA structures by centering and \n scaling each structure independently, with robust handling\n to avoid numerical issues.\n \n Parameters:\n -----------\n coords: Numpy array with shape (batch_size, seq_length, 3)\n 3D coordinates to normalize\n \n Returns:\n --------\n normalized: Numpy array with shape (batch_size, seq_length, 3)\n Normalized coordinates in the range [-1, 1] \n \"\"\"\n # Create copy to avoid modifying the original\n normalized = np.copy(coords)\n \n # Check for problematic values upfront\n if np.isnan(coords).any():\n print(\"WARNING: NaN values detected in input coordinates. They will be ignored during normalization.\")\n if np.isinf(coords).any():\n print(\"WARNING: Infinite values detected in input coordinates. They will be ignored during normalization.\")\n \n # Handle each structure in the batch separately\n for i in range(coords.shape[0]):\n # Identify valid positions (non-zero, non-NaN, non-Inf)\n valid_mask = ~np.all(coords[i] == 0, axis=-1) \n valid_mask = valid_mask & ~np.any(np.isnan(coords[i]), axis=-1)\n valid_mask = valid_mask & ~np.any(np.isinf(coords[i]), axis=-1)\n \n # Extract only valid coordinates\n valid_coords = coords[i][valid_mask]\n \n if len(valid_coords) > 0:\n try:\n # 1. Center at the geometric center\n center = np.nanmean(valid_coords, axis=0)\n \n # Check if the calculated center contains valid values \n if np.isnan(center).any() or np.isinf(center).any():\n print(f\"WARNING: Invalid center calculated for structure {i}. Using [0,0,0].\")\n center = np.zeros(3)\n \n # Apply translation to the center\n centered = valid_coords - center\n \n # 2. Determine appropriate scale factor\n # Calculate maximum distance from the center\n dist_from_center = np.sqrt(np.sum(centered * centered, axis=1))\n \n # Exclude NaN or infinite values for scale_factor calculation\n valid_dists = dist_from_center[~np.isnan(dist_from_center) & ~np.isinf(dist_from_center)]\n \n if len(valid_dists) > 0:\n scale_factor = np.max(valid_dists)\n # Protect against very small scale_factor\n if scale_factor < 1e-10:\n scale_factor = 1.0\n else:\n scale_factor = 1.0\n \n # 3. Normalize coordinates to [-1, 1] range\n normalized_valid = centered / scale_factor\n \n # 4. Replace values in the normalized array\n normalized[i][valid_mask] = normalized_valid\n \n # Debug info\n # print(f\"Structure {i}: center={center}, scale_factor={scale_factor}, \" \n # f\"min={np.min(normalized_valid)}, max={np.max(normalized_valid)}\")\n \n except Exception as e:\n print(f\"ERROR during normalization of structure {i}: {str(e)}\")\n print(\"Keeping original values for this structure.\")\n else:\n print(f\"WARNING: No valid coordinates found for structure {i}.\")\n \n # Final check to detect any issues\n if np.isnan(normalized).any():\n print(\"WARNING: NaN values present after normalization. Replacing with zeros.\")\n normalized = np.nan_to_num(normalized, nan=0.0)\n \n if np.isinf(normalized).any():\n print(\"WARNING: Infinite values present after normalization. Replacing with zeros.\") \n normalized = np.nan_to_num(normalized, posinf=0.0, neginf=0.0)\n \n return normalized\n\ndef check_structure_validity(coords, min_distance=0.8, max_distance=7.0, allow_clashes=0.05):\n \"\"\"\n More refined and realistic biophysical validation.\n \"\"\"\n valid = True\n valid_mask = ~np.all(coords == 0, axis=1)\n valid_coords = coords[valid_mask]\n \n if len(valid_coords) < 3:\n return True\n \n # Check distances between consecutive residues\n invalid_bonds = 0\n for i in range(1, len(valid_coords)):\n dist = np.linalg.norm(valid_coords[i] - valid_coords[i-1])\n if dist < min_distance or dist > max_distance:\n invalid_bonds += 1\n \n # Allow a small percentage of invalid bonds\n if invalid_bonds / len(valid_coords) > 0.1: # More than 10% invalid bonds\n valid = False\n \n # Check for clashes, allowing some\n clashes = 0\n total_pairs = 0\n for i in range(len(valid_coords)):\n for j in range(i+3, len(valid_coords)): # Skip adjacent\n total_pairs += 1\n dist = np.linalg.norm(valid_coords[i] - valid_coords[j])\n if dist < min_distance:\n clashes += 1\n \n # Allow a small percentage of clashes\n if total_pairs > 0 and clashes / total_pairs > allow_clashes:\n valid = False\n \n return valid\n\ndef refine_with_distance_geometry(initial_coords, target_distances, weights, max_iterations=200):\n \"\"\"\n Optimizes 3D coordinates to better satisfy distance constraints.\n\n Parameters:\n -----------\n initial_coords: Initial coordinates (array of shape (n, 3))\n target_distances: Target distance matrix (array of shape (n, n)) \n weights: Matrix of weights for each constraint (array of shape (n, n))\n max_iterations: Maximum number of iterations\n\n Returns:\n --------\n Refined coordinates (array of shape (n, 3))\n \"\"\"\n coords = initial_coords.copy()\n n = coords.shape[0]\n learning_rate = 0.01\n\n for iteration in range(max_iterations):\n # Calculate current distance matrix\n current_distances = np.zeros((n, n))\n for i in range(n):\n for j in range(i+1, n):\n dist = np.linalg.norm(coords[i] - coords[j])\n current_distances[i, j] = dist\n current_distances[j, i] = dist\n\n # Calculate gradients \n grad = np.zeros_like(coords)\n for i in range(n):\n for j in range(i+1, n):\n if weights[i, j] > 0:\n # Vector from i to j\n direction = coords[j] - coords[i]\n current_dist = np.linalg.norm(direction)\n \n # Avoid division by zero\n if current_dist < 1e-10:\n continue\n\n direction = direction / current_dist\n \n # Difference between current and target distance\n diff = current_distances[i, j] - target_distances[i, j]\n \n # Update gradients\n grad_ij = weights[i, j] * diff * direction\n grad[i] += grad_ij\n grad[j] -= grad_ij\n \n # Update coordinates\n coords = coords - learning_rate * grad\n \n # Gradually reduce learning rate\n learning_rate *= 0.995\n\n return coords","metadata":{"execution":{"iopub.status.busy":"2025-04-09T15:13:37.956754Z","iopub.execute_input":"2025-04-09T15:13:37.957035Z","iopub.status.idle":"2025-04-09T15:13:37.989005Z","shell.execute_reply.started":"2025-04-09T15:13:37.957011Z","shell.execute_reply":"2025-04-09T15:13:37.988208Z"},"papermill":{"duration":0.033629,"end_time":"2025-03-26T03:43:07.921623","exception":false,"start_time":"2025-03-26T03:43:07.887994","status":"completed"},"tags":[],"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"def reference_based_approach(X_valid, y_valid, geometric_sampling=True, noise_level=0.21, correlation=0.83):\n \"\"\"\n Placeholder for reference-based RNA structure prediction approach.\n \n Parameters:\n -----------\n X_valid : array-like\n Validation input features\n y_valid : array-like\n Validation target structures\n geometric_sampling : bool, optional\n Whether to use geometric sampling (default True)\n noise_level : float, optional\n Level of noise to add to the structure (default 0.21)\n correlation : float, optional\n Correlation parameter for structure generation (default 0.83)\n \n Returns:\n --------\n model : object\n A placeholder model object with a predict method\n \"\"\"\n class ReferenceModel:\n def __init__(self, noise_level, correlation):\n self.noise_level = noise_level\n self.correlation = correlation\n \n def predict(self, X):\n \"\"\"\n Generate placeholder predictions based on input features.\n \n Parameters:\n -----------\n X : array-like\n Input features for prediction\n \n Returns:\n --------\n predictions : numpy.ndarray\n Generated 3D coordinates\n \"\"\"\n # Create placeholder predictions with some randomness\n predictions = []\n for seq_features in X:\n # Generate a simple 3D structure \n # Assume sequence length based on input features\n seq_length = seq_features.shape[0]\n \n # Create a simple linear structure with some noise\n structure = np.zeros((seq_length, 3))\n for i in range(1, seq_length):\n # Simple linear progression with small random variations\n structure[i] = structure[i-1] + np.array([3.8, 0, 0]) + \\\n np.random.normal(0, self.noise_level, 3)\n \n predictions.append(structure)\n \n return np.array(predictions)\n \n # Create and return the reference model\n return ReferenceModel(noise_level, correlation)\n\ndef sample_structural_variation(coords, noise_level=0.5, preserve_distance=True, \n use_global_movement=False, correlation=0.7):\n \"\"\"\n Enhanced version of structural variation sampling with better\n handling of large RNAs and improved noise distribution.\n \"\"\"\n new_coords = coords.copy()\n valid_mask = ~np.all(coords == 0, axis=1)\n valid_indices = np.where(valid_mask)[0]\n \n if len(valid_indices) < 3:\n return new_coords\n \n # Parameters optimized for RNA structure\n typical_bond_length = 3.8 # Angstroms - typical RNA backbone distance\n \n # Add global domain movements if requested\n if use_global_movement and len(valid_indices) > 20:\n # More natural domain identification - try to find natural hinge points\n # For RNA, these often occur at junctions between helices\n \n # Calculate distance between consecutive residues as a heuristic\n # for finding potential hinge points (larger distances often indicate junctions)\n distances = []\n for i in range(1, len(valid_indices)):\n idx1 = valid_indices[i-1]\n idx2 = valid_indices[i]\n dist = np.linalg.norm(coords[idx1] - coords[idx2])\n distances.append((i, dist))\n \n # Sort by distance to find potential hinges\n distances.sort(key=lambda x: x[1], reverse=True)\n \n # Take top 2 potential hinge points (if we have enough points)\n num_hinges = min(2, len(distances)//3)\n \n for h in range(num_hinges):\n if h < len(distances):\n hinge_point = distances[h][0]\n if hinge_point < 5 or hinge_point > len(valid_indices) - 5:\n continue\n \n hinge_idx = valid_indices[hinge_point]\n \n # Angle of rotation with natural distribution\n # More small movements than large ones\n angle = np.random.exponential(0.2) # Mostly small angles with occasional larger ones\n if np.random.random() < 0.5:\n angle = -angle # Allow both directions\n \n # Create a more natural rotation matrix with slight 3D component\n # RNAs often bend and twist in 3D\n sin_a, cos_a = np.sin(angle), np.cos(angle)\n tilt = np.random.normal(0, 0.1) # Small tilt in 3D\n rotation_matrix = np.array([\n [cos_a, -sin_a, 0],\n [sin_a, cos_a, tilt],\n [0, -tilt, 1]\n ])\n \n # Apply rotation around hinge point\n ref_point = new_coords[hinge_idx]\n for i in valid_indices[hinge_point+1:]:\n vector = new_coords[i] - ref_point\n rotated = np.dot(vector, rotation_matrix)\n new_coords[i] = ref_point + rotated\n \n # Propagate variation residue by residue, with correlation\n # RNA has strong local correlations in structure\n prev_noise = np.zeros(3)\n \n correlation = 0.5 # High correlation for smoother variations\n \n for i in range(1, len(coords)):\n if not valid_mask[i] or not valid_mask[i-1]:\n continue\n \n vec = new_coords[i-1] - new_coords[i]\n vec_length = np.linalg.norm(vec)\n \n # Generate correlated noise (smoother transitions)\n new_noise = np.random.normal(0, noise_level, size=3)\n noise_vec = correlation * prev_noise + (1 - correlation) * new_noise\n prev_noise = noise_vec.copy()\n \n noise_norm = np.linalg.norm(noise_vec)\n if noise_norm > 0:\n # Scale noise proportionally\n noise_vec = noise_vec / noise_norm * (noise_level * vec_length)\n \n # Add noise to the direction\n new_vec = vec + noise_vec\n \n # Preserve distance if requested\n if preserve_distance:\n current_length = np.linalg.norm(new_vec)\n if current_length > 0:\n # Allow slight variation in bond length (RNA is not rigid)\n target_length = typical_bond_length * (1 + np.random.normal(0, 0.05))\n new_vec = new_vec / current_length * target_length\n \n new_coords[i] = new_coords[i-1] - new_vec\n \n return new_coords\n\ndef get_rotation_matrix(axis, theta):\n \"\"\"\n Return the rotation matrix for rotation around an arbitrary axis.\n \n Parameters:\n -----------\n axis: Unit vector defining the rotation axis\n theta: Rotation angle in radians\n \n Returns:\n --------\n 3x3 rotation matrix\n \"\"\"\n # Ensure axis is a unit vector\n axis = axis / np.linalg.norm(axis)\n \n a = np.cos(theta / 2.0)\n b, c, d = -axis * np.sin(theta / 2.0)\n \n return np.array([\n [a*a + b*b - c*c - d*d, 2*(b*c - a*d), 2*(b*d + a*c)],\n [2*(b*c + a*d), a*a + c*c - b*b - d*d, 2*(c*d - a*b)],\n [2*(b*d - a*c), 2*(c*d + a*b), a*a + d*d - b*b - c*c]\n ])\n\n# Auxiliary function to calculate the dihedral angle (in degrees)\ndef calculate_dihedral(p0, p1, p2, p3):\n \"\"\"\n Calculates the dihedral angle (in degrees) defined by the points p0, p1, p2, and p3.\n \"\"\"\n b0 = p1 - p0\n b1 = p2 - p1\n b2 = p3 - p2\n\n # Normalize b1 so its length does not influence the calculation\n b1 /= np.linalg.norm(b1) + 1e-8\n\n # Normal vectors to the planes formed by (p0,p1,p2) and (p1,p2,p3)\n v = b0 - np.dot(b0, b1) * b1\n w = b2 - np.dot(b2, b1) * b1\n\n x = np.dot(v, w)\n y = np.dot(np.cross(b1, v), w)\n angle = np.degrees(np.arctan2(y, x))\n return angle\n\n# Auxiliary function to generate a rotation matrix\ndef get_rotation_matrix(axis, theta):\n \"\"\"\n Returns the 3x3 rotation matrix for a rotation of theta radians around the given 'axis'.\n \"\"\"\n a = np.cos(theta / 2.0)\n b, c, d = -axis * np.sin(theta / 2.0)\n aa, bb, cc, dd = a*a, b*b, c*c, d*d\n bc, ad, ac, ab, bd, cd = b*c, a*d, a*c, a*b, b*d, c*d\n return np.array([\n [aa+bb-cc-dd, 2*(bc+ad), 2*(bd-ac)],\n [2*(bc-ad), aa+cc-bb-dd, 2*(cd+ab)],\n [2*(bd+ac), 2*(cd-ab), aa+dd-bb-cc]\n ])\n\n# New function to refine the RNA backbone with dihedral angle adjustment\ndef refine_rna_backbone_with_dihedrals(structure, ideal_dihedral=180.0):\n \"\"\"\n Refines the geometry of the RNA backbone by adjusting the dihedral angles to an ideal value.\n\n Parameters:\n structure: np.array of shape (seq_length, 3) containing the coordinates.\n ideal_dihedral: Ideal angle (in degrees) for the backbone segments (e.g., 180°).\n\n Returns:\n np.array with the refined structure.\n \"\"\"\n refined = structure.copy()\n n = len(refined)\n if n < 4:\n return refined # There are no dihedral angles to correct\n\n for i in range(n - 3):\n p0, p1, p2, p3 = refined[i], refined[i+1], refined[i+2], refined[i+3]\n current_angle = calculate_dihedral(p0, p1, p2, p3)\n # Compute the difference (in radians) between the ideal angle and the current one\n angle_diff = np.radians(ideal_dihedral - current_angle)\n \n # Define the rotation axis as the direction of the segment (p3 - p2)\n axis = p3 - p2\n norm_axis = np.linalg.norm(axis)\n if norm_axis < 1e-6:\n continue\n axis /= norm_axis\n \n # Get the rotation matrix for the correction angle\n R = get_rotation_matrix(axis, angle_diff)\n \n # Apply the rotation to all points starting from p3\n for j in range(i+3, n):\n vec = refined[j] - p2\n refined[j] = p2 + np.dot(vec, R.T)\n \n return refined","metadata":{"execution":{"iopub.status.busy":"2025-04-09T15:13:37.990122Z","iopub.execute_input":"2025-04-09T15:13:37.990372Z","iopub.status.idle":"2025-04-09T15:13:38.011151Z","shell.execute_reply.started":"2025-04-09T15:13:37.990351Z","shell.execute_reply":"2025-04-09T15:13:38.010297Z"},"papermill":{"duration":0.039629,"end_time":"2025-03-26T03:43:07.980725","exception":false,"start_time":"2025-03-26T03:43:07.941096","status":"completed"},"tags":[],"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"def repair_invalid_structure(structure):\n \"\"\"\n Attempt to repair an invalid RNA structure.\n \n Parameters:\n -----------\n structure: Potentially invalid RNA structure\n \n Returns:\n --------\n Repaired structure\n \"\"\"\n # Create a copy to repair\n repaired = structure.copy()\n \n # Check for valid residues\n valid_mask = ~np.all(repaired == 0, axis=1)\n \n # Fix bond lengths\n for i in range(1, len(repaired)):\n if valid_mask[i] and valid_mask[i-1]:\n # Get current bond\n bond_vector = repaired[i] - repaired[i-1]\n bond_length = np.linalg.norm(bond_vector)\n \n # Check if bond is too short or too long\n if bond_length < 1.0 or bond_length > 7.0:\n # Fix bond to ideal length\n ideal_length = 3.8\n if bond_length > 0:\n repaired[i] = repaired[i-1] + (bond_vector / bond_length) * ideal_length\n else:\n # Generate a random direction if bond length is zero\n random_direction = np.random.randn(3)\n random_direction = random_direction / np.linalg.norm(random_direction)\n repaired[i] = repaired[i-1] + random_direction * ideal_length\n \n # Check for clashes (atoms too close to each other)\n for i in range(len(repaired)):\n if valid_mask[i]:\n for j in range(i+3, len(repaired)): # Skip adjacent residues\n if valid_mask[j]:\n # Calculate distance\n distance = np.linalg.norm(repaired[j] - repaired[i])\n \n # If atoms are too close\n if distance < 1.0:\n # Move one atom away slightly in a random direction\n random_direction = np.random.randn(3)\n random_direction = random_direction / np.linalg.norm(random_direction)\n repaired[j] = repaired[i] + random_direction * 4.0 # Place at safe distance\n \n # Final normalization\n repaired = normalize_structure(repaired)\n \n return repaired\n\ndef create_emergency_structure(seq_length):\n \"\"\"\n Create an emergency structure when all else fails.\n Generates a physically plausible RNA structure.\n \n Parameters:\n -----------\n seq_length: Length of the RNA sequence\n \n Returns:\n --------\n Basic RNA structure\n \"\"\"\n # Create a simple linear structure as fallback\n emergency_structure = np.zeros((seq_length, 3))\n \n # Define canonical nucleotide step (3.8Å)\n step = np.array([3.8, 0.0, 0.0])\n \n # Generate a straight chain with some randomness\n for i in range(seq_length):\n if i == 0:\n emergency_structure[i] = np.zeros(3)\n else:\n # Add slight random deviation to prevent perfect linearity\n random_noise = np.random.normal(0, 0.2, 3)\n emergency_structure[i] = emergency_structure[i-1] + step + random_noise\n \n # Add a slight curve to make it more RNA-like\n # Apply a gentle curve in the y-z plane\n for i in range(seq_length):\n angle = i * 0.1 # Gradual rotation\n emergency_structure[i, 1] += 2 * np.sin(angle) # Y-component\n emergency_structure[i, 2] += 2 * np.cos(angle) # Z-component\n \n # Normalize\n emergency_structure = normalize_structure(emergency_structure)\n \n return emergency_structure","metadata":{"execution":{"iopub.status.busy":"2025-04-09T15:13:38.011988Z","iopub.execute_input":"2025-04-09T15:13:38.012228Z","iopub.status.idle":"2025-04-09T15:13:38.025709Z","shell.execute_reply.started":"2025-04-09T15:13:38.012194Z","shell.execute_reply":"2025-04-09T15:13:38.025043Z"},"papermill":{"duration":0.031523,"end_time":"2025-03-26T03:43:08.034120","exception":false,"start_time":"2025-03-26T03:43:08.002597","status":"completed"},"tags":[],"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"def calculate_tm_score(pred_coords, true_coords, d0_scale=1.24):\n \"\"\"\n Calculates a robust approximation of the TM-score between predicted and true coordinates.\n Adds protections against division by zero and NaN.\n \"\"\"\n # Remove padding (rows with zeros) from the true structures\n mask = ~np.all(true_coords == 0, axis=1)\n pred = pred_coords[mask]\n true = true_coords[mask]\n \n L = len(true)\n if L < 3:\n return 0.0\n \n # Define d0 based on L (values adapted for RNA)\n if L >= 30:\n d0 = 0.6 * np.sqrt(L - 0.5) - 2.5\n d0 = max(0.1, d0)\n elif L >= 24:\n d0 = 0.7\n elif L >= 20:\n d0 = 0.6\n elif L >= 16:\n d0 = 0.5\n elif L >= 12:\n d0 = 0.4\n else:\n d0 = 0.3\n \n distances = np.sqrt(np.sum((pred - true) ** 2, axis=1))\n tm_terms = 1.0 / (1.0 + (distances / (d0 + 1e-8)) ** 2)\n tm_score = np.sum(tm_terms) / L\n return float(tm_score)\n\ndef calculate_tm_score_exact(pred_coords, true_coords):\n \"\"\"\n Implementation more closely matching US-align with sequence-independent alignment.\n Includes multiple rotation schemes to find the optimal structural alignment.\n \"\"\"\n # Remove padding\n mask = ~np.all(true_coords == 0, axis=1)\n pred = pred_coords[mask]\n true = true_coords[mask]\n \n Lref = len(true)\n if Lref < 3:\n return 0.0\n \n # Define d0 exactly as in the evaluation formula\n if Lref >= 30:\n d0 = 0.6 * np.sqrt(Lref - 0.5) - 2.5\n elif Lref >= 24:\n d0 = 0.7\n elif Lref >= 20:\n d0 = 0.6\n elif Lref >= 16:\n d0 = 0.5\n elif Lref >= 12:\n d0 = 0.4\n else:\n d0 = 0.3\n \n # Normalize structures\n pred_centered = pred - np.mean(pred, axis=0)\n true_centered = true - np.mean(true, axis=0)\n \n # Try multiple fragment lengths for sequence-independent alignment\n # This mimics US-align's approach to find the best fragment alignment\n best_tm_score = 0.0\n fragment_lengths = [Lref, max(5, Lref//2), max(5, Lref//4)]\n \n for frag_len in fragment_lengths:\n # Try different fragment start positions\n for i in range(0, Lref - frag_len + 1, max(1, frag_len//2)):\n pred_frag = pred_centered[i:i+frag_len]\n \n # Try aligning with different parts of the true structure\n for j in range(0, Lref - frag_len + 1, max(1, frag_len//2)):\n true_frag = true_centered[j:j+frag_len]\n \n # Covariance matrix for optimal rotation\n covariance = np.dot(pred_frag.T, true_frag)\n U, S, Vt = np.linalg.svd(covariance)\n rotation = np.dot(U, Vt)\n \n # Try different rotation schemes - this is the new part\n rotations_to_try = [\n rotation, # Original rotation from SVD\n np.dot(rotation, np.array([[0, 1, 0], [-1, 0, 0], [0, 0, 1]])), # 90 degree Z rotation\n np.dot(rotation, np.array([[-1, 0, 0], [0, -1, 0], [0, 0, 1]])) # 180 degree Z rotation\n ]\n \n for rot in rotations_to_try:\n # Apply rotation to the full structure\n pred_aligned = np.dot(pred_centered, rot)\n \n # Calculate distances\n distances = np.sqrt(np.sum((pred_aligned - true_centered) ** 2, axis=1))\n \n # Calculate TM-score terms\n tm_terms = 1.0 / (1.0 + (distances / d0) ** 2)\n tm_score = np.sum(tm_terms) / Lref\n \n best_tm_score = max(best_tm_score, tm_score)\n \n return float(best_tm_score)\n\ndef load_processed_data():\n \"\"\"\n Loads processed data for training.\n \"\"\"\n X_train = np.load(os.path.join(OUTPUT_DIR, 'X_train.npy'))\n y_train = np.load(os.path.join(OUTPUT_DIR, 'y_train.npy'))\n X_valid = np.load(os.path.join(OUTPUT_DIR, 'X_valid.npy'))\n y_valid = np.load(os.path.join(OUTPUT_DIR, 'y_valid.npy'))\n \n print(f\"Data loaded - X_train: {X_train.shape}, y_train: {y_train.shape}\")\n print(f\"Data loaded - X_valid: {X_valid.shape}, y_valid: {y_valid.shape}\")\n \n return X_train, y_train, X_valid, y_valid","metadata":{"execution":{"iopub.status.busy":"2025-04-09T15:13:38.026575Z","iopub.execute_input":"2025-04-09T15:13:38.026851Z","iopub.status.idle":"2025-04-09T15:13:38.042713Z","shell.execute_reply.started":"2025-04-09T15:13:38.026821Z","shell.execute_reply":"2025-04-09T15:13:38.041943Z"},"papermill":{"duration":0.031127,"end_time":"2025-03-26T03:43:08.084534","exception":false,"start_time":"2025-03-26T03:43:08.053407","status":"completed"},"tags":[],"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"def prepare_test_features(test_seq_df, max_length=720):\n \"\"\"\n Prepares test features (one-hot encoding of the sequence).\n \"\"\"\n X_test = []\n for _, row in test_seq_df.iterrows():\n seq = row['sequence']\n features = []\n for nucleotide in seq:\n if nucleotide == 'A':\n features.append([1, 0, 0, 0, 0])\n elif nucleotide == 'C':\n features.append([0, 1, 0, 0, 0])\n elif nucleotide == 'G':\n features.append([0, 0, 1, 0, 0])\n elif nucleotide == 'U':\n features.append([0, 0, 0, 1, 0])\n else:\n features.append([0, 0, 0, 0, 1])\n if len(features) < max_length:\n padding = [[0, 0, 0, 0, 0]] * (max_length - len(features))\n features.extend(padding)\n else:\n features = features[:max_length]\n X_test.append(features)\n return np.array(X_test)\n\ndef extract_sequence_features(seq_features):\n \"\"\"\n Extract relevant sequence features from one-hot encoding.\n \"\"\"\n # Get valid rows (non-padding)\n valid_mask = ~np.all(seq_features == 0, axis=1)\n valid_features = seq_features[valid_mask]\n \n # Calculate nucleotide composition\n a_content = np.mean(valid_features[:, 0])\n c_content = np.mean(valid_features[:, 1])\n g_content = np.mean(valid_features[:, 2])\n u_content = np.mean(valid_features[:, 3])\n gc_content = c_content + g_content\n \n return {\n 'length': np.sum(valid_mask),\n 'a_content': a_content,\n 'c_content': c_content,\n 'g_content': g_content, \n 'u_content': u_content,\n 'gc_content': gc_content,\n 'au_content': a_content + u_content\n }\n\ndef visualize_3d_structure(true_coords, pred_coords, sample_idx=0, title=\"3D Structure Comparison\", show_plot=False):\n \"\"\"\n Visualizes the true and predicted 3D structures for a sample.\n Only shows the plot if explicitly requested.\n \"\"\"\n true = true_coords[sample_idx]\n pred = pred_coords[sample_idx]\n mask = ~np.all(true == 0, axis=1)\n true = true[mask]\n pred = pred[mask]\n \n fig = plt.figure(figsize=(15, 7))\n ax1 = fig.add_subplot(121, projection='3d')\n ax1.plot(true[:, 0], true[:, 1], true[:, 2], 'b-', label='True')\n ax1.scatter(true[:, 0], true[:, 1], true[:, 2], c='b', s=20, alpha=0.5)\n ax1.set_title('True Structure')\n ax1.set_xlabel('X')\n ax1.set_ylabel('Y')\n ax1.set_zlabel('Z')\n ax1.grid(True)\n \n ax2 = fig.add_subplot(122, projection='3d')\n ax2.plot(pred[:, 0], pred[:, 1], pred[:, 2], 'r-', label='Predicted')\n ax2.scatter(pred[:, 0], pred[:, 1], pred[:, 2], c='r', s=20, alpha=0.5)\n ax2.set_title('Predicted Structure')\n ax2.set_xlabel('X')\n ax2.set_ylabel('Y')\n ax2.set_zlabel('Z')\n ax2.grid(True)\n \n plt.suptitle(title)\n plt.tight_layout()\n \n # Always save the figure\n filename = f'structure_comparison_{sample_idx}.png'\n plt.savefig(os.path.join(OUTPUT_DIR, filename))\n \n # Only show the plot if requested\n if show_plot:\n plt.show()\n else:\n plt.close(fig)\n \n return filename # Return the filename for reference","metadata":{"execution":{"iopub.status.busy":"2025-04-09T15:13:38.043541Z","iopub.execute_input":"2025-04-09T15:13:38.043854Z","iopub.status.idle":"2025-04-09T15:13:38.060579Z","shell.execute_reply.started":"2025-04-09T15:13:38.043823Z","shell.execute_reply":"2025-04-09T15:13:38.059711Z"},"papermill":{"duration":0.029709,"end_time":"2025-03-26T03:43:08.132009","exception":false,"start_time":"2025-03-26T03:43:08.102300","status":"completed"},"tags":[],"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"> ## RNA-specific functions","metadata":{"papermill":{"duration":0.017717,"end_time":"2025-03-26T03:43:08.167465","exception":false,"start_time":"2025-03-26T03:43:08.149748","status":"completed"},"tags":[]}},{"cell_type":"code","source":"def identify_stem_loops(sequence):\n \"\"\"\n Simple function to identify potential stem-loop regions in RNA.\n \n Parameters:\n -----------\n sequence: RNA sequence\n \n Returns:\n --------\n List of (start, end) indices for potential stem loops\n \"\"\"\n # This is a simplified implementation\n # A real implementation would use a more sophisticated algorithm\n \n stem_loops = []\n min_stem_length = 3\n \n # Look for complementary regions that could form stems\n for i in range(len(sequence) - 2*min_stem_length - 3):\n for j in range(i + min_stem_length + 3, len(sequence) - min_stem_length):\n # Check if regions could form a stem\n potential_stem = True\n for k in range(min_stem_length):\n if not are_complementary(sequence[i+k], sequence[j+min_stem_length-1-k]):\n potential_stem = False\n break\n \n if potential_stem:\n # Potential stem-loop found\n stem_loops.append((i, j + min_stem_length))\n break\n \n return stem_loops\n\ndef are_complementary(base1, base2):\n \"\"\"Check if two bases are complementary in RNA.\"\"\"\n return (base1 == 'A' and base2 == 'U') or \\\n (base1 == 'U' and base2 == 'A') or \\\n (base1 == 'G' and base2 == 'C') or \\\n (base1 == 'C' and base2 == 'G') or \\\n (base1 == 'G' and base2 == 'U') or \\\n (base1 == 'U' and base2 == 'G') # G-U wobble pairs are valid in RNA\n\ndef apply_stem_loop_template(structure, start, end):\n \"\"\"\n Apply a stem-loop template to a specific region of the structure.\n \n Parameters:\n -----------\n structure: RNA 3D structure\n start, end: Indices of the stem-loop region\n \n Returns:\n --------\n Modified structure with stem-loop template applied\n \"\"\"\n # Create a copy to modify\n result = structure.copy()\n \n # Length of the region\n region_length = end - start + 1\n \n # Not enough residues to form a proper stem-loop\n if region_length < 7:\n return result\n \n # Calculate stem length (approximately 1/3 of the region on each side)\n stem_length = max(2, region_length // 6)\n loop_start = start + stem_length\n loop_end = end - stem_length\n \n # Loop length\n loop_length = loop_end - loop_start + 1\n \n # Apply stem template (roughly parallel strands)\n for i in range(stem_length):\n # Base positions in the two stems\n pos1 = start + i\n pos2 = end - i\n \n if pos1 < len(result) and pos2 < len(result):\n # Create roughly parallel strands\n if i > 0:\n # Base the position on the previous nucleotide in the strand\n result[pos1] = result[pos1-1] + np.array([0.0, 3.8, 0.0])\n result[pos2] = result[pos2+1] + np.array([0.0, -3.8, 0.0])\n \n # Apply loop template (roughly circular)\n if loop_length > 0:\n # Calculate center of the loop\n if loop_start < len(result) and loop_end < len(result):\n center = (result[loop_start-1] + result[loop_end+1]) / 2\n center[1] += 4.0 # Offset in y direction\n \n # Create a circular loop\n radius = 3.8 # approximately nucleotide distance\n for i in range(loop_length):\n idx = loop_start + i\n if idx < len(result):\n angle = np.pi * i / (loop_length - 1)\n result[idx] = center + np.array([\n radius * np.cos(angle),\n 0.0,\n radius * np.sin(angle)\n ])\n \n return result\n\ndef post_process_rna_structure(structure, sequence, gc_content, use_global_movement=True):\n \"\"\"\n Apply RNA-specific post-processing to refine a structure.\n \n Parameters:\n -----------\n structure: Predicted 3D coordinates\n sequence: RNA sequence\n gc_content: GC content of the sequence\n use_global_movement: Whether to apply global movement transformations\n \n Returns:\n --------\n Refined structure\n \"\"\"\n # Create a new structure for modifications\n result = structure.copy()\n \n # 1. Apply mild refinement based on sequence composition\n noise_level = 0.1\n if gc_content > 0.6:\n # GC-rich regions tend to form more stable structures\n noise_level = 0.05 # Lower noise for more stable structures\n elif gc_content < 0.4:\n # AT-rich regions tend to be more flexible\n noise_level = 0.15 # Higher noise for more flexible regions\n \n # Apply noise proportional to sequence characteristics\n result = sample_structural_variation(\n result,\n noise_level=noise_level,\n preserve_distance=True, # Always preserve distances for realistic structures\n use_global_movement=use_global_movement,\n correlation=0.85 # High correlation for smoother changes\n )\n \n # 2. Look for motifs in the sequence and apply structure templates\n # This is a simplified example - a complete implementation would include more motifs\n stem_loops = identify_stem_loops(sequence)\n if stem_loops:\n for start, end in stem_loops:\n # Apply stem-loop template to these regions\n result = apply_stem_loop_template(result, start, end)\n \n # 3. Normalize bond lengths to ideal values for RNA\n valid_mask = ~np.all(result == 0, axis=1)\n for i in range(1, len(result)):\n if valid_mask[i] and valid_mask[i-1]:\n # Get the current bond vector\n bond_vector = result[i] - result[i-1]\n bond_length = np.linalg.norm(bond_vector)\n \n if bond_length > 0:\n # Normalize to ideal RNA backbone distance with small variation\n ideal_length = 3.8 * (1 + np.random.normal(0, 0.03))\n result[i] = result[i-1] + (bond_vector / bond_length) * ideal_length\n \n return result","metadata":{"execution":{"iopub.status.busy":"2025-04-09T15:13:38.061322Z","iopub.execute_input":"2025-04-09T15:13:38.061613Z","iopub.status.idle":"2025-04-09T15:13:38.075842Z","shell.execute_reply.started":"2025-04-09T15:13:38.061592Z","shell.execute_reply":"2025-04-09T15:13:38.074939Z"},"papermill":{"duration":0.031324,"end_time":"2025-03-26T03:43:08.217678","exception":false,"start_time":"2025-03-26T03:43:08.186354","status":"completed"},"tags":[],"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"> ## Definition of Graph Model Classes and Functions 🧬🔬","metadata":{}},{"cell_type":"code","source":"class GraphAttentionLayer(nn.Module):\n \"\"\"\n Graph attention layer for processing node and edge features\n \"\"\"\n def __init__(self, in_features, out_features, heads=8, dropout=0.1):\n super(GraphAttentionLayer, self).__init__()\n self.in_features = in_features\n self.out_features = out_features\n self.heads = heads\n self.dropout = dropout\n \n # Linear transformations for query, key and value\n self.query = nn.Linear(in_features, out_features * heads)\n self.key = nn.Linear(in_features, out_features * heads)\n self.value = nn.Linear(in_features, out_features * heads)\n \n # Transformation for edge type\n self.edge_attn = nn.Linear(in_features, heads)\n \n # Dropout and output layer\n self.dropout_layer = nn.Dropout(dropout)\n self.output_transform = nn.Linear(out_features * heads, out_features)\n \n def forward(self, x, edge_index, edge_attr=None):\n # Shape: x = [num_nodes, in_features]\n # edge_index = [2, num_edges]\n # edge_attr = [num_edges, edge_features]\n \n # Calculate queries, keys and values\n queries = self.query(x).view(-1, self.heads, self.out_features)\n keys = self.key(x).view(-1, self.heads, self.out_features)\n values = self.value(x).view(-1, self.heads, self.out_features)\n \n # Extract source and destination nodes for each edge\n src, dst = edge_index\n \n # Calculate attention scores\n q_i = queries[dst]\n k_j = keys[src]\n \n # Dot product attention\n alpha = torch.sum(q_i * k_j, dim=-1) / np.sqrt(self.out_features)\n \n # Add attention based on edge type, if available\n if edge_attr is not None:\n edge_attention = self.edge_attn(edge_attr).view(-1, self.heads)\n alpha = alpha + edge_attention\n \n # Softmax to normalize attention between neighbors\n alpha = torch.softmax(alpha, dim=0)\n alpha = self.dropout_layer(alpha)\n \n # Apply attention to values\n v_j = values[src].view(-1, self.heads, self.out_features)\n weighted_values = v_j * alpha.unsqueeze(-1)\n \n # Aggregate weighted values\n output = torch.zeros_like(queries)\n for i in range(dst.max() + 1):\n mask = (dst == i)\n if mask.any():\n output[i] = weighted_values[mask].sum(dim=0)\n \n # Concatenate/transform attention heads for final dimensionality\n output = output.reshape(-1, self.heads * self.out_features)\n output = self.output_transform(output)\n \n return output\n\nclass RNAGraphTransformer(nn.Module):\n def __init__(self, node_features, edge_features, hidden_dim=128, n_layers=6):\n super(RNAGraphTransformer, self).__init__()\n \n # Embedding layers\n self.node_embedding = nn.Linear(node_features, hidden_dim)\n self.edge_embedding = nn.Linear(edge_features, hidden_dim)\n \n # Graph attention layers\n self.graph_layers = nn.ModuleList([\n GraphAttentionLayer(hidden_dim, hidden_dim, heads=8, dropout=0.1)\n for _ in range(n_layers)\n ])\n \n # MLPs for processing by interaction type\n self.covalent_mlp = nn.Sequential(\n nn.Linear(hidden_dim*2, hidden_dim),\n nn.ReLU(),\n nn.Linear(hidden_dim, hidden_dim)\n )\n \n self.basepair_mlp = nn.Sequential(\n nn.Linear(hidden_dim*2, hidden_dim),\n nn.ReLU(),\n nn.Linear(hidden_dim, hidden_dim)\n )\n \n self.tertiary_mlp = nn.Sequential(\n nn.Linear(hidden_dim*2, hidden_dim),\n nn.ReLU(),\n nn.Linear(hidden_dim, hidden_dim)\n )\n \n # Output layer for 3D coordinates\n self.coords_output = nn.Linear(hidden_dim, 3)\n \n # Output layer for distances between residues\n self.distance_output = nn.Linear(hidden_dim*2, 1)\n \n def forward(self, G):\n # Get node and edge features\n node_features = torch.stack([G.nodes[i]['features'] for i in G.nodes])\n edge_index, edge_features = self._get_edge_data(G)\n \n # Initial embedding\n h_nodes = self.node_embedding(node_features)\n \n # Pass through graph layers\n for layer in self.graph_layers:\n h_nodes = layer(h_nodes, edge_index, edge_features)\n \n # Predict coordinates for each node\n coords = self.coords_output(h_nodes)\n \n # Predict distances between connected pairs\n distances = {}\n for (i, j) in G.edges():\n # Concatenate features of the two nodes\n edge_repr = torch.cat([h_nodes[i], h_nodes[j]], dim=0)\n \n # Process with specific MLP based on edge type\n edge_type = G[i][j]['type']\n if edge_type == 'covalent':\n edge_repr = self.covalent_mlp(edge_repr)\n elif edge_type == 'basepair':\n edge_repr = self.basepair_mlp(edge_repr)\n elif edge_type == 'tertiary':\n edge_repr = self.tertiary_mlp(edge_repr)\n \n # Predict distance\n distances[(i, j)] = self.distance_output(edge_repr)\n \n return coords, distances\n \n def _get_edge_data(self, G):\n # Convert edge data to the format expected by PyTorch Geometric\n edge_indices = []\n edge_features_list = []\n \n for i, j in G.edges():\n edge_indices.append([i, j])\n \n # Convert edge attributes to feature vector\n edge_attr = G[i][j]\n edge_type = edge_attr['type']\n \n # One-hot encoding of edge type\n edge_type_vector = [0, 0, 0, 0] # [covalent, basepair, stacking, tertiary]\n if edge_type == 'covalent':\n edge_type_vector[0] = 1\n elif edge_type == 'basepair':\n edge_type_vector[1] = 1\n elif edge_type == 'stacking':\n edge_type_vector[2] = 1\n elif edge_type == 'tertiary':\n edge_type_vector[3] = 1\n \n # Add edge weight\n weight = edge_attr.get('weight', 1.0)\n \n # Combine into a single feature vector\n combined_features = edge_type_vector + [weight]\n edge_features_list.append(combined_features)\n \n # Convert to PyTorch tensors\n edge_index = torch.tensor(edge_indices).t().contiguous() # Transpose to shape [2, num_edges]\n edge_features = torch.tensor(edge_features_list, dtype=torch.float)\n \n return edge_index, edge_features\n\n# ==== 2. FUNCTIONS FOR RNA GRAPH PROCESSING ====\n\ndef get_base_features(base):\n \"\"\"\n Converts a nucleotide into a feature vector.\n \"\"\"\n # One-hot encoding for base type\n base_encoding = {\n 'A': [1, 0, 0, 0, 0], # Adenine\n 'C': [0, 1, 0, 0, 0], # Cytosine\n 'G': [0, 0, 1, 0, 0], # Guanine\n 'U': [0, 0, 0, 1, 0], # Uracil\n 'T': [0, 0, 0, 1, 0], # Treats T as U\n 'N': [0, 0, 0, 0, 1] # Unknown base\n }\n \n # Basic base features\n encoding = base_encoding.get(base, [0, 0, 0, 0, 1]) # Default to N if not recognized\n \n # Additional properties\n is_purine = 1.0 if base in ['A', 'G'] else 0.0\n is_pyrimidine = 1.0 if base in ['C', 'U', 'T'] else 0.0\n \n # Pairing features\n can_pair_with_A = 1.0 if base in ['U', 'T'] else 0.0\n can_pair_with_C = 1.0 if base in ['G'] else 0.0\n can_pair_with_G = 1.0 if base in ['C', 'U', 'T'] else 0.0 # G-U wobble\n can_pair_with_U = 1.0 if base in ['A', 'G'] else 0.0 # G-U wobble\n \n # Combine all features\n features = encoding + [is_purine, is_pyrimidine, \n can_pair_with_A, can_pair_with_C, \n can_pair_with_G, can_pair_with_U]\n \n return torch.tensor(features, dtype=torch.float)\n\ndef predict_pair_type(base1, base2):\n \"\"\"\n Determines the type of base pair between two nucleotides.\n \"\"\"\n if (base1 == 'A' and base2 == 'U') or (base1 == 'U' and base2 == 'A'):\n return 'AU'\n elif (base1 == 'G' and base2 == 'C') or (base1 == 'C' and base2 == 'G'):\n return 'GC'\n elif (base1 == 'G' and base2 == 'U') or (base1 == 'U' and base2 == 'G'):\n return 'GU' # wobble pair\n else:\n return 'noncanonical'\n\ndef predict_tertiary_interactions(sequence):\n \"\"\"\n Predicts possible tertiary interactions based on the sequence.\n This is a simplified implementation and should be replaced by a trained model.\n \"\"\"\n tertiary_interactions = []\n \n # Simple rules for potential tertiary interactions\n for i in range(len(sequence)):\n for j in range(i + 4, len(sequence)): # At least 4 bases apart\n base_i = sequence[i]\n base_j = sequence[j]\n \n # Heuristic rules for tertiary interactions\n if (base_i == 'A' and base_j == 'G') or (base_i == 'G' and base_j == 'A'):\n interaction_type = 'A-minor'\n prob = 0.3\n tertiary_interactions.append((i, j, interaction_type, prob))\n elif (base_i == 'G' and base_j == 'G'):\n interaction_type = 'G-quadruplex'\n prob = 0.2\n tertiary_interactions.append((i, j, interaction_type, prob))\n \n return tertiary_interactions\n\ndef create_rna_graph(sequence, predicted_contacts=None):\n \"\"\"\n Creates a graph representing the RNA molecule.\n \"\"\"\n G = nx.Graph()\n \n # 1. Add nodes (nucleotides)\n for i, base in enumerate(sequence):\n G.add_node(i, base=base, position=None, features=get_base_features(base))\n \n # 2. Add covalent bond edges (backbone)\n for i in range(len(sequence)-1):\n G.add_edge(i, i+1, type='covalent', weight=1.0)\n \n # 3. Add base pairing edges (if predicted)\n if predicted_contacts is not None:\n for i, j, prob in predicted_contacts:\n if i < j-3: # Avoid trivial contacts\n G.add_edge(i, j, type='basepair', weight=prob, \n pair_type=predict_pair_type(sequence[i], sequence[j]))\n \n # 4. Add possible stacking interactions\n for i in range(len(sequence)-1):\n G.add_edge(i, i+1, type='stacking', weight=0.8) # Adjacent stacking\n \n # 5. Add possible tertiary interactions\n tertiary_interactions = predict_tertiary_interactions(sequence)\n for i, j, interaction_type, prob in tertiary_interactions:\n G.add_edge(i, j, type='tertiary', interaction=interaction_type, weight=prob)\n \n return G\n\n# ==== 3. GEOMETRY AND REFINEMENT FUNCTIONS ====\n\ndef compute_distance_matrix(coords):\n \"\"\"\n Calculates the distance matrix from 3D coordinates.\n \"\"\"\n n = coords.shape[0]\n dist_matrix = np.zeros((n, n))\n \n for i in range(n):\n for j in range(i+1, n):\n dist = np.linalg.norm(coords[i] - coords[j])\n dist_matrix[i, j] = dist\n dist_matrix[j, i] = dist\n \n return dist_matrix\n\ndef distance_geometry_optimization(initial_coords, target_distances, weights, max_iterations=200):\n \"\"\"\n Optimizes 3D coordinates to better satisfy distance constraints.\n \"\"\"\n coords = initial_coords.copy()\n n = coords.shape[0]\n learning_rate = 0.01\n \n for iteration in range(max_iterations):\n # Calculate current distance matrix\n current_distances = compute_distance_matrix(coords)\n \n # Calculate gradients\n grad = np.zeros_like(coords)\n for i in range(n):\n for j in range(i+1, n):\n if weights[i, j] > 0:\n # Unit vector from i to j\n direction = coords[j] - coords[i]\n current_dist = np.linalg.norm(direction)\n \n # Avoid division by zero\n if current_dist < 1e-10:\n continue\n \n direction = direction / current_dist\n \n # Difference between current and target distance\n diff = current_distances[i, j] - target_distances[i, j]\n \n # Update gradients\n grad_ij = weights[i, j] * diff * direction\n grad[i] += grad_ij\n grad[j] -= grad_ij\n \n # Update coordinates\n coords = coords - learning_rate * grad\n \n # Gradually reduce learning rate\n learning_rate *= 0.995\n \n return coords\n\ndef weighted_coordinate_average(coord_weight_pairs):\n \"\"\"\n Calculates a weighted average of coordinate sets.\n \"\"\"\n total_weight = sum(weight for _, weight in coord_weight_pairs)\n avg_coords = np.zeros_like(coord_weight_pairs[0][0])\n \n for coords, weight in coord_weight_pairs:\n avg_coords += (coords * weight / total_weight)\n \n return avg_coords\n\ndef load_fragment_library():\n \"\"\"\n Loads library of RNA structural fragments.\n This is a simplified implementation and should be replaced by a real library.\n \"\"\"\n # Create a simple example library\n fragment_library = {\n 'stem': [np.random.randn(10, 3) for _ in range(5)], # 5 examples of stems with 10 nucleotides\n 'loop': [np.random.randn(5, 3) for _ in range(3)], # 3 examples of loops with 5 nucleotides\n 'bulge': [np.random.randn(3, 3) for _ in range(2)], # 2 examples of bulges with 3 nucleotides\n 'junction': [np.random.randn(8, 3) for _ in range(2)], # 2 examples of junctions with 8 nucleotides\n }\n \n print(\"Loaded simplified structural fragment library\")\n return fragment_library\n\ndef predict_secondary_structure(sequence):\n \"\"\"\n Predicts RNA secondary structure from sequence.\n This is a simplified implementation.\n \"\"\"\n # Example implementation - in a real system, you would use\n # secondary structure prediction methods like ViennaRNA\n structure_elements = []\n \n # Simplification: treat everything as stem or loop\n i = 0\n while i < len(sequence):\n if i < len(sequence) - 10:\n # Check for possible stem\n stem_length = min(5, (len(sequence) - i) // 2)\n structure_elements.append({\n 'id': len(structure_elements),\n 'type': 'stem',\n 'start': i,\n 'end': i + 2*stem_length - 1,\n 'sequence': sequence[i:i + 2*stem_length]\n })\n i += 2*stem_length\n else:\n # Remainder as loop\n structure_elements.append({\n 'id': len(structure_elements),\n 'type': 'loop',\n 'start': i,\n 'end': len(sequence) - 1,\n 'sequence': sequence[i:]\n })\n i = len(sequence)\n \n return structure_elements\n\ndef find_best_fragment(sequence, fragment_library):\n \"\"\"\n Finds the best fragment in the library for a given sequence.\n \"\"\"\n # Simple implementation - just returns the first fragment\n # In a real system, you would do a comparison based on sequence/geometry\n if len(fragment_library) > 0:\n return fragment_library[0]\n return np.zeros((len(sequence), 3)) # Return zeros if library is empty\n\ndef assemble_fragments(selected_fragments, structure_elements):\n \"\"\"\n Assembles an initial structure from selected fragments.\n \"\"\"\n # Determine total sequence length\n max_pos = max([elem['end'] for elem in structure_elements]) + 1\n \n # Initialize structure\n assembled_structure = np.zeros((max_pos, 3))\n \n # Current position for assembly\n current_pos = np.zeros(3)\n \n # Assemble each structural element\n for elem in sorted(structure_elements, key=lambda x: x['start']):\n elem_id = elem['id']\n fragment = selected_fragments[elem_id]\n \n # Place fragment at current position\n length = elem['end'] - elem['start'] + 1\n fragment_resized = fragment\n \n # Resize or truncate fragment if necessary\n if len(fragment) != length:\n if len(fragment) > length:\n fragment_resized = fragment[:length]\n else:\n # Extend fragment by repeating last coordinate\n fragment_resized = np.vstack([fragment, np.tile(fragment[-1], (length - len(fragment), 1))])\n \n # Position in 3D space (simplified)\n fragment_centered = fragment_resized - fragment_resized[0] + current_pos\n \n # Add to assembled structure\n assembled_structure[elem['start']:elem['end']+1] = fragment_centered\n \n # Update current position\n current_pos = fragment_centered[-1] + np.array([3.8, 0, 0]) # Approximate bond distance\n \n return assembled_structure\n\ndef generate_model_ensemble(structure, sequence, num_models=5, perturbation_scale=0.2):\n \"\"\"\n Generates an ensemble of structural models based on controlled perturbations.\n \"\"\"\n ensemble = [structure] # Include original model\n \n for i in range(1, num_models):\n # Apply perturbation with gradually increasing scale\n scale = perturbation_scale * i / num_models\n \n # Use existing structural variation function\n perturbed = sample_structural_variation(\n structure,\n noise_level=scale,\n preserve_distance=True,\n use_global_movement=(i % 2 == 0),\n correlation=0.9 - (i * 0.1 / num_models)\n )\n \n # Normalize and add to ensemble\n perturbed = normalize_structure(perturbed)\n ensemble.append(perturbed)\n \n return ensemble[:num_models]\n\n# ==== 4. MAIN PREDICTION PIPELINE ====\n\ndef predict_rna_contacts(sequence):\n \"\"\"\n Predicts contacts to form base pairs in the RNA molecule.\n Simplified rule-based implementation; in production, replace with trained model.\n \"\"\"\n contacts = []\n \n # Simple implementation based on Watson-Crick pairing rules\n for i in range(len(sequence)):\n for j in range(i+4, len(sequence)): # At least 4 bases separation\n base_i = sequence[i]\n base_j = sequence[j]\n \n # Check pairing compatibility\n if (base_i == 'A' and base_j == 'U') or (base_i == 'U' and base_j == 'A'):\n contacts.append((i, j, 0.95)) # High confidence\n elif (base_i == 'G' and base_j == 'C') or (base_i == 'C' and base_j == 'G'):\n contacts.append((i, j, 0.98)) # Highest confidence\n elif (base_i == 'G' and base_j == 'U') or (base_i == 'U' and base_j == 'G'):\n contacts.append((i, j, 0.85)) # Wobble pair, lower confidence\n \n # Filter redundant or mutually exclusive contacts\n filtered_contacts = []\n used_positions = set()\n \n for i, j, prob in sorted(contacts, key=lambda x: x[2], reverse=True):\n if i not in used_positions and j not in used_positions:\n filtered_contacts.append((i, j, prob))\n used_positions.add(i)\n used_positions.add(j)\n \n return filtered_contacts\n\ndef refine_with_constraints(structure, predicted_contacts, physical_constraints):\n \"\"\"\n Refines the structure to satisfy contact and physical constraints.\n \"\"\"\n # Simplified implementation\n refined = structure.copy()\n \n # Apply distance constraints for predicted contacts\n for i, j, prob in predicted_contacts:\n if i >= len(structure) or j >= len(structure):\n continue\n \n # Current coordinates\n pos_i = structure[i]\n pos_j = structure[j]\n \n # Current distance\n current_dist = np.linalg.norm(pos_i - pos_j)\n \n # Target distance for a base pair (~5-6Å)\n target_dist = 5.5\n \n # Move nucleotides closer or further as needed\n if current_dist > 0:\n direction = (pos_j - pos_i) / current_dist\n adjustment = (current_dist - target_dist) * prob * 0.5\n \n refined[i] = pos_i + direction * adjustment\n refined[j] = pos_j - direction * adjustment\n \n # Apply constraints for backbone bond lengths\n bond_length = physical_constraints.get('bond_lengths', {}).get('backbone', 3.8)\n \n for i in range(1, len(structure)):\n pos_prev = refined[i-1]\n pos_curr = refined[i]\n \n current_bond = np.linalg.norm(pos_curr - pos_prev)\n \n if current_bond > 0:\n direction = (pos_curr - pos_prev) / current_bond\n refined[i] = pos_prev + direction * bond_length\n \n return refined\n\ndef geometric_assembly(sequence, predicted_contacts, fragment_library):\n \"\"\"\n Assembles a 3D structure using fragments and satisfying physical constraints.\n \"\"\"\n # 1. Decomposition of sequence into likely structural elements\n structure_elements = predict_secondary_structure(sequence)\n \n # 2. Select appropriate fragments from the library\n selected_fragments = {}\n for element in structure_elements:\n element_type = element['type']\n element_seq = element['sequence']\n \n # Find the most compatible fragment\n if element_type in fragment_library:\n best_fragment = find_best_fragment(element_seq, fragment_library[element_type])\n selected_fragments[element['id']] = best_fragment\n \n # 3. Initial assembly with geometric superposition\n initial_structure = assemble_fragments(selected_fragments, structure_elements)\n \n # 4. Refinement to satisfy physical constraints\n refined_structure = refine_with_constraints(\n initial_structure, \n predicted_contacts,\n physical_constraints={\n 'bond_lengths': {'backbone': 3.8},\n 'bond_angles': {'backbone': 110},\n 'stacking_geometry': {'distance': 3.4},\n 'nonbonded_distances': {'min': 3.0}\n }\n )\n \n return refined_structure\n\ndef energy_based_refinement(structure, sequence):\n \"\"\"\n Simplified version that uses the existing function for refinement.\n \"\"\"\n # Extract GC content to use with the post_process_rna_structure function\n gc_count = sum(1 for base in sequence if base in ['G', 'C'])\n gc_content = gc_count / len(sequence) if len(sequence) > 0 else 0.5\n \n # Use existing functions from original code\n refined = refine_rna_backbone(structure)\n refined = post_process_rna_structure(refined, sequence, gc_content, use_global_movement=True)\n \n return refined\n\ndef advanced_rna_structure_prediction(sequence, X_valid, y_valid, optimal_params):\n \"\"\"\n Advanced 3D RNA structure prediction pipeline integrated with existing code.\n \n Parameters:\n -----------\n sequence: str\n RNA sequence to predict the structure for\n X_valid, y_valid: \n Validation data from existing code (required for compatibility)\n optimal_params: dict\n Optimal parameters from existing code\n \n Returns:\n --------\n ensemble: list\n List of predicted 3D structures (ensemble of models)\n \"\"\"\n print(f\"Starting advanced 3D structure prediction for RNA sequence of size {len(sequence)}\")\n \n # 1. Predict RNA-RNA contact map\n print(\"Predicting contact map...\")\n contact_map = predict_rna_contacts(sequence)\n \n # 2. Build initial graph\n print(\"Building RNA graph...\")\n rna_graph = create_rna_graph(sequence, contact_map)\n \n # 3. Create simplified GNN model (or load if trained)\n print(\"Initializing GNN model...\")\n # In a complete implementation, you would load a pre-trained model\n # Here, we create a simple model for demonstration\n node_features = 11 # 5 for base one-hot encoding + 6 for additional features\n edge_features = 5 # 4 for edge type (one-hot) + 1 for weight\n \n # Initializing the model (note: untrained, just for example)\n gnn_model = RNAGraphTransformer(node_features, edge_features)\n \n # 4. Generate initial coordinates using existing reference model\n print(\"Generating initial coordinates with reference model...\")\n np.random.seed(optimal_params.get('seed', 42))\n reference_model = reference_based_approach(\n X_valid, \n y_valid,\n geometric_sampling=True,\n noise_level=optimal_params.get('noise', 0.21),\n correlation=optimal_params.get('corr', 0.83)\n )\n \n # Convert sequence to input format of reference model\n print(\"Preparing sequence for prediction...\")\n seq_features = np.zeros((1, 720, 5)) # One-hot encoding with padding to 720\n for i, base in enumerate(sequence):\n if i >= 720: # Limit to avoid out of bounds index\n break\n if base == 'A':\n seq_features[0, i, 0] = 1\n elif base == 'C':\n seq_features[0, i, 1] = 1\n elif base == 'G':\n seq_features[0, i, 2] = 1\n elif base == 'U' or base == 'T':\n seq_features[0, i, 3] = 1\n else:\n seq_features[0, i, 4] = 1 # Unknown base\n \n # FIX: Ensure correct dimensions for initial coordinates\n try:\n # Generate initial prediction with reference model\n predictions = reference_model.predict(seq_features)\n \n # Extract coordinates for current sequence and ensure correct dimensions\n raw_coords = predictions[0]\n \n # Ensure coordinates have exactly the length of the sequence\n initial_coords = np.zeros((len(sequence), 3))\n \n # Copy only the necessary length\n min_length = min(len(raw_coords), len(sequence))\n initial_coords[:min_length] = raw_coords[:min_length]\n \n print(f\"Generated initial coordinates with shape: {initial_coords.shape}\")\n except Exception as e:\n print(f\"Error generating initial coordinates: {str(e)}\")\n # Emergency coordinates in case of failure\n initial_coords = np.zeros((len(sequence), 3))\n # Create a simple linear structure\n for i in range(1, len(sequence)):\n initial_coords[i] = initial_coords[i-1] + np.array([3.8, 0, 0])\n print(\"Using emergency coordinates due to error.\")\n \n # 5. Generate fragment-based structure\n print(\"Generating fragment-based structure...\")\n fragment_library = load_fragment_library()\n try:\n fragment_based_coords = geometric_assembly(sequence, contact_map, fragment_library)\n \n # Ensure fragment-based structure has correct size\n if len(fragment_based_coords) != len(sequence):\n temp_coords = np.zeros((len(sequence), 3))\n min_length = min(len(fragment_based_coords), len(sequence))\n temp_coords[:min_length] = fragment_based_coords[:min_length]\n fragment_based_coords = temp_coords\n \n # Normalize fragment-based structure\n fragment_based_coords = normalize_structure(fragment_based_coords)\n except Exception as e:\n print(f\"Error in fragment-based assembly: {str(e)}\")\n # Fallback: use only reference model\n fragment_based_coords = initial_coords.copy()\n \n # 6. Combine evidence (weighted average of coordinates)\n print(\"Combining structures...\")\n combined_coords = weighted_coordinate_average([\n (initial_coords, 0.6), # Higher weight for reference model prediction\n (fragment_based_coords, 0.4)\n ])\n \n # 7. Refinement with geometric distance and structure\n print(\"Refining geometry...\")\n # Convert predicted contacts to target distance matrix\n distances = {}\n for i, j, prob in contact_map:\n if i < len(sequence) and j < len(sequence) and prob > 0.5: # Validate indices and use only high probability contacts\n # Typical base pair distance (~5-6Å)\n distances[(i, j)] = 5.5\n distances[(j, i)] = 5.5\n \n # Add backbone distances (consecutive bonds)\n for i in range(len(sequence) - 1):\n distances[(i, i+1)] = 3.8 # Typical RNA backbone distance\n distances[(i+1, i)] = 3.8\n \n # Refine with distance geometry\n try:\n distance_matrix = np.zeros((len(sequence), len(sequence)))\n for i in range(len(sequence)):\n for j in range(i+1, len(sequence)):\n dist = np.linalg.norm(combined_coords[i] - combined_coords[j])\n distance_matrix[i, j] = dist\n distance_matrix[j, i] = dist\n \n weights = np.zeros((len(sequence), len(sequence)))\n \n for (i, j), _ in distances.items():\n if i < len(weights) and j < len(weights): # Validate indices\n weights[i, j] = 1.0\n \n target_distances = distance_matrix.copy()\n for (i, j), dist in distances.items():\n if i < len(target_distances) and j < len(target_distances): # Validate indices\n target_distances[i, j] = dist\n target_distances[j, i] = dist\n \n refined_coords = refine_with_distance_geometry(\n combined_coords, target_distances, weights\n )\n except Exception as e:\n print(f\"Error in geometric refinement: {str(e)}\")\n refined_coords = combined_coords\n \n # 8. Final refinement with energy minimization\n print(\"Applying final refinement...\")\n try:\n # Calculate GC content for post_process_rna_structure\n gc_count = sequence.count('G') + sequence.count('C')\n gc_content = gc_count / len(sequence) if len(sequence) > 0 else 0.5\n \n # Use post_process_rna_structure instead of energy_based_refinement\n final_structure = post_process_rna_structure(\n refined_coords, \n sequence, \n gc_content, \n use_global_movement=True\n )\n \n # Also apply backbone structure refinement\n if 'refine_rna_backbone_with_dihedrals' in globals():\n final_structure = refine_rna_backbone_with_dihedrals(final_structure)\n elif 'refine_rna_backbone' in globals():\n final_structure = refine_rna_backbone(final_structure)\n except Exception as e:\n print(f\"Error in energy refinement: {str(e)}\")\n final_structure = refined_coords\n \n # 9. Quality assessment and ensemble generation\n print(\"Generating final model ensemble...\")\n ensemble = []\n \n # Add the main model\n ensemble.append(normalize_structure(final_structure))\n \n # Generate variations to complete the ensemble\n for i in range(1, 5): # Generate 4 additional variations\n try:\n # Use progressively larger perturbation scale\n perturbation = 0.05 * i\n variation = sample_structural_variation(\n final_structure,\n noise_level=perturbation,\n preserve_distance=True,\n use_global_movement=(i % 2 == 0),\n correlation=0.9 - (i * 0.05)\n )\n ensemble.append(normalize_structure(variation))\n except Exception as e:\n print(f\"Error generating variation {i}: {str(e)}\")\n # In case of error, duplicate the main model\n ensemble.append(normalize_structure(final_structure))\n \n print(\"Advanced 3D structure prediction completed.\")\n return ensemble\n\n# ==== 5. PREDICTION PIPELINE EXECUTION 🚀🧬 ====\n\ndef run_graph_based_pipeline(X_valid, y_valid, test_seq_df, sample_submission_df, output_dir, optimal_params):\n \"\"\"\n Pipeline that uses the graph-based approach for predicting 3D RNA structures.\n This function replaces or complements the original run_hybrid_pipeline function.\n \n Parameters:\n -----------\n X_valid, y_valid: Validation data\n test_seq_df: DataFrame with test sequences\n sample_submission_df: Example submission format\n output_dir: Output directory for files\n optimal_params: Optimized parameters for the model\n \n Returns:\n --------\n submission_df: DataFrame\n Submission file with predicted structures\n status_dict: dict\n Dictionary with pipeline execution status\n \"\"\"\n print(\"=\" * 80)\n print(\"ADVANCED PIPELINE: GRAPH-BASED APPROACH\".center(80))\n print(\"=\" * 80)\n \n status = {\n 'success': False,\n 'model_type': 'graph_based',\n 'structures_generated': 0,\n 'error': None\n }\n \n try:\n # Dictionary to store structures for each sequence\n seq_to_structures = {}\n \n # Process each test sequence\n for i, (_, row) in enumerate(test_seq_df.iterrows()):\n target_id = row['target_id']\n sequence = row['sequence']\n seq_length = len(sequence)\n \n print(f\"Processing sequence {i+1}/{len(test_seq_df)}, ID: {target_id}, length: {seq_length}\")\n \n # Generate 3D structures using advanced pipeline\n structures = advanced_rna_structure_prediction(\n sequence, \n X_valid, \n y_valid, \n optimal_params\n )\n \n # Store structures\n seq_to_structures[target_id] = structures\n status['structures_generated'] += 1\n \n # Create submission DataFrame\n print(\"\\nCreating submission file...\")\n submission_df = sample_submission_df.copy()\n \n for i, row in submission_df.iterrows():\n if i % 1000 == 0:\n print(f\"Processing row {i}/{len(submission_df)}\")\n \n # Parse ID to get sequence ID and residue index\n id_parts = row['ID'].split('_')\n seq_id = id_parts[0]\n residue_idx = int(id_parts[1]) - 1 # Convert to zero-based indexing\n \n # Check if we have structures for this sequence\n if seq_id in seq_to_structures:\n structures = seq_to_structures[seq_id]\n \n # Check if residue index is valid\n if residue_idx < len(structures[0]):\n # Fill coordinates for all 5 structures\n for struct_idx in range(5):\n if struct_idx < len(structures):\n submission_df.at[i, f'x_{struct_idx+1}'] = structures[struct_idx][residue_idx][0]\n submission_df.at[i, f'y_{struct_idx+1}'] = structures[struct_idx][residue_idx][1]\n submission_df.at[i, f'z_{struct_idx+1}'] = structures[struct_idx][residue_idx][2]\n else:\n # If we have fewer than 5 structures, duplicate the last one\n last_idx = len(structures) - 1\n submission_df.at[i, f'x_{struct_idx+1}'] = structures[last_idx][residue_idx][0]\n submission_df.at[i, f'y_{struct_idx+1}'] = structures[last_idx][residue_idx][1]\n submission_df.at[i, f'z_{struct_idx+1}'] = structures[last_idx][residue_idx][2]\n \n # Save submission\n graph_file = os.path.join(output_dir, 'submission_graph_based.csv')\n submission_df.to_csv(graph_file, index=False)\n print(f\"Graph-based submission saved at {graph_file}\")\n \n # Save as standard submission\n standard_file = os.path.join(output_dir, 'submission.csv')\n submission_df.to_csv(standard_file, index=False)\n \n # Mark as success\n status['success'] = True\n \n return submission_df, status\n \n except Exception as e:\n print(f\"ERROR in graph-based pipeline: {str(e)}\")\n import traceback\n traceback.print_exc()\n status['error'] = str(e)\n return None, status","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-04-09T15:13:38.076780Z","iopub.execute_input":"2025-04-09T15:13:38.077060Z","iopub.status.idle":"2025-04-09T15:13:38.139481Z","shell.execute_reply.started":"2025-04-09T15:13:38.077029Z","shell.execute_reply":"2025-04-09T15:13:38.138744Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"> ## Submission Creation","metadata":{"papermill":{"duration":0.017635,"end_time":"2025-03-26T03:43:08.747861","exception":false,"start_time":"2025-03-26T03:43:08.730226","status":"completed"},"tags":[]}},{"cell_type":"code","source":"def create_submission_dataframe(seq_to_coords, sample_submission_df):\n \"\"\"\n Create a submission DataFrame from the final structures.\n \n Parameters:\n -----------\n seq_to_coords: Dictionary mapping sequence IDs to lists of structures\n sample_submission_df: Sample submission format\n \n Returns:\n --------\n Submission DataFrame\n \"\"\"\n # Create a copy of the sample submission\n submission_df = sample_submission_df.copy()\n \n # Fill in the coordinates for each structure\n for i, row in submission_df.iterrows():\n if i % 1000 == 0:\n print(f\"Processing row {i}/{len(submission_df)}\")\n \n # Parse the ID to get sequence ID and residue index\n id_parts = row['ID'].split('_')\n seq_id = id_parts[0]\n residue_idx = int(id_parts[1]) - 1 # Convert to 0-based indexing\n \n # Check if we have structures for this sequence\n if seq_id in seq_to_coords:\n structures = seq_to_coords[seq_id]\n \n # Check if the residue index is valid\n if residue_idx < len(structures[0]):\n # Fill in coordinates for all 5 structures\n for struct_idx in range(5):\n if struct_idx < len(structures):\n submission_df.at[i, f'x_{struct_idx+1}'] = structures[struct_idx][residue_idx][0]\n submission_df.at[i, f'y_{struct_idx+1}'] = structures[struct_idx][residue_idx][1]\n submission_df.at[i, f'z_{struct_idx+1}'] = structures[struct_idx][residue_idx][2]\n else:\n # If we have fewer than 5 structures, duplicate the last one\n last_idx = len(structures) - 1\n submission_df.at[i, f'x_{struct_idx+1}'] = structures[last_idx][residue_idx][0]\n submission_df.at[i, f'y_{struct_idx+1}'] = structures[last_idx][residue_idx][1]\n submission_df.at[i, f'z_{struct_idx+1}'] = structures[last_idx][residue_idx][2]\n \n return submission_df\n\ndef generate_nn_pruned_submission(model, quality_model, test_seq_df, sample_submission_df):\n \"\"\"\n Enhanced submission generation that uses NN-based pruning for structure selection.\n \"\"\"\n print(\"Generating submission with Neural Network pruning...\")\n \n # Prepare test features\n X_test = prepare_test_features(test_seq_df)\n \n # Generate multiple predictions for ensemble diversity\n print(\"Generating base predictions...\")\n base_predictions = model.predict(X_test)\n \n seq_to_coords = {}\n for i, (_, row) in enumerate(test_seq_df.iterrows()):\n target_id = row['target_id']\n seq = row['sequence']\n seq_length = len(seq)\n \n print(f\"Processing sequence {i+1}/{len(test_seq_df)}, ID: {target_id}, length: {seq_length}\")\n \n # Get base coordinates\n base_coords = base_predictions[i][:seq_length]\n \n # Extract sequence features for this RNA\n seq_features = X_test[i][:seq_length]\n \n # Generate and prune structures using the NN model\n structures = generate_and_prune_structures(\n base_coords, \n seq_features, \n quality_model,\n num_candidates=30, # Generate more candidates\n top_k=5 # Keep top 5 for submission\n )\n \n # Store the structures\n seq_to_coords[target_id] = structures\n \n # Create submission DataFrame\n print(\"Creating submission file...\")\n submission_df = sample_submission_df.copy()\n \n for i, row in submission_df.iterrows():\n id_parts = row['ID'].split('_')\n seq_id = id_parts[0]\n residue_idx = int(id_parts[1]) - 1\n \n if seq_id in seq_to_coords:\n structures = seq_to_coords[seq_id]\n if residue_idx < len(structures[0]):\n for struct_idx in range(5):\n submission_df.at[i, f'x_{struct_idx+1}'] = structures[struct_idx][residue_idx][0]\n submission_df.at[i, f'y_{struct_idx+1}'] = structures[struct_idx][residue_idx][1]\n submission_df.at[i, f'z_{struct_idx+1}'] = structures[struct_idx][residue_idx][2]\n \n submission_file = os.path.join(OUTPUT_DIR, 'submission_nn_pruned.csv')\n submission_df.to_csv(submission_file, index=False)\n print(f\"NN-pruned submission file saved to {submission_file}\")\n \n # Also save as standard submission\n standard_file = os.path.join(OUTPUT_DIR, 'submission.csv')\n submission_df.to_csv(standard_file, index=False)\n \n return submission_df","metadata":{"execution":{"iopub.status.busy":"2025-04-09T15:13:38.140328Z","iopub.execute_input":"2025-04-09T15:13:38.140547Z","iopub.status.idle":"2025-04-09T15:13:38.151875Z","shell.execute_reply.started":"2025-04-09T15:13:38.140528Z","shell.execute_reply":"2025-04-09T15:13:38.151255Z"},"papermill":{"duration":0.029933,"end_time":"2025-03-26T03:43:08.795620","exception":false,"start_time":"2025-03-26T03:43:08.765687","status":"completed"},"tags":[],"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"> ## Main Execution Block","metadata":{"papermill":{"duration":0.053971,"end_time":"2025-03-26T03:43:09.046403","exception":false,"start_time":"2025-03-26T03:43:08.992432","status":"completed"},"tags":[]}},{"cell_type":"code","source":"if __name__ == \"__main__\":\n # Print startup banner\n print(\"=\" * 80)\n print(\"RNA 3D STRUCTURE PREDICTION PIPELINE\".center(80))\n print(\"=\" * 80)\n \n # Print selected mode\n mode_description = \"Advanced Graph-Based Pipeline: Modeling Long-Range Interactions\"\n print(f\"Selected mode: {mode_description}\")\n print(\"-\" * 80)\n \n try:\n # Start the graph-based pipeline\n start_time = time.time()\n \n print(\"Loading processed data...\")\n X_train, y_train, X_valid, y_valid = load_processed_data()\n \n print(\"\\nLoading test data...\")\n test_seq_df = pd.read_csv(os.path.join(DATA_DIR, \"test_sequences.csv\"))\n sample_submission_df = pd.read_csv(os.path.join(DATA_DIR, \"sample_submission.csv\"))\n \n # Optimal parameters (can be adjusted)\n optimal_params = {\n 'noise': 0.21,\n 'corr': 0.83,\n 'seed': 42\n }\n \n # Run the graph-based pipeline\n submission_df, status = run_graph_based_pipeline(\n X_valid, y_valid,\n test_seq_df, sample_submission_df,\n OUTPUT_DIR,\n optimal_params=optimal_params\n )\n \n # Calculate total runtime\n runtime = time.time() - start_time\n hours, remainder = divmod(runtime, 3600)\n minutes, seconds = divmod(remainder, 60)\n \n # Display results summary\n print(\"\\n\" + \"=\" * 80)\n print(\"GRAPH-BASED PIPELINE RESULTS SUMMARY\".center(80))\n print(\"=\" * 80)\n print(f\"Total runtime: {int(hours)}h {int(minutes)}m {int(seconds)}s\")\n \n if status['success']:\n print(\"\\nGRAPH-BASED PIPELINE STATISTICS:\")\n print(f\" - Model type: {status['model_type']}\")\n print(f\" - Structures generated: {status['structures_generated']}\")\n else:\n print(f\"\\nPipeline failed with error: {status['error']}\")\n \n # Display output file information\n print(\"\\nOUTPUT FILES:\")\n submission_file = os.path.join(OUTPUT_DIR, 'submission_graph_based.csv')\n if os.path.exists(submission_file):\n try:\n file_size = os.path.getsize(submission_file)\n print(f\" - Graph-based submission: {submission_file} ({file_size/1024/1024:.2f} MB)\")\n except:\n print(f\" - Graph-based submission: {submission_file}\")\n \n standard_file = os.path.join(OUTPUT_DIR, 'submission.csv')\n if os.path.exists(standard_file):\n try:\n file_size = os.path.getsize(standard_file)\n print(f\" - Standard submission: {standard_file} ({file_size/1024/1024:.2f} MB)\")\n except:\n print(f\" - Standard submission: {standard_file}\")\n \n print(\"=\" * 80)\n \n print(\"\\nProcess completed.\")\n \n except Exception as e:\n print(\"\\n\" + \"=\" * 80)\n print(\"ERROR IN MAIN EXECUTION\".center(80))\n print(\"=\" * 80)\n print(f\"Critical error: {str(e)}\")\n traceback.print_exc()\n print(\"=\" * 80)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-04-09T15:13:38.152665Z","iopub.execute_input":"2025-04-09T15:13:38.152897Z","iopub.status.idle":"2025-04-09T15:20:21.812752Z","shell.execute_reply.started":"2025-04-09T15:13:38.152866Z","shell.execute_reply":"2025-04-09T15:20:21.812029Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"submission_df = pd.read_csv('/kaggle/working/submission.csv')\nprint(\"Overview of the DataFrame:\")\nprint(submission_df.shape) # Print the shape (rows, columns)\nprint(submission_df.head()) # Display the first 5 rows","metadata":{"execution":{"iopub.status.busy":"2025-04-09T15:20:21.813567Z","iopub.execute_input":"2025-04-09T15:20:21.813823Z","iopub.status.idle":"2025-04-09T15:20:21.833678Z","shell.execute_reply.started":"2025-04-09T15:20:21.813803Z","shell.execute_reply":"2025-04-09T15:20:21.833025Z"},"papermill":{"duration":0.068636,"end_time":"2025-03-26T03:44:52.857355","exception":false,"start_time":"2025-03-26T03:44:52.788719","status":"completed"},"tags":[],"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## 🧬 RNA Structure: Real vs Predicted 🔍","metadata":{}},{"cell_type":"code","source":"def normalize_for_visualization(coords):\n \"\"\"\n Normaliza e centraliza coordenadas para visualização consistente.\n Lida com valores inválidos e padroniza a escala.\n \n Parameters:\n -----------\n coords : numpy.ndarray\n Coordenadas 3D para normalizar\n \n Returns:\n --------\n numpy.ndarray\n Coordenadas normalizadas\n \"\"\"\n import numpy as np\n \n # Cria uma cópia para evitar modificar o original\n normalized = coords.copy()\n \n # Identifica coordenadas válidas (não-zero e não-NaN)\n valid_mask = ~np.all(normalized == 0, axis=1) & ~np.any(np.isnan(normalized), axis=1)\n \n # Se não houver coordenadas válidas, retornar as originais\n if not np.any(valid_mask):\n print(\"AVISO: Nenhuma coordenada válida encontrada para normalização\")\n return normalized\n \n # Extrair apenas coordenadas válidas\n valid_coords = normalized[valid_mask]\n \n # 1. Centralizar na origem\n center = np.mean(valid_coords, axis=0)\n valid_coords = valid_coords - center\n \n # 2. Normalizar para escala padrão (valores máximos entre -50 e 50)\n max_dist = np.max(np.abs(valid_coords))\n if max_dist > 0:\n scale_factor = 50.0 / max_dist\n valid_coords = valid_coords * scale_factor\n \n # Aplicar transformações apenas às coordenadas válidas\n normalized[valid_mask] = valid_coords\n \n return normalized\n\ndef visualize_rna_structure_comparison(sequence_str, real_structure, X_valid, y_valid, optimal_params, title=None):\n \"\"\"\n Visualize comparison between real and predicted RNA 3D structures using the graph-based approach.\n Improved with better normalization and error handling.\n \n Parameters:\n -----------\n sequence_str : str\n RNA sequence as string\n real_structure : array\n True structure coordinates \n X_valid, y_valid : array\n Validation data for the model\n optimal_params : dict\n Parameters for the model\n title : str, optional\n Title for the visualization\n \"\"\"\n import matplotlib.pyplot as plt\n import numpy as np\n \n # Verificar se real_structure contém dados válidos\n if real_structure is None or np.all(np.isnan(real_structure)) or len(real_structure) == 0:\n print(\"ERRO: Estrutura real inválida ou vazia\")\n return None\n \n # Generate prediction using the graph-based approach\n try:\n structures = advanced_rna_structure_prediction(\n sequence_str, \n X_valid, \n y_valid, \n optimal_params\n )\n \n # Take the best structure (first one in the ensemble)\n predicted_structure = structures[0]\n except Exception as e:\n print(f\"ERRO ao gerar previsão: {str(e)}\")\n return None\n \n # Verificar se predicted_structure contém dados válidos\n if predicted_structure is None or np.all(np.isnan(predicted_structure)) or len(predicted_structure) == 0:\n print(\"ERRO: Estrutura prevista inválida ou vazia\")\n return None\n \n # Ensure structures have the same length\n min_length = min(len(real_structure), len(predicted_structure))\n real_structure = real_structure[:min_length].copy()\n predicted_structure = predicted_structure[:min_length].copy()\n \n # Normalizar ambas estruturas para visualização consistente\n real_structure_viz = normalize_for_visualization(real_structure)\n predicted_structure_viz = normalize_for_visualization(predicted_structure)\n \n # Criar uma figura com dois subplots lado a lado\n fig = plt.figure(figsize=(16, 6))\n \n # Plot real structure\n ax1 = fig.add_subplot(121, projection='3d')\n ax1.set_title('Real RNA Structure', fontsize=12)\n \n # Verificar se há dados válidos para plotar\n valid_mask_real = ~np.all(real_structure_viz == 0, axis=1) & ~np.any(np.isnan(real_structure_viz), axis=1)\n if np.any(valid_mask_real):\n # Plot points\n ax1.scatter(real_structure_viz[valid_mask_real, 0], \n real_structure_viz[valid_mask_real, 1],\n real_structure_viz[valid_mask_real, 2],\n c=np.arange(np.sum(valid_mask_real)),\n cmap='viridis',\n s=50)\n \n # Connect consecutive points to show backbone\n ax1.plot(real_structure_viz[valid_mask_real, 0],\n real_structure_viz[valid_mask_real, 1], \n real_structure_viz[valid_mask_real, 2],\n color='gray',\n alpha=0.5,\n linewidth=2)\n else:\n ax1.text(0, 0, 0, \"No valid data\", ha='center', va='center', fontsize=14)\n \n ax1.set_xlabel('X')\n ax1.set_ylabel('Y')\n ax1.set_zlabel('Z')\n \n # Definir limites consistentes\n ax1.set_xlim([-60, 60])\n ax1.set_ylim([-60, 60])\n ax1.set_zlim([-60, 60])\n \n # Plot predicted structure\n ax2 = fig.add_subplot(122, projection='3d')\n ax2.set_title('Predicted RNA Structure (Graph-Based)', fontsize=12)\n \n # Verificar se há dados válidos para plotar\n valid_mask_pred = ~np.all(predicted_structure_viz == 0, axis=1) & ~np.any(np.isnan(predicted_structure_viz), axis=1)\n if np.any(valid_mask_pred):\n # Plot points\n ax2.scatter(predicted_structure_viz[valid_mask_pred, 0],\n predicted_structure_viz[valid_mask_pred, 1], \n predicted_structure_viz[valid_mask_pred, 2],\n c=np.arange(np.sum(valid_mask_pred)),\n cmap='plasma',\n s=50)\n \n # Connect consecutive points to show backbone \n ax2.plot(predicted_structure_viz[valid_mask_pred, 0],\n predicted_structure_viz[valid_mask_pred, 1],\n predicted_structure_viz[valid_mask_pred, 2], \n color='red',\n alpha=0.5,\n linewidth=2)\n else:\n ax2.text(0, 0, 0, \"No valid data\", ha='center', va='center', fontsize=14)\n \n ax2.set_xlabel('X') \n ax2.set_ylabel('Y')\n ax2.set_zlabel('Z')\n \n # Definir limites consistentes\n ax2.set_xlim([-60, 60])\n ax2.set_ylim([-60, 60])\n ax2.set_zlim([-60, 60])\n \n # Overall title if provided\n if title:\n fig.suptitle(title, fontsize=16)\n \n fig.tight_layout()\n plt.show()\n \n return predicted_structure # Return for metrics calculation\n\ndef visualize_ensemble_structures(sequence_str, X_valid, y_valid, optimal_params, title=None):\n \"\"\"\n Visualize all 5 structures in the ensemble generated by the graph-based approach.\n Enhanced with better normalization and error handling.\n \n Parameters:\n -----------\n sequence_str : str\n RNA sequence as string\n X_valid, y_valid : array\n Validation data for the model\n optimal_params : dict\n Parameters for the model\n title : str, optional\n Title for the visualization\n \"\"\"\n import matplotlib.pyplot as plt\n import numpy as np\n \n # Generate prediction using the graph-based approach\n try:\n structures = advanced_rna_structure_prediction(\n sequence_str, \n X_valid, \n y_valid, \n optimal_params\n )\n except Exception as e:\n print(f\"ERRO ao gerar previsões: {str(e)}\")\n return\n \n # Normalize all structures for visualization\n normalized_structures = []\n for structure in structures:\n if structure is not None and len(structure) > 0:\n norm_struct = normalize_for_visualization(structure)\n normalized_structures.append(norm_struct)\n \n if not normalized_structures:\n print(\"ERRO: Nenhuma estrutura válida para visualizar\")\n return\n \n # Create a figure with subplots for all structures\n fig = plt.figure(figsize=(20, 10))\n \n # Plot each structure in the ensemble\n for i, structure in enumerate(normalized_structures):\n ax = fig.add_subplot(1, len(normalized_structures), i+1, projection='3d')\n ax.set_title(f'Ensemble Structure {i+1}', fontsize=12)\n \n # Verificar dados válidos\n valid_mask = ~np.all(structure == 0, axis=1) & ~np.any(np.isnan(structure), axis=1)\n if np.any(valid_mask):\n # Plot points\n ax.scatter(structure[valid_mask, 0], \n structure[valid_mask, 1],\n structure[valid_mask, 2],\n c=np.arange(np.sum(valid_mask)),\n cmap='plasma',\n s=40)\n \n # Connect consecutive points to show backbone\n ax.plot(structure[valid_mask, 0],\n structure[valid_mask, 1], \n structure[valid_mask, 2],\n color='red',\n alpha=0.5,\n linewidth=2)\n else:\n ax.text(0, 0, 0, \"No valid data\", ha='center', va='center', fontsize=14)\n \n ax.set_xlabel('X')\n ax.set_ylabel('Y')\n ax.set_zlabel('Z')\n \n # Definir limites consistentes\n ax.set_xlim([-60, 60])\n ax.set_ylim([-60, 60])\n ax.set_zlim([-60, 60])\n \n # Overall title if provided\n if title:\n fig.suptitle(title, fontsize=16)\n \n fig.tight_layout()\n plt.show()\n\ndef calculate_structure_metrics(real_structure, predicted_structure):\n \"\"\"\n Calculate key metrics to compare real and predicted structures.\n Enhanced with better handling of invalid data.\n \n Parameters:\n -----------\n real_structure : numpy.ndarray\n Original 3D structure coordinates\n predicted_structure : numpy.ndarray\n Predicted 3D structure coordinates\n \n Returns:\n --------\n metrics : dict\n Dictionary of comparison metrics\n \"\"\"\n import numpy as np\n \n # Verificar se as estruturas são válidas\n if (real_structure is None or predicted_structure is None or\n len(real_structure) == 0 or len(predicted_structure) == 0):\n print(\"ERRO: Estruturas vazias ou inválidas\")\n return {\n 'Distance MAE': 0.0,\n 'Coordinate RMSE': 0.0,\n 'Structural Similarity': 0.0,\n 'TM-Score': 0.0\n }\n \n # Ensure structures are the same length\n min_length = min(len(real_structure), len(predicted_structure))\n real_structure = real_structure[:min_length].copy()\n predicted_structure = predicted_structure[:min_length].copy()\n \n # Identificar coordenadas válidas em ambas estruturas\n valid_mask = (~np.all(real_structure == 0, axis=1) & \n ~np.any(np.isnan(real_structure), axis=1) &\n ~np.all(predicted_structure == 0, axis=1) &\n ~np.any(np.isnan(predicted_structure), axis=1))\n \n # Verificar se há coordenadas válidas suficientes\n if np.sum(valid_mask) < 3:\n print(\"AVISO: Menos de 3 coordenadas válidas para comparação\")\n return {\n 'Distance MAE': 0.0,\n 'Coordinate RMSE': 0.0,\n 'Structural Similarity': 0.0,\n 'TM-Score': 0.0\n }\n \n # Usar apenas coordenadas válidas\n real_valid = real_structure[valid_mask]\n pred_valid = predicted_structure[valid_mask]\n \n try:\n # Calculate pairwise distances\n real_dist_matrix = np.linalg.norm(\n real_valid[:, np.newaxis] - real_valid, \n axis=2\n )\n pred_dist_matrix = np.linalg.norm(\n pred_valid[:, np.newaxis] - pred_valid, \n axis=2\n )\n \n # Mean absolute error of distances\n distance_mae = np.mean(np.abs(real_dist_matrix - pred_dist_matrix))\n \n # Root Mean Squared Error (RMSE) of coordinates\n rmse = np.sqrt(np.mean((real_valid - pred_valid)**2))\n \n # Structural similarity (cosine similarity of distance matrices)\n try:\n similarity = np.corrcoef(\n real_dist_matrix.ravel(), \n pred_dist_matrix.ravel()\n )[0, 1]\n except:\n print(\"AVISO: Erro ao calcular similaridade estrutural\")\n similarity = 0.0\n \n # Calculate TM-score\n try:\n tm_score = calculate_tm_score(pred_valid, real_valid)\n except:\n print(\"AVISO: Erro ao calcular TM-score\")\n tm_score = 0.0\n \n except Exception as e:\n print(f\"ERRO ao calcular métricas: {str(e)}\")\n return {\n 'Distance MAE': 0.0,\n 'Coordinate RMSE': 0.0,\n 'Structural Similarity': 0.0,\n 'TM-Score': 0.0\n }\n \n return {\n 'Distance MAE': distance_mae,\n 'Coordinate RMSE': rmse, \n 'Structural Similarity': similarity,\n 'TM-Score': tm_score\n }\n\ndef plot_structure_metrics(metrics):\n \"\"\"\n Visualize structure comparison metrics.\n Enhanced with better error handling.\n \n Parameters:\n -----------\n metrics : dict \n Dictionary of comparison metrics\n \"\"\"\n import matplotlib.pyplot as plt\n import numpy as np\n \n # Verificar se há métricas válidas\n if not metrics:\n print(\"ERRO: Métricas vazias ou inválidas\")\n return\n \n # Garantir que os valores são números válidos\n for k, v in list(metrics.items()):\n if v is None or np.isnan(v) or np.isinf(v):\n print(f\"AVISO: Valor inválido para {k}, substituindo por 0.0\")\n metrics[k] = 0.0\n \n fig, ax = plt.subplots(figsize=(12, 6))\n metrics_names = list(metrics.keys())\n metrics_values = list(metrics.values())\n \n bars = ax.bar(metrics_names, metrics_values, color=['#1f77b4', '#2ca02c', '#d62728', '#9467bd'])\n ax.set_title('🧬 RNA Structure Prediction Metrics (Graph-Based Approach) 🔬', fontsize=14)\n ax.set_ylabel('Metric Value', fontsize=12)\n ax.tick_params(axis='x', labelrotation=45)\n \n # Add value labels on top of each bar\n for bar, v in zip(bars, metrics_values):\n ax.text(bar.get_x() + bar.get_width()/2., v, \n f'{v:.4f}', ha='center', va='bottom')\n \n # Ajustar limites do eixo y para garantir que todos os valores sejam visíveis\n y_min, y_max = ax.get_ylim()\n ax.set_ylim(min(y_min, -0.05), max(y_max, max(metrics_values) * 1.1))\n \n fig.tight_layout()\n plt.show()\n\ndef main_graph_visualization():\n \"\"\"\n Main function to visualize RNA structures predicted with the graph-based approach.\n Enhanced with better error handling.\n \"\"\"\n import numpy as np\n \n # Load processed data \n try:\n X_train, y_train, X_valid, y_valid = load_processed_data()\n except Exception as e:\n print(f\"ERRO ao carregar dados: {str(e)}\")\n return\n \n # Verificar se os dados são válidos\n if X_valid is None or y_valid is None or len(X_valid) == 0 or len(y_valid) == 0:\n print(\"ERRO: Dados de validação vazios ou inválidos\")\n return\n \n # Set optimal parameters\n optimal_params = {\n 'noise': 0.21,\n 'corr': 0.83,\n 'seed': 42\n }\n \n # Visualize multiple sequences\n num_sequences = min(3, len(X_valid))\n \n for i in range(num_sequences):\n try:\n print(f\"\\nVisualizando sequência {i+1}\")\n \n # Converter one-hot para sequência de bases\n try:\n sequence = np.argmax(X_valid[i], axis=-1)\n base_map = {0: 'A', 1: 'C', 2: 'G', 3: 'U', 4: 'N'}\n valid_indices = np.where(X_valid[i].sum(axis=1) > 0)[0]\n \n if len(valid_indices) == 0:\n print(\"AVISO: Nenhum índice válido na sequência\")\n sequence_str = \"N\" * 10 # Sequência padrão\n else:\n sequence_str = ''.join(base_map[sequence[j]] for j in valid_indices)\n \n print(f\"RNA Sequence: {sequence_str[:10]}... (length: {len(sequence_str)})\")\n except Exception as e:\n print(f\"ERRO ao extrair sequência: {str(e)}\")\n sequence_str = \"N\" * 10 # Sequência padrão em caso de erro\n \n # Verificar se a estrutura real é válida\n if np.all(np.isnan(y_valid[i])):\n print(\"AVISO: Estrutura real contém apenas valores NaN\")\n continue\n \n # Visualize comparison\n predicted_structure = visualize_rna_structure_comparison(\n sequence_str, \n y_valid[i], \n X_valid, \n y_valid, \n optimal_params,\n title=f'🧬 RNA Structure Comparison - Sequence {i+1} (Graph-Based Method) 🔍'\n )\n \n if predicted_structure is not None:\n # Calculate and plot metrics\n metrics = calculate_structure_metrics(y_valid[i], predicted_structure)\n plot_structure_metrics(metrics)\n \n # Visualize all structures in the ensemble\n visualize_ensemble_structures(\n sequence_str,\n X_valid, \n y_valid, \n optimal_params,\n title=f'🧬 RNA Structure Ensemble - Sequence {i+1} 🔍'\n )\n \n except Exception as e:\n print(f\"ERRO ao processar sequência {i+1}: {str(e)}\")\n import traceback\n traceback.print_exc()\n\n# Para executar a visualização\nmain_graph_visualization()","metadata":{"papermill":{"duration":0.037766,"end_time":"2025-03-26T03:44:52.935922","exception":false,"start_time":"2025-03-26T03:44:52.898156","status":"completed"},"tags":[],"trusted":true,"execution":{"iopub.status.busy":"2025-04-09T15:28:01.190072Z","iopub.execute_input":"2025-04-09T15:28:01.190395Z","iopub.status.idle":"2025-04-09T15:30:34.041495Z","shell.execute_reply.started":"2025-04-09T15:28:01.190368Z","shell.execute_reply":"2025-04-09T15:30:34.040639Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"","metadata":{"trusted":true},"outputs":[],"execution_count":null}]}