{"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":31259,"isInternetEnabled":false,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"#!/usr/bin/env python3\n\"\"\"\nStanford RNA 3D Folding Challenge — Phase Transition-Aware RNA Predictor\nOptimized Implementation (SciSpace Edition)\n\nConcept:\n--------\nThis model treats RNA folding as a phase-transition system. Sequences near\ncritical points (balanced GC, high entropy, symmetric purine-pyrimidine content)\nexhibit structural bistability. We adapt ensemble sampling accordingly.\n\nKey improvements:\n1. Integrates base-pair probability map (pseudo-thermodynamic constraint)\n2. Uses energy-weighted ensemble sampling instead of pure random noise\n3. Dynamically tunes temperature based on criticality\n4. Smooths coordinates with local backbone regularization\n\nAuthor: SciSpace Research Assistant\n\"\"\"\n\n# ============================================================\n# Setup\n# ============================================================\nimport os\nimport numpy as np\nimport pandas as pd\nfrom pathlib import Path\nfrom collections import Counter\n\nnp.random.seed(42)\n\n# Canonical base map\nNUCLEOTIDE_MAPPING = {'A': 'A', 'U': 'U', 'G': 'G', 'C': 'C', 'T': 'U', 'I': 'A'}\n\n# ============================================================\n# Feature Computation\n# ============================================================\n\ndef clean_sequence(seq: str) -> str:\n    \"\"\"Convert modified nucleotides to canonical A/U/G/C.\"\"\"\n    return \"\".join(NUCLEOTIDE_MAPPING.get(b, 'A') for b in seq)\n\ndef compute_entropy(seq: str) -> float:\n    \"\"\"Compute normalized Shannon entropy of nucleotide composition.\"\"\"\n    counts = Counter(seq)\n    total = len(seq)\n    entropy = -sum((c/total) * np.log2(c/total) for c in counts.values() if c > 0)\n    return entropy / 2.0  # max entropy for A,C,G,U = log2(4)=2\n\ndef compute_gc_criticality(seq: str) -> float:\n    \"\"\"Measure proximity to 50% GC (optimal folding flexibility).\"\"\"\n    gc_frac = (seq.count('G') + seq.count('C')) / len(seq)\n    return 1 - 2 * abs(gc_frac - 0.5)\n\ndef compute_purine_pyrimidine_balance(seq: str) -> float:\n    \"\"\"Purine (A,G) vs Pyrimidine (C,U) balance.\"\"\"\n    pur = seq.count('A') + seq.count('G')\n    pyr = seq.count('C') + seq.count('U')\n    return 1 - abs(pur - pyr) / len(seq)\n\ndef compute_criticality_score(seq: str) -> float:\n    \"\"\"Weighted phase-transition criticality.\"\"\"\n    seq = clean_sequence(seq)\n    gc_c = compute_gc_criticality(seq)\n    ent = compute_entropy(seq)\n    bal = compute_purine_pyrimidine_balance(seq)\n    run_penalty = 1.0 - min(max_run_length(seq) / len(seq), 1.0)\n    return np.clip(0.4*gc_c + 0.3*ent + 0.2*bal + 0.1*run_penalty, 0, 1)\n\ndef max_run_length(seq: str) -> int:\n    \"\"\"Longest homopolymer run.\"\"\"\n    max_run = 1\n    for base in 'ACGU':\n        run, longest = 0, 1\n        for b in seq:\n            if b == base:\n                run += 1\n                longest = max(longest, run)\n            else:\n                run = 0\n        max_run = max(max_run, longest)\n    return max_run\n\n# ============================================================\n# Pseudo Base-Pair Probability Model (Thermodynamic Prior)\n# ============================================================\n\ndef base_pair_probabilities(seq: str) -> np.ndarray:\n    \"\"\"\n    Compute a heuristic base-pair probability matrix (NxN).\n    Uses simple complementarity weighting + distance penalty.\n    \"\"\"\n    n = len(seq)\n    bp_matrix = np.zeros((n, n))\n    pairs = {('A', 'U'), ('U', 'A'), ('G', 'C'), ('C', 'G'), ('G', 'U'), ('U', 'G')}\n    \n    for i in range(n):\n        for j in range(i+3, n):  # minimum loop length constraint\n            if (seq[i], seq[j]) in pairs:\n                # distance-penalized complementarity\n                bp_matrix[i, j] = 1.0 / (1.0 + 0.05*(j - i))\n    \n    # Normalize row-wise probabilities\n    bp_matrix = bp_matrix / (bp_matrix.sum(axis=1, keepdims=True) + 1e-6)\n    return bp_matrix\n\n# ============================================================\n# Structure Generator\n# ============================================================\n\ndef generate_structure(seq: str, temperature: float, bp_map: np.ndarray) -> np.ndarray:\n    \"\"\"\n    Generate 3D coordinates guided by base-pair probabilities and temperature.\n    \"\"\"\n    n = len(seq)\n    coords = np.zeros((n, 3))\n    angle_step = 32.7 * np.pi / 180\n    radius, rise = 10.0, 2.8\n\n    for i in range(n):\n        angle = i * angle_step\n        x, y, z = radius * np.cos(angle), radius * np.sin(angle), i * rise\n        \n        # Attraction toward base-paired partners (pseudo-hydrogen bonds)\n        partner_force = np.zeros(3)\n        for j in range(n):\n            if bp_map[i, j] > 0:\n                ideal_dist = 10.0\n                partner_vec = np.array([\n                    ideal_dist * np.cos(j * angle_step),\n                    ideal_dist * np.sin(j * angle_step),\n                    j * rise\n                ]) - np.array([x, y, z])\n                partner_force += bp_map[i, j] * partner_vec\n        \n        # Apply energy-scaled noise (higher T = more disorder)\n        noise = np.random.normal(0, temperature * 1.2, size=3)\n        coords[i] = np.array([x, y, z]) + partner_force * 0.05 + noise\n\n    # Backbone smoothing\n    for i in range(1, n-1):\n        coords[i] = 0.25*coords[i-1] + 0.5*coords[i] + 0.25*coords[i+1]\n    \n    return coords\n\n# ============================================================\n# Ensemble Generator\n# ============================================================\n\ndef generate_adaptive_ensemble(seq: str, crit: float, n=5):\n    bp_map = base_pair_probabilities(seq)\n    \n    # Temperature scaling by criticality\n    if crit < 0.3:\n        temps = np.linspace(0.4, 0.8, n)\n    elif crit > 0.7:\n        temps = np.linspace(1.0, 3.0, n)\n    else:\n        temps = np.linspace(0.6, 2.0, n)\n    \n    return [generate_structure(seq, t, bp_map) for t in temps]\n\n# ============================================================\n# Main Pipeline\n# ============================================================\n\ndef predict_rna_structures(df_test, df_template=None):\n    all_preds = []\n    print(f\"\\nProcessing {len(df_test)} sequences...\\n\")\n    \n    for idx, row in df_test.iterrows():\n        seq_id = row['ID']\n        seq = clean_sequence(row['sequence'])\n        crit = compute_criticality_score(seq)\n        structures = generate_adaptive_ensemble(seq, crit, n=5)\n        \n        for i, base in enumerate(seq):\n            row_dict = {'ID': f\"{seq_id}_{i+1}\", 'resname': base, 'resid': i+1}\n            for k in range(5):\n                row_dict[f\"x_{k+1}\"], row_dict[f\"y_{k+1}\"], row_dict[f\"z_{k+1}\"] = structures[k][i]\n            all_preds.append(row_dict)\n        \n        print(f\"[{idx+1}] {seq_id} | Criticality={crit:.3f}\")\n    \n    pred_df = pd.DataFrame(all_preds)\n    if df_template is not None:\n        pred_df = pred_df.reindex(columns=df_template.columns, fill_value=0.0)\n    \n    return pred_df\n\n# ============================================================\n# Execution\n# ============================================================\n\ndef run_pipeline():\n    test_path = Path(\"test_sequences.csv\")\n    sub_path = Path(\"sample_submission.csv\")\n    \n    if not test_path.exists():\n        print(\"⚠️ No test set found. Using dummy dataset.\")\n        df_test = pd.DataFrame({\n            'ID': ['RNA1', 'RNA2'],\n            'sequence': ['ACGUACGUACGU', 'GGCCAUUAGCGA']\n        })\n    else:\n        df_test = pd.read_csv(test_path)\n    \n    df_template = pd.read_csv(sub_path) if sub_path.exists() else None\n    df_pred = predict_rna_structures(df_test, df_template)\n    \n    out_file = \"submission_optimized.csv\"\n    df_pred.to_csv(out_file, index=False)\n    print(f\"\\n✅ Saved optimized predictions to {out_file}\")\n\nif __name__ == \"__main__\":\n    run_pipeline()","metadata":{"trusted":true},"outputs":[],"execution_count":null}]}