{"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":"# Setup","metadata":{}},{"cell_type":"code","source":"!pip install kaggle\n!pip install hdf5plugin~=2.0","metadata":{"jupyter":{"source_hidden":true},"execution":{"iopub.status.busy":"2022-09-13T05:44:31.472147Z","iopub.execute_input":"2022-09-13T05:44:31.473107Z","iopub.status.idle":"2022-09-13T05:44:55.857251Z","shell.execute_reply.started":"2022-09-13T05:44:31.473039Z","shell.execute_reply":"2022-09-13T05:44:55.855410Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Import all the core libraries\nimport numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\nimport seaborn as sns\nimport os\nimport gc\n\nimport h5py\nimport hdf5plugin\n\n# Change pandas display setting\npd.set_option('display.max_columns', None)\npd.set_option('display.max_rows', None)\n\n# Import sklearn and all the machine-learning related libraries\nfrom sklearn.metrics import mean_absolute_error\nfrom sklearn.metrics import confusion_matrix\nfrom sklearn.model_selection import train_test_split\nfrom sklearn.metrics import mean_absolute_error\nimport lightgbm as lgbm\nfrom lightgbm import LGBMClassifier","metadata":{"execution":{"iopub.status.busy":"2022-09-13T05:44:55.860819Z","iopub.execute_input":"2022-09-13T05:44:55.861455Z","iopub.status.idle":"2022-09-13T05:44:58.448917Z","shell.execute_reply.started":"2022-09-13T05:44:55.861391Z","shell.execute_reply":"2022-09-13T05:44:58.447438Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# calculate file size in KB, MB, GB\ndef convert_bytes(size):\n    \"\"\" Convert bytes to KB, or MB or GB\"\"\"\n    for x in ['bytes', 'KB', 'MB', 'GB', 'TB']:\n        if size < 1024.0:\n            return \"%3.1f %s\" % (size, x)\n        size /= 1024.0","metadata":{"jupyter":{"source_hidden":true},"execution":{"iopub.status.busy":"2022-09-13T05:44:58.450609Z","iopub.execute_input":"2022-09-13T05:44:58.451021Z","iopub.status.idle":"2022-09-13T05:44:58.458760Z","shell.execute_reply.started":"2022-09-13T05:44:58.450983Z","shell.execute_reply":"2022-09-13T05:44:58.457451Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Competition Data Check\nFirst of all, check the size of the available files.\nAs only 16GB RAM and 73.1GB Disk are assigned to a Kaggle notebook, I need to always think about the utilization of those resources.\nThus, firstly check the size of each file in the competition.","metadata":{}},{"cell_type":"code","source":"file_list = [\n    'train_cite_inputs.h5', \n    'train_cite_targets.h5', \n    'train_multi_inputs.h5', \n    'train_multi_targets.h5',\n    'test_cite_inputs.h5', \n    'test_multi_inputs.h5'\n    ]\nfor f in file_list:\n  f_path = f'../input/open-problems-multimodal/{f}' \n  f_size = os.path.getsize(f_path)\n  f_size_converted = convert_bytes(f_size)\n  print(f'{f} :', f_size_converted)","metadata":{"execution":{"iopub.status.busy":"2022-09-13T05:44:58.461510Z","iopub.execute_input":"2022-09-13T05:44:58.462418Z","iopub.status.idle":"2022-09-13T05:44:58.477723Z","shell.execute_reply.started":"2022-09-13T05:44:58.462345Z","shell.execute_reply":"2022-09-13T05:44:58.476363Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Given the assigned resources, I should not try to load those files into Pandas Dataframe. Particularly, train_cite_targets.h5 and \ntrain_multi_inputs.h5 are too massive to expand in this notebook.","metadata":{}},{"cell_type":"markdown","source":"# Goal of This Notebook\nIn this, I attempt to read one of the huge training file, train_multi_inputs.h5 and visualize the multiome data through a heatmap from the viewpoint of a non-expert not familiar with multiome data, so that more people like in me can participate in and contribute to this competition.","metadata":{}},{"cell_type":"markdown","source":"# Understanding Multiome Data from The Description\nData description mentions the Multiome data as follows:  \n> For the Multiome samples: given chromatin accessibility, predict gene expression.\n\n>train/test_multi_inputs.h5 - ATAC-seq peak counts transformed with TF-IDF using the default log(TF) * log(IDF) output (chromatin accessibility), with rows corresponding to cells and columns corresponding to the location of the genome whose level of accessibility is measured, here identified by the genomic coordinates on reference genome GRCh38 provided in the 10x References - 2020-A (July 7, 2020).\n\n> train_multi_targets.h5 - RNA gene expression levels as library-size normalized and log1p transformed counts for the same cells.\n\nWhat I understand after a few reads of those sentences is the followings:\nMultiome training input data (train_multi_inputs.h5) has rows of cells. Every column represents a location of the genome and its value shows the chromatin accessibility. <br>\nMultiome training target data (train_multi_targets.h5) has rows of cells. Every column represents a location of the genome and its value shows the gene expression level. <br>","metadata":{}},{"cell_type":"markdown","source":"# Peeking at The Competition Data\nI am skeptical about what I learned through the reading.   \nSo, in the following cells, I try to look at a fragment of data.\n\nFirst, load the file using h5py.","metadata":{}},{"cell_type":"code","source":"train_multi_inputs = 'train_multi_inputs.h5'\nf_path = f'../input/open-problems-multimodal/{train_multi_inputs}'\nf = h5py.File(f_path, 'r')\ngroup = f['train_multi_inputs']\nprint(list(group.keys()))\n\n# Check out dataset inside the group.\nfor group_key in group.keys():\n    print('Scanning ' + group_key + '...')\n    print('Shape: '+str(group[group_key].shape))  \n# Inside a hdf5 file, the data has a tree structure; each node is classified as either Dataset or Group.\n# Thus, if the structure is not known, the data should be processed as per the data type using isinstance(group[group_key], h5py.Dataset).","metadata":{"execution":{"iopub.status.busy":"2022-09-13T05:44:58.479175Z","iopub.execute_input":"2022-09-13T05:44:58.479602Z","iopub.status.idle":"2022-09-13T05:44:58.735825Z","shell.execute_reply.started":"2022-09-13T05:44:58.479564Z","shell.execute_reply":"2022-09-13T05:44:58.734075Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"A bit of tweaking and checking of the h5 file gives the following ideas.\n* train_multi_inputs.h5 consists 4 Dataset, namely axis0, axis1, block0_items, and block0_values. \n* block0_values contains the core data. The row length is 105,942 whilst the column length is 228,942. So many columns!\n* axis1 has the same length as the rows. It is likely that this represents cells.\n* axis0/block0_items has the same length as the columns. Either one shall have the identifiers of genome locations.\n\n","metadata":{}},{"cell_type":"markdown","source":"Now, let's print some of the data for each Dataset to validate those ideas. <br>\nLet's start from axis1.","metadata":{}},{"cell_type":"code","source":"# Print 10 records of axis1 as numpy array\nnp.char.decode(group['axis1'][:10])","metadata":{"execution":{"iopub.status.busy":"2022-09-13T05:44:58.738424Z","iopub.execute_input":"2022-09-13T05:44:58.739707Z","iopub.status.idle":"2022-09-13T05:44:58.754560Z","shell.execute_reply.started":"2022-09-13T05:44:58.739653Z","shell.execute_reply":"2022-09-13T05:44:58.752842Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"As speculated, each value contains 12 bytes of characters and numbers. \nIt is pretty much the same as cell_id in the Data section of this competition.<br>\nSo, I keep the data as cell_ids.","metadata":{}},{"cell_type":"code","source":"cell_ids = np.char.decode(group['axis1'])","metadata":{"execution":{"iopub.status.busy":"2022-09-13T05:44:58.756467Z","iopub.execute_input":"2022-09-13T05:44:58.757161Z","iopub.status.idle":"2022-09-13T05:44:58.838863Z","shell.execute_reply.started":"2022-09-13T05:44:58.757121Z","shell.execute_reply":"2022-09-13T05:44:58.837403Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"I have a guess that axis0 and block0_items are essentially the same. \nLet's check their data and confirm it.","metadata":{}},{"cell_type":"code","source":"# Print 10 records of axis0 as numpy array\nnp.char.decode(group['axis0'][:10])","metadata":{"execution":{"iopub.status.busy":"2022-09-13T05:44:58.840809Z","iopub.execute_input":"2022-09-13T05:44:58.841828Z","iopub.status.idle":"2022-09-13T05:44:58.855167Z","shell.execute_reply.started":"2022-09-13T05:44:58.841769Z","shell.execute_reply":"2022-09-13T05:44:58.854099Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Print 10 records of block0_items as numpy array\nnp.char.decode(group['block0_items'][:10])","metadata":{"execution":{"iopub.status.busy":"2022-09-13T05:44:58.856766Z","iopub.execute_input":"2022-09-13T05:44:58.858103Z","iopub.status.idle":"2022-09-13T05:44:58.873163Z","shell.execute_reply.started":"2022-09-13T05:44:58.858060Z","shell.execute_reply":"2022-09-13T05:44:58.871876Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Bingo. They shall be the same; I retain only axis0 as genome_locs.","metadata":{}},{"cell_type":"code","source":"genome_locs = np.char.decode(group['axis0'])","metadata":{"execution":{"iopub.status.busy":"2022-09-13T05:44:58.876586Z","iopub.execute_input":"2022-09-13T05:44:58.876980Z","iopub.status.idle":"2022-09-13T05:44:59.075576Z","shell.execute_reply.started":"2022-09-13T05:44:58.876947Z","shell.execute_reply":"2022-09-13T05:44:59.074195Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Just in case, I google about the value to make sure if the value really suggests the genome location.\nI could validate it in the link:\nhttps://www.ncbi.nlm.nih.gov/nuccore/GL000194.1","metadata":{}},{"cell_type":"markdown","source":"By now, I became certain that block0_values stores the actual dataset that we are going to train, and each value represents the chromatine accessibility. ","metadata":{}},{"cell_type":"code","source":"train_multi_inputs_data = group['block0_values']","metadata":{"execution":{"iopub.status.busy":"2022-09-13T05:44:59.077238Z","iopub.execute_input":"2022-09-13T05:44:59.077663Z","iopub.status.idle":"2022-09-13T05:44:59.084717Z","shell.execute_reply.started":"2022-09-13T05:44:59.077630Z","shell.execute_reply":"2022-09-13T05:44:59.083096Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Our ultimate mission is to predict the gene expression levels. <br>\nThat said, my notebook is only to visualize this chromatine accessibility data in the available genome locations.","metadata":{}},{"cell_type":"markdown","source":"# Creating Heatmap of Chromatine Accessability\nLoading the whole data and visualizing it will not be possible with the limited CPU avilability. But still, I do not want to cut any records. <br>\nSo this time, I slice the dataset horizontally to get the first 100 columns and save it into a file.<br>\nDesclaimer: I decided to keep it into a file in my local as I did not want to keep the data in the memory only as the CPU was already reaching 60%. ","metadata":{}},{"cell_type":"code","source":"%%time\nchunksize = 100;\nchunk_data = train_multi_inputs_data[:, 0:chunksize]\nnp.savetxt(f'train_multi_input_100.csv', chunk_data, delimiter=',', fmt='%1.8f')","metadata":{"execution":{"iopub.status.busy":"2022-09-13T05:44:59.086693Z","iopub.execute_input":"2022-09-13T05:44:59.087237Z","iopub.status.idle":"2022-09-13T05:48:11.176236Z","shell.execute_reply.started":"2022-09-13T05:44:59.087199Z","shell.execute_reply":"2022-09-13T05:48:11.174817Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":" ","metadata":{}},{"cell_type":"code","source":"gc.collect()\n# Read the data stored in the temporary file \nchunk_data = np.loadtxt(\"train_multi_input_100.csv\", delimiter=',')","metadata":{"execution":{"iopub.status.busy":"2022-09-13T05:48:11.178002Z","iopub.execute_input":"2022-09-13T05:48:11.178547Z","iopub.status.idle":"2022-09-13T05:48:19.267256Z","shell.execute_reply.started":"2022-09-13T05:48:11.178488Z","shell.execute_reply":"2022-09-13T05:48:19.265814Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Visualize it using plt.imshow\nplt.imshow(chunk_data, cmap='Blues', aspect='auto', vmin=0, vmax=1)","metadata":{"execution":{"iopub.status.busy":"2022-09-13T05:48:30.265490Z","iopub.execute_input":"2022-09-13T05:48:30.266386Z","iopub.status.idle":"2022-09-13T05:48:32.848220Z","shell.execute_reply.started":"2022-09-13T05:48:30.266318Z","shell.execute_reply":"2022-09-13T05:48:32.846844Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The heatmap drew some light-and-dark bands, indicating chromatines are more accessible to certain regions of the genome.","metadata":{}},{"cell_type":"markdown","source":"# Convert to Sparse Matrix\nAt this point, we see some patterns in the training data.\nHowever, the heatmap visualizes only the first 100 columns.\nWhen I increase the column number, the CPU usage goes up higher. \nGood news is that it is apparent that there are so many 0s in the data.\nSo we can convert the data to sparse matrix so that we can reduce the data size to be loaded.","metadata":{}},{"cell_type":"markdown","source":"Check the original data size of the first 100 columns.","metadata":{}},{"cell_type":"code","source":"convert_bytes(chunk_data.nbytes)","metadata":{"execution":{"iopub.status.busy":"2022-09-13T05:48:36.944536Z","iopub.execute_input":"2022-09-13T05:48:36.944990Z","iopub.status.idle":"2022-09-13T05:48:36.953042Z","shell.execute_reply.started":"2022-09-13T05:48:36.944954Z","shell.execute_reply":"2022-09-13T05:48:36.951456Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Convert the data to sparse matrix and check the data size after the conversion.","metadata":{}},{"cell_type":"code","source":"from scipy.sparse import csr_matrix\nchunk_data_s_matrix = csr_matrix(chunk_data)\nconvert_bytes(chunk_data_s_matrix.data.nbytes)","metadata":{"execution":{"iopub.status.busy":"2022-09-13T05:48:39.292608Z","iopub.execute_input":"2022-09-13T05:48:39.293483Z","iopub.status.idle":"2022-09-13T05:48:39.395473Z","shell.execute_reply.started":"2022-09-13T05:48:39.293432Z","shell.execute_reply":"2022-09-13T05:48:39.394265Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The sparse matrix conversion could downsize the data.\nUsing plt.spy, I am drawing the same heatmap below.","metadata":{}},{"cell_type":"code","source":"# visualize the sparse matrix with Spy\nplt.spy(chunk_data_s_matrix, aspect='auto', markersize=0.1)","metadata":{"execution":{"iopub.status.busy":"2022-09-13T05:48:56.778393Z","iopub.execute_input":"2022-09-13T05:48:56.779388Z","iopub.status.idle":"2022-09-13T05:48:57.284659Z","shell.execute_reply.started":"2022-09-13T05:48:56.779310Z","shell.execute_reply":"2022-09-13T05:48:57.283237Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Lastly, as I managed to successfully find a way to downsize the data and project it in the heatmap without crashing the CPU, I am going to increase the column numbers to be loaded from 100 to 5000.","metadata":{}},{"cell_type":"code","source":"# train_multi_inputs_data_s_matrix = csr_matrix(train_multi_inputs_data)\n# convert_bytes(train_multi_inputs_data_s_matrix.data.nbytes)","metadata":{"execution":{"iopub.status.busy":"2022-09-13T05:48:19.371078Z","iopub.status.idle":"2022-09-13T05:48:19.371735Z","shell.execute_reply.started":"2022-09-13T05:48:19.371415Z","shell.execute_reply":"2022-09-13T05:48:19.371446Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\ngc.collect()\nchunksize = 5000;\nchunk_data = train_multi_inputs_data[:, 0:chunksize]\nprint(\"original data size:\", convert_bytes(chunk_data.nbytes))\n# Convert to sparse matrix\nchunk_data_s_matrix = csr_matrix(chunk_data)\nprint(\"data size after conversion:\", convert_bytes(chunk_data_s_matrix.data.nbytes))","metadata":{"execution":{"iopub.status.busy":"2022-09-13T05:51:30.675181Z","iopub.execute_input":"2022-09-13T05:51:30.676195Z","iopub.status.idle":"2022-09-13T05:53:35.933811Z","shell.execute_reply.started":"2022-09-13T05:51:30.676136Z","shell.execute_reply":"2022-09-13T05:53:35.931864Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# visualize the sparse matrix with Spy\nplt.figure(figsize=(20,12))\nplt.spy(chunk_data_s_matrix, aspect='auto', markersize=0.01)","metadata":{"execution":{"iopub.status.busy":"2022-09-13T05:55:39.582045Z","iopub.execute_input":"2022-09-13T05:55:39.582608Z","iopub.status.idle":"2022-09-13T05:55:49.740750Z","shell.execute_reply.started":"2022-09-13T05:55:39.582548Z","shell.execute_reply":"2022-09-13T05:55:49.737622Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"This makes clear that there are bands observed vertically too.\nMeaning, it is likely that the chromatine accessibility also depends on the cell type.","metadata":{}},{"cell_type":"markdown","source":"Lastly, I am not an expert of biology or any -ome studies. What I know about central dogma is pretty much what you learned through biology classes in high school and what is written in the available [learning resources](https://openproblems.bio/neurips_docs/data/about_multimodal/). <br> \nBut with this heatmap, I can imagine that the nucleosome, the basic repeating subunit of chromatin, has some areas where histons are very/little crowded with each other. If an area is less crowded, transcription shall be more active that the gene \nexpression level may be high.\n\nI hope that having this image may help you build the right model.","metadata":{}}]}