{"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":"What we do here:\n- violin/ridgeline plots of expression levels and correlation (normalized data) of several proteins, important for the prediction of cell type (according to [this notebok](http://www.kaggle.com/code/andreylalaley/cell-type-by-cite-important-features/notebook) or [ChatGPT](http://chat.openai.com/auth/login)), and their RNA","metadata":{}},{"cell_type":"markdown","source":"**Data:** RDS files with [sparse matrices of normalised counts](http://www.kaggle.com/datasets/stautxie/sparse-measurement-data-open-problems-multimodal) data 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.","metadata":{}},{"cell_type":"code","source":"# Libraries\nsuppressMessages(library(Matrix))\nsuppressMessages(library(openxlsx))\nsuppressMessages(library(dplyr)) \nsuppressMessages(library(tidyr)) #pivot_longer/wider\nsuppressMessages(library(ggplot2))\nsuppressMessages(library(tictoc))\nsuppressMessages(library(ggridges))\nsuppressMessages(library(ggExtra))\nsuppressMessages(library(ggpubr))\n\n# Function for figure size adjusment\nfig <- function(width, heigth) {\n    options(repr.plot.width = width, repr.plot.height = heigth) }","metadata":{"_uuid":"051d70d956493feee0c6d64651c6a088724dca2a","_execution_state":"idle","_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-12-19T20:45:03.004808Z","iopub.execute_input":"2022-12-19T20:45:03.042110Z","iopub.status.idle":"2022-12-19T20:45:03.076449Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Load data**","metadata":{}},{"cell_type":"markdown","source":"Targets: for the surface protein levels, each row corresponds to a cell (e.g. \"45006fe3e4c8\") and each column to a protein (e.g. \"CD86\").\n\nInputs: 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":"markdown","source":"Cell types:  \nMasP = Mast Cell Progenitor  \nMkP = Megakaryocyte Progenitor  \nNeuP = Neutrophil Progenitor  \nMoP = Monocyte Progenitor  \nEryP = Erythrocyte Progenitor  \nHSC = Hematoploetic Stem Cell  \nBP = B-Cell Progenitor","metadata":{}},{"cell_type":"code","source":"# Load Metadata (day, donor, cell type, and technology)\nmetadata <- read.csv('../input/open-problems-multimodal/metadata.csv',row.names=1)\n\nmetadata <- \n  metadata %>% \n  filter(technology == \"citeseq\" &\n         donor %in% c(\"13176\", \"31800\", \"32606\") &\n         day %in% c(\"2\",\"3\",\"4\")) %>%\n  mutate_all(as.character) %>%\n  mutate(\"Row.names\" = row.names(.))\n\n# Load RNA normalized data\n#path <- \"/kaggle/input/sparse-raw-counts-data-open-problems-multimodal/citeseq/sp_train_cite_inputs_raw.rds\"\npath <- \"/kaggle/input/sparse-measurement-data-open-problems-multimodal/sp_train_cite_inputs.rds\"\nmat_RNA <- readRDS(path)\n\n# dgCMatrix to matrix\nmat_RNA <- as.matrix(mat_RNA)\n\n# Load protein normalized data\n#path <- \"/kaggle/input/sparse-raw-counts-data-open-problems-multimodal/citeseq/sp_train_cite_targets_raw.rds\"\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)\ngc()","metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"execution":{"iopub.status.busy":"2022-12-19T20:45:20.839194Z","iopub.execute_input":"2022-12-19T20:45:20.841171Z","iopub.status.idle":"2022-12-19T20:46:21.413287Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Export RNA Ensemble IDs for all proteins in the data set**","metadata":{}},{"cell_type":"markdown","source":"To find needed RNA in data set we use Ensemble IDs collected previously in the [Research Project 01 Around Multimodal Single-Cell](http://www.kaggle.com/datasets/alexandervc/research-project-01-around-multimodal-singlecell).","metadata":{}},{"cell_type":"code","source":"path <- \"/kaggle/input/research-project-01-around-multimodal-singlecell/TotalSeq_B_Universal_Cocktail_v1_140_Antibodies_399904_Barcodes.xlsx\"\nensemblid <- read.xlsx(path, sheet = 1)\n#head(ensemblid) - a bit messy table, contain NAN at TCR\nno_ID <- ensemblid[is.na(ensemblid$Ensemble.ID), ]\nno_ID <- rbind(no_ID, filter(ensemblid, Ensemble.ID == \"NaN\"))\nensemblid <- filter(ensemblid, !Description %in% no_ID$Description)\nensemblid$Description <- gsub(\"\\\\s*\\\\([^\\\\)]+\\\\)\", \"\", ensemblid$Description)\nensemblid$Description <- gsub(\"anti-human \", \"\", ensemblid$Description)\nensemblid$Description <- gsub(\" Recombinant\", \"\", ensemblid$Description)\n\n#assign EnsemblIDs for each protein (names of some proteins differs from RNA names)\nfor (n in 1:ncol(mat_prot)) {\n    \n    id <- unique(ensemblid[grepl(paste0(colnames(mat_prot)[n],\"$\"), ensemblid$Description), \"Ensemble.ID\"])\n    colnames(mat_prot)[n] <- \n      paste0(colnames(mat_prot)[n],\"_\", id)\n    #print(paste(\"interation:\", n,\"found:\",\n    #ensemblid[grepl(paste0(colnames(mat_prot)[n],\"$\"),\n    #ensemblid$Description), \"Description\"],\n    #\"looking for:\", colnames(mat_prot)[n]))\n}\n\n#have to set manually several IDs (now ends with \"_\")\n#colnames(mat_prot)[grep('\\\\_$', colnames(mat_prot))]\ncolnames(mat_prot)[colnames(mat_prot) == \"FceRIa_\"] <- \"FceRIa_ENSG00000179639\"\ncolnames(mat_prot)[colnames(mat_prot) == \"HLA-A-B-C_\"] <- \"HLA-A-B-C_ENSG00000206503\"\ncolnames(mat_prot)[colnames(mat_prot) == \"integrinB7_\"] <- \"integrinB7_ENSG00000139626\"\n\n\ncat(\"Found EnsemblIDs for\", length(colnames(mat_prot)[!grepl('\\\\_$', colnames(mat_prot))]), \"of\", length(colnames(mat_prot)), \"proteins\")\ncat(\"\\nFor\", length(colnames(mat_prot)[grepl('\\\\_$', colnames(mat_prot))]),\n\"following proteins EnsemblIDs were not found.\")\ncolnames(mat_prot)[grep('\\\\_$', colnames(mat_prot))]","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-12-19T20:46:29.651280Z","iopub.execute_input":"2022-12-19T20:46:29.653041Z","iopub.status.idle":"2022-12-19T20:46:29.840859Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# a little work on renaming RNA and combining data\nprep_dat <- function(prot_name, RNA_name, excl_dropouts = TRUE) {\n    \n    my_prots <- mat_prot[, prot_name] %>%\n      as.data.frame\n    \n    if (ncol(my_prots) == 1) names(my_prots) <- prot_name\n    my_prots <- my_prots %>%\n      mutate(Row.names = rownames(.)) %>%\n      pivot_longer(cols = 1:(ncol(.)-1), values_to = \"Expression level (normalized)\") %>%\n      separate(col = \"name\", into = c(\"name\", \"id\"))\n    my_prots$group <- \"Protein\"\n    \n    my_prots$old_name <- NA\n\n    my_RNA <- mat_RNA[, RNA_name] %>%\n      as.data.frame \n    \n    if (ncol(my_RNA) == 1) names(my_RNA) <- RNA_name\n    my_RNA <- my_RNA %>%\n      mutate(Row.names = rownames(.)) %>%\n      pivot_longer(cols = 1:(ncol(.)-1), values_to = \"Expression level (normalized)\") %>%\n      separate(col = \"name\", into = c(\"id\", \"old_name\"))\n\n    my_RNA$name <- NA\n    for (id in unique(my_RNA$id)) {\n\n        try(my_RNA[my_RNA$id == id, ]$name <-\n          paste0(unique(my_prots[my_prots$id == id, ]$name), collapse = \"_\"))\n    }\n    my_RNA$group <- \"RNA\"\n    \n    if (excl_dropouts == TRUE) {\n        my_RNA <- filter(my_RNA, `Expression level (normalized)` != 0) #exclude dropouts \n        my_prots <- filter(my_prots, interaction(name, Row.names) %in% interaction(my_RNA$name, my_RNA$Row.names) )\n    }\n    all_dat <- rbind(my_RNA, my_prots) %>%\n      mutate_if(is.character, as.factor)\n\n    # add metadata\n    all_dat <- merge(all_dat, metadata, by = \"Row.names\") %>%\n      mutate_if(is.character, as.factor)\n    \n    # Add sample sizes to the data\n    sample_size <- all_dat %>%\n      group_by(group, name, cell_type) %>%\n      summarize(num=sum(!is.na(`Expression level (normalized)`)),.groups = \"keep\" )\n\n    all_dat <- all_dat %>%\n      left_join(sample_size, by = c(\"name\", \"group\", \"cell_type\")) %>%\n      mutate(sample_size_name = paste0(name, \"\\n\", \"n=\", num),\n             sample_size_ct = paste0(cell_type, \"\\n\", \"n=\", num))\n    return(all_dat)\n\n}","metadata":{"execution":{"iopub.status.busy":"2022-12-19T20:46:37.584342Z","iopub.execute_input":"2022-12-19T20:46:37.587038Z","iopub.status.idle":"2022-12-19T20:46:37.606722Z"},"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# EryP cell type","metadata":{}},{"cell_type":"markdown","source":"- Predictors according to [this notebok](http://www.kaggle.com/code/andreylalaley/cell-type-by-cite-important-features/notebook):  \n\"CD36\", \"CD71\",\"CD244\", \"CD38\", \"CD115\", \"CD81\", \"CD142\", \"CD88\", \"CD112\"\n- Predictors according to [ChatGPT](http://chat.openai.com/auth/login): CD71, CD36, CD235a (CD235a is not present in the data) !","metadata":{}},{"cell_type":"markdown","source":"**Select the protein/RNA of interest** ","metadata":{}},{"cell_type":"code","source":"# select the protein if interest\nprot_name <- c(\"CD36\", \"CD71\",\"CD244\", \"CD38\", \"CD115\", \"CD81\", \"CD142\", \"CD88\", \"CD112\")\n\n#find respective RNA\nRNA_name <- ensemblid[grepl(paste0(prot_name, collapse = \"|\"),\n                            ensemblid$Description,\n                            ignore.case = T),]$Ensemble.ID\n\nRNA_name <- colnames(mat_RNA)[which(grepl(paste0(RNA_name, collapse = \"|\"),\n                                     colnames(mat_RNA),\n                                     ignore.case = T))]\n\n#add EnsemblIDs to protein names \nprot_name <- colnames(mat_prot)[which(grepl(paste0(prot_name, collapse = \"|\"),\n                                     colnames(mat_prot),\n                                     ignore.case = T))]","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-12-19T21:12:36.704317Z","iopub.execute_input":"2022-12-19T21:12:36.705912Z","iopub.status.idle":"2022-12-19T21:12:36.809551Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Prepare the data for plotting**","metadata":{}},{"cell_type":"code","source":"all_dat <- prep_dat(prot_name, RNA_name, excl_dropouts = TRUE)\n\ncat(\"\\nRenamed RNA\")\nselect(all_dat, old_name, name, id) %>% filter(as.character(old_name) != as.character(name)) %>% na.omit %>% unique\n\ncat(\"\\nHead of the resulting data frame for plotting\")\nhead(all_dat[order(all_dat$Row.names, all_dat$name, all_dat$group, all_dat$cell_type),], n = 3)\ncat(\"\\nOveral number of cells:\", all_dat$Row.names %>% n_distinct)\ncat(\"\\nMean expression levels by cell type\")\nall_dat %>% group_by(name,group, cell_type) %>% summarise (mean_expression = mean(`Expression level (normalized)`),\n                                                          .groups = \"keep\") %>%\npivot_wider(names_from = \"cell_type\", values_from = c(\"mean_expression\"))","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-12-19T21:12:41.676039Z","iopub.execute_input":"2022-12-19T21:12:41.678783Z","iopub.status.idle":"2022-12-19T21:13:08.243432Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Plots to compare the expression of proteins and their RNA within cell type**  \n(Dropouts are excluded)","metadata":{}},{"cell_type":"code","source":"fig(20,12)\n\nall_dat[all_dat$cell_type %in% c(\"EryP\", \"HSC\"),] %>%\n  ggplot(aes(x = sample_size_name, y = `Expression level (normalized)`, fill = group)) +\n  geom_violin() +\n  xlab(\"\")+\n  geom_hline(yintercept=0, linetype=\"dashed\", color = \"black\") +\n  facet_wrap(~group+cell_type, scales = \"free\") +\n  geom_boxplot(width=0.1, color=\"black\", alpha=0.2,\n              outlier.shape = NA) +\n  theme_bw(base_size = 20) +\n  theme(legend.position = \"none\")","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-12-19T20:47:22.072878Z","iopub.execute_input":"2022-12-19T20:47:22.074350Z","iopub.status.idle":"2022-12-19T20:47:26.820550Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"all_dat[all_dat$cell_type %in% c(\"MkP\", \"MoP\"),] %>%\n  ggplot(aes(x = sample_size_name, y = `Expression level (normalized)`, fill = group)) +\n  geom_violin() +\n  xlab(\"\")+\n  geom_hline(yintercept=0, linetype=\"dashed\", color = \"black\") +\n  facet_wrap(~group+cell_type, scales = \"free\") +\n  geom_boxplot(width=0.1, color=\"black\", alpha=0.2,\n              outlier.shape = NA) +\n  theme_bw(base_size = 20) +\n  theme(legend.position = \"none\")","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-12-19T20:47:30.707581Z","iopub.execute_input":"2022-12-19T20:47:30.709126Z","iopub.status.idle":"2022-12-19T20:47:32.788407Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig(27,12)\nall_dat[all_dat$cell_type %in% c(\"BP\", \"NeuP\", \"MasP\"),] %>%\n  ggplot(aes(x=sample_size_name, y=`Expression level (normalized)`, fill=group)) +\n  geom_violin() +\n  xlab(\"\")+\n  geom_hline(yintercept=0, linetype=\"dashed\", color = \"black\") +\n  facet_wrap(~group+cell_type, scales = \"free\") +\n  geom_boxplot(width=0.1, color=\"black\", alpha=0.2,\n              outlier.shape = NA) +\n  theme_bw(base_size = 20) +\n  theme(legend.position = \"none\")","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-12-19T20:47:36.883135Z","iopub.execute_input":"2022-12-19T20:47:36.884681Z","iopub.status.idle":"2022-12-19T20:47:40.358383Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Plots to compare the expression of proteins and their RNA between cell types**  \n(Dropouts are excluded)","metadata":{}},{"cell_type":"code","source":"fig(27,12)\n\nall_dat[all_dat$name %in% c(\"CD36\", \"CD71\", \"CD88\"),] %>%\n  ggplot(aes(x = sample_size_ct, y = `Expression level (normalized)`, fill = group)) +\n  geom_violin() +\n  xlab(\"\")+\n  geom_hline(yintercept=0, linetype=\"dashed\", color = \"black\") +\n  facet_wrap(~group+name, scales = \"free\", ncol = 3) +\n  geom_boxplot(width=0.1, color=\"black\", alpha=0.2,\n              outlier.shape = NA) +\n  theme_bw(base_size = 20) +\n  theme(legend.position = \"none\")","metadata":{"execution":{"iopub.status.busy":"2022-12-19T20:47:51.420242Z","iopub.execute_input":"2022-12-19T20:47:51.421712Z","iopub.status.idle":"2022-12-19T20:47:54.892532Z"},"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig(27,12)\n\nall_dat[all_dat$name %in% c(\"CD112\", \"CD115\", \"CD142\"),] %>%\n  ggplot(aes(x = sample_size_ct, y = `Expression level (normalized)`, fill = group)) +\n  geom_violin() +\n  xlab(\"\")+\n  geom_hline(yintercept=0, linetype=\"dashed\", color = \"black\") +\n  facet_wrap(~group+name, scales = \"free\", ncol = 3) +\n  geom_boxplot(width=0.1, color=\"black\", alpha=0.2,\n              outlier.shape = NA) +\n  theme_bw(base_size = 20) +\n  theme(legend.position = \"none\")","metadata":{"execution":{"iopub.status.busy":"2022-12-19T20:48:04.278980Z","iopub.execute_input":"2022-12-19T20:48:04.280684Z","iopub.status.idle":"2022-12-19T20:48:07.079361Z"},"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig(27,12)\n\nall_dat[all_dat$name %in% c(\"CD244\", \"CD38\", \"CD81\"),] %>%\n  ggplot(aes(x = sample_size_ct, y = `Expression level (normalized)`, fill = group)) +\n  geom_violin() +\n  xlab(\"\")+\n  geom_hline(yintercept=0, linetype=\"dashed\", color = \"black\") +\n  facet_wrap(~group+name, scales = \"free\", ncol = 3) +\n  geom_boxplot(width=0.1, color=\"black\", alpha=0.2,\n              outlier.shape = NA) +\n  theme_bw(base_size = 20) +\n  theme(legend.position = \"none\")","metadata":{"execution":{"iopub.status.busy":"2022-12-19T20:48:11.070638Z","iopub.execute_input":"2022-12-19T20:48:11.072174Z","iopub.status.idle":"2022-12-19T20:48:15.058125Z"},"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Plots to compare the density of the expression of proteins and their RNA between cell types**  \n(Dropouts are excluded)","metadata":{}},{"cell_type":"code","source":"fig(27,12)\n\nall_dat[all_dat$name %in% c(\"CD36\", \"CD71\", \"CD88\", \"CD112\", \"CD115\"),] %>%\n  ggplot(aes(x = `Expression level (normalized)`, y = sample_size_ct, fill = cell_type)) +\n  geom_density_ridges() +\n  ylab(\"Cell type & sample size\")+\n  theme_ridges() +\n  #geom_hline(yintercept=0, linetype=\"dashed\", color = \"black\") +\n  facet_wrap(~group+name, scales = \"free\", ncol = 5) +\n  theme_classic(base_size = 20) +\n  theme(legend.position = \"none\")","metadata":{"execution":{"iopub.status.busy":"2022-12-19T20:48:19.220879Z","iopub.execute_input":"2022-12-19T20:48:19.222476Z","iopub.status.idle":"2022-12-19T20:48:22.921586Z"},"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig(22,12)\n\nall_dat[all_dat$name %in% c(\"CD142\", \"CD244\", \"CD38\", \"CD81\"),] %>%\n  ggplot(aes(x = `Expression level (normalized)`, y = sample_size_ct, fill = cell_type)) +\n  geom_density_ridges() +\n  ylab(\"Cell type & sample size\")+\n  theme_ridges() +\n  #geom_hline(yintercept=0, linetype=\"dashed\", color = \"black\") +\n  facet_wrap(~group+name, scales = \"free\", ncol = 4) +\n  theme_classic(base_size = 20) +\n  theme(legend.position = \"none\")","metadata":{"execution":{"iopub.status.busy":"2022-12-19T20:48:27.052713Z","iopub.execute_input":"2022-12-19T20:48:27.054159Z","iopub.status.idle":"2022-12-19T20:48:29.958959Z"},"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Plot protein expression vs its RNA**  \n(Dropouts are excluded)","metadata":{}},{"cell_type":"code","source":"all_dat <- pivot_wider(select(all_dat, -c(old_name)), names_from = \"group\", values_from = \"Expression level (normalized)\")\ncat(\"The head of the data frame for plotting\")\nhead(all_dat)","metadata":{"execution":{"iopub.status.busy":"2022-12-19T21:13:29.899896Z","iopub.execute_input":"2022-12-19T21:13:29.901594Z","iopub.status.idle":"2022-12-19T21:13:30.052232Z"},"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig(22,22)\n\nggscatter(all_dat[all_dat$cell_type == \"EryP\",], x = \"Protein\", y = \"RNA\",\n          shape = 20, alpha = 0.5, size = 3,\n          add = \"reg.line\", conf.int = TRUE, \n          cor.coef = TRUE, cor.coef.size = 8,\n          cor.coeff.args = list(method = \"spearman\", label.sep = \"\\n\"),\n          add.params = list(fill = \"lightgray\"),\n          facet.by = c(\"sample_size_name\"),\n          title = \"Spearman correlation, HSC cells, dropouts are excluded\",\n          ggtheme = theme_bw(base_size = 26) )","metadata":{"execution":{"iopub.status.busy":"2022-12-19T21:14:56.454782Z","iopub.execute_input":"2022-12-19T21:14:56.457398Z","iopub.status.idle":"2022-12-19T21:15:01.761311Z"},"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# HSC cell type","metadata":{}},{"cell_type":"markdown","source":"\n- Predictors according to [this notebok](http://www.kaggle.com/code/andreylalaley/cell-type-by-cite-important-features/notebook):\n  \"CD36\", \"CD47\", \"CD38\", \"CD41\", \"CD48\", \"CD62L\", \"CD49b\", \"CD63\", \"CD49d\", \"CD11a\", \"CD71\", \"CD123\", \"CD101\"\n- Predictors according to [ChatGPT](http://chat.openai.com/auth/login):  \nCD34, CD38, CD133, CD90, CD45 (CD34, CD133, CD90 - not present in the data)","metadata":{}},{"cell_type":"markdown","source":"**Select the protein/RNA of interest**  ","metadata":{}},{"cell_type":"code","source":"# select the protein if interest\nprot_name <- c(\"CD36\", \"CD47\", \"CD38\", \"CD41\", \"CD48\", \"CD62L\", \"CD49b\", \"CD63\", \"CD49d\", \"CD11a\", \"CD71\", \"CD123\", \"CD101\")\n\n#find respective RNA\nRNA_name <- ensemblid[grepl(paste0(prot_name, collapse = \"|\"),\n                            ensemblid$Description,\n                            ignore.case = T),]$Ensemble.ID\n\nRNA_name <- colnames(mat_RNA)[which(grepl(paste0(RNA_name, collapse = \"|\"),\n                                     colnames(mat_RNA),\n                                     ignore.case = T))]\n\n#add EnsemblIDs to protein names \nprot_name <- colnames(mat_prot)[which(grepl(paste0(prot_name, collapse = \"|\"),\n                                     colnames(mat_prot),\n                                     ignore.case = T))]","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-12-19T20:49:04.097099Z","iopub.execute_input":"2022-12-19T20:49:04.098642Z","iopub.status.idle":"2022-12-19T20:49:04.226419Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Prepare the data for plotting**","metadata":{}},{"cell_type":"code","source":"all_dat <- prep_dat(prot_name, RNA_name, excl_dropouts = TRUE)\n\ncat(\"\\nRenamed RNA\")\nselect(all_dat, old_name, name, id) %>% filter(as.character(old_name) != as.character(name)) %>% na.omit %>% unique\n\ncat(\"\\nHead of the resulting data frame for plotting\")\nhead(all_dat[order(all_dat$Row.names, all_dat$name, all_dat$group, all_dat$cell_type),], n = 3)\ncat(\"\\nOveral number of cells:\", all_dat$Row.names %>% n_distinct)\ncat(\"\\nMean expression levels by cell type\")\nall_dat %>% group_by(name,group, cell_type) %>% summarise (mean_expression = mean(`Expression level (normalized)`),\n                                                          .groups = \"keep\") %>%\npivot_wider(names_from = \"cell_type\", values_from = c(\"mean_expression\"))","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-12-19T20:49:10.764802Z","iopub.execute_input":"2022-12-19T20:49:10.766509Z","iopub.status.idle":"2022-12-19T20:49:53.519437Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Plots to compare the expression of proteins and their RNA within cell type**  \n(Dropouts are excluded)","metadata":{}},{"cell_type":"code","source":"fig(25,12)\n\nall_dat[all_dat$cell_type %in% c(\"EryP\", \"HSC\"),] %>%\n  ggplot(aes(x = sample_size_name, y = `Expression level (normalized)`, fill = group)) +\n  geom_violin() +\n  xlab(\"\")+\n  geom_hline(yintercept=0, linetype=\"dashed\", color = \"black\") +\n  facet_wrap(~group+cell_type, scales = \"free\") +\n  geom_boxplot(width=0.1, color=\"black\", alpha=0.2,\n              outlier.shape = NA) +\n  theme_bw(base_size = 20) +\n  theme(legend.position = \"none\")","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-12-19T20:50:20.425055Z","iopub.execute_input":"2022-12-19T20:50:20.426552Z","iopub.status.idle":"2022-12-19T20:50:27.559463Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig(25,12)\nall_dat[all_dat$cell_type %in% c(\"MkP\", \"MoP\"),] %>%\n  ggplot(aes(x = sample_size_name, y = `Expression level (normalized)`, fill = group)) +\n  geom_violin() +\n  xlab(\"\")+\n  geom_hline(yintercept=0, linetype=\"dashed\", color = \"black\") +\n  facet_wrap(~group+cell_type, scales = \"free\") +\n  geom_boxplot(width=0.1, color=\"black\", alpha=0.2,\n              outlier.shape = NA) +\n  theme_bw(base_size = 20) +\n  theme(legend.position = \"none\")","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-12-19T20:50:33.381028Z","iopub.execute_input":"2022-12-19T20:50:33.382516Z","iopub.status.idle":"2022-12-19T20:50:36.350490Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig(30,12)\nall_dat[all_dat$cell_type %in% c(\"BP\", \"NeuP\", \"MasP\"),] %>%\n  ggplot(aes(x=sample_size_name, y=`Expression level (normalized)`, fill=group)) +\n  geom_violin() +\n  xlab(\"\")+\n  geom_hline(yintercept=0, linetype=\"dashed\", color = \"black\") +\n  facet_wrap(~group+cell_type, scales = \"free\") +\n  geom_boxplot(width=0.1, color=\"black\", alpha=0.2,\n              outlier.shape = NA) +\n  theme_bw(base_size = 16) +\n  theme(legend.position = \"none\")","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-12-19T20:50:40.673738Z","iopub.execute_input":"2022-12-19T20:50:40.675289Z","iopub.status.idle":"2022-12-19T20:50:46.446293Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Plots to compare the expression of proteins and their RNA between cell types**  \n(Dropouts are excluded)","metadata":{}},{"cell_type":"code","source":"fig(27,12)\n\nall_dat[all_dat$name %in% c(\"CD11a\", \"CD36\", \"CD38\", \"CD41\"),] %>%\n  ggplot(aes(x = sample_size_ct, y = `Expression level (normalized)`, fill = group)) +\n  geom_violin() +\n  xlab(\"\")+\n  geom_hline(yintercept=0, linetype=\"dashed\", color = \"black\") +\n  facet_wrap(~group+name, scales = \"free\", ncol = 4) +\n  geom_boxplot(width=0.1, color=\"black\", alpha=0.2,\n              outlier.shape = NA) +\n  theme_bw(base_size = 20) +\n  theme(legend.position = \"none\")","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-12-19T20:50:52.222409Z","iopub.execute_input":"2022-12-19T20:50:52.224162Z","iopub.status.idle":"2022-12-19T20:50:56.444049Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig(27,12)\n\nall_dat[all_dat$name %in% c(\"CD47\", \"CD48\",\"CD49b\", \"CD49d\"),] %>%\n  ggplot(aes(x = sample_size_ct, y = `Expression level (normalized)`, fill = group)) +\n  geom_violin() +\n  xlab(\"\")+\n  geom_hline(yintercept=0, linetype=\"dashed\", color = \"black\") +\n  facet_wrap(~group+name, scales = \"free\", ncol = 4) +\n  geom_boxplot(width=0.1, color=\"black\", alpha=0.2,\n              outlier.shape = NA) +\n  theme_bw(base_size = 20) +\n  theme(legend.position = \"none\")","metadata":{"execution":{"iopub.status.busy":"2022-12-19T20:51:00.224212Z","iopub.execute_input":"2022-12-19T20:51:00.225682Z","iopub.status.idle":"2022-12-19T20:51:05.742946Z"},"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig(27,12) #\"CD101\" - small simple sizes (n = 2-115)\n\nall_dat[all_dat$name %in% c(\"CD62L\", \"CD63\", \"CD71\", \"CD123\"),] %>%\n  ggplot(aes(x = sample_size_ct, y = `Expression level (normalized)`, fill = group)) +\n  geom_violin() +\n  xlab(\"\")+\n  geom_hline(yintercept=0, linetype=\"dashed\", color = \"black\") +\n  facet_wrap(~group+name, scales = \"free\", ncol = 4) +\n  geom_boxplot(width=0.1, color=\"black\", alpha=0.2,\n              outlier.shape = NA) +\n  theme_bw(base_size = 20) +\n  theme(legend.position = \"none\")","metadata":{"execution":{"iopub.status.busy":"2022-12-18T22:12:38.316371Z","iopub.execute_input":"2022-12-18T22:12:38.318180Z","iopub.status.idle":"2022-12-18T22:12:44.549036Z"},"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Plots to compare the density of the expression of proteins and their RNA between cell types**  \n(Dropouts are excluded)","metadata":{}},{"cell_type":"code","source":"fig(27,12)\n\nall_dat[all_dat$name %in% c(\"CD11a\", \"CD36\", \"CD38\", \"CD41\",\"CD47\", \"CD48\"),] %>%\n  ggplot(aes(x = `Expression level (normalized)`, y = sample_size_ct, fill = cell_type)) +\n  geom_density_ridges() +\n  ylab(\"Cell type & sample size\")+\n  theme_ridges() +\n  #geom_hline(yintercept=0, linetype=\"dashed\", color = \"black\") +\n  facet_wrap(~group+name, scales = \"free\", ncol = 6) +\n  theme_classic(base_size = 20) +\n  theme(legend.position = \"none\")","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-12-18T22:16:25.277251Z","iopub.execute_input":"2022-12-18T22:16:25.282721Z","iopub.status.idle":"2022-12-18T22:16:31.006198Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig(27,12)\n\nall_dat[all_dat$name %in% c(\"CD49b\", \"CD49d\", \"CD62L\", \"CD63\", \"CD71\", \"CD123\"),] %>%\n  ggplot(aes(x = `Expression level (normalized)`, y = sample_size_ct, fill = cell_type)) +\n  geom_density_ridges() +\n  ylab(\"Cell type & sample size\")+\n  theme_ridges() +\n  #geom_hline(yintercept=0, linetype=\"dashed\", color = \"black\") +\n  facet_wrap(~group+name, scales = \"free\", ncol = 6) +\n  theme_classic(base_size = 20) +\n  theme(legend.position = \"none\")","metadata":{"execution":{"iopub.status.busy":"2022-12-18T22:17:41.089411Z","iopub.execute_input":"2022-12-18T22:17:41.091016Z","iopub.status.idle":"2022-12-18T22:17:46.470747Z"},"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Plot protein expression vs its RNA in HSC cells**  \n(Dropouts are excluded)","metadata":{}},{"cell_type":"code","source":"all_dat <- pivot_wider(select(all_dat, -c(old_name)), names_from = \"group\", values_from = \"Expression level (normalized)\")\ncat(\"The head of the data frame for plotting\")\nhead(all_dat)","metadata":{"execution":{"iopub.status.busy":"2022-12-19T20:53:29.096549Z","iopub.execute_input":"2022-12-19T20:53:29.098549Z","iopub.status.idle":"2022-12-19T20:53:29.385231Z"},"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig(27,27)\n\nggscatter(all_dat[all_dat$cell_type == \"HSC\",], x = \"Protein\", y = \"RNA\",\n          shape = 20, alpha = 0.5, size = 3,\n          add = \"reg.line\", conf.int = TRUE, \n          cor.coef = TRUE, cor.coef.size = 8,\n          cor.coeff.args = list(method = \"spearman\", label.sep = \"\\n\"),\n          add.params = list(fill = \"lightgray\"),\n          facet.by = c(\"sample_size_name\"),\n          title = \"Spearman correlation, HSC cells, dropouts are excluded\",\n          ggtheme = theme_bw(base_size = 26) )","metadata":{"execution":{"iopub.status.busy":"2022-12-19T21:10:35.855311Z","iopub.execute_input":"2022-12-19T21:10:35.858068Z","iopub.status.idle":"2022-12-19T21:10:54.439365Z"},"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Mouse/rat proteins","metadata":{}},{"cell_type":"markdown","source":"**Select the protein/RNA of interest**","metadata":{}},{"cell_type":"code","source":"# select the protein if interest\nprot_name <- c(\"Mouse-IgG1\", \"Mouse-IgG2a\", \"Mouse-IgG2b\", \"Rat-IgG2b\", \"Rat-IgG1\", \"Rat-IgG2a\")\n\n#add EnsemblIDs to protein names \nprot_name <- colnames(mat_prot)[which(grepl(paste0(prot_name, collapse = \"|\"),\n                                     colnames(mat_prot),\n                                     ignore.case = T))]","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-12-19T21:16:13.363544Z","iopub.execute_input":"2022-12-19T21:16:13.366149Z","iopub.status.idle":"2022-12-19T21:16:13.391235Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Prepare the data for plotting**","metadata":{}},{"cell_type":"code","source":"# a little work on renaming\nmy_prots <- mat_prot[, prot_name] %>%\n  as.data.frame %>%\n  mutate(Row.names = rownames(.)) %>%\n  pivot_longer(cols = 1:(ncol(.)-1), values_to = \"Expression level (normalized)\")# %>%\n  #separate(col = \"name\", into = c(\"name\", \"id\"))\nmy_prots$group <- \"Protein\"\n\nall_dat <- my_prots %>%\n  mutate_if(is.character, as.factor)\n\n# add metadata\nall_dat <- merge(all_dat, metadata, by = \"Row.names\") %>%\n  mutate_if(is.character, as.factor)\n\n# Add sample sizes to the data\nsample_size <- all_dat %>%\n  group_by(group, name, cell_type) %>%\n  summarize(num=sum(!is.na(`Expression level (normalized)`)),.groups = \"keep\" )\n\nall_dat <- all_dat %>%\n  left_join(sample_size, by = c(\"name\", \"group\", \"cell_type\")) %>%\n  mutate(sample_size = paste0(name, \"\\n\", \"n=\", num)) %>%\n  select(-num)\n\ncat(\"\\nHead of the resulting data frame for plotting\")\nhead(all_dat[order(all_dat$Row.names, all_dat$name, all_dat$group, all_dat$cell_type),], n = 3)\ncat(\"\\nOveral number of cells:\", all_dat$Row.names %>% n_distinct)\ncat(\"\\nMean expression levels by cell type\")\nall_dat %>% group_by(name, group, cell_type) %>% summarise (mean_expression = mean(`Expression level (normalized)`),\n                                                          .groups = \"keep\") %>%\npivot_wider(names_from = \"cell_type\", values_from = c(\"mean_expression\"))\n","metadata":{"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Plots**  \n(No RNA for them in the data set)","metadata":{}},{"cell_type":"code","source":"fig(20,7)\n\nall_dat[all_dat$cell_type %in% c(\"EryP\", \"HSC\"),] %>%\n  ggplot(aes(x = sample_size, y = `Expression level (normalized)`, fill = group)) +\n  geom_violin() +\n  xlab(\"\")+\n  geom_hline(yintercept=0, linetype=\"dashed\", color = \"black\") +\n  facet_wrap(~group+cell_type, scales = \"free\") +\n  geom_boxplot(width=0.1, color=\"black\", alpha=0.2,\n              outlier.shape = NA) +\n  theme_classic(base_size = 20) +\n  theme(legend.position = \"none\")","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-12-19T21:19:02.604253Z","iopub.execute_input":"2022-12-19T21:19:02.605938Z","iopub.status.idle":"2022-12-19T21:19:07.250780Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig(20,12)\n\nall_dat %>%\n  ggplot(aes(x = `Expression level (normalized)`, y = cell_type, fill = cell_type)) +\n  geom_density_ridges() +\n  #xlab(\"\")+\n  theme_ridges() +\n  #geom_hline(yintercept=0, linetype=\"dashed\", color = \"black\") +\n  facet_wrap(~group+name, scales = \"free\") +\n  theme_classic(base_size = 20) +\n  theme(legend.position = \"none\")","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-12-19T21:19:10.995641Z","iopub.execute_input":"2022-12-19T21:19:10.997168Z","iopub.status.idle":"2022-12-19T21:19:14.610137Z"},"trusted":true},"execution_count":null,"outputs":[]}]}