{"metadata":{"kernelspec":{"name":"ir","display_name":"R","language":"R"},"language_info":{"name":"R","codemirror_mode":"r","pygments_lexer":"r","mimetype":"text/x-r-source","file_extension":".r","version":"4.0.5"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# Reduce Memory Footprint of Competition Data","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19"}},{"cell_type":"markdown","source":"UPDATE: \n\n1. We transpose the output measurement matrices so that the rows correspond to cell_id. For multi ATAC measurement, the columns correspond to genome coordinates. For CITE, the columns correspond to genes.\n\n2. We also output the measure matrices in RDS file so that R users can use directly.\n\n3. The associated dataset is public and can be found [here](https://www.kaggle.com/datasets/stautxie/sparse-measurement-data-open-problems-multimodal).","metadata":{}},{"cell_type":"markdown","source":"The [Open Problems - Multimodal Single-Cell Integration](https://www.kaggle.com/competitions/open-problems-multimodal/) presents a challenge on building effective models to predict interactions of DNA, RNA and proteins in single cells. Even before stepping on modeling, the first challenge is to manage the memory footprint. Particularly, the input files from CITE and Multiome ATAC measures. A naive approach to simply load the underlying data set will consume more than 200 GB RAM, which easily exceeds all possible RAM of most of the computers (mine included).\n\nThe key here is to exploit the sparsity inherent in the measurement. In another word, the interactions between DNA, RNA and proteins are fairly narrow. An individual gene typically regulates a small fraction of the proteins, while a single protein tends to impact the expression levels of a small fraction of genes if any. Therefore, we will be able to save significant amount of memory footprint to save the measurement in sparse matrix format instead of the dense matrix format provided by the competition host.\n\nWhile this [notebook](https://www.kaggle.com/code/sbunzini/reduce-memory-usage-by-95-with-sparse-matrices) presents similar idea, our work presents new contributions:\n\n1. We package the result into reusable h5 format that can be readily used by other users. By doing so, this notebook saves fellow Kagglers' time and carbon output.\n\n2. This work is implemented in R. Since there are much fewer R users in this competition than Python, hopefully this will attract more participation.","metadata":{}},{"cell_type":"code","source":"library(rhdf5)\nlibrary(Matrix)\nlibrary(tictoc)","metadata":{"execution":{"iopub.status.busy":"2022-08-29T19:51:11.346328Z","iopub.execute_input":"2022-08-29T19:51:11.349205Z","iopub.status.idle":"2022-08-29T19:51:12.682510Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The key functionality of this notebook resides within the following function as we convert a dense matrix into a sparse matrix (compressed sparse column array in [scipy](https://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.csc_array.html#scipy.sparse.csc_array) or dgCMatrix in [R](https://cran.r-project.org/web/packages/Matrix/index.html)) in batches. We collect these submatrix through loop. ","metadata":{}},{"cell_type":"code","source":"\nmatrix_to_dgC <- function(h5_file, h5_dat_name, row_index, col_index, chunk_size=5000) {\n  # divide col_index into list of chunk_size and then call selected_matrix_to_dgC and combine results\n  chunk_list <- 0:floor(length(col_index) / chunk_size)\n  chunk_range_list <- lapply(chunk_list, \n                             function(x) {(chunk_size * x + 1):min(length(col_index), (chunk_size * (x+1)))})\n  \n  full_sp_matrix <- Matrix(nrow=length(row_index), ncol=0, sparse=TRUE)\n  for (i in 1:length(chunk_range_list)) {\n    cat(paste(\"iter = \", i, \"\\n\", sep=\"\"))\n    cat(\"Sparsifying ...\\n\")\n    selected_mat <- h5read(h5_file, h5_dat_name, index=list(NULL, chunk_range_list[[i]]))\n    sp_mat <- Matrix(data=selected_mat, sparse=TRUE)\n    rm(selected_mat); gc()\n    \n    cat(\"Binding columns ... \\n\")\n    full_sp_matrix <- cbind(full_sp_matrix, sp_mat)\n    rm(sp_mat); gc()\n  }\n  return (full_sp_matrix)\n}","metadata":{"execution":{"iopub.status.busy":"2022-08-29T19:51:12.684729Z","iopub.execute_input":"2022-08-29T19:51:12.715709Z","iopub.status.idle":"2022-08-29T19:51:12.728832Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The next function is a wrapper so that we collect column and row indices within the original h5 measure file and save them in a new h5 file so that the resultant file is self-sufficient.","metadata":{}},{"cell_type":"code","source":"\nreformat_h5_file <- function(input_scenario) {\n  cat(paste(\"Processing: \", input_scenario, \"\\n\", sep=\"\"))\n  h5ls(file.path(DAT_DIR, paste(input_scenario, \".h5\", sep=\"\")))\n  \n  mu_tr_axis0 <- h5read(file.path(DAT_DIR, paste(input_scenario, \".h5\", sep=\"\")),\n                        paste(\"/\", input_scenario, \"/axis0\", sep=\"\"))\n  \n  mu_tr_axis1 <- h5read(file.path(DAT_DIR, paste(input_scenario, \".h5\", sep=\"\")),\n                        paste(\"/\", input_scenario, \"/axis1\", sep=\"\"))\n  \n  mu_tr_val <- h5read(file.path(DAT_DIR, paste(input_scenario, \".h5\", sep=\"\")),\n                      paste(\"/\", input_scenario, \"/block0_values\", sep=\"\"),\n                      index=list(NULL, 1:100))\n  \n  mu_tr_dgC <- Matrix(data=mu_tr_val, sparse=TRUE) \n  \n  cat(paste(\"Before sparsify: \", format(object.size(mu_tr_val), \"MB\"), \" MB\\n\", sep=\"\"))\n  cat(paste(\"After sparsify: \", format(object.size(mu_tr_dgC), \"MB\"), \" MB\\n\", sep=\"\"))\n  cat(paste(\"Size reduction: \", round(100 *(1- object.size(mu_tr_dgC) / object.size(mu_tr_val)), 2), \"%\\n\", sep=\"\"))\n    \n  rm(mu_tr_dgC, mu_tr_val); gc()\n  \n  \n  full_sp_matrix <- matrix_to_dgC(file.path(DAT_DIR, paste(input_scenario, \".h5\", sep=\"\")),\n                                  paste(\"/\", input_scenario, \"/block0_values\", sep=\"\"),\n                                  mu_tr_axis0,\n                                  mu_tr_axis1,\n                                  chunk_size=5000)\n    \n  full_sp_matrix <- t(full_sp_matrix)\n  gc()\n  dimnames(full_sp_matrix) <- list(mu_tr_axis1, mu_tr_axis0)\n  \n  #~ write into a new h5 file\n  \n  output_h5_file <- paste(\"sp_\", input_scenario, \".h5\", sep=\"\")\n  h5write(mu_tr_axis0, output_h5_file, name=\"axis0\")\n  h5write(mu_tr_axis1, output_h5_file, name=\"axis1\")\n  h5write(full_sp_matrix@i, output_h5_file, name=\"value_i\")\n  h5write(full_sp_matrix@p, output_h5_file, name=\"value_p\")\n  h5write(full_sp_matrix@x, output_h5_file, name=\"value_x\")\n  \n  h5ls(output_h5_file)  \n    \n  #~ save into rds file\n  \n  saveRDS(full_sp_matrix, paste(\"sp_\", input_scenario, \".rds\", sep=\"\"))\n}","metadata":{"execution":{"iopub.status.busy":"2022-08-29T20:00:08.857040Z","iopub.execute_input":"2022-08-29T20:00:08.858743Z","iopub.status.idle":"2022-08-29T20:00:08.874169Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Finally, we loop through all six measurement files and generate the resultant sparse files.","metadata":{}},{"cell_type":"code","source":"\nDAT_DIR <- \"../input/open-problems-multimodal\"\n\n\n#all_scenarios <- c(\"train_multi_inputs\", \"train_multi_targets\",\n#                   \"train_cite_inputs\", \"train_cite_targets\",\n#                   \"test_multi_inputs\", \"test_cite_inputs\")\n\nall_scenarios <- c(\"train_cite_inputs\")\n\n\nfor (sc in all_scenarios) {\n  reformat_h5_file(input_scenario = sc)\n}","metadata":{"execution":{"iopub.status.busy":"2022-08-29T20:00:12.800539Z","iopub.execute_input":"2022-08-29T20:00:12.802167Z","iopub.status.idle":"2022-08-29T20:06:28.313182Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Due to inherent sparsity in the measurement files, sparse measurement files are much smaller:\n\n1. the size of the sparse measurement file is only 4% of the original measurement file in dense matrix format for multiome ATAC experiments and 36% of the CITE measurement. Therefore, it becomes possible to use them on Kaggle server as well as other normal computers, together with statistical modeling approaches that accept sparse matrix. \n\n2. The output h5 file is less than 50% of the original h5 file. ","metadata":{}},{"cell_type":"markdown","source":"Finally a note how to use the data set, which can be found [here](https://www.kaggle.com/datasets/stautxie/sparse-measurement-data-open-problems-multimodal). The new h5 file contains 5 arrays:\n\n1. axis0 (row index from the original h5 file)\n2. axis1 (column index from the original h5 file)\n3. value_i (attribute i in dgCMatrix in R or index indices in csc_array in scipy.sparse\n4. value_p (attribute p in dgCMatrix in R or index indptr in csc_array in scipy.sparse\n5. value_x (attribute x in dgcMatrix in R or index data in csc_array in scipy.sparse.","metadata":{}},{"cell_type":"markdown","source":"Hope this helps. Thanks.","metadata":{}}]}