{"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":"# TOC \n### (This notebook is under construction !!!)\n\n---\n\n#### Version 1\n- Copy the code from openproblems-bio repo\n\n#### Version 2,3\n- code for limma\n\n#### Version 4\n- save adata.h5ad\n\n#### Version 5,6,7,8,9,10\n- Add scatter plot\n- Debug for the bulk data mismatch","metadata":{}},{"cell_type":"markdown","source":"## The original code: https://github.com/openproblems-bio/neurips-2023-scripts/blob/main/compute_de.ipynb\n### I made some small changes to make the code can be running on Kaggle (except Limma, which is using R)","metadata":{}},{"cell_type":"code","source":"!pip install scanpy","metadata":{"_kg_hide-output":true,"scrolled":true,"_kg_hide-input":false,"execution":{"iopub.status.busy":"2023-10-18T09:52:29.898473Z","iopub.execute_input":"2023-10-18T09:52:29.898878Z","iopub.status.idle":"2023-10-18T09:52:41.878384Z","shell.execute_reply.started":"2023-10-18T09:52:29.898847Z","shell.execute_reply":"2023-10-18T09:52:41.876608Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# This Python 3 environment comes with many helpful analytics libraries installed\n# It is defined by the kaggle/python Docker image: https://github.com/kaggle/docker-python\n# For example, here's several helpful packages to load\n\nimport numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\n\n# Input data files are available in the read-only \"../input/\" directory\n# For example, running this (by clicking run or pressing Shift+Enter) will list all files under the input directory\n\nimport os\nfor dirname, _, filenames in os.walk('/kaggle/input'):\n    for filename in filenames:\n        print(os.path.join(dirname, filename))\n\n# You can write up to 20GB to the current directory (/kaggle/working/) that gets preserved as output when you create a version using \"Save & Run All\" \n# You can also write temporary files to /kaggle/temp/, but they won't be saved outside of the current session","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2023-10-18T09:52:41.881151Z","iopub.execute_input":"2023-10-18T09:52:41.881571Z","iopub.status.idle":"2023-10-18T09:52:41.898154Z","shell.execute_reply.started":"2023-10-18T09:52:41.881538Z","shell.execute_reply":"2023-10-18T09:52:41.896508Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"---","metadata":{}},{"cell_type":"code","source":"import pandas as pd\nimport numpy as np\nimport sklearn\nimport scipy\n\nimport anndata as ad\nimport scanpy as sc\n\nfrom dask import delayed\nfrom dask.distributed import Client, LocalCluster\n\nimport os, binascii, gc\n\nfrom tqdm.auto import tqdm","metadata":{"execution":{"iopub.status.busy":"2023-10-18T09:52:41.899682Z","iopub.execute_input":"2023-10-18T09:52:41.900049Z","iopub.status.idle":"2023-10-18T09:52:41.909007Z","shell.execute_reply.started":"2023-10-18T09:52:41.900021Z","shell.execute_reply":"2023-10-18T09:52:41.907747Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Loading expression data\n\nHere we load expression data (long format) and converting it into an AnnData object (wide sparse format).\n\nYou'll need to increase your instance RAM to at least 64 GB.","metadata":{}},{"cell_type":"code","source":"data_dir = '/kaggle/input/open-problems-single-cell-perturbations'\nadata_train_df = pd.read_parquet(os.path.join(data_dir, 'adata_train.parquet'))\nadata_obs_meta_df = pd.read_csv(os.path.join(data_dir, 'adata_obs_meta.csv'))\n\nadata_train_df['obs_id'] = adata_train_df['obs_id'].astype('category')\nadata_train_df['gene'] = adata_train_df['gene'].astype('category')\n\nobs_ids = adata_train_df['obs_id'].unique()\nobs_id_map = dict(zip(obs_ids, range(len(obs_ids))))\n\ngenes = adata_train_df['gene'].unique()\ngene_map = dict(zip(genes, range(len(genes))))\n\nadata_train_df['obs_index'] = adata_train_df['obs_id'].map(obs_id_map)\nadata_train_df['gene_index'] = adata_train_df['gene'].map(gene_map)\n\nnormalized_counts_values = adata_train_df['normalized_count'].to_numpy()\ncounts_values = adata_train_df['count'].to_numpy()\n\nrow_indices = adata_train_df['obs_index'].to_numpy()\ncol_indices = adata_train_df['gene_index'].to_numpy()\n\ncounts = scipy.sparse.csr_matrix((counts_values, (row_indices, col_indices)))\n\nobs_df = pd.Series(obs_ids, name='obs_id').to_frame()\nvar_df = pd.Series(genes, name='gene').to_frame()\n\nobs_df = obs_df.set_index('obs_id')\nvar_df = var_df.set_index('gene')\n\nobs_df.index = obs_df.index.astype('str')\nvar_df.index = var_df.index.astype('str')\n\ncounts_adata = ad.AnnData(\n    X=counts,\n    obs=obs_df,\n    var=var_df,\n    dtype=np.uint32,\n)\n\nindex_ordering_before_join = counts_adata.obs.index\ncounts_adata.obs = counts_adata.obs.join(adata_obs_meta_df.set_index('obs_id'))\nindex_ordering_after_join = counts_adata.obs.index\nassert (index_ordering_before_join == index_ordering_after_join).all()","metadata":{"execution":{"iopub.status.busy":"2023-10-18T08:30:03.090646Z","iopub.execute_input":"2023-10-18T08:30:03.091088Z","iopub.status.idle":"2023-10-18T08:33:41.310195Z","shell.execute_reply.started":"2023-10-18T08:30:03.091052Z","shell.execute_reply":"2023-10-18T08:33:41.308972Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"counts_adata","metadata":{"execution":{"iopub.status.busy":"2023-10-18T08:33:41.311968Z","iopub.execute_input":"2023-10-18T08:33:41.312505Z","iopub.status.idle":"2023-10-18T08:33:41.322819Z","shell.execute_reply.started":"2023-10-18T08:33:41.312460Z","shell.execute_reply":"2023-10-18T08:33:41.321118Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"counts_adata.write_h5ad('singlecell_adata.h5ad');","metadata":{"execution":{"iopub.status.busy":"2023-10-18T08:33:41.325411Z","iopub.execute_input":"2023-10-18T08:33:41.326003Z","iopub.status.idle":"2023-10-18T08:33:49.344936Z","shell.execute_reply.started":"2023-10-18T08:33:41.325951Z","shell.execute_reply":"2023-10-18T08:33:49.343334Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Pseudobulking counts by cell type","metadata":{}},{"cell_type":"code","source":"from scipy import sparse\n\ndef sum_by(adata: ad.AnnData, col: str) -> 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    return sum_adata","metadata":{"execution":{"iopub.status.busy":"2023-10-18T08:33:49.346953Z","iopub.execute_input":"2023-10-18T08:33:49.347358Z","iopub.status.idle":"2023-10-18T08:33:49.363806Z","shell.execute_reply.started":"2023-10-18T08:33:49.347323Z","shell.execute_reply":"2023-10-18T08:33:49.362876Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"counts_adata.obs['plate_well_cell_type'] = counts_adata.obs['plate_name'].astype('str') \\\n    + '_' + counts_adata.obs['well'].astype('str') \\\n    + '_' + counts_adata.obs['cell_type'].astype('str')\ncounts_adata.obs['plate_well_cell_type'] = counts_adata.obs['plate_well_cell_type'].astype('category')\n\nbulk_adata = sum_by(counts_adata, 'plate_well_cell_type')\nbulk_adata.obs = bulk_adata.obs.drop(columns=['plate_well_cell_type'])\nbulk_adata.X = np.array(bulk_adata.X.todense())\nbulk_adata.X = bulk_adata.X.astype('float64')\nbulk_adata = bulk_adata.copy()","metadata":{"execution":{"iopub.status.busy":"2023-10-18T08:33:49.365094Z","iopub.execute_input":"2023-10-18T08:33:49.365930Z","iopub.status.idle":"2023-10-18T08:33:55.106804Z","shell.execute_reply.started":"2023-10-18T08:33:49.365896Z","shell.execute_reply":"2023-10-18T08:33:55.105631Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"bulk_adata.write_h5ad('raw_bulk_adata.h5ad');","metadata":{"execution":{"iopub.status.busy":"2023-10-18T08:33:55.108338Z","iopub.execute_input":"2023-10-18T08:33:55.108727Z","iopub.status.idle":"2023-10-18T08:33:55.793560Z","shell.execute_reply.started":"2023-10-18T08:33:55.108680Z","shell.execute_reply":"2023-10-18T08:33:55.792280Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"bulk_adata = sc.read_h5ad('raw_bulk_adata.h5ad');\nbulk_adata","metadata":{"execution":{"iopub.status.busy":"2023-10-18T10:40:09.614760Z","iopub.execute_input":"2023-10-18T10:40:09.615597Z","iopub.status.idle":"2023-10-18T10:40:11.633776Z","shell.execute_reply.started":"2023-10-18T10:40:09.615550Z","shell.execute_reply":"2023-10-18T10:40:11.632209Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# stop reordering to debug\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# }","metadata":{"execution":{"iopub.status.busy":"2023-10-18T10:40:11.636366Z","iopub.execute_input":"2023-10-18T10:40:11.638286Z","iopub.status.idle":"2023-10-18T10:40:11.645558Z","shell.execute_reply.started":"2023-10-18T10:40:11.638142Z","shell.execute_reply":"2023-10-18T10:40:11.643770Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Loading pseudobulked counts from correctly filtered AnnData (@Daniel: delete this section once the fixed RNA expression AnnData is uploaded to Kaggle)","metadata":{}},{"cell_type":"code","source":"fixed_bulk_adata = sc.read_h5ad('/kaggle/input/open-problems-single-cell-perturbations-optional/train_or_control_bulk_by_cell_type_adata.h5ad')\nfixed_bulk_adata.X = fixed_bulk_adata.layers['counts']","metadata":{"execution":{"iopub.status.busy":"2023-10-18T10:40:12.003300Z","iopub.execute_input":"2023-10-18T10:40:12.003833Z","iopub.status.idle":"2023-10-18T10:40:18.635211Z","shell.execute_reply.started":"2023-10-18T10:40:12.003796Z","shell.execute_reply":"2023-10-18T10:40:18.633832Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"new_plate_names = bulk_adata.obs['plate_name'].sort_values().unique()\noriginal_plate_names = fixed_bulk_adata.obs['plate_name'].sort_values().unique()\nplate_name_map = dict(zip(original_plate_names, new_plate_names))\n\nfixed_bulk_adata.obs['plate_name'] = fixed_bulk_adata.obs['plate_name'].map(plate_name_map)\n# fixed_bulk_adata.obs['plate_name'] = fixed_bulk_adata.obs['plate_name'].map(plate_reordering)","metadata":{"execution":{"iopub.status.busy":"2023-10-18T10:40:18.638452Z","iopub.execute_input":"2023-10-18T10:40:18.639379Z","iopub.status.idle":"2023-10-18T10:40:18.651589Z","shell.execute_reply.started":"2023-10-18T10:40:18.639328Z","shell.execute_reply":"2023-10-18T10:40:18.650508Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"new_donor_ids = bulk_adata.obs['donor_id'].sort_values().unique()\noriginal_donor_ids = fixed_bulk_adata.obs['raw_cell_id'].sort_values().unique()\ndonor_id_map = dict(zip(original_donor_ids, new_donor_ids))\n\nfixed_bulk_adata.obs['donor_id'] = fixed_bulk_adata.obs['raw_cell_id'].map(donor_id_map)","metadata":{"execution":{"iopub.status.busy":"2023-10-18T10:40:18.653186Z","iopub.execute_input":"2023-10-18T10:40:18.654759Z","iopub.status.idle":"2023-10-18T10:40:18.670097Z","shell.execute_reply.started":"2023-10-18T10:40:18.654700Z","shell.execute_reply":"2023-10-18T10:40:18.668400Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"lincs_id_mapping_df = pd.read_parquet('/kaggle/input/open-problems-single-cell-perturbations-optional/lincs_id_compound_mapping.parquet')","metadata":{"execution":{"iopub.status.busy":"2023-10-18T10:40:18.673134Z","iopub.execute_input":"2023-10-18T10:40:18.673811Z","iopub.status.idle":"2023-10-18T10:40:18.703343Z","shell.execute_reply.started":"2023-10-18T10:40:18.673771Z","shell.execute_reply":"2023-10-18T10:40:18.701892Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"compound_id_to_sm_lincs_id = lincs_id_mapping_df.set_index('compound_id')['sm_lincs_id'].to_dict()\ncompound_id_to_sm_name = lincs_id_mapping_df.set_index('compound_id')['sm_name'].to_dict()\ncompound_id_to_smiles = lincs_id_mapping_df.set_index('compound_id')['smiles'].to_dict()\n\nfixed_bulk_adata.obs['sm_lincs_id'] = \\\n    fixed_bulk_adata.obs['compound_id'].map(compound_id_to_sm_lincs_id)\nfixed_bulk_adata.obs['sm_name'] = \\\n    fixed_bulk_adata.obs['compound_id'].map(compound_id_to_sm_name)\nfixed_bulk_adata.obs['SMILES'] = \\\n    fixed_bulk_adata.obs['compound_id'].map(compound_id_to_smiles)","metadata":{"execution":{"iopub.status.busy":"2023-10-18T10:40:18.705008Z","iopub.execute_input":"2023-10-18T10:40:18.705394Z","iopub.status.idle":"2023-10-18T10:40:18.721262Z","shell.execute_reply.started":"2023-10-18T10:40:18.705361Z","shell.execute_reply":"2023-10-18T10:40:18.720221Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sorted_index = fixed_bulk_adata.obs.sort_values(['plate_name', 'sm_name', 'cell_type']).index\nfixed_bulk_adata = fixed_bulk_adata[sorted_index].copy()\nsorted_index = bulk_adata.obs.sort_values(['plate_name', 'sm_name', 'cell_type']).index\nbulk_adata = bulk_adata[sorted_index].copy()","metadata":{"execution":{"iopub.status.busy":"2023-10-18T10:40:18.722755Z","iopub.execute_input":"2023-10-18T10:40:18.723737Z","iopub.status.idle":"2023-10-18T10:40:19.726652Z","shell.execute_reply.started":"2023-10-18T10:40:18.723704Z","shell.execute_reply":"2023-10-18T10:40:19.724888Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"If the new expression data is uploaded to Kaggle (or placed in '/home/jovyan/kaggle/input/open-problems-multimodal-2023' in the correct format), the assertion below should evaluate to True.","metadata":{}},{"cell_type":"code","source":"# bulk_adata.obs['plate_name'] = bulk_adata.obs['plate_name'].map(plate_reordering)","metadata":{"execution":{"iopub.status.busy":"2023-10-18T10:40:19.728284Z","iopub.execute_input":"2023-10-18T10:40:19.729417Z","iopub.status.idle":"2023-10-18T10:40:19.733932Z","shell.execute_reply.started":"2023-10-18T10:40:19.729386Z","shell.execute_reply":"2023-10-18T10:40:19.732803Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"bulk_adata","metadata":{"execution":{"iopub.status.busy":"2023-10-18T10:40:19.735711Z","iopub.execute_input":"2023-10-18T10:40:19.736109Z","iopub.status.idle":"2023-10-18T10:40:19.755753Z","shell.execute_reply.started":"2023-10-18T10:40:19.736082Z","shell.execute_reply":"2023-10-18T10:40:19.753924Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fixed_bulk_adata","metadata":{"execution":{"iopub.status.busy":"2023-10-18T10:40:19.757791Z","iopub.execute_input":"2023-10-18T10:40:19.759178Z","iopub.status.idle":"2023-10-18T10:40:19.770497Z","shell.execute_reply.started":"2023-10-18T10:40:19.759133Z","shell.execute_reply":"2023-10-18T10:40:19.768908Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# assert np.allclose(fixed_bulk_adata.X, bulk_adata.X)\nnp.allclose(fixed_bulk_adata.X, bulk_adata[:, fixed_bulk_adata.var.index].X)","metadata":{"execution":{"iopub.status.busy":"2023-10-18T10:40:19.774206Z","iopub.execute_input":"2023-10-18T10:40:19.775136Z","iopub.status.idle":"2023-10-18T10:40:21.955796Z","shell.execute_reply.started":"2023-10-18T10:40:19.775099Z","shell.execute_reply":"2023-10-18T10:40:21.954525Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"---","metadata":{}},{"cell_type":"markdown","source":"# Metadata mismatch","metadata":{}},{"cell_type":"code","source":"clean_bulk_adata = bulk_adata[:, fixed_bulk_adata.var.index]","metadata":{"execution":{"iopub.status.busy":"2023-10-18T10:40:23.137105Z","iopub.execute_input":"2023-10-18T10:40:23.137607Z","iopub.status.idle":"2023-10-18T10:40:23.161597Z","shell.execute_reply.started":"2023-10-18T10:40:23.137571Z","shell.execute_reply":"2023-10-18T10:40:23.159879Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"pd.merge(\n    clean_bulk_adata.obs[['cell_type', 'plate_name', 'donor_id', 'sm_lincs_id']].drop_duplicates().reset_index(),\n    fixed_bulk_adata.obs[['cell_type', 'plate_name', 'donor_id', 'sm_lincs_id']].drop_duplicates().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-18T10:40:23.811823Z","iopub.execute_input":"2023-10-18T10:40:23.812319Z","iopub.status.idle":"2023-10-18T10:40:23.885251Z","shell.execute_reply.started":"2023-10-18T10:40:23.812276Z","shell.execute_reply":"2023-10-18T10:40:23.883924Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fixed_bulk_adata.obs[['cell_type', 'plate_name', 'donor_id', 'sm_lincs_id']].drop_duplicates().reset_index().shape","metadata":{"execution":{"iopub.status.busy":"2023-10-18T10:40:27.101484Z","iopub.execute_input":"2023-10-18T10:40:27.101939Z","iopub.status.idle":"2023-10-18T10:40:27.115275Z","shell.execute_reply.started":"2023-10-18T10:40:27.101907Z","shell.execute_reply":"2023-10-18T10:40:27.114043Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"---","metadata":{}},{"cell_type":"markdown","source":"# Debug mismatch","metadata":{}},{"cell_type":"code","source":"(\n    clean_bulk_adata.obs.query('cell_type==\"NK cells\" and col>3')[['cell_type', 'plate_name', 'donor_id', 'sm_lincs_id']]\n    .astype(str)\n    .pivot(index='sm_lincs_id', columns='donor_id', values='plate_name')\n).drop_duplicates()","metadata":{"execution":{"iopub.status.busy":"2023-10-18T10:40:28.632930Z","iopub.execute_input":"2023-10-18T10:40:28.633471Z","iopub.status.idle":"2023-10-18T10:40:28.663718Z","shell.execute_reply.started":"2023-10-18T10:40:28.633403Z","shell.execute_reply":"2023-10-18T10:40:28.662247Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"(\n    fixed_bulk_adata.obs.query('cell_type==\"NK cells\" and col>3')[['cell_type', 'plate_name', 'donor_id', 'sm_lincs_id']]\n    .astype(str)\n    .pivot(index='sm_lincs_id', columns='donor_id', values='plate_name')\n).drop_duplicates()","metadata":{"execution":{"iopub.status.busy":"2023-10-18T10:40:29.353367Z","iopub.execute_input":"2023-10-18T10:40:29.353885Z","iopub.status.idle":"2023-10-18T10:40:29.379945Z","shell.execute_reply.started":"2023-10-18T10:40:29.353852Z","shell.execute_reply":"2023-10-18T10:40:29.379047Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Fix mismatch","metadata":{}},{"cell_type":"code","source":"donor_rename_dict = {\n    'donor_0':'donor_0',\n    'donor_1':'donor_2',\n    'donor_2':'donor_1',\n}\n\nplate_rename_dict = {\n    'plate_0' : 'plate_2',\n    'plate_1' : 'plate_0',\n    'plate_2' : 'plate_1',\n    'plate_3' : 'plate_4',\n    'plate_4' : 'plate_3',\n    'plate_5' : 'plate_5',\n}","metadata":{"execution":{"iopub.status.busy":"2023-10-18T10:40:30.208782Z","iopub.execute_input":"2023-10-18T10:40:30.209855Z","iopub.status.idle":"2023-10-18T10:40:30.216859Z","shell.execute_reply.started":"2023-10-18T10:40:30.209808Z","shell.execute_reply":"2023-10-18T10:40:30.215205Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fixed_bulk_adata.obs['donor_id'] = fixed_bulk_adata.obs['donor_id'].map(donor_rename_dict)\nfixed_bulk_adata.obs['plate_name'] = fixed_bulk_adata.obs['plate_name'].map(plate_rename_dict)","metadata":{"execution":{"iopub.status.busy":"2023-10-18T10:40:30.606197Z","iopub.execute_input":"2023-10-18T10:40:30.606630Z","iopub.status.idle":"2023-10-18T10:40:30.615914Z","shell.execute_reply.started":"2023-10-18T10:40:30.606599Z","shell.execute_reply":"2023-10-18T10:40:30.614574Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# check after fix","metadata":{}},{"cell_type":"code","source":"(\n    fixed_bulk_adata.obs.query('cell_type==\"NK cells\" and col>3')[['cell_type', 'plate_name', 'donor_id', 'sm_lincs_id']]\n    .astype(str)\n    .pivot(index='sm_lincs_id', columns='donor_id', values='plate_name')\n).drop_duplicates()","metadata":{"execution":{"iopub.status.busy":"2023-10-18T10:40:31.437068Z","iopub.execute_input":"2023-10-18T10:40:31.437885Z","iopub.status.idle":"2023-10-18T10:40:31.463603Z","shell.execute_reply.started":"2023-10-18T10:40:31.437843Z","shell.execute_reply":"2023-10-18T10:40:31.461627Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"(\n    clean_bulk_adata.obs.query('cell_type==\"NK cells\" and col>3')[['cell_type', 'plate_name', 'donor_id', 'sm_lincs_id']]\n    .astype(str)\n    .pivot(index='sm_lincs_id', columns='donor_id', values='plate_name')\n).drop_duplicates()","metadata":{"execution":{"iopub.status.busy":"2023-10-18T10:40:32.057250Z","iopub.execute_input":"2023-10-18T10:40:32.058753Z","iopub.status.idle":"2023-10-18T10:40:32.082210Z","shell.execute_reply.started":"2023-10-18T10:40:32.058694Z","shell.execute_reply":"2023-10-18T10:40:32.080853Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"pd.merge(\n    clean_bulk_adata.obs[['cell_type', 'plate_name', 'donor_id', 'sm_lincs_id']].drop_duplicates().reset_index(),\n    fixed_bulk_adata.obs[['cell_type', 'plate_name', 'donor_id', 'sm_lincs_id']].drop_duplicates().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-18T10:40:33.460224Z","iopub.execute_input":"2023-10-18T10:40:33.461629Z","iopub.status.idle":"2023-10-18T10:40:33.482211Z","shell.execute_reply.started":"2023-10-18T10:40:33.461585Z","shell.execute_reply":"2023-10-18T10:40:33.480542Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fixed_bulk_adata.obs[['cell_type', 'plate_name', 'donor_id', 'sm_lincs_id']].drop_duplicates().reset_index().shape","metadata":{"execution":{"iopub.status.busy":"2023-10-18T10:40:34.076499Z","iopub.execute_input":"2023-10-18T10:40:34.076993Z","iopub.status.idle":"2023-10-18T10:40:34.092383Z","shell.execute_reply.started":"2023-10-18T10:40:34.076960Z","shell.execute_reply":"2023-10-18T10:40:34.090888Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"---","metadata":{}},{"cell_type":"markdown","source":"# Value mismatch, with sort (some genes were zero? I'm still working on it ...)","metadata":{}},{"cell_type":"code","source":"sorted_index = fixed_bulk_adata.obs.astype(str).sort_values(['plate_name', 'sm_name', 'cell_type']).index\nfixed_bulk_adata = fixed_bulk_adata[sorted_index].copy()\nsorted_index = clean_bulk_adata.obs.astype(str).sort_values(['plate_name', 'sm_name', 'cell_type']).index\nclean_bulk_adata = clean_bulk_adata[sorted_index].copy()","metadata":{"execution":{"iopub.status.busy":"2023-10-18T10:41:21.399855Z","iopub.execute_input":"2023-10-18T10:41:21.400476Z","iopub.status.idle":"2023-10-18T10:41:22.649390Z","shell.execute_reply.started":"2023-10-18T10:41:21.400405Z","shell.execute_reply":"2023-10-18T10:41:22.647870Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"np.allclose(fixed_bulk_adata.X, clean_bulk_adata.X)","metadata":{"execution":{"iopub.status.busy":"2023-10-18T10:41:22.651722Z","iopub.execute_input":"2023-10-18T10:41:22.652131Z","iopub.status.idle":"2023-10-18T10:41:23.505727Z","shell.execute_reply.started":"2023-10-18T10:41:22.652098Z","shell.execute_reply":"2023-10-18T10:41:23.504160Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"mdf = pd.merge(\n    clean_bulk_adata.obs[['cell_type', 'plate_name', 'donor_id', 'sm_lincs_id']].drop_duplicates().reset_index(),\n    fixed_bulk_adata.obs[['cell_type', 'plate_name', 'donor_id', 'sm_lincs_id']].drop_duplicates().reset_index(),\n    on=['cell_type', 'plate_name', 'donor_id', 'sm_lincs_id'],\n    how='inner'\n)","metadata":{"execution":{"iopub.status.busy":"2023-10-18T10:41:24.350295Z","iopub.execute_input":"2023-10-18T10:41:24.351648Z","iopub.status.idle":"2023-10-18T10:41:24.378797Z","shell.execute_reply.started":"2023-10-18T10:41:24.351601Z","shell.execute_reply":"2023-10-18T10:41:24.376807Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import matplotlib.pyplot as plt\nimport seaborn as sns","metadata":{"execution":{"iopub.status.busy":"2023-10-18T10:41:25.626478Z","iopub.execute_input":"2023-10-18T10:41:25.626949Z","iopub.status.idle":"2023-10-18T10:41:25.633958Z","shell.execute_reply.started":"2023-10-18T10:41:25.626912Z","shell.execute_reply":"2023-10-18T10:41:25.632097Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig, axs = plt.subplots(ncols=5, nrows=2, figsize=(14, 6), gridspec_kw={'hspace':0.6, 'wspace':0.3})\naxs = axs.flatten()\n\nfor i, idx in enumerate(np.random.choice(mdf.index, 10, replace=False)):\n    x = clean_bulk_adata[mdf.loc[idx, 'index_x']].X\n    y = fixed_bulk_adata[mdf.loc[idx, 'index_y']].X\n    axs[i].scatter(x, y)\n    \n    tag = (x == y).all()\n\n    axs[i].set_title('{} - {}\\n{} - {}'.format(mdf.loc[idx, 'cell_type'], mdf.loc[idx, 'sm_lincs_id'],  mdf.loc[idx, 'plate_name'], tag))\n\n    axs[i].set_xscale('log')\n    axs[i].set_yscale('log')\n    \n    x_min, x_max = axs[i].get_xlim()\n    y_min, y_max = axs[i].get_ylim()\n    axs[i].set_xlim(min(x_min, y_min), max(x_max, y_max))\n    axs[i].set_ylim(min(x_min, y_min), max(x_max, y_max))\n\n    ticks = np.power(10, np.arange(np.log(max(x_max, y_max)) / np.log(10)))\n    axs[i].set_xticks(ticks)\n    axs[i].set_yticks(ticks)\n\nplt.show()\n","metadata":{"execution":{"iopub.status.busy":"2023-10-18T10:41:26.380095Z","iopub.execute_input":"2023-10-18T10:41:26.381088Z","iopub.status.idle":"2023-10-18T10:41:31.583806Z","shell.execute_reply.started":"2023-10-18T10:41:26.381047Z","shell.execute_reply":"2023-10-18T10:41:31.582902Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig, axs = plt.subplots(ncols=5, nrows=4, figsize=(14, 12), gridspec_kw={'hspace':0.6, 'wspace':0.3})\naxs = axs.flatten()\n\nfor i, gene in enumerate(np.random.choice(clean_bulk_adata.var.index, 20, replace=False)):\n    x = clean_bulk_adata[:, gene].X\n    y = fixed_bulk_adata[:, gene].X\n    axs[i].scatter(x, y)\n    \n    tag = (x == y).all()\n\n    axs[i].set_title('{} - {}'.format(gene, tag))\n\n    x_min, x_max = axs[i].get_xlim()\n    y_min, y_max = axs[i].get_ylim()\n    axs[i].set_xlim(min(x_min, y_min), max(x_max, y_max))\n    axs[i].set_ylim(min(x_min, y_min), max(x_max, y_max))\n\n    ticks = np.linspace(0, max(x_max, y_max), 4).round()\n    axs[i].set_xticks(ticks)\n    axs[i].set_yticks(ticks)\n\nplt.show()\n","metadata":{"execution":{"iopub.status.busy":"2023-10-18T10:41:31.585798Z","iopub.execute_input":"2023-10-18T10:41:31.586237Z","iopub.status.idle":"2023-10-18T10:41:34.253300Z","shell.execute_reply.started":"2023-10-18T10:41:31.586202Z","shell.execute_reply":"2023-10-18T10:41:34.251849Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"---","metadata":{}},{"cell_type":"code","source":"bulk_adata = fixed_bulk_adata.copy()","metadata":{"execution":{"iopub.status.busy":"2023-10-09T02:17:56.685759Z","iopub.execute_input":"2023-10-09T02:17:56.686231Z","iopub.status.idle":"2023-10-09T02:17:57.050883Z","shell.execute_reply.started":"2023-10-09T02:17:56.686197Z","shell.execute_reply":"2023-10-09T02:17:57.049600Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"bulk_adata.write_h5ad('bulk_adata.h5ad');\ndel fixed_bulk_adata; gc.collect();","metadata":{"execution":{"iopub.status.busy":"2023-10-09T02:17:57.410391Z","iopub.execute_input":"2023-10-09T02:17:57.410739Z","iopub.status.idle":"2023-10-09T02:17:59.564995Z","shell.execute_reply.started":"2023-10-09T02:17:57.410712Z","shell.execute_reply":"2023-10-09T02:17:59.563880Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"---","metadata":{}},{"cell_type":"code","source":"de_pert_cols = [\n    'sm_name',\n    'sm_lincs_id',\n    'SMILES',\n    'dose_uM',\n    'timepoint_hr',\n    'cell_type',\n]\n\ncontrol_compound = 'Dimethyl Sulfoxide'","metadata":{"execution":{"iopub.status.busy":"2023-10-09T02:18:08.708208Z","iopub.execute_input":"2023-10-09T02:18:08.708619Z","iopub.status.idle":"2023-10-09T02:18:08.714359Z","shell.execute_reply.started":"2023-10-09T02:18:08.708590Z","shell.execute_reply":"2023-10-09T02:18:08.712987Z"},"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}","metadata":{"execution":{"iopub.status.busy":"2023-10-09T02:19:23.082874Z","iopub.execute_input":"2023-10-09T02:19:23.083329Z","iopub.status.idle":"2023-10-09T02:19:23.089409Z","shell.execute_reply.started":"2023-10-09T02:19:23.083296Z","shell.execute_reply":"2023-10-09T02:19:23.087987Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"---\n---\n---","metadata":{}},{"cell_type":"markdown","source":"## Running Limma (using R, not running in this kernal)\n\nhttps://github.com/openproblems-bio/neurips-2023-scripts","metadata":{}},{"cell_type":"markdown","source":"# I write some glue code to run Limma in Kaggle","metadata":{}},{"cell_type":"code","source":"compound_name_col = de_pert_cols[0]\ncell_types = bulk_adata.obs['cell_type'].unique()\n\nwith open('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 = 'cell_type_{}'.format(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        \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\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","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"! head r_code.sh","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"---\n---\n---","metadata":{}},{"cell_type":"markdown","source":"## The original code","metadata":{}},{"cell_type":"markdown","source":"```python\nimport limma_utils\nimport importlib\nimportlib.reload(limma_utils)\n\ndef _run_limma_for_cell_type(bulk_adata):\n    import limma_utils\n    bulk_adata = bulk_adata.copy()\n    \n    compound_name_col = de_pert_cols[0]\n    \n    # limma doesn't like dashes etc. in the compound names\n    rpert_mapping = 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    bulk_adata.obs['Rpert'] = 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 = bulk_adata.obs.set_index(compound_name_col)['Rpert'].to_dict()\n    ref_pert = compound_name_to_Rpert[control_compound]\n            \n    random_string = binascii.b2a_hex(os.urandom(15)).decode()\n    \n    \n    limma_utils.limma_fit(\n        bulk_adata, \n        design='~0+Rpert+donor_id+plate_name+row',\n        output_path=f'output/{random_string}_limma.rds',\n        plot_output_path=f'output/{random_string}_voom',\n        exec_path='limma_fit.r',\n        verbose=True,\n    )\n\n    pert_de_dfs = []\n    \n\n\n    for pert in bulk_adata.obs['Rpert'].unique():\n        if pert == ref_pert:\n            continue\n\n        pert_de_df = limma_utils.limma_contrast(\n            fit_path=f'output/{random_string}_limma.rds',\n            contrast='Rpert'+pert+'-Rpert'+ref_pert,\n            exec_path='limma_contrast.r',\n        )\n\n        pert_de_df['Rpert'] = pert\n\n        pert_obs = bulk_adata.obs[bulk_adata.obs['Rpert'].eq(pert)]\n        for col in de_pert_cols:\n            pert_de_df[col] = pert_obs[col].unique()[0]\n        pert_de_dfs.append(pert_de_df)\n\n    de_df = pd.concat(pert_de_dfs, axis=0)\n\n    try:\n        os.remove(f'output/{random_string}_limma.rds')\n        os.remove(f'output/{random_string}_voom')\n    except FileNotFoundError:\n        pass\n    \n    return de_df\n\nrun_limma_for_cell_type = delayed(_run_limma_for_cell_type)\n\n```","metadata":{"execution":{"iopub.status.busy":"2023-10-08T06:55:12.725215Z","iopub.execute_input":"2023-10-08T06:55:12.726345Z","iopub.status.idle":"2023-10-08T06:55:12.763982Z","shell.execute_reply.started":"2023-10-08T06:55:12.726309Z","shell.execute_reply":"2023-10-08T06:55:12.762548Z"}}},{"cell_type":"markdown","source":"---","metadata":{}},{"cell_type":"markdown","source":"```python\n%%capture\n\ncluster = LocalCluster(\n    n_workers=6,\n    processes=True,\n    threads_per_worker=1,\n    memory_limit='20GB',\n)\n\nc = Client(cluster)\n\n```","metadata":{"execution":{"iopub.status.busy":"2023-10-08T07:06:51.961065Z","iopub.execute_input":"2023-10-08T07:06:51.961439Z","iopub.status.idle":"2023-10-08T07:06:51.970275Z","shell.execute_reply.started":"2023-10-08T07:06:51.961411Z","shell.execute_reply":"2023-10-08T07:06:51.968625Z"}}},{"cell_type":"markdown","source":"---","metadata":{}},{"cell_type":"markdown","source":"```python\n%%capture\n\ncell_types = bulk_adata.obs['cell_type'].unique()\nde_dfs = []\n\nfor cell_type in cell_types:\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    de_df = run_limma_for_cell_type(cell_type_bulk_adata)\n    \n    de_dfs.append(de_df)\n\nde_dfs = c.compute(de_dfs, sync=True)\nde_df = pd.concat(de_dfs)\n```","metadata":{"execution":{"iopub.status.busy":"2023-10-08T07:06:36.395041Z","iopub.execute_input":"2023-10-08T07:06:36.395482Z","iopub.status.idle":"2023-10-08T07:06:36.405811Z","shell.execute_reply.started":"2023-10-08T07:06:36.395454Z","shell.execute_reply":"2023-10-08T07:06:36.404791Z"}}},{"cell_type":"markdown","source":"---","metadata":{}},{"cell_type":"markdown","source":"```python\nde_df.to_csv('de_df.csv.gz', index=False)\n```","metadata":{}},{"cell_type":"markdown","source":"---\n---\n---","metadata":{}},{"cell_type":"code","source":"de_df = pd.read_csv('/kaggle/input/open-problems-single-cell-perturbations-optional/de_df.csv')","metadata":{"execution":{"iopub.status.busy":"2023-10-08T08:04:33.214079Z","iopub.execute_input":"2023-10-08T08:04:33.214496Z","iopub.status.idle":"2023-10-08T08:05:48.421811Z","shell.execute_reply.started":"2023-10-08T08:04:33.214465Z","shell.execute_reply":"2023-10-08T08:05:48.420449Z"},"trusted":true},"execution_count":null,"outputs":[]},{"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(tdf.values, obs=pd.DataFrame(tdf.index.to_frame().values, columns=tdf.index.names), dtype=np.float64)\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":{"execution":{"iopub.status.busy":"2023-10-08T08:06:14.992119Z","iopub.execute_input":"2023-10-08T08:06:14.992533Z","iopub.status.idle":"2023-10-08T08:06:15.003664Z","shell.execute_reply.started":"2023-10-08T08:06:14.992496Z","shell.execute_reply":"2023-10-08T08:06:15.002283Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"de_adata = convert_de_df_to_anndata(de_df, de_pert_cols, 0.05)","metadata":{"execution":{"iopub.status.busy":"2023-10-08T08:06:15.847709Z","iopub.execute_input":"2023-10-08T08:06:15.848146Z","iopub.status.idle":"2023-10-08T08:08:39.198742Z","shell.execute_reply.started":"2023-10-08T08:06:15.848111Z","shell.execute_reply":"2023-10-08T08:08:39.197746Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"kaggle_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(kaggle_train_de_df.values, obs=pd.DataFrame(kaggle_train_de_df.index.to_frame().values, columns=kaggle_train_de_df.index.names))","metadata":{"execution":{"iopub.status.busy":"2023-10-08T08:10:27.503627Z","iopub.execute_input":"2023-10-08T08:10:27.504083Z","iopub.status.idle":"2023-10-08T08:10:28.960914Z","shell.execute_reply.started":"2023-10-08T08:10:27.504052Z","shell.execute_reply":"2023-10-08T08:10:28.959768Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sorting_index = kaggle_train_de_adata.obs.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-08T08:10:30.606948Z","iopub.execute_input":"2023-10-08T08:10:30.607596Z","iopub.status.idle":"2023-10-08T08:10:30.848183Z","shell.execute_reply.started":"2023-10-08T08:10:30.607565Z","shell.execute_reply":"2023-10-08T08:10:30.846943Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"de_adata.obs.index = de_adata.obs.index.astype('str')\n\nsorting_index = de_adata.obs.sort_values(['sm_name', 'cell_type']).index\nde_adata = de_adata[sorting_index].copy()","metadata":{"execution":{"iopub.status.busy":"2023-10-08T08:10:32.355364Z","iopub.execute_input":"2023-10-08T08:10:32.355815Z","iopub.status.idle":"2023-10-08T08:10:32.842918Z","shell.execute_reply.started":"2023-10-08T08:10:32.355780Z","shell.execute_reply":"2023-10-08T08:10:32.841854Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"(de_adata.var.index == kaggle_train_de_adata.var.index).all()","metadata":{"execution":{"iopub.status.busy":"2023-10-08T08:10:33.388465Z","iopub.execute_input":"2023-10-08T08:10:33.388878Z","iopub.status.idle":"2023-10-08T08:10:33.399582Z","shell.execute_reply.started":"2023-10-08T08:10:33.388847Z","shell.execute_reply":"2023-10-08T08:10:33.398245Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import matplotlib.pyplot as plt\nimport seaborn as sns","metadata":{"execution":{"iopub.status.busy":"2023-10-08T08:10:35.286781Z","iopub.execute_input":"2023-10-08T08:10:35.287199Z","iopub.status.idle":"2023-10-08T08:10:35.641831Z","shell.execute_reply.started":"2023-10-08T08:10:35.287171Z","shell.execute_reply":"2023-10-08T08:10:35.640297Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sns.scatterplot(\n    x=kaggle_train_de_adata.X[100],\n    y=de_adata.X[100],\n)","metadata":{"execution":{"iopub.status.busy":"2023-10-08T08:10:35.644061Z","iopub.execute_input":"2023-10-08T08:10:35.644420Z","iopub.status.idle":"2023-10-08T08:10:36.041865Z","shell.execute_reply.started":"2023-10-08T08:10:35.644390Z","shell.execute_reply":"2023-10-08T08:10:36.040522Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"---\n---\n---","metadata":{}},{"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":"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.shape[0]), 10)):\n    sns.scatterplot(\n        x=kaggle_train_de_adata.X[idx],\n        y=de_adata.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\ndiff = mrrmse(\n    pd.DataFrame(de_adata.X),\n    pd.DataFrame(kaggle_train_de_adata.X),\n)\nfig.suptitle('MRRMSE={:.3f}'.format(diff), fontsize=16)\nplt.show()","metadata":{},"execution_count":null,"outputs":[]}]}