{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","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":31234,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"!pip install -U plotly\n!pip install py3Dmol\n!pip install biopython","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-09T05:20:29.851491Z","iopub.execute_input":"2026-01-09T05:20:29.851865Z","iopub.status.idle":"2026-01-09T05:21:04.639472Z","shell.execute_reply.started":"2026-01-09T05:20:29.851832Z","shell.execute_reply":"2026-01-09T05:21:04.638038Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import pandas as pd\nimport numpy as np\nimport matplotlib.pyplot as plt\nimport seaborn as sns\nfrom collections import Counter\nimport plotly.express as px\nimport plotly.graph_objects as go\nfrom plotly.subplots import make_subplots\nimport warnings\nwarnings.filterwarnings('ignore')\n\n# Set style for better visualizations\nplt.style.use('default')\nsns.set_palette(\"husl\")\n","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true,"execution":{"iopub.status.busy":"2026-01-09T05:21:04.641964Z","iopub.execute_input":"2026-01-09T05:21:04.642424Z","iopub.status.idle":"2026-01-09T05:21:07.378849Z","shell.execute_reply.started":"2026-01-09T05:21:04.64239Z","shell.execute_reply":"2026-01-09T05:21:07.377865Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Load the main datasets\ntrain_seq = pd.read_csv('/kaggle/input/stanford-rna-3d-folding-2/train_sequences.csv')\ntrain_labels = pd.read_csv('/kaggle/input/stanford-rna-3d-folding-2/train_labels.csv')\nval_seq = pd.read_csv('/kaggle/input/stanford-rna-3d-folding-2/validation_sequences.csv')\n\nprint(\"Dataset Shapes:\")\nprint(f\"Train Sequences: {train_seq.shape}\")\nprint(f\"Train Labels: {train_labels.shape}\")\nprint(f\"Validation Sequences: {val_seq.shape}\")\nprint(\"\\nTrain Sequences Info:\")\nprint(train_seq.info())\nprint(\"\\nTrain Labels Info:\")\nprint(train_labels.info())\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-09T05:21:07.380122Z","iopub.execute_input":"2026-01-09T05:21:07.380763Z","iopub.status.idle":"2026-01-09T05:21:19.069573Z","shell.execute_reply.started":"2026-01-09T05:21:07.380724Z","shell.execute_reply":"2026-01-09T05:21:19.068644Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## Sequence Length Distribution Analysis","metadata":{}},{"cell_type":"code","source":"# Sequence length distribution\ntrain_seq['sequence_length'] = train_seq['sequence'].str.len()\nval_seq['sequence_length'] = val_seq['sequence'].str.len()\n\nfig = make_subplots(rows=2, cols=2, \n                   subplot_titles=('Sequence Length Distribution', 'Length Histogram', \n                                 'Length Box Plot', 'Cumulative Distribution'),\n                   specs=[[{\"secondary_y\": False}, {\"secondary_y\": False}],\n                         [{\"secondary_y\": False}, {\"secondary_y\": False}]])\n\n# Length distribution\nfig.add_trace(go.Histogram(x=train_seq['sequence_length'], name='Train', \n                          nbinsx=50, opacity=0.7), row=1, col=1)\nfig.add_trace(go.Histogram(x=val_seq['sequence_length'], name='Validation', \n                          nbinsx=50, opacity=0.7), row=1, col=1)\n\n# Length histogram\nfig.add_trace(go.Histogram(x=train_seq['sequence_length'], name='Train Hist', \n                          nbinsx=100), row=1, col=2)\nfig.add_trace(go.Histogram(x=val_seq['sequence_length'], name='Val Hist', \n                          nbinsx=100), row=1, col=2)\n\n# Box plot\nfig.add_trace(go.Box(y=train_seq['sequence_length'], name='Train Box'), row=2, col=1)\nfig.add_trace(go.Box(y=val_seq['sequence_length'], name='Validation Box'), row=2, col=1)\n\n# CDF\ntrain_sorted = np.sort(train_seq['sequence_length'])\nval_sorted = np.sort(val_seq['sequence_length'])\nfig.add_trace(go.Scatter(x=np.arange(len(train_sorted))/len(train_sorted), \n                        y=train_sorted, name='Train CDF', line=dict(color='blue')), row=2, col=2)\nfig.add_trace(go.Scatter(x=np.arange(len(val_sorted))/len(val_sorted), \n                        y=val_sorted, name='Val CDF', line=dict(color='red')), row=2, col=2)\n\nfig.update_layout(height=800, title_text=\"Sequence Length Analysis\", showlegend=True)\nfig.show()\n\nprint(f\"Sequence Length Stats:\")\nprint(train_seq['sequence_length'].describe())\nprint(f\"\\nOutliers (>Q3+1.5*IQR): {len(train_seq[train_seq['sequence_length'] > train_seq['sequence_length'].quantile(0.75) + 1.5*(train_seq['sequence_length'].quantile(0.75) - train_seq['sequence_length'].quantile(0.25))])}\")\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-09T05:21:19.071456Z","iopub.execute_input":"2026-01-09T05:21:19.07173Z","iopub.status.idle":"2026-01-09T05:21:20.99467Z","shell.execute_reply.started":"2026-01-09T05:21:19.071707Z","shell.execute_reply":"2026-01-09T05:21:20.993853Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## Nucleotide Composition Heatmap / Analysis","metadata":{}},{"cell_type":"code","source":"def nucleotide_composition(seq_df):\n    \"\"\"Calculate nucleotide frequencies for each sequence\"\"\"\n    comp = []\n    for seq in seq_df['sequence']:\n        counts = Counter(seq)\n        total = len(seq)\n        comp.append({\n            'A_freq': counts.get('A', 0) / total,\n            'C_freq': counts.get('C', 0) / total,\n            'G_freq': counts.get('G', 0) / total,\n            'U_freq': counts.get('U', 0) / total,\n            'GC_content': (counts.get('G', 0) + counts.get('C', 0)) / total,\n            'AU_content': (counts.get('A', 0) + counts.get('U', 0)) / total\n        })\n    return pd.DataFrame(comp)\n\n# Calculate composition\ntrain_comp = nucleotide_composition(train_seq)\nval_comp = nucleotide_composition(val_seq)\n\n# Combine for correlation analysis\ncomp_df = pd.concat([train_comp.assign(set='train'), val_comp.assign(set='validation')])\n\n# Composition heatmap\nfig = px.imshow(comp_df[['A_freq', 'C_freq', 'G_freq', 'U_freq', 'GC_content']].corr(),\n               title=\"Nucleotide Composition Correlation Matrix\",\n               color_continuous_scale='RdBu_r', aspect=\"auto\")\nfig.show()\n\n# Composition distribution\nfig, axes = plt.subplots(2, 3, figsize=(18, 12))\naxes = axes.ravel()\n\nfor i, col in enumerate(['A_freq', 'C_freq', 'G_freq', 'U_freq', 'GC_content', 'AU_content']):\n    sns.histplot(data=comp_df, x=col, hue='set', kde=True, ax=axes[i], alpha=0.6)\n    axes[i].set_title(f'{col.replace(\"_\", \" \").title()} Distribution')\n    axes[i].grid(True, alpha=0.3)\n\nplt.tight_layout()\nplt.show()\n\nprint(\"Composition Statistics:\")\nprint(comp_df[['A_freq', 'C_freq', 'G_freq', 'U_freq', 'GC_content']].describe())\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-09T05:21:20.995835Z","iopub.execute_input":"2026-01-09T05:21:20.996156Z","iopub.status.idle":"2026-01-09T05:21:24.945716Z","shell.execute_reply.started":"2026-01-09T05:21:20.996123Z","shell.execute_reply":"2026-01-09T05:21:24.944765Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## Stoichiometric Analysis","metadata":{}},{"cell_type":"code","source":"# Stoichiometry parsing (ensure these columns exist first)\ntrain_seq['num_chains'] = train_seq['stoichiometry'].apply(\n    lambda x: len(x.split(';')) if pd.notna(x) and x else 1\n)\ntrain_seq['chain_pattern'] = train_seq['stoichiometry'].fillna('single').apply(\n    lambda x: x.split(';')[0] if x and ';' in x else 'single'\n)\n\ntrain_seq['has_ligands'] = train_seq['ligand_ids'].notna() & (train_seq['ligand_ids'] != '')\ntrain_seq['num_ligands'] = train_seq['ligand_ids'].apply(\n    lambda x: len(x.split(';')) if pd.notna(x) and x and x != '' else 0\n)\n\nfig = make_subplots(\n    rows=2, cols=2, \n    subplot_titles=('Chain Distribution', 'Ligand Presence', \n                    'Chains vs Length', 'Ligands vs GC Content'),\n    specs=[[{\"type\": \"histogram\"}, {\"type\": \"pie\"}],\n           [{\"type\": \"scatter\"}, {\"type\": \"scatter\"}]]\n)\n\n# Chain distribution histogram\nfig.add_trace(\n    go.Histogram(x=train_seq['num_chains'], name='Num Chains', \n                nbinsx=10, marker_color='skyblue'),\n    row=1, col=1\n)\n\nligand_counts = train_seq['has_ligands'].value_counts()\nfig.add_trace(\n    go.Pie(\n        values=ligand_counts.values,\n        labels=['No Ligands', 'Has Ligands'],\n        marker_colors=['lightgray', '#FF6B6B'],\n        textinfo='label+percent',\n        hole=0.4,  # Donut style\n        showlegend=True\n    ),\n    row=1, col=2\n)\n\n# Chains vs Length scatter (with GC content coloring)\nfig.add_trace(\n    go.Scatter(\n        x=train_seq['num_chains'], \n        y=train_seq['sequence_length'],\n        mode='markers',\n        marker=dict(\n            color=train_comp['GC_content']*100, \n            colorscale='Viridis',\n            size=8,\n            opacity=0.7,\n            colorbar=dict(title=\"GC Content (%)\")\n        ),\n        name='Chains vs Length',\n        text=[f\"ID: {tid[:10]}...\" for tid in train_seq['target_id']],\n        hovertemplate='<b>%{text}</b><br>Chains: %{x}<br>Length: %{y}<br>GC: %{marker.color:.1f}%<extra></extra>'\n    ),\n    row=2, col=1\n)\n\n# Ligands vs GC scatter\nfig.add_trace(\n    go.Scatter(\n        x=train_seq['num_ligands'], \n        y=train_comp['GC_content'],\n        mode='markers',\n        marker=dict(color='coral', size=10, line=dict(color='white', width=1)),\n        name='Ligands vs GC',\n        hovertemplate='<b>Ligands: %{x}</b><br>GC Content: %{y:.1%}<extra></extra>'\n    ),\n    row=2, col=2\n)\n\nfig.update_layout(\n    height=800, \n    title_text=\"🔬 RNA Structural Context Analysis\",\n    showlegend=True,\n    font=dict(size=12)\n)\nfig.show()\n\n# Print statistics\nprint(\"🧬 Structural Features Summary:\")\nprint(f\"• Multi-chain structures: {sum(train_seq['num_chains'] > 1)} ({100*sum(train_seq['num_chains'] > 1)/len(train_seq):.1f}%)\")\nprint(f\"• Structures with ligands: {train_seq['has_ligands'].sum()} ({100*train_seq['has_ligands'].mean():.1f}%)\")\nprint(f\"• Max ligands per structure: {train_seq['num_ligands'].max()}\")\nprint(f\"• Chain distribution: {train_seq['num_chains'].value_counts().to_dict()}\")\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-09T05:21:24.947011Z","iopub.execute_input":"2026-01-09T05:21:24.947382Z","iopub.status.idle":"2026-01-09T05:21:25.089671Z","shell.execute_reply.started":"2026-01-09T05:21:24.947355Z","shell.execute_reply":"2026-01-09T05:21:25.088496Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## Publication Trend Analysis","metadata":{}},{"cell_type":"code","source":"# Parse temporal cutoff dates\ntrain_seq['temporal_cutoff'] = pd.to_datetime(train_seq['temporal_cutoff'])\nval_seq['temporal_cutoff'] = pd.to_datetime(val_seq['temporal_cutoff'])\n\nfig, (ax1, ax2) = plt.subplots(2, 1, figsize=(15, 10))\n\n# Publication timeline\ntrain_seq['year'] = train_seq['temporal_cutoff'].dt.year\nyear_counts = train_seq['year'].value_counts().sort_index()\nax1.bar(year_counts.index, year_counts.values, alpha=0.7, color='skyblue')\nax1.set_title('RNA Structure Publications Over Time')\nax1.set_xlabel('Year')\nax1.set_ylabel('Number of Structures')\nax1.grid(True, alpha=0.3)\n\n# Length trends over time\ntrend_data = train_seq.groupby('year')['sequence_length'].agg(['mean', 'median', 'count']).reset_index()\nax2.plot(trend_data['year'], trend_data['mean'], marker='o', label='Mean Length', linewidth=2)\nax2.plot(trend_data['year'], trend_data['median'], marker='s', label='Median Length', linewidth=2)\nax2.set_title('RNA Sequence Length Trends Over Time')\nax2.set_xlabel('Year')\nax2.set_ylabel('Sequence Length')\nax2.legend()\nax2.grid(True, alpha=0.3)\n\nplt.tight_layout()\nplt.show()\n\nprint(\"Publication Trends:\")\nprint(train_seq['year'].value_counts().sort_index())\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-09T05:21:25.090917Z","iopub.execute_input":"2026-01-09T05:21:25.091273Z","iopub.status.idle":"2026-01-09T05:21:25.58493Z","shell.execute_reply.started":"2026-01-09T05:21:25.091246Z","shell.execute_reply":"2026-01-09T05:21:25.583709Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"def dinucleotide_freq(seq):\n    \"\"\"Calculate dinucleotide frequencies\"\"\"\n    di_counts = Counter()\n    for i in range(len(seq)-1):\n        di_counts[seq[i:i+2]] += 1\n    total = sum(di_counts.values())\n    return {k: v/total for k, v in di_counts.items()}\n\n# Top dinucleotides\nall_di = []\nfor seq in train_seq['sequence'][:1000]:  # Sample for speed\n    all_di.extend(list(dinucleotide_freq(seq).keys()))\n\ndi_counter = Counter(all_di)\ntop_di = di_counter.most_common(16)\n\nfig = go.Figure(data=[\n    go.Bar(x=[x[0] for x in top_di], y=[x[1] for x in top_di],\n           marker_color=['#1f77b4', '#ff7f0e', '#2ca02c', '#d62728', '#9467bd',\n                        '#8c564b', '#e377c2', '#7f7f7f', '#bcbd22', '#17becf']*2)\n])\nfig.update_layout(title=\"Top 16 Dinucleotide Frequencies\", \n                 xaxis_title=\"Dinucleotide\", yaxis_title=\"Frequency\")\nfig.show()\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-09T05:21:25.586132Z","iopub.execute_input":"2026-01-09T05:21:25.586731Z","iopub.status.idle":"2026-01-09T05:21:25.766932Z","shell.execute_reply.started":"2026-01-09T05:21:25.586702Z","shell.execute_reply":"2026-01-09T05:21:25.766141Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## Chain Distribution Analysis","metadata":{}},{"cell_type":"code","source":"chain_counts = train_labels['chain'].value_counts()\ncopy_counts = train_labels['copy'].value_counts()\n\nfig = make_subplots(rows=1, cols=2, subplot_titles=('Chain Distribution', 'Copy Number Distribution'))\n\nfig.add_trace(go.Bar(x=chain_counts.index, y=chain_counts.values, name='Chains'), row=1, col=1)\nfig.add_trace(go.Bar(x=copy_counts.index, y=copy_counts.values, name='Copies'), row=1, col=2)\n\nfig.update_layout(height=500, title_text=\"Chain and Copy Number Distributions\")\nfig.show()\n\nprint(\"Chain Statistics:\")\nprint(train_labels.groupby('chain').size().describe())\nprint(\"\\nCopy Statistics:\")\nprint(train_labels['copy'].describe())\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-09T05:21:25.768077Z","iopub.execute_input":"2026-01-09T05:21:25.768338Z","iopub.status.idle":"2026-01-09T05:21:27.145423Z","shell.execute_reply.started":"2026-01-09T05:21:25.768316Z","shell.execute_reply":"2026-01-09T05:21:27.144621Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## Coordinate Analysis","metadata":{}},{"cell_type":"code","source":"# 1. FIRST: Check data size and clean efficiently\nprint(\"Dataset Overview:\")\nprint(f\"Total rows: {len(train_labels):,}\")\nprint(f\"Coordinate columns: {[col for col in train_labels.columns if col.startswith(('x_', 'y_', 'z_'))][:5]}...\")\n\n# 2. Smart sampling - ONLY 5000 points max for 3D viz\nSAMPLE_SIZE = 5000\nvalid_coords = train_labels[['x_1', 'y_1', 'z_1', 'resid', 'resname', 'chain']].dropna()\nsample_coords = valid_coords.sample(n=min(SAMPLE_SIZE, len(valid_coords)), random_state=42)\nprint(f\" Using {len(sample_coords):,} valid points for 3D plot\")\n\n# 3. ULTRA-FAST 3D scatter (SIMPLE VERSION)\nfig = go.Figure()\n\nfig.add_trace(go.Scatter3d(\n    x=sample_coords['x_1'],\n    y=sample_coords['y_1'], \n    z=sample_coords['z_1'],\n    mode='markers',\n    marker=dict(\n        size=4,\n        color=sample_coords['resid'],\n        colorscale='Viridis',\n        opacity=0.6,\n        colorbar=dict(title=\"Residue #\")\n    ),\n    name='C1\\' Atoms',\n    hovertemplate='<b>Res %{customdata[0]}</b><br>Chain: %{customdata[1]}<br>X: %{x:.1f}Å<br>Y: %{y:.1f}Å<br>Z: %{z:.1f}Å<extra></extra>',\n    customdata=sample_coords[['resid', 'chain']].values\n))\n\nfig.update_layout(\n    title=f\"🧬 RNA C1' Backbone (Sample of {len(sample_coords):,} points)\",\n    scene=dict(\n        xaxis_title='X (Å)', yaxis_title='Y (Å)', zaxis_title='Z (Å)',\n        aspectmode='cube',\n        camera=dict(eye=dict(x=1.5, y=1.5, z=1.5))\n    ),\n    height=700,\n    showlegend=True\n)\nfig.show()\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-09T05:21:27.147773Z","iopub.execute_input":"2026-01-09T05:21:27.14813Z","iopub.status.idle":"2026-01-09T05:21:29.279124Z","shell.execute_reply.started":"2026-01-09T05:21:27.148099Z","shell.execute_reply":"2026-01-09T05:21:29.278086Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## RNA Structure Visualization ","metadata":{}},{"cell_type":"code","source":"import py3Dmol\nimport ipywidgets\n\ndef view_cif_kaggle(cif_path):\n    with open(cif_path, 'r') as f:\n        cif_data = f.read()\n    \n    view = py3Dmol.view(width=800, height=600)\n    view.addModel(cif_data, 'cif')\n    view.setStyle({}, {'cartoon': {'color': 'spectrum'}})\n    view.zoomTo()\n    return view\n\nview = view_cif_kaggle('/kaggle/input/stanford-rna-3d-folding-2/PDB_RNA/124d.cif') #3rd file in PDB_RNA/\nview","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-09T05:21:29.280263Z","iopub.execute_input":"2026-01-09T05:21:29.280562Z","iopub.status.idle":"2026-01-09T05:21:29.304822Z","shell.execute_reply.started":"2026-01-09T05:21:29.280528Z","shell.execute_reply":"2026-01-09T05:21:29.304142Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## MSA Conservation Score \n\n**Multiple Sequence Alignment (MSA)** aligns your target RNA sequence with evolutionarily related sequences. \n\n**Conservation scores** show how often each position has the **same nucleotide** across all sequences\n\nFor the given file below these are the stats:\n- Position 16: ~0.85 conservation → Highly stable structural element\n- Position 17: ~0.80 conservation → Base-paired region\n- Position 13: ~0.75 conservation → Helix core\n","metadata":{}},{"cell_type":"code","source":"from Bio import AlignIO\n\ndef plot_msa_conservation(fasta_path):\n    \"\"\"Conservation plot for MSA\"\"\"\n    alignment = AlignIO.read(fasta_path, \"fasta\")\n    length = alignment.get_alignment_length()\n    \n    # Calculate conservation\n    conservation = []\n    for i in range(length):\n        col = [str(rec.seq[i]) for rec in alignment]\n        cons = 1 - len(set(col))/4  # Simple conservation score\n        conservation.append(cons)\n    \n    plt.figure(figsize=(15, 4))\n    plt.plot(conservation, linewidth=2)\n    plt.fill_between(range(length), conservation, alpha=0.3)\n    plt.title('MSA Conservation Score by Position')\n    plt.xlabel('Alignment Position')\n    plt.ylabel('Conservation (0-1)')\n    plt.grid(True, alpha=0.3)\n    plt.show()\n    \n    print(f\"MSA: {len(alignment)} sequences, {length} positions\")\n    print(\"Most conserved regions:\", np.argsort(conservation)[-5:][::-1])\n\nplot_msa_conservation('/kaggle/input/stanford-rna-3d-folding-2/MSA/17RA.MSA.fasta') #3rd file in MSA/\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-09T05:21:29.305759Z","iopub.execute_input":"2026-01-09T05:21:29.305992Z","iopub.status.idle":"2026-01-09T05:21:29.575158Z","shell.execute_reply.started":"2026-01-09T05:21:29.305971Z","shell.execute_reply":"2026-01-09T05:21:29.574157Z"}},"outputs":[],"execution_count":null}]}