{"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":"code","source":"! pip install scanpy\n! git clone https://github.com/openproblems-bio/neurips-2023-scripts","metadata":{"_kg_hide-output":true,"execution":{"iopub.status.busy":"2023-10-17T06:18:30.300942Z","iopub.execute_input":"2023-10-17T06:18:30.302055Z","iopub.status.idle":"2023-10-17T06:18:46.831774Z","shell.execute_reply.started":"2023-10-17T06:18:30.302022Z","shell.execute_reply":"2023-10-17T06:18:46.830715Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import os, gc, binascii\n\nfrom glob import glob\nfrom tqdm.auto import tqdm\nfrom itertools import product\n\nimport numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\n\nfrom scipy import sparse\nfrom sklearn.decomposition import PCA, TruncatedSVD\n\nimport matplotlib.pyplot as plt\nimport seaborn as sns\n\nimport anndata as ad\nimport scanpy as sc","metadata":{"execution":{"iopub.status.busy":"2023-10-17T06:18:46.833817Z","iopub.execute_input":"2023-10-17T06:18:46.834142Z","iopub.status.idle":"2023-10-17T06:18:53.968925Z","shell.execute_reply.started":"2023-10-17T06:18:46.834109Z","shell.execute_reply":"2023-10-17T06:18:53.967985Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"---","metadata":{}},{"cell_type":"markdown","source":"## The metadata of psudobulk were NOT match with the ground truth psudobulk\n## (still working on it ...)","metadata":{}},{"cell_type":"code","source":"# fixed_bulk_adata = sc.read_h5ad('/kaggle/input/op2-02-pseudobulk-v2/bulk_adata.h5ad')\n# bulk_adata = sc.read_h5ad('/kaggle/input/op2-02-pseudobulk-v2/raw_bulk_adata.h5ad')\n# bulk_adata = bulk_adata[:, fixed_bulk_adata.var.index]","metadata":{"execution":{"iopub.status.busy":"2023-10-17T07:10:37.784270Z","iopub.execute_input":"2023-10-17T07:10:37.785069Z","iopub.status.idle":"2023-10-17T07:10:37.789957Z","shell.execute_reply.started":"2023-10-17T07:10:37.785031Z","shell.execute_reply":"2023-10-17T07:10:37.788651Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# idx1 = pd.DataFrame(fixed_bulk_adata.X, columns=fixed_bulk_adata.var.index).sort_values(by=fixed_bulk_adata.var.index.tolist()).index\n# idx2 = pd.DataFrame(bulk_adata.X, columns=bulk_adata.var.index).sort_values(by=bulk_adata.var.index.tolist()).index","metadata":{"execution":{"iopub.status.busy":"2023-10-17T06:35:53.309993Z","iopub.execute_input":"2023-10-17T06:35:53.310339Z","iopub.status.idle":"2023-10-17T06:36:04.283673Z","shell.execute_reply.started":"2023-10-17T06:35:53.310313Z","shell.execute_reply":"2023-10-17T06:36:04.282746Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# fixed_bulk_adata = fixed_bulk_adata[fixed_bulk_adata.obs.index[idx1]]\n# bulk_adata = bulk_adata[bulk_adata.obs.index[idx2]]","metadata":{"execution":{"iopub.status.busy":"2023-10-17T07:10:34.849907Z","iopub.execute_input":"2023-10-17T07:10:34.850687Z","iopub.status.idle":"2023-10-17T07:10:34.855119Z","shell.execute_reply.started":"2023-10-17T07:10:34.850638Z","shell.execute_reply":"2023-10-17T07:10:34.853954Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# pd.merge(\n#     fixed_bulk_adata.obs[['cell_type', 'plate_name', 'donor_id', 'sm_lincs_id']].reset_index(),\n#     bulk_adata.obs[['cell_type', 'plate_name', 'donor_id', 'sm_lincs_id']].reset_index(),\n#     on=['cell_type', 'plate_name', 'donor_id', 'sm_lincs_id'],\n#     how='inner'\n# ).shape","metadata":{"execution":{"iopub.status.busy":"2023-10-17T07:10:47.203714Z","iopub.execute_input":"2023-10-17T07:10:47.204544Z","iopub.status.idle":"2023-10-17T07:10:47.222768Z","shell.execute_reply.started":"2023-10-17T07:10:47.204510Z","shell.execute_reply":"2023-10-17T07:10:47.221855Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"---","metadata":{}},{"cell_type":"code","source":"def preprocessing(adata_path, preprocess_name, reduction_name):\n    adata = sc.read(adata_path)\n\n    if preprocess_name == 'none':\n        pass\n    elif preprocess_name == 'logCP10K':\n        sc.pp.normalize_total(adata=adata, target_sum=1e6, copy=False)\n        sc.pp.log1p(adata, copy=False)\n    elif preprocess_name == 'size':\n        sc.pp.normalize_total(adata=adata, target_sum=1e6, copy=False)\n    elif preprocess_name == 'scale':\n        sc.pp.scale(adata, copy=False)\n    elif preprocess_name == 'median':\n        adata.X /= np.median(adata.X.todense(), axis=1).clip(min=1).reshape(-1, 1)\n    elif preprocess_name == 'mean':\n        adata.X /= adata.X.mean(axis=1).clip(min=1).reshape(-1, 1)\n    else:\n        raise NotImplementedError\n    \n    if reduction_name == 'none':\n        pass\n    else:\n        if reduction_name == 'pca128':\n            rd_obj = PCA(n_components=128, svd_solver='arpack', random_state=42)\n        elif reduction_name == 'pca512':\n            rd_obj = PCA(n_components=512, svd_solver='arpack', random_state=42)\n        elif reduction_name == 'svd128':\n            rd_obj = TruncatedSVD(n_components=128, algorithm='arpack', n_iter=7, random_state=42)\n        elif reduction_name == 'svd512':\n            rd_obj = TruncatedSVD(n_components=512, algorithm='arpack', n_iter=7, random_state=42)\n        else:\n            raise NotImplementedError\n        \n        if 'pca' in reduction_name:\n            X_reduced = rd_obj.fit_transform(adata.X.todense().astype('float64'))\n        else:\n            X_reduced = rd_obj.fit_transform(adata.X.astype(float))\n        adata.X = rd_obj.inverse_transform(X_reduced)\n\n    adata.obs['plate_well_cell_type'] = adata.obs['plate_name'].astype('str') \\\n        + '_' + adata.obs['well'].astype('str') \\\n        + '_' + adata.obs['cell_type'].astype('str')\n    adata.obs['plate_well_cell_type'] = adata.obs['plate_well_cell_type'].astype('category')\n    \n    return adata","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Prepare for Limma","metadata":{}},{"cell_type":"code","source":"# The original code: https://github.com/openproblems-bio/neurips-2023-scripts/blob/main/compute_de.ipynb\ndef sum_by(adata: ad.AnnData, col: str='plate_well_cell_type') -> ad.AnnData:\n    \"\"\"\n    Adapted from this forum post: \n    https://discourse.scverse.org/t/group-sum-rows-based-on-jobs-feature/371/4\n    \"\"\"\n    \n    assert pd.api.types.is_categorical_dtype(adata.obs[col])\n\n    # sum `.X` entries for each unique value in `col`\n    cat = adata.obs[col].values\n    indicator = sparse.coo_matrix(\n        (\n            np.broadcast_to(True, adata.n_obs),\n            (cat.codes, np.arange(adata.n_obs))\n        ),\n        shape=(len(cat.categories), adata.n_obs),\n    )\n    sum_adata = ad.AnnData(\n        indicator @ adata.X,\n        var=adata.var,\n        obs=pd.DataFrame(index=cat.categories),\n        dtype=adata.X.dtype,\n    )\n    \n    # copy over `.obs` values that have a one-to-one-mapping with `.obs[col]`\n    obs_cols = adata.obs.columns\n    obs_cols = list(set(adata.obs.columns) - set([col]))\n    \n    one_to_one_mapped_obs_cols = []\n    nunique_in_col = adata.obs[col].nunique()\n    for other_col in obs_cols:\n        if len(adata.obs[[col, other_col]].drop_duplicates()) == nunique_in_col:\n            one_to_one_mapped_obs_cols.append(other_col)\n\n    joining_df = adata.obs[[col] + one_to_one_mapped_obs_cols].drop_duplicates().set_index(col)\n    assert (sum_adata.obs.index == sum_adata.obs.join(joining_df).index).all()\n    sum_adata.obs = sum_adata.obs.join(joining_df)\n    sum_adata.obs.index.name = col\n    sum_adata.obs = sum_adata.obs.reset_index()\n    sum_adata.obs.index = sum_adata.obs.index.astype('str')\n\n    sum_adata.obs = sum_adata.obs.drop(columns=['plate_well_cell_type'])\n    sum_adata.X = np.array(sum_adata.X.todense())\n    sum_adata.X = sum_adata.X.astype('float64')\n    \n#     plate_reordering = {\n#         'plate_0': 'plate_1',\n#         'plate_1': 'plate_2',\n#         'plate_2': 'plate_3',\n#         'plate_3': 'plate_0',\n#         'plate_4': 'plate_4',\n#         'plate_5': 'plate_5',\n#     }\n\n#     sum_adata.obs['plate_name'] = sum_adata.obs['plate_name'].map(plate_reordering)\n    return sum_adata","metadata":{"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def limma_prep_worker(bulk_adata, exp_id):\n    de_pert_cols = [\n        'sm_name',\n        'sm_lincs_id',\n        'SMILES',\n        'dose_uM',\n        'timepoint_hr',\n        'cell_type',\n    ]\n\n    control_compound = 'Dimethyl Sulfoxide'\n\n    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    compound_name_col = de_pert_cols[0]\n    cell_types = bulk_adata.obs['cell_type'].unique()\n\n    with open(f'exp{exp_id}_r_code.sh', 'w') as f:\n\n        for cell_type in tqdm(cell_types, position=1, desc='cell_type'):\n            rename_cell_type = cell_rename_dict[cell_type]\n            cell_type_working_dir = 'exp{}_cell_type_{}'.format(exp_id, rename_cell_type)\n\n            os.makedirs(cell_type_working_dir, exist_ok=True)\n\n            cell_type_selection = bulk_adata.obs['cell_type'].eq(cell_type)\n            cell_type_bulk_adata = bulk_adata[cell_type_selection].copy()\n\n            rpert_mapping = cell_type_bulk_adata.obs[compound_name_col].drop_duplicates() \\\n            .reset_index(drop=True).reset_index() \\\n            .set_index(compound_name_col)['index'].to_dict()\n\n            cell_type_bulk_adata.obs['Rpert'] = cell_type_bulk_adata.obs.apply(\n                lambda row: rpert_mapping[row[compound_name_col]], \n                axis='columns',\n            ).astype('str')\n\n            compound_name_to_Rpert = cell_type_bulk_adata.obs.set_index(compound_name_col)['Rpert'].to_dict()\n            ref_pert = compound_name_to_Rpert[control_compound]\n\n            # save h5ad for each cell type\n            cell_type_bulk_adata.write_h5ad(os.path.join(cell_type_working_dir, 'input.h5ad'))\n\n            random_string = binascii.b2a_hex(os.urandom(15)).decode()\n\n            gc.collect();\n\n            Rscript = '`which Rscript`'\n\n            # for limma fit\n            exec_path_fit = '/kaggle/working/neurips-2023-scripts/limma_fit.r'\n            design_fit = '~0+Rpert+donor_id+plate_name+row'\n\n            input_path = os.path.join(cell_type_working_dir, 'input.h5ad')\n            fit_output_path = os.path.join(cell_type_working_dir, 'limma.rds')\n            fit_plot_output_path = os.path.join(cell_type_working_dir, 'voom.pdf')\n\n            r_code_for_fit_in_bash = f'''\necho '# limma fit for cell type {rename_cell_type}'\n{Rscript} \\\n{exec_path_fit} \\\n--input_h5ad \\\n{input_path} \\\n--design \\\n'{design_fit}' \\\n--fit_output_path \\\n{fit_output_path} \\\n--plot_output_path \\\n{fit_plot_output_path} \n'''\n\n            print (r_code_for_fit_in_bash, file=f)\n\n            # for limma contrast\n            exec_path_contrast = '/kaggle/working/neurips-2023-scripts/limma_contrast.r'\n\n            for pert in cell_type_bulk_adata.obs['Rpert'].unique():\n                if pert == ref_pert:\n                    continue\n                else:\n                    contrast_output_path = os.path.join(cell_type_working_dir, 'pert_{}_contrast_result.csv'.format(pert))\n                    contrast = 'Rpert'+pert+'-Rpert'+ref_pert\n                    r_code_for_contrast_in_bash = f'''\necho '# limma contrast for cell type {rename_cell_type} with Rpert={pert}'\n{Rscript} \\\n{exec_path_contrast} \\\n--input_fit \\\n{fit_output_path} \\\n--contrast \\\n{contrast} \\\n--contrast_output_path \\\n{contrast_output_path}\n'''\n\n                print (r_code_for_contrast_in_bash, file=f)\n\n    return","metadata":{"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"---","metadata":{}},{"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}","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# adata_path = '/kaggle/input/op2-02-pseudobulk-v2/singlecell_adata.h5ad'\n# adata = sc.read_h5ad(adata_path)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"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']\nreduction_names = ['none']\n\nadata_path = '/kaggle/input/op2-02-pseudobulk-v2/singlecell_adata.h5ad'\n\nfor exp_id, (preprocess_name, reduction_name) in enumerate(tqdm(list(product(preprocess_names, reduction_names)))):\n    print ('[DEBUG]', exp_id, preprocess_name, reduction_name)\n    for rename_cell_type in cell_rename_dict.values():\n        cell_type_working_dir = 'exp{}_cell_type_{}'.format(exp_id, rename_cell_type)\n        os.makedirs(cell_type_working_dir, exist_ok=True)\n\n    print ('[DEBUG] preprocessing')\n    adata = preprocessing(adata_path, preprocess_name, reduction_name)\n    gc.collect();\n    print ('[DEBUG] pseudobulk')\n    bulk_adata = sum_by(adata)\n    gc.collect();\n    print ('[DEBUG] limma preprocessing')\n    limma_prep_worker(bulk_adata, exp_id)\n    gc.collect();\n    ","metadata":{"_kg_hide-output":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"---\n---\n---","metadata":{}},{"cell_type":"markdown","source":"# Setup R environment","metadata":{}},{"cell_type":"code","source":"r_dependencies = '''\n- r-base=4.3.1=h29c4799_3\n- r-base64enc=0.1_3=r43h57805ef_1006\n- r-cli=3.6.1=r43ha503ecb_1\n- r-cpp11=0.4.6=r43hc72bb7e_0\n- r-crayon=1.5.2=r43hc72bb7e_2\n- r-digest=0.6.33=r43ha503ecb_0\n- r-ellipsis=0.3.2=r43h57805ef_2\n- r-evaluate=0.21=r43hc72bb7e_1\n- r-fansi=1.0.4=r43h57805ef_1\n- r-fastmap=1.1.1=r43ha503ecb_1\n- r-glue=1.6.2=r43h57805ef_2\n- r-htmltools=0.5.6=r43ha503ecb_0\n- r-irdisplay=1.1=r43hd8ed1ab_2\n- r-irkernel=1.3.2=r43h785f33e_1\n- r-jsonlite=1.8.7=r43h57805ef_0\n- r-lifecycle=1.0.3=r43hc72bb7e_2\n- r-pbdzmq=0.3_9=r43ha81a24b_1\n- r-pillar=1.9.0=r43hc72bb7e_1\n- r-ragg=1.2.5=r43h85cdef0_2\n- r-repr=1.1.6=r43h785f33e_1\n- r-rlang=1.1.1=r43ha503ecb_1\n- r-systemfonts=1.0.4=r43haf97adc_2\n- r-textshaping=0.3.6=r43h24cd192_6\n- r-utf8=1.2.3=r43h57805ef_1\n- r-uuid=1.1_1=r43h57805ef_0\n- r-vctrs=0.6.3=r43ha503ecb_0\n'''\n\n\npackages = ' '.join(\n    pd.Series(\n        r_dependencies.replace('- ', '').split('\\n')[1:-1]\n    ).apply(lambda x: x.split('=')[0]).to_list()\n) + ' r-optparse r-reticulate r-rcpp gcc_impl_linux-64=13'\n\nos.system('conda install -q -c conda-forge -y {}'.format(packages))","metadata":{"_kg_hide-output":true,"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"r_install = '''\ninstall.packages(\"remotes\", repos = \"http://cran.us.r-project.org\", quiet=TRUE)\ninstall.packages(\"BiocManager\", repos = \"http://cran.us.r-project.org\", quiet=TRUE)\nBiocManager::install(\"edgeR\")\nremotes::install_cran(c(\"anndata\"), repos = \"https://cran.rstudio.com\")\n# install.packages(\"optparse\", repos = \"http://cran.us.r-project.org\", quiet=TRUE)   # package ‘optparse’ is not available for this version of R\nremotes::install_version(\"reticulate\", version = \"1.22\", repos = \"http://cran.us.r-project.org\")\n'''\n\nwith open('r_install.r', 'w') as f:\n    print (r_install, file=f)","metadata":{"_kg_hide-input":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"! `which Rscript` r_install.r","metadata":{"_kg_hide-output":true,"_kg_hide-input":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# RUN Limma","metadata":{}},{"cell_type":"code","source":"!for file in exp*_r_code.sh; do bash $file; done;","metadata":{"_kg_hide-output":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"!rm -rf /kaggle/working/neurips-2023-scripts","metadata":{"_kg_hide-output":true,"_kg_hide-input":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"---\n---\n---","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    \n    return de_adata","metadata":{"_kg_hide-input":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Compared with the `de_train`","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.sort_values(['sm_name', 'cell_type']).index\nkaggle_train_de_adata = kaggle_train_de_adata[sorting_index].copy()\n\ndel (kaggle_train_de_df); gc.collect();","metadata":{"execution":{"iopub.status.busy":"2023-10-16T06:50:31.394968Z","iopub.execute_input":"2023-10-16T06:50:31.395716Z","iopub.status.idle":"2023-10-16T06:50:35.597880Z","shell.execute_reply.started":"2023-10-16T06:50:31.395673Z","shell.execute_reply":"2023-10-16T06:50:35.596491Z"},"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_count":null,"outputs":[]},{"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_count":null,"outputs":[]},{"cell_type":"code","source":"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 = '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[:, kaggle_train_de_adata.var.index].X),\n        pd.DataFrame(kaggle_train_de_adata.X),\n    )\n    fig.suptitle('{}: MRRMSE={:.3f}'.format(preprocess_name, diff), fontsize=16)\n    plt.show()\n\n    gc.collect();\n    print ()\n    print ()\n","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}