{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.10.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":87793,"databundleVersionId":11403143,"sourceType":"competition"},{"sourceId":10855324,"sourceType":"datasetVersion","datasetId":6742586}],"isInternetEnabled":false,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nfrom sklearn import *\nfrom collections import Counter\nfrom IPython.display import display, HTML\nimport xgboost as xgb\nfrom xgboost import XGBRegressor\nimport warnings, pickle, tqdm\nwarnings.filterwarnings(\"ignore\")\n#import hdbscan\n\np = '/kaggle/input/stanford-rna-3d-folding/'\n\ntrains = pd.read_csv(p+'train_sequences.csv')\ntrainl = pd.read_csv(p+'train_labels.csv')\n\nvals = pd.read_csv(p+'validation_sequences.csv')\nvall = pd.read_csv(p+'validation_labels.csv')\n\ntests = pd.read_csv(p+'test_sequences.csv')\nsubs = pd.read_csv(p+'sample_submission.csv')","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true,"execution":{"iopub.status.busy":"2025-03-16T18:33:38.681984Z","iopub.execute_input":"2025-03-16T18:33:38.682456Z","iopub.status.idle":"2025-03-16T18:33:41.935115Z","shell.execute_reply.started":"2025-03-16T18:33:38.682416Z","shell.execute_reply":"2025-03-16T18:33:41.934062Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Temporal Cutoff & Null Review","metadata":{}},{"cell_type":"code","source":"print('Train Sequence Before: ', len(trains))\ntrains = trains[trains['temporal_cutoff']<'2022-05-27']\nprint('Train Sequence After: ', len(trains))\n\n#Remove Nulls in XYZ, whole sequence\nremove_ids = trainl[trainl.isna().any(axis=1)]['ID'].values\nremove_ids = list(set(['_'.join(i.split('_')[:-1]) for i in remove_ids]))\nprint('Removed IDs: ', len(remove_ids))\n\ntrains = trains[~(trains['target_id'].isin(remove_ids))]\nprint('Train Sequence After: ', len(trains))\n\ntrains = trains.reset_index(drop=True)\n\ncutoff_train_ids = []\nfor i in range(len(trains)):\n    target = trains['target_id'][i]\n    seq = [s for s in trains['sequence'][i]]\n    for j, s in enumerate(seq):\n        res = str(target) + '_' + str(j+1)\n        cutoff_train_ids.append(res)\n\nprint('Train Labels Before: ', len(trainl))\ntrainl = trainl[trainl['ID'].isin(cutoff_train_ids)]\nprint('Train Labels After: ', len(trainl))","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-03-16T18:33:41.936162Z","iopub.execute_input":"2025-03-16T18:33:41.936451Z","iopub.status.idle":"2025-03-16T18:33:42.103834Z","shell.execute_reply.started":"2025-03-16T18:33:41.936417Z","shell.execute_reply":"2025-03-16T18:33:42.102729Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Exploring possibilities in Sequences","metadata":{}},{"cell_type":"code","source":"RNA_SEQ_TRANS = {}\nfrom itertools import product\nli = ['A', 'C', 'G', 'U']\n\nfor i in range(1,12):\n    comb = product(li, repeat=i)\n    RNA_SEQ_TRANS[i] = {''.join(l):i for i,l in enumerate(comb)}\n\nfor i in RNA_SEQ_TRANS:  #4**i\n    print(i, len(RNA_SEQ_TRANS[i]))","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-03-16T18:33:42.105015Z","iopub.execute_input":"2025-03-16T18:33:42.105373Z","iopub.status.idle":"2025-03-16T18:33:46.044216Z","shell.execute_reply.started":"2025-03-16T18:33:42.105341Z","shell.execute_reply":"2025-03-16T18:33:46.043062Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Metric\nhttps://www.kaggle.com/code/metric/ribonanza-tm-score","metadata":{}},{"cell_type":"code","source":"#####| default_exp core #HAVE TO RESOLVE PATH ISSUES","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-03-16T11:23:10.889593Z","iopub.execute_input":"2025-03-16T11:23:10.889865Z","iopub.status.idle":"2025-03-16T11:23:10.893817Z","shell.execute_reply.started":"2025-03-16T11:23:10.889840Z","shell.execute_reply":"2025-03-16T11:23:10.892898Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"#####| export #HAVE TO RESOLVE PATH ISSUES\n\nimport pandas.api.types\nimport pandas as pd\nimport os, re, shutil\nimport kagglehub\n\ntemp_path = '/kaggle/tempfile/'\nif not os.path.exists(temp_path):\n    os.makedirs(temp_path)\nusalign_path = kagglehub.dataset_download('metric/usalign')\nif not os.path.exists(temp_path + '/USalign'):\n    shutil.copy(usalign_path + '/USalign', temp_path + 'USalign')\n\n!chmod +rwx {temp_path}/USalign\n\nclass RNA3DFMetric:\n    def __init__(self):\n        pass\n    \n    def parse_tmscore_output(self, output):\n        tm_score_match = re.findall(r'TM-score=\\s+([\\d.]+)', output)[1]\n        if not tm_score_match:\n            raise ValueError('No TM score found')\n        return float(tm_score_match)\n\n    def write_target_line(self, asr, an, rn, rnum, xc, yc, zc, atom_type='P') -> str:\n        return f'ATOM{asr:>7d}  {an:<6s}{rn:<4s}{rnum:>3d}{xc:>12.3f}{yc:>8.3f}{zc:>8.3f}  1.00  0.00{atom_type:>12s}\\n'\n    \n    def write2pdb(self, df: pd.DataFrame, xyz_id: str, target_path: str) -> int:\n        resolved_cnt = 0\n        with open(target_path, 'w') as target_file:\n            for _, row in df.iterrows():\n                x_coord = row[f'x_{xyz_id}']\n                y_coord = row[f'y_{xyz_id}']\n                z_coord = row[f'z_{xyz_id}']\n                if x_coord > -1e17 and y_coord > -1e17 and z_coord > -1e17:\n                    resolved_cnt += 1\n                    target_line = self.write_target_line(int(row['resid']), \"C1'\", row['resname'], int(row['resid']),x_coord,y_coord,z_coord,'C')\n                    target_file.write(target_line)\n        return resolved_cnt\n\n    def score(self, solution: pd.DataFrame, submission: pd.DataFrame, verbose=False) -> float:\n        solution['target_id'] = solution['ID'].apply(lambda x: x.split('_')[0])\n        submission['target_id'] = submission['ID'].apply(lambda x: x.split('_')[0])\n        results = []\n        for target_id, group_native in solution.groupby('target_id'):\n            group_predicted = submission[submission['target_id'] == target_id]\n            native_pdb = 'native.pdb'\n            predicted_pdb = 'predicted.pdb'\n            target_id_scores = []\n            for pred_cnt in range(1, 6):\n                prediction_scores = []\n                for native_cnt in range(1, 41):\n                    resolved_cnt = self.write2pdb(group_native, native_cnt, native_pdb)\n                    _ = self.write2pdb(group_predicted, pred_cnt, predicted_pdb)\n                    if resolved_cnt > 0:\n                        command = f'{temp_path}/USalign {predicted_pdb} {native_pdb} -atom \" C1\\'\"'\n                        usalign_output = os.popen(command).read()\n                        if verbose: print(usalign_output)\n                        prediction_scores.append(self.parse_tmscore_output(usalign_output))\n                        print(prediction_scores[-1])\n                target_id_scores.append(max(prediction_scores))\n            results.append(max(target_id_scores))\n        return float(sum(results) / len(results))","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-03-16T18:33:46.045403Z","iopub.execute_input":"2025-03-16T18:33:46.045798Z","iopub.status.idle":"2025-03-16T18:33:46.668437Z","shell.execute_reply.started":"2025-03-16T18:33:46.045761Z","shell.execute_reply":"2025-03-16T18:33:46.666805Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Develop Training Set Features","metadata":{}},{"cell_type":"code","source":"bases = {'A':{'id': 0, 'x':0, 'y':0, 'z':0},\n         'C':{'id': 1, 'x':0, 'y':0, 'z':0},\n         'G':{'id': 2, 'x':0, 'y':0, 'z':0}, \n         'U':{'id': 3, 'x':0, 'y':0, 'z':0}}\n\nfor b in ['A','C','G','U']:\n    bases[b]['x'] =  trainl[((trainl['resid']==1) & (trainl['resname']==b))]['x_1'].mean()\n    bases[b]['y'] =  trainl[((trainl['resid']==1) & (trainl['resname']==b))]['y_1'].mean()\n    bases[b]['z'] =  trainl[((trainl['resid']==1) & (trainl['resname']==b))]['z_1'].mean()\nbases","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-03-16T18:33:46.672880Z","iopub.execute_input":"2025-03-16T18:33:46.673282Z","iopub.status.idle":"2025-03-16T18:33:46.785666Z","shell.execute_reply.started":"2025-03-16T18:33:46.673236Z","shell.execute_reply":"2025-03-16T18:33:46.784609Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"#ADDITIONAL FEATURE IDEAS:\n#  Normalize starting points for training and validation data (Min/Max scale to 0 by IDs)\n#  Test averaging the previous points from the different model outputs\n#  Try blending outputs with the NN models\n#  Add +/- Sequence Bases\n\ndef ngrams(seq, n):\n    ngrams = zip(*[seq[i:] for i in range(n)])\n    return {''.join(k):v for k,v in Counter(ngrams).items()}\n    \ndef getBasePositionFeatures(df, dflabels, labels=True, models=[]):\n    if labels:\n        dflabels = {l: [x,y,z] for l, x, y, z in trainl[['ID', 'x_1', 'y_1', 'z_1']].values}\n    df['seq_len'] = df['sequence'].map(len)\n    ngramc = ['AA', 'AC', 'AG', 'AU', 'CA', 'CC', 'CG', 'CU', 'GA', 'GC', 'GG', 'GU', 'UA', 'UC', 'UG', 'UU', 'AAA', 'AAC', 'AAG', 'AAU', 'ACA', 'ACC', 'ACG', 'ACU', 'AGA', 'AGC', 'AGG', 'AGU', 'AUA', 'AUC', 'AUG', 'AUU', 'CAA', 'CAC', 'CAG', 'CAU', 'CCA', 'CCC', 'CCG', 'CCU', 'CGA', 'CGC', 'CGG', 'CGU', 'CUA', 'CUC', 'CUG', 'CUU', 'GAA', 'GAC', 'GAG', 'GAU', 'GCA', 'GCC', 'GCG', 'GCU', 'GGA', 'GGC', 'GGG', 'GGU', 'GUA', 'GUC', 'GUG', 'GUU', 'UAA', 'UAC', 'UAG', 'UAU', 'UCA', 'UCC', 'UCG', 'UCU', 'UGA', 'UGC', 'UGG', 'UGU', 'UUA', 'UUC', 'UUG', 'UUU']\n    cols = ['ID', 'resname', 'base', 'resid', 'seq_len', 'prior_a_count', 'prior_c_count', 'prior_g_count', 'prior_u_count', 'A', 'C', 'G', 'U'] + ngramc + ['prev_x', 'prev_y', 'prev_z', 'x_1', 'y_1', 'z_1']\n    tcols = [c for c in cols if c not in ['ID', 'resname', 'x_1', 'y_1', 'z_1']]\n\n    sub = []\n    for i in range(len(df)):\n        lstart = [0] * len(cols)\n        lstart[cols.index('ID')] = str(df['target_id'][i])\n        lstart[cols.index('seq_len')] = df['seq_len'][i]\n        total_bases = dict(Counter(df['sequence'][i]))\n        for k in total_bases:\n            lstart[cols.index(k)] = total_bases[k]\n        ngrams2 = ngrams(df['sequence'][i], 2)\n        for k in ngrams2:\n            lstart[cols.index(k)] = ngrams2[k]\n        ngrams3 = ngrams(df['sequence'][i], 3)\n        for k in ngrams3:\n            lstart[cols.index(k)] = ngrams3[k]\n        seq = [s for s in df['sequence'][i]]\n        for j, s in enumerate(seq):\n            l_item = lstart[:]\n            l_item[cols.index('ID')] += '_' + str(j+1)\n            ID_ = l_item[cols.index('ID')]\n            l_item[cols.index('resname')] = s\n            l_item[cols.index('base')] = bases[s]['id']\n            l_item[cols.index('resid')] = j+1\n            \n            c = dict(Counter(df['sequence'][i][:j+1]))\n            for k in c:\n                l_item[cols.index('prior_' + str(k).lower() + '_count')] = c[k]\n\n            #Last Labels\n            if j+1 > 1:\n                l_item[cols.index('prev_x')] = sub[-1][cols.index('x_1')]\n                l_item[cols.index('prev_y')] = sub[-1][cols.index('y_1')]\n                l_item[cols.index('prev_z')] = sub[-1][cols.index('z_1')]\n                \n            if labels: #Training\n                l_item[cols.index('x_1')] = dflabels[ID_][0]\n                l_item[cols.index('y_1')] = dflabels[ID_][1]\n                l_item[cols.index('z_1')] = dflabels[ID_][2]\n            else: #Prediction\n                l_item = l_item[:-3] #remove x_1, y_1, z_1\n                resp = [l_item[cols.index(k)] for k in tcols]\n                for mid, model in enumerate(models):\n                    dfp = pd.DataFrame([resp], columns=tcols)\n                    if j+1 > 1:\n                        dfp['prev_x'] = sub[-1][cols.index('x_1') + (mid*3)]\n                        dfp['prev_y'] = sub[-1][cols.index('y_1') + (mid*3)]\n                        dfp['prev_z'] = sub[-1][cols.index('z_1') + (mid*3)]\n                    l_item +=  list(model.predict(dfp)[0])\n            sub.append(l_item)\n    if labels == False:\n        cols += ['x_2', 'y_2', 'z_2', 'x_3', 'y_3', 'z_3', 'x_4', 'y_4', 'z_4', 'x_5', 'y_5', 'z_5']\n    sub = pd.DataFrame(sub, columns=cols)\n    return sub","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-03-16T11:23:11.430506Z","iopub.execute_input":"2025-03-16T11:23:11.430782Z","iopub.status.idle":"2025-03-16T11:23:11.449146Z","shell.execute_reply.started":"2025-03-16T11:23:11.430757Z","shell.execute_reply":"2025-03-16T11:23:11.448017Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"%%time\ntrain = getBasePositionFeatures(trains, trainl)  #Extend to 5 sets for use with 5 models\ntrain.head()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-03-16T11:23:11.450189Z","iopub.execute_input":"2025-03-16T11:23:11.450429Z","iopub.status.idle":"2025-03-16T11:23:20.079054Z","shell.execute_reply.started":"2025-03-16T11:23:11.450406Z","shell.execute_reply":"2025-03-16T11:23:20.077982Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"%%time\n\ntcols = ['base', 'resid', 'seq_len', 'prior_a_count', 'prior_c_count', 'prior_g_count', 'prior_u_count', 'A', 'C', 'G', 'U', 'AA', 'AC', 'AG', 'AU', 'CA', 'CC', 'CG', 'CU', 'GA', 'GC', 'GG', 'GU', 'UA', 'UC', 'UG', 'UU', 'AAA', 'AAC', 'AAG', 'AAU', 'ACA', 'ACC', 'ACG', 'ACU', 'AGA', 'AGC', 'AGG', 'AGU', 'AUA', 'AUC', 'AUG', 'AUU', 'CAA', 'CAC', 'CAG', 'CAU', 'CCA', 'CCC', 'CCG', 'CCU', 'CGA', 'CGC', 'CGG', 'CGU', 'CUA', 'CUC', 'CUG', 'CUU', 'GAA', 'GAC', 'GAG', 'GAU', 'GCA', 'GCC', 'GCG', 'GCU', 'GGA', 'GGC', 'GGG', 'GGU', 'GUA', 'GUC', 'GUG', 'GUU', 'UAA', 'UAC', 'UAG', 'UAU', 'UCA', 'UCC', 'UCG', 'UCU', 'UGA', 'UGC', 'UGG', 'UGU', 'UUA', 'UUC', 'UUG', 'UUU', 'prev_x', 'prev_y', 'prev_z']\npcols = ['x_1', 'y_1', 'z_1']\n\nmodels = []\nfor i in tqdm.tqdm(range(5)): \n    model = XGBRegressor(n_estimators=7000, max_depth=6+i, learning_rate=0.2, tree_method='hist', n_jobs=-1, random_state=27)\n    model.fit(train[tcols], train[pcols])\n    models.append(model)\n    with open('model'+str(i)+'.pkl', \"wb\") as f:\n        pickle.dump(model, f)\n#models = [models[0]] * 5","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-03-16T11:36:16.876345Z","iopub.execute_input":"2025-03-16T11:36:16.876726Z","iopub.status.idle":"2025-03-16T11:50:49.732021Z","shell.execute_reply.started":"2025-03-16T11:36:16.876698Z","shell.execute_reply":"2025-03-16T11:50:49.730529Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"#models = []\n#for i in range(5):\n    #with open('model'+str(i)+'.pkl', \"rb\") as f:\n    #    model = pickle.load(f)\n    #medels.append(model)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-03-16T11:23:39.059229Z","iopub.execute_input":"2025-03-16T11:23:39.059493Z","iopub.status.idle":"2025-03-16T11:23:39.064927Z","shell.execute_reply.started":"2025-03-16T11:23:39.059466Z","shell.execute_reply":"2025-03-16T11:23:39.063143Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"%%time\nval = getBasePositionFeatures(vals, vall, False, models)\nval.head()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-03-16T11:58:18.840718Z","iopub.execute_input":"2025-03-16T11:58:18.841162Z","iopub.status.idle":"2025-03-16T12:00:18.175394Z","shell.execute_reply.started":"2025-03-16T11:58:18.841131Z","shell.execute_reply":"2025-03-16T12:00:18.173766Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"scols = subs.columns\nsub = val[:]\nsub = sub[scols]","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-03-16T12:00:18.176139Z","iopub.execute_input":"2025-03-16T12:00:18.176396Z","iopub.status.idle":"2025-03-16T12:00:18.181425Z","shell.execute_reply.started":"2025-03-16T12:00:18.176371Z","shell.execute_reply":"2025-03-16T12:00:18.180317Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"%%time\nRNA3DFM = RNA3DFMetric()\nRNA3DFM.score(vall, sub, False)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-03-16T12:00:18.182421Z","iopub.execute_input":"2025-03-16T12:00:18.182681Z","iopub.status.idle":"2025-03-16T12:00:30.258771Z","shell.execute_reply.started":"2025-03-16T12:00:18.182656Z","shell.execute_reply":"2025-03-16T12:00:30.257395Z"},"_kg_hide-output":false},"outputs":[],"execution_count":null},{"cell_type":"code","source":"%%time\ntest = getBasePositionFeatures(tests, tests, False, models)\ntest.head()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-03-16T11:32:58.003859Z","iopub.status.idle":"2025-03-16T11:32:58.004644Z","shell.execute_reply.started":"2025-03-16T11:32:58.004031Z","shell.execute_reply":"2025-03-16T11:32:58.004077Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"scols = subs.columns\nsub = test[:]\nsub = sub[scols]\nsub.to_csv('submission.csv', index=False)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-03-16T11:32:58.005428Z","iopub.status.idle":"2025-03-16T11:32:58.006030Z","shell.execute_reply.started":"2025-03-16T11:32:58.005619Z","shell.execute_reply":"2025-03-16T11:32:58.005659Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"#https://zhanggroup.org/US-align/help/\n#Reference Visualizations","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-03-16T11:31:09.164310Z","iopub.status.idle":"2025-03-16T11:31:09.164749Z","shell.execute_reply.started":"2025-03-16T11:31:09.164467Z","shell.execute_reply":"2025-03-16T11:31:09.164506Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Ｈ𝐀𝑷𝑷𝓎 🇰𝗮𝘨𝘨🇱𝖎Ｎɢ 💯","metadata":{}}]}