{"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.12.12"},"kaggle":{"accelerator":"none","dataSources":[{"sourceId":118765,"databundleVersionId":15231210,"isSourceIdPinned":false,"sourceType":"competition"},{"sourceId":11775065,"sourceType":"datasetVersion","datasetId":7392749},{"sourceId":14441699,"sourceType":"datasetVersion","datasetId":9224635}],"dockerImageVersionId":31234,"isInternetEnabled":false,"language":"python","sourceType":"notebook","isGpuEnabled":false},"papermill":{"default_parameters":{},"duration":217.135968,"end_time":"2026-01-09T06:43:37.639371","environment_variables":{},"exception":null,"input_path":"__notebook__.ipynb","output_path":"__notebook__.ipynb","parameters":{},"start_time":"2026-01-09T06:40:00.503403","version":"2.6.0"}},"nbformat_minor":4,"nbformat":4,"cells":[{"id":"b66340c4","cell_type":"markdown","source":"# 🏛️ Stanford RNA 3D Folding: TBM Pipeline (MMseqs2) 🧬\n\nThis notebook implements **Template-Based Modeling (TBM)** using **MMseqs2**.\nIt searches for similar RNA structures in the PDB database and maps their coordinates to the target sequences.\nThis is often the strongest component of a structure prediction pipeline.\n\n### Workflow\n1.  **Setup**: Install MMseqs2 and Biopython (Offline Mode).\n2.  **Database**: Create a searchable DB from `PDB_RNA` sequences.\n3.  **Search**: Align test sequences against the DB.\n4.  **Extract**: Parse `.cif` files for matching templates and extract 3D coordinates.\n5.  **Submit**: Generate `submission.csv` using the best templates found.","metadata":{"papermill":{"duration":0.003645,"end_time":"2026-01-09T06:40:04.312008","exception":false,"start_time":"2026-01-09T06:40:04.308363","status":"completed"},"tags":[]}},{"id":"175eeea1","cell_type":"code","source":"import os\nimport sys\nimport glob\nimport csv\nimport gzip\nimport shutil\nimport warnings\nimport pandas as pd\nimport numpy as np\nfrom datetime import datetime\n\n# 📦 Offline Installation (Hardcoded Paths)\nNUMPY_WHL = '/kaggle/input/biopython-cp312/numpy-2.2.6-cp312-cp312-manylinux_2_17_x86_64.manylinux2014_x86_64.whl'\nBIO_WHL = '/kaggle/input/biopython-cp312/biopython-1.86-cp312-cp312-manylinux_2_17_x86_64.manylinux2014_x86_64.whl'\n\n# Check & Install Numpy\ntry:\n    import numpy\n    print(\"✅ Numpy already installed.\")\nexcept ImportError:\n    if os.path.exists(NUMPY_WHL):\n        print(f\"Installing Numpy from {NUMPY_WHL}...\")\n        !pip install {NUMPY_WHL}\n    else:\n        print(f\"❌ Numpy wheel not found at {NUMPY_WHL}\")\n\n# Check & Install Biopython\ntry:\n    import Bio\n    print(\"✅ Biopython already installed.\")\nexcept ImportError:\n    if os.path.exists(BIO_WHL):\n        print(f\"Installing Biopython from {BIO_WHL}...\")\n        !pip install {BIO_WHL}\n    else:\n        sys.exit(f\"❌ Biopython wheel not found at {BIO_WHL}\")\n\n# Import Biopython Modules\ntry:\n    from Bio import SeqIO, PDB, Align, BiopythonWarning\n    from Bio.PDB.MMCIF2Dict import MMCIF2Dict\n    from Bio.PDB import MMCIFParser\n    print(\"✅ Biopython import successful (including Align).\")\nexcept ImportError as e:\n    sys.exit(f\"❌ Import failed: {e}\")\n\n# Suppress Warnings (Safe now because BiopythonWarning is imported)\nwarnings.simplefilter('ignore', BiopythonWarning)\n\n# ⚙️ Config\nclass Config:\n    MAX_TEMPLATES = 5\n    NULL_VALUE = 0.0\n    CHECK_CUTOFF = True\n    BASE_DIR = '/kaggle/input/stanford-rna-3d-folding-2'\n    PDB_DIR = f'{BASE_DIR}/PDB_RNA'\n    WORK_DIR = '/kaggle/working'\n\nprint(\"Config Loaded.\")","metadata":{"execution":{"iopub.status.busy":"2026-01-09T11:40:23.764808Z","iopub.execute_input":"2026-01-09T11:40:23.765365Z","iopub.status.idle":"2026-01-09T11:40:31.861376Z","shell.execute_reply.started":"2026-01-09T11:40:23.765201Z","shell.execute_reply":"2026-01-09T11:40:31.858615Z"},"papermill":{"duration":40.37759,"end_time":"2026-01-09T06:40:44.692266","exception":false,"start_time":"2026-01-09T06:40:04.314676","status":"completed"},"tags":[],"trusted":true},"outputs":[],"execution_count":null},{"id":"01fc6c49","cell_type":"markdown","source":"## 1. Setup MMseqs2 🛠️\nWe use the pre-uploaded Kaggle Dataset (`/kaggle/input/mmseqs2/mmseqs`).","metadata":{"papermill":{"duration":0.003107,"end_time":"2026-01-09T06:40:44.698699","exception":false,"start_time":"2026-01-09T06:40:44.695592","status":"completed"},"tags":[]}},{"id":"33108fb1","cell_type":"code","source":"# ⚙️ MMseqs2 Setup (Hardcoded Path)\nMMSEQS_SRC = '/kaggle/input/mmseqs2/mmseqs/bin/mmseqs'\nMMSEQS_BIN = f\"{Config.WORK_DIR}/mmseqs/bin/mmseqs\"\n\nif os.path.exists(MMSEQS_SRC):\n    print(f\"✅ Found MMseqs2 at {MMSEQS_SRC}\")\n    !mkdir -p {Config.WORK_DIR}/mmseqs/bin\n    !cp {MMSEQS_SRC} {MMSEQS_BIN}\n    !chmod +x {MMSEQS_BIN}\n    print(\"MMseqs2 setup complete.\")\nelse:\n    sys.exit(f\"❌ MMseqs2 not found at {MMSEQS_SRC}\")\n","metadata":{"execution":{"iopub.status.busy":"2026-01-09T11:40:31.864108Z","iopub.execute_input":"2026-01-09T11:40:31.864922Z","iopub.status.idle":"2026-01-09T11:40:32.39446Z","shell.execute_reply.started":"2026-01-09T11:40:31.864875Z","shell.execute_reply":"2026-01-09T11:40:32.393273Z"},"papermill":{"duration":0.508669,"end_time":"2026-01-09T06:40:45.210468","exception":false,"start_time":"2026-01-09T06:40:44.701799","status":"completed"},"tags":[],"trusted":true},"outputs":[],"execution_count":null},{"id":"3726939e","cell_type":"markdown","source":"## 2. Create Database & Search 🔎\n1. Convert `test_sequences.csv` to FASTA.\n2. Create MMseqs2 DB from `PDB_RNA/pdb_seqres_NA.fasta`.\n3. Run `easy-search`.","metadata":{"papermill":{"duration":0.003097,"end_time":"2026-01-09T06:40:45.217064","exception":false,"start_time":"2026-01-09T06:40:45.213967","status":"completed"},"tags":[]}},{"id":"b7818076","cell_type":"code","source":"# 1. Convert Test CSV to FASTA\ntest_csv = f\"{Config.BASE_DIR}/test_sequences.csv\"\ntest_fasta = f\"{Config.WORK_DIR}/test_sequences.fasta\"\n\nif not os.path.exists(test_fasta):\n    df = pd.read_csv(test_csv)\n    with open(test_fasta, 'w') as f:\n        for idx, row in df.iterrows():\n            f.write(f\">{row['target_id']}\\n{row['sequence']}\\n\")\n    print(f\"Converted {len(df)} sequences to FASTA.\")\n\n# 2. Database & Search\npdb_fasta = f\"{Config.PDB_DIR}/pdb_seqres_NA.fasta\"\nsearch_result = f\"{Config.WORK_DIR}/testResult.txt\"\ntmp_dir = f\"{Config.WORK_DIR}/tmp_mmseqs\"\n\ncmd = f\"{MMSEQS_BIN} easy-search {test_fasta} {pdb_fasta} {search_result} {tmp_dir} \" \\\n      f\"--search-type 3 --format-output query,target,evalue,qstart,qend,tstart,tend,qaln,taln\"\n\nprint(\"Running Search... (This may take a few minutes)\")\n!{cmd}\nprint(\"Search Complete.\")","metadata":{"execution":{"iopub.status.busy":"2026-01-09T11:40:32.396104Z","iopub.execute_input":"2026-01-09T11:40:32.396646Z","iopub.status.idle":"2026-01-09T11:40:50.582039Z","shell.execute_reply.started":"2026-01-09T11:40:32.39659Z","shell.execute_reply":"2026-01-09T11:40:50.581079Z"},"papermill":{"duration":18.250112,"end_time":"2026-01-09T06:41:03.47053","exception":false,"start_time":"2026-01-09T06:40:45.220418","status":"completed"},"tags":[],"trusted":true},"outputs":[],"execution_count":null},{"id":"2ede1724","cell_type":"markdown","source":"## 3. Helpers for Coordinate Extraction 🧩\nFunctions to parse `.cif` files and extract C1' atoms.","metadata":{"papermill":{"duration":0.008504,"end_time":"2026-01-09T06:41:03.48808","exception":false,"start_time":"2026-01-09T06:41:03.479576","status":"completed"},"tags":[]}},{"id":"61d44391","cell_type":"code","source":"def clean_res_name(res_name):\n    if res_name in ['A', 'C', 'G', 'U']:\n        return res_name\n    return 'X' # Modified residue\n\ndef is_before_or_on(cutoff_str, release_str):\n    # Returns True if cutoff <= release (i.e. template is NEWER or EQUAL to cutoff)\n    try:\n        d_cut = pd.to_datetime(cutoff_str)\n        d_rel = pd.to_datetime(release_str)\n        return d_cut <= d_rel\n    except:\n        return False\n\ndef extract_rna_sequence(cif_path, chain_id):\n    if cif_path.endswith('.gz'):\n        open_func = gzip.open\n        mode = 'rt'\n    else:\n        open_func = open\n        mode = 'r'\n        \n    with open_func(cif_path, mode) as f:\n        mmcif_dict = MMCIF2Dict(f)\n    \n    strand_ids = mmcif_dict.get('_pdbx_poly_seq_scheme.pdb_strand_id', [])\n    mon_ids = mmcif_dict.get('_pdbx_poly_seq_scheme.mon_id', [])\n    pdb_seq_nums = mmcif_dict.get('_pdbx_poly_seq_scheme.pdb_seq_num', [])\n    \n    full_seq = []\n    seq_nums = []\n    \n    for s_id, m_id, num in zip(strand_ids, mon_ids, pdb_seq_nums):\n        if s_id == chain_id:\n            full_seq.append(clean_res_name(m_id))\n            seq_nums.append(num)\n            \n    return ''.join(full_seq), seq_nums\n\ndef get_c1prime_labels(cif_path, chain_id, alignment, chain_seq_nums):\n    # alignment[0]: Query (A,C,G,-)\n    # alignment[1]: Template (A,C,G,-)\n    \n    parser = MMCIFParser(QUIET=True)\n    try:\n        if cif_path.endswith('.gz'):\n            with gzip.open(cif_path, 'rt') as f:\n                structure = parser.get_structure('RNA', f)\n        else:\n            structure = parser.get_structure('RNA', cif_path)\n            \n        chain = structure[0][chain_id]\n        residues = {str(r.id[1]): r for r in chain}\n        \n        result = []\n        \n        ref_resid_idx = 0 # 1-based index for Query\n        chain_idx = 0     # 0-based index for Template Sequence List\n        \n        for q_char, t_char in zip(alignment[0], alignment[1]):\n            if t_char != '-': \n                chain_idx += 1\n            \n            if q_char != '-':\n                ref_resid_idx += 1\n                \n                # If template has gap or missing residue, fill with NULL\n                if t_char == '-': \n                    result.append((q_char, ref_resid_idx, Config.NULL_VALUE, Config.NULL_VALUE, Config.NULL_VALUE))\n                else:\n                    try:\n                        # Get PDB residue number from our extracted list\n                        pdb_num = chain_seq_nums[chain_idx-1]\n                        residue = residues[str(pdb_num)]\n                        \n                        # Extract C1'\n                        if \"C1'\" in residue:\n                            atom = residue[\"C1'\"]\n                            c = atom.get_coord()\n                            result.append((q_char, ref_resid_idx, c[0], c[1], c[2]))\n                        else:\n                            result.append((q_char, ref_resid_idx, Config.NULL_VALUE, Config.NULL_VALUE, Config.NULL_VALUE))\n                    except KeyError:\n                        result.append((q_char, ref_resid_idx, Config.NULL_VALUE, Config.NULL_VALUE, Config.NULL_VALUE))\n        \n        return result\n        \n    except Exception as e:\n        return []","metadata":{"execution":{"iopub.status.busy":"2026-01-09T11:40:50.583569Z","iopub.execute_input":"2026-01-09T11:40:50.583906Z","iopub.status.idle":"2026-01-09T11:40:50.599207Z","shell.execute_reply.started":"2026-01-09T11:40:50.583871Z","shell.execute_reply":"2026-01-09T11:40:50.598195Z"},"papermill":{"duration":0.0298,"end_time":"2026-01-09T06:41:03.526434","exception":false,"start_time":"2026-01-09T06:41:03.496634","status":"completed"},"tags":[],"trusted":true},"outputs":[],"execution_count":null},{"id":"05c92b0c","cell_type":"markdown","source":"## 4. Processing Results & Generating Submission 📝\n\nIterate through targets, find best templates, and write to CSV.","metadata":{"papermill":{"duration":0.008228,"end_time":"2026-01-09T06:41:03.543517","exception":false,"start_time":"2026-01-09T06:41:03.535289","status":"completed"},"tags":[]}},{"id":"7c61082e","cell_type":"code","source":"# Load Metadata\ndf_test = pd.read_csv(test_csv)\ntemporal_map = dict(zip(df_test['target_id'], df_test['temporal_cutoff']))\nsequence_map = dict(zip(df_test['target_id'], df_test['sequence']))\n\n# Load Release Dates\nrelease_date_file = f\"{Config.PDB_DIR}/pdb_release_dates_NA.csv\"\nrelease_dates = {}\nif os.path.exists(release_date_file):\n    df_dates = pd.read_csv(release_date_file)\n    release_dates = dict(zip(df_dates['Entry ID'], df_dates['Release Date']))\n\n# Parse Search Results\nfrom collections import defaultdict\nhits = defaultdict(list)\nif os.path.exists(search_result):\n    with open(search_result, 'r') as f:\n        for line in f:\n            # query,target,evalue,qstart,qend,tstart,tend,qaln,taln\n            parts = line.strip().split()\n            if len(parts) >= 9:\n                hits[parts[0]].append(parts)\n\n# Process\noutput_rows = []\nprocessed_count = 0\n\nprint(f\"Processing {len(df_test)} targets...\")\n\nfor t_id in df_test['target_id']:\n    processed_count += 1\n    target_seq = sequence_map[t_id]\n    cutoff_date = temporal_map[t_id]\n    L = len(target_seq)\n    \n    templates_found = []\n    \n    # Get Hits for this target\n    target_hits = hits.get(t_id, [])\n    \n    for hit in target_hits:\n        if len(templates_found) >= Config.MAX_TEMPLATES: break\n        \n        _, tmpl_name, evalue, qstart, qend, tstart, tend, qaln, taln = hit\n        pdb_id, chain_id = tmpl_name.split('_')\n        pdb_id = pdb_id.lower()\n        \n        # Check Temporal Cutoff\n        rel_date = release_dates.get(pdb_id.upper(), '2099-01-01')\n        if Config.CHECK_CUTOFF and is_before_or_on(cutoff_date, rel_date):\n            continue # Leakage prevention\n            \n        # Path check\n        cif_path = f\"{Config.PDB_DIR}/{pdb_id}.cif\"\n        if not os.path.exists(cif_path): continue\n        \n        # Extract Sequence Numbers from CIF\n        try:\n            full_chain_seq, chain_seq_nums = extract_rna_sequence(cif_path, chain_id)\n        except Exception as e:\n            continue\n        \n        # Build Alignment Strings with proper padding\n        # Global Alignment Upgrade (Bio.Align.PairwiseAligner)\n        aligner = Align.PairwiseAligner()\n        aligner.mode = 'global'\n        aligner.match_score = 2\n        aligner.mismatch_score = -1\n        aligner.open_gap_score = -2.0\n        aligner.extend_gap_score = -0.5\n        \n        # Align full sequences\n        alignments = aligner.align(target_seq, full_chain_seq)\n        best_aln = alignments[0]\n        query_aln_str = str(best_aln[0])\n        tmpl_aln_str = str(best_aln[1])\n        \n        # Extract Coordinates\n        coords = get_c1prime_labels(cif_path, chain_id, [query_aln_str, tmpl_aln_str], chain_seq_nums)\n        \n        # Validation\n        if len(coords) == L:\n            templates_found.append(coords)\n    \n    # Generate Output Rows\n    for i in range(L):\n        row = {\n            'ID': f\"{t_id}_{i+1}\",\n            'resname': target_seq[i],\n            'resid': i+1\n        }\n        \n        # Fill templates\n        for n, template in enumerate(templates_found):\n            # template[i] is (res, resid, x, y, z)\n            c_res, c_resid, x, y, z = template[i]\n            row[f'x_{n+1}'] = x\n            row[f'y_{n+1}'] = y\n            row[f'z_{n+1}'] = z\n            \n        # Fill remaining with NULL\n        for n in range(len(templates_found), Config.MAX_TEMPLATES):\n             row[f'x_{n+1}'] = Config.NULL_VALUE\n             row[f'y_{n+1}'] = Config.NULL_VALUE\n             row[f'z_{n+1}'] = Config.NULL_VALUE\n             \n        output_rows.append(row)\n        \n    if processed_count % 10 == 0:\n        print(f\"Processed {processed_count}/{len(df_test)}\", end='\\r')\n\nprint(\"\\nSaving submission...\")\ndf_out = pd.DataFrame(output_rows)\ndf_out.to_csv('submission.csv', index=False)\nprint(\"TBM submission.csv saved!\")","metadata":{"execution":{"iopub.status.busy":"2026-01-09T11:40:50.601844Z","iopub.execute_input":"2026-01-09T11:40:50.602162Z"},"papermill":{"duration":153.456473,"end_time":"2026-01-09T06:43:37.008873","exception":false,"start_time":"2026-01-09T06:41:03.5524","status":"completed"},"tags":[],"trusted":true},"outputs":[],"execution_count":null}]}