{"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":"# EDA for the Multimodal Single-Cell Integration Competition","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19"}},{"cell_type":"code","source":"import os, gc, scipy.sparse\nimport pandas as pd\nimport matplotlib.pyplot as plt\nimport seaborn as sns\nimport numpy as np\nfrom colorama import Fore, Back, Style\nfrom sklearn.decomposition import TruncatedSVD\nfrom matplotlib.ticker import PercentFormatter\n\nDATA_DIR = \"/kaggle/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,"execution":{"iopub.status.busy":"2022-10-30T15:33:11.556073Z","iopub.execute_input":"2022-10-30T15:33:11.557438Z","iopub.status.idle":"2022-10-30T15:33:11.569332Z","shell.execute_reply.started":"2022-10-30T15:33:11.557369Z","shell.execute_reply":"2022-10-30T15:33:11.568047Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"A little trick to save time with pip: On the first run of the notebook, we need to install the `tables` module with pip. If the module is already installed (after a restart of the notebook, for instance), pip wastes 10 seconds by checking whether a newer version exists. We can skip this check by testing for the presence of the module in a simple if statement.","metadata":{}},{"cell_type":"code","source":"%%time\n# If you see a warning \"Failed to establish a new connection\" running this cell,\n# go to \"Settings\" on the right hand side, \n# and turn on internet. Note, you need to be phone verified.\n# We need this library to read HDF files.\nif not os.path.exists('/opt/conda/lib/python3.7/site-packages/tables'):\n    !pip install --quiet tables\n","metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"execution":{"iopub.status.busy":"2022-10-30T15:32:28.717410Z","iopub.execute_input":"2022-10-30T15:32:28.718229Z","iopub.status.idle":"2022-10-30T15:32:45.588955Z","shell.execute_reply.started":"2022-10-30T15:32:28.718178Z","shell.execute_reply":"2022-10-30T15:32:45.587797Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# The metadata table\n\nThe metadata table (which describes training and test data) shows us:\n- There is data about 281528 unique cells.\n- The cells belong to five days, four donors, eight cell types (including one type named 'hidden'), and two technologies.\n- The metadata table has no missing values.\n\n**Insight:** \n- Every cell is used only on a single day and then discarded. There are no time series over single cells.\n- The two technologies do not share cells. It looks like we may create two completely independent models, one per technology, even if they share the same four donors. It's two Kaggle competitions in one!\n- As the models are independent, it is a good idea to work with two separate notebooks, one for CITEseq, the other one for Multiome.\n- Donor and cell_type are categorical features, which can be one-hot encoded. ","metadata":{}},{"cell_type":"code","source":"df_meta = pd.read_csv(FP_CELL_METADATA, index_col='cell_id')\ndisplay(df_meta)\nif not df_meta.index.duplicated().any(): print('All cell_ids are unique.')\nif not df_meta.isna().any().any(): print('There are no missing values.')\n    ","metadata":{"execution":{"iopub.status.busy":"2022-10-30T15:32:45.590555Z","iopub.execute_input":"2022-10-30T15:32:45.590936Z","iopub.status.idle":"2022-10-30T15:32:46.180807Z","shell.execute_reply.started":"2022-10-30T15:32:45.590899Z","shell.execute_reply":"2022-10-30T15:32:46.179569Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"_, axs = plt.subplots(2, 2, figsize=(11, 6))\nfor col, ax in zip(['day', 'donor', 'cell_type', 'technology'], axs.ravel()):\n    vc = df_meta[col].astype(str).value_counts()\n    if col == 'day':\n        vc.sort_index(key = lambda x : x.astype(int), ascending=False, inplace=True)\n    else:\n        vc.sort_index(ascending=False, inplace=True)\n    ax.barh(vc.index, vc, color=['MediumSeaGreen'])\n    ax.set_ylabel(col)\n    ax.set_xlabel('# cells')\nplt.tight_layout(h_pad=4, w_pad=4)\nplt.suptitle('Metadata distribution', y=1.04, fontsize=20)\nplt.show()\n","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-10-30T15:32:46.183334Z","iopub.execute_input":"2022-10-30T15:32:46.183734Z","iopub.status.idle":"2022-10-30T15:32:47.257894Z","shell.execute_reply.started":"2022-10-30T15:32:46.183701Z","shell.execute_reply":"2022-10-30T15:32:47.256715Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The CITEseq measurements took place on four days, the Multiome measurements on five (except that there are no measurements for donor 27678 on day 4. For every combination of day, donor and technology, there are around 8000 cells:","metadata":{}},{"cell_type":"code","source":"# From https://www.kaggle.com/code/peterholderrieth/getting-started-data-loading\ndf_meta_cite = df_meta[df_meta.technology==\"citeseq\"]\ndf_meta_multi = df_meta[df_meta.technology==\"multiome\"]\n\nfig, axs = plt.subplots(1,2,figsize=(12,6))\ndf_cite_cell_dist = df_meta_cite[[\"day\",\"donor\"]].value_counts().to_frame()\\\n                .sort_values(\"day\").reset_index()\\\n                .rename(columns={0:\"# cells\"})\nsns.barplot(data=df_cite_cell_dist, x=\"day\",hue=\"donor\",y=\"# cells\", ax=axs[0])\naxs[0].set_title(f\"{len(df_meta_cite)} cells measured with CITEseq\")\n\ndf_multi_cell_dist = df_meta_multi[[\"day\",\"donor\"]].value_counts().to_frame()\\\n                .sort_values(\"day\").reset_index()\\\n                .rename(columns={0:\"# cells\"})\nsns.barplot(data=df_multi_cell_dist, x=\"day\",hue=\"donor\",y=\"# cells\", ax=axs[1])\naxs[1].set_title(f\"{len(df_meta_multi)} cells measured with Multiome\")\nplt.suptitle('# Cells per day, donor and technology', y=1.04, fontsize=20)\nplt.show()\nprint('Average:', round(len(df_meta) / 35))","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-10-30T15:32:47.259621Z","iopub.execute_input":"2022-10-30T15:32:47.259970Z","iopub.status.idle":"2022-10-30T15:32:47.863667Z","shell.execute_reply.started":"2022-10-30T15:32:47.259938Z","shell.execute_reply":"2022-10-30T15:32:47.862414Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"A diagram (taken from the [competition homepage](https://www.kaggle.com/competitions/open-problems-multimodal/data)) illustrates the relationships:\n- In the CITEseq competition, we train on data for three donors and three days. For the public leaderboard, we predict the fourth donor; for the private leaderboard, we predict another day for all donors.\n- In the Multiome competition, we train on three donors and four days. For the public leaderboard, we predict the fourth donor; for the private leaderboard, we predict a fifth day.\n\n**Insight:** As the data is grouped and we have to predict unseen groups (the last donor or the last day), we might choose `GroupKFold` as cross-validation scheme.\n\n![Diagram](https://www.googleapis.com/download/storage/v1/b/kaggle-user-content/o/inbox%2F4308072%2F23e8c1f6faea1453998544cdc116a20e%2FNeurIPS%202022%20-%20Frame%204.jpg?generation=1660755395301873&alt=media)","metadata":{}},{"cell_type":"markdown","source":"# The time series of the cell types\n\nIn this competition, we are working with a time series. The cells in the given data can be classified into seven cell types. If we plot the ratio of cell types by day, we see that during the course of the experiment the hematopoietic stem cells (HSC, green in the diagram) slowly change into other cell types. We begin with 50 % HSC, and on day 7 only 20 % HSC remain. We don't get the cell types for day 10, but we may guess that there will be even fewer hematopoietic stem cells left.\n","metadata":{}},{"cell_type":"code","source":"daily_cell_types = df_meta_cite.groupby(['day', 'cell_type']).size().unstack()\ndaily_cell_types[daily_cell_types.columns] = daily_cell_types.values / daily_cell_types.values.sum(axis=1).reshape(-1, 1)\n_, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 4), sharey=True)\nfor cell_type in daily_cell_types.columns:\n    ax1.plot(daily_cell_types.index,\n             daily_cell_types[cell_type],\n             label=cell_type,\n             lw = 7 if cell_type == 'HSC' else 3)\nax1.yaxis.set_major_formatter(PercentFormatter(xmax=1, decimals=0))\nax1.set_xticks(daily_cell_types.index)\nax1.legend()\nax1.set_xlabel('Day')\nax1.set_ylabel('Ratio')\nax1.set_title('CITEseq ratio of cell types by day')\n\ncolors = plt.rcParams['axes.prop_cycle'].by_key()['color']\ndaily_cell_types = df_meta_multi.groupby(['day', 'cell_type']).size().unstack()\ndaily_cell_types.drop(columns=['hidden'], inplace=True)\ndaily_cell_types.drop(index=[10], inplace=True)\ndaily_cell_types[daily_cell_types.columns] = daily_cell_types.values / daily_cell_types.values.sum(axis=1).reshape(-1, 1)\ndaily_cell_types.loc[10] = [0.002, 0.128, 0.03, 0.33, 0.19, 0.05, 0.27] # estimate by visual extrapolation\nfor i, cell_type in enumerate(daily_cell_types.columns):\n    ax2.plot(daily_cell_types.index[:-1],\n             daily_cell_types[cell_type].iloc[:-1],\n             label=cell_type,\n             color=colors[i],\n             lw = 7 if cell_type == 'HSC' else 3)\n    ax2.plot(daily_cell_types.index[-2:],\n             daily_cell_types[cell_type].iloc[-2:],\n             color=colors[i],\n             linestyle='dotted',\n             lw = 7 if cell_type == 'HSC' else 3)\nax2.yaxis.set_major_formatter(PercentFormatter(xmax=1, decimals=0))\nax2.set_xticks([2,3,4,7,10])\n#ax2.legend()\nax2.set_xlabel('Day')\nax2.set_ylabel('Ratio')\nax2.set_title('Multiome ratio of cell types by day')\nplt.show()\n\ndel daily_cell_types\n","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-10-30T15:51:15.717104Z","iopub.execute_input":"2022-10-30T15:51:15.717536Z","iopub.status.idle":"2022-10-30T15:51:16.139235Z","shell.execute_reply.started":"2022-10-30T15:51:15.717502Z","shell.execute_reply":"2022-10-30T15:51:16.138234Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# CITEseq inputs\n\nWe start by looking at CITEseq, which has the smaller datafiles than Multiome and is more tractable.\n\nThe CITEseq input files contain 70988 samples (i.e., cells) for train and 48663 samples for test. 70988 + 48663 = 119651, which matches the number of rows in the CITEseq metadata table. No values are missing.\n\nThe input data corresponds to RNA expression levels for 22050 genes (there are 22050 columns).\n\nThe data have dtype float32, which means we need 119651 * 22050 * 4 = 10553218200 bytes = 10.6 GByte of RAM just for the features (train and test) without the targets. \n\nOriginally, these RNA expression levels were counts (i.e., nonnegative integers), but they have been normalized and log1p-transformed. With the log1p transformation, the data remain nonnegative.\n\nMost columns have a minimum of zero, which means that for most genes there are cells which didn't express this gene (the count was 0). In fact, 78 % of all table entries are zero, and for some columns, the count is always zero.\n\n**Insight:**\n- This is big data. Make sure you don't waste RAM and use efficient algorithms.\n- Perhaps we should first load only the training data into RAM, fit one or more models, and then delete the training data from RAM before loading the test data.\n- The columns which are zero for every cell should be dropped before modeling.","metadata":{}},{"cell_type":"code","source":"%%time\n# Analyze train features\ndf_cite_train_x = pd.read_hdf(FP_CITE_TRAIN_INPUTS)\ndisplay(df_cite_train_x.head())\nprint('Shape:', df_cite_train_x.shape)\nprint(\"Missing values:\", df_cite_train_x.isna().sum().sum())\nprint(\"Genes which never occur in train:\", (df_cite_train_x == 0).all(axis=0).sum())\nprint(f\"Zero entries in train: {(df_cite_train_x == 0).sum().sum() / df_cite_train_x.size:.0%}\")\ncite_gene_names = list(df_cite_train_x.columns)","metadata":{"execution":{"iopub.status.busy":"2022-09-14T07:36:01.166165Z","iopub.execute_input":"2022-09-14T07:36:01.166640Z","iopub.status.idle":"2022-09-14T07:37:18.427962Z","shell.execute_reply.started":"2022-09-14T07:36:01.166606Z","shell.execute_reply":"2022-09-14T07:37:18.426654Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The distribution of the zeros can be visualized with the pyplot `spy()` function. This function plots a black dot for every nonzero entry of the array. The resulting image shows us that the differences between the columns are substantial: some columns are almost white (i.e., they contain mostly zeros), others are dark (i.e., they contain a lot of nonzero values).\n\nThe rows of the matrix look homogeneous.\n\n**Insight:** Maybe we can exploit the column differences for feature selection.","metadata":{}},{"cell_type":"code","source":"plt.figure(figsize=(10, 4))\nplt.spy(df_cite_train_x[:5000])\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2022-09-14T07:37:18.429409Z","iopub.execute_input":"2022-09-14T07:37:18.430351Z","iopub.status.idle":"2022-09-14T07:37:23.255834Z","shell.execute_reply.started":"2022-09-14T07:37:18.430313Z","shell.execute_reply":"2022-09-14T07:37:23.254956Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The histogram shows some artefacts because the data originally were integers. We don't show the zeros in the histogram, because with 78 % zeros, the histogram would have such a high peak at zero that we couldn't see anything else.\n\n**Insight:** The feature values are either 0 or between 2.9 and 12 - the distribution is far away from normal. If we apply statistical tests, this fact has to be taken into consideration.","metadata":{}},{"cell_type":"code","source":"%%time\nnonzeros = df_cite_train_x.values.ravel()\nnonzeros = nonzeros[nonzeros != 0] # comment this line if you want to see the peak at zero\nplt.figure(figsize=(16, 4))\nplt.gca().set_facecolor('#0057b8')\nplt.hist(nonzeros, bins=500, density=True, color='#ffd700')\nprint('Minimum nonzero value:', nonzeros.min())\ndel nonzeros\nplt.title(\"Histogram of nonzero RNA expression levels in train\")\nplt.xlabel(\"log1p-transformed expression count\")\nplt.ylabel(\"density\")\nplt.show()","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-09-14T07:37:23.257165Z","iopub.execute_input":"2022-09-14T07:37:23.258094Z","iopub.status.idle":"2022-09-14T07:38:03.910174Z","shell.execute_reply.started":"2022-09-14T07:37:23.258057Z","shell.execute_reply":"2022-09-14T07:38:03.908998Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"_, axs = plt.subplots(5, 4, figsize=(16, 16))\nfor col, ax in zip(df_cite_train_x.columns[:20], axs.ravel()):\n    nonzeros = df_cite_train_x[col].values\n    nonzeros = nonzeros[nonzeros != 0] # comment this line if you want to see the peak at zero\n    ax.hist(nonzeros, bins=100, density=True)\n    ax.set_title(col)\nplt.tight_layout(h_pad=2)\nplt.suptitle('Histograms of nonzero RNA expression levels for selected features', fontsize=20, y=1.04)\nplt.show()\ndel nonzeros","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-09-14T07:38:03.912057Z","iopub.execute_input":"2022-09-14T07:38:03.912526Z","iopub.status.idle":"2022-09-14T07:38:10.580860Z","shell.execute_reply.started":"2022-09-14T07:38:03.912479Z","shell.execute_reply":"2022-09-14T07:38:10.579508Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"As we have seen that most entries of df_cite_train_x are zero and memory is scarce, we convert the numpy array to a [compressed sparse row](https://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.csr_matrix.html) (CSR) matrix:","metadata":{}},{"cell_type":"code","source":"cell_index = df_cite_train_x.index\nmeta = df_meta_cite.reindex(cell_index)\ngc.collect()\ndf_cite_train_x = scipy.sparse.csr_matrix(df_cite_train_x.values)\n","metadata":{"execution":{"iopub.status.busy":"2022-09-14T07:38:10.584595Z","iopub.execute_input":"2022-09-14T07:38:10.585060Z","iopub.status.idle":"2022-09-14T07:39:07.949894Z","shell.execute_reply.started":"2022-09-14T07:38:10.585017Z","shell.execute_reply":"2022-09-14T07:39:07.948444Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We have freed enough memory to analyze the test data; afterwards we convert the test data to a CSR matrix as well.","metadata":{}},{"cell_type":"code","source":"df_cite_test_x = pd.read_hdf(FP_CITE_TEST_INPUTS)\nprint('Shape of CITEseq test:', df_cite_test_x.shape)\nprint(\"Missing values:\", df_cite_test_x.isna().sum().sum())\nprint(\"Genes which never occur in test: \", (df_cite_test_x == 0).all(axis=0).sum())\nprint(f\"Zero entries in test:  {(df_cite_test_x == 0).sum().sum() / df_cite_test_x.size:.0%}\")\n\n\ngc.collect()\ncell_index_test = df_cite_test_x.index\nmeta_test = df_meta_cite.reindex(cell_index_test)\ndf_cite_test_x = scipy.sparse.csr_matrix(df_cite_test_x.values)\n","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-09-14T07:39:07.951607Z","iopub.execute_input":"2022-09-14T07:39:07.952040Z","iopub.status.idle":"2022-09-14T07:40:38.356097Z","shell.execute_reply.started":"2022-09-14T07:39:07.951999Z","shell.execute_reply":"2022-09-14T07:40:38.355134Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# The data leak\n\nIt has been pointed out in several discussion posts (the first one was [CITEseq data: same RNA expression matrices from different donors in day2?](https://www.kaggle.com/competitions/open-problems-multimodal/discussion/349867) by @gwentea) that the first 7476 rows of test (day 2, donor 27678) are identical to the first 7476 rows of train (day 2, donor 32606):\n","metadata":{}},{"cell_type":"code","source":"print('Data leak:', (df_cite_train_x[:7476] == df_cite_test_x[:7476]).toarray().all())","metadata":{"execution":{"iopub.status.busy":"2022-09-14T07:40:38.358147Z","iopub.execute_input":"2022-09-14T07:40:38.359428Z","iopub.status.idle":"2022-09-14T07:40:47.739939Z","shell.execute_reply.started":"2022-09-14T07:40:38.359385Z","shell.execute_reply":"2022-09-14T07:40:47.738642Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Insight:**\n- Some mistake happened when the data was prepared; it can't be that 7476 cells of one donor return exactly the same measurements as 7476 cells of another donor.\n- These rows belong to the public test set; the private test set is not affected.\n- Kaggle has decided to exclude the 7476 wrong cells from scoring, see [this discussion](https://www.kaggle.com/competitions/open-problems-multimodal/discussion/350933).\n- (Before that decision, one had to copy the 7476 rows from the training targets into the test predictions to get a good score.)","metadata":{}},{"cell_type":"markdown","source":"# The distributions of train and test in feature space\n\nFor the next diagrams, we need to project the data (train and test together) to two dimensions. Usually I'd do that with a PCA, but the scikit-learn PCA implementation needs too much memory. Fortunately, TruncatedSVD does a similar projection and uses much less memory:","metadata":{}},{"cell_type":"code","source":"# Concatenate train and test for the SVD\nboth = scipy.sparse.vstack([df_cite_train_x, df_cite_test_x])\nprint(f\"Shape of both before SVD: {both.shape}\")\n\n# Project to two dimensions\nsvd = TruncatedSVD(n_components=2, random_state=1)\nboth = svd.fit_transform(both)\nprint(f\"Shape of both after SVD:  {both.shape}\")\n\n# Separate train and test\nX = both[:df_cite_train_x.shape[0]]\nXt = both[df_cite_train_x.shape[0]:]\n","metadata":{"execution":{"iopub.status.busy":"2022-09-14T07:40:47.741684Z","iopub.execute_input":"2022-09-14T07:40:47.742757Z","iopub.status.idle":"2022-09-14T07:41:44.328693Z","shell.execute_reply.started":"2022-09-14T07:40:47.742704Z","shell.execute_reply":"2022-09-14T07:41:44.327224Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The scatterplots below show the extent of the data for every day and every donor (SVD projection to two dimensions). The nine diagrams with orange dots make up the training data. The three diagrams with orange-red dots below are the public test set and the four dark red diagrams at the right are the private test set.\n\nThe gray area marks the complete training data (union of the nine orange diagrams), and the black area behind marks the test data.\n\nWe see that the distributions differ. In particular the day 2 distribution (left column of diagrams) is much less wide than the others.\n\n**Insight:**\n- For the public test we predict a previously unseen donor, for the private test a previously unseen day. This situation suggests that we use a GroupKFold for cross-validation, although its unclear whether we should group on day, donor, or both. \n- Every small diagram consists of only 8000 samples and the data is noisy. It will be difficult not to overfit.\n- As the public test set is three times smaller than the training dataset, we must not rely on the public leaderboard to evaluate models. \n- Day 7 (private test) covers areas of the feature space which occur neither in train nor in public test. We'll need a model which can extrapolate to this area, and the public leaderboard will give us no clue about the quality of this extrapolation.\n- The diagrams confirm the data leak described above: The lower two diagrams of day 2 are identical.","metadata":{}},{"cell_type":"code","source":"# Scatterplot for every day and donor\n_, axs = plt.subplots(4, 4, sharex=True, sharey=True, figsize=(12, 11))\nfor donor, axrow in zip([13176, 31800, 32606, 27678], axs):\n    for day, ax in zip([2, 3, 4, 7], axrow):\n        ax.scatter(Xt[:,0], Xt[:,1], s=1, c='k')\n        ax.scatter(X[:,0], X[:,1], s=1, c='lightgray')\n        if day != 7 and donor != 27678: # train\n            temp = X[(meta.donor == donor) & (meta.day == day)]\n            ax.scatter(temp[:,0], temp[:,1], s=1, c='orange')\n        else: # test\n            temp = Xt[(meta_test.donor == donor) & (meta_test.day == day)]\n            ax.scatter(temp[:,0], temp[:,1], s=1, c='darkred' if day == 7 else 'orangered')\n        ax.set_title(f'Donor {donor} day {day}')\n        ax.set_aspect('equal')\nplt.suptitle('CITEseq features, projected to the first two SVD components', y=0.95, fontsize=20)\nplt.show()\n","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-09-14T07:41:44.330344Z","iopub.execute_input":"2022-09-14T07:41:44.330729Z","iopub.status.idle":"2022-09-14T07:41:47.378927Z","shell.execute_reply.started":"2022-09-14T07:41:44.330695Z","shell.execute_reply":"2022-09-14T07:41:47.377765Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_cite_train_x, df_cite_test_x, X, Xt = None, None, None, None # release the memory","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-09-14T07:41:47.380884Z","iopub.execute_input":"2022-09-14T07:41:47.381558Z","iopub.status.idle":"2022-09-14T07:41:47.400810Z","shell.execute_reply.started":"2022-09-14T07:41:47.381519Z","shell.execute_reply":"2022-09-14T07:41:47.399457Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# CITEseq targets (surface protein levels)\n\nThe CITEseq output (target) file is much smaller - it has 70988 rows like the training input file, but only 140 columns. The 140 columns correspond to 140 proteins.\n\nThe targets are dsb-normalized surface protein levels. We plot the histograms of a few selected columns and see that the distributions vary: some columns are normally distributed, some columns are multimodal, some have other shapes, and there seem to be outliers.\n\n**Insight:**\n- This is a multi-output regression problem with 140 outputs. We won't be able to apply some standard methods of single-output regression - e.g., we can't compute the correlation between every feature and the target. \n- As the targets are so diverse, a one-size-fits-all approach might not give the best results.","metadata":{}},{"cell_type":"code","source":"df_cite_train_y = pd.read_hdf(FP_CITE_TRAIN_TARGETS)\ndisplay(df_cite_train_y.head())\nprint('Output shape:', df_cite_train_y.shape)\n\n_, axs = plt.subplots(5, 4, figsize=(16, 16))\nfor col, ax in zip(['CD86', 'CD270', 'CD48', 'CD8', 'CD7', 'CD14', 'CD62L', 'CD54', 'CD42b', 'CD2', 'CD18', 'CD36', 'CD328', 'CD224', 'CD35', 'CD57', 'TCRVd2', 'HLA-E', 'CD82', 'CD101'], axs.ravel()):\n    ax.hist(df_cite_train_y[col], bins=100, density=True)\n    ax.set_title(col)\nplt.tight_layout(h_pad=2)\nplt.suptitle('Selected target histograms (surface protein levels)', fontsize=20, y=1.04)\nplt.show()\n\ncite_protein_names = list(df_cite_train_y.columns)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-09-14T07:41:47.402424Z","iopub.execute_input":"2022-09-14T07:41:47.402776Z","iopub.status.idle":"2022-09-14T07:41:54.604893Z","shell.execute_reply.started":"2022-09-14T07:41:47.402744Z","shell.execute_reply":"2022-09-14T07:41:54.603717Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"A projection of the CITEseq targets to two dimensions again shows that the groups have different distributions.","metadata":{}},{"cell_type":"code","source":"svd = TruncatedSVD(n_components=2, random_state=1)\nX = svd.fit_transform(df_cite_train_y)\n\n# Scatterplot for every day and donor\n_, axs = plt.subplots(4, 4, sharex=True, sharey=True, figsize=(12, 11))\nfor donor, axrow in zip([13176, 31800, 32606, 27678], axs):\n    for day, ax in zip([2, 3, 4, 7], axrow):\n        if day != 7 and donor != 27678: # train\n            ax.scatter(X[:,0], X[:,1], s=1, c='lightgray')\n            temp = X[(meta.donor == donor) & (meta.day == day)]\n            ax.scatter(temp[:,0], temp[:,1], s=1, c='orange')\n        else: # test\n            ax.text(50, -25, '?', fontsize=100, color='gray', ha='center')\n        ax.set_title(f'Donor {donor} day {day}')\n        ax.set_aspect('equal')\nplt.suptitle('CITEseq target, projected to the first two SVD components', y=0.95, fontsize=20)\nplt.show()","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-09-14T07:41:54.606526Z","iopub.execute_input":"2022-09-14T07:41:54.607150Z","iopub.status.idle":"2022-09-14T07:41:57.417688Z","shell.execute_reply.started":"2022-09-14T07:41:54.607115Z","shell.execute_reply":"2022-09-14T07:41:57.416586Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_cite_train_x, df_cite_train_y, X, svd = None, None, None, None # release the memory","metadata":{"execution":{"iopub.status.busy":"2022-09-14T07:41:57.419281Z","iopub.execute_input":"2022-09-14T07:41:57.420097Z","iopub.status.idle":"2022-09-14T07:41:57.427099Z","shell.execute_reply.started":"2022-09-14T07:41:57.420060Z","shell.execute_reply":"2022-09-14T07:41:57.425980Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Name matching\n\nThe CITEseq task has genes as input and proteins as output. Genes encode proteins, and it is more or less known which genes encode which proteins. This information is encoded in the column names: The input dataframe has the genes as column names, and the target dataframe has the proteins as column names. According to the naming convention, the gene names contain the protein name as suffix after a '_'. \n\nIf we match the input column names with the target column names, we find 151 genes which encode a target protein (see the table below). It doesn't matter that some proteins are encoded by more than one gene (e.g., rows 146 and 147 of the table). We may assume that these 151 features will have a high feature importance in our models. \n\n**Insight:** If we apply dimensionality reduction (PCA, SVD, whatever) to the 22050 features, we should make sure that we don't reduce away the 151 features which encode the proteins we want to predict.","metadata":{}},{"cell_type":"code","source":"matching_names = []\nfor protein in cite_protein_names:\n    matching_names += [(gene, protein) for gene in cite_gene_names if protein in gene]\npd.DataFrame(matching_names, columns=['Gene', 'Protein'])","metadata":{"execution":{"iopub.status.busy":"2022-09-14T07:41:57.429048Z","iopub.execute_input":"2022-09-14T07:41:57.429496Z","iopub.status.idle":"2022-09-14T07:41:57.666013Z","shell.execute_reply.started":"2022-09-14T07:41:57.429461Z","shell.execute_reply":"2022-09-14T07:41:57.664530Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Multiome input\n\nThe Multiome dataset is much larger than the CITEseq part and way too large to fit into 16 GByte RAM:\n- train inputs:  105942 * 228942 float32 values (97 GByte)\n- train targets: 105942 *  23418 float32 values (10 GByte)\n- test inputs:    55935 * 228942 float32 values (13 GByte)\n\nFor this EDA, we read all the data to check for missing, zero and negative values, we plot a histogram and we look at the Y chromosome, but we don't analyze much more.\n\nNo values are missing.\n\nThe data consists of ATAC-seq peak counts transformed with [TF-IDF](https://scikit-learn.org/stable/modules/generated/sklearn.feature_extraction.text.TfidfTransformer.html). They are all nonnegative. In the sample we are looking at, 98 % of the entries are zero.\n\n**Insight:**\n- The Multiome data is even bigger data than the CITEseq data. Most of us can't afford notebooks with 97 GByte RAM.\n- Dimensionality reduction and feature selection will be key to get successful models.\n- Maybe we can do something with sparse array data structures.\n- Maybe we can train on a subset of the rows (and sacrifice some precision).\n- Maybe we can convert the data to float16.\n- We should never have training and test data in memory at the same time.\n- We may even want to look for algorithms which don't need all the training data in RAM at the same time. Neural networks are a candidate.\n- The columns which are zero for every cell should be dropped before modeling.\n- The training set has n_features > n_samples. We'll use algorithms which can deal with more features than samples.\n","metadata":{}},{"cell_type":"code","source":"%%time\nbins = 100\ncell_summary = pd.DataFrame()\n\ndef analyze_multiome_x(filename):\n    global cell_summary\n    start = 0\n    chunksize = 5000\n    total_rows = 0\n    maximum_x = 0\n\n    while True: # read the next chunk of the file\n        X = pd.read_hdf(filename, start=start, stop=start+chunksize)\n        if X.isna().any().any(): print('There are missing values.')\n        if (X < 0).any().any(): print('There are negative values.')\n        total_rows += len(X)\n        print(total_rows, 'rows read')\n\n        donors = df_meta_multi.donor.reindex(X.index) # metadata: donor of cell\n        days = df_meta_multi.day.reindex(X.index) # metadata: day of cell\n        chrY_cols = [f for f in X.columns if 'chrY' in f]\n        maximum_x = max(maximum_x, X[chrY_cols].values.ravel().max())\n        for donor in [13176, 31800, 32606, 27678]:\n            hist, _ = np.histogram(X[chrY_cols][donors == donor].values.ravel(), bins=bins, range=(0, 15))\n            chrY_histo[donor] += hist\n\n        cell_summary = pd.concat([cell_summary,\n                                  pd.DataFrame({'donor': donors,\n                                                'day': days,\n                                                'total': X.sum(axis=1),\n                                                'total_nonzero': (X != 0).sum(axis=1)})])\n        if len(X) < chunksize: break\n        start += chunksize\n\n    display(X.head(3))\n    print(f\"Zero entries in {filename}: {(X == 0).sum().sum() / X.size:.0%}\")\n\nchrY_histo = dict()\nfor donor in [13176, 31800, 32606, 27678]:\n    chrY_histo[donor] = np.zeros((bins, ), int)\n\n# Look at the training data\nanalyze_multiome_x(FP_MULTIOME_TRAIN_INPUTS)","metadata":{"execution":{"iopub.status.busy":"2022-09-14T07:41:57.667758Z","iopub.execute_input":"2022-09-14T07:41:57.668537Z","iopub.status.idle":"2022-09-14T07:54:13.804929Z","shell.execute_reply.started":"2022-09-14T07:41:57.668485Z","shell.execute_reply":"2022-09-14T07:54:13.803273Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"In the histogram, we again hide the peak for the 98 % zero values:","metadata":{}},{"cell_type":"code","source":"df_multi_train_x = pd.read_hdf(FP_MULTIOME_TRAIN_INPUTS, start=0, stop=5000)\nnonzeros = df_multi_train_x.values.ravel()\nnonzeros = nonzeros[nonzeros != 0] # comment this line if you want to see the peak at zero\nplt.figure(figsize=(16, 4))\nplt.gca().set_facecolor('#0057b8')\nplt.hist(nonzeros, bins=500, density=True, color='#ffd700')\ndel nonzeros\nplt.title(\"Histogram of nonzero feature values (subset)\")\nplt.xlabel(\"TFIDF-transformed peak count\")\nplt.ylabel(\"density\")\nplt.show()\n\ndel df_multi_train_x # free the memory","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-09-14T07:54:13.806929Z","iopub.execute_input":"2022-09-14T07:54:13.807485Z","iopub.status.idle":"2022-09-14T07:55:05.194245Z","shell.execute_reply.started":"2022-09-14T07:54:13.807438Z","shell.execute_reply":"2022-09-14T07:55:05.193041Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We do the same checks for the test data:","metadata":{}},{"cell_type":"code","source":"%%time\n# Look at the test data\nanalyze_multiome_x(FP_MULTIOME_TEST_INPUTS)\n","metadata":{"execution":{"iopub.status.busy":"2022-09-14T07:55:05.197002Z","iopub.execute_input":"2022-09-14T07:55:05.197525Z","iopub.status.idle":"2022-09-14T08:02:08.899396Z","shell.execute_reply.started":"2022-09-14T07:55:05.197471Z","shell.execute_reply":"2022-09-14T08:02:08.898096Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Batch effects\n\nBiological experiments are notorious for batch effects: If you repeat the same experiment several times, you get systematic differences in the measurements. We can show these batch effects for the ATACseq measurements. Because of the huge data size, we look only at the row sums of the data (i.e. the sum of all measured values of a cell). In the diagram, we plot these sums for every row of the train and test datasets, and we color the dots by the day of the experiment (day 2 is red, day 3 is green and so on).\n\nThe diagram clearly shows that day 3 (green) gave the highest measurements, day 7 (yellow) the lowest. There are differences among donors as well, albeit smaller ones.\n\n**Insight:**\n- We should try to normalize the data in a preprocessing step so that the batch effects don't affect the predictions.\n- Because the batch effects depend more on the day than on the donor, we should use a `GroupKFold` on days for cross-validation to validate that the differences between the days don't break the model.","metadata":{}},{"cell_type":"code","source":"def plot_batch_effects():\n    _, axs = plt.subplots(1, 2, sharey=True, figsize=(12, 10))\n    color = cell_summary.day.map({2: 'r', 3: 'g', 4: 'b', 7: 'y', 10: 'gray'})\n    axs[0].scatter(cell_summary.total_nonzero, np.arange(len(cell_summary)), s=0.1, c=color)\n    axs[0].set_xlabel('cell total')\n    axs[1].scatter(cell_summary.total, np.arange(len(cell_summary)), s=0.1, c=color)\n    axs[1].set_xlabel('cell total nonzeros')\n    axs[0].set_ylabel('cell')\n    axs[0].invert_yaxis()\n    plt.suptitle('Row totals colored by day', y=0.94, fontsize=20)\n    plt.show()\n    \nplot_batch_effects()","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-09-14T08:02:08.900934Z","iopub.execute_input":"2022-09-14T08:02:08.901266Z","iopub.status.idle":"2022-09-14T08:02:13.710035Z","shell.execute_reply.started":"2022-09-14T08:02:08.901235Z","shell.execute_reply":"2022-09-14T08:02:13.708830Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# The Y chromosome\n\nLet's plot a histogram of the nonzero values of the Y chromosome for every donor. ","metadata":{}},{"cell_type":"code","source":"plt.rcParams['savefig.facecolor'] = \"1.0\"\n_, axs = plt.subplots(1, 4, sharex=True, sharey=True, figsize=(14, 4))\nfor donor, ax in zip([13176, 31800, 32606, 27678], axs):\n    ax.set_title(f\"Donor {donor} {'(test)' if donor == 27678 else ''}\", fontsize=16)\n    total = chrY_histo[donor].sum()\n    ax.fill_between(range(bins-1), chrY_histo[donor][1:] / total, color='limegreen')\n    ax.get_xaxis().set_visible(False)\n    ax.get_yaxis().set_visible(False)\nplt.suptitle(\"Histogram of nonzero Y chromosome accessibility\", y=0.95, fontsize=20)\nplt.tight_layout()\nplt.show()","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-09-14T08:02:13.711407Z","iopub.execute_input":"2022-09-14T08:02:13.712340Z","iopub.status.idle":"2022-09-14T08:02:14.044840Z","shell.execute_reply.started":"2022-09-14T08:02:13.712300Z","shell.execute_reply":"2022-09-14T08:02:14.043456Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The histograms of the Y chromosomes illustrate the diversity of the donors: Donor 13176 seems to have (almost) no Y chromosome (maybe the few nonzero values are measuring errors). It appears that the donors are one woman and three men.\n\n**Insight:**\n- If we assume that the true values of donor 13176 are all zero, these data show us the magnitude of the measurement error.\n- Maybe we can create a new (binary) feature for the presence of the Y chromosome.\n- The diagram reminds us that the donors are different and that our model should be robust to these differences.","metadata":{}},{"cell_type":"markdown","source":"# Multiome target\n\nThe Multiome targets (RNA count data) are in similar shape as the CITEseq inputs: They have 105942 rows and 23418 columns. All targets are nonnegative and no values are missing.\n","metadata":{}},{"cell_type":"code","source":"%%time\ncell_summary = pd.DataFrame()\nstart = 0\nchunksize = 10000\ntotal_rows = 0\nwhile True:\n    df_multi_train_y = pd.read_hdf(FP_MULTIOME_TRAIN_TARGETS, start=start, stop=start+chunksize)\n    if df_multi_train_y.isna().any().any(): print('There are missing values.')\n    if (df_multi_train_y < 0).any().any(): print('There are negative values.')\n    total_rows += len(df_multi_train_y)\n    print(total_rows, 'rows read')\n\n    donors = df_meta_multi.donor.reindex(df_multi_train_y.index) # metadata: donor of cell\n    days = df_meta_multi.day.reindex(df_multi_train_y.index) # metadata: day of cell\n    cell_summary = pd.concat([cell_summary,\n                              pd.DataFrame({'donor': donors,\n                                            'day': days,\n                                            'total': df_multi_train_y.sum(axis=1),\n                                            'total_nonzero': (df_multi_train_y != 0).sum(axis=1)})])\n    \n    if len(df_multi_train_y) < chunksize: break\n    start += chunksize\n    \ndisplay(df_multi_train_y.head())\n","metadata":{"execution":{"iopub.status.busy":"2022-09-14T08:02:14.046727Z","iopub.execute_input":"2022-09-14T08:02:14.047081Z","iopub.status.idle":"2022-09-14T08:03:51.594802Z","shell.execute_reply.started":"2022-09-14T08:02:14.047049Z","shell.execute_reply":"2022-09-14T08:03:51.592846Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nnonzeros = df_multi_train_y.values.ravel()\nnonzeros = nonzeros[nonzeros != 0]\nplt.figure(figsize=(16, 4))\nplt.gca().set_facecolor('#0057b8')\nplt.hist(nonzeros, bins=500, density=True, color='#ffd700')\ndel nonzeros\nplt.title(\"Histogram of nonzero target values (based on a subset of the rows)\")\nplt.xlabel(\"log1p-transformed expression count\")\nplt.ylabel(\"density\")\nplt.show()\n\ndf_multi_train_y = None # release the memory","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-09-14T08:03:51.596984Z","iopub.execute_input":"2022-09-14T08:03:51.599902Z","iopub.status.idle":"2022-09-14T08:03:55.002297Z","shell.execute_reply.started":"2022-09-14T08:03:51.599763Z","shell.execute_reply":"2022-09-14T08:03:55.000993Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Above, we saw the batch effects in the ATACseq data; below we see that the gene expression dataset has its own, independent, batch effects: For the gene expression data, day 7 (yellow) clearly has the lowest totals. Day 3 (green), which deviated the most from the other days for the ATACseq data, is not far from the mean for the gene expressions. But again, differences among days seem to exceed differences among donors.\n\nThe good news here is that the competition metric does not depend on shifts in the target values: Even if for day 7 (and maybe for day 10) all measurements are shifted by a constant, the predictions can still achieve a good correlation with the true values. With a mean-squared-error metric, this would be different.","metadata":{}},{"cell_type":"code","source":"plot_batch_effects()","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-09-14T08:03:55.003972Z","iopub.execute_input":"2022-09-14T08:03:55.004356Z","iopub.status.idle":"2022-09-14T08:03:58.696332Z","shell.execute_reply.started":"2022-09-14T08:03:55.004313Z","shell.execute_reply":"2022-09-14T08:03:58.694907Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Summary\n\n- This is a big data competition. The more RAM you have, the better. In any case, you'll have to treat RAM as a scarce resource and write memory-efficient code.\n- As there are two competitions in one, you can choose to participate in only one of them (and merge your predictions with somebody else's predictions for the other part). CITEseq is easier to begin with.\n- Make sure you decide for a good cross-validation scheme so that you avoid the shakedown at the end of the competition.\n\nAfter this EDA, you may want to look at some models which generate predictions in spite of having only 16 GB RAM:\n- [MSCI CITEseq Keras Quickstart](https://www.kaggle.com/code/ambrosm/msci-citeseq-keras-quickstart)\n- [MSCI CITEseq Quickstart](https://www.kaggle.com/ambrosm/msci-citeseq-quickstart) (using LightGBM)\n- [MSCI Multiome Quickstart](https://www.kaggle.com/ambrosm/msci-multiome-quickstart) (using ridge regression)","metadata":{}},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}