{"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":"%%capture\n!pip3 install scanpy[leiden]","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2022-12-18T03:18:01.582287Z","iopub.execute_input":"2022-12-18T03:18:01.583164Z","iopub.status.idle":"2022-12-18T03:18:18.345167Z","shell.execute_reply.started":"2022-12-18T03:18:01.583065Z","shell.execute_reply":"2022-12-18T03:18:18.343604Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import gc\nimport os\nimport sys\n\nimport anndata as ad\nimport matplotlib.image as mpimg\nimport matplotlib.pyplot as plt\nimport numpy as np\nimport pandas as pd\nimport scipy.stats as stats\nimport seaborn as sns\nimport scanpy as sc\nfrom tqdm.notebook import tqdm\n\n\nsc.settings.verbosity = 3\nsc.set_figure_params(dpi=150)\nsns.set_style(\"ticks\")\n\nDATA_DIR = \"/kaggle/input/openproblemsmultimodalanndata\"","metadata":{"execution":{"iopub.status.busy":"2022-12-18T03:18:18.347893Z","iopub.execute_input":"2022-12-18T03:18:18.348419Z","iopub.status.idle":"2022-12-18T03:18:20.594968Z","shell.execute_reply.started":"2022-12-18T03:18:18.348371Z","shell.execute_reply":"2022-12-18T03:18:20.593753Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Cite\nFP_CITE_TRAIN_INPUTS_ann = os.path.join(DATA_DIR,\"train_cite_inputs_values.sparse.ann.h5\")\nFP_CITE_TRAIN_TARGETS_ann = os.path.join(DATA_DIR,\"train_cite_targets_values.sparse.ann.h5\")\nFP_CITE_TEST_INPUTS_ann = os.path.join(DATA_DIR,\"test_cite_inputs_values.sparse.ann.h5\")\n\nDATA_DIR_RAW = \"/kaggle/input/msci-raw-counts-anndata\"\n# Cite\nFP_CITE_RAW_TRAIN_INPUTS = os.path.join(DATA_DIR_RAW,\"train_cite_inputs_raw.ann.h5\")\nFP_CITE_RAW_TRAIN_TARGETS = os.path.join(DATA_DIR_RAW,\"train_cite_targets_raw.ann.h5\")\nFP_CITE_RAW_TEST_INPUTS = os.path.join(DATA_DIR_RAW,\"test_cite_inputs_raw.ann.h5\")","metadata":{"execution":{"iopub.status.busy":"2022-12-18T03:18:20.596502Z","iopub.execute_input":"2022-12-18T03:18:20.596816Z","iopub.status.idle":"2022-12-18T03:18:20.603083Z","shell.execute_reply.started":"2022-12-18T03:18:20.596787Z","shell.execute_reply":"2022-12-18T03:18:20.602005Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"metadata_fp = \"/kaggle/input/open-problems-multimodal/metadata.csv\"\nadd_metadata_fp = \"/kaggle/input/open-problems-multimodal/metadata_cite_day_2_donor_27678.csv\"\nmetadata = pd.read_csv(metadata_fp)\nadd_metadata = pd.read_csv(add_metadata_fp)\nmetadata = pd.concat([metadata, add_metadata])\nmetadata.set_index(\"cell_id\", inplace=True)\nmetadata.head()","metadata":{"execution":{"iopub.status.busy":"2022-12-18T03:18:20.605439Z","iopub.execute_input":"2022-12-18T03:18:20.605826Z","iopub.status.idle":"2022-12-18T03:18:21.085596Z","shell.execute_reply.started":"2022-12-18T03:18:20.605796Z","shell.execute_reply":"2022-12-18T03:18:21.084613Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def annotate(adata: ad.AnnData, metadata: pd.DataFrame, f_donors: [str]):\n    adata.obs[\"cell_type\"] = metadata[\"cell_type\"][adata.obs_names]\n    adata.obs[\"day\"] = metadata[\"day\"][adata.obs_names].astype(str)\n    adata.obs[\"donor\"] = metadata[\"donor\"][adata.obs_names].astype(str)\n    adata.obs[\"gender\"] = \"m\"\n    adata.obs.loc[adata.obs[\"donor\"].isin(f_donors), \"gender\"] = \"f\"\n    adata.obs[\"n_genes_counts\"] = (adata.X > 0).sum(axis=1).A1\n    adata.obs[\"read_depth\"] = adata.X.sum(axis=1).A1","metadata":{"execution":{"iopub.status.busy":"2022-12-18T03:18:21.087164Z","iopub.execute_input":"2022-12-18T03:18:21.087606Z","iopub.status.idle":"2022-12-18T03:18:21.095212Z","shell.execute_reply.started":"2022-12-18T03:18:21.087559Z","shell.execute_reply":"2022-12-18T03:18:21.093918Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"cite_target = sc.read_h5ad(FP_CITE_TRAIN_TARGETS_ann)\nannotate(cite_target, metadata=metadata, f_donors=[\"13176\", ])\ncite_target","metadata":{"execution":{"iopub.status.busy":"2022-12-18T03:18:21.096575Z","iopub.execute_input":"2022-12-18T03:18:21.096867Z","iopub.status.idle":"2022-12-18T03:18:22.756663Z","shell.execute_reply.started":"2022-12-18T03:18:21.096841Z","shell.execute_reply":"2022-12-18T03:18:22.755596Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"cite_target_df = pd.DataFrame(cite_target.X.A, columns=cite_target.var_names, index=cite_target.obs_names)\ncite_target_df = pd.concat([cite_target_df, cite_target.obs], axis=1)\ncite_target_df.head()","metadata":{"execution":{"iopub.status.busy":"2022-12-18T03:18:22.757924Z","iopub.execute_input":"2022-12-18T03:18:22.758282Z","iopub.status.idle":"2022-12-18T03:18:22.947531Z","shell.execute_reply.started":"2022-12-18T03:18:22.758251Z","shell.execute_reply":"2022-12-18T03:18:22.946364Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"cite = sc.read_h5ad(FP_CITE_RAW_TRAIN_INPUTS)\nannotate(cite, metadata, f_donors = [\"13176\", ])","metadata":{"execution":{"iopub.status.busy":"2022-12-18T03:18:22.948880Z","iopub.execute_input":"2022-12-18T03:18:22.949255Z","iopub.status.idle":"2022-12-18T03:18:55.304293Z","shell.execute_reply.started":"2022-12-18T03:18:22.949222Z","shell.execute_reply":"2022-12-18T03:18:55.303200Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sc.pp.normalize_total(cite)","metadata":{"execution":{"iopub.status.busy":"2022-12-18T03:18:55.305859Z","iopub.execute_input":"2022-12-18T03:18:55.306288Z","iopub.status.idle":"2022-12-18T03:18:57.326311Z","shell.execute_reply.started":"2022-12-18T03:18:55.306259Z","shell.execute_reply":"2022-12-18T03:18:57.325262Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"CD_info = pd.read_csv(\"/kaggle/input/research-project-01-around-multimodal-singlecell/CD_info_alt.csv\", index_col=0).sort_values(\"cd_name\")\nCD_info.head()","metadata":{"execution":{"iopub.status.busy":"2022-12-18T03:18:57.330569Z","iopub.execute_input":"2022-12-18T03:18:57.331397Z","iopub.status.idle":"2022-12-18T03:18:57.355721Z","shell.execute_reply.started":"2022-12-18T03:18:57.331349Z","shell.execute_reply":"2022-12-18T03:18:57.354473Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"CD_info.shape[0]","metadata":{"execution":{"iopub.status.busy":"2022-12-18T03:18:57.356927Z","iopub.execute_input":"2022-12-18T03:18:57.357387Z","iopub.status.idle":"2022-12-18T03:18:57.364647Z","shell.execute_reply.started":"2022-12-18T03:18:57.357358Z","shell.execute_reply":"2022-12-18T03:18:57.363559Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# jointplot on subplots\nfrom https://stackoverflow.com/a/61330067","metadata":{}},{"cell_type":"code","source":"for i in tqdm(range(CD_info.shape[0])):\n    \n    cd = CD_info.iloc[i,:].cd_name\n    rna_id = CD_info.iloc[i,:].cite_input\n    \n    rna_counts = cite[:, rna_id].X.A.T[0]\n    if any(rna_counts > 0):\n        rna_counts /= (rna_counts.sum() / sum(rna_counts > 0))**0.5\n    cd_counts = cite_target_df.loc[:, cd]\n    if any(cd_counts > 0):\n        cd_counts /= (cd_counts.sum() / len(cd_counts))**0.5\n    \n    cd_df = pd.DataFrame({rna_id: rna_counts, cd: cite_target_df.loc[:, cd], \"cell_type\": cite.obs[\"cell_type\"]})\n    \n    g = sns.jointplot(data=cd_df, x=cd, y=rna_id, hue=\"cell_type\")\n    plt.suptitle(f\"{cd}, non-zero/zero: {sum(rna_counts > 0)}/{sum(rna_counts == 0)}\")\n\n    g.savefig(f'{cd}.png')\n    plt.close(g.fig)","metadata":{"execution":{"iopub.status.busy":"2022-12-18T03:28:46.173643Z","iopub.execute_input":"2022-12-18T03:28:46.174079Z","iopub.status.idle":"2022-12-18T03:39:09.758844Z","shell.execute_reply.started":"2022-12-18T03:28:46.174017Z","shell.execute_reply":"2022-12-18T03:39:09.757741Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"cols = 9\nrows = CD_info.shape[0] // cols + int((CD_info.shape[0] % cols) > 0)\n\nfig, axs = plt.subplots(rows, cols, figsize=(4 * cols, 4 * rows))\n\nfor i in tqdm(range(CD_info.shape[0])):\n    j = i // cols\n    k = i % cols\n    cd = CD_info.iloc[i,:].cd_name\n    axs[j, k].imshow(mpimg.imread(f'{cd}.png'))\n\n# turn off x and y axis\n[ax.set_axis_off() for ax in axs.ravel()]\n\nplt.tight_layout()","metadata":{"execution":{"iopub.status.busy":"2022-12-18T03:39:09.760869Z","iopub.execute_input":"2022-12-18T03:39:09.761227Z","iopub.status.idle":"2022-12-18T03:40:23.407362Z","shell.execute_reply.started":"2022-12-18T03:39:09.761195Z","shell.execute_reply":"2022-12-18T03:40:23.405259Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Variance RNA vs variance protein","metadata":{}},{"cell_type":"code","source":"rna_counts = cite[:, CD_info.cite_input].X.A\nmean_rna = rna_counts.mean(axis=0)\nvar_rna = np.var(rna_counts, axis=0)\n\ncd_counts = cite_target[:, CD_info.cd_name].X.A\nmean_cd = cd_counts.mean(axis=0)\nvar_cd = np.var(cd_counts, axis=0)","metadata":{"execution":{"iopub.status.busy":"2022-12-18T03:40:23.409366Z","iopub.execute_input":"2022-12-18T03:40:23.409717Z","iopub.status.idle":"2022-12-18T03:40:24.574409Z","shell.execute_reply.started":"2022-12-18T03:40:23.409686Z","shell.execute_reply":"2022-12-18T03:40:24.572959Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sns.scatterplot(x=var_cd, y=var_rna)\nplt.suptitle(\"Var(RNA) vs var(CD)\")\nplt.xlabel(\"var(CD)\")\nplt.ylabel(\"var(RNA)\")","metadata":{"execution":{"iopub.status.busy":"2022-12-18T03:40:24.576808Z","iopub.execute_input":"2022-12-18T03:40:24.577144Z","iopub.status.idle":"2022-12-18T03:40:24.926188Z","shell.execute_reply.started":"2022-12-18T03:40:24.577114Z","shell.execute_reply":"2022-12-18T03:40:24.925090Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sns.scatterplot(x=mean_rna, y=var_rna)\nplt.suptitle(\"Var(RNA) vs mean(RNA)\")\nplt.xlabel(\"mean(RNA)\")\nplt.ylabel(\"var(RNA)\")","metadata":{"execution":{"iopub.status.busy":"2022-12-18T03:40:24.927616Z","iopub.execute_input":"2022-12-18T03:40:24.927981Z","iopub.status.idle":"2022-12-18T03:40:25.247353Z","shell.execute_reply.started":"2022-12-18T03:40:24.927920Z","shell.execute_reply":"2022-12-18T03:40:25.245697Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sns.scatterplot(x=mean_cd, y=var_cd)\nplt.suptitle(\"Var(CD) vs mean(CD)\")\nplt.xlabel(\"mean(CD)\")\nplt.ylabel(\"var(CD)\")","metadata":{"execution":{"iopub.status.busy":"2022-12-18T03:40:25.249373Z","iopub.execute_input":"2022-12-18T03:40:25.249851Z","iopub.status.idle":"2022-12-18T03:40:25.682109Z","shell.execute_reply.started":"2022-12-18T03:40:25.249806Z","shell.execute_reply":"2022-12-18T03:40:25.680969Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"zero = (mean_cd == 0) | (mean_rna == 0)\ncd_ratio = var_cd[~zero] / mean_cd[~zero]\nrna_ratio = var_rna[~zero] / (mean_rna[~zero] + 0.14*mean_rna[~zero]**2)\n\nsns.scatterplot(x=cd_ratio, y=rna_ratio)\n# plt.suptitle(\"$Var(RNA)/(mean(RNA)+0.14*mean(RNA)^2) vs Var(CD)/mean(CD)$\")\nplt.xlabel(\"Var(CD)/mean(CD)\")\nplt.ylabel(\"$Var(RNA)/(mean(RNA)+0.14*mean(RNA)^2)$\")","metadata":{"execution":{"iopub.status.busy":"2022-12-18T03:40:25.683809Z","iopub.execute_input":"2022-12-18T03:40:25.684462Z","iopub.status.idle":"2022-12-18T03:40:26.327409Z","shell.execute_reply.started":"2022-12-18T03:40:25.684424Z","shell.execute_reply":"2022-12-18T03:40:26.326312Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"CD_info.cd_name.loc[~zero][np.argsort(-rna_ratio)][:10].values","metadata":{"execution":{"iopub.status.busy":"2022-12-18T03:40:26.328654Z","iopub.execute_input":"2022-12-18T03:40:26.328973Z","iopub.status.idle":"2022-12-18T03:40:26.338004Z","shell.execute_reply.started":"2022-12-18T03:40:26.328922Z","shell.execute_reply":"2022-12-18T03:40:26.336754Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"rna_ratio[np.argsort(-rna_ratio)][:10]","metadata":{"execution":{"iopub.status.busy":"2022-12-18T03:40:26.339646Z","iopub.execute_input":"2022-12-18T03:40:26.339988Z","iopub.status.idle":"2022-12-18T03:40:26.350300Z","shell.execute_reply.started":"2022-12-18T03:40:26.339949Z","shell.execute_reply":"2022-12-18T03:40:26.349211Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"CD_info.cd_name.loc[~zero][np.argsort(-cd_ratio)][:10].values","metadata":{"execution":{"iopub.status.busy":"2022-12-18T03:40:26.353101Z","iopub.execute_input":"2022-12-18T03:40:26.353402Z","iopub.status.idle":"2022-12-18T03:40:26.361633Z","shell.execute_reply.started":"2022-12-18T03:40:26.353375Z","shell.execute_reply":"2022-12-18T03:40:26.360491Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# draw plts only for suspicious CDs","metadata":{}},{"cell_type":"code","source":"cds_suspicious = CD_info.cd_name.loc[~zero][np.argsort(-rna_ratio)][:12].values\n\ncols = 4\nrows = len(cds_suspicious) // cols + int((len(cds_suspicious) % cols) > 0)\n\nfig, axs = plt.subplots(rows, cols, figsize=(4 * cols, 4 * rows))\n\nfor i in tqdm(range(len(cds_suspicious))):\n    j = i // cols\n    k = i % cols\n    cd = cds_suspicious[i]\n    axs[j, k].imshow(mpimg.imread(f'{cd}.png'))\n\n# turn off x and y axis\n[ax.set_axis_off() for ax in axs.ravel()]\n\nplt.tight_layout()","metadata":{"execution":{"iopub.status.busy":"2022-12-18T03:41:10.234278Z","iopub.execute_input":"2022-12-18T03:41:10.234999Z","iopub.status.idle":"2022-12-18T03:41:17.863903Z","shell.execute_reply.started":"2022-12-18T03:41:10.234962Z","shell.execute_reply":"2022-12-18T03:41:17.862910Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"cds_suspicious = CD_info.cd_name.loc[~zero][np.argsort(-cd_ratio)][:12].values\n\ncols = 4\nrows = len(cds_suspicious) // cols + int((len(cds_suspicious) % cols) > 0)\n\nfig, axs = plt.subplots(rows, cols, figsize=(4 * cols, 4 * rows))\n\nfor i in tqdm(range(len(cds_suspicious))):\n    j = i // cols\n    k = i % cols\n    cd = cds_suspicious[i]\n    axs[j, k].imshow(mpimg.imread(f'{cd}.png'))\n\n# turn off x and y axis\n[ax.set_axis_off() for ax in axs.ravel()]\n\nplt.tight_layout()","metadata":{"execution":{"iopub.status.busy":"2022-12-18T03:42:18.989744Z","iopub.execute_input":"2022-12-18T03:42:18.990801Z","iopub.status.idle":"2022-12-18T03:42:26.638749Z","shell.execute_reply.started":"2022-12-18T03:42:18.990750Z","shell.execute_reply":"2022-12-18T03:42:26.637345Z"},"trusted":true},"execution_count":null,"outputs":[]}]}