{"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##### - adata -> pseudobulk\n##### - log10, genes mean vs std\n##### - design matrix for limma","metadata":{"execution":{"iopub.status.busy":"2023-09-25T06:50:27.716285Z","iopub.execute_input":"2023-09-25T06:50:27.716645Z","iopub.status.idle":"2023-09-25T06:50:27.751224Z","shell.execute_reply.started":"2023-09-25T06:50:27.716614Z","shell.execute_reply":"2023-09-25T06:50:27.74964Z"}}},{"cell_type":"code","source":"!pip install pydeseq2","metadata":{},"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\n\nimport matplotlib.pyplot as plt\n\nimport statsmodels.api as sm","metadata":{"execution":{"iopub.status.busy":"2023-09-28T02:40:16.994668Z","iopub.execute_input":"2023-09-28T02:40:16.995041Z","iopub.status.idle":"2023-09-28T02:40:17.00513Z","shell.execute_reply.started":"2023-09-28T02:40:16.995011Z","shell.execute_reply":"2023-09-28T02:40:17.003883Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# adata -> pseudobulk","metadata":{}},{"cell_type":"markdown","source":"```python\nadata_train = pd.read_parquet('/kaggle/input/open-problems-single-cell-perturbations/adata_train.parquet')\nadata_obs_meta = pd.read_csv('/kaggle/input/open-problems-single-cell-perturbations/adata_obs_meta.csv')\n\npseudobulk = (\n    # merge adata with obs metadata\n    pd.merge(\n        adata_train,\n        adata_obs_meta,\n        on='obs_id',\n        how='left'\n    )\n    # pseudobulk on [library] x [plate] x [donor] x [cell type] x [compound] x [gene]\n    .groupby(['library_id', 'plate_name', 'donor_id', 'cell_type', 'sm_lincs_id', 'gene'])[['count', 'normalized_count']]\n    # averaging \n    .mean()\n    # convert (pivot) to [sample] x [gene] table\n    .reset_index('gene')\n    .pivot(columns='gene', values='count')\n    .fillna(0)\n    .reset_index()\n)\n\npseudobulk.to_csv('pseudobulk.csv', index=False)\n```","metadata":{"execution":{"iopub.status.busy":"2023-09-22T10:45:50.601591Z","iopub.execute_input":"2023-09-22T10:45:50.602189Z","iopub.status.idle":"2023-09-22T10:45:50.61292Z","shell.execute_reply.started":"2023-09-22T10:45:50.602153Z","shell.execute_reply":"2023-09-22T10:45:50.611672Z"}}},{"cell_type":"code","source":"# pseudobulk with mean counts, convert it to raw (sum) counts\n\nadata_obs_meta = pd.read_csv('/kaggle/input/open-problems-single-cell-perturbations/adata_obs_meta.csv')\ncount_df = adata_obs_meta.groupby(['library_id', 'plate_name', 'donor_id', 'cell_type', 'sm_lincs_id'])['obs_id'].count()\n\npseudobulk = pd.read_csv('/kaggle/input/op2-pseudobulk/pseudobulk.csv').set_index(['library_id', 'plate_name', 'donor_id', 'cell_type', 'sm_lincs_id'])\nsum_count_df = (pseudobulk.T*count_df.loc[pseudobulk.index]).astype(int).T\n\npseudobulk = pseudobulk.reset_index()\npseudobulk['group'] = pseudobulk.apply(lambda row: '{}x{}'.format(row['cell_type'], row['sm_lincs_id']), axis=1).values\npseudobulk = pseudobulk.set_index(['library_id', 'plate_name', 'donor_id', 'cell_type', 'sm_lincs_id', 'group'])","metadata":{"execution":{"iopub.status.busy":"2023-09-28T03:02:52.897257Z","iopub.execute_input":"2023-09-28T03:02:52.897829Z","iopub.status.idle":"2023-09-28T03:03:48.031559Z","shell.execute_reply.started":"2023-09-28T03:02:52.897789Z","shell.execute_reply":"2023-09-28T03:03:48.030247Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"pseudobulk","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"---","metadata":{}},{"cell_type":"markdown","source":"# DESeq2 (ongoing)","metadata":{}},{"cell_type":"markdown","source":"```python\nfrom pydeseq2.ds import DeseqStats, DeseqDataSet\n\ndds = DeseqDataSet(\n    counts=sum_count_df.reset_index().iloc[:, 5:],\n    metadata=pd.DataFrame(pseudobulk.index.tolist(), columns=pseudobulk.index.names),\n    design_factors=[\"library_id\", 'group'],\n#     design_factors=[\"library_id\", 'cell_type', 'sm_lincs_id'],\n#     design_factors=[\"library_id\", \"plate_name\", 'donor_id','cell_type', 'sm_lincs_id'],\n    refit_cooks=True,\n    n_cpus=12,\n)\n\ndds.deseq2()\n\n\nimport pickle as pkl\n\nwith open(\"dds.pkl\", \"wb\") as f:\n    pkl.dump(dds, f)\n```","metadata":{}},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"---\n---\n---","metadata":{}},{"cell_type":"markdown","source":"# For limma (draft code with disclaimer)","metadata":{}},{"cell_type":"code","source":"# rename for R matrix\n\ncell_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' : 'M',\n}\n\nsm_lincs_id = pseudobulk.reset_index('sm_lincs_id')['sm_lincs_id'].unique()\n\nsm_rename_dict = dict(\n    zip(\n        sm_lincs_id,\n        [sm_id.split(';')[0].replace('-','_') for sm_id in sm_lincs_id],\n    )\n)\n\npseudobulk_rename = (\n    pseudobulk\n    .reset_index()\n    .replace(sm_rename_dict)\n    .replace(cell_rename_dict)\n    .set_index(['library_id', 'plate_name', 'donor_id', 'cell_type', 'sm_lincs_id'])\n    .drop('group', axis=1)\n)","metadata":{"execution":{"iopub.status.busy":"2023-09-28T03:21:03.587757Z","iopub.execute_input":"2023-09-28T03:21:03.588248Z","iopub.status.idle":"2023-09-28T03:21:05.325144Z","shell.execute_reply.started":"2023-09-28T03:21:03.588215Z","shell.execute_reply":"2023-09-28T03:21:05.323827Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"---","metadata":{}},{"cell_type":"markdown","source":"## Some QC / Normalization / Correction (TBD)","metadata":{"execution":{"iopub.status.busy":"2023-09-25T07:32:12.258165Z","iopub.execute_input":"2023-09-25T07:32:12.259634Z","iopub.status.idle":"2023-09-25T07:32:12.265182Z","shell.execute_reply.started":"2023-09-25T07:32:12.259569Z","shell.execute_reply":"2023-09-25T07:32:12.263942Z"}}},{"cell_type":"code","source":"log_df = (\n    (\n        (\n            pseudobulk_rename.T \n            / \n            pseudobulk_rename.sum(axis=1).values\n        ).T\n    )*1e6\n).apply(lambda x: np.log10(x+1)).reset_index()\n\nlog_df['group'] = log_df.apply(lambda row: '{}x{}'.format(row['cell_type'], row['sm_lincs_id']), axis=1)\nlog_df = log_df.set_index(['library_id', 'plate_name', 'donor_id', 'cell_type', 'sm_lincs_id', 'group'])\nlog_df.to_csv('preprocessing.csv')","metadata":{"execution":{"iopub.status.busy":"2023-09-28T03:21:05.327753Z","iopub.execute_input":"2023-09-28T03:21:05.328256Z","iopub.status.idle":"2023-09-28T03:23:30.945266Z","shell.execute_reply.started":"2023-09-28T03:21:05.328215Z","shell.execute_reply":"2023-09-28T03:23:30.943758Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# genes_mean = pseudobulk.mean()\n# genes_std = pseudobulk.std()\n\ngenes_mean = log_df.mean()\ngenes_std = log_df.std()\n\n# LOWESS: locally weighted regression\n# gene_std_sm = sm.nonparametric.lowess(genes_std, genes_mean, return_sorted=False)","metadata":{"execution":{"iopub.status.busy":"2023-09-28T03:23:30.950134Z","iopub.execute_input":"2023-09-28T03:23:30.950534Z","iopub.status.idle":"2023-09-28T03:23:31.715751Z","shell.execute_reply.started":"2023-09-28T03:23:30.9505Z","shell.execute_reply":"2023-09-28T03:23:31.71467Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# plt.scatter(genes_mean, genes_std)\nplt.scatter(genes_mean, np.sqrt(genes_std))\n\n# plt.scatter(genes_mean, gene_std_sm)\n\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-09-28T03:23:31.718961Z","iopub.execute_input":"2023-09-28T03:23:31.71962Z","iopub.status.idle":"2023-09-28T03:23:32.466243Z","shell.execute_reply.started":"2023-09-28T03:23:31.719577Z","shell.execute_reply":"2023-09-28T03:23:32.465355Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"---","metadata":{}},{"cell_type":"markdown","source":"# limma\n### refs \n#### https://www.ncbi.nlm.nih.gov/pmc/articles/PMC7873980/\n#### https://bioconductor.org/packages/release/bioc/vignettes/limma/inst/doc/usersguide.pdf","metadata":{"execution":{"iopub.status.busy":"2023-09-25T08:11:44.734018Z","iopub.execute_input":"2023-09-25T08:11:44.734387Z","iopub.status.idle":"2023-09-25T08:11:44.739623Z","shell.execute_reply.started":"2023-09-25T08:11:44.734357Z","shell.execute_reply":"2023-09-25T08:11:44.73831Z"}}},{"cell_type":"code","source":"design_columns = ['intercept']\n\nfor level in (0, 1, 2, 5):\n    level_label = log_df.index.levels[level]\n    # others\n    if level in (0, 1, 2):  \n        level_label = sorted(level_label, key=lambda x: int(x.split('_')[1]))\n    elif level == 5:\n        level_label = sorted(level_label, key=lambda x: int(x.split('_')[1]))\n    else:\n        raise\n        \n    design_columns += level_label\n\n    \ndesign_matrix = pd.DataFrame(columns=design_columns, index=log_df.index)\n\nfor idx in log_df.index:\n    group_idx = pd.Series(idx)[[0, 1, 2, 5]].tolist()\n    design_matrix.loc[idx, group_idx] = 1\n    \ndesign_matrix['intercept'] = 1\ndesign_matrix = design_matrix.fillna(0)\n\n# save design matrix\ndesign_matrix.to_csv('design_matrix.csv', index=False)\ndesign_matrix","metadata":{"execution":{"iopub.status.busy":"2023-09-28T03:23:32.527077Z","iopub.execute_input":"2023-09-28T03:23:32.527745Z","iopub.status.idle":"2023-09-28T03:23:39.51937Z","shell.execute_reply.started":"2023-09-28T03:23:32.527708Z","shell.execute_reply":"2023-09-28T03:23:39.517763Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig, ax = plt.subplots(figsize=(12, 4))\nsns.heatmap(design_matrix, ax=ax, cbar=False)\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-09-28T03:01:08.85252Z","iopub.execute_input":"2023-09-28T03:01:08.852981Z","iopub.status.idle":"2023-09-28T03:01:11.631193Z","shell.execute_reply.started":"2023-09-28T03:01:08.852925Z","shell.execute_reply":"2023-09-28T03:01:11.629822Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"---","metadata":{}},{"cell_type":"code","source":"de_train_exp = (\n    pd.read_parquet('/kaggle/input/open-problems-single-cell-perturbations/de_train.parquet')\n    .groupby(['cell_type', 'sm_lincs_id'])\n    .size()\n    .reset_index()\n    .drop(0, axis=1)\n    .replace(sm_rename_dict)\n    .replace(cell_rename_dict)\n)\n\ngroup_name = de_train_exp.apply(lambda x: '{}x{}-{}x{}'.format(x['cell_type'], x['sm_lincs_id'], x['cell_type'], 'LSM_36361'), axis=1)\n# group_name = de_train_exp.apply(lambda x: 'group{}x{}-group{}x{}'.format(x['cell_type'], x['sm_lincs_id'], x['cell_type'], 'LSM_36361'), axis=1)\n\nde_train_exp['cont'] = group_name\n\ncontrast_matrix = pd.DataFrame(index=design_matrix.columns, columns=de_train_exp['cont']).fillna(0)\n\nfor col_idx in contrast_matrix.columns:\n    sample, neg_ctrl = col_idx.split('-')\n    contrast_matrix.loc[sample, col_idx] = 1\n    contrast_matrix.loc[neg_ctrl, col_idx] = -1\n\n    \n# save contrast matrix\ncontrast_matrix.to_csv('contrast_matrix.csv', index=False)\ncontrast_matrix","metadata":{"execution":{"iopub.status.busy":"2023-09-28T03:31:48.50227Z","iopub.execute_input":"2023-09-28T03:31:48.502762Z","iopub.status.idle":"2023-09-28T03:31:50.817327Z","shell.execute_reply.started":"2023-09-28T03:31:48.502726Z","shell.execute_reply":"2023-09-28T03:31:50.815903Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"---","metadata":{}},{"cell_type":"markdown","source":"# limma code (NOT Working !!!)","metadata":{}},{"cell_type":"code","source":"limma_code = '''\nlibrary(limma)\n\ndata <- read.csv('preprocessing.csv')\nvdata <- data[, 7:ncol(data)]\nvdata <- t(vdata)\n\nplate_name = factor(data$plate_name)\nlibrary_id = factor(data$library_id)\ndonor_id = factor(data$donor_id)\ncell_type = factor(data$cell_type)\nsm_lincs_id = factor(data$sm_lincs_id)\ngroup=factor(data$group)\n\ndesign = read.csv('design_matrix.csv')\n# design <- model.matrix(~plate_name+library_id+donor_id+group)\n# design <- model.matrix(~plate_name+library_id+donor_id+cell_type+sm_lincs_id)\n\nfit <- lmFit(vdata, design)\n'''","metadata":{"execution":{"iopub.status.busy":"2023-09-28T03:38:16.025141Z","iopub.execute_input":"2023-09-28T03:38:16.025632Z","iopub.status.idle":"2023-09-28T03:38:16.03231Z","shell.execute_reply.started":"2023-09-28T03:38:16.025598Z","shell.execute_reply":"2023-09-28T03:38:16.031002Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{"trusted":true},"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":[]}]}