{"metadata":{"kaggle":{"accelerator":"none","dataSources":[{"sourceId":51294,"databundleVersionId":6923401,"sourceType":"competition"},{"sourceId":6900449,"sourceType":"datasetVersion","datasetId":3963795}],"dockerImageVersionId":30587,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":false},"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.10.13"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# Experiments on data preparation and cleaning\n\n### The data set used is taken from here: [Generating a Dataset with structure field](https://www.kaggle.com/code/konstantinboyko/generating-a-dataset-with-structure-field)\n\n## This example describes various attempts to prepare and clean data for further experiments.\n\n### `The construction of heat maps requires a large amount of RAM; for this, a local computer with 64 GB of RAM was used.`","metadata":{}},{"cell_type":"code","source":"import pandas as pd\nimport numpy as np\nimport seaborn as sns\nfrom pathlib import Path\n\nimport matplotlib.pyplot as plt\nimport matplotlib.mlab as mlab\nimport matplotlib\nplt.style.use('ggplot')\nfrom matplotlib.pyplot import figure\n\n%matplotlib inline\nmatplotlib.rcParams['figure.figsize'] = (12,8)\n\npd.options.mode.chained_assignment = None\n\n#%load_ext cudf.pandas\n\nclass APP:\n    path_src = Path('/kaggle/input/stanford-ribonanza-rna-folding-prepare')\n    path_trg = Path('/kaggle/working')\n    file_src = path_src / 'train_data_struct_ext.parquet'\n    file_trg = path_trg / 'train_struct_ext_clear_v01.parquet'\n    colours = ['#000099', '#ffff00']  # define colors: yellow - missing data, blue - not missing\n    show_heatmap = True\n    show_hist_skip_val = True\n    show_hist_seq_len = True\n    show_hist_react_error = True\n\nclass CFG:\n    drop_col_activity = True\n    drop_col_act_max = True\n    remove_react_nan = True\n    remove_error_nan = True\n    cut_sequence = True\n    remove_seq_not_177 = False\n    remove_duplicates = True\n    dropna = True\n    merge_react_in_one_col = True\n    merge_by_columns = True\n    skip_seq_beg = 26\n    skip_seq_end = 21\n    drop_col_end_min = 135\n    drop_col_end_max = 126\n\n# reading data\ndf = pd.read_parquet(APP.file_src)\nprint('df=',df.shape, '2A3_MaP=', len(df[df['experiment_type']=='2A3_MaP']), 'DMS_MaP=',len(df[df['experiment_type']=='DMS_MaP']))\n\n# shape and data types of the data\nprint(df.shape)\nprint(df.dtypes)\n\n# selection of non-numeric columns\ndf_non_numeric = df.select_dtypes(exclude=[np.number])\nnon_numeric_cols = df_non_numeric.columns.values\nprint(non_numeric_cols)\n\n# selection of numeric columns\ndf_numeric = df.select_dtypes(include=[np.number])\nnumeric_cols = df_numeric.columns.values\nprint(numeric_cols)","metadata":{"execution":{"iopub.execute_input":"2023-11-30T09:39:06.114747Z","iopub.status.busy":"2023-11-30T09:39:06.114329Z","iopub.status.idle":"2023-11-30T09:39:30.390307Z","shell.execute_reply":"2023-11-30T09:39:30.388914Z","shell.execute_reply.started":"2023-11-30T09:39:06.114711Z"},"collapsed":true,"jupyter":{"outputs_hidden":true},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Heat map of missing values\n\nWhen there are not very many features in the set, visualize the missing values using a heat map.\n\nThe map below shows the pattern of missing values for the set features. Features are located along the horizontal axis, and the number of records/rows is located along the vertical axis. Yellow color indicates missing data.","metadata":{}},{"cell_type":"code","source":"react_cols = df.filter(like='reactivity_0').columns.to_list()\nif APP.show_heatmap:\n    sns.heatmap(df[react_cols].isnull(), cmap=sns.color_palette(APP.colours))","metadata":{"execution":{"iopub.execute_input":"2023-11-30T09:39:30.393138Z","iopub.status.busy":"2023-11-30T09:39:30.392624Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Percentage list of missing data\n\nIf there are many features in the set and rendering takes a long time, you can list the proportion of missing records for each feature.\n\nThe reactivity_0027 attribute is missing 5% of the values, and reactivity_0136 is missing 100%.\n\nThis list is a useful summary that can perfectly complement a heat map visualization.","metadata":{}},{"cell_type":"code","source":"for col in df.columns:\n    pct_missing = np.mean(df[col].isnull())\n    print('{} - {}%'.format(col, round(pct_missing*100)))","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Histogram of missing data\n\nAnother good visualization technique for sets with a large number of features is to plot a histogram of the number of missing values in a record.\n\nFrom this it is clear that from? thousand records more than ? thousand rows do not have a single missing value, and about ? thousand - only one. Such rows can be used as “reference” lines to test various hypotheses for data addition.","metadata":{}},{"cell_type":"code","source":"if APP.show_hist_skip_val:\n   # first create an indicator for features with missing data\n    hist_cols=[]\n    dict = {}\n    for col in df.columns:\n        missing = df[col].isnull()\n        num_missing = np.sum(missing)\n        if num_missing > 0:  \n            hist_cols.append(col)\n            dict[col] = missing\n\n    # then we build a histogram based on the indicator\n    df_hist = pd.DataFrame(dict)\n    df['num_missing'] = df_hist[hist_cols].sum(axis=1)\n    df['num_missing'].value_counts().reset_index(drop=True).sort_index().plot.bar(x='index', y='num_missing')\n    df.drop(['num_missing'], axis=1,inplace=True)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## What to do with missing values?\n\nThere are no general solutions to the problem of missing data. For each specific set, one has to look for the most suitable methods or combinations thereof. Let's look at common techniques. They will help in simple situations, but most likely you will have to be creative and look for non-trivial solutions, for example, modeling gaps.","metadata":{}},{"cell_type":"markdown","source":"## Dropping features\n\nFeature rejection can only be used for uninformative features.\n\nIn the percentage list built earlier, we saw that some features have a high percentage of missing values: 97% - 100%. We can completely abandon these signs.","metadata":{}},{"cell_type":"code","source":"if CFG.drop_col_activity:\n    # Removing columns 0001 - 0026 and 0136 - 0206 or 0127-0206\n    col_end = CFG.drop_col_end_max if CFG.drop_col_act_max else CFG.drop_col_end_min\n    react_cols = df.filter(like='reactivity_0').columns.to_list()\n    drop_cols = []\n    drop_cols.extend(react_cols[0:CFG.skip_seq_beg])\n    drop_cols.extend(react_cols[col_end:])\n\n    react_error_cols = df.filter(like='reactivity_error_0').columns.to_list()\n    drop_cols.extend(react_error_cols[0:CFG.skip_seq_beg])\n    drop_cols.extend(react_error_cols[col_end:])\n\n    df.drop(drop_cols, axis=1,inplace=True)\n\n    react_cols = df.filter(like='reactivity_0').columns.to_list()\n    print(react_cols)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"if APP.show_heatmap:\n    sns.heatmap(df[react_cols].isnull(), cmap=sns.color_palette(APP.colours))","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Dropping records\n\nThe first technique in statistics is called listwise deletion and simply discards records containing missing values. This solution is only suitable if the missing data is not informative.\n\nOther criteria can be used for discarding.","metadata":{}},{"cell_type":"code","source":"react_cols = df.filter(like='reactivity_0').columns.to_list()\nfor col in react_cols:\n    pct_missing = np.mean(df[col].isnull())\n    print('{} - {}%'.format(col, round(pct_missing*100)))","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"if CFG.remove_react_nan:\n    hist_cols=[]\n    dict = {}\n    for col in react_cols:\n        missing = df[col].isnull()\n        num_missing = np.sum(missing)  # Number of null values in the column\n        if num_missing > 0:  \n            hist_cols.append(col)\n            dict[col] = missing\n\n    # then we build a histogram based on the indicator\n    df_hist = pd.DataFrame(dict)\n    df['num_miss_react'] = df_hist[hist_cols].sum(axis=1)\n\n    # discard lines with a large number of gaps\n    ind_missing = df[df['num_miss_react'] >= 10].index\n    df = df.drop(ind_missing, axis=0)\n    df.drop(['num_miss_react'], axis=1,inplace=True)\n    print(df.shape)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"react_cols = df.filter(like='reactivity_0').columns.to_list()\nif APP.show_heatmap:\n    sns.heatmap(df[react_cols].isnull(), cmap=sns.color_palette(APP.colours))\nprint(len(react_cols))","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"if CFG.remove_error_nan:\n    react_error_cols = df.filter(like='reactivity_error_0').columns.to_list()\n    hist_cols=[]\n    dict = {}\n    for col in react_error_cols:\n        missing = df[col].isnull()  # Number of null values in the column\n        num_missing = np.sum(missing)\n        if num_missing > 0:  \n            hist_cols.append(col)\n            dict[col] = missing\n\n    # then we build a histogram based on the indicator\n    df_hist = pd.DataFrame(dict)\n    df['num_miss_error_react'] = df_hist[hist_cols].sum(axis=1)\n\n    # discard lines with a large number of gaps\n    ind_missing = df[df['num_miss_error_react'] >= 10].index\n    df = df.drop(ind_missing, axis=0)\n    df.drop(['num_miss_error_react'], axis=1,inplace=True)\n    print(df.shape)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"react_error_cols = df.filter(like='reactivity_error_0').columns.to_list()\nif APP.show_heatmap:\n    sns.heatmap(df[react_error_cols].isnull(), cmap=sns.color_palette(APP.colours))","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(df[\"sequence\"].apply(lambda str: str[0:26] ).nunique())\nprint(df[\"sequence\"].apply(lambda str: str[:-22:-1] ).nunique())","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"if APP.show_hist_seq_len:\n    plt.hist(df[\"sequence\"].apply(len), bins=100)\n    plt.xlabel('Sequence length')\n    plt.ylabel('Quantity')\n    plt.title('Sequence length distribution')\n    plt.show()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"if CFG.remove_seq_not_177:\n    # discard lines with a lot of gaps\n    ind_missing = df[df['sequence'].apply(len) != 177].index\n    df = df.drop(ind_missing, axis=0)\n    print(df.shape)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"if APP.show_hist_seq_len:\n    plt.hist(df[\"sequence\"].apply(len), bins=100)\n    plt.xlabel('Sequence length')\n    plt.ylabel('Quantity')\n    plt.title('Sequence length distribution')\n    plt.show()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"if CFG.cut_sequence:\n    df[\"sequence\"] = df[\"sequence\"].apply(lambda str: str[CFG.skip_seq_beg:len(str)-CFG.skip_seq_end])\n    df[\"structure\"] = df[\"structure\"].apply(lambda str: str[CFG.skip_seq_beg:len(str)-CFG.skip_seq_end])\n    df[\"sequence_ext\"] = df[\"sequence_ext\"].apply(lambda str: str[CFG.skip_seq_beg:len(str)-CFG.skip_seq_end])\n    print(len(df.iloc[0].sequence),len(df.iloc[0].structure),len(df.iloc[0].sequence_ext))","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"if APP.show_hist_seq_len:\n    plt.hist(df[\"sequence\"].apply(len), bins=100)\n    plt.xlabel('Sequence length')\n    plt.ylabel('Quantity')\n    plt.title('Sequence length distribution')\n    plt.show()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Descriptive statistics\n\nDeviations in numerical features may be too clear to not be visualized by a boxplot. Instead, you can analyze their descriptive statistics.","metadata":{}},{"cell_type":"code","source":"if APP.show_hist_react_error:\n    df['reactivity_error_0027'].describe()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Column chart\n\nFor categorical features, you can build a bar chart to visualize data about the categories and their distribution.\n\nFor example, the distribution of the ecology attribute is quite uniform and acceptable. But if there is a category with only one value \"other\", then it will be an outlier.","metadata":{}},{"cell_type":"code","source":"if APP.show_hist_react_error:\n    df['reactivity_error_0027'].value_counts().plot.bar()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Merge reactivity columns into one column","metadata":{}},{"cell_type":"code","source":"if CFG.merge_react_in_one_col:\n    react_cols = df.filter(like='reactivity_0').columns.to_list()\n    react_error_cols = df.filter(like='reactivity_error_0').columns.to_list()\n    df['reactivity'] = df[react_cols].values.tolist()\n    df['react_error'] = df[react_error_cols].values.tolist()\n    \n    drop_cols = []\n    drop_cols.extend(react_cols)\n    drop_cols.extend(react_error_cols)\n    df.drop(drop_cols, axis=1,inplace=True)\n    print(df.shape)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df.head(1)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Split datasets","metadata":{}},{"cell_type":"code","source":"df_2A3 = df.loc[df.experiment_type == '2A3_MaP']  # Shape=(821840, 420)\ndf_DMS = df.loc[df.experiment_type == 'DMS_MaP']  # Shape=(821840, 420)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## What to do with duplicates?\n\nObviously, we don’t need duplicate records, which means they need to be excluded from the set.","metadata":{}},{"cell_type":"code","source":"if CFG.remove_duplicates:\n    df_2A3.sort_values(by='signal_to_noise',ascending=False,inplace=True)\n    df_DMS.sort_values(by='signal_to_noise',ascending=False,inplace=True)\n    df_2A3.drop_duplicates(subset=['sequence'],keep='first',ignore_index=True,inplace=True)\n    df_DMS.drop_duplicates(subset=['sequence'],keep='first',ignore_index=True,inplace=True)\n    df_2A3.sort_values(by='sequence_id',ascending=False,inplace=True)\n    df_DMS.sort_values(by='sequence_id',ascending=False,inplace=True)\n    df_2A3.reset_index(drop=True,inplace=True)\n    df_DMS.reset_index(drop=True,inplace=True)\n    print(df_2A3.loc[df_2A3.sequence_id == '3e5377bfc2c2'].index)\n    print(df_DMS.loc[df_DMS.sequence_id == '3e5377bfc2c2'].index)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Merge sets\n\nColumns:\n- General: sequence_id sequence structure sequence_ext experiment_type\n- Unnecessary: dataset_name\n- Individual: reads signal_to_noise SN_filter reactivity react_error","metadata":{}},{"cell_type":"code","source":"if CFG.merge_by_columns:\n    # Merge datasets by columns\n    df_2A3.drop(['dataset_name','experiment_type'], axis=1,inplace=True)\n    df_DMS.drop(['dataset_name','experiment_type','sequence','structure','sequence_ext'], axis=1,inplace=True)\n    rename_2A3_cols = {}\n    rename_DMS_cols = {}\n    for col in df_DMS.columns:\n        rename_2A3_cols[col] = '2A3_'+col\n        rename_DMS_cols[col] = 'DMS_'+col\n    df_2A3.rename(columns=rename_2A3_cols,inplace=True)    \n    df_DMS.rename(columns=rename_DMS_cols,inplace=True)\n    df_2A3['sequence_id'] = df_2A3['2A3_sequence_id']\n\n    # Merge datasets\n    df_2A3.set_index(keys='2A3_sequence_id', inplace = True)\n    df_DMS.set_index(keys='DMS_sequence_id', inplace = True)\n    df = pd.concat([df_2A3, df_DMS], axis=1)\nelse:\n    # Merge datasets by row\n    df = pd.concat([df_2A3, df_DMS], axis=0)\n    df.sort_index(inplace=True)\ndf.reset_index(drop=True,inplace=True)\nprint(df.loc[df.sequence_id == '3e5377bfc2c2'].index)\nprint(df.shape)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df.head(1)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Remove missing values","metadata":{}},{"cell_type":"code","source":"if CFG.dropna:\n    for col in df.columns:\n        pct_missing = np.mean(df[col].isnull())\n        print('{} - {}%'.format(col, round(pct_missing*100)))\n    \n    df.dropna(inplace=True) \n    df.reset_index(drop=True,inplace=True)\n    print(df.shape)\n    \n    for col in df.columns:\n        pct_missing = np.mean(df[col].isnull())\n        print('{} - {}%'.format(col, round(pct_missing*100)))","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Save data","metadata":{}},{"cell_type":"code","source":"df.reset_index(drop=True,inplace=True)\ndf.to_parquet(APP.file_trg)\n\ndf=pd.read_parquet(APP.file_trg)\nprint(df.shape)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df.head(-5)","metadata":{"trusted":true},"execution_count":null,"outputs":[]}]}