{"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":"# Protein vs top 100 correlated RNA: intersections & expression","metadata":{}},{"cell_type":"markdown","source":"This is the analysis of \n- within (fixed day and donor) and between cell type intersections of top RNAs correlated with the protein of interest (tables with correlations are the output of [this notebook](http://https://www.kaggle.com/code/antoninadolgorukova/mmscel-protein-analysis))\n- expression levels of several RNA that are \"stable\" within EryP\n- example of dropouts analysis\n\nData: RDS files with sparse matrices of [normalised counts data](http://www.kaggle.com/datasets/stautxie/sparse-measurement-data-open-problems-multimodal) for [Open Problems - Multimodal Single-Cell Integration](http://www.kaggle.com/competitions/open-problems-multimodal). The dataset for this competition comprises single-cell multiomics data collected from mobilized peripheral CD34+ hematopoietic stem and progenitor cells (HSPCs) isolated from four healthy human donors. The train data we use here consists of both gene expression (RNA) and surface protein data for days 2,3,4 for donors 1-3 (donor IDs: 32606,13176, and 31800).","metadata":{}},{"cell_type":"code","source":"# Libraries and some functions\nsuppressMessages(library(dplyr)) \nsuppressMessages(library(tidyr))\nsuppressMessages(library(gtools))\nsuppressMessages(library(rstatix))\nsuppressMessages(library(purrr))\nsuppressMessages(library(Matrix))\nsuppressMessages(library(tictoc))\nsuppressMessages(library(ggplot2))\nsuppressMessages(library(openxlsx))\nsuppressMessages(library(RVenn)) #package for set operations on multiple sets\nlibrary(UpSetR) #intersection visualisation\nsuppressMessages(remotes::install_github(\"jokergoo/ComplexHeatmap\")) #takes time to install\nlibrary(ComplexHeatmap)\nlibrary(superheat)\n\n# function for figure size adjusment\nfig <- function(width, heigth) {\n    options(repr.plot.width = width, repr.plot.height = heigth)\n}","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-01-16T08:51:55.421743Z","iopub.execute_input":"2023-01-16T08:51:55.424181Z","iopub.status.idle":"2023-01-16T08:55:26.303980Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Load proteins and RNA data**","metadata":{}},{"cell_type":"markdown","source":"Proteins (Targets dataset): for the surface protein levels, each row corresponds to a cell (e.g. \"45006fe3e4c8\") and each column to a protein (e.g. \"CD86\").\n\nRNA (Inputs dataset): For the RNA counts, each row corresponds to a cell (e.g. \"45006fe3e4c8\") and each column to a gene. The column format for a gene is given by {EnsemblID}_{GeneName} where EnsemblID refers to the Ensembl Gene ID and GeneName to the gene name (e.g. \"ENSG00000159840_ZYX\").\n\nMetadata: Donor and cell types. The train data consists of both gene expression (RNA) and surface protein data for days 2,3,4 for donors 1-3 (donor IDs: 32606,13176, and 31800), the public test data consists of RNA for days 2,3,4 for donor 4 (donor ID: 27678) and the private test data consists data from day 7 from all donors.","metadata":{}},{"cell_type":"code","source":"tic() # ~ 1 min\n# Load Metadata (day, donor, cell type, and technology)\nmetadata <- read.csv('../input/open-problems-multimodal/metadata.csv',row.names=1)\nmetadata <- filter(metadata, technology == \"citeseq\" &\n                   donor %in% c(\"13176\", \"31800\", \"32606\") &\n                   day %in% c(\"2\",\"3\",\"4\")) %>%\n  mutate(\"Row.names\" = row.names(.))\n# Load protein data\n# raw data\n#path <- \"/kaggle/input/sparse-raw-counts-data-open-problems-multimodal/citeseq/sp_train_cite_targets_raw.rds\"\n# normalized data\npath <- \"/kaggle/input/sparse-measurement-data-open-problems-multimodal/sp_train_cite_targets.rds\"\nmat_prot <- readRDS(path)\n\n#dgCMatrix to matrix\nmat_prot <- as.matrix(mat_prot)\n\ntoc()\ngc()","metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"execution":{"iopub.status.busy":"2023-01-16T08:56:33.364758Z","iopub.execute_input":"2023-01-16T08:56:33.415853Z","iopub.status.idle":"2023-01-16T08:56:38.728770Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Load Spearman correlations**","metadata":{}},{"cell_type":"code","source":"# Load correlations data\npath <- \"/kaggle/input/mmscel-protein-analysis/Prot-RNA_corr_all.csv\"\nall_dat <- read.csv(path)\nall_dat <- all_dat %>% select(RNA, everything())\n\ncat(\"The head of the data frame with Spearman correlations between\",\n    ncol(all_dat)-1, \"proteins and\",nrow(all_dat),\n    \"RNA in all cells\")\nhead(all_dat)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-01-16T09:05:40.016977Z","iopub.execute_input":"2023-01-16T09:05:40.018874Z","iopub.status.idle":"2023-01-16T09:05:47.450971Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Select the protein of interest**","metadata":{}},{"cell_type":"code","source":"prot_of_interest <- \"CD36\"\nmy_prot <- mat_prot[, prot_of_interest]","metadata":{"execution":{"iopub.status.busy":"2023-01-16T09:05:51.055503Z","iopub.execute_input":"2023-01-16T09:05:51.057021Z","iopub.status.idle":"2023-01-16T09:05:51.075649Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# function for the selection of most pronounced correlations by cutoff or the top with n correlations\n\n# data - a df with correlation coefficients for protein-protein|RNA pairs,\n# calculated for subgroups by donor, day or cell type (names in group).\n# factor - \"donor\", \"day\" or \"cell type\",\n# cut_off if a cutoff for correlation coefficients\n# n - max number of pairs to show\n\nprep_df <- function(data, fact, cut_off = -1, n) {\n    dat <- filter(data, factor == fact & abs(corr.coeff) > cut_off)\n\n    all_tmp <- NULL\n    for (f in unique(dat$group)) {\n\n        tmp <- dat[dat$group == f, ]\n        tmp <- tmp[order(abs(tmp$corr.coeff), decreasing = TRUE), ]\n        tmp <- head(tmp, n = n)\n        all_tmp <- rbind(all_tmp, tmp)\n    }\n    dat <- all_tmp\n    return(dat)\n}","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-01-16T09:05:53.644317Z","iopub.execute_input":"2023-01-16T09:05:53.645891Z","iopub.status.idle":"2023-01-16T09:05:53.660724Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# 1. Between cell types","metadata":{}},{"cell_type":"markdown","source":"Let's first visualize intersections between cell types.","metadata":{}},{"cell_type":"code","source":"# Load correlations data by cell type\ntic() # ~2 min\npath <- \"/kaggle/input/mmscel-protein-analysis/Prot-RNA_corr_13gr.csv\"\nall_dat <- read.csv(path)\nall_dat <- all_dat %>% select(factor, group, sample_size, RNA, everything())\ntoc()\n\ncat(\"The head of the data frame with Spearman correlations between\",\n    ncol(all_dat)-4, \"proteins and\",nrow(all_dat)/13,\n    \"RNA in each donor, day, cell type (13 groups)\")\nhead(all_dat)\n\n# select the protein of interest\ndat <- all_dat %>% select(factor, group, sample_size, RNA, all_of(prot_of_interest))\nnames(dat)[names(dat) == prot_of_interest] <- \"corr.coeff\"","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-01-16T09:05:57.253093Z","iopub.execute_input":"2023-01-16T09:05:57.254743Z","iopub.status.idle":"2023-01-16T09:10:12.973071Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#select top 100 RNA correlated with the protein of interest\nfact = \"cell_type\"\ncut_off = -1\nn = 100\n\ndat <- prep_df(dat, fact, cut_off, n)\n#cat(\"Head of the resulting data frame\")\n#head(dat)\n\ncat(\"\\nSample sizes (number of cells) in the datasets with top correlations by cell type\")\ndat %>%\n  group_by(group, sample_size) %>%\n  summarise(n.of.corrs = sum(!is.na(corr.coeff)), .groups = \"keep\") %>% \n  pivot_wider(names_from = \"group\", values_from = \"sample_size\")","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-01-16T09:13:51.641983Z","iopub.execute_input":"2023-01-16T09:13:51.643483Z","iopub.status.idle":"2023-01-16T09:13:57.634083Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#calculate expression levels by cell_type\n\nout_tab_prot <- NULL\nfor (ct in unique(metadata$cell_type)) {\n    \n    meta <- filter(metadata, cell_type == ct)\n\n    prot <- my_prot[names(my_prot) %in% meta$Row.names]\n    prot_expr <- mean(prot)\n    \n    tmp_prot <- data.frame(\"cell_type\" = ct,\n                           \"mean\" = prot_expr)\n    out_tab_prot <- rbind(out_tab_prot, tmp_prot)\n}\n\nout_tab_prot <- out_tab_prot %>%\n  pivot_wider(names_from = \"cell_type\", values_from = \"mean\")\n\nout_tab_prot <- cbind(\"ALL\" = mean(my_prot),out_tab_prot)\n\ncat(\"Mean expression levels of\",prot_of_interest, \": in all data and  by cell type\")\nout_tab_prot","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-01-16T09:14:01.228790Z","iopub.execute_input":"2023-01-16T09:14:01.230387Z","iopub.status.idle":"2023-01-16T09:14:01.379796Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"res <- \n  dat[order(dat$corr.coeff, decreasing = TRUE), ] %>%\n  group_by(group) %>%\n  summarize(RNA = paste0(gsub(\".*_\",\"\",RNA),\" (\", round(corr.coeff,2),\")\"),\n           .groups = \"keep\") %>%\n  pivot_wider(names_from = \"group\", values_from = RNA, values_fn = list)\n\nt1 <- gsub(\"\\\\s*\\\\([^\\\\)]+\\\\)\",\"\",unlist(res$`BP`))\nt2 <- gsub(\"\\\\s*\\\\([^\\\\)]+\\\\)\",\"\",unlist(res$`EryP`))\nt3 <- gsub(\"\\\\s*\\\\([^\\\\)]+\\\\)\",\"\",unlist(res$`HSC`))\nt4 <- gsub(\"\\\\s*\\\\([^\\\\)]+\\\\)\",\"\",unlist(res$`MasP`))\nt5 <- gsub(\"\\\\s*\\\\([^\\\\)]+\\\\)\",\"\",unlist(res$`MkP`))\nt6 <- gsub(\"\\\\s*\\\\([^\\\\)]+\\\\)\",\"\",unlist(res$`MoP`))\nt7 <- gsub(\"\\\\s*\\\\([^\\\\)]+\\\\)\",\"\",unlist(res$`NeuP`))\n\nres$`In all cells` <- t1[t1 %in% t2 & t1 %in% t3\n                       & t1 %in% t4\n                       & t1 %in% t5\n                       & t1 %in% t6\n                       & t1 %in% t7] %>% list\ncat(\"Lists with top\", n, \"RNA correlated with\", prot_of_interest, \"by cell type,\",\n   \"sorted by the value of Spearman correlation coefficient (in brackets)\")\nres\ncat(length(unlist(res$`In all cells`)), \"RNA are in the top 100 in each cell type\")","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-01-16T09:14:05.437146Z","iopub.execute_input":"2023-01-16T09:14:05.439147Z","iopub.status.idle":"2023-01-16T09:14:05.562332Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Let's try different visualisation techniques to get an idea about intersections of top correlated RNA in the data sets: \n- [RVenn](https://cran.r-project.org/web/packages/RVenn/vignettes/vignette.html) - a package for dealing with multiple sets\n- [UpSetR](https://cran.r-project.org/web/packages/UpSetR/vignettes/basic.usage.html) - a More Scalable Alternative to Venn and Euler Diagrams for Visualizing Intersecting Sets or [ComplexHeatmap package](https://jokergoo.github.io/ComplexHeatmap-reference/book/upset-plot.html) - to visualize associations between different sources of data sets and reveal potential patterns.\n- [superheat](http://rlbarter.github.io/superheat/index.html) - customizable and extendable heatmaps which act as a tool for the visual exploration of complex datasets.","metadata":{}},{"cell_type":"markdown","source":"**Setmap (presence/absence of RNAs in cell types)**","metadata":{}},{"cell_type":"markdown","source":"The union of the 7 data sets (top 100 RNAs for each cell type = 700 RNA): 508, thus, we take top 20 here.","metadata":{}},{"cell_type":"code","source":"fig(12,22)\nintr = map(unique(dat$group),\n          function(x) {\n              dat %>% filter(group == x) %>% \n              arrange(desc(abs(corr.coeff))) %>% \n              head(n = 20) %>%\n              select(RNA) %>% unlist\n              } )\nnames(intr) <- unique(dat$group)   \n#str(intr)\nintrVenn <- Venn(intr)\n\ncat(\"The union of the top 20 RNAs correlated with\",prot_of_interest, \"in the 7 cell types:\", unite(intrVenn) %>% length)\nsetmap(intrVenn, title = \"A clustered heatmap showing presence/absence of the RNA in cell types\",\n      set_fontsize = 18, element_fontsize = 14)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-01-16T09:14:10.970387Z","iopub.execute_input":"2023-01-16T09:14:10.972373Z","iopub.status.idle":"2023-01-16T09:14:11.996078Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Upset plot (intersection sizes)**","metadata":{}},{"cell_type":"markdown","source":"With upset plot we can see intersections in complete top 100 lists.","metadata":{}},{"cell_type":"code","source":"intr = map(unique(dat$group),\n          function(x) {\n              dat %>% filter(group == x) %>% \n              arrange(desc(corr.coeff)) %>% \n              head(n = 100) %>%\n              select(RNA) %>% unlist\n              } )\nnames(intr) <- unique(dat$group)\n\n#upset(fromList(intr), nsets = 7, text.scale = 2.5, nintersects = NA) # UpSetR package\nintr <- make_comb_mat(intr)\nfig(15,5)\nUpSet(intr, top_annotation = upset_top_annotation(intr, add_numbers = TRUE, \n                                                  numbers_gp = gpar(fontsize = 14),\n                                                  #axis_param = list(gp = gpar(fontsize = 16)),\n                                                  ),\n     pt_size = unit(5, \"mm\"), lwd = 2,\n      row_names_gp = gpar(fontsize = 18)) ","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-01-16T09:14:16.128682Z","iopub.execute_input":"2023-01-16T09:14:16.130572Z","iopub.status.idle":"2023-01-16T09:14:17.725595Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Heatmap**","metadata":{}},{"cell_type":"markdown","source":"To see more closely RNA that have > 2 intersections across cell types we plot a heatmap.","metadata":{}},{"cell_type":"code","source":"# Find intersections\nintr <- dat %>% \n  select(RNA, c(group, corr.coeff)) %>%\n  pivot_wider(names_from = \"group\", values_from = \"corr.coeff\") %>%\n  mutate(intersections = rowSums(!is.na(.[-1])) ) %>%\n  filter(intersections > 2)\nRNA_res <- as.matrix(select(intr, -c(1, ncol(intr))) )\nrownames(RNA_res) <- intr$RNA\n\n#to add corr coefficients in the plot\nRNA_txt <- round(RNA_res, 2)\nRNA_txt[is.na(RNA_txt)] <- \"\"\n\n#plot heatmap with superheat package https://rlbarter.github.io/superheat/index.html\nfig(18,12)\nsuperheat(RNA_res, heat.na.col = \"white\", heat.pal = c(\"#1E90FF\", \"#FF4500\"), #heat.col.scheme = \"red\", \n         left.label.col = \"white\", bottom.label.col = \"white\",\n         legend.height = 0.1, legend.text.size = 16,\n         X.text = RNA_txt, \n         title = paste(\"RNA correlated with\",prot_of_interest, \"with > 2 intersections in top 100 lists for cell types\"),title.size = 8 )","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-01-16T09:14:21.264515Z","iopub.execute_input":"2023-01-16T09:14:21.266637Z","iopub.status.idle":"2023-01-16T09:14:38.247280Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"for (thing in ls()) { message(thing); print(object.size(get(thing)), units='auto') }\nrm(all_dat,dat, ct,n,tmp_prot,out_tab_prot, RNA_res, RNA_txt, intr,intrVenn, fact, meta, cut_off,\n  prot, prot_expr,res, t1,t2,t3,t4,t5,t6,t7)\ngc()","metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Stable all 3 days","metadata":{}},{"cell_type":"markdown","source":"Next, let's analyze the lists of the top 100 RNA for each cell type at each day (300 RNA per cell type).","metadata":{}},{"cell_type":"code","source":"# Load correlations data\n# by cell type at each day\ntic() # ~3 min\npath <- \"/kaggle/input/mmscel-protein-analysis/Prot-RNA_corr_21gr.csv\"\nall_dat <- read.csv(path)\nall_dat <- all_dat %>% select(cell_type, day, sample_size, RNA, everything())\ntoc()\n\ncat(\"The head of the data frame with Spearman correlations between\",\n    ncol(all_dat)-4, \"proteins and\",nrow(all_dat)/21,\n    \"RNA in each cell type by donor, day (21 groups)\")\nhead(all_dat)\n\n# select the protein of interest\ndat <- all_dat %>% select(cell_type, day, sample_size, RNA, all_of(prot_of_interest))\nnames(dat)[names(dat) == prot_of_interest] <- \"corr.coeff\"\ndat$factor <- \"cell_type\"\ndat$group <- interaction(dat$cell_type, dat$day)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-01-16T09:16:13.156573Z","iopub.execute_input":"2023-01-16T09:16:13.158010Z","iopub.status.idle":"2023-01-16T09:20:59.258972Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#select top 100 RNA correlated with the protein of interest\nfact = \"cell_type\"\ncut_off = -1\nn = 100\n\ndat <- prep_df(dat, fact, cut_off, n)\n#cat(\"Head of the resulting data frame\")\n#head(dat)\ncat(\"\\nSample sizes (number of cells) in the datasets with top\", n, \"correlations by cell type at each day\")\ndat %>%\n  group_by(cell_type, day, sample_size) %>%\n  summarise(n.of.corrs = sum(!is.na(corr.coeff)), .groups = \"keep\") %>% \n  pivot_wider(names_from = \"cell_type\", values_from = \"sample_size\")","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-01-16T09:27:46.482120Z","iopub.execute_input":"2023-01-16T09:27:46.483803Z","iopub.status.idle":"2023-01-16T09:27:51.903369Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#calculate expression levels by cell_type at each day\n\nout_tab_prot <- NULL\nfor (ct in unique(metadata$cell_type)) {\n    for (da in unique(metadata$day)) {\n        \n        meta <- filter(metadata, cell_type == ct & day == da)\n\n        prot <- my_prot[names(my_prot) %in% meta$Row.names]\n        prot_expr <- mean(prot)\n\n        tmp_prot <- data.frame(\"cell_type\" = ct,\n                               \"day\" = da,\n                               \"mean\" = prot_expr)\n        out_tab_prot <- rbind(out_tab_prot, tmp_prot)\n        \n    }\n}\n\nout_tab_prot <- out_tab_prot %>%\n  pivot_wider(names_from = \"cell_type\", values_from = \"mean\")\n\n#calculate expression levels for all data at each day\nout_prot <- NULL\nfor(gr in unique(metadata$day)) {\n    meta <- filter(metadata, day == gr)\n\n    prot <- my_prot[names(my_prot) %in% meta$Row.names]\n    prot_expr <- mean(prot)\n    \n    tmp_prot <- data.frame(\"day\" = gr,\n                           \"ALL\" = prot_expr)\n    out_prot <- rbind(out_prot, tmp_prot)  \n}\n\nout_tab_prot <- cbind(out_prot, out_tab_prot[-1])\ncat(\"Mean expression levels of\",prot_of_interest, \"at each day: in all data and  by cell type\")\nout_tab_prot","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-01-16T09:27:55.559697Z","iopub.execute_input":"2023-01-16T09:27:55.561616Z","iopub.status.idle":"2023-01-16T09:28:01.306636Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#filter only RNA that ae in the top 100 lists for all 3 days\nstable_3_days <- group_by(dat, cell_type, RNA) %>%\n  select(-c(sample_size, group) ) %>%\n  pivot_wider(names_from = \"day\", values_from = \"corr.coeff\") %>% ungroup\n\nstable_3_days$total <- rowSums(!is.na(stable_3_days[4:6])) \nstable_3_days <- filter(stable_3_days, total == 3 )\nstable_3_days$mean <- rowMeans(stable_3_days[4:6])\nstable_3_days <- stable_3_days[order(stable_3_days$mean, decreasing = TRUE), ]\n\nres <- stable_3_days %>% \n  mutate(RNA_corr = paste0(gsub(\".*_\",\"\",RNA),\" (\", round(mean,2),\")\")) %>% \n  select(cell_type, RNA_corr) %>%\n  pivot_wider(names_from = \"cell_type\", values_from = RNA_corr, values_fn = list)\n\nt1 <- gsub(\"\\\\s*\\\\([^\\\\)]+\\\\)\",\"\",unlist(res$`MoP`))\nt2 <- gsub(\"\\\\s*\\\\([^\\\\)]+\\\\)\",\"\",unlist(res$`EryP`))\nt3 <- gsub(\"\\\\s*\\\\([^\\\\)]+\\\\)\",\"\",unlist(res$`HSC`))\nt4 <- gsub(\"\\\\s*\\\\([^\\\\)]+\\\\)\",\"\",unlist(res$`MasP`))\nt5 <- gsub(\"\\\\s*\\\\([^\\\\)]+\\\\)\",\"\",unlist(res$`MkP`))\n#t6 <- gsub(\"\\\\s*\\\\([^\\\\)]+\\\\)\",\"\",unlist(res$`BP`))\nt7 <- gsub(\"\\\\s*\\\\([^\\\\)]+\\\\)\",\"\",unlist(res$`NeuP`))\n\nres$`In all cells` <- t1[t1 %in% t2 & t1 %in% t3\n                       & t1 %in% t4\n                       & t1 %in% t5\n                       #& t1 %in% t6\n                       & t1 %in% t7] %>% list\ncat(\"RNA form the top\", n, \"lists that correlates with\", prot_of_interest, \"all 3 days by cell type,\",\n   \"sorted by the value of the mean Spearman correlation coefficient (in brackets)\")\nres\ncat(length(unlist(res$`In all cells`)), \"RNA are in the top 100 in each cell type\")","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-01-16T09:28:14.854404Z","iopub.execute_input":"2023-01-16T09:28:14.856392Z","iopub.status.idle":"2023-01-16T09:28:14.998290Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"intr <- list()\nfor (gr in unique(stable_3_days$cell_type) ) {\n    tmp <- stable_3_days %>% filter(cell_type == gr) %>% \n              arrange(desc(mean)) %>% \n              #head(n = 10) %>%\n              select(RNA) %>% as.list\n    names(tmp) = gr\n    intr <- append(intr, tmp)\n}\ncat(\"The lists of the RNA that correlates with CD36 all 3 days by cell type:\\n\")\nstr(intr)\nintr_Venn <- Venn(intr)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-01-16T09:28:21.444525Z","iopub.execute_input":"2023-01-16T09:28:21.446680Z","iopub.status.idle":"2023-01-16T09:28:21.583609Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"cat(\"The union of the RNAs that correlate with\", prot_of_interest, \"all 3 days:\", unite(intr_Venn) %>% length)\nfig(12,22)\nsetmap(intr_Venn)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-01-16T09:28:24.912082Z","iopub.execute_input":"2023-01-16T09:28:24.914175Z","iopub.status.idle":"2023-01-16T09:28:25.323828Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"intr_UpSet <- make_comb_mat(intr)\nfig(15,5)\nUpSet(intr_UpSet, top_annotation = upset_top_annotation(intr_UpSet, add_numbers = TRUE,\n                                                 numbers_gp = gpar(fontsize = 14)),\n     pt_size = unit(5, \"mm\"), lwd = 2,\n     row_names_gp = gpar(fontsize = 18))","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-01-16T09:28:29.583077Z","iopub.execute_input":"2023-01-16T09:28:29.584753Z","iopub.status.idle":"2023-01-16T09:28:32.892181Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Find intersections\nintr_heat <- pivot_wider(select(stable_3_days, c(RNA, cell_type, mean)), names_from = \"cell_type\", values_from = \"mean\") %>%\nmutate(intersections = rowSums(!is.na(.[-(1)])) ) %>% filter(intersections > 1)\n\nres <- as.matrix(select(intr_heat, -c(1, ncol(intr_heat))) )\nrownames(res) <- intr_heat$RNA\n\n#to add corr coefficients in the plot\nres_txt <- round(res, 2)\nres_txt[is.na(res_txt)] <- \"\"\n\n#plot heatmap with superheat package https://rlbarter.github.io/superheat/index.html\nfig(18,5)\ncat(\"Mean Spearman correlation coefficients for RNAs correlated with\", prot_of_interest,\n    \"\\nthat are in the top 100 for 3 days and overlap between > 1 of 7 cell types.\")\nsuperheat(res, heat.na.col = \"white\", heat.pal = c(\"#1E90FF\", \"#FF4500\"), #, heat.col.scheme = \"red\",\n         left.label.col = \"white\", bottom.label.col = \"white\",\n         legend.height = 0.15, legend.text.size = 16,\n         X.text = res_txt,\n         title.size = 8 )","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-01-16T09:28:36.253106Z","iopub.execute_input":"2023-01-16T09:28:36.254950Z","iopub.status.idle":"2023-01-16T09:28:37.971137Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"for (thing in ls()) { message(thing); print(object.size(get(thing)), units='auto') }\nrm(all_dat,dat, ct,da,gr,n,tmp_prot,tmp,out_prot, out_tab_prot,intr,intr_heat,intr_UpSet,intr_Venn, fact,meta, cut_off,\n  stable_3_days,prot, prot_expr,res,res_txt, t1,t2,t3,t4,t5,t6,t7)\ngc()","metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# 2. Within cell types: fixed day and donor","metadata":{}},{"cell_type":"markdown","source":"Next, we analyze the lists of the top 100 RNA for each cell type at each donor and day (900 RNA per 7 cell types = 63 groups).","metadata":{}},{"cell_type":"code","source":"# Load correlations data by cell type at each day\ntic() # ~20 min\npath <- \"/kaggle/input/mmscel-protein-analysis/Prot-RNA_corr_63gr.csv\"\nall_dat <- read.csv(path)\nall_dat <- all_dat %>% select(cell_type, day, donor, sample_size, RNA, everything())\ntoc()\n\ncat(\"The head of the data frame with Spearman correlations between\",\n    ncol(all_dat)-5, \"proteins and\",nrow(all_dat)/62,\n    \"RNA in each cell type by donor, day (63 groups)\")\nhead(all_dat)\n\n# select the protein of interest\ndat <- all_dat %>% select(cell_type, day, donor, sample_size, RNA, all_of(prot_of_interest))\nnames(dat)[names(dat) == prot_of_interest] <- \"corr.coeff\"\ndat$factor <- \"cell_type\"\ndat$group <- interaction(dat$cell_type, dat$day, dat$donor)\ndat$Protein <- prot_of_interest","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-01-16T09:32:06.757092Z","iopub.execute_input":"2023-01-16T09:32:06.758567Z","iopub.status.idle":"2023-01-16T09:49:09.165994Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#select top 100 RNA correlated with the protein of interest\nfact = \"cell_type\"\ncut_off = -1\nn = 100\n\ndat <- prep_df(dat, fact, cut_off, n)\n#cat(\"Head of the resulting data frame\")\n#head(dat)\ncat(\"\\nSample sizes (number of cells) in the datasets with top\", n, \"correlations by cell type at each day in each donor\")\ndat %>%\n  group_by(cell_type, day, donor, sample_size) %>%\n  summarise(n.of.corrs = sum(!is.na(corr.coeff)), .groups = \"keep\") %>% \n  pivot_wider(names_from = \"cell_type\", values_from = \"sample_size\") %>%\n  arrange(day)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-01-16T10:08:08.772884Z","iopub.execute_input":"2023-01-16T10:08:08.774487Z","iopub.status.idle":"2023-01-16T10:09:05.941377Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#calculate expression levels by cell_type at each day in each donor\nout_tab_prot <- NULL\nfor (ct in unique(metadata$cell_type)) {\n    for (da in unique(metadata$day)) {\n        for (do in unique(metadata$donor)) {\n        \n        meta <- filter(metadata, cell_type == ct & day == da & donor == do)\n\n        prot <- my_prot[names(my_prot) %in% meta$Row.names]\n        prot_expr <- mean(prot)\n\n        tmp_prot <- data.frame(\"cell_type\" = ct,\n                               \"day\" = da,\n                               \"donor\" = do,\n                               \"mean\" = prot_expr)\n        out_tab_prot <- rbind(out_tab_prot, tmp_prot)\n        }\n    }\n}\n\nout_tab_prot <- out_tab_prot %>%\n  pivot_wider(names_from = \"cell_type\", values_from = \"mean\")\n\n#calculate expression levels for all data at each day\nout_prot <- NULL\nfor(da in unique(metadata$day)) {\n    for(do in unique(metadata$donor)) {\n        meta <- filter(metadata, day == da & donor == do)\n\n        prot <- my_prot[names(my_prot) %in% meta$Row.names]\n        prot_expr <- mean(prot)\n\n        tmp_prot <- data.frame(\"day\" = da,\n                               \"donor\" = do,\n                               \"ALL\" = prot_expr)\n        out_prot <- rbind(out_prot, tmp_prot)\n        }\n}\n\nout_tab_prot <- merge(out_prot, out_tab_prot, by = c(\"day\", \"donor\"))\ncat(\"Mean expression levels of\",prot_of_interest, \"at each day in each donor: in all data and  by cell type\")\nout_tab_prot","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-01-16T10:09:16.940294Z","iopub.execute_input":"2023-01-16T10:09:16.941780Z","iopub.status.idle":"2023-01-16T10:09:18.471853Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Here we look for intersections of the top 100 RNAs correlated with the protein of interest within cell types (same RNA is in the top 100 for each donor and day).","metadata":{}},{"cell_type":"markdown","source":"## EryP","metadata":{}},{"cell_type":"code","source":"eryp <- filter(dat, cell_type == \"EryP\") %>%\n      mutate(set = paste0(donor,\"_\", day)) %>%\n      select(-c(sample_size, donor, day, cell_type, factor, group)) %>%\n      pivot_wider(names_from = \"set\", values_from = \"corr.coeff\") %>%\n      as.data.frame %>%\n      mutate(intersections = rowSums(!is.na(.[-(1:2)])) )\n\ncat(\"In the top 100 lists for EryP cells there are:\", nrow(eryp), \"unique RNA\\n\",\n    \"RNA with more than 1 intersection:\", nrow(filter(eryp, intersections > 1)),\"\\n\",\n   \"RNA with at least 7 intersections:\", nrow(filter(eryp, intersections > 6)) )\n \neryp[order(eryp$intersections, decreasing = T), ] %>%\n filter(intersections > 6)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-01-16T10:09:40.176060Z","iopub.execute_input":"2023-01-16T10:09:40.177739Z","iopub.status.idle":"2023-01-16T10:09:40.315047Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#eryp[order(eryp$intersections, decreasing = T), ] %>%\n# filter(intersections > 6) %>% select(RNA)","metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## MoP","metadata":{}},{"cell_type":"code","source":"mop <- filter(dat, cell_type == \"MoP\") %>%\n      mutate(set = paste0(donor,\"_\", day)) %>%\n      select(-c(sample_size, donor, day, cell_type, factor, group)) %>%\n      pivot_wider(names_from = \"set\", values_from = \"corr.coeff\") %>%\n      as.data.frame %>%\n      mutate(intersections = rowSums(!is.na(.[-(1:2)])) )\n\ncat(\"In the top 100 lists for MkP cells there are:\", nrow(mop), \"unique RNA\\n\",\n    \"RNA with more than 1 intersection:\", nrow(filter(mop, intersections > 1)),\"\\n\",\n   \"RNA with at least 7 intersections:\", nrow(filter(mop, intersections > 6)) )\n\nmop[order(mop$intersections, decreasing = T), ] %>%\n filter(intersections > 6)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-01-16T10:09:47.055990Z","iopub.execute_input":"2023-01-16T10:09:47.057684Z","iopub.status.idle":"2023-01-16T10:09:47.138901Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## MkP","metadata":{}},{"cell_type":"code","source":"mkp <- filter(dat, cell_type == \"MkP\") %>%\n      mutate(set = paste0(donor,\"_\", day)) %>%\n      select(-c(sample_size, donor, day, cell_type, factor, group)) %>%\n      pivot_wider(names_from = \"set\", values_from = \"corr.coeff\") %>%\n      as.data.frame %>%\n      mutate(intersections = rowSums(!is.na(.[-(1:2)])) )\n\ncat(\"In the top 100 lists for MkP cells there are:\", nrow(mkp), \"unique RNA\\n\",\n    \"RNA with more than 1 intersection:\", nrow(filter(mkp, intersections > 1)),\"\\n\",\n   \"RNA with at least 7 intersections:\", nrow(filter(mkp, intersections > 6)) )\n\nmkp[order(mkp$intersections, decreasing = T), ] %>%\n filter(intersections > 6)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-01-16T10:09:50.507118Z","iopub.execute_input":"2023-01-16T10:09:50.508676Z","iopub.status.idle":"2023-01-16T10:09:50.586803Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## HSC","metadata":{}},{"cell_type":"code","source":"hsc <- filter(dat, cell_type == \"HSC\") %>%\n      mutate(set = paste0(donor,\"_\", day)) %>%\n      select(-c(sample_size, donor, day, cell_type, factor, group)) %>%\n      pivot_wider(names_from = \"set\", values_from = \"corr.coeff\") %>%\n      as.data.frame %>%\n      mutate(intersections = rowSums(!is.na(.[-(1:2)])) )\n\ncat(\"In the top 100 lists for HSC cells there are:\", nrow(hsc), \"unique RNA\\n\",\n    \"RNA with more than 1 intersection:\", nrow(filter(hsc, intersections > 1)),\"\\n\",\n   \"RNA with at least 7 intersections:\", nrow(filter(hsc, intersections > 6)) )\n\nhsc[order(hsc$intersections, decreasing = T), ] %>%\n filter(intersections > 6)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-01-16T10:09:53.841401Z","iopub.execute_input":"2023-01-16T10:09:53.842927Z","iopub.status.idle":"2023-01-16T10:09:53.922071Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"XIST appearance","metadata":{}},{"cell_type":"code","source":"hsc[grepl(\"XIST\", hsc$RNA), ]","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-01-16T10:09:59.561820Z","iopub.execute_input":"2023-01-16T10:09:59.563431Z","iopub.status.idle":"2023-01-16T10:09:59.589046Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## MasP","metadata":{}},{"cell_type":"code","source":"masp <- filter(dat, cell_type == \"MasP\") %>%\n      mutate(set = paste0(donor,\"_\", day)) %>%\n      select(-c(sample_size, donor, day, cell_type, factor, group)) %>%\n      pivot_wider(names_from = \"set\", values_from = \"corr.coeff\") %>%\n      as.data.frame %>%\n      mutate(intersections = rowSums(!is.na(.[-(1:2)])) )\n\ncat(\"In the top 100 lists for MasP cells there are:\", nrow(masp), \"unique RNA\\n\",\n    \"RNA with more than 1 intersection:\", nrow(filter(masp, intersections > 1)),\"\\n\",\n   \"RNA with at least 7 intersections:\", nrow(filter(masp, intersections > 6)) )\n\nmasp[order(masp$intersections, decreasing = T), ] %>%\n filter(intersections > 6)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-01-16T10:10:02.334874Z","iopub.execute_input":"2023-01-16T10:10:02.336650Z","iopub.status.idle":"2023-01-16T10:10:02.411261Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## NeuP","metadata":{}},{"cell_type":"code","source":"sbgr <- filter(dat, cell_type == \"NeuP\") %>%\n      mutate(set = paste0(donor,\"_\", day)) %>%\n      select(-c(sample_size, donor, day, cell_type, factor, group)) %>%\n      pivot_wider(names_from = \"set\", values_from = \"corr.coeff\") %>%\n      as.data.frame %>%\n      mutate(intersections = rowSums(!is.na(.[-(1:2)])) )\n\ncat(\"In the top 100 lists for NeuP cells there are:\", nrow(sbgr), \"unique RNA\\n\",\n    \"RNA with more than 1 intersection:\", nrow(filter(sbgr, intersections > 1)),\"\\n\",\n   \"RNA with at least 7 intersections:\", nrow(filter(sbgr, intersections > 6)) )\n\nsbgr[order(sbgr$intersections, decreasing = T), ] %>%\n filter(intersections > 6)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-01-16T10:10:13.831602Z","iopub.execute_input":"2023-01-16T10:10:13.833132Z","iopub.status.idle":"2023-01-16T10:10:13.905738Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"XIST appearance","metadata":{}},{"cell_type":"code","source":"sbgr[grepl(\"XIST\", sbgr$RNA), ]","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-01-16T10:10:19.053252Z","iopub.execute_input":"2023-01-16T10:10:19.054801Z","iopub.status.idle":"2023-01-16T10:10:19.079172Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# 3. Different subsets for cell types","metadata":{}},{"cell_type":"markdown","source":"Well, we have found 49 RNA with at least 7 intersections in EryP. Here we look for intersections of these RNAs with RNA in the top 100 lists for other cell types (100 RNA per each donor and day, 900 RNA per cell type).","metadata":{}},{"cell_type":"code","source":"# EryP and MkP\ntmp <- eryp[eryp$intersections %in% c(7,8,9), ] %>%\n select(RNA) %>%\n filter(RNA %in% mkp$RNA)\ncat(\"Intersections between\", nrow(filter(eryp, intersections > 6)),\n    \"RNA with at least 7 intersections in EryP and all\", nrow(mkp), \"unique RNA in the top 100 lists for MkP:\")\ntmp\n\n# EryP and MasP\ntmp <- eryp[eryp$intersections %in% c(7,8,9), ] %>%\n select(RNA) %>%\n filter(RNA  %in% masp$RNA)\ncat(\"Intersections between\", nrow(filter(eryp, intersections > 6)),\n    \"RNA with at least 7 intersections in EryP and all\", nrow(masp), \"unique RNA in the top 100 lists for MasP:\")\ntmp\n\n# EryP and MoP\ntmp <- eryp[eryp$intersections %in% c(7,8,9), ] %>%\n select(RNA) %>%\n filter(RNA  %in% mop$RNA)\ncat(\"Intersections between\", nrow(filter(eryp, intersections > 6)),\n    \"RNA with at least 7 intersections in EryP and all\", nrow(mop), \"unique RNA in the top 100 lists for MoP:\")\ntmp\n\n# EryP and HSC\ntmp <- eryp[eryp$intersections %in% c(7,8,9), ] %>%\n select(RNA) %>%\n filter(RNA  %in% hsc$RNA)\ncat(\"Intersections between\", nrow(filter(eryp, intersections > 6)),\n    \"RNA with at least 7 intersections in EryP and all\", nrow(hsc), \"unique RNA in the top 100 lists for HSC:\")\ntmp","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-01-16T10:10:23.705568Z","iopub.execute_input":"2023-01-16T10:10:23.707460Z","iopub.status.idle":"2023-01-16T10:10:23.868503Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#XIST <- all_dat %>% filter(grepl(\"XIST\", RNA))\n#XIST[is.na(XIST[6:ncol(XIST)]),]$donor %>% unique\n#XIST[!is.na(XIST[6:ncol(XIST)]),]$donor %>% unique","metadata":{"execution":{"iopub.status.busy":"2023-01-04T09:41:29.061470Z","iopub.execute_input":"2023-01-04T09:41:29.063622Z","iopub.status.idle":"2023-01-04T09:41:30.501853Z"},"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"for (thing in ls()) { message(thing); print(object.size(get(thing)), units='auto') }\nrm(all_dat,dat, n,tmp,fact,cut_off,\n  prot,sbgr, hsc,masp,mkp,mop)\ngc()","metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"execution":{"iopub.status.busy":"2023-01-16T10:10:45.947643Z","iopub.execute_input":"2023-01-16T10:10:45.949299Z","iopub.status.idle":"2023-01-16T10:10:49.541524Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Expression levels by day","metadata":{}},{"cell_type":"code","source":"tic()\ndat_prot <- as.data.frame(mat_prot)\n#Merge with Metadata (day, donor, cell type, and technology)\nmetadata <- read.csv('../input/open-problems-multimodal/metadata.csv',row.names=1)\ndat_prot <- merge(metadata, dat_prot,  by=0, all.y = TRUE)\n\n# select one protein of interest\nprot_of_interest <- \"CD36\"\nmy_prot <- select(dat_prot, c(\"Row.names\",\"day\",\"donor\",\"cell_type\", all_of(prot_of_interest)))\nmy_prot$Protein <- prot_of_interest\nmy_prot$day <- factor(my_prot$day, levels = c('2','3','4'))\nnames(my_prot)[names(my_prot) == prot_of_interest] <- \"exp_lvl\"\n\n#Load RNA for the protein of interest\n#raw data\n#path <- \"/kaggle/input/sparse-raw-counts-data-open-problems-multimodal/citeseq/sp_train_cite_inputs_raw.rds\"\n#normalized data\npath <- \"/kaggle/input/sparse-measurement-data-open-problems-multimodal/sp_train_cite_inputs.rds\"\nmat_RNA <- readRDS(path)\n\n#dgCMatrix to dataframe\nmat_RNA <- as.matrix(mat_RNA)\n\n#Set names of the RNA of interest\nRNA <- eryp[eryp$intersections == 9, \"RNA\"]\nmy_RNA <- mat_RNA[, RNA] #select the needed columns\nrm(mat_RNA, dat_prot, mat_prot)\ngc()\n\n#Merge with Metadata (day, donor, cell type, and technology)\nmy_RNA <- merge(metadata, my_RNA,  by=0, all.y = TRUE)\nmy_RNA <- pivot_longer(my_RNA, cols = (6:ncol(my_RNA)), names_to = \"RNA\", values_to = \"exp_lvl\")\nmy_RNA$day <- factor(my_RNA$day, levels = c('2','3','4'))\ntoc()","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-01-03T18:26:35.858793Z","iopub.execute_input":"2023-01-03T18:26:35.861607Z","iopub.status.idle":"2023-01-03T18:27:49.489389Z"},"_kg_hide-output":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"cat(\"The protein data in the long format:\", nrow(my_prot), \"rows (\",length(prot_of_interest),\"protein x 70988 cells)\" )\nhead(my_prot)\n\ncat(\"The RNA data in the long format:\", nrow(my_RNA), \"rows (\",length(RNA),\"RNA x 70988 cells)\" )\nhead(my_RNA)","metadata":{"execution":{"iopub.status.busy":"2023-01-03T18:28:30.344428Z","iopub.execute_input":"2023-01-03T18:28:30.349767Z","iopub.status.idle":"2023-01-03T18:28:30.447183Z"},"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Expression of the protein of interest","metadata":{}},{"cell_type":"code","source":"fig(15,5)\nsbst <- my_prot\n\nggplot(data = sbst, aes(x = day, y = exp_lvl, colour = Protein))+\n geom_jitter() +\n ggtitle(paste(prot_of_interest, \"expression level by day\")) +\n theme_minimal(base_size = 18)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-01-03T18:32:01.578271Z","iopub.execute_input":"2023-01-03T18:32:01.579793Z","iopub.status.idle":"2023-01-03T18:32:05.873688Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig(22,10)\nsbst <- my_prot\n\nggplot(data = sbst, aes(x = day, y = exp_lvl, colour = Protein))+\n geom_jitter() +\n ggtitle(paste(prot_of_interest, \"expression level by day and cell type\")) +\n facet_wrap(~cell_type) +\n theme_minimal(base_size = 20)\n ","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-01-03T18:33:23.588530Z","iopub.execute_input":"2023-01-03T18:33:23.593543Z","iopub.status.idle":"2023-01-03T18:33:28.614756Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## TOP RNA with 9 intersections within EryP cell type: expression levels by day","metadata":{}},{"cell_type":"markdown","source":"## Plots","metadata":{}},{"cell_type":"code","source":"fig(22,10)\n\nfor (rna in RNA) {\n    sbst <- filter(my_RNA, RNA == rna)\n\n    plt <- \n     ggplot(data = sbst, aes(x = day, y = exp_lvl, colour = RNA))+\n     geom_jitter() +\n     ggtitle(paste(rna, \"expression level by day and cell type\")) +\n     facet_wrap(~cell_type) +\n     theme_minimal(base_size = 20)\n    print(plt)    \n}","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-01-03T18:33:43.122176Z","iopub.execute_input":"2023-01-03T18:33:43.123855Z","iopub.status.idle":"2023-01-03T18:35:28.837075Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Tables","metadata":{}},{"cell_type":"markdown","source":"### EryP","metadata":{}},{"cell_type":"code","source":"my_RNA %>%\n      filter(cell_type == \"EryP\") %>%\n      group_by(RNA, day) %>% #, cell_type\n      get_summary_stats(exp_lvl, type = \"common\") %>%\n      select(day, RNA, n, min,median, max, mean,sd) %>%\n      pivot_longer(cols = c(n, min,median, max, mean,sd), names_to = \"stats\", values_to = \"value\") %>%\n      pivot_wider(names_from = day, values_from = \"value\")","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-01-03T18:38:44.109625Z","iopub.execute_input":"2023-01-03T18:38:44.111449Z","iopub.status.idle":"2023-01-03T18:38:47.461553Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### MoP","metadata":{}},{"cell_type":"code","source":"my_RNA %>%\n      filter(cell_type == \"MoP\") %>%\n      group_by(RNA, day) %>% #, cell_type\n      get_summary_stats(exp_lvl, type = \"common\") %>%\n      select(day, RNA, n, min,median, max, mean,sd) %>%\n      pivot_longer(cols = c(n, min,median, max, mean,sd), names_to = \"stats\", values_to = \"value\") %>%\n      pivot_wider(names_from = day, values_from = \"value\")","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-01-03T18:38:54.181834Z","iopub.execute_input":"2023-01-03T18:38:54.183583Z","iopub.status.idle":"2023-01-03T18:38:58.210266Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### HSC","metadata":{}},{"cell_type":"code","source":"my_RNA %>%\n      filter(cell_type == \"MoP\") %>%\n      group_by(RNA, day) %>% #, cell_type\n      get_summary_stats(exp_lvl, type = \"common\") %>%\n      select(day, RNA, n, min,median, max, mean,sd) %>%\n      pivot_longer(cols = c(n, min,median, max, mean,sd), names_to = \"stats\", values_to = \"value\") %>%\n      pivot_wider(names_from = day, values_from = \"value\")","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-01-03T18:39:06.170483Z","iopub.execute_input":"2023-01-03T18:39:06.172292Z","iopub.status.idle":"2023-01-03T18:39:09.391771Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Number of dropouts for CD36 and its RNA ","metadata":{}},{"cell_type":"code","source":"dropouts_RNA <- filter(my_RNA, exp_lvl == 0 & RNA == \"ENSG00000135218_CD36\")\n\n# merge with protein data\ndropouts_prot <- my_prot %>%\n  rename(CD36_lvl = exp_lvl) %>%\n  select(Row.names, CD36_lvl) \ndropouts_RNA <- merge(dropouts_RNA, dropouts_prot, by = \"Row.names\", all.x = TRUE)\n\ncat(\"There are\", nrow(dropouts_RNA), \"rows with 0 values of CD36 RNA\")\ndropouts_RNA %>%\n group_by(day, cell_type) %>%\n summarise(n = n(), .groups = \"keep\") %>%\n pivot_wider(names_from = \"cell_type\", values_from = \"n\")","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-01-03T18:39:22.879967Z","iopub.execute_input":"2023-01-03T18:39:22.882638Z","iopub.status.idle":"2023-01-03T18:39:23.313565Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"ggplot(data = dropouts_RNA, aes(x = day, y = CD36_lvl, colour = CD36_lvl))+\n     geom_jitter() +\n     ggtitle(paste(gsub(\".*\\\\_\",\"\", unique(dropouts_RNA$RNA)),\n                   \"protein expression level by day and cell type. Only cells with zero RNA of that protein.\")) +\n     facet_wrap(~cell_type) +\n     theme_minimal(base_size = 18)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-01-03T18:39:27.011487Z","iopub.execute_input":"2023-01-03T18:39:27.013114Z","iopub.status.idle":"2023-01-03T18:39:30.996317Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#dat_prot %>%\n#  group_by(donor, day, cell_type) %>%\n#  summarise(n = n(), .groups = \"keep\") %>%\n#  pivot_wider(names_from = \"cell_type\", values_from = \"n\")","metadata":{"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#my_RNA %>%\n#  group_by(RNA, donor, day, cell_type) %>%\n#  summarise(n = n(), .groups = \"keep\") %>%\n#  pivot_wider(names_from = \"cell_type\", values_from = \"n\")","metadata":{"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]}]}