{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"### - Add volcano plot\n### - Add comparison `pseudobulk` vs `ground truth pseudobulk`\n\n\n### some code from: https://github.com/openproblems-bio/neurips-2023-scripts/blob/main/compute_de.ipynb","metadata":{"execution":{"iopub.status.busy":"2023-10-09T09:23:11.908667Z","iopub.execute_input":"2023-10-09T09:23:11.909766Z","iopub.status.idle":"2023-10-09T09:23:11.914591Z","shell.execute_reply.started":"2023-10-09T09:23:11.909729Z","shell.execute_reply":"2023-10-09T09:23:11.913422Z"}}},{"cell_type":"code","source":"! pip install scanpy","metadata":{"_kg_hide-output":true,"execution":{"iopub.status.busy":"2023-10-20T09:13:08.257024Z","iopub.execute_input":"2023-10-20T09:13:08.258190Z","iopub.status.idle":"2023-10-20T09:13:19.647327Z","shell.execute_reply.started":"2023-10-20T09:13:08.258134Z","shell.execute_reply":"2023-10-20T09:13:19.645810Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import os, gc\n\nfrom glob import glob\nfrom tqdm.auto import tqdm\nfrom itertools import product\n\n\nimport numpy as np\nimport pandas as pd\n\nimport matplotlib.pyplot as plt\nimport seaborn as sns\n\nimport anndata as ad\nimport scanpy as sc","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2023-10-20T09:13:19.649552Z","iopub.execute_input":"2023-10-20T09:13:19.649907Z","iopub.status.idle":"2023-10-20T09:13:19.657748Z","shell.execute_reply.started":"2023-10-20T09:13:19.649876Z","shell.execute_reply":"2023-10-20T09:13:19.656329Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"cell_rename_dict = {\n    'NK cells' : 'NK',\n    'T cells CD4+' : 'T4',\n    'T cells CD8+' : 'T8',\n    'T regulatory cells' : 'Treg',\n    'B cells' : 'B',\n    'Myeloid cells' : 'Myeloid',\n}\n\n\nde_pert_cols = [\n    'sm_name',\n    'sm_lincs_id',\n    'SMILES',\n    'dose_uM',\n    'timepoint_hr',\n    'cell_type',\n]","metadata":{"execution":{"iopub.status.busy":"2023-10-20T09:13:19.659330Z","iopub.execute_input":"2023-10-20T09:13:19.659650Z","iopub.status.idle":"2023-10-20T09:13:19.671202Z","shell.execute_reply.started":"2023-10-20T09:13:19.659624Z","shell.execute_reply":"2023-10-20T09:13:19.670303Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"pert_de_dfs = []\n\nfor rename_cell_type in tqdm(cell_rename_dict.values(), position=1, desc='cell_type', leave=True):\n    cell_type_working_dir = 'cell_type_{}'.format(rename_cell_type)\n\n    cell_type_bulk_adata = sc.read_h5ad(f'/kaggle/input/op2-02-pseudobulk-v2/{cell_type_working_dir}/input.h5ad')\n    for file in tqdm(glob(f'/kaggle/input/op2-03-limma/v2/{cell_type_working_dir}/pert_*_contrast_result.csv'), position=2, desc='pert', leave=False):\n        pert_de_df = pd.read_csv(file)\n        pert_de_df = pert_de_df.rename({pert_de_df.columns[0]: 'gene'}, axis=1)\n\n        pert = os.path.basename(file).split('_')[1]\n        pert_de_df['Rpert'] = pert\n\n        pert_obs = cell_type_bulk_adata.obs[cell_type_bulk_adata.obs['Rpert'].eq(pert)]\n        for col in de_pert_cols:\n            pert_de_df[col] = pert_obs[col].unique()[0]\n\n        pert_de_dfs.append(pert_de_df)\n\nde_df_v2 = pd.concat(pert_de_dfs)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"pert_de_dfs = []\n\nfor rename_cell_type in tqdm(cell_rename_dict.values(), position=1, desc='cell_type', leave=True):\n    cell_type_working_dir = 'cell_type_{}'.format(rename_cell_type)\n\n    cell_type_bulk_adata = sc.read_h5ad(f'/kaggle/input/op2-02-pseudobulk-v3/{cell_type_working_dir}/input.h5ad')\n    for file in tqdm(glob(f'/kaggle/input/op2-03-limma/v3/{cell_type_working_dir}/pert_*_contrast_result.csv'), position=2, desc='pert', leave=False):\n        pert_de_df = pd.read_csv(file)\n        pert_de_df = pert_de_df.rename({pert_de_df.columns[0]: 'gene'}, axis=1)\n\n        pert = os.path.basename(file).split('_')[1]\n        pert_de_df['Rpert'] = pert\n\n        pert_obs = cell_type_bulk_adata.obs[cell_type_bulk_adata.obs['Rpert'].eq(pert)]\n        for col in de_pert_cols:\n            pert_de_df[col] = pert_obs[col].unique()[0]\n\n        pert_de_dfs.append(pert_de_df)\n\nde_df_v3 = pd.concat(pert_de_dfs)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"---","metadata":{}},{"cell_type":"markdown","source":"# Volcano plot","metadata":{}},{"cell_type":"code","source":"de_df_v2.head()","metadata":{"execution":{"iopub.status.busy":"2023-10-20T09:22:50.531533Z","iopub.execute_input":"2023-10-20T09:22:50.531935Z","iopub.status.idle":"2023-10-20T09:22:50.552262Z","shell.execute_reply.started":"2023-10-20T09:22:50.531903Z","shell.execute_reply":"2023-10-20T09:22:50.550855Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig_df = de_df_v2.sample(int(1e4))\nfig_df['log10_p'] = fig_df['P.Value'].apply(np.log10)\n\nsns.scatterplot(\n    data=fig_df,\n    x='logFC', y='log10_p',\n    hue='cell_type',\n    edgecolor='k', lw=0.1, alpha=0.5\n)\n\nplt.xlim(-4, 4)\nplt.ylim(0.2, -5.2)\n\nplt.show()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"---","metadata":{}},{"cell_type":"markdown","source":"### https://www.kaggle.com/competitions/open-problems-single-cell-perturbations/discussion/448556\n\n> What is the functional connection between DE(A) and p-value(A)?\n>\n> Is it possible for Gene A and Gene B to have identical p-values while having drastically different differential expressions like +2 vs +10?","metadata":{}},{"cell_type":"code","source":"de_df_v2['logFC_round4'] = de_df_v2['logFC'].apply(lambda x: round(x, 4))\nde_df_v2['logFC_round4_count'] = de_df_v2.groupby(['sm_name', 'cell_type','logFC_round4'])['logFC'].transform('count')\n\nde_df_v2['log10_p'] = de_df_v2['P.Value'].apply(np.log10)\nde_df_v2['log10_p_round4'] = de_df_v2['log10_p'].apply(lambda x: round(x, 4))\nde_df_v2['log10_p_round4_count'] = de_df_v2.groupby(['sm_name', 'cell_type','log10_p_round4'])['log10_p'].transform('count')","metadata":{"execution":{"iopub.status.busy":"2023-10-20T09:45:03.132129Z","iopub.execute_input":"2023-10-20T09:45:03.132980Z","iopub.status.idle":"2023-10-20T09:45:32.245647Z","shell.execute_reply.started":"2023-10-20T09:45:03.132931Z","shell.execute_reply":"2023-10-20T09:45:32.244686Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"logFC_round4 = de_df_v2.query('10>logFC_round4_count>=3')['logFC_round4'].unique()\n# pd.Series(logFC_round4[(2<logFC_round4) & (logFC_round4<3)]).sort_values()","metadata":{"execution":{"iopub.status.busy":"2023-10-20T09:46:23.734633Z","iopub.execute_input":"2023-10-20T09:46:23.735601Z","iopub.status.idle":"2023-10-20T09:46:24.864255Z","shell.execute_reply.started":"2023-10-20T09:46:23.735559Z","shell.execute_reply":"2023-10-20T09:46:24.863448Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"de_df_v2.query('logFC_round4_count>=3 and logFC_round4==2.0002')[['gene', 'cell_type', 'sm_lincs_id', 'logFC', 'P.Value', 'log10_p']]","metadata":{"execution":{"iopub.status.busy":"2023-10-20T09:48:40.894190Z","iopub.execute_input":"2023-10-20T09:48:40.894652Z","iopub.status.idle":"2023-10-20T09:48:41.030504Z","shell.execute_reply.started":"2023-10-20T09:48:40.894609Z","shell.execute_reply":"2023-10-20T09:48:41.029206Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"de_df_v2.query('logFC_round4_count>=3 and logFC_round4==5.0437')[['gene', 'cell_type', 'sm_lincs_id', 'logFC', 'P.Value', 'log10_p']]","metadata":{"execution":{"iopub.status.busy":"2023-10-20T10:08:18.144210Z","iopub.execute_input":"2023-10-20T10:08:18.144780Z","iopub.status.idle":"2023-10-20T10:08:18.472848Z","shell.execute_reply.started":"2023-10-20T10:08:18.144743Z","shell.execute_reply":"2023-10-20T10:08:18.471882Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"logp_round4 = de_df_v2.query('5>log10_p_round4_count>=3')['log10_p_round4'].unique()\n# pd.Series(logp_round4[(-5<logp_round4) & (logp_round4<-4)]).sort_values()","metadata":{"execution":{"iopub.status.busy":"2023-10-20T09:54:23.098307Z","iopub.execute_input":"2023-10-20T09:54:23.098668Z","iopub.status.idle":"2023-10-20T09:54:24.081687Z","shell.execute_reply.started":"2023-10-20T09:54:23.098642Z","shell.execute_reply":"2023-10-20T09:54:24.080481Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"de_df_v2.query('log10_p_round4_count>=3 and log10_p_round4==-4.0010')[['gene', 'cell_type', 'sm_lincs_id', 'logFC', 'P.Value', 'log10_p']]","metadata":{"execution":{"iopub.status.busy":"2023-10-20T09:49:26.981704Z","iopub.execute_input":"2023-10-20T09:49:26.982958Z","iopub.status.idle":"2023-10-20T09:49:27.118596Z","shell.execute_reply.started":"2023-10-20T09:49:26.982909Z","shell.execute_reply":"2023-10-20T09:49:27.117317Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"de_df_v2.query('log10_p_round4_count>=3 and log10_p_round4==-9.7344')[['gene', 'cell_type', 'sm_lincs_id', 'logFC', 'P.Value', 'log10_p']]","metadata":{"execution":{"iopub.status.busy":"2023-10-20T09:49:52.253782Z","iopub.execute_input":"2023-10-20T09:49:52.254267Z","iopub.status.idle":"2023-10-20T09:49:52.438491Z","shell.execute_reply.started":"2023-10-20T09:49:52.254215Z","shell.execute_reply":"2023-10-20T09:49:52.437322Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"de_df_v2.query('log10_p_round4_count>=3 and log10_p_round4==-18.9565')[['gene', 'cell_type', 'sm_lincs_id', 'logFC', 'P.Value', 'log10_p']]","metadata":{"execution":{"iopub.status.busy":"2023-10-20T09:52:04.408412Z","iopub.execute_input":"2023-10-20T09:52:04.409495Z","iopub.status.idle":"2023-10-20T09:52:04.593780Z","shell.execute_reply.started":"2023-10-20T09:52:04.409447Z","shell.execute_reply":"2023-10-20T09:52:04.592986Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"de_df_v2.query('log10_p_round4_count>=3 and log10_p_round4==-30.9375')[['gene', 'cell_type', 'sm_lincs_id', 'logFC', 'P.Value', 'log10_p']]","metadata":{"execution":{"iopub.status.busy":"2023-10-20T09:54:30.223495Z","iopub.execute_input":"2023-10-20T09:54:30.224184Z","iopub.status.idle":"2023-10-20T09:54:30.349322Z","shell.execute_reply.started":"2023-10-20T09:54:30.224151Z","shell.execute_reply":"2023-10-20T09:54:30.348520Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"---\n---\n---","metadata":{}},{"cell_type":"markdown","source":"## Converting DataFrame to Anndata","metadata":{}},{"cell_type":"code","source":"def convert_de_df_to_anndata(de_df, pert_cols, de_sig_cutoff):\n    de_df = de_df.copy()\n    zero_pval_selection = de_df['P.Value'].eq(0)\n    de_df.loc[zero_pval_selection, 'P.Value'] = np.finfo(np.float64).eps\n\n    de_df['sign_log10_pval'] = np.sign(de_df['logFC']) * -np.log10(de_df['P.Value'])\n    de_df['is_de'] = de_df['P.Value'].lt(de_sig_cutoff)\n    de_df['is_de_adj'] = de_df['adj.P.Val'].lt(de_sig_cutoff)\n\n    de_feature_dfs = {}\n    for feature in tqdm(['is_de', 'is_de_adj', 'sign_log10_pval', 'logFC', 'P.Value', 'adj.P.Val']):\n        df = de_df.reset_index().pivot_table(\n            index=['gene'], \n            columns=pert_cols,\n            values=[feature],\n            dropna=True,\n        )\n        de_feature_dfs[feature] = df\n\n    # de_adata = ad.AnnData(de_feature_dfs['sign_log10_pval'].T, dtype=np.float64)\n    # de_adata.obs = de_adata.obs.reset_index()\n    # de_adata.obs = de_adata.obs.drop(columns=['level_0'])\n    # de_adata.obs.index = de_adata.obs.index.astype('string')\n\n    tdf = de_feature_dfs['sign_log10_pval'].T.reset_index().drop(columns=['level_0'])\n    tdf = tdf.astype({'dose_uM':str, 'timepoint_hr':str}).set_index(['sm_name', 'sm_lincs_id', 'SMILES', 'dose_uM', 'timepoint_hr', 'cell_type'])\n    de_adata = ad.AnnData(\n        tdf.values,\n        obs=pd.DataFrame(tdf.index.to_frame().values, columns=tdf.index.names),\n        var=tdf.columns.to_frame(),\n        dtype=np.float64\n    )\n\n    de_adata.layers['is_de'] = de_feature_dfs['is_de'].to_numpy().T\n    de_adata.layers['is_de_adj'] = de_feature_dfs['is_de_adj'].to_numpy().T\n    de_adata.layers['logFC'] = de_feature_dfs['logFC'].to_numpy().T\n    de_adata.layers['P.Value'] = de_feature_dfs['P.Value'].to_numpy().T\n    de_adata.layers['adj.P.Val'] = de_feature_dfs['adj.P.Val'].to_numpy().T\n    de_adata.obs.index = de_adata.obs.index.astype('str')\n    \n    return de_adata","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"de_adata_v2 = convert_de_df_to_anndata(de_df_v2, de_pert_cols, 0.05)\nde_adata_v3 = convert_de_df_to_anndata(de_df_v3, de_pert_cols, 0.05)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sorting_index = de_adata_v2.obs.astype(str).sort_values(['sm_name', 'cell_type']).index\nde_adata_v2 = de_adata_v2[sorting_index].copy()\n\nsorting_index = de_adata_v3.obs.astype(str).sort_values(['sm_name', 'cell_type']).index\nde_adata_v3 = de_adata_v3[sorting_index].copy()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"de_adata_v2.write_h5ad('de_adata_v2.h5ad')\nde_adata_v3.write_h5ad('de_adata_v3.h5ad')","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"---","metadata":{}},{"cell_type":"markdown","source":"# Volcano plot","metadata":{}},{"cell_type":"code","source":"genes = np.random.choice(de_adata_v2.var.index, 1000, replace=False)\n\nfig, axs = plt.subplots(ncols=2, nrows=2, figsize=(10, 10))\naxs = axs.flatten()\n\nfor i, sm_lincs_id in enumerate(np.random.choice(de_adata_v2.obs['sm_lincs_id'].unique(), 4, replace=False)):\n\n    idx = de_adata_v2.obs.query('sm_lincs_id==@sm_lincs_id').index\n\n    fig_df_x = (\n        pd.DataFrame(\n            de_adata_v2[idx, genes].layers['logFC'].copy(),\n            index=pd.MultiIndex.from_frame(de_adata_v2.obs.loc[idx, ['sm_lincs_id', 'cell_type']]),\n            columns=genes,\n        )\n        .reset_index()\n        .melt(\n            id_vars=['sm_lincs_id', 'cell_type'],\n            var_name='gene',\n            value_name='logFC'\n        )\n    )\n\n    fig_df_y = (\n        pd.DataFrame(\n            de_adata_v2[idx, genes].layers['P.Value'].copy(),\n            index=pd.MultiIndex.from_frame(de_adata_v2.obs.loc[idx, ['sm_lincs_id', 'cell_type']]),\n            columns=genes,\n        )\n        .reset_index()\n        .melt(\n            id_vars=['sm_lincs_id', 'cell_type'],\n            var_name='gene',\n            value_name='log10_p'\n        )\n    )\n\n    fig_df_y['log10_p'] = fig_df_y['log10_p'].apply(np.log10)\n    \n    sns.scatterplot(\n        data=pd.merge(fig_df_x, fig_df_y, on=['sm_lincs_id', 'cell_type', 'gene']),\n        x='logFC', y='log10_p',\n        hue='cell_type',\n        edgecolor='k', lw=0.1, alpha=0.5,\n        ax=axs[i]\n    )\n\n    axs[i].set_title(sm_lincs_id)\n    axs[i].set_xlim(-4, 4)\n    axs[i].set_ylim(0.2, -5.2)\n\nplt.show()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"---","metadata":{}},{"cell_type":"markdown","source":"## Checking the results","metadata":{}},{"cell_type":"code","source":"data_dir = '/kaggle/input/open-problems-single-cell-perturbations'\n\nkaggle_train_de_df = pd.read_parquet(os.path.join(data_dir, 'de_train.parquet'))\nkaggle_train_de_df = kaggle_train_de_df.set_index(list(kaggle_train_de_df.columns[:5]))\n\n# kaggle_train_de_adata = ad.AnnData(kaggle_train_de_df)\n# kaggle_train_de_adata.obs = kaggle_train_de_adata.obs.reset_index()\n# kaggle_train_de_adata.obs.index = kaggle_train_de_adata.obs.index.astype('str')\n\nkaggle_train_de_adata = ad.AnnData(\n    kaggle_train_de_df.values,\n    obs=pd.DataFrame(kaggle_train_de_df.index.to_frame().values, columns=kaggle_train_de_df.index.names),\n    var=kaggle_train_de_df.columns.rename('gene').to_frame(),\n)\n\nsorting_index = kaggle_train_de_adata.obs.astype(str).sort_values(['sm_name', 'cell_type']).index\nkaggle_train_de_adata = kaggle_train_de_adata[sorting_index].copy()","metadata":{"execution":{"iopub.status.busy":"2023-10-19T08:39:43.226310Z","iopub.execute_input":"2023-10-19T08:39:43.227912Z","iopub.status.idle":"2023-10-19T08:39:45.159299Z","shell.execute_reply.started":"2023-10-19T08:39:43.227865Z","shell.execute_reply":"2023-10-19T08:39:45.157623Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"---\n---\n---","metadata":{"execution":{"iopub.status.busy":"2023-10-13T10:22:43.876030Z","iopub.execute_input":"2023-10-13T10:22:43.876489Z","iopub.status.idle":"2023-10-13T10:22:44.134588Z","shell.execute_reply.started":"2023-10-13T10:22:43.876448Z","shell.execute_reply":"2023-10-13T10:22:44.133430Z"}}},{"cell_type":"code","source":"def mrrmse(y_pred: pd.DataFrame, y_true: pd.DataFrame):\n    return ((y_pred - y_true)**2).mean(axis=1).apply(np.sqrt).mean()","metadata":{"execution":{"iopub.status.busy":"2023-10-19T08:45:31.010546Z","iopub.execute_input":"2023-10-19T08:45:31.011077Z","iopub.status.idle":"2023-10-19T08:45:31.018329Z","shell.execute_reply.started":"2023-10-19T08:45:31.011045Z","shell.execute_reply":"2023-10-19T08:45:31.016947Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# kaggle `de_train` vs `ground truth pseudobulk` (v2)","metadata":{}},{"cell_type":"code","source":"fig, axs = plt.subplots(figsize=(14, 6), ncols=5, nrows=2, gridspec_kw={'hspace':0.5})\naxs = axs.flatten()\n\nfor i, idx in enumerate(np.random.choice(np.arange(de_adata_v2.shape[0]), 10)):\n    sns.scatterplot(\n        x=kaggle_train_de_adata.X[idx],\n        y=de_adata_v2[:, kaggle_train_de_adata.var.index].X[idx],\n        ax=axs[i],\n    )\n    cell_type, sm_lincs_id = de_adata_v2.obs.loc[str(idx), ['cell_type', 'sm_lincs_id']]\n    axs[i].set_title(f'{cell_type} / {sm_lincs_id}')\n\n\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-10-19T08:45:58.978656Z","iopub.execute_input":"2023-10-19T08:45:58.979204Z","iopub.status.idle":"2023-10-19T08:46:01.438419Z","shell.execute_reply.started":"2023-10-19T08:45:58.979163Z","shell.execute_reply":"2023-10-19T08:46:01.437203Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"mrrmse(\n    pd.DataFrame(de_adata_v2.X),\n    pd.DataFrame(kaggle_train_de_adata.X),\n)","metadata":{"execution":{"iopub.status.busy":"2023-10-19T08:46:01.440501Z","iopub.execute_input":"2023-10-19T08:46:01.441112Z","iopub.status.idle":"2023-10-19T08:46:01.658394Z","shell.execute_reply.started":"2023-10-19T08:46:01.441077Z","shell.execute_reply":"2023-10-19T08:46:01.657137Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# kaggle `de_train` vs `pseudobulk` (v3)","metadata":{}},{"cell_type":"code","source":"fig, axs = plt.subplots(figsize=(14, 6), ncols=5, nrows=2, gridspec_kw={'hspace':0.5})\naxs = axs.flatten()\n\nfor i, idx in enumerate(np.random.choice(np.arange(de_adata_v3.shape[0]), 10)):\n    sns.scatterplot(\n        x=kaggle_train_de_adata.X[idx],\n        y=de_adata_v3[:, kaggle_train_de_adata.var.index].X[idx],\n        ax=axs[i],\n    )\n    cell_type, sm_lincs_id = de_adata_v3.obs.loc[str(idx), ['cell_type', 'sm_lincs_id']]\n    axs[i].set_title(f'{cell_type} / {sm_lincs_id}')\n\n\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-10-19T08:46:01.660368Z","iopub.execute_input":"2023-10-19T08:46:01.661229Z","iopub.status.idle":"2023-10-19T08:46:04.253058Z","shell.execute_reply.started":"2023-10-19T08:46:01.661183Z","shell.execute_reply":"2023-10-19T08:46:04.251784Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"mrrmse(\n    pd.DataFrame(de_adata_v3.X),\n    pd.DataFrame(kaggle_train_de_adata.X),\n)","metadata":{"execution":{"iopub.status.busy":"2023-10-19T08:46:04.256153Z","iopub.execute_input":"2023-10-19T08:46:04.257232Z","iopub.status.idle":"2023-10-19T08:46:04.464368Z","shell.execute_reply.started":"2023-10-19T08:46:04.257186Z","shell.execute_reply":"2023-10-19T08:46:04.463091Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# `ground truth pseudobulk` (v2) vs `pseudobulk` (v3)","metadata":{}},{"cell_type":"code","source":"fig, axs = plt.subplots(figsize=(14, 6), ncols=5, nrows=2, gridspec_kw={'hspace':0.5})\naxs = axs.flatten()\n\nfor i, idx in enumerate(np.random.choice(np.arange(de_adata_v2.shape[0]), 10)):\n    sns.scatterplot(\n        x=de_adata_v2.X[idx],\n        y=de_adata_v3[:, de_adata_v2.var.index].X[idx],\n        ax=axs[i],\n    )\n    cell_type, sm_lincs_id = de_adata_v2.obs.loc[str(idx), ['cell_type', 'sm_lincs_id']]\n    axs[i].set_title(f'{cell_type} / {sm_lincs_id}')\n\n\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-10-19T08:46:04.466299Z","iopub.execute_input":"2023-10-19T08:46:04.467574Z","iopub.status.idle":"2023-10-19T08:46:07.746403Z","shell.execute_reply.started":"2023-10-19T08:46:04.467528Z","shell.execute_reply":"2023-10-19T08:46:07.745105Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"mrrmse(\n    pd.DataFrame(de_adata_v2.X),\n    pd.DataFrame(de_adata_v3.X),\n)","metadata":{"execution":{"iopub.status.busy":"2023-10-19T08:46:07.748321Z","iopub.execute_input":"2023-10-19T08:46:07.748796Z","iopub.status.idle":"2023-10-19T08:46:07.952261Z","shell.execute_reply.started":"2023-10-19T08:46:07.748752Z","shell.execute_reply":"2023-10-19T08:46:07.951015Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"---\n---\n---","metadata":{}},{"cell_type":"markdown","source":"# Visualize the results from https://www.kaggle.com/code/awater1223/op2-05-preprocessing-compare/notebook","metadata":{}},{"cell_type":"code","source":"# preprocess_names = [\n#     'logCP10K',\n# #     'scale',\n# #     'median',\n#     'mean',\n#     'size',\n#     'none',\n# ]\n\n# # reduction_names = ['none', 'pca128', 'pca512', 'svd128', 'svd512']\n# # reduction_names = ['none', 'svd128', 'svd512']\n# reduction_names = ['none']\n\n# for exp_id, (preprocess_name, reduction_name) in enumerate(tqdm(list(product(preprocess_names, reduction_names)), position=1, desc='exp_id')):\n#     print ('[DEBUG]', exp_id, preprocess_name, reduction_name)\n#     pert_de_dfs = []\n\n#     for rename_cell_type in tqdm(cell_rename_dict.values(), position=2, desc='cell_type', leave=False):\n#         cell_type_working_dir = '/kaggle/input/op2-05-preprocessing-compare/exp{}_cell_type_{}'.format(exp_id, rename_cell_type)\n\n#         cell_type_bulk_adata = sc.read_h5ad(f'{cell_type_working_dir}/input.h5ad')\n#         for file in tqdm(glob(f'{cell_type_working_dir}/pert_*_contrast_result.csv'), position=3, desc='pert', leave=False):\n#             pert_de_df = pd.read_csv(file)\n#             pert_de_df = pert_de_df.rename({pert_de_df.columns[0]: 'gene'}, axis=1)\n\n#             pert = os.path.basename(file).split('_')[1]\n#             pert_de_df['Rpert'] = pert\n\n#             pert_obs = cell_type_bulk_adata.obs[cell_type_bulk_adata.obs['Rpert'].eq(pert)]\n#             for col in de_pert_cols:\n#                 pert_de_df[col] = pert_obs[col].unique()[0]\n\n#             pert_de_dfs.append(pert_de_df)\n\n#     gc.collect();\n#     de_df = pd.concat(pert_de_dfs)\n#     de_adata = convert_de_df_to_anndata(de_df, de_pert_cols, 0.05)\n\n#     de_adata.obs.index = de_adata.obs.index.astype('str')\n\n#     sorting_index = de_adata.obs.sort_values(['sm_name', 'cell_type']).index\n#     de_adata = de_adata[sorting_index].copy()\n#     de_adata.write(f'exp{exp_id}_bulk.h5ad')\n\n#     fig, axs = plt.subplots(figsize=(14, 6), ncols=5, nrows=2, gridspec_kw={'hspace':0.5})\n#     axs = axs.flatten()\n\n#     for i, idx in enumerate(np.random.choice(np.arange(de_adata.shape[0]), 10)):\n#         sns.scatterplot(\n#             x=kaggle_train_de_adata.X[idx],\n#             y=de_adata[:, kaggle_train_de_adata.var.index].X[idx],\n#             ax=axs[i],\n#         )\n#         cell_type, sm_lincs_id = de_adata.obs.loc[str(idx), ['cell_type', 'sm_lincs_id']]\n#         axs[i].set_title(f'{cell_type} / {sm_lincs_id}')\n\n#     diff = mrrmse(\n#         pd.DataFrame(de_adata.X),\n#         pd.DataFrame(kaggle_train_de_adata.X),\n#     )\n#     fig.suptitle('{}-{:.3f}'.format(preprocess_name, diff), fontsize=16)\n#     plt.show()\n\n#     gc.collect()\n#     print ()\n","metadata":{"execution":{"iopub.status.busy":"2023-10-19T08:46:15.018492Z","iopub.execute_input":"2023-10-19T08:46:15.018876Z","iopub.status.idle":"2023-10-19T08:46:15.029060Z","shell.execute_reply.started":"2023-10-19T08:46:15.018847Z","shell.execute_reply":"2023-10-19T08:46:15.027766Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}