{"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 scvi-tools\n!pip install scanpy\n!pip install scikit-misc","metadata":{"scrolled":true,"_kg_hide-output":true,"execution":{"iopub.status.busy":"2023-10-23T10:57:51.700376Z","iopub.execute_input":"2023-10-23T10:57:51.700769Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import gc\n\nimport numpy as np\nimport pandas as pd\n\nimport scanpy as sc\nimport scvi\n\nimport matplotlib.pyplot as plt\nimport seaborn as sns","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"adata = sc.read_h5ad('/kaggle/input/op2-02-pseudobulk-v2/singlecell_adata.h5ad')\nbulk_adata = sc.read_h5ad('/kaggle/input/op2-02-pseudobulk-v2/bulk_adata.h5ad')\n\nadata = adata[:, bulk_adata.var.index].copy()\ndel (bulk_adata); gc.collect();","metadata":{"execution":{"iopub.status.busy":"2023-10-18T03:11:37.099572Z","iopub.execute_input":"2023-10-18T03:11:37.099936Z","iopub.status.idle":"2023-10-18T03:12:24.101092Z","shell.execute_reply.started":"2023-10-18T03:11:37.099908Z","shell.execute_reply":"2023-10-18T03:12:24.099887Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sc.pp.filter_genes(adata, min_counts=3)\n# sc.pp.normalize_total(adata, target_sum=1e4)\n# sc.pp.log1p(adata)","metadata":{"execution":{"iopub.status.busy":"2023-10-18T03:11:29.090170Z","iopub.execute_input":"2023-10-18T03:11:29.090844Z","iopub.status.idle":"2023-10-18T03:11:37.092477Z","shell.execute_reply.started":"2023-10-18T03:11:29.090805Z","shell.execute_reply":"2023-10-18T03:11:37.091261Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# sc.pp.highly_variable_genes(\n#     adata,\n#     n_top_genes=1200,\n#     subset=True,\n#     flavor=\"seurat_v3\",\n#     batch_key=\"donor_id\",\n# )","metadata":{"execution":{"iopub.status.busy":"2023-10-18T03:11:37.093680Z","iopub.execute_input":"2023-10-18T03:11:37.094000Z","iopub.status.idle":"2023-10-18T03:11:37.098162Z","shell.execute_reply.started":"2023-10-18T03:11:37.093974Z","shell.execute_reply":"2023-10-18T03:11:37.097446Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"scvi.model.SCVI.setup_anndata(\n    adata,\n    batch_key='donor_id',\n    labels_key='sm_lincs_id',\n    categorical_covariate_keys=[\"plate_name\", \"library_id\"],\n)\n\nmodel = scvi.model.SCVI(adata)","metadata":{"execution":{"iopub.status.busy":"2023-10-17T10:12:03.574087Z","iopub.execute_input":"2023-10-17T10:12:03.574523Z","iopub.status.idle":"2023-10-17T10:12:03.946156Z","shell.execute_reply.started":"2023-10-17T10:12:03.574493Z","shell.execute_reply":"2023-10-17T10:12:03.945053Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"model.module","metadata":{"execution":{"iopub.status.busy":"2023-10-17T10:44:27.466067Z","iopub.execute_input":"2023-10-17T10:44:27.466526Z","iopub.status.idle":"2023-10-17T10:44:27.474240Z","shell.execute_reply.started":"2023-10-17T10:44:27.466481Z","shell.execute_reply":"2023-10-17T10:44:27.473111Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"model.train()","metadata":{"execution":{"iopub.status.busy":"2023-10-17T10:27:49.411522Z","iopub.execute_input":"2023-10-17T10:27:49.411947Z","iopub.status.idle":"2023-10-17T10:27:53.070442Z","shell.execute_reply.started":"2023-10-17T10:27:49.411914Z","shell.execute_reply":"2023-10-17T10:27:53.069221Z"},"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":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"model.module.eval()\n\nde_dfs = []\nfor cell_type in cell_rename_dict.keys():\n    \n    de_df = model.differential_expression(\n        adata[adata.obs.query('cell_type==@cell_type').index].copy(),\n        groupby='sm_lincs_id', group2='LSM-36361',\n        batch_size=128\n    ).assign(cell_type=cell_type)\n    de_dfs.append(de_df)\n\n    \nde_dfs = pd.concat(de_dfs)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"de_dfs.to_csv('de_dfs.csv')","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"---","metadata":{}},{"cell_type":"code","source":"de_train = pd.read_parquet('/kaggle/input/open-problems-single-cell-perturbations/de_train.parquet').set_index(['cell_type', 'sm_lincs_id']).drop(['sm_name', 'SMILES', 'control'], axis=1)","metadata":{"execution":{"iopub.status.busy":"2023-10-18T03:24:53.930016Z","iopub.execute_input":"2023-10-18T03:24:53.930421Z","iopub.status.idle":"2023-10-18T03:24:56.809732Z","shell.execute_reply.started":"2023-10-18T03:24:53.930393Z","shell.execute_reply":"2023-10-18T03:24:56.808564Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"de_df = de_dfs.drop(['pseudocounts', 'delta', 'group2', 'comparison'], axis=1).reset_index()\nde_df[\"log10_p\"] = -np.log10(de_df[\"proba_not_de\"]+1e-3)\nde_df[\"lfc_neg_log10_p\"] = -np.log10(de_df[\"proba_not_de\"]+1e-3) * de_df['lfc_mean'].apply(np.sign)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"comp_df = de_df.pivot(index=['cell_type', 'group1'], columns='gene', values='lfc_neg_log10_p').loc[de_train.index]","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sns.scatterplot(\n    data=de_df,\n    x='lfc_mean', y='log10_p',\n    hue='cell_type',\n    edgecolor='k', lw=0.1, alpha=0.5\n)\n\nplt.show()","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"share_gene = list(set(comp_df.columns) & set(de_train.columns))\n\ndef 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":"len(share_gene)","metadata":{"execution":{"iopub.status.busy":"2023-10-18T11:11:01.582539Z","iopub.execute_input":"2023-10-18T11:11:01.582977Z","iopub.status.idle":"2023-10-18T11:11:01.605784Z","shell.execute_reply.started":"2023-10-18T11:11:01.582946Z","shell.execute_reply":"2023-10-18T11:11:01.604055Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"mrrmse(comp_df.loc[:, share_gene], de_train.loc[:, share_gene])","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":[]}]}