{"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 --quiet tables\n\nimport os\nimport pandas as pd\nimport matplotlib.pyplot as plt\nimport numpy as np\n\nDATA_DIR = \"../input/open-problems-multimodal/\"\nFP_CELL_METADATA = os.path.join(DATA_DIR,\"metadata.csv\")\n\nFP_CITE_TRAIN_INPUTS = os.path.join(DATA_DIR,\"train_cite_inputs.h5\")\nFP_CITE_TRAIN_TARGETS = os.path.join(DATA_DIR,\"train_cite_targets.h5\")\nFP_CITE_TEST_INPUTS = os.path.join(DATA_DIR,\"test_cite_inputs.h5\")\n\nFP_MULTIOME_TRAIN_INPUTS = os.path.join(DATA_DIR,\"train_multi_inputs.h5\")\nFP_MULTIOME_TRAIN_TARGETS = os.path.join(DATA_DIR,\"train_multi_targets.h5\")\nFP_MULTIOME_TEST_INPUTS = os.path.join(DATA_DIR,\"test_multi_inputs.h5\")\n\nFP_SUBMISSION = os.path.join(DATA_DIR,\"sample_submission.csv\")\nFP_EVALUATION_IDS = os.path.join(DATA_DIR,\"evaluation_ids.csv\")","metadata":{"execution":{"iopub.status.busy":"2022-08-25T16:25:21.261208Z","iopub.execute_input":"2022-08-25T16:25:21.262007Z","iopub.status.idle":"2022-08-25T16:25:30.724140Z","shell.execute_reply.started":"2022-08-25T16:25:21.261965Z","shell.execute_reply":"2022-08-25T16:25:30.722981Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"START = int(1e4)\nSTOP = START+1000\n\ndf_multi_train_x = pd.read_hdf(FP_MULTIOME_TRAIN_INPUTS,start=START,stop=STOP)\ndf_multi_train_x.head()","metadata":{"execution":{"iopub.status.busy":"2022-08-25T16:25:30.726632Z","iopub.execute_input":"2022-08-25T16:25:30.727289Z","iopub.status.idle":"2022-08-25T16:25:34.895485Z","shell.execute_reply.started":"2022-08-25T16:25:30.727246Z","shell.execute_reply":"2022-08-25T16:25:34.894414Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(sorted(list({i[:i.find(':')] for i in df_multi_train_x.columns})))","metadata":{"execution":{"iopub.status.busy":"2022-08-25T16:25:34.897154Z","iopub.execute_input":"2022-08-25T16:25:34.897536Z","iopub.status.idle":"2022-08-25T16:25:34.979939Z","shell.execute_reply.started":"2022-08-25T16:25:34.897499Z","shell.execute_reply":"2022-08-25T16:25:34.978972Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(f\"Values range from {df_multi_train_x.min().min()} to {df_multi_train_x.max().max()}\")","metadata":{"execution":{"iopub.status.busy":"2022-08-25T16:25:34.982754Z","iopub.execute_input":"2022-08-25T16:25:34.983159Z","iopub.status.idle":"2022-08-25T16:25:37.056170Z","shell.execute_reply.started":"2022-08-25T16:25:34.983119Z","shell.execute_reply":"2022-08-25T16:25:37.055092Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"min_cells_non_zero = df_multi_train_x.gt(0).sum(axis=1).min()\nmax_cells_non_zero = df_multi_train_x.gt(0).sum(axis=1).max()\nsum_non_zero_values = df_multi_train_x.sum().sum()\ncount_non_zero_values = df_multi_train_x.gt(0).sum().sum()\naverage_non_zero_per_gene = df_multi_train_x[df_multi_train_x.gt(0)].count(axis = 1).mean()\n\nprint(f\"Each cell has at least {min_cells_non_zero} genes with non-zero accessibility values and a maximum of {max_cells_non_zero}.\")\nprint(f\"On average there are {round(average_non_zero_per_gene)} genes with non-zero accessibility values in each cell.\")\nprint(f\"The average non-zero value is about {sum_non_zero_values / count_non_zero_values:.2f}.\")\n\n# investigate standard deviation of features\nstd_dev_of_genes = df_multi_train_x.std()\n\n# ignore genes that are only accessible in a single cell\nstd_dev_of_genes_without_singles = std_dev_of_genes[df_multi_train_x.gt(0).sum().gt(1)]\nprint(f\"The standard deviation of our features is between {std_dev_of_genes_without_singles.min():.2f} and {std_dev_of_genes_without_singles.max():.2f}.\\nThe average standard deviation is {std_dev_of_genes_without_singles.mean():.2f}\")","metadata":{"execution":{"iopub.status.busy":"2022-08-25T16:25:37.057565Z","iopub.execute_input":"2022-08-25T16:25:37.058708Z","iopub.status.idle":"2022-08-25T16:25:47.381685Z","shell.execute_reply.started":"2022-08-25T16:25:37.058668Z","shell.execute_reply":"2022-08-25T16:25:47.380625Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"s = df_multi_train_x.gt(0).sum(axis = 1)\ncounts = s.groupby(lambda x: s[x] // 300).count()\ncounts.index = counts.index * 300\n\nfig, ax = plt.subplots()\nax.plot(counts.index, counts.values)\nax.set_xlabel('number of accessible genes')\nax.set_ylabel('number of cells')\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2022-08-25T16:25:47.385133Z","iopub.execute_input":"2022-08-25T16:25:47.386053Z","iopub.status.idle":"2022-08-25T16:25:48.543362Z","shell.execute_reply.started":"2022-08-25T16:25:47.386014Z","shell.execute_reply":"2022-08-25T16:25:48.542426Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_multi_train_y = pd.read_hdf(FP_MULTIOME_TRAIN_TARGETS, start=START, stop=STOP)\ndf_multi_train_y.head()","metadata":{"execution":{"iopub.status.busy":"2022-08-25T16:25:48.545043Z","iopub.execute_input":"2022-08-25T16:25:48.545386Z","iopub.status.idle":"2022-08-25T16:25:49.337051Z","shell.execute_reply.started":"2022-08-25T16:25:48.545350Z","shell.execute_reply":"2022-08-25T16:25:49.336091Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(sorted(list({i[:10] for i in df_multi_train_y.columns})))\nprint(df_multi_train_y.columns.str.len().unique().item())","metadata":{"execution":{"iopub.status.busy":"2022-08-25T16:25:49.338414Z","iopub.execute_input":"2022-08-25T16:25:49.338815Z","iopub.status.idle":"2022-08-25T16:25:49.361123Z","shell.execute_reply.started":"2022-08-25T16:25:49.338779Z","shell.execute_reply":"2022-08-25T16:25:49.360089Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(f\"Values for gene expression range from {df_multi_train_y.min().min():.2f} to {df_multi_train_y.max().max():.2f}\")\n\n# get data about non-zero values\nmin_cells_non_zero_y = df_multi_train_y.gt(0).sum(axis=1).min()\nmax_cells_non_zero_y = df_multi_train_y.gt(0).sum(axis=1).max()\nsum_non_zero_values_y = df_multi_train_y.sum().sum()\ncount_non_zero_values_y = df_multi_train_y.gt(0).sum().sum()\naverage_non_zero_per_gene_y = df_multi_train_y[df_multi_train_y.gt(0)].count(axis = 1).mean()\n\nprint(f\"Each cell has at least {min_cells_non_zero_y} genes with non-zero gene expression values and a maximum of {max_cells_non_zero_y}.\")\nprint(f\"On average there are {round(average_non_zero_per_gene_y)} genes with non-zero gene expression values in each cell.\")\nprint(f\"The average non-zero value for gene expression is about {sum_non_zero_values_y / count_non_zero_values_y:.2f}.\")\n\n# investigate standard deviation of features\nstd_dev_of_genes_y = df_multi_train_y.std()\n\n# ignore genes that are only accessible in a single cell\nstd_dev_of_genes_without_singles_y = std_dev_of_genes_y[df_multi_train_y.gt(0).sum().gt(1)]\nprint(f\"The standard deviation of gene expression values is between {std_dev_of_genes_without_singles_y.min():.2f} and {std_dev_of_genes_without_singles_y.max():.2f}.\\nThe average standard deviation is {std_dev_of_genes_without_singles_y.mean():.2f}\")","metadata":{"execution":{"iopub.status.busy":"2022-08-25T16:25:49.364435Z","iopub.execute_input":"2022-08-25T16:25:49.364736Z","iopub.status.idle":"2022-08-25T16:25:50.614207Z","shell.execute_reply.started":"2022-08-25T16:25:49.364709Z","shell.execute_reply":"2022-08-25T16:25:50.612783Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"s = df_multi_train_y.gt(0).sum(axis = 1)\ncounts = s.groupby(lambda x: s[x] // 100).count()\ncounts.index = counts.index * 100\n\nfig, ax = plt.subplots()\nax.plot(counts.index, counts.values)\nax.set_xlabel('number of genes expressed')\nax.set_ylabel('number of cells')\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2022-08-25T16:25:50.618156Z","iopub.execute_input":"2022-08-25T16:25:50.618461Z","iopub.status.idle":"2022-08-25T16:25:50.906864Z","shell.execute_reply.started":"2022-08-25T16:25:50.618417Z","shell.execute_reply":"2022-08-25T16:25:50.905944Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_cite_train_x = pd.read_hdf(FP_CITE_TRAIN_INPUTS,start=START,stop=STOP)\ndf_cite_train_x.head()","metadata":{"execution":{"iopub.status.busy":"2022-08-25T16:25:50.908295Z","iopub.execute_input":"2022-08-25T16:25:50.908662Z","iopub.status.idle":"2022-08-25T16:25:51.709464Z","shell.execute_reply.started":"2022-08-25T16:25:50.908627Z","shell.execute_reply":"2022-08-25T16:25:51.708450Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"gene_ids_multiome = set(df_multi_train_y.columns)\nprint(f\"Different Gene IDs in multiome: {len(gene_ids_multiome)}\")\n#for now we just keep the stem of the gene_id\ngene_ids_citeseq = set([i[:i.find(\"_\")] for i in df_cite_train_x.columns])\nprint(f\"Different Gene IDs in citeseq: {len(gene_ids_citeseq)}\")","metadata":{"execution":{"iopub.status.busy":"2022-08-25T16:25:51.710916Z","iopub.execute_input":"2022-08-25T16:25:51.711346Z","iopub.status.idle":"2022-08-25T16:25:51.743901Z","shell.execute_reply.started":"2022-08-25T16:25:51.711301Z","shell.execute_reply":"2022-08-25T16:25:51.742284Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(f\"Elements in Set Union: {len(gene_ids_citeseq | gene_ids_multiome)}\")\nprint(f\"Elements in Set Intersection: {len(gene_ids_citeseq & gene_ids_multiome)}\")\nprint(f\"multiome has {len(gene_ids_multiome - gene_ids_citeseq)} unique gene ids.\")\nprint(f\"Citeseq has {len(gene_ids_citeseq - gene_ids_multiome)} unique gene ids.\")","metadata":{"execution":{"iopub.status.busy":"2022-08-25T16:25:51.747015Z","iopub.execute_input":"2022-08-25T16:25:51.747782Z","iopub.status.idle":"2022-08-25T16:25:51.768350Z","shell.execute_reply.started":"2022-08-25T16:25:51.747737Z","shell.execute_reply":"2022-08-25T16:25:51.766361Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(f\"Values for gene expression range from {df_cite_train_x.min().min():.2f} to {df_cite_train_x.max().max():.2f}\")\n\n# get data about non-zero values\nmin_cells_non_zero_y = df_cite_train_x.gt(0).sum(axis=1).min()\nmax_cells_non_zero_y = df_cite_train_x.gt(0).sum(axis=1).max()\nsum_non_zero_values_y = df_cite_train_x.sum().sum()\ncount_non_zero_values_y = df_cite_train_x.gt(0).sum().sum()\naverage_non_zero_per_gene_y = df_cite_train_x[df_cite_train_x.gt(0)].count(axis = 1).mean()\n\nprint(f\"Each cell has at least {min_cells_non_zero_y} genes with non-zero gene expression values and a maximum of {max_cells_non_zero_y}.\")\nprint(f\"On average there are {round(average_non_zero_per_gene_y)} genes with non-zero gene expression values in each cell.\")\nprint(f\"The average non-zero value for gene expression is about {sum_non_zero_values_y / count_non_zero_values_y:.2f}.\")\n\n# investigate standard deviation of features\nstd_dev_of_genes_y = df_cite_train_x.std()\n\n# ignore genes that are only accessible in a single cell\nstd_dev_of_genes_without_singles_y = std_dev_of_genes_y[df_cite_train_x.gt(0).sum().gt(1)]\nprint(f\"The standard deviation of gene expression values is between {std_dev_of_genes_without_singles_y.min():.2f} and {std_dev_of_genes_without_singles_y.max():.2f}.\\nThe average standard deviation is {std_dev_of_genes_without_singles_y.mean():.2f}\")","metadata":{"execution":{"iopub.status.busy":"2022-08-25T16:25:51.770311Z","iopub.execute_input":"2022-08-25T16:25:51.771663Z","iopub.status.idle":"2022-08-25T16:25:52.943995Z","shell.execute_reply.started":"2022-08-25T16:25:51.771621Z","shell.execute_reply":"2022-08-25T16:25:52.942723Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"s = df_cite_train_x.gt(0).sum(axis = 1)\ncounts = s.groupby(lambda x: s[x] // 300).count()\ncounts.index = counts.index * 300\n\nfig, ax = plt.subplots()\nax.plot(counts.index, counts.values)\nax.set_xlabel('number of genes expressed')\nax.set_ylabel('number of cells')\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2022-08-25T16:25:52.945521Z","iopub.execute_input":"2022-08-25T16:25:52.946053Z","iopub.status.idle":"2022-08-25T16:25:53.238465Z","shell.execute_reply.started":"2022-08-25T16:25:52.946013Z","shell.execute_reply":"2022-08-25T16:25:53.237389Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_cite_train_y = pd.read_hdf(FP_CITE_TRAIN_TARGETS)\ndf_cite_train_y.head()","metadata":{"execution":{"iopub.status.busy":"2022-08-25T16:25:53.239768Z","iopub.execute_input":"2022-08-25T16:25:53.240367Z","iopub.status.idle":"2022-08-25T16:25:53.971605Z","shell.execute_reply.started":"2022-08-25T16:25:53.240326Z","shell.execute_reply":"2022-08-25T16:25:53.970475Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(f\"Measurements of surface protein levels range from {df_cite_train_y.min().min():.2f} to {df_cite_train_y.max().max():.2f}.\")\nprint(f\"The average value is {df_cite_train_y.mean().mean():.2f}.\")\nprint(f\"The standard deviation of surface protein levels is between {df_cite_train_y.std().min():.2f} and {df_cite_train_y.std().max():.2f}.\")\nprint(f\"The average standard deviation is {df_cite_train_y.std().mean():.2f}.\")","metadata":{"execution":{"iopub.status.busy":"2022-08-25T16:25:53.973212Z","iopub.execute_input":"2022-08-25T16:25:53.973609Z","iopub.status.idle":"2022-08-25T16:25:54.310928Z","shell.execute_reply.started":"2022-08-25T16:25:53.973573Z","shell.execute_reply":"2022-08-25T16:25:54.309636Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"threshold = 0.1\ndf_cite_train_y.applymap(lambda x: abs(x)).gt(threshold).sum(axis = 1)\nprint(f\"Each cell has between {df_cite_train_y.applymap(lambda x: abs(x)).gt(threshold).sum(axis = 1).min()} and {df_cite_train_y.applymap(lambda x: abs(x)).gt(threshold).sum(axis = 1).max()} measurements with absolute values over {threshold}.\")","metadata":{"execution":{"iopub.status.busy":"2022-08-25T16:25:54.312368Z","iopub.execute_input":"2022-08-25T16:25:54.312863Z","iopub.status.idle":"2022-08-25T16:26:03.381903Z","shell.execute_reply.started":"2022-08-25T16:25:54.312823Z","shell.execute_reply":"2022-08-25T16:26:03.380840Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_meta = pd.read_csv(FP_CELL_METADATA).set_index(\"cell_id\")\ndf_meta","metadata":{"execution":{"iopub.status.busy":"2022-08-25T16:26:03.383269Z","iopub.execute_input":"2022-08-25T16:26:03.384054Z","iopub.status.idle":"2022-08-25T16:26:03.620432Z","shell.execute_reply.started":"2022-08-25T16:26:03.384016Z","shell.execute_reply":"2022-08-25T16:26:03.619407Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig, axarr = plt.subplots(nrows=1, ncols=3, figsize=(12, 5))\nfor i, col in enumerate([\"donor\", \"day\", \"technology\"]):\n    _= df_meta[[col]].value_counts().plot.pie(ax=axarr[i], autopct='%1.1f%%', ylabel=col)","metadata":{"execution":{"iopub.status.busy":"2022-08-25T16:26:03.621753Z","iopub.execute_input":"2022-08-25T16:26:03.622601Z","iopub.status.idle":"2022-08-25T16:26:03.930659Z","shell.execute_reply.started":"2022-08-25T16:26:03.622563Z","shell.execute_reply":"2022-08-25T16:26:03.929662Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}