{"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":11512973,"sourceType":"competition"}],"dockerImageVersionId":30919,"isInternetEnabled":false,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# Overview\n\nSo, our aims to predict the three-dimensional (3D) structure of RNA molecules from their nucleotide sequences.\n\n[Visualisation](https://www.kaggle.com/code/asarvazyan/interactive-3d-sequence-visualization) \n\n[Features ](https://www.kaggle.com/code/dantheshark/rna-3d-folding-understand-data)\n\nhttps://www.kaggle.com/code/fernandosr85/rna-3d-structure/notebook#RNA-3D-Structure-Prediction-with-Reference-Based-Approach-and-Structural-Sampling","metadata":{}},{"cell_type":"code","source":"# packages\n\n# standard\nimport numpy as np\nimport pandas as pd\nimport time\n\nimport nltk\nfrom sklearn.feature_extraction.text import CountVectorizer,TfidfVectorizer\nfrom sklearn.decomposition import TruncatedSVD\n\n# plots\nimport matplotlib.pyplot as plt\nimport plotly.express as px\nimport seaborn as sns\n\nimport plotly.graph_objects as go\n\n# warning handling\nimport warnings\nwarnings.filterwarnings('ignore', category=FutureWarning)  \nwarnings.filterwarnings('ignore', category=RuntimeWarning)\n\n# configs\npd.set_option('display.max_columns', 300)\npd.set_option('display.max_rows', 150)\n\ndefault_color_1 = 'darkblue'\nset_plot = False\n\n\n\nfrom colorama import Style, Fore\nblk = Style.BRIGHT + Fore.BLACK\nred = Style.BRIGHT + Fore.RED\nblu = Style.BRIGHT + Fore.BLUE\nclr = Style.RESET_ALL\n\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-03-25T10:20:33.670488Z","iopub.execute_input":"2025-03-25T10:20:33.670826Z","iopub.status.idle":"2025-03-25T10:20:36.741375Z","shell.execute_reply.started":"2025-03-25T10:20:33.670798Z","shell.execute_reply":"2025-03-25T10:20:36.740675Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Data preprocessing","metadata":{}},{"cell_type":"code","source":"# load data\ndf_train = pd.read_csv('../input/stanford-rna-3d-folding/train_labels.csv')\ndf_valid = pd.read_csv('../input/stanford-rna-3d-folding/validation_labels.csv')\ndf_train_seq = pd.read_csv('../input/stanford-rna-3d-folding/train_sequences.csv')\ndf_valid_seq = pd.read_csv('../input/stanford-rna-3d-folding/validation_sequences.csv')\ndf_test = pd.read_csv('../input/stanford-rna-3d-folding/test_sequences.csv')\n# ------------------\n# Train data\n# ------------------\n\n#df_train_seq[\"all_sequences_str\"] = df_train_seq[\"all_sequences\"].astype(str)\ndf_train_seq['sequence_id'] = df_train_seq[\"target_id\"]\ndf_train_seq[\"sequence_num\"] = df_train_seq[\"all_sequences\"].astype(str).apply(lambda x: \"_\".join(x.split(\"|\")[0:1])[1:])\ndf_train_seq[\"sequence_group\"] = df_train_seq[\"all_sequences\"].astype(str).apply(lambda x: \"_\".join(x.split(\"_\")[0:1])[1:])\ndf_train_seq[\"sequence_group_num\"] = df_train_seq[\"sequence_num\"].apply(lambda x: \"_\".join(x.split(\"_\")[1:2]))\n\n# for nan in all_sequences handle  sequence_group    df_train_seq[df_train_seq[\"all_sequences\"].isna()]\n\ndf_train_seq.loc[df_train_seq.target_id == '2ZJQ_Y', 'sequence_num'] = '2ZJQ_1'\ndf_train_seq.loc[df_train_seq.target_id == '2ZJQ_X', 'sequence_num'] = '2ZJQ_1'\ndf_train_seq.loc[df_train_seq.target_id == '4V65_A1', 'sequence_num'] = '4V65_1'\ndf_train_seq.loc[df_train_seq.target_id == '4V65_BB', 'sequence_num'] = '4V65_1'\ndf_train_seq.loc[df_train_seq.target_id == '4V5F_CA', 'sequence_num'] = '4V5F_1'\n\ndf_train_seq.loc[df_train_seq.target_id == '2ZJQ_Y', 'sequence_group'] = '2ZJQ_Y'\ndf_train_seq.loc[df_train_seq.target_id == '2ZJQ_X', 'sequence_group'] = '2ZJQ_X'\ndf_train_seq.loc[df_train_seq.target_id == '4V65_A1', 'sequence_group'] = '4V65_A1'\ndf_train_seq.loc[df_train_seq.target_id == '4V65_BB', 'sequence_group'] = '4V65_BB'\ndf_train_seq.loc[df_train_seq.target_id == '4V5F_CA', 'sequence_group'] = '4V5F_CA'\n\ndf_train_seq.loc[df_train_seq.target_id == '2ZJQ_Y', 'sequence_group_num'] = 1\ndf_train_seq.loc[df_train_seq.target_id == '2ZJQ_X', 'sequence_group_num'] = 1\ndf_train_seq.loc[df_train_seq.target_id == '4V65_A1', 'sequence_group_num'] = 1\ndf_train_seq.loc[df_train_seq.target_id == '4V65_BB', 'sequence_group_num'] = 1\ndf_train_seq.loc[df_train_seq.target_id == '4V5F_CA', 'sequence_group_num'] = 1\n\ndf_train = df_train.fillna(0)\n\ndf_train[\"sequence_id\"] = df_train[\"ID\"].apply(lambda x: \"_\".join(x.split(\"_\")[:-1]))\ndf_train[\"sequence_target\"] = df_train[\"ID\"].apply(lambda x: \"\".join(x.split(\"_\")[0]))\ndf_train[\"sequence_class\"] = df_train[\"ID\"].apply(lambda x: \"\".join(x.split(\"_\")[1])) #\"_\".join(    #df_train[\"ID\"][1].split(\"_\")[0:2]\n\n\ntrain = pd.merge(df_train, df_train_seq[['sequence','temporal_cutoff',\t'description',\t'all_sequences',\t'sequence_id',\t'sequence_num',\t'sequence_group',\t'sequence_group_num']], on=\"sequence_id\")\n\n# ------------------\n# Valid data\n# ------------------\n# replace extreme values by NaN\ndf_valid.replace(to_replace=-1E18, value=np.nan, inplace=True);\n\ndf_valid_seq['sequence_id'] = df_valid_seq[\"target_id\"]\ndf_valid_seq[\"sequence_num\"] = df_valid_seq[\"all_sequences\"].astype(str).apply(lambda x: \"_\".join(x.split(\"|\")[0:1])[1:])\ndf_valid_seq[\"sequence_group\"] = df_valid_seq[\"all_sequences\"].astype(str).apply(lambda x: \"_\".join(x.split(\"_\")[0:1])[1:])\ndf_valid_seq[\"sequence_group_num\"] = df_valid_seq[\"sequence_num\"].apply(lambda x: \"_\".join(x.split(\"_\")[1:2]))\n\ndf_valid = df_valid.fillna(0)\n\ndf_valid[\"sequence_id\"] = df_valid[\"ID\"].apply(lambda x: \"_\".join(x.split(\"_\")[:-1]))\ndf_valid[\"sequence_target\"] = df_valid[\"ID\"].apply(lambda x: \"\".join(x.split(\"_\")[0]))\ndf_valid[\"sequence_class\"] = df_valid[\"ID\"].apply(lambda x: \"\".join(x.split(\"_\")[1])) #\"_\".join(    #df_train[\"ID\"][1].split(\"_\")[0:2]\n\nvalid = pd.merge(df_valid[['ID',\t'resname',\t'resid',\t'x_1',\t'y_1',\t'z_1', 'sequence_id',\t'sequence_target',\t'sequence_class']], \\\n                 df_valid_seq[['sequence','temporal_cutoff',\t'description',\t'all_sequences',\t'sequence_id',\t'sequence_num',\t'sequence_group',\t'sequence_group_num']], on=\"sequence_id\")\n\n# ------------------\n# Test data\n# ------------------\n\ndf_test['sequence_id'] = df_test[\"target_id\"]\ndf_test[\"sequence_num\"] = df_test[\"all_sequences\"].astype(str).apply(lambda x: \"_\".join(x.split(\"|\")[0:1])[1:])\ndf_test[\"sequence_group\"] = df_test[\"all_sequences\"].astype(str).apply(lambda x: \"_\".join(x.split(\"_\")[0:1])[1:])\ndf_test[\"sequence_group_num\"] = df_test[\"sequence_num\"].apply(lambda x: \"\".join(x.split(\"_\")[1:2]))\n\n\ndf_test_seq = df_test.copy()\n\ndf_test_seq['resname'] = df_test_seq['sequence'].apply(list)  # Преобразуем каждую строку в список букв\ndf_test_seq['target_id'+'sequence'] = df_test_seq['target_id']+df_test_seq['sequence']\n\n# create new dataset for every resid by row\ntest = df_test_seq.explode('resname')  \ntest['resid'] = test.groupby('target_idsequence').cumcount() + 1\ntest['x_1'] = 0\ntest['y_1'\n] = 0\ntest['z_1'] = 0\n\ntest[\"sequence_target\"] = test[\"target_id\"].apply(lambda x: \"\".join(x.split(\"_\")[0]))\ntest[\"sequence_class\"] = 1\ntest[\"ID\"] = test[\"target_id\"]+\"_\"+test[\"resid\"].astype(str)\n\ntest = test[['ID','resname'\t,'resid',\t'x_1',\t'y_1',\t'z_1',\t'sequence_id',\t'sequence_target',\t'sequence_class',\t'sequence',\t'temporal_cutoff',\t'description'\t,'all_sequences', 'sequence_num'\t,'sequence_group'\t,'sequence_group_num']]\n\ndel df_train ,df_valid ,df_train_seq ,df_valid_seq ,df_test \n\ndisplay(train[1:2])\ndisplay(valid[1:2])\ndisplay(test[1:2])","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-03-25T10:20:39.170214Z","iopub.execute_input":"2025-03-25T10:20:39.170738Z","iopub.status.idle":"2025-03-25T10:20:40.179274Z","shell.execute_reply.started":"2025-03-25T10:20:39.170712Z","shell.execute_reply":"2025-03-25T10:20:40.178224Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"target_col = ['x_1', 'y_1', 'z_1']\n\ntrain[target_col] = train[target_col].fillna(train[target_col].mean())","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-03-24T19:26:43.602647Z","iopub.execute_input":"2025-03-24T19:26:43.602969Z","iopub.status.idle":"2025-03-24T19:26:43.616515Z","shell.execute_reply.started":"2025-03-24T19:26:43.602946Z","shell.execute_reply":"2025-03-24T19:26:43.615575Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"#find group from test in train dataset\nprint(\"sequence_group find in train test:\")\nfor i in train['sequence_group'].unique():\n    for j in test['sequence_group'].unique():\n        if i == j:\n            print(i)\nprint(\"sequence_group find in train valid:\")\nfor i in train['sequence_group'].unique():\n    for j in valid['sequence_group'].unique():\n        if i == j:\n            print(i)\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-03-25T10:26:23.011565Z","iopub.execute_input":"2025-03-25T10:26:23.011934Z","iopub.status.idle":"2025-03-25T10:26:23.233507Z","shell.execute_reply.started":"2025-03-25T10:26:23.011896Z","shell.execute_reply":"2025-03-25T10:26:23.232494Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"valid[valid['sequence_group'] == '8UYS'][1:4]","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-03-25T10:30:11.344097Z","iopub.execute_input":"2025-03-25T10:30:11.344440Z","iopub.status.idle":"2025-03-25T10:30:11.360006Z","shell.execute_reply.started":"2025-03-25T10:30:11.344411Z","shell.execute_reply":"2025-03-25T10:30:11.358943Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"from sklearn.base import BaseEstimator, TransformerMixin\nclass AggFeatureExtractor(BaseEstimator, TransformerMixin):\n    \n    def __init__(self, group_col, agg_col, agg_func):\n        self.group_col = group_col\n        self.group_col_name = ''\n        for col in group_col:\n            self.group_col_name += col\n        self.agg_col = agg_col\n        self.agg_func = agg_func\n        self.agg_df = None\n        self.medians = None\n        \n    def fit(self, X, y=None):\n        group_col = self.group_col\n        agg_col = self.agg_col\n        agg_func = self.agg_func\n        \n        self.agg_df = X.groupby(group_col)[agg_col].agg(agg_func)\n        self.agg_df.columns = [f'{self.group_col_name}_{agg}_{_agg_col}' for _agg_col in agg_col for agg in agg_func]\n        self.medians = X[agg_col].median()\n        \n        return self\n    \n    def transform(self, X):\n        group_col = self.group_col\n        agg_col = self.agg_col\n        agg_func = self.agg_func\n        agg_df = self.agg_df\n        medians = self.medians\n        \n        X_merged = pd.merge(X, agg_df, left_on=group_col, right_index=True, how='left')\n        X_merged.fillna(medians, inplace=True)\n        X_agg = X_merged.loc[:, [f'{self.group_col_name}_{agg}_{_agg_col}' for _agg_col in agg_col for agg in agg_func]]\n        \n        return X_agg\n    \n    def fit_transform(self, X, y=None):\n        self.fit(X, y)\n        X_agg = self.transform(X)\n        return X_agg","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-03-25T10:30:24.067432Z","iopub.execute_input":"2025-03-25T10:30:24.067756Z","iopub.status.idle":"2025-03-25T10:30:24.075780Z","shell.execute_reply.started":"2025-03-25T10:30:24.067730Z","shell.execute_reply":"2025-03-25T10:30:24.074728Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"def find_atoms_in_region(x, y, z, x_target, y_target, z_target, tolerance=100):\n    \"\"\"\n    Filter all RNA atoms that are within a given tolerance range around a target coordinate.\n    \"\"\"\n    filtered_atoms =  (min(x_target - tolerance, x_target + tolerance) < x < max(x_target - tolerance, x_target + tolerance)) \\\n        & (min(y_target - tolerance, y_target + tolerance) < y < max(y_target - tolerance, y_target + tolerance)) \\\n        & (min(z_target - tolerance, z_target + tolerance) < z < max(z_target - tolerance, z_target + tolerance))\n\n    return filtered_atoms\n    \ndef euclidean_distance(x1,y1,z1,x2,y2,z2):\n    return np.sqrt((x2 - x1)**2 +\n                   (y2 - y1)**2 +\n                   (z2 - z1)**2)\n\ndef count_nucleotides_map(sequence):\n    return {\n        'A': sequence.count('A'),\n        'C': sequence.count('C'),\n        'G': sequence.count('G'),\n        'U': sequence.count('U')\n    }\ndef count_nucleotides(sequence):\n      return    f\"A{sequence.count('A')}C{sequence.count('C')}G{sequence.count('G')}U{sequence.count('U')}\"\n\ndef get_nucleotide_combinations(seq, start_pos, length):\n        if start_pos + length <= len(seq):\n            result = seq[start_pos:start_pos + length]\n        else:\n            result = '-'  \n        return result\n    \ndef FE (df):\n    df['temporal_cutoff'] = pd.to_datetime(df['temporal_cutoff']).astype('int64') // 10**9\n    df[\"res\"] = df[\"resname\"] + df[\"resid\"].astype(str)\n                                    \n    df[\"length\"] = df[\"sequence\"].str.len()\n    \n    df['A_cnt'] = df['sequence'].astype(str).str.count(\"A\")\n    df['C_cnt'] = df['sequence'].astype(str).str.count(\"C\")\n    df['U_cnt'] = df['sequence'].astype(str).str.count(\"U\")\n    df['G_cnt'] = df['sequence'].astype(str).str.count(\"G\")\n    df['AC_cnt'] = df['sequence'].astype(str).str.count(\"AC\")\n    df['AU_cnt'] = df['sequence'].astype(str).str.count(\"AU\")\n    df['AG_cnt'] = df['sequence'].astype(str).str.count(\"AG\")\n    df['CA_cnt'] = df['sequence'].astype(str).str.count(\"CA\")\n    df['CU_cnt'] = df['sequence'].astype(str).str.count(\"CU\")\n    df['CG_cnt'] = df['sequence'].astype(str).str.count(\"CG\")\n    df['UA_cnt'] = df['sequence'].astype(str).str.count(\"UA\")\n    df['UC_cnt'] = df['sequence'].astype(str).str.count(\"UC\")\n    df['UG_cnt'] = df['sequence'].astype(str).str.count(\"UG\")\n    df['GA_cnt'] = df['sequence'].astype(str).str.count(\"GA\")\n    df['GC_cnt'] = df['sequence'].astype(str).str.count(\"GC\")\n    df['GU_cnt'] = df['sequence'].astype(str).str.count(\"GU\")\n    df['AA_cnt'] = df['sequence'].astype(str).str.count(\"AA\")\n    df['CC_cnt'] = df['sequence'].astype(str).str.count(\"CC\")\n    df['UU_cnt'] = df['sequence'].astype(str).str.count(\"UU\")\n    df['GG_cnt'] = df['sequence'].astype(str).str.count(\"GG\")\n    \n    df['begin_seq'] = df['sequence'].astype(str).str[0]\n    df['end_seq'] = df['sequence'].astype(str).str[-1]\n\n    df['nucleotide_ngram2'] = df.apply(lambda row: get_nucleotide_combinations( row['sequence'], row['resid']-1, 2),   axis=1)\n    df['nucleotide_ngram3'] = df.apply(lambda row: get_nucleotide_combinations( row['sequence'], row['resid']-1, 3),   axis=1)\n    df['nucleotide_ngram4'] = df.apply(lambda row: get_nucleotide_combinations( row['sequence'], row['resid']-1, 4),   axis=1)\n    df['nucleotide_ngram2_prev'] = df['nucleotide_ngram2'].shift(1)\n    df['nucleotide_ngram2_prev'].replace(to_replace=np.nan, value=\"<\", inplace=True)\n  \n    df['nucleotide_counts'] = df['sequence'].apply(count_nucleotides)\n    df['counts_a'] = df['sequence'].apply(lambda x: x.count('A'))\n    df['counts_c'] = df['sequence'].apply(lambda x: x.count('C'))\n    df['counts_g'] = df['sequence'].apply(lambda x: x.count('G'))\n    df['counts_u'] = df['sequence'].apply(lambda x: x.count('U'))\n\n    df[\"GC_content\"] = df[\"sequence\"].apply(\n        lambda seq: (seq.count(\"G\") + seq.count(\"C\")) / len(seq) )\n    df[\"GA_content\"] = df[\"sequence\"].apply(\n        lambda seq: (seq.count(\"G\") + seq.count(\"A\")) / len(seq) )\n    df[\"GU_content\"] = df[\"sequence\"].apply(\n        lambda seq: (seq.count(\"G\") + seq.count(\"U\")) / len(seq) )\n    df[\"CA_content\"] = df[\"sequence\"].apply(\n        lambda seq: (seq.count(\"C\") + seq.count(\"A\")) / len(seq) )\n    df[\"CU_content\"] = df[\"sequence\"].apply(\n        lambda seq: (seq.count(\"C\") + seq.count(\"U\")) / len(seq) )\n    df[\"AU_content\"] = df[\"sequence\"].apply(\n        lambda seq: (seq.count(\"A\") + seq.count(\"U\")) / len(seq) )\n\n    df['resname_prev'] = df['resname'].shift(1)\n    df['resname_next'] = df['resname'].shift(-1)\n    df['resname_prev'].replace(to_replace=np.nan, value=\"<\", inplace=True)\n    df['resname_next'].replace(to_replace=np.nan, value=\">\", inplace=True)\n\n    #df['prev_x'] = df['x_1'].shift(1)\n    #df['prev_y'] = df['y_1'].shift(1)\n    #df['prev_z'] = df['z_1'].shift(1)\n    #df['next_x'] = df['x_1'].shift(-1)\n    #df['next_y'] = df['y_1'].shift(-1)\n    #df['next_z'] = df['z_1'].shift(-1)\n\n    #df[\"distance_from_origin\"] = np.sqrt(df[\"x_1\"]**2 + df[\"y_1\"]**2 + df[\"z_1\"]**2 )\n\n   # df['distance_prev'] = df.apply(lambda row: euclidean_distance( row[\"x_1\"], row[\"y_1\"], row[\"x_1\"],row['prev_x'], row['prev_y'], row['prev_z']), axis=1)\n    #df['distance_next'] = df.apply(lambda row: euclidean_distance( row[\"x_1\"], row[\"y_1\"], row[\"x_1\"],row['next_x'], row['next_y'], row['next_z']), axis=1)\n\n    df[\"chain\"] = df[\"all_sequences\"].astype(str).apply(lambda x: \"_\".join(x.split(\"|\")[1:2])[6:].replace('auth', '').replace(' ', '').strip())\n    df[\"RNA\"] = df[\"all_sequences\"].astype(str).apply(lambda x: \"_\".join(x.split(\"|\")[2:3])[0:])  \n    df[\"RNA_chain\"] = df[\"RNA\"].astype(str).apply(lambda x: \"_\".join(x.split(\"(5'-R(\")[1:])[0:].replace(\")-3')\",\"\").replace(\") -3')\",\"\"))  \n    df[\"add_chain\"] = df[\"all_sequences\"].astype(str).apply(lambda x: \"_\".join(x.split(\"|\")[3:][4:5] )[6:].replace('auth', '').replace(' ', '').strip())\n    df[\"add_RNA\"] = df[\"all_sequences\"].astype(str).apply(lambda x: \"_\".join(x.split(\"|\")[3:][5:6] ))\n    df[\"type_RNA\"] = df[\"all_sequences\"].astype(str).apply(lambda x: \"_\".join(x.split(\"|\")[3:][6:7] ))  \n    df[\"class\"] = df[\"type_RNA\"].astype(str).apply(lambda x: \"_\".join(x.split(\"\\n\")[0:1] ))\n    df[\"class_chain\"] = df[\"type_RNA\"].astype(str).apply(lambda x: \"_\".join(x.split(\"\\n\")[1:] ))\n    \n    df = df.fillna(0)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-03-25T10:31:31.454822Z","iopub.execute_input":"2025-03-25T10:31:31.455221Z","iopub.status.idle":"2025-03-25T10:31:31.480512Z","shell.execute_reply.started":"2025-03-25T10:31:31.455192Z","shell.execute_reply":"2025-03-25T10:31:31.479449Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ------------------\n# Prepare dataset for model \n# ------------------\nprint(f\"Prepare datasets for model \")\n\n#train['set'] = 'train'\n#valid['set'] = 'valid'\n#test['set'] = 'test'\n\n\n# ------------------\n# Create Features\n# ------------------\n\nFE(train)\nFE(valid)\nFE(test)\n\n\nFEATURES = train.columns.tolist()\nprint(f\"train shape :{blu}{train.shape}{clr}, \", f\"valid shape :{blu}{valid.shape}{clr}, \", f\"test shape :{blu}{test.shape}{clr}\")\n\nprint(f\"features shape :{blu}{len(FEATURES)}{clr} \")\nprint(f'{\"-\" * 100}')\n\none_hot_encoded = pd.get_dummies(train['resname'], prefix='nucleotide')\ntrain = pd.concat([train, one_hot_encoded], axis=1)\n\none_hot_encoded = pd.get_dummies(valid['resname'], prefix='nucleotide')\nvalid = pd.concat([valid, one_hot_encoded], axis=1)\n\none_hot_encoded = pd.get_dummies(test['resname'], prefix='nucleotide')\ntest = pd.concat([test, one_hot_encoded], axis=1)\n\n\nprint(f\"train shape :{blu}{train.shape}{clr}, \", f\"valid shape :{blu}{valid.shape}{clr}, \", f\"test shape :{blu}{test.shape}{clr}\")\n\nprint(f'{\"-\" * 100}')\n\n\n\ntrain.head(3)\n\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-03-25T10:39:19.424809Z","iopub.execute_input":"2025-03-25T10:39:19.425206Z","iopub.status.idle":"2025-03-25T10:40:19.656413Z","shell.execute_reply.started":"2025-03-25T10:39:19.425177Z","shell.execute_reply":"2025-03-25T10:40:19.655446Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"For a text description, we will form a bag of words, then reduce the dimension of the new set of features","metadata":{}},{"cell_type":"code","source":"# ------------------\n# Create Features with bag word for Descripction and RNA columns\n# ------------------\ndef bag_of_word(df):\n    count_vectorizer = CountVectorizer(\n        analyzer=\"word\", tokenizer=nltk.word_tokenize,\n        preprocessor=None, stop_words='english', max_features=None)    \n    bag_of_words_combine = count_vectorizer.fit_transform(df['description'])\n    bag_of_words_combine1 = count_vectorizer.fit_transform(df['RNA'])\n    \n    svd = TruncatedSVD(n_components=50, n_iter=30, random_state=12)\n    truncated_bag_of_words_combine = svd.fit_transform(bag_of_words_combine)\n    svd1 = TruncatedSVD(n_components=15, n_iter=25, random_state=12)\n    truncated_bag_of_words_combine1 = svd1.fit_transform(bag_of_words_combine1)\n    \n    add_col = [f'desc_{i}' for i in svd.get_feature_names_out()]\n    add_col1 = [f'rna_{i}' for i in svd1.get_feature_names_out()]\n    \n    bag_df = pd.DataFrame(truncated_bag_of_words_combine, columns = add_col )\n    bag_df1 = pd.DataFrame(truncated_bag_of_words_combine1, columns = add_col1)\n    \n    df = df.reset_index(drop=True)\n    df = pd.concat([df, bag_df1], axis = 1)\n    \n    df = X_combine.reset_index(drop=True)\n    df = pd.concat([df, bag_df], axis = 1)\n    \n    del bag_df,bag_df1\n\n#train = bag_of_word(train)\n#valid = bag_of_word(valid)\n#test = bag_of_word(test)\n\ntrain.head(5)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-03-25T10:50:17.298741Z","iopub.execute_input":"2025-03-25T10:50:17.299152Z","iopub.status.idle":"2025-03-25T10:50:17.339631Z","shell.execute_reply.started":"2025-03-25T10:50:17.299120Z","shell.execute_reply":"2025-03-25T10:50:17.338413Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ------------------\n# Label encode features \n# ------------------\nRMV = ['x_1',\t'y_1',\t'z_1']    \n# from the final dataset we remove all columns calculated based on coordinates and leave those calculated by groups\n\nFEATURES = [c for c in train.columns if not c in RMV]\n\nfor_delete = ['ID','all_sequences',  'description',\"type_RNA\",'sequence_id',\t'sequence_target'] \n\ncat_features = train.drop(RMV, axis=1).select_dtypes(include=['object']).columns.tolist()\nnum_features = train.drop(RMV, axis=1).select_dtypes(exclude=['object']).columns.tolist()\n\nfor_label_encode = [c for c in cat_features if not c in for_delete] \n\nprint(f\"Category features will be encoded: \",for_label_encode)\nprint(f\"\\nNumerical: \",num_features)\n\n#for i,c in enumerate(for_label_encode):\n#    combine = X_combine[c]\n#    combine,_ = pd.factorize(combine)\n#    X_combine[c] = combine.astype(\"float32\")\n  \n# ------------------\n# Label encode features \n# ------------------\n\n\nfrom sklearn.preprocessing import LabelEncoder\n\n# Create dictionary to store label encoders\nlabel_encoders = {}\n\nfor col in for_label_encode:\n    # Handle potential NaN values before encoding\n    train[col] = train[col].astype(str).fillna('missing')\n    \n    # Create and fit label encoder\n    le = LabelEncoder()\n    train[col] = le.fit_transform(train[col])\n    label_encoders[col] = le\n\nfor col in for_label_encode:\n    # Handle potential NaN values before encoding\n    valid[col] = train[col].astype(str).fillna('missing')\n    \n    # Create and fit label encoder\n    le = LabelEncoder()\n    valid[col] = le.fit_transform(valid[col])\n    label_encoders[col] = le\n\nfor col in for_label_encode:\n    # Handle potential NaN values before encoding\n    test[col] = test[col].astype(str).fillna('missing')\n    \n    # Create and fit label encoder\n    le = LabelEncoder()\n    test[col] = le.fit_transform(test[col])\n    label_encoders[col] = le\n\n\n\n# ------------------\n# Result datasets for model \n# ------------------\n\nX_train = train.drop(for_delete, axis =1, errors='ignore')\nX_valid = valid.drop(for_delete, axis =1, errors='ignore')\nX_test = test.drop(for_delete, axis =1, errors='ignore')\n\nFEATURES = [c for c in X_train.columns if not c in RMV]\ntarget_col = ['x_1', 'y_1', 'z_1']\n\nprint(f\"train shape :{blu}{X_train.shape}{clr}, \", f\"valid shape :{blu}{X_valid.shape}{clr}, \", f\"test shape :{blu}{X_test.shape}{clr}\")\n\nprint(f\"X_train ->  isnull :{X_train.isnull().values.sum()}\")\nprint(f\"X_valid ->  isnull :{X_valid.isnull().values.sum()}\")\nprint(f\"X_test -> isnull :{X_test.isnull().values.sum()}\")\n\nX_train.head(5)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-03-25T10:55:34.400181Z","iopub.execute_input":"2025-03-25T10:55:34.400513Z","iopub.status.idle":"2025-03-25T10:55:36.880861Z","shell.execute_reply.started":"2025-03-25T10:55:34.400488Z","shell.execute_reply":"2025-03-25T10:55:36.879977Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"def reduce_mem_usage(df, verbose=True):\n    numerics = ['int16', 'int32', 'int64', 'float16', 'float32', 'float64']\n    start_mem = df.memory_usage().sum() / 1024**2\n    for col in df.columns:\n        col_type = df[col].dtypes\n        if col_type in numerics:\n            c_min = df[col].min()\n            c_max = df[col].max()\n            if str(col_type)[:3] == 'int':\n                if c_min > np.iinfo(np.int8).min and c_max < np.iinfo(np.int8).max:\n                    df[col] = df[col].astype(np.int8)\n                elif c_min > np.iinfo(np.int16).min and c_max < np.iinfo(np.int16).max:\n                    df[col] = df[col].astype(np.int16)\n                elif c_min > np.iinfo(np.int32).min and c_max < np.iinfo(np.int32).max:\n                    df[col] = df[col].astype(np.int32)\n                elif c_min > np.iinfo(np.int64).min and c_max < np.iinfo(np.int64).max:\n                    df[col] = df[col].astype(np.int64)\n            else:\n                if c_min > np.finfo(np.float32).min and c_max < np.finfo(np.float32).max:\n                    df[col] = df[col].astype(np.float32)\n                else:\n                    df[col] = df[col].astype(np.float64)\n\n    end_mem = df.memory_usage().sum() / 1024**2\n    print('Memory usage after optimization is: {:.2f} MB'.format(end_mem))\n    print('Decreased by {:.1f}%'.format(100 * (start_mem - end_mem) / start_mem))\n\n    return df\n\n\nreduce_mem_usage(X_train)\nreduce_mem_usage(X_valid)\nreduce_mem_usage(X_test)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-03-25T10:55:50.284700Z","iopub.execute_input":"2025-03-25T10:55:50.285021Z","iopub.status.idle":"2025-03-25T10:55:50.424866Z","shell.execute_reply.started":"2025-03-25T10:55:50.284996Z","shell.execute_reply":"2025-03-25T10:55:50.423827Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"X_train = X_train.drop('nucleotide_-', axis =1)\nX_train = X_train.drop('nucleotide_X', axis =1)\nFEATURES = [c for c in X_train.columns if not c in RMV]\nX_train","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-03-25T11:04:08.764897Z","iopub.execute_input":"2025-03-25T11:04:08.765372Z","iopub.status.idle":"2025-03-25T11:04:08.905752Z","shell.execute_reply.started":"2025-03-25T11:04:08.765330Z","shell.execute_reply":"2025-03-25T11:04:08.904148Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"X_train.info()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-03-25T10:55:55.511329Z","iopub.execute_input":"2025-03-25T10:55:55.511644Z","iopub.status.idle":"2025-03-25T10:55:55.560615Z","shell.execute_reply.started":"2025-03-25T10:55:55.511620Z","shell.execute_reply":"2025-03-25T10:55:55.559608Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"#FEATURES = ['resname','resid'] + list_f + add_col + add_col1\nprint(len(FEATURES))\n\nFEATURES","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-03-25T10:56:12.429426Z","iopub.execute_input":"2025-03-25T10:56:12.429779Z","iopub.status.idle":"2025-03-25T10:56:12.436410Z","shell.execute_reply.started":"2025-03-25T10:56:12.429747Z","shell.execute_reply":"2025-03-25T10:56:12.435522Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"#  Data Modeling","metadata":{}},{"cell_type":"markdown","source":"# XGBoost","metadata":{}},{"cell_type":"code","source":"from sklearn.model_selection import KFold,cross_val_score\nfrom sklearn.linear_model import LinearRegression\nfrom sklearn.multioutput import RegressorChain\nfrom sklearn.model_selection import RepeatedKFold\nfrom sklearn.metrics import mean_absolute_error\nfrom xgboost import plot_importance\nfrom xgboost import XGBRegressor\nimport xgboost as xgb\nprint(f\"XGBoost version\",xgb.__version__)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-03-25T10:56:19.018185Z","iopub.execute_input":"2025-03-25T10:56:19.018523Z","iopub.status.idle":"2025-03-25T10:56:19.249199Z","shell.execute_reply.started":"2025-03-25T10:56:19.018496Z","shell.execute_reply":"2025-03-25T10:56:19.248164Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"clf = 'xgb'\n#clf = 'lnr'","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-03-25T10:56:21.121360Z","iopub.execute_input":"2025-03-25T10:56:21.121706Z","iopub.status.idle":"2025-03-25T10:56:21.125715Z","shell.execute_reply.started":"2025-03-25T10:56:21.121676Z","shell.execute_reply":"2025-03-25T10:56:21.124666Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"\nmodels = {}\nfor coord in target_col:\n    print(f'Training model for {coord}...')\n    model = xgb.XGBRegressor(\n        objective='reg:squarederror',\n        n_estimators=1000,\n        max_depth=7,\n        learning_rate=0.1,\n        subsample=0.8,\n        colsample_bytree=0.8\n    )\n    \n    # Filter out rows with missing coordinates\n    mean_val = X_train[coord].mean()\n    X_train[coord].fillna(mean_val, inplace=True)\n    valid_idx = X_train[coord].notna()  # Now just indicates all rows are valid\n    model.fit(X_train[valid_idx][FEATURES], X_train[valid_idx][coord])\n    models[coord] = model\n\n\n# Predict coordinates for test set\npred = X_test[FEATURES].copy()\nfor coord in coord_targets:\n    # Use trained model to make predictions\n    pred[coord] = models[coord].predict(X_test[FEATURES])\n    \n# Display sample predictions\nprint(\"Test predictions sample:\")\npred.head()","metadata":{"execution":{"iopub.status.busy":"2025-03-24T19:27:46.284535Z","iopub.status.idle":"2025-03-24T19:27:46.284943Z","shell.execute_reply":"2025-03-24T19:27:46.284782Z"}}},{"cell_type":"code","source":"%%time\nclf = 'xgb'\n#clf = 'lnr'\nxgb_params = { 'learning_rate': 0.001,\n                        #'min_child_weight': 0.3,\n                        'early_stopping_rounds': 300,\n                        #multi_strategy=\"multi_output_tree\",\n                        #'num_target':  3,\n                        'booster': 'gbtree',\n                        'tree_method': \"hist\",\n                        'n_estimators': 3000,\n                        'subsample': 0.8,\n                        'colsample_bytree': 0.8,\n                        'n_jobs': 10,\n                        'max_depth': 15,\n                        'eval_metric':  'rmse' }\n\nFOLDS = 5\nkf = KFold(n_splits=FOLDS, shuffle=True, random_state=42)\n\noof = np.zeros((len(X_train[FEATURES]),3))  # x,y,z\npred = np.zeros((len(X_test[FEATURES]),3))\nX1_test = X_test[FEATURES]\n\n#  K FOLD x\nfor i, (train_index, valid_index) in enumerate(kf.split(X_train[FEATURES])):\n    print(f\"### X Fold {i+1} ###\")\n    \n    X1_train = X_train.loc[train_index,FEATURES].reset_index(drop=True).copy()\n    X1_valid = X_train.loc[valid_index,FEATURES].reset_index(drop=True).copy()\n    y1_train = X_train.loc[train_index,['x_1']]\n    y1_valid = X_train.loc[valid_index,['x_1']]\n    X1_test = X1_test[FEATURES].reset_index(drop=True).copy()\n    print(f\"X train shape :{blu}{X1_train.shape}{clr}, \", f\"X valid shape :{blu}{X1_valid.shape}{clr}, \", f\"X test shape :{blu}{X1_test.shape}{clr}\")\n    print(f\"x train shape :{blu}{y1_train.shape}{clr}, \", f\"x valid shape :{blu}{y1_valid.shape}{clr}, \")\n\n    \n    xgb_x = XGBRegressor(**xgb_params)\n    xgb_x.fit(X1_train, y1_train, eval_set=[(X1_valid, y1_valid)],  verbose=300 )\n    \n    oof[valid_index,0] = xgb_x.predict(X1_valid)\n    pred[:,0] += xgb_x.predict(X1_test)\n    print('MAE x valid: %s' % mean_absolute_error(oof[valid_index,0], y1_valid))\n    \npred[:,0] /= FOLDS\nprint(\"~result:\")\ndisplay(pred)\n\nX1_test['x_1'] = pred[:,0]\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-03-25T11:04:14.285751Z","iopub.execute_input":"2025-03-25T11:04:14.286045Z","iopub.status.idle":"2025-03-25T11:10:15.187460Z","shell.execute_reply.started":"2025-03-25T11:04:14.286021Z","shell.execute_reply":"2025-03-25T11:10:15.185667Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"#  K FOLD y\nfor i, (train_index, valid_index) in enumerate(kf.split(X_train[FEATURES])):\n    print(f\"### Y Fold {i+1} ###\")\n    \n    X2_train = X_train.loc[train_index,FEATURES + ['x_1']].reset_index(drop=True).copy()\n    \n    X2_valid = X_train.loc[valid_index,FEATURES + ['x_1']].reset_index(drop=True).copy()\n    y2_train = X_train.loc[train_index,['y_1']]\n    y2_valid = X_train.loc[valid_index,['y_1']]\n    X2_test = X1_test[FEATURES + ['x_1']].reset_index(drop=True).copy()\n    print(f\"X2 train shape :{blu}{X2_train.shape}{clr}, \", f\"X2 valid shape :{blu}{X2_valid.shape}{clr}, \", f\"X2 test shape :{blu}{X2_test.shape}{clr}\")\n    print(f\"y2 train shape :{blu}{y2_train.shape}{clr}, \", f\"y2 valid shape :{blu}{y2_valid.shape}{clr}, \")\n\n    xgb_y = XGBRegressor(**xgb_params)\n    xgb_y.fit(X2_train, y2_train, eval_set=[(X2_valid, y2_valid)],  verbose=300 )\n    oof[valid_index,1] = xgb_y.predict(X2_valid)\n    pred[:,1] += xgb_y.predict(X2_test)\n    print('MAE y valid: %s' % mean_absolute_error(oof[valid_index,1], y2_valid))\n    \npred[:,1] /= FOLDS\nprint(\"~result:\")\ndisplay(pred)\n\nX1_test['y_1'] = pred[:,1]\n\n\n#  K FOLD z\nfor i, (train_index, valid_index) in enumerate(kf.split(X_train[FEATURES])):\n    print(f\"### Z Fold {i+1} ###\")\n    \n    X3_train = X_train.loc[train_index,FEATURES + ['x_1','y_1']].reset_index(drop=True).copy()\n    X3_valid = X_train.loc[valid_index,FEATURES + ['x_1','y_1']].reset_index(drop=True).copy()\n    y3_train = X_train.loc[train_index,['z_1']]\n    y3_valid = X_train.loc[valid_index,['z_1']]\n    X3_test = X1_test[FEATURES + ['x_1','y_1']].reset_index(drop=True).copy()\n    print(f\"X3 train shape :{blu}{X3_train.shape}{clr}, \", f\"X3 valid shape :{blu}{X3_valid.shape}{clr}, \", f\"X3 test shape :{blu}{X3_test.shape}{clr}\")\n    print(f\"y3 train shape :{blu}{y3_train.shape}{clr}, \", f\"y3 valid shape :{blu}{y3_valid.shape}{clr}, \")\n    xgb_z = XGBRegressor(**xgb_params)\n    xgb_z.fit(X3_train, y3_train, eval_set=[(X3_valid, y3_valid)],  verbose=300 )\n    oof[valid_index,2] = xgb_z.predict(X3_valid)\n    pred[:,2] += xgb_z.predict(X3_test)\n    print('MAE z valid: %s' % mean_absolute_error(oof[valid_index,2], y3_valid))\n    \n    \npred[:,2] /= FOLDS\nprint(\"~result:\")\ndisplay(pred)\n\nX1_test['z_1'] = pred[:,2]\n\n\n    \n    # CLEAR MEMORY\n#del X1_train, X1_valid, X2_train, X2_valid,X3_train, X3_valid\n#del y1_train, y1_valid, y2_train, y2_valid, y3_train, y3_valid\n#if i != FOLDS-1: del xgb_x,xgb_y,xgb_z\n    \n","metadata":{"execution":{"iopub.status.busy":"2025-03-25T11:10:43.201338Z","iopub.execute_input":"2025-03-25T11:10:43.201658Z","iopub.status.idle":"2025-03-25T11:23:29.797843Z","shell.execute_reply.started":"2025-03-25T11:10:43.201634Z","shell.execute_reply":"2025-03-25T11:23:29.796955Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"<pre>from sklearn.inspection import permutation_importance\nperm_importance = permutation_importance(xgb_z, X3_valid, y3_valid)\nsorted_idx = perm_importance.importances_mean.argsort()[:50]\n#plt.barh(X3_valid.columns[sorted_idx], perm_importance.importances_mean[sorted_idx], )\n#plt.xlabel(\"Permutation Importance\")\nprint(\"done\")\n# Figure Size\nfig, ax = plt.subplots(figsize=(16, 20))\n \n# Horizontal Bar Plot\nax.barh(X3_valid.columns[sorted_idx], perm_importance.importances_mean[sorted_idx])\n\nax = plt.gca() #Getting the current axis\nax.spines['bottom'].set_visible(False) \nax.spines['top'].set_visible(False) \nax.spines['left'].set_visible(False) \nax.spines['right'].set_visible(False) \n\nplt.setp(ax.spines.values(), visible=False) \n# Remove x, y Ticks\nax.xaxis.set_ticks_position('none')\nax.yaxis.set_ticks_position('none')\n \n# Add padding between axes and labels\nax.xaxis.set_tick_params(pad=5)\nax.yaxis.set_tick_params(pad=10)\n \n# Add x, y gridlines\nax.grid( color='grey',\n        linestyle='-.', linewidth=0.5,\n        alpha=0.2)\n \n# Show top values\nax.invert_yaxis()\n \n# Add annotation to bars\nfor i in ax.patches:\n    plt.text(i.get_width(), i.get_y(),\n             str(round((i.get_width()), 3)),\n             fontsize=8, \n             color='grey')\n \n# Add Plot Title\nax.set_title('Permutation Importance',\n             loc='left', )\n \n# Add Text watermark\nfig.text(0.9, 0.15, 'xgb_z', fontsize=12,\n         color='grey', ha='right', va='top',\n         alpha=0.7)\n \n# Show Plot\nplt.show()\n</pre>","metadata":{"execution":{"iopub.status.busy":"2025-03-23T10:51:35.372226Z","iopub.execute_input":"2025-03-23T10:51:35.372637Z","iopub.status.idle":"2025-03-23T10:53:03.360229Z","shell.execute_reply.started":"2025-03-23T10:51:35.372603Z","shell.execute_reply":"2025-03-23T10:53:03.358554Z"}}},{"cell_type":"code","source":"sorted_idx = np.argsort(xgb_x.feature_importances_)[::-1]\nfor index in sorted_idx:\n    print([X1_test.columns[index], xgb_x.feature_importances_[index]]) \nfig, ax = plt.subplots(figsize=(16, 20))\n \n# Horizontal Bar Plot\nax.barh(X1_test.columns[sorted_idx][:30], xgb_x.feature_importances_[sorted_idx][:30])\n\nax = plt.gca() #Getting the current axis\nax.spines['bottom'].set_visible(False) \nax.spines['top'].set_visible(False) \nax.spines['left'].set_visible(False) \nax.spines['right'].set_visible(False) \n\nplt.setp(ax.spines.values(), visible=False) \n# Remove x, y Ticks\nax.xaxis.set_ticks_position('none')\nax.yaxis.set_ticks_position('none')\n \n# Add padding between axes and labels\nax.xaxis.set_tick_params(pad=5)\nax.yaxis.set_tick_params(pad=10)\n \n# Add x, y gridlines\nax.grid( color='grey',\n        linestyle='-.', linewidth=0.5,\n        alpha=0.2)\n \n# Show top values\nax.invert_yaxis()\n \n# Add annotation to bars\nfor i in ax.patches:\n    plt.text(i.get_width(), i.get_y(),\n             str(round((i.get_width()), 3)),\n             fontsize=8, \n             color='grey')\n \n# Add Plot Title\nax.set_title('Permutation Importance',\n             loc='left', )\n \n# Add Text watermark\nfig.text(0.9, 0.15, 'xgb_z', fontsize=12,\n         color='grey', ha='right', va='top',\n         alpha=0.7)\n \n# Show Plot\nplt.show()\n#plot_importance(xgb_x, max_num_features = 15)\n#plt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-03-25T11:23:45.860025Z","iopub.execute_input":"2025-03-25T11:23:45.860367Z","iopub.status.idle":"2025-03-25T11:23:47.268179Z","shell.execute_reply.started":"2025-03-25T11:23:45.860342Z","shell.execute_reply":"2025-03-25T11:23:47.267183Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"feature_important = xgb_x.get_booster().get_score(importance_type='weight')\nkeys = list(feature_important.keys())\nvalues = list(feature_important.values())\n\ndata = pd.DataFrame(data=values, index=keys, columns=[\"score\"]).sort_values(by = \"score\", ascending=True)\ndata.nlargest(40, columns=\"score\").plot(kind='barh', figsize = (20,10)) ## plot top 40 features","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-03-25T11:24:00.983975Z","iopub.execute_input":"2025-03-25T11:24:00.984358Z","iopub.status.idle":"2025-03-25T11:24:01.541537Z","shell.execute_reply.started":"2025-03-25T11:24:00.984331Z","shell.execute_reply":"2025-03-25T11:24:01.540668Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"y_valid = X_valid[target_col].copy()\ny_valid = y_valid.reset_index(drop = True)\ny_pred = pd.DataFrame(pred[:,0:3], columns = target_col)\n\ndef calculate_tm_score_exact(pred_coords, true_coords):\n    \"\"\"\n    Implementation closer to the official method used by US-align.\n    \"\"\"\n    # Remove padding\n    mask = ~np.all(true_coords == 0)\n    pred = pred_coords[mask]\n    true = true_coords[mask]\n    \n    Lref = len(true_coords)\n    \n    # Define d0 exactly as in the evaluation formula\n    if Lref >= 30:\n        d0 = 0.6 * np.sqrt(Lref - 0.5) - 2.5\n    elif Lref >= 24:\n        d0 = 0.7\n    elif Lref >= 20:\n        d0 = 0.6\n    elif Lref >= 16:\n        d0 = 0.5\n    elif Lref >= 12:\n        d0 = 0.4\n    else:\n        d0 = 0.3\n    \n    # Normalize structures\n    pred_centered = pred - np.mean(pred, axis=0)\n    true_centered = true - np.mean(true, axis=0)\n    \n    # Covariance matrix for optimal rotation\n    covariance = np.dot(pred_centered.T, true_centered)\n    U, S, Vt = np.linalg.svd(covariance)\n    rotation = np.dot(U, Vt)\n    \n    # Apply rotation\n    pred_aligned = np.dot(pred_centered, rotation)\n    \n    # Calculate distances\n    distances = np.sqrt(np.sum((pred_aligned - true_centered) ** 2, axis=1))\n    \n    # Calculate TM-score terms\n    tm_terms = 1.0 / (1.0 + (distances / d0) ** 2)\n    tm_score = np.sum(tm_terms) / Lref\n    \n    return float(tm_score)\n    \ndef calculate_tm_score(pred_coords, true_coords, d0_scale=1.24):\n    \"\"\"\n    Calculates a robust approximation of the TM-score between predicted and true coordinates.\n    Adds protection against division by zero and NaN.\n    \"\"\"\n    # Remove padding (rows with zeros) from true structures\n    mask = ~np.all(true_coords == 0)\n    pred = pred_coords[mask]\n    true = true_coords[mask]\n    \n    L = len(true_coords)\n\n    if L < 3:\n        return 0.0\n    \n    # Define d0 based on L (values adapted for RNA)\n    if L >= 30:\n        d0 = 0.6 * np.sqrt(L - 0.5) - 2.5\n        d0 = max(0.1, d0)\n    elif L >= 24:\n        d0 = 0.7\n    elif L >= 20:\n        d0 = 0.6\n    elif L >= 16:\n        d0 = 0.5\n    elif L >= 12:\n        d0 = 0.4\n    else:\n        d0 = 0.3\n    \n    distances = np.sqrt(np.sum((pred - true) ** 2, axis=1))\n    tm_terms = 1.0 / (1.0 + (distances / (d0 + 1e-8)) ** 2)\n    tm_score = np.sum(tm_terms) / L\n    return float(tm_score)\n\ny_valid = np.nan_to_num(y_valid, nan=0.0)\ny_pred = np.nan_to_num(y_pred, nan=0.0)\ntm_scores = []\ntm_scores1 = []\nfor i in range(len(X_valid[target_col])):\n       \n        tm = calculate_tm_score(y_pred[i], y_valid[i])\n        tm1 = calculate_tm_score_exact(y_pred[i], y_valid[i])\n        tm_scores.append(tm)\n        tm_scores1.append(tm1)\navg_tm_score = np.mean(tm_scores)\nprint(f\"Approximate average TM-score: {avg_tm_score:.8f}\")\navg_tm_score1 = np.mean(tm_scores1)\nprint(f\"Approximate average TM-score: {avg_tm_score1:.8f}\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-03-25T11:24:05.474603Z","iopub.execute_input":"2025-03-25T11:24:05.475107Z","iopub.status.idle":"2025-03-25T11:24:05.846515Z","shell.execute_reply.started":"2025-03-25T11:24:05.475042Z","shell.execute_reply":"2025-03-25T11:24:05.845709Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"print(pred[:, 0].shape)\npred[:, 0]","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-03-25T11:24:14.493538Z","iopub.execute_input":"2025-03-25T11:24:14.493856Z","iopub.status.idle":"2025-03-25T11:24:14.500669Z","shell.execute_reply.started":"2025-03-25T11:24:14.493832Z","shell.execute_reply":"2025-03-25T11:24:14.499677Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"\nsub = pd.read_csv('../input/stanford-rna-3d-folding/sample_submission.csv')\ncol_name = ['x_1',\t'y_1',\t'z_1',\t'x_2',\t'y_2',\t'z_2',\t'x_3',\t'y_3',\t'z_3',\t'x_4',\t'y_4',\t'z_4',\t'x_5',\t'y_5',\t'z_5']\nsub = sub.drop(col_name, axis =1)\nsub['x_1'] = pred[:, 0]\nsub['y_1'] = pred[:, 1]\nsub['z_1'] = pred[:, 2]\n\nsub['x_2'] = 0.0\nsub['y_2'] = 0.0\nsub['z_2'] = 0.0\n\nsub['x_3'] = 0.0\nsub['y_3'] = 0.0\nsub['z_3'] = 0.0\n\nsub['x_4'] = 0.0\nsub['y_4'] = 0.0\nsub['z_4'] = 0.0\n\nsub['x_5'] = 0.0\nsub['y_5'] = 0.0\nsub['z_5'] = 0.0\n\n#submission[['x_1','y_1','z_1','x_2','y_2','z_2','x_3','y_3','z_3','x_4','y_4','z_4','x_5','y_5','z_5']] = test_predictions[['x_1','y_1','z_1','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\n#for i in range(1,6):\n #   columns+=[f\"x_{i}\"]\n#    columns+=[f\"y_{i}\"]\n #   columns+=[f\"z_{i}\"]\n\nVER = 1\nsub.to_csv('submission.csv', index=False)\nsub","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-03-25T11:25:10.532593Z","iopub.execute_input":"2025-03-25T11:25:10.532895Z","iopub.status.idle":"2025-03-25T11:25:10.603881Z","shell.execute_reply.started":"2025-03-25T11:25:10.532871Z","shell.execute_reply":"2025-03-25T11:25:10.602810Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Visualize predictions","metadata":{}},{"cell_type":"code","source":"# Plot the results\nplt.figure(figsize=(15, 5))\ns = 100\na = 0.4\n\n\nplt.subplot(1, 6, 1)  \nplt.scatter(valid['x_1'],valid['y_1'], edgecolor=\"k\",c=\"cornflowerblue\", s=s,alpha=a)\nplt.title('X Y projection valid')\n\nplt.subplot(1, 6, 2)  \nplt.scatter(valid['y_1'], valid['z_1'], edgecolor=\"k\",c=\"cornflowerblue\", s=s,alpha=a)\nplt.title('Y Z projection valid')\n\nplt.subplot(1, 6, 3) \nplt.scatter(valid['x_1'],valid['z_1'], edgecolor=\"k\",c=\"cornflowerblue\", s=s,alpha=a)\nplt.title('X Z projection valid')\n\nplt.subplot(1, 6, 4)  \nplt.scatter(pred[:, 0],pred[:, 1], edgecolor=\"k\",c=\"cornflowerblue\", s=s,alpha=a)\nplt.title('X Y projection')\n\nplt.subplot(1, 6, 5)  \nplt.scatter(pred[:, 1],pred[:, 2], edgecolor=\"k\",c=\"cornflowerblue\", s=s,alpha=a)\nplt.title('Y Z projection')\n\nplt.subplot(1, 6, 6) \nplt.scatter(pred[:, 0],pred[:, 1], edgecolor=\"k\",c=\"cornflowerblue\", s=s,alpha=a)\nplt.title('X Z projection')\n\nplt.tight_layout()  \n\nplt.legend()\nplt.show()\n\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-03-25T11:25:20.331950Z","iopub.execute_input":"2025-03-25T11:25:20.332416Z","iopub.status.idle":"2025-03-25T11:25:21.540523Z","shell.execute_reply.started":"2025-03-25T11:25:20.332381Z","shell.execute_reply":"2025-03-25T11:25:21.539484Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"def plot_structure(df: pd.DataFrame, sequence_id: str) -> None:\n    sequence_df = df[df[\"sequence_id\"] == sequence_id]\n    sequence_points = sequence_df[[\"x_1\", \"y_1\", \"z_1\", \"resname\"]]\n    seq_lst = sequence_df['resname'].to_list()\n    seq_str = ''.join(seq_lst)\n    print(seq_str)\n    #print(sequence_points)\n    \n    colors = {\"A\": \"red\", \"G\": \"blue\", \"C\": \"green\", \"U\": \"orange\"}\n    fig = go.Figure()\n    \n    for resname, color in colors.items():\n        subset = sequence_df[sequence_df[\"resname\"] == resname]\n        fig.add_trace(go.Scatter3d(\n            x=subset[\"x_1\"], y=subset[\"y_1\"], z=subset[\"z_1\"],\n            mode='markers',\n            marker=dict(size=5, color=color),\n            name=resname,\n            opacity=0.8\n        ))\n    \n    fig.add_trace(go.Scatter3d(\n        x=sequence_df[\"x_1\"], y=sequence_df[\"y_1\"], z=sequence_df[\"z_1\"],\n        mode='lines',\n        line=dict(color='gray', width=2),\n        name='RNA Backbone'\n    ))\n    \n    fig.update_layout(\n            scene=dict(xaxis_title='X', yaxis_title='Y', zaxis_title='Z'),\n            title=f'3D RNA Structure of sequence {sequence_id}',\n        )\n            \n    fig.show(renderer=\"iframe\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-03-25T11:25:30.144241Z","iopub.execute_input":"2025-03-25T11:25:30.144550Z","iopub.status.idle":"2025-03-25T11:25:30.151831Z","shell.execute_reply.started":"2025-03-25T11:25:30.144526Z","shell.execute_reply":"2025-03-25T11:25:30.150939Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"sequence_id = \"8UYS\"\nsequence_df = valid[valid[\"sequence_group\"] == sequence_id].copy()\n\nsequence_points = sequence_df[[\"x_1\", \"y_1\", \"z_1\", \"resname\"]]\nseq_lst = sequence_df['resname'].to_list()\nseq_str = ''.join(seq_lst)\nsequence_df[0:1]\nplot_structure(sequence_df, sequence_df[\"sequence_id\"].max())","metadata":{"execution":{"iopub.status.busy":"2025-03-25T11:25:46.031394Z","iopub.execute_input":"2025-03-25T11:25:46.031695Z","iopub.status.idle":"2025-03-25T11:25:46.088634Z","shell.execute_reply.started":"2025-03-25T11:25:46.031672Z","shell.execute_reply":"2025-03-25T11:25:46.087782Z"}}},{"cell_type":"code","source":"sub_id = sub[\"ID\"].apply(lambda x: \"_\".join(x.split(\"_\")[:-1])).unique()\n#sub[\"sequence_group\"] = sub[\"ID\"].apply(lambda x: \"_\".join(x.split(\"_\")[:-1]))\n#sub[\"sequence_id\"] = sub[\"ID\"].apply(lambda x: \"_\".join(x.split(\"_\")[:-1]))\n\n#sequence_id = \"R1149\"\n#sequence_df = sub[sub[\"sequence_group\"] == sequence_id].copy()\n\ndef print_plot_str(c):\n    #for c in sub[\"sequence_group\"].unique():\n    print(\"----------------------------------------------------------------\")\n    print(c)\n    sequence_df = sub[sub[\"sequence_group\"] == c].copy()\n    #display(sequence_df.head(3))\n    plot_structure(sequence_df, c)\n\nsub[\"sequence_group\"] = sub[\"ID\"].apply(lambda x: \"_\".join(x.split(\"_\")[:-1]))\nsub[\"sequence_id\"] = sub[\"ID\"].apply(lambda x: \"_\".join(x.split(\"_\")[:-1]))\nprint_plot_str('R1107')\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-03-25T11:25:49.440757Z","iopub.execute_input":"2025-03-25T11:25:49.441083Z","iopub.status.idle":"2025-03-25T11:25:49.492923Z","shell.execute_reply.started":"2025-03-25T11:25:49.441037Z","shell.execute_reply":"2025-03-25T11:25:49.491992Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"print_plot_str('R1108')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-03-25T11:25:57.715535Z","iopub.execute_input":"2025-03-25T11:25:57.715902Z","iopub.status.idle":"2025-03-25T11:25:57.758837Z","shell.execute_reply.started":"2025-03-25T11:25:57.715874Z","shell.execute_reply":"2025-03-25T11:25:57.758031Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"print_plot_str('R1156')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-03-25T11:26:07.447547Z","iopub.execute_input":"2025-03-25T11:26:07.447924Z","iopub.status.idle":"2025-03-25T11:26:07.497562Z","shell.execute_reply.started":"2025-03-25T11:26:07.447882Z","shell.execute_reply":"2025-03-25T11:26:07.496721Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"print_plot_str('R1136')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-03-24T19:22:23.894998Z","iopub.status.idle":"2025-03-24T19:22:23.895284Z","shell.execute_reply":"2025-03-24T19:22:23.895167Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"print_plot_str('R1126')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-03-24T19:22:23.895926Z","iopub.status.idle":"2025-03-24T19:22:23.896281Z","shell.execute_reply":"2025-03-24T19:22:23.896159Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"print_plot_str('R1116')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-03-24T19:22:23.897484Z","iopub.status.idle":"2025-03-24T19:22:23.897835Z","shell.execute_reply":"2025-03-24T19:22:23.897666Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"print_plot_str('R1138')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-03-24T19:22:23.898769Z","iopub.status.idle":"2025-03-24T19:22:23.899143Z","shell.execute_reply":"2025-03-24T19:22:23.899009Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"print_plot_str('R1117v2')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-03-24T19:22:23.900309Z","iopub.status.idle":"2025-03-24T19:22:23.900635Z","shell.execute_reply":"2025-03-24T19:22:23.900488Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"print_plot_str('R1128')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-03-24T19:22:23.901779Z","iopub.status.idle":"2025-03-24T19:22:23.902139Z","shell.execute_reply":"2025-03-24T19:22:23.902003Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"print_plot_str('R1190')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-03-24T19:22:23.903266Z","iopub.status.idle":"2025-03-24T19:22:23.903625Z","shell.execute_reply":"2025-03-24T19:22:23.903501Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"print_plot_str('R1149')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-03-24T19:22:23.904278Z","iopub.status.idle":"2025-03-24T19:22:23.904576Z","shell.execute_reply":"2025-03-24T19:22:23.904457Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"print_plot_str('R1189')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-03-24T19:22:23.905442Z","iopub.status.idle":"2025-03-24T19:22:23.905814Z","shell.execute_reply":"2025-03-24T19:22:23.905623Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"sequence_df = sub[sub[\"sequence_group\"] == 'R1138'].copy()\nX = sequence_df['x_1']\nY = sequence_df['y_1']\nZ = sequence_df['z_1']\n\n\nplt.figure(figsize=(15, 5))\n\nplt.subplot(1, 3, 1)  \nsns.scatterplot(data=sequence_df, x=X, y=Y, s=Z*Z.max()/X.max(),hue = 'resname')\nplt.title('X Y projection')\n\nplt.subplot(1, 3, 2)  \nsns.scatterplot(data=sequence_df, x=Y, y=Z, s=X*X.max()/Y.max(),hue = 'resname')\nplt.title('Y Z projection')\n\nplt.subplot(1, 3, 3) \nsns.scatterplot(data=sequence_df, x=X, y=Z, s=Y*Y.max()/X.max(),hue = 'resname')\nplt.title('X Z projection')\n\nplt.tight_layout()  \nplt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-03-25T11:26:16.392016Z","iopub.execute_input":"2025-03-25T11:26:16.392408Z","iopub.status.idle":"2025-03-25T11:26:17.391020Z","shell.execute_reply.started":"2025-03-25T11:26:16.392366Z","shell.execute_reply":"2025-03-25T11:26:17.389926Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"g = sns.pairplot(data=sequence_df[['x_1', 'y_1', 'z_1', 'resname']],\n             hue = 'resname', height=2.5,diag_kind=\"kde\",\n             diag_kws = {'color' : default_color_1},\n             plot_kws = {'s' : 15, \n                         'alpha' : 0.5,\n                         'color' : default_color_1})\ng.map_lower(sns.kdeplot, levels=1, color=\".2\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-03-25T11:26:19.716721Z","iopub.execute_input":"2025-03-25T11:26:19.717083Z","iopub.status.idle":"2025-03-25T11:26:24.877279Z","shell.execute_reply.started":"2025-03-25T11:26:19.717030Z","shell.execute_reply":"2025-03-25T11:26:24.876276Z"}},"outputs":[],"execution_count":null}]}