Main train data
","metadata":{}},{"cell_type":"markdown","source":"From the competition [description](https://www.kaggle.com/competitions/open-problems-single-cell-perturbations/data):\nThe **de_train.parquet** file comprises the main competition data. It contains values for a number of cell_type / sm_name pairs.\n\n- cell_type - The annotated cell type of each cell based on RNA expression.\n- sm_name - The primary name for the (parent) compound (in a standardized representation) as chosen by LINCS. This is provided to map the data in this experiment to the LINCS Connectivity Map data.\n- sm_lincs_id - The global LINCS ID (parent) compound (in a standardized representation). This is provided to map the data in this experiment to the [LINCS Connectivity Map data](https://maayanlab.cloud/Harmonizome/resource/LINCS+L1000+Connectivity+Map).\n- SMILES - Simplified molecular-input line-entry system (SMILES) representations of the compounds used in the experiment. This is a 1D representation of molecular structure. These SMILES are provided by Cellarity based on the specific compounds ordered for this experiment.\n- control - Boolean indicating whether this instance was used as a control.\n- genes A1BG, A1BG-AS1, …, ZZEF1 (numbering 18,211 in all) - Differential expression value **(-log10(p-value) * sign(LFC))** for each gene.","metadata":{}},{"cell_type":"code","source":"de_train <- read_parquet('/kaggle/input/open-problems-single-cell-perturbations/de_train.parquet')\nde_train <- setDT(de_train)\n\nhead(de_train)\ndata.summary(de_train)","metadata":{"execution":{"iopub.status.busy":"2023-11-05T05:45:00.286190Z","iopub.execute_input":"2023-11-05T05:45:00.287221Z","iopub.status.idle":"2023-11-05T05:45:07.707605Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Number of compounds per cell type:","metadata":{}},{"cell_type":"code","source":"de_train[, .N, by = \"cell_type\"]","metadata":{"execution":{"iopub.status.busy":"2023-11-05T05:45:07.709468Z","iopub.execute_input":"2023-11-05T05:45:07.710465Z","iopub.status.idle":"2023-11-05T05:45:07.735676Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Number of rows with data (= number of cell types) per compounds:","metadata":{}},{"cell_type":"code","source":"de_train[, .(num_cell_types = .N), by = \"sm_name\"][\n , .(num_compounds = .N), by = \"num_cell_types\"]","metadata":{"execution":{"iopub.status.busy":"2023-11-05T05:45:07.737934Z","iopub.execute_input":"2023-11-05T05:45:07.738887Z","iopub.status.idle":"2023-11-05T05:45:07.759012Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"cat(\"Drugs with data for 3 cell types only\")\nde_train[, .N, by = \"sm_name\"][N == 3]\ncat(\"\\nDrugs with data for 5 cell types\")\nde_train[, .N, by = \"sm_name\"][N == 5]\ncat(\"\\nDrugs with data for all 6 cell types\")\nde_train[, .N, by = \"sm_name\"][N == 6]","metadata":{"execution":{"iopub.status.busy":"2023-11-05T05:45:07.761966Z","iopub.execute_input":"2023-11-05T05:45:07.763004Z","iopub.status.idle":"2023-11-05T05:45:07.832859Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# - Differential expression (DE)","metadata":{}},{"cell_type":"markdown","source":"Differential expression (DE) values in the data, as indicated in the description, are calculated as (-log10(p-value) * sign(LFC)), where:\n- -log10(p-value) represents the transformation of the p-value into a more interpretable scale\n- sign(LFC) represents the direction of the expression change (positive or negative) \nThus, for each DE value (x) we can back-calculate the p-value (10^x) of the significance of the expression change and determine the direction of the change (sign(x)). Let's see the range of the given DE values.","metadata":{}},{"cell_type":"code","source":"cat(\"For all compounds, cell types, genes:\\n\")\nval_range <- de_train[,\n .(min_value = apply(de_train[, 6:ncol(de_train)], 1, min),\n max_value = apply(de_train[, 6:ncol(de_train)], 1, max))]\n\ncat(\"Min DE value = \", min(val_range[, min_value]), \"\\n\",\n \"Max DE value = \", max(val_range[, max_value]), sep = \"\")","metadata":{"execution":{"iopub.status.busy":"2023-11-05T05:45:07.835395Z","iopub.execute_input":"2023-11-05T05:45:07.836655Z","iopub.status.idle":"2023-11-05T05:45:09.772982Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Multidimensional Scaling (MDS)**","metadata":{}},{"cell_type":"markdown","source":" NMDS (Non-Metric Multidimensional Scaling) is a technique used for reducing the dimensionality (here - to 2 coordinates) and visualizing the dissimilarity between samples.","metadata":{}},{"cell_type":"code","source":"dist_matrix <- dist(as.matrix(de_train[, 6:ncol(de_train)]),\n method = \"euclidean\")\nnmds_result <- metaMDS(dist_matrix)","metadata":{"execution":{"iopub.status.busy":"2023-11-05T05:45:09.774737Z","iopub.execute_input":"2023-11-05T05:45:09.775830Z","iopub.status.idle":"2023-11-05T05:46:02.919811Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig(20, 10)\n\nplot_data <- data.frame(\n NMDS1 = nmds_result$points[, 1],\n NMDS2 = nmds_result$points[, 2],\n Cell_Type = de_train$cell_type,\n Color = factor(de_train$cell_type)\n)\n\ncolors <- rainbow(length(unique(de_train$cell_type)))\n\nggplot(plot_data, aes(x = NMDS1, y = NMDS2)) +\n geom_point(aes(fill = Color), shape = 21, size = 4) +\n ggtitle(\"Multidimensional Scaling (MDS) of DE values using Euclidean distance\") +\n scale_fill_manual(values = colors) +\n theme_bw(base_size = 22) +\n theme(legend.position = \"right\")","metadata":{"execution":{"iopub.status.busy":"2023-11-05T05:46:02.921503Z","iopub.execute_input":"2023-11-05T05:46:02.922553Z","iopub.status.idle":"2023-11-05T05:46:03.622375Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Unaggregated counts
","metadata":{}},{"cell_type":"markdown","source":"The **adata_train.parquet** file - is a supplement to de_train, it conntains unaggregated count and normalized data in COO sparse-array format. Columns:\n- obs_id - This is a unique identifier assigned to each cell in the raw dataset.\n- gene - Corresponds to the columns of de_train.\n- count - The raw molecular counts for the gene expression data measured in the experiment as output by 10x CellRanger.\n- normalized_count - These counts have been library size normalized and log(X+1) transformed.\n\n","metadata":{}},{"cell_type":"code","source":"adata_train <- read_parquet(\"/kaggle/input/open-problems-single-cell-perturbations/adata_train.parquet\",\n as_data_frame = FALSE) %>%\n slice_head(n = 10) %>%\n collect()\nadata_train\nrm(adata_train) ; g <- gc()","metadata":{"execution":{"iopub.status.busy":"2023-10-09T10:04:04.145023Z","iopub.execute_input":"2023-10-09T10:04:04.146441Z","iopub.status.idle":"2023-10-09T10:06:11.949336Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"This dataset is a superset to de_train and contains additional genes. In my PC (the RAM limits on Kagge do not allow to handle this) I read this file, filtered only genes from the de_train, and saved it in qs format. I used this exact code below. ","metadata":{}},{"cell_type":"code","source":"# for(col in c(\"obs_id\", \"gene\", \"count\", \"normalized_count\") ) {\n \n# read_parquet('data/adata_train.parquet',\n# col_select = all_of(col),\n# as_data_frame = FALSE) %>% \n# collect() %>%\n# qsave(paste0(\"adata_train_\", col, \".qs\"))\n# g <- gc()\n# }\n\n# lapply(c(\"obs_id\", \"gene\", \"count\", \"normalized_count\"),\n# function(col) {\n# read_parquet('data/adata_train.parquet',\n# col_select = all_of(col),\n# as_data_frame = FALSE) %>%\n# collect() %>%\n# qsave(paste0(\"data/adata_train_\", col, \".qs\"))\n# g <- gc()\n \n# })\n\n# adata_train <- data.table(\n# \"obs_id\" = qread(\"data/adata_train_obs_id.qs\"),\n# \"gene\" = qread(\"data/adata_train_gene.qs\"),\n# \"count\" = qread(\"data/adata_train_count.qs\")\n# )\n# head(adata_train)\n# setnames(adata_train, names(adata_train), c(\"obs_id\", \"gene\", \"count\"))\n\n# adata_train <- adata_train[gene %chin% names(de_train)]\n# qsave(adata_train, \"adata_train_genes_counts.qs\")\n\n# rm(adata_train) ; g <- gc()\n# # may be helpful: https://hbs-rcs.github.io/large_data_in_R/#read-and-count-lyft-records-with-arrow","metadata":{"execution":{"iopub.status.busy":"2023-10-09T10:06:11.951853Z","iopub.execute_input":"2023-10-09T10:06:11.953161Z","iopub.status.idle":"2023-10-09T10:06:11.964437Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"There are 416,442,312 rows, 240,090 observation ids (~1200-3300 rows per id, not shown), 21255 genes.","metadata":{}},{"cell_type":"markdown","source":"There is also a **adata_obs_meta.csv** file with observation metadata for adata_train:\n- library_id - A unique identifier for each library, which is a measurement made on pooled samples from each row of the plate. All cells from wells on the same row of the same plate will share a - library_id.\n- plate_name - A unique ID for all samples from the same plate.\n- well - The well location of the sample on each plate (this is standard across 96 well plate experiments). It is a concatenation of row and col.\n- row - Which row on the plate the sample came from.\n- col - Which column on the plate the sample came from.\n- donor_id - Identifies the donor source of the sample, one of three.\n- cell_type - The annotated cell type of each cell based on RNA expression. This matches the cell_type in the de_train.parquet.\n- cell_id - This is included for consistency with LINCS Connectivity Map metadata, which denotes a cell_id for each cell line.\n- sm_name - The primary name for the (parent) compound (in a standardized representation) as chosen by LINCS. This is provided to map the data in this experiment to the LINCS Connectivity Map data.\n- sm_lincs_id - The global LINCS ID (parent) compound (in a standardized representation). This is provided to map the data in this experiment to the LINCS Connectivity Map data.\n- SMILES - Simplified molecular-input line-entry system (SMILES) representations of the compounds used in the experiment. This is a 1D representation of molecular structure. These SMILES are provided by Cellarity based on the specific compounds ordered for this experiment.\n- dose_uM - Dose of the compound in on a micro-molar scale. This maps to the pert_idose field in LINCS.\n- timepoint_hr - Duration of treatment in hours. This maps to the pert_itime field in LINCS.\n- control - Whether this observation was used as a control, True or False.","metadata":{}},{"cell_type":"code","source":"adata_obs_meta <- fread(\"/kaggle/input/open-problems-single-cell-perturbations/adata_obs_meta.csv\")\n\nhead(adata_obs_meta)\ndata.summary(adata_obs_meta)","metadata":{"execution":{"iopub.status.busy":"2023-11-05T07:06:24.601097Z","iopub.execute_input":"2023-11-05T07:06:24.603790Z","iopub.status.idle":"2023-11-05T07:06:25.613336Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Plot the experimental setup**","metadata":{}},{"cell_type":"code","source":"plt <- copy(adata_obs_meta)\nplt[, sm_name := stringr::str_wrap(sm_name, width = 10)]\n#plt[nchar(sm_name) > 10 & !grepl(\"\\\\\\n\", sm_name), unique(sm_name)]\n\nplt[sm_name == \"5-(9-Isopropyl-8-methyl-2-morpholino-9H-purin-6-yl)pyrimidin-2-amine\",\n sm_name := \"5-(9-Isopropyl-8-methyl\\n-2-morpholino-9H-purin\\n-6-yl)pyrimidin-2-amine\"]\n\nplt[, n_cell_types := uniqueN(cell_type), by = c(\"well\", \"plate_name\")]\nplt[, n_cells := uniqueN(obs_id), by = c(\"well\", \"plate_name\")]\n\nplt[, info := paste0(sm_name,\n \"\\nn cell types = \", n_cell_types,\n \"\\nn cells = \", n_cells\n )]\nplt <- unique(plt[, .(plate_name, col, row, sm_name, info)])","metadata":{"execution":{"iopub.status.busy":"2023-11-05T07:40:47.349271Z","iopub.execute_input":"2023-11-05T07:40:47.351602Z","iopub.status.idle":"2023-11-05T07:41:00.074087Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig(30, 70)\n\nggplot(plt, aes(x = factor(col), y = row, label = info,\n fill = sm_name)) +\n geom_tile(alpha = 0.3) +\n geom_text(size = 6) +\n facet_wrap(~ plate_name, ncol = 1, scales = \"free_x\") +\n labs(x = \"Column\", y = \"Row\") +\n scale_colour_hue(l = 70, c = 30) + \n theme_bw(base_size = 22) +\n theme(panel.grid = element_blank()) +\n guides(fill = guide_legend(ncol = 1)) #\"none\"","metadata":{"execution":{"iopub.status.busy":"2023-11-05T08:11:54.562099Z","iopub.execute_input":"2023-11-05T08:11:54.564513Z","iopub.status.idle":"2023-11-05T08:12:01.456480Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Number of cells from which the data were collected by cell type**","metadata":{}},{"cell_type":"code","source":"sum_by_ct <- adata_obs_meta[, .N, by = \"cell_type\"][order(-N)]\ncat(\"Number of obs_ids ( = cells) in each cell type\")\nsum_by_ct\n\nsum_by_ct <- adata_obs_meta[, .N, by = \"cell_type\"][order(-N)]\ncat(\"Number of obs_ids ( = cells) in each cell type (cells treated with controls are excluded)\")\n\nadata_obs_meta[!sm_name %in% c(\"Dabrafenib\", \"Belinostat\", \"Dimethyl Sulfoxide\"),\n .N, by = \"cell_type\"][order(-N)] ","metadata":{"execution":{"iopub.status.busy":"2023-10-31T05:43:05.202254Z","iopub.execute_input":"2023-10-31T05:43:05.205563Z","iopub.status.idle":"2023-10-31T05:43:05.800696Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Number of cells from which the data were collected by compounds**","metadata":{}},{"cell_type":"code","source":"sum_by_drug <- adata_obs_meta[, .N, by = \"sm_name\"][order(-N)]\ncat(\"Top componds by number of obs_ids ( = cells)\")\nhead(sum_by_drug, 10)\ncat(\"\\nBottom componds by number of obs_ids ( = cells)\")\ntail(sum_by_drug, 10)","metadata":{"execution":{"iopub.status.busy":"2023-10-15T16:49:26.656889Z","iopub.execute_input":"2023-10-15T16:49:26.658562Z","iopub.status.idle":"2023-10-15T16:49:26.722178Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig(30, 10)\n\nggplot(sum_by_drug, aes(x = sm_name, y = N) ) +\n geom_bar(stat = 'identity', col=\"grey\", fill = \"#BF9039\") + \n xlab(\"\") +\n #facet_wrap(~variable, scales = \"free\", ncol = 4) +\n ggtitle(paste0(\"Number of cell per compond\")) +\n theme_bw(base_size = 14) +\n theme(axis.text.x = element_text(angle = 45, hjust=1),\n panel.grid.major = element_blank()) ","metadata":{"execution":{"iopub.status.busy":"2023-10-15T16:49:30.697166Z","iopub.execute_input":"2023-10-15T16:49:30.698800Z","iopub.status.idle":"2023-10-15T16:49:31.477551Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Multiome (in progress)
","metadata":{}},{"cell_type":"markdown","source":"**multiome_train.parquet** contains optional additional 10x Multiome data for each sample at baseline.\n\n- obs_id - Unique identifier for each observation. (Distinct from identifiers used in adata.)\n- location - This is a feature ID. If the feature_type in multiome_var_meta.csv is Gene Expression then this is a gene symbol. If feature_type is Peaks, then this is the genomic interval of the peak.\n- count - This is the raw molecular counts of the transcript for accessible DNA measurement as output by Cellranger-Arc.\n- normalized_count - If the feature_type in multiome_var_meta.csv is Gene Expression then this is library size normalized and log(X+1) transformed counts. If feature_type is Peaks, then this is ATAC-seq peak counts transformed with TF-IDF using the default log(TF) * log(IDF).","metadata":{}},{"cell_type":"code","source":"multiome_train <- read_parquet('/kaggle/input/open-problems-single-cell-perturbations/multiome_train.parquet')\nmultiome_train <- setDT(multiome_train)\n\nhead(multiome_train)\ndata.summary(multiome_train)\nrm(multiome_train) ; gc()","metadata":{"execution":{"iopub.status.busy":"2023-10-09T10:31:40.374063Z","iopub.execute_input":"2023-10-09T10:31:40.376598Z","iopub.status.idle":"2023-10-09T10:34:57.832584Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**multiome_obs_meta.csv**:\n\n- obs_id - Identifier corresponding to that in multiome_train.parquet.\n- cell_type - The annotated cell type of each cell based on RNA expression.\n- donor_id - Identifies the donor source of the sample, one of three.","metadata":{}},{"cell_type":"code","source":"multiome_obs_meta <- fread(\"/kaggle/input/open-problems-single-cell-perturbations/multiome_obs_meta.csv\")\n\nhead(multiome_obs_meta)\ndata.summary(multiome_obs_meta)","metadata":{"execution":{"iopub.status.busy":"2023-10-09T10:36:41.096904Z","iopub.execute_input":"2023-10-09T10:36:41.098239Z","iopub.status.idle":"2023-10-09T10:36:41.177817Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**multiome_var_meta.csv**:\n\n- location - This is a feature ID. If the feature_type is Gene Expression then this is a gene symbol. If feature_type is Peaks, then this is the genomic interval of the peak.\n- gene_id - This is an alternative unique feature ID. If the feature_type is Gene Expression then this is an Ensembl Stable Gene ID. If feature_type is Peaks, then this is the genomic interval of the peak.\n- feature_type - Denotes whether the feature is an RNA expression measurement or a Chromatin Accessibility measurement.\n- genome - The genome version used when running CellRanger-Arc\n- interval - The genomic coordinates of each feature on reference genome GRCh38. Genomic coordinates are directly related to the reference genome and include the chromosome name, start position, - and end position in the following format: chr1:1234570-1234870.","metadata":{}},{"cell_type":"code","source":"multiome_var_meta <- fread(\"/kaggle/input/open-problems-single-cell-perturbations/multiome_var_meta.csv\")\n\nhead(multiome_var_meta)\ndata.summary(multiome_var_meta)","metadata":{"execution":{"iopub.status.busy":"2023-10-09T10:36:41.222037Z","iopub.execute_input":"2023-10-09T10:36:41.223492Z","iopub.status.idle":"2023-10-09T10:36:41.796043Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# multiome_var_meta <- multiome_var_meta[feature_type == \"Gene Expression\" &\n# location %chin% names(de_train)]\n# multiome_var_meta[, uniqueN(location)]\n\n# multiome_train <- multiome_train[location %chin% names(de_train)]\n# multiome_train[, uniqueN(location)]","metadata":{"execution":{"iopub.status.busy":"2023-10-09T10:37:13.423446Z","iopub.execute_input":"2023-10-09T10:37:13.425444Z","iopub.status.idle":"2023-10-09T10:37:13.481885Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# length(setdiff(names(de_train),\n# multiome_train[, unique(location)]))","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#Test data
","metadata":{}},{"cell_type":"markdown","source":"Our goal is to predict corresponding differential expression value (-log10(p-value) * sign(LFC)) for the cell_type / sm_name pairs given in id_map.csv. The p-value is a score for the perturbation effect output from a linear model used by competition hosts (Limma), that takes into account various factors such as gene expression, experimental covariates, and Bayesian priors. So, it apppears, we have to predict the statistical significance of the change, and its direction, not the DE itself.","metadata":{}},{"cell_type":"code","source":"id_map <- fread(\"/kaggle/input/open-problems-single-cell-perturbations/id_map.csv\")\n\nhead(id_map)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"data.summary(id_map)\n\ncat(\"\\nCell types in the test set:\", paste0(unique(id_map[, cell_type]), collapse = \", \"))","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"So, there are 129 different compounds in the test set. We have to predict their impact on gene expression in cells of 2 types: B cells and Myeloid cells.\n- Public Test includes 50 randomly selected compounds\n- Private Test - 79 randomly selected compounds","metadata":{}},{"cell_type":"markdown","source":"Number of compounds per cell type:","metadata":{}},{"cell_type":"code","source":"id_map[, .N, by = \"cell_type\"]","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The input to your model will be a tuple of cell_type and sm_name and the output of your model will be predicted signed -log10(p-values) for all 18211 genes.","metadata":{}},{"cell_type":"code","source":"#","metadata":{"trusted":true},"execution_count":null,"outputs":[]}]}