{"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"name":"python","version":"3.12.12","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"none","dataSources":[{"sourceId":118765,"databundleVersionId":15231210,"sourceType":"competition"}],"dockerImageVersionId":31259,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# Stanford RNA 3D Folding Part 2 - Complete EDA\n\n\n\n## 📊 Analysis Sections\n1. Setup and Data Loading\n2. Train/Validation/Test Split Analysis\n3. Sequence Analysis (RNA sequences)\n4. Structure Analysis (PDB files)\n5. MSA Analysis (Multiple Sequence Alignments)\n6. Metadata Analysis\n7. Label Distribution\n8. Data Quality Assessment\n9. Feature Engineering\n10. Correlations and Patterns\n11. Summary and Insights\n\n---","metadata":{}},{"cell_type":"markdown","source":"## 1. Setup and Configuration","metadata":{}},{"cell_type":"code","source":"# Core libraries\nimport pandas as pd\nimport numpy as np\nimport os\nimport gc\nimport warnings\nfrom pathlib import Path\nfrom collections import Counter, defaultdict\nimport json\nimport glob\n\n# Visualization\nimport matplotlib.pyplot as plt\nimport seaborn as sns\nimport plotly.express as px\nimport plotly.graph_objects as go\nfrom plotly.subplots import make_subplots\n\n# Bioinformatics (optional but recommended)\ntry:\n    from Bio import SeqIO\n    from Bio.Seq import Seq\n    from Bio.PDB import PDBParser, PDBIO, Select\n    BIOPYTHON_AVAILABLE = True\n    print(\"✅ BioPython available\")\nexcept ImportError:\n    BIOPYTHON_AVAILABLE = False\n    print(\"⚠️ BioPython not available. Install with: pip install biopython\")\n\n# Statistical analysis\nfrom scipy import stats\nfrom scipy.spatial.distance import pdist, squareform\nfrom sklearn.preprocessing import StandardScaler\nfrom sklearn.decomposition import PCA\n\n# Progress bars\nfrom tqdm.notebook import tqdm\n\n# Configure\nwarnings.filterwarnings('ignore')\nplt.style.use('seaborn-v0_8-darkgrid')\nsns.set_palette(\"husl\")\n%matplotlib inline\n\npd.set_option('display.max_columns', None)\npd.set_option('display.max_rows', 100)\npd.set_option('display.float_format', lambda x: f'{x:.4f}')\n\nprint(\"✅ All libraries imported successfully!\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-02-06T18:32:03.172306Z","iopub.execute_input":"2026-02-06T18:32:03.172669Z","iopub.status.idle":"2026-02-06T18:32:09.638775Z","shell.execute_reply.started":"2026-02-06T18:32:03.172627Z","shell.execute_reply":"2026-02-06T18:32:09.637578Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Configuration - UPDATE THIS PATH!\nCONFIG = {\n    'DATA_PATH': '/kaggle/input/stanford-rna-3d-folding-2/',\n    'CACHE_DIR': './cache/',\n    'OUTPUT_DIR': './outputs/',\n    'RANDOM_SEED': 42,\n    'FIGURE_SIZE': (14, 8),\n    'DPI': 100,\n    'NUCLEOTIDES': ['A', 'C', 'G', 'U'],\n}\n\n# Create directories\nfor dir_path in [CONFIG['CACHE_DIR'], CONFIG['OUTPUT_DIR']]:\n    os.makedirs(dir_path, exist_ok=True)\n\nnp.random.seed(CONFIG['RANDOM_SEED'])\n\nprint(f\"📂 Data path: {CONFIG['DATA_PATH']}\")\nprint(f\"💾 Cache: {CONFIG['CACHE_DIR']}\")\nprint(f\"📊 Output: {CONFIG['OUTPUT_DIR']}\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-02-06T18:32:09.641643Z","iopub.execute_input":"2026-02-06T18:32:09.643043Z","iopub.status.idle":"2026-02-06T18:32:09.653799Z","shell.execute_reply.started":"2026-02-06T18:32:09.642977Z","shell.execute_reply":"2026-02-06T18:32:09.652416Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Utility functions\ndef print_section(title, char='=', length=80):\n    print(f\"\\n{char * length}\")\n    print(f\"{title.center(length)}\")\n    print(f\"{char * length}\\n\")\n\ndef get_file_size(filepath):\n    \"\"\"Get file size in MB\"\"\"\n    return os.path.getsize(filepath) / (1024**2)\n\ndef calculate_gc_content(sequence):\n    \"\"\"Calculate GC content\"\"\"\n    gc = sequence.count('G') + sequence.count('C')\n    return (gc / len(sequence)) * 100 if len(sequence) > 0 else 0\n\ndef get_nucleotide_composition(sequence):\n    \"\"\"Get nucleotide percentages\"\"\"\n    total = len(sequence)\n    return {nuc: (sequence.count(nuc) / total) * 100 for nuc in CONFIG['NUCLEOTIDES']}\n\ndef parse_pdb_c1_coords(pdb_file):\n    \"\"\"Parse C1' coordinates from PDB file\"\"\"\n    coords = []\n    try:\n        with open(pdb_file, 'r') as f:\n            for line in f:\n                if line.startswith('ATOM') and \"C1'\" in line:\n                    x = float(line[30:38].strip())\n                    y = float(line[38:46].strip())\n                    z = float(line[46:54].strip())\n                    coords.append([x, y, z])\n    except:\n        pass\n    return np.array(coords)\n\ndef calculate_radius_of_gyration(coords):\n    \"\"\"Calculate Rg\"\"\"\n    if len(coords) == 0:\n        return 0\n    centroid = np.mean(coords, axis=0)\n    return np.sqrt(np.mean(np.sum((coords - centroid)**2, axis=1)))\n\nprint(\"✅ Utility functions loaded\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-02-06T18:32:09.65553Z","iopub.execute_input":"2026-02-06T18:32:09.655875Z","iopub.status.idle":"2026-02-06T18:32:09.748222Z","shell.execute_reply.started":"2026-02-06T18:32:09.655847Z","shell.execute_reply":"2026-02-06T18:32:09.74656Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"---\n## 2. Data Loading and File Inventory","metadata":{}},{"cell_type":"code","source":"print_section(\"📂 FILE INVENTORY\")\n\ndata_path = Path(CONFIG['DATA_PATH'])\n\n# Check main directory\nif data_path.exists():\n    print(\"Main directory structure:\")\n    print(\"=\"*60)\n    \n    all_files = []\n    for item in sorted(data_path.iterdir()):\n        if item.is_file():\n            size = get_file_size(item)\n            all_files.append({'name': item.name, 'size_mb': size, 'type': 'file'})\n            print(f\"📄 {item.name:<40} {size:>10.2f} MB\")\n        elif item.is_dir():\n            num_files = len(list(item.glob('*')))\n            all_files.append({'name': item.name, 'files': num_files, 'type': 'dir'})\n            print(f\"📁 {item.name:<40} ({num_files} files)\")\n    \n    # Check subdirectories\n    print(\"\\n\" + \"=\"*60)\n    print(\"Subdirectory details:\")\n    print(\"=\"*60)\n    \n    # MSA directory\n    msa_dir = data_path / 'MSA'\n    if msa_dir.exists():\n        msa_files = list(msa_dir.glob('*'))\n        print(f\"\\n📁 MSA/ - {len(msa_files)} files\")\n        if msa_files:\n            extensions = Counter([f.suffix for f in msa_files])\n            for ext, count in extensions.most_common():\n                print(f\"   • {ext if ext else 'no extension'}: {count} files\")\n    \n    # PDB directory\n    pdb_dir = data_path / 'PDB_RNA'\n    if pdb_dir.exists():\n        pdb_files = list(pdb_dir.glob('*.pdb'))\n        print(f\"\\n📁 PDB_RNA/ - {len(pdb_files)} PDB files\")\n        if pdb_files:\n            total_size = sum([get_file_size(f) for f in pdb_files])\n            print(f\"   • Total size: {total_size:.2f} MB\")\n            print(f\"   • Average size: {total_size/len(pdb_files):.3f} MB\")\n    \n    # Extra directory\n    extra_dir = data_path / 'extra'\n    if extra_dir.exists():\n        extra_files = list(extra_dir.glob('*'))\n        print(f\"\\n📁 extra/ - {len(extra_files)} files\")\n        for f in extra_files:\n            if f.is_file():\n                size = get_file_size(f)\n                print(f\"   • {f.name:<35} {size:>8.2f} MB\")\n                \nelse:\n    print(f\"⚠️ Data path not found: {CONFIG['DATA_PATH']}\")\n    print(\"Please update CONFIG['DATA_PATH'] in the configuration cell.\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-02-06T18:32:09.751294Z","iopub.execute_input":"2026-02-06T18:32:09.751708Z","iopub.status.idle":"2026-02-06T18:32:10.114662Z","shell.execute_reply.started":"2026-02-06T18:32:09.751669Z","shell.execute_reply":"2026-02-06T18:32:10.113341Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"print_section(\"📊 LOADING CSV FILES\")\n\n# Load all CSV files\ncsv_files = {\n    'train_sequences': 'train_sequences.csv',\n    'train_labels': 'train_labels.csv',\n    'validation_sequences': 'validation_sequences.csv',\n    'validation_labels': 'validation_labels.csv',\n    'test_sequences': 'test_sequences.csv',\n    'sample_submission': 'sample_submission.csv',\n    'rna_metadata': 'extra/rna_metadata.csv',\n}\n\ndataframes = {}\n\nfor name, filepath in csv_files.items():\n    full_path = data_path / filepath\n    if full_path.exists():\n        df = pd.read_csv(full_path)\n        dataframes[name] = df\n        print(f\"✅ {name:<25} Shape: {str(df.shape):<15} Memory: {df.memory_usage(deep=True).sum()/(1024**2):.2f} MB\")\n    else:\n        print(f\"❌ {name:<25} NOT FOUND: {filepath}\")\n\nprint(f\"\\n📈 Total datasets loaded: {len(dataframes)}\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-02-06T18:32:10.115875Z","iopub.execute_input":"2026-02-06T18:32:10.116188Z","iopub.status.idle":"2026-02-06T18:32:30.139561Z","shell.execute_reply.started":"2026-02-06T18:32:10.11616Z","shell.execute_reply":"2026-02-06T18:32:30.138279Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Display quick preview of each dataset\nprint_section(\"👀 DATASET PREVIEWS\")\n\nfor name, df in dataframes.items():\n    print(f\"\\n{'='*70}\")\n    print(f\"{name.upper()}\")\n    print(f\"{'='*70}\")\n    print(f\"Shape: {df.shape}\")\n    print(f\"\\nColumns: {', '.join(df.columns.tolist())}\")\n    print(f\"\\nFirst 3 rows:\")\n    display(df.head(3))\n    print(f\"\\nData types:\")\n    print(df.dtypes)\n    print(f\"\\nMissing values:\")\n    missing = df.isnull().sum()\n    if missing.sum() > 0:\n        print(missing[missing > 0])\n    else:\n        print(\"None\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-02-06T18:32:30.140717Z","iopub.execute_input":"2026-02-06T18:32:30.140994Z","iopub.status.idle":"2026-02-06T18:32:31.722955Z","shell.execute_reply.started":"2026-02-06T18:32:30.140969Z","shell.execute_reply":"2026-02-06T18:32:31.721859Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"---\n## 3. Train/Validation/Test Split Analysis","metadata":{}},{"cell_type":"code","source":"print_section(\"🔍 DATA SPLIT ANALYSIS\")\n\n# Get split sizes\nsplit_info = {\n    'Train': len(dataframes.get('train_sequences', pd.DataFrame())),\n    'Validation': len(dataframes.get('validation_sequences', pd.DataFrame())),\n    'Test': len(dataframes.get('test_sequences', pd.DataFrame())),\n}\n\nprint(\"Dataset Split:\")\nprint(\"=\"*50)\ntotal = sum(split_info.values())\nfor split, count in split_info.items():\n    pct = (count / total * 100) if total > 0 else 0\n    print(f\"{split:<15} {count:>6} samples ({pct:>5.1f}%)\")\nprint(f\"{'Total':<15} {total:>6} samples\")\n\n# Visualize split\nfig, axes = plt.subplots(1, 2, figsize=(14, 5))\n\n# Pie chart\ncolors = ['#FF6B6B', '#4ECDC4', '#95E1D3']\naxes[0].pie(split_info.values(), labels=split_info.keys(), autopct='%1.1f%%',\n           colors=colors, startangle=90)\naxes[0].set_title('Dataset Split Distribution', fontsize=14, fontweight='bold')\n\n# Bar chart\naxes[1].bar(split_info.keys(), split_info.values(), color=colors, alpha=0.7, edgecolor='black')\naxes[1].set_ylabel('Number of Samples', fontsize=12)\naxes[1].set_title('Dataset Split Counts', fontsize=14, fontweight='bold')\nfor i, (k, v) in enumerate(split_info.items()):\n    axes[1].text(i, v + max(split_info.values())*0.02, str(v), \n                ha='center', va='bottom', fontsize=11, fontweight='bold')\n\nplt.tight_layout()\nplt.savefig(f\"{CONFIG['OUTPUT_DIR']}/data_split.png\", dpi=CONFIG['DPI'], bbox_inches='tight')\nplt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-02-06T18:32:31.72641Z","iopub.execute_input":"2026-02-06T18:32:31.726747Z","iopub.status.idle":"2026-02-06T18:32:32.344596Z","shell.execute_reply.started":"2026-02-06T18:32:31.726719Z","shell.execute_reply":"2026-02-06T18:32:32.343165Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"---\n## 4. Sequence Analysis","metadata":{}},{"cell_type":"code","source":"print_section(\"🧬 SEQUENCE ANALYSIS\")\n\n# Combine sequences for analysis\nseq_datasets = {\n    'train': dataframes.get('train_sequences'),\n    'validation': dataframes.get('validation_sequences'),\n    'test': dataframes.get('test_sequences'),\n}\n\n# Analyze each dataset\nfor name, df in seq_datasets.items():\n    if df is not None and 'sequence' in df.columns:\n        print(f\"\\n{name.upper()} Sequences:\")\n        print(\"=\"*50)\n        \n        # Add sequence length\n        df['seq_length'] = df['sequence'].str.len()\n        \n        # Calculate features\n        print(\"Calculating sequence features...\")\n        df['gc_content'] = df['sequence'].apply(calculate_gc_content)\n        \n        # Nucleotide composition\n        comp_list = []\n        for seq in tqdm(df['sequence'].head(1000), desc=f\"Processing {name}\"):\n            comp_list.append(get_nucleotide_composition(seq))\n        \n        comp_df = pd.DataFrame(comp_list)\n        for nuc in CONFIG['NUCLEOTIDES']:\n            df.loc[:len(comp_df)-1, nuc] = comp_df[nuc].values\n        \n        # Statistics\n        print(f\"\\nSequence Length Statistics:\")\n        print(df['seq_length'].describe())\n        \n        print(f\"\\nGC Content Statistics:\")\n        print(df['gc_content'].describe())\n        \n        # Update dataframe\n        seq_datasets[name] = df","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-02-06T18:32:32.345661Z","iopub.execute_input":"2026-02-06T18:32:32.345957Z","iopub.status.idle":"2026-02-06T18:32:32.538323Z","shell.execute_reply.started":"2026-02-06T18:32:32.345932Z","shell.execute_reply":"2026-02-06T18:32:32.537005Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Visualize sequence features across datasets\nfig, axes = plt.subplots(2, 3, figsize=(18, 10))\n\ncolors_map = {'train': '#FF6B6B', 'validation': '#4ECDC4', 'test': '#95E1D3'}\n\n# Row 1: Sequence length\nfor idx, (name, df) in enumerate(seq_datasets.items()):\n    if df is not None and 'seq_length' in df.columns:\n        axes[0, idx].hist(df['seq_length'], bins=50, alpha=0.7, \n                         color=colors_map[name], edgecolor='black')\n        axes[0, idx].set_xlabel('Sequence Length (nt)', fontsize=10)\n        axes[0, idx].set_ylabel('Frequency', fontsize=10)\n        axes[0, idx].set_title(f'{name.capitalize()} - Length Distribution', \n                              fontsize=11, fontweight='bold')\n        axes[0, idx].axvline(df['seq_length'].mean(), color='red', \n                            linestyle='--', linewidth=2, \n                            label=f\"Mean: {df['seq_length'].mean():.0f}\")\n        axes[0, idx].legend()\n\n# Row 2: GC content\nfor idx, (name, df) in enumerate(seq_datasets.items()):\n    if df is not None and 'gc_content' in df.columns:\n        axes[1, idx].hist(df['gc_content'], bins=50, alpha=0.7, \n                         color=colors_map[name], edgecolor='black')\n        axes[1, idx].set_xlabel('GC Content (%)', fontsize=10)\n        axes[1, idx].set_ylabel('Frequency', fontsize=10)\n        axes[1, idx].set_title(f'{name.capitalize()} - GC Content', \n                              fontsize=11, fontweight='bold')\n        axes[1, idx].axvline(df['gc_content'].mean(), color='red', \n                            linestyle='--', linewidth=2,\n                            label=f\"Mean: {df['gc_content'].mean():.1f}%\")\n        axes[1, idx].legend()\n\nplt.tight_layout()\nplt.savefig(f\"{CONFIG['OUTPUT_DIR']}/sequence_distributions.png\", \n           dpi=CONFIG['DPI'], bbox_inches='tight')\nplt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-02-06T18:32:32.539775Z","iopub.execute_input":"2026-02-06T18:32:32.540158Z","iopub.status.idle":"2026-02-06T18:32:35.484435Z","shell.execute_reply.started":"2026-02-06T18:32:32.540121Z","shell.execute_reply":"2026-02-06T18:32:35.482666Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Compare distributions between splits\nprint(\"\\n📊 Statistical Comparison of Splits:\")\nprint(\"=\"*70)\n\nif 'train' in seq_datasets and 'validation' in seq_datasets:\n    train_df = seq_datasets['train']\n    val_df = seq_datasets['validation']\n    \n    # KS test for sequence length\n    from scipy.stats import ks_2samp\n    \n    stat, pval = ks_2samp(train_df['seq_length'], val_df['seq_length'])\n    print(f\"\\nSequence Length Distribution:\")\n    print(f\"  KS statistic: {stat:.4f}\")\n    print(f\"  P-value: {pval:.4e}\")\n    print(f\"  Distributions are {'SIMILAR' if pval > 0.05 else 'DIFFERENT'} (α=0.05)\")\n    \n    # KS test for GC content\n    stat, pval = ks_2samp(train_df['gc_content'].dropna(), val_df['gc_content'].dropna())\n    print(f\"\\nGC Content Distribution:\")\n    print(f\"  KS statistic: {stat:.4f}\")\n    print(f\"  P-value: {pval:.4e}\")\n    print(f\"  Distributions are {'SIMILAR' if pval > 0.05 else 'DIFFERENT'} (α=0.05)\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-02-06T18:32:35.485732Z","iopub.execute_input":"2026-02-06T18:32:35.486052Z","iopub.status.idle":"2026-02-06T18:32:35.501756Z","shell.execute_reply.started":"2026-02-06T18:32:35.486023Z","shell.execute_reply":"2026-02-06T18:32:35.500358Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Nucleotide composition comparison\nprint(\"\\n🔬 Nucleotide Composition Analysis:\")\nprint(\"=\"*70)\n\nfig, axes = plt.subplots(1, len(seq_datasets), figsize=(16, 5))\n\nfor idx, (name, df) in enumerate(seq_datasets.items()):\n    if df is not None and all(nuc in df.columns for nuc in CONFIG['NUCLEOTIDES']):\n        avg_comp = df[CONFIG['NUCLEOTIDES']].mean()\n        \n        ax = axes[idx] if len(seq_datasets) > 1 else axes\n        bars = ax.bar(CONFIG['NUCLEOTIDES'], avg_comp.values, \n                     color=['#FF6B6B', '#4ECDC4', '#95E1D3', '#FFE66D'],\n                     alpha=0.8, edgecolor='black', linewidth=1.5)\n        ax.set_ylabel('Average Percentage (%)', fontsize=11)\n        ax.set_title(f'{name.capitalize()} - Nucleotide Composition', \n                    fontsize=12, fontweight='bold')\n        ax.set_ylim([0, max(avg_comp.values) * 1.2])\n        \n        # Add value labels\n        for i, (bar, val) in enumerate(zip(bars, avg_comp.values)):\n            height = bar.get_height()\n            ax.text(bar.get_x() + bar.get_width()/2., height,\n                   f'{val:.1f}%', ha='center', va='bottom', \n                   fontsize=10, fontweight='bold')\n        \n        print(f\"\\n{name.upper()}:\")\n        for nuc, pct in avg_comp.items():\n            print(f\"  {nuc}: {pct:.2f}%\")\n\nplt.tight_layout()\nplt.savefig(f\"{CONFIG['OUTPUT_DIR']}/nucleotide_composition.png\", \n           dpi=CONFIG['DPI'], bbox_inches='tight')\nplt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-02-06T18:32:35.50313Z","iopub.execute_input":"2026-02-06T18:32:35.503595Z","iopub.status.idle":"2026-02-06T18:32:36.459592Z","shell.execute_reply.started":"2026-02-06T18:32:35.503557Z","shell.execute_reply":"2026-02-06T18:32:36.458323Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"---\n## 5. Label Analysis (Training Structures)","metadata":{}},{"cell_type":"code","source":"print_section(\"🎯 LABEL ANALYSIS\")\n\n# Analyze train labels\nif 'train_labels' in dataframes:\n    labels_df = dataframes['train_labels']\n    \n    print(\"Train Labels Information:\")\n    print(\"=\"*60)\n    print(f\"Shape: {labels_df.shape}\")\n    print(f\"\\nColumns: {labels_df.columns.tolist()}\")\n    print(f\"\\nFirst few rows:\")\n    display(labels_df.head())\n    \n    # Check if labels are 3D coordinates or references\n    print(f\"\\nLabel format analysis:\")\n    for col in labels_df.columns:\n        print(f\"  {col}: {labels_df[col].dtype}\")\n        if labels_df[col].dtype in ['float64', 'int64']:\n            print(f\"    Range: [{labels_df[col].min():.2f}, {labels_df[col].max():.2f}]\")\n        else:\n            print(f\"    Unique values: {labels_df[col].nunique()}\")\n            if labels_df[col].nunique() < 10:\n                print(f\"    Values: {labels_df[col].unique().tolist()}\")\nelse:\n    print(\"⚠️ No train_labels data found\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-02-06T18:32:36.460897Z","iopub.execute_input":"2026-02-06T18:32:36.461228Z","iopub.status.idle":"2026-02-06T18:32:48.098517Z","shell.execute_reply.started":"2026-02-06T18:32:36.461193Z","shell.execute_reply":"2026-02-06T18:32:48.097414Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"---\n## 6. PDB Structure Analysis","metadata":{}},{"cell_type":"code","source":"!pip install -q biopython","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# --- Section 6: PDB Structure Analysis ---\nfrom Bio.PDB.MMCIFParser import MMCIFParser\nimport numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\nfrom tqdm.notebook import tqdm\n\nprint_section(\"🏗️ PDB STRUCTURE ANALYSIS\")\n\n# 1. Improved Parser for .cif files\ndef parse_cif_c1_coords(cif_file):\n    \"\"\"Accurately extract C1' coordinates using BioPython MMCIFParser\"\"\"\n    parser = MMCIFParser(QUIET=True)\n    coords = []\n    try:\n        structure = parser.get_structure('RNA', str(cif_file))\n        for model in structure:\n            for chain in model:\n                for residue in chain:\n                    if \"C1'\" in residue:\n                        coords.append(residue[\"C1'\"].get_coord())\n    except:\n        pass\n    return np.array(coords)\n\n# 2. File Inventory\npdb_dir = data_path / 'PDB_RNA'\npdb_files = sorted(list(pdb_dir.glob('*.cif')))\nprint(f\"Found {len(pdb_files)} CIF files in directory.\")\n\nstructure_stats = []\n\n# 3. Processing (Setting to 500 for a representative sample, change to len(pdb_files) for full run)\n#process_limit = 500 \nfor pdb_file in tqdm(pdb_files, desc=\"Processing Structures\"):\n    coords = parse_cif_c1_coords(pdb_file)\n    \n    if len(coords) > 0:\n        # Calculate metrics for visualization\n        rg = calculate_radius_of_gyration(coords)\n        centroid = np.mean(coords, axis=0)\n        \n        # Bounding Box calculations\n        bbox_min = np.min(coords, axis=0)\n        bbox_max = np.max(coords, axis=0)\n        bbox_size = bbox_max - bbox_min\n        \n        # Max pairwise distance (Approximate for speed)\n        max_dist = np.linalg.norm(bbox_max - bbox_min)\n        \n        structure_stats.append({\n            'pdb_id': pdb_file.stem,\n            'num_residues': len(coords),\n            'radius_of_gyration': rg,\n            'max_distance': max_dist,\n            'bbox_x': bbox_size[0],\n            'bbox_y': bbox_size[1],\n            'bbox_z': bbox_size[2],\n            'compactness': len(coords) / (rg + 1e-6)\n        })\n\n# 4. Visualization and Data Export\nif structure_stats:\n    # Important: define struct_df so Section 10 can find it\n    struct_df = pd.DataFrame(structure_stats)\n    \n    # Create 2x3 Plotting Area\n    fig, axes = plt.subplots(2, 3, figsize=(18, 12))\n    \n    # Plot 1: Size Distribution\n    axes[0, 0].hist(struct_df['num_residues'], bins=30, color='skyblue', edgecolor='black', alpha=0.7)\n    axes[0, 0].set_title('Structure Size Distribution', fontweight='bold')\n    axes[0, 0].set_xlabel('Number of Residues')\n    axes[0, 0].set_ylabel('Frequency')\n\n    # Plot 2: Radius of Gyration\n    axes[0, 1].hist(struct_df['radius_of_gyration'], bins=30, color='coral', edgecolor='black', alpha=0.7)\n    axes[0, 1].set_title('Radius of Gyration (Rg) Distribution', fontweight='bold')\n    axes[0, 1].set_xlabel('Rg (Å)')\n    axes[0, 1].set_ylabel('Frequency')\n\n    # Plot 3: Compactness vs Size\n    axes[0, 2].scatter(struct_df['num_residues'], struct_df['compactness'], color='green', alpha=0.5)\n    axes[0, 2].set_title('Compactness vs Structure Size', fontweight='bold')\n    axes[0, 2].set_xlabel('Number of Residues')\n    axes[0, 2].set_ylabel('Compactness (N/Rg)')\n\n    # Plot 4: Size vs Rg (Scaling Law)\n    axes[1, 0].scatter(struct_df['num_residues'], struct_df['radius_of_gyration'], color='purple', alpha=0.5)\n    axes[1, 0].set_title('Size vs Radius of Gyration', fontweight='bold')\n    axes[1, 0].set_xlabel('Number of Residues')\n    axes[1, 0].set_ylabel('Rg (Å)')\n\n    # Plot 5: Bounding Box dimensions\n    bbox_data = [struct_df['bbox_x'], struct_df['bbox_y'], struct_df['bbox_z']]\n    axes[1, 1].boxplot(bbox_data, labels=['X', 'Y', 'Z'], patch_artist=True)\n    axes[1, 1].set_title('Bounding Box Dimensions', fontweight='bold')\n    axes[1, 1].set_ylabel('Size (Å)')\n\n    # Plot 6: Summary Stats Table (Simplified)\n    axes[1, 2].axis('off')\n    summary_text = (\n        f\"Analyzed Samples: {len(struct_df)}\\n\\n\"\n        f\"Avg residues: {struct_df['num_residues'].mean():.1f}\\n\"\n        f\"Max residues: {struct_df['num_residues'].max()}\\n\"\n        f\"Avg Rg: {struct_df['radius_of_gyration'].mean():.2f} Å\\n\"\n        f\"Avg Compactness: {struct_df['compactness'].mean():.2f}\"\n    )\n    axes[1, 2].text(0.1, 0.5, summary_text, fontsize=12, family='monospace', bbox=dict(facecolor='white', alpha=0.5))\n\n    plt.tight_layout()\n    plt.savefig(f\"{CONFIG['OUTPUT_DIR']}/pdb_structure_analysis_v2.png\", dpi=CONFIG['DPI'])\n    plt.show()\n\n    # Save to Cache for Section 10\n    struct_df.to_csv(f\"{CONFIG['CACHE_DIR']}/structure_stats.csv\", index=False)\n    print(f\"✅ Successfully analyzed {len(struct_df)} structures and saved plots.\")\n    display(struct_df.head())\nelse:\n    print(\"⚠️ No valid PDB structures could be analyzed. Check file formats.\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-02-06T18:50:59.765897Z","iopub.execute_input":"2026-02-06T18:50:59.766225Z","iopub.status.idle":"2026-02-06T18:56:05.077379Z","shell.execute_reply.started":"2026-02-06T18:50:59.7662Z","shell.execute_reply":"2026-02-06T18:56:05.07653Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"---\n## 7. MSA (Multiple Sequence Alignment) Analysis","metadata":{}},{"cell_type":"code","source":"# --- Section 7: MSA (Multiple Sequence Alignment) Analysis ---\nfrom Bio import SeqIO\nimport numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\nfrom tqdm.notebook import tqdm\n\nprint_section(\"🧬 MSA (Multiple Sequence Alignment) ANALYSIS\")\n\n# 1. Directory and File Setup\nmsa_dir = data_path / 'MSA'\nif msa_dir.exists():\n    msa_files = sorted(list(msa_dir.glob('*.fasta')))\n    print(f\"Found {len(msa_files)} MSA files in {msa_dir}\")\n\n    # 2. Analyzing Sample MSA files\n   \n   \n    msa_stats = []\n    \n    print(f\"Analyzing MSA files...\")\n    for msa_file in tqdm(msa_files, desc=\"Processing MSAs\"):\n        try:\n            # Parse FASTA alignment\n            alignment = list(SeqIO.parse(msa_file, \"fasta\"))\n            \n            if len(alignment) > 0:\n                num_seqs = len(alignment)\n                align_len = len(alignment[0].seq)\n                \n                # Calculate gap statistics\n                all_seqs_str = [str(record.seq) for record in alignment]\n                gap_counts = [s.count('-') + s.count('.') for s in all_seqs_str]\n                gap_rates = [c / align_len for c in gap_counts]\n                \n                msa_stats.append({\n                    'msa_id': msa_file.stem,\n                    'num_sequences': num_seqs, # Depth of MSA\n                    'alignment_length': align_len,\n                    'avg_gap_rate': np.mean(gap_rates),\n                    'max_gap_rate': np.max(gap_rates)\n                })\n        except Exception as e:\n            # print(f\"Error parsing {msa_file.name}: {e}\")\n            pass\n\n    # 3. Visualization\n    if msa_stats:\n        msa_df = pd.DataFrame(msa_stats)\n        \n        fig, axes = plt.subplots(1, 3, figsize=(18, 5))\n        \n        # Plot 1: MSA Depth (Number of Sequences)\n        axes[0].hist(msa_df['num_sequences'], bins=20, color='skyblue', edgecolor='black')\n        axes[0].set_title('MSA Depth Distribution', fontweight='bold')\n        axes[0].set_xlabel('Number of Sequences')\n        axes[0].set_ylabel('Frequency')\n        axes[0].set_yscale('log') # Log scale as depth varies greatly\n\n        # Plot 2: Alignment Length\n        axes[1].hist(msa_df['alignment_length'], bins=20, color='lightgreen', edgecolor='black')\n        axes[1].set_title('Alignment Length Distribution', fontweight='bold')\n        axes[1].set_xlabel('Length (nt)')\n        axes[1].set_ylabel('Frequency')\n\n        # Plot 3: Gap Rate\n        axes[2].hist(msa_df['avg_gap_rate'], bins=20, color='salmon', edgecolor='black')\n        axes[2].set_title('Average Gap Rate Distribution', fontweight='bold')\n        axes[2].set_xlabel('Gap Fraction')\n        axes[2].set_ylabel('Frequency')\n\n        plt.tight_layout()\n        plt.show()\n\n        # Save to Cache for final summary\n        msa_df.to_csv(f\"{CONFIG['CACHE_DIR']}/msa_stats.csv\", index=False)\n        print(f\"✅ MSA Analysis complete. Sampled average depth: {msa_df['num_sequences'].mean():.1f}\")\n        display(msa_df.head())\n    else:\n        print(\"❌ Could not extract stats from MSA files. Check file integrity.\")\nelse:\n    print(\"⚠️ MSA directory not found. Please check your data path.\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-02-06T18:57:04.599587Z","iopub.execute_input":"2026-02-06T18:57:04.599927Z","iopub.status.idle":"2026-02-06T18:58:37.009456Z","shell.execute_reply.started":"2026-02-06T18:57:04.599899Z","shell.execute_reply":"2026-02-06T18:58:37.00838Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"---\n## 8. Metadata Analysis","metadata":{}},{"cell_type":"code","source":"print_section(\"📋 METADATA ANALYSIS\")\n\nif 'rna_metadata' in dataframes:\n    metadata_df = dataframes['rna_metadata']\n    \n    print(\"RNA Metadata Overview:\")\n    print(\"=\"*70)\n    print(f\"Shape: {metadata_df.shape}\")\n    print(f\"\\nColumns: {metadata_df.columns.tolist()}\")\n    print(f\"\\nFirst few rows:\")\n    display(metadata_df.head(10))\n    \n    print(f\"\\nColumn Statistics:\")\n    for col in metadata_df.columns:\n        print(f\"\\n{col}:\")\n        print(f\"  Type: {metadata_df[col].dtype}\")\n        print(f\"  Missing: {metadata_df[col].isnull().sum()} ({metadata_df[col].isnull().sum()/len(metadata_df)*100:.1f}%)\")\n        \n        if metadata_df[col].dtype in ['float64', 'int64']:\n            print(f\"  Range: [{metadata_df[col].min()}, {metadata_df[col].max()}]\")\n            print(f\"  Mean: {metadata_df[col].mean():.2f}\")\n        else:\n            nunique = metadata_df[col].nunique()\n            print(f\"  Unique values: {nunique}\")\n            if nunique < 20:\n                print(f\"  Values: {metadata_df[col].value_counts().to_dict()}\")\nelse:\n    print(\"⚠️ No metadata file found\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-02-06T18:32:48.361873Z","iopub.execute_input":"2026-02-06T18:32:48.362527Z","iopub.status.idle":"2026-02-06T18:32:48.888746Z","shell.execute_reply.started":"2026-02-06T18:32:48.362489Z","shell.execute_reply":"2026-02-06T18:32:48.886632Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"---\n## 9. Sample Submission Analysis","metadata":{}},{"cell_type":"code","source":"print_section(\"📝 SAMPLE SUBMISSION FORMAT\")\n\nif 'sample_submission' in dataframes:\n    submission_df = dataframes['sample_submission']\n    \n    print(\"Sample Submission Structure:\")\n    print(\"=\"*70)\n    print(f\"Shape: {submission_df.shape}\")\n    print(f\"Columns: {submission_df.columns.tolist()}\")\n    print(f\"\\nFirst few rows:\")\n    display(submission_df.head())\n    \n    print(\"\\n⚠️ IMPORTANT SUBMISSION NOTES:\")\n    print(\"=\"*70)\n    print(\"• Competition requires predicting 5 models per sequence\")\n    print(\"• Output format: C1' coordinates for each residue\")\n    print(\"• Evaluation metric: TM-score (Template Modeling score)\")\n    print(\"• Best of 5 predictions will be used for scoring\")\nelse:\n    print(\"⚠️ No sample submission file found\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-02-06T18:32:48.890718Z","iopub.execute_input":"2026-02-06T18:32:48.891339Z","iopub.status.idle":"2026-02-06T18:32:48.912197Z","shell.execute_reply.started":"2026-02-06T18:32:48.891148Z","shell.execute_reply":"2026-02-06T18:32:48.910111Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"---\n## 10. Summary and Key Insights","metadata":{}},{"cell_type":"code","source":"print_section(\"📝 COMPREHENSIVE SUMMARY\", char='*')\n\nsummary = {\n    'Dataset Overview': {\n        'Total Samples': sum(split_info.values()),\n        'Training Samples': split_info['Train'],\n        'Validation Samples': split_info['Validation'],\n        'Test Samples': split_info['Test'],\n    },\n    'Data Files': {\n        'CSV Files': len([k for k in dataframes.keys()]),\n        'PDB Structures': len(list((data_path / 'PDB_RNA').glob('*.pdb'))) if (data_path / 'PDB_RNA').exists() else 0,\n        'MSA Files': len(list((data_path / 'MSA').glob('*'))) if (data_path / 'MSA').exists() else 0,\n    },\n}\n\n# Add sequence statistics if available\nif 'train' in seq_datasets and seq_datasets['train'] is not None:\n    train_df = seq_datasets['train']\n    summary['Sequence Characteristics'] = {\n        'Avg Length (train)': f\"{train_df['seq_length'].mean():.1f} nt\",\n        'Length Range (train)': f\"[{train_df['seq_length'].min()}, {train_df['seq_length'].max()}] nt\",\n        'Avg GC Content (train)': f\"{train_df['gc_content'].mean():.2f}%\",\n    }\n\n# Add structure statistics if available\nif 'struct_df' in locals() and struct_df is not None:\n    summary['Structure Properties'] = {\n        'Structures Analyzed': len(struct_df),\n        'Avg Radius of Gyration': f\"{struct_df['radius_of_gyration'].mean():.2f} Å\",\n        'Avg Residues per Structure': f\"{struct_df['num_residues'].mean():.1f}\",\n    }\n\nprint(json.dumps(summary, indent=2))\n\nprint(\"\\n\" + \"=\"*80)\nprint(\"KEY INSIGHTS\")\nprint(\"=\"*80)\n\ninsights = [\n    \"\\n1. DATA STRUCTURE:\",\n    \"   • Well-organized train/validation/test split\",\n    \"   • Separate sequence and label files\",\n    \"   • Rich auxiliary data (MSA, PDB structures, metadata)\",\n    \n    \"\\n2. SEQUENCE DIVERSITY:\",\n    \"   • Wide range of sequence lengths\",\n    \"   • Varied GC content across samples\",\n    \"   • Balanced nucleotide composition\",\n    \n    \"\\n3. STRUCTURAL DATA:\",\n    \"   • PDB files contain C1' coordinates for backbone\",\n    \"   • Structures show varying compactness\",\n    \"   • Need to predict 5 models per sequence\",\n    \n    \"\\n4. MSA AVAILABILITY:\",\n    \"   • Multiple sequence alignments available\",\n    \"   • Critical for evolutionary information\",\n    \"   • Variable depth across samples\",\n    \n    \"\\n5. MODELING CONSIDERATIONS:\",\n    \"   • TM-score evaluation (structure similarity)\",\n    \"   • Need ensemble approach (5 predictions)\",\n    \"   • Can leverage MSA features\",\n    \"   • May benefit from transfer learning\",\n]\n\nfor insight in insights:\n    print(insight)\n\nprint(\"\\n\" + \"=\"*80)\nprint(\"NEXT STEPS\")\nprint(\"=\"*80)\n\nnext_steps = [\n    \"\\n1. DATA PREPROCESSING:\",\n    \"   → Parse PDB files to extract C1' coordinates\",\n    \"   → Process MSA files for coevolution features\",\n    \"   → Normalize sequences and coordinates\",\n    \"   → Create efficient data loaders\",\n    \n    \"\\n2. FEATURE ENGINEERING:\",\n    \"   → MSA features (coevolution, conservation)\",\n    \"   → Secondary structure prediction\",\n    \"   → Distance maps and contact predictions\",\n    \"   → RNA language model embeddings\",\n    \n    \"\\n3. MODEL DEVELOPMENT:\",\n    \"   → Baseline: Template-based methods\",\n    \"   → Deep learning: Transformer/GNN architectures\",\n    \"   → Ensemble: Multiple model averaging\",\n    \"   → Refinement: Energy minimization\",\n    \n    \"\\n4. EVALUATION:\",\n    \"   → Implement TM-score calculator\",\n    \"   → Cross-validation strategy\",\n    \"   → Visualization tools (PyMOL, py3Dmol)\",\n    \"   → Track multiple metrics (RMSD, GDT, etc.)\",\n]\n\nfor step in next_steps:\n    print(step)\n\nprint(\"\\n\" + \"*\"*80)\nprint(\"EDA COMPLETE! 🎉\")\nprint(\"*\"*80)\nprint(f\"\\nOutputs saved to: {CONFIG['OUTPUT_DIR']}\")\nprint(f\"Cache directory: {CONFIG['CACHE_DIR']}\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-02-06T18:32:48.913879Z","iopub.execute_input":"2026-02-06T18:32:48.914272Z","iopub.status.idle":"2026-02-06T18:32:48.974596Z","shell.execute_reply.started":"2026-02-06T18:32:48.914211Z","shell.execute_reply":"2026-02-06T18:32:48.9735Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Save all enhanced datasets\nprint(\"💾 Saving processed data...\\n\")\n\nfor name, df in seq_datasets.items():\n    if df is not None:\n        output_file = f\"{CONFIG['CACHE_DIR']}/{name}_sequences_enhanced.csv\"\n        df.to_csv(output_file, index=False)\n        print(f\"✅ Saved: {output_file}\")\n\n# Save summary\nwith open(f\"{CONFIG['OUTPUT_DIR']}/eda_summary.json\", 'w') as f:\n    json.dump(summary, f, indent=2)\nprint(f\"✅ Saved: {CONFIG['OUTPUT_DIR']}/eda_summary.json\")\n\nprint(\"\\n🎉 All analysis complete!\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-02-06T19:03:37.641432Z","iopub.execute_input":"2026-02-06T19:03:37.642314Z","iopub.status.idle":"2026-02-06T19:03:38.677359Z","shell.execute_reply.started":"2026-02-06T19:03:37.64223Z","shell.execute_reply":"2026-02-06T19:03:38.676324Z"}},"outputs":[],"execution_count":null}]}