{"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":"# 1. Overview\nThe goal of this competition is to better understand the relationship between different modalities in cells. The goal of this notebook is to gain a better understanding of the associated data. This equips us with the knowledge needed to make good decisions about model design and data layout.\n\n**This is a work in progress. If any aspect needs clarification, please let me know. My understanding of genetics is very limited. Feel free to point out anything that is false.**\n\n
\n

1.1 What do we want to learn?

\n
\n\nDuring transcription in cells, there is a known flow of information. DNA must be accessible to produce RNA. Produced RNA is used as a template to build proteins. Therefore, one could assume that we can use knowledge about the accessibility of DNA to predict future states of RNA and that we could use knowledge about RNA to predict the concentration of proteins in the future. In this challenge, we want to learn more about this relationship between DNA, RNA, and proteins. We thus need to capture information about three distinct properties of a cell:\n* chromatin accessibility\n* gene expression\n* surface protein levels\n\n
\n

1.2 How are those three properties of a cell presented?

\n
\n\nBefore we have a look at how the information about those properties of a cell is laid out, we must note that the methods used to obtain the data do not capture all properties at once. We have two distinct methods for testing. The first one is the \"10x Chromium Single Cell Multiome ATAC + Gene Expression\" short \"multiome\" test. The second one is the \"10x Genomics Single Cell Gene Expression with Feature Barcoding technology\" short \"citeseq\" test.\n\nWith the multiome test, we can measure **chromatin accessibility and gene expression**. With the citeseq test, we can measure **gene expression and surface protein levels**.\n\nTherefore, we will have data about chromatin accessibility and surface protein levels once (from multiome and citeseq, respectively). And we will have data about gene expression two times, one from each test. With that out of the way, let's dive into how the data is actually presented.\n\n
\n

1.3 Imports

\n
","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19"}},{"cell_type":"code","source":"# installs\n!pip install --quiet tables\n\n# imports\nimport os\nimport pandas as pd\nimport matplotlib.pyplot as plt\nimport numpy as np\n\n\n# set paths\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":{"_kg_hide-input":true,"_kg_hide-output":true,"execution":{"iopub.status.busy":"2022-08-19T09:48:52.243766Z","iopub.execute_input":"2022-08-19T09:48:52.244313Z","iopub.status.idle":"2022-08-19T09:49:07.053007Z","shell.execute_reply.started":"2022-08-19T09:48:52.244215Z","shell.execute_reply":"2022-08-19T09:49:07.051547Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# 2. Data\n\n
\n

2.1 Chromatin accessibility data

\n
\n\nFirst we will have a look at the data about chromatin accessibility. Inspecting the corresponding HDF5 file, we see that the data is stored in one single table having dimensions (228942, 105942). Each value is stored as a 32bit float. Thus, to load the full table into memory while preserving the 32bit accuracy, we will need about 90 GB of RAM. This is quite a lot. Therefore, we will only look at a chunk of the data here and also don't load the whole dataset into memory while training or doing transformations. The other HDF5 tables also store values as 32bit floats.","metadata":{}},{"cell_type":"code","source":"# Loading the whole dataset into pandas exceeds the memory,\n# therefore we define start and stop values\nSTART = 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-19T09:49:07.056483Z","iopub.execute_input":"2022-08-19T09:49:07.057004Z","iopub.status.idle":"2022-08-19T09:49:11.568880Z","shell.execute_reply.started":"2022-08-19T09:49:07.056954Z","shell.execute_reply":"2022-08-19T09:49:11.567644Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"As we can see, each individual cell is identified by a cell_id. We then have 228942 columns that are named something like \"STUFF:NUMBER-NUMBER\". STUFF is actually the name of a chromosome, while the numbers are a range indicating where the gene starts and ends. Let's have a look at what kind of chromosomes we have:","metadata":{}},{"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-19T09:49:11.570459Z","iopub.execute_input":"2022-08-19T09:49:11.570846Z","iopub.status.idle":"2022-08-19T09:49:11.663303Z","shell.execute_reply.started":"2022-08-19T09:49:11.570814Z","shell.execute_reply":"2022-08-19T09:49:11.662053Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We actually find the chromosomes we expect, namely chr1-chr22, the 22 chromosomes humans have (called autosomes), and also chrX and chrY, being the gender-specific chromosomes. What about the ones starting with KI and GL? According to a quick internet search, those are unplaced genes. They most likely are part of the human genome, but we don't know yet on which chromosome they are. Noteworthy at this point is that the number of protein-coding genes in humans is estimated to be between 19.9k and 21.3k. Therefore, it looks like we have measurements of much more than just the protein-coding genes.\n\nNext, we check the range of the values we have.","metadata":{}},{"cell_type":"code","source":"# first call to min/max gives us the min/max in each column. \n# Than we min/max again to get total min/max\nprint(f\"Values range from {df_multi_train_x.min().min()} to {df_multi_train_x.max().max()}\")","metadata":{"execution":{"iopub.status.busy":"2022-08-19T09:49:11.664716Z","iopub.execute_input":"2022-08-19T09:49:11.665069Z","iopub.status.idle":"2022-08-19T09:49:13.751572Z","shell.execute_reply.started":"2022-08-19T09:49:11.665037Z","shell.execute_reply":"2022-08-19T09:49:13.750294Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"So let's summarize what we have learned about the data corresponding to the accessibility of DNA so far:\n\n* We have chromatin accessibility measurements for approximately 106k cells in total.\n* We measure how accessible certain genes are in each cell, approximately 229k genes per cell.\n* Accessibility is given in numbers from 0.0 to ~18. We do not know the upper bound because we have not looked at all the data yet.\n* The values in our dataset use 32bit precision floats","metadata":{}},{"cell_type":"markdown","source":"What else do we want to know about chromatin accessibility data?\n\n* How many values are non-zero for each cell?\n* What is the standard deviation and what is the average non-zero value?\n\nFirst, we will have a closer look at the chromatin accessibility values of each cell.","metadata":{}},{"cell_type":"code","source":"# get data about non-zero values\nmin_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-19T09:49:13.754451Z","iopub.execute_input":"2022-08-19T09:49:13.754822Z","iopub.status.idle":"2022-08-19T09:49:24.174983Z","shell.execute_reply.started":"2022-08-19T09:49:13.754789Z","shell.execute_reply":"2022-08-19T09:49:24.173728Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"That's already good information about what we can expect from our features for the first problem. To even better understand how many features we have for each sample, we will plot the number of cells per feature count.","metadata":{}},{"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-19T09:52:42.303926Z","iopub.execute_input":"2022-08-19T09:52:42.304886Z","iopub.status.idle":"2022-08-19T09:52:43.407985Z","shell.execute_reply.started":"2022-08-19T09:52:42.304843Z","shell.execute_reply":"2022-08-19T09:52:43.407036Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"As we can see, the majority of our cells have between 2K and 7K accessible genes.\n\nWe now have quite a good understanding of the chromatin accessibility measured with the multiome test. We continue with investigating gene expression features.","metadata":{}},{"cell_type":"markdown","source":"
\n

2.2 Gene expression data

\n
\n\nAs mentioned before we have two datasets containing gene expression data. We will first look at the data from the multiome test.\n\n### 2.2.1 Gene expression from multiome\nWe would actually be able to load the whole dataset at once, but we will only look at the part corresponding to the already seen X values for now.","metadata":{}},{"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-19T09:49:25.381363Z","iopub.execute_input":"2022-08-19T09:49:25.382142Z","iopub.status.idle":"2022-08-19T09:49:26.150298Z","shell.execute_reply.started":"2022-08-19T09:49:25.382094Z","shell.execute_reply":"2022-08-19T09:49:26.149145Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"As we can see, we have 23418 values that our model will need to predict. But what exactly are those values?","metadata":{}},{"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-19T09:55:20.700649Z","iopub.execute_input":"2022-08-19T09:55:20.701184Z","iopub.status.idle":"2022-08-19T09:55:20.726863Z","shell.execute_reply.started":"2022-08-19T09:55:20.701090Z","shell.execute_reply":"2022-08-19T09:55:20.725406Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"As we can see, all of the features start with ENSG and then 5 zeroes. What we have here is called the Ensambl ID. The general form is ENS(species)(object type)(identifier).(version).\nENS tells us that we are looking at an ensembl ID. The species part is empty for human genes by convention. The object type is G for gene. It looks like the identifier is always 11 decimals long. And it looks like we don't have any version specifications in our data.\n\nWe will now check for similar properties than before.","metadata":{}},{"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-19T09:49:26.183500Z","iopub.execute_input":"2022-08-19T09:49:26.183837Z","iopub.status.idle":"2022-08-19T09:49:27.473562Z","shell.execute_reply.started":"2022-08-19T09:49:26.183798Z","shell.execute_reply":"2022-08-19T09:49:27.472263Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We can see that the range of gene expression values is smaller than that for chromatin accessibility, but the standard deviation is higher. This might be important for the design of our model.\n\nEven though this information will probably not influence the design of our model, let's still have a look at how many genes are expressed in cells. Just because it is interesting.","metadata":{}},{"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-19T09:49:27.474977Z","iopub.execute_input":"2022-08-19T09:49:27.475333Z","iopub.status.idle":"2022-08-19T09:49:27.761735Z","shell.execute_reply.started":"2022-08-19T09:49:27.475294Z","shell.execute_reply":"2022-08-19T09:49:27.760660Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Having an overview of gene expression data obtained by the multiome test, we will now compare it to those obtained by the citeseq test.\n\n### 2.2.2 Gene expression from citeseq","metadata":{}},{"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-19T09:49:27.763196Z","iopub.execute_input":"2022-08-19T09:49:27.763536Z","iopub.status.idle":"2022-08-19T09:49:28.602416Z","shell.execute_reply.started":"2022-08-19T09:49:27.763506Z","shell.execute_reply":"2022-08-19T09:49:28.601196Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The first thing we notice is that the start of our gene_id looks much like what we have seen in the multiome data, but there is a new suffix. So what is it about the suffix?\n\nChecking the Ensembl ID of gene_id on [ensembl.org](https://www.ensembl.org/Homo_sapiens/Gene/Summary?db=core;g=ENSG00000121410;r=19:58345178-58353492) (in this case for ENSG00000121410) we see that the suffix is actually the name of the gene. As we will see in the next code cell, the gene_id is unique even without this suffix, so it looks like redundant information for now.","metadata":{}},{"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-19T09:49:28.603997Z","iopub.execute_input":"2022-08-19T09:49:28.604719Z","iopub.status.idle":"2022-08-19T09:49:28.626628Z","shell.execute_reply.started":"2022-08-19T09:49:28.604675Z","shell.execute_reply":"2022-08-19T09:49:28.625369Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"As mentioned, stripping of the suffix still produces unique names.\n\nLet's check for overlap in both datasets about gene expression.","metadata":{}},{"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-19T09:49:28.628199Z","iopub.execute_input":"2022-08-19T09:49:28.629094Z","iopub.status.idle":"2022-08-19T09:49:28.649523Z","shell.execute_reply.started":"2022-08-19T09:49:28.629048Z","shell.execute_reply":"2022-08-19T09:49:28.648070Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Even though we have a huge intersection, there are quite a few genes unique to each test.\n\nWe will now again get information about distribution of values:","metadata":{}},{"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-19T09:49:28.653094Z","iopub.execute_input":"2022-08-19T09:49:28.653739Z","iopub.status.idle":"2022-08-19T09:49:29.861258Z","shell.execute_reply.started":"2022-08-19T09:49:28.653703Z","shell.execute_reply":"2022-08-19T09:49:29.860009Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We can see that the gene expression features obtained by the citeseq and the multiome test are quite similar. Values are in the same range and also the standard deviation and average non-zero values are of comparable size. For now, we only had a look at part of the data, so it will be interesting to see if this holds for all the data. We will get to that later.\n\nFor now, let's have a look at how many genes are expressed in individual cells.","metadata":{}},{"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-19T09:49:29.862598Z","iopub.execute_input":"2022-08-19T09:49:29.862935Z","iopub.status.idle":"2022-08-19T09:49:30.159536Z","shell.execute_reply.started":"2022-08-19T09:49:29.862904Z","shell.execute_reply":"2022-08-19T09:49:30.158257Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"This sums up our investigation of the gene expression data obtained by the citeseq test. We saw that the data is comparable to that obtained by multiome. One key difference is that both have unique genes that are only measured in one test. Later, we will also investigate if the comparability of the data is due to preceding normalization of the raw data obtained by the tests and also address the question of what an expression value of e. g. 2.4 actually means.","metadata":{}},{"cell_type":"markdown","source":"
\n

2.3 Surface protein level data

\n
\n\nLastly, we will have a look at the surface protein levels data gathered by citeseq.","metadata":{}},{"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-19T09:49:30.161148Z","iopub.execute_input":"2022-08-19T09:49:30.161985Z","iopub.status.idle":"2022-08-19T09:49:30.831174Z","shell.execute_reply.started":"2022-08-19T09:49:30.161948Z","shell.execute_reply":"2022-08-19T09:49:30.829860Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Compared to what we have seen so far, the number of columns in this data set is quite small. We have measurements of 140 features per cell. Most of the names start with CD, which is short for \"Cluster of differentiation\". CDs are used to classify surface molecules a cell expresses. This information can then be used to get an idea of what kind of cell is present, or what function this cell is supposed to serve in the body (I am not sure if my understanding here is even remotely accurate).\n\nLet's forget about the biological view for a moment and focus on the data science centric view. As we will see in the next cell, we have no zero values in this dataset, and thus much of the computation we did before (counting non-zero values, for example) is pointless. We will look at other features:","metadata":{}},{"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-19T09:49:30.832681Z","iopub.execute_input":"2022-08-19T09:49:30.833007Z","iopub.status.idle":"2022-08-19T09:49:31.214535Z","shell.execute_reply.started":"2022-08-19T09:49:30.832978Z","shell.execute_reply":"2022-08-19T09:49:31.213381Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We would also like to know if the absence of zero values could be due to inaccuracy in measurements. The following code checks for that. ATTENTION: As of right now, the threshold is completely arbitrary and will be revised when I know how accurate the test is!","metadata":{}},{"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-19T09:49:31.216064Z","iopub.execute_input":"2022-08-19T09:49:31.216427Z","iopub.status.idle":"2022-08-19T09:49:41.729428Z","shell.execute_reply.started":"2022-08-19T09:49:31.216395Z","shell.execute_reply":"2022-08-19T09:49:41.728338Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"It does not look like inaccuracy in measurements is the reason for the absence of 0 values, but as said, this needs to be checked!\n\n# 3 Other Data\nOn top of the data we already saw, there is also a file containing metadata and two files needed to make a submission for the competition. We will first have a look at the metadata and then at the competition-related data.\n\n
\n

3.1 Inspection of metadata.csv

\n
\n\nThe Metadata is stored in the file metadata.csv. The file contains one row for each cell in the dataset and provides some additional information about that cell.","metadata":{}},{"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-19T10:06:59.028843Z","iopub.execute_input":"2022-08-19T10:06:59.029324Z","iopub.status.idle":"2022-08-19T10:06:59.364839Z","shell.execute_reply.started":"2022-08-19T10:06:59.029287Z","shell.execute_reply":"2022-08-19T10:06:59.363365Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"For each cell, we see the **day** column tells us on which day the test was performed. The test that was actually performed is shown in the **technology** column. Note that experiments started on day 1, therefore the first tests were performed one day after the cells were injected with Neupogen for the first time ([compare here](https://allcells.com/research-grade-tissue-products/mobilized-leukopak/)). For each of the four donors, we have a **donor** ID giving us information about the origin of the cell. Lastly, there is a **cell_type** column. The cell types are labels assigned by humans. They might be imprecise and this information is not available for test data since it would be possible to draw conclusions about surface protein levels, for example. It is not clear if we can make use of that later (maybe we can use it for creating balanced splits, but we will see).\n\n[jirkaborovec](https://www.kaggle.com/jirkaborovec) already composed a great analysis of the metadata in this [notebook](https://www.kaggle.com/code/jirkaborovec/mmscel-inst-eda-stat-predictions). I will copy some parts here for readability and add some comments myself.","metadata":{}},{"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-19T10:40:28.209414Z","iopub.execute_input":"2022-08-19T10:40:28.209801Z","iopub.status.idle":"2022-08-19T10:40:28.564472Z","shell.execute_reply.started":"2022-08-19T10:40:28.209769Z","shell.execute_reply":"2022-08-19T10:40:28.562379Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"As we can see, the cell data is pretty balanced. We have almost an equal number of cells from each donor (the big numbers in the first picture are the donor ids). Also, the days of the experiment were fairly balanced. The last day, day 10, only receives an 11% share of the cells. Day 10 is also the only day not present in the train data at all! Also, the train set does not contain any data from donor 27678!\n\nWe have slightly more data available for the multinome test. It is not the worst since our model also has to predict many more features for that test.\n\nThe distribution of data in general is pretty well balanced (e. g., the number of tests taken / test / day is well distributed). For a more in-depth analysis of the metadata, I highly recommend the already mentioned [notebook](https://www.kaggle.com/code/jirkaborovec/mmscel-inst-eda-stat-predictions).\n\n
\n

3.2 Inspection of submission helpers

\n
\n\nThere are two files related to result submission. evaluation_ids.csv specifies the data that needs to be submitted to the competition, and sample_submission.csv is meant as a guide for formatting the submitted data.","metadata":{}},{"cell_type":"markdown","source":"# X Notes\n\nThere are still some open questions in the text we need to address. Also, we want to get an understanding of the accuracy of the values and thus how many bits we will take for storage of data in the final data format.","metadata":{}}]}