{"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":"# The protein of interest vs its RNA: correlation & expression with and without dropouts","metadata":{}},{"cell_type":"markdown","source":"What we do here:\n- Calculate Spearman correlations between the protein of interest and its RNA in all cells and by cell type, donor, day\n- Summarise their expression levels in all cells and by cell type, donor, day\n- Plot expression levels and correlations by day and cell type\n- Explore dropouts in the RNA data\n- Do the same exploratory analysis after dropouts removal\n- Plot expression levels and correlations by cell type for several other proteins","metadata":{}},{"cell_type":"markdown","source":"**Versions:**  \n8: prot_of_interest = CD36  \n9: prot_of_interest = CD115  \n10: prot_of_interest = CD71  \n11: prot_of_interest = CD88  \n12: prot_of_interest = CD44  \n13: prot_of_interest = CD45  \n14: prot_of_interest = CD45RA  ","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(gridExtra))\nsuppressMessages(library(grid))\nsuppressMessages(library(ggExtra))\nsuppressMessages(library(plotly)) #3D plot\n\nlibrary(ggpubr) # ggscatter\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,"_kg_hide-output":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Load 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\n**RNA (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\n**Metadata**: 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":"# 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)","metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Select the protein and RNA of interest**","metadata":{}},{"cell_type":"code","source":"prot_of_interest <- \"CD45RA\"","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# extract Ensembl IDs to identify the RNA of the protein of interest\nensembl <- read.csv('/kaggle/input/research-project-01-around-multimodal-singlecell/CD_to_EnsemblSymbol_correspondence_NIPS2022.csv',row.names = 1)\nensembl <- ensembl[ensembl$CD_name == prot_of_interest, ]\nensembl","metadata":{"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Find the RNA of interest and bind cell IDs\nRNA_name <- ensembl$ID.in.RNA.dataset\nmy_RNA <- mat_RNA[, RNA_name]\nmy_prot <- mat_prot[, prot_of_interest]","metadata":{"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# 1. All cells","metadata":{}},{"cell_type":"markdown","source":"**Calculate correlations**","metadata":{}},{"cell_type":"code","source":"#calculate correlations by cell_type\nout <- NULL\nfor (ct in unique(metadata$cell_type)) {\n    meta <- filter(metadata, cell_type == ct)\n\n    prot <- my_prot[names(my_prot) %in% meta$Row.names]\n    RNA <- my_RNA[names(my_RNA) %in% meta$Row.names]\n    corr <- cor(prot, RNA, method = \"spearman\") \n\n    tmp <- data.frame(\"group\" = \"All data\",\n                      \"cell_type\" = ct,\n                     \"corr.coeff\" = corr)\n    out <- rbind(out, tmp)            \n}\nout_tab <- out %>%\n  pivot_wider(names_from = \"cell_type\", values_from = \"corr.coeff\")\n#out_tab\n\n#calculate correlations by cell_type for each day and donor\nfactor = c(\"donor\", \"day\")\nout <- NULL\nfor (ct in unique(metadata$cell_type)) {\n    for(fa in (factor)) {\n        metadata$filt <- metadata[, fa]\n        for(gr in unique(metadata$filt)) {\n            meta <- filter(metadata, cell_type == ct & filt == gr)\n\n            prot <- my_prot[names(my_prot) %in% meta$Row.names]\n            RNA <- my_RNA[names(my_RNA) %in% meta$Row.names]\n            corr <- cor(prot, RNA, method = \"spearman\") \n\n            tmp <- data.frame(\"group\" = gr,\n                              \"cell_type\" = ct,\n                             \"corr.coeff\" = corr)\n            out <- rbind(out, tmp)\n        }            \n    }  \n}\n\nout <- pivot_wider(out, names_from = \"cell_type\", values_from = \"corr.coeff\")\nout_tab <- rbind(out_tab, out)\n#out_tab\n\n#calculate correlations for all data, each day, donor\nfactor = c(\"donor\", \"day\")\nout <- NULL\nfor(fa in (factor)) {\n    metadata$filt <- metadata[, fa]\n    for(gr in unique(metadata$filt)) {\n        meta <- filter(metadata, filt == gr)\n\n        prot <- my_prot[names(my_prot) %in% meta$Row.names]\n        RNA <- my_RNA[names(my_RNA) %in% meta$Row.names]\n        corr <- cor(prot, RNA, method = \"spearman\") \n\n        tmp <- data.frame(\"group\" = gr,\n                          \"cell_type\" = \"ALL\",\n                         \"corr.coeff\" = corr)\n        out <- rbind(out, tmp)    \n    }  \n}\nout <- pivot_wider(out, names_from = \"cell_type\", values_from = \"corr.coeff\")\nout <- rbind(data.frame(\n    \"group\" = \"All data\",\n    \"ALL\" = cor(my_prot, my_RNA, method = \"spearman\") ),\n             out)\n#out\n#combine all\nout_tab <- full_join(out, out_tab, by = \"group\")\ncat(\"Spearman correlations between\",prot_of_interest,\"and\", RNA_name,\n    \": in all data and  by cell type + day, donor\")\nout_tab","metadata":{"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Calculate expression levels**","metadata":{}},{"cell_type":"code","source":"#calculate expression levels by cell_type\nout_prot <- NULL\nout_RNA <- 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    RNA <- my_RNA[names(my_RNA) %in% meta$Row.names]\n    prot_expr <- mean(prot)\n    RNA_expr <- mean(RNA)\n    \n    tmp_prot <- data.frame(\"group\" = \"All data\",\n                           \"cell_type\" = ct,\n                           \"mean\" = prot_expr)\n    out_prot <- rbind(out_prot, tmp_prot)\n\n    tmp_RNA <- data.frame(\"group\" = \"All data\",\n                           \"cell_type\" = ct,\n                           \"mean\" = RNA_expr)\n    out_RNA <- rbind(out_RNA, tmp_RNA)\n}\n\nout_tab_prot <- out_prot %>%\n  pivot_wider(names_from = \"cell_type\", values_from = \"mean\")\nout_tab_RNA <- out_RNA %>%\n  pivot_wider(names_from = \"cell_type\", values_from = \"mean\")\n#out_tab_prot\n#out_tab_RNA\n\n#calculate expression levels by cell_type for each day and donor\nfactor = c(\"donor\", \"day\")\nout_prot <- NULL\nout_RNA <- NULL\nfor (ct in unique(metadata$cell_type)) {\n    for(fa in (factor)) {\n        metadata$filt <- metadata[, fa]\n        for(gr in unique(metadata$filt)) {\n            meta <- filter(metadata, cell_type == ct & filt == gr)\n\n            prot <- my_prot[names(my_prot) %in% meta$Row.names]\n            RNA <- my_RNA[names(my_RNA) %in% meta$Row.names]\n            prot_expr <- mean(prot)\n            RNA_expr <- mean(RNA)\n\n            tmp_prot <- data.frame(\"group\" = gr,\n                                   \"cell_type\" = ct,\n                                   \"mean\" = prot_expr)\n            out_prot <- rbind(out_prot, tmp_prot)\n\n            tmp_RNA <- data.frame(\"group\" = gr,\n                                   \"cell_type\" = ct,\n                                   \"mean\" = RNA_expr)\n            out_RNA <- rbind(out_RNA, tmp_RNA)\n        }            \n    }  \n}\n\nout_prot <- pivot_wider(out_prot, names_from = \"cell_type\", values_from = \"mean\")\nout_tab_prot <- rbind(out_tab_prot, out_prot)\n\nout_RNA <- pivot_wider(out_RNA, names_from = \"cell_type\", values_from = \"mean\")\nout_tab_RNA <- rbind(out_tab_RNA, out_RNA)\n#out_tab_prot\n#out_tab_RNA\n\n#calculate expression levels for all data, each day, donor\nfactor = c(\"donor\", \"day\")\nout_prot <- NULL\nout_RNA <- NULL\nfor(fa in (factor)) {\n    metadata$filt <- metadata[, fa]\n    for(gr in unique(metadata$filt)) {\n        meta <- filter(metadata, filt == gr)\n\n        prot <- my_prot[names(my_prot) %in% meta$Row.names]\n        RNA <- my_RNA[names(my_RNA) %in% meta$Row.names]        \n        prot_expr <- mean(prot)\n        RNA_expr <- mean(RNA)\n\n        tmp_prot <- data.frame(\"group\" = gr,\n                               \"cell_type\" = \"ALL\",\n                               \"mean\" = prot_expr)\n        out_prot <- rbind(out_prot, tmp_prot)\n\n        tmp_RNA <- data.frame(\"group\" = gr,\n                               \"cell_type\" = \"ALL\",\n                               \"mean\" = RNA_expr)\n        out_RNA <- rbind(out_RNA, tmp_RNA)    \n    }  \n}\nout_prot <- pivot_wider(out_prot, names_from = \"cell_type\", values_from = \"mean\")\nout_prot <- rbind(data.frame(\n    \"group\" = \"All data\",\n    \"ALL\" = mean(my_prot)), out_prot)\n#out_prot\n#combine all\nout_tab_prot <- full_join(out_prot, out_tab_prot, by = \"group\")\ncat(\"Mean expression levels of\",prot_of_interest, \": in all data and  by cell type + day, donor\")\nout_tab_prot\n\nout_RNA <- pivot_wider(out_RNA, names_from = \"cell_type\", values_from = \"mean\")\nout_RNA <- rbind(data.frame(\n    \"group\" = \"All data\",\n    \"ALL\" = mean(my_RNA)), out_RNA)\n#out_RNA\n#combine all\nout_tab_RNA <- full_join(out_RNA, out_tab_RNA, by = \"group\")\ncat(\"\\nMean expression levels of\",RNA_name, \": in all data and  by cell type + day, donor\")\nout_tab_RNA","metadata":{"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Plot the protein of interest expression levels by days**","metadata":{}},{"cell_type":"code","source":"options(scipen = 100) # suppress scientific notations\nmy_colors <- RColorBrewer::brewer.pal(9, 'Greens')[c(3,5,7)]\n\ndat <- data.frame(\"RNA\" = my_RNA, \"protein\" = my_prot)\ndat <- merge(dat,select(metadata, c(day,donor,cell_type)) , by = 0)\n\n# add sample sizes\n# the number of cells with protein and RNA data is the same, so calculate by any column\nsample_size <- dat %>%\n  group_by(cell_type) %>%\n  summarize(num_prot=sum(!is.na(protein)), .groups = \"keep\" )\n\ndat <- dat %>%\n  left_join(sample_size, by = c(\"cell_type\")) %>%\n  mutate(sample_size_ct = paste0(cell_type, \"\\n\", \"n=\", num_prot)) %>%\n  select(-num_prot)\ncat(\"The head of the data frame made for plotting.\")\nhead(dat)","metadata":{"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig(18,6)\nplt_prot <- ggplot(dat[dat$RNA != 0, ], aes(x = protein)) +\n    geom_histogram(binwidth = 2, color = \"black\", fill = my_colors[2]) +\n    xlab(paste(prot_of_interest, \"expression level (normalized)\")) +\n    theme_bw(base_size = 18) + theme(aspect.ratio = 1)\n\nplt_rna <- ggplot(dat[dat$RNA != 0, ], aes(x = RNA)) +\n    geom_histogram(binwidth = 0.1, color = \"black\", fill = my_colors[2]) +\n    xlab(paste(RNA_name, \"expression level (normalized)\")) +\n    theme_bw(base_size = 18) + theme(aspect.ratio = 1)\n\ngrid.arrange(plt_prot, plt_rna, ncol = 2,\n             top = textGrob(paste0(prot_of_interest,\n                                   \" protein and RNA expression levels, cells with RNA dropouts are removed (n = \", nrow(dat[dat$RNA != 0, ]), \")\" ),\n                            gp=gpar(fontsize=20))\n              )","metadata":{"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig(11,6)\nggplot(dat[dat$RNA == 0, ], aes(x = protein)) +\n    geom_histogram(binwidth = 2, color = \"black\", fill = my_colors[2]) +\n    xlab(paste(prot_of_interest, \"expression level (normalized)\")) +\n    theme_bw(base_size = 18) + theme(aspect.ratio = 1) +\n    ggtitle(paste0(prot_of_interest, \" expression level\\nin cells without its RNA (n = \", nrow(dat[dat$RNA == 0, ]), \")\" ))","metadata":{"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig(18,10)\nggplot(dat, aes(x = protein)) + #,fill=day, color= day #[dat$cell_type %in% c(\"EryP\", \"HSC\"),]\ngeom_histogram(binwidth=5, position=\"dodge\",colour=\"black\", fill = my_colors[2]) +\nxlab(paste(\"Expression level (normalized)\")) +\nfacet_wrap(~sample_size_ct, scales = \"fixed\", ncol = 4) +\ntheme_bw(base_size = 20) +\ntheme(aspect.ratio = 1) +\nggtitle(paste(prot_of_interest, \"expression levels\", \"by cell type (all cells)\")) +\nscale_fill_manual(values = my_colors)","metadata":{"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig(18,10)\nggplot(dat, aes(x = RNA)) + #[dat$cell_type %in% c(\"EryP\", \"HSC\"),]\ngeom_histogram(binwidth=0.4, position=\"dodge\",colour=\"black\", fill = my_colors[2]) +\nxlab(paste(RNA_name, \"Expression level (normalized)\")) +\nfacet_wrap(~sample_size_ct, scales = \"fixed\", ncol = 4) +\ntheme_bw(base_size = 20) +\ntheme(aspect.ratio = 1) +\nggtitle(paste(RNA_name, \"expression levels\", \"by cell type (all cells)\")) +\nscale_fill_manual(values = my_colors)","metadata":{"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**PLot Correlations**","metadata":{}},{"cell_type":"code","source":"fig(12,12)\n\nplt <- ggscatter(dat, x = \"protein\", y = \"RNA\",\n          shape = 20, alpha = 0.5, size = 5, color = my_colors[3],\n          add = \"reg.line\", conf.int = TRUE,\n          add.params = list(color = \"black\", fill = \"lightgray\"),\n          cor.coef = TRUE, cor.coef.size = 6,\n          cor.coeff.args = list(method = \"spearman\", label.sep = \"\\n\"),\n          xlab = prot_of_interest, ylab = RNA_name,\n          title = paste0(\"R - Spearman correlation between expression levels,\\nall data (n = \", n_distinct(dat$Row.names), \")\"),\n          ggtheme = theme_bw(base_size = 18)\n          ) + theme(aspect.ratio = 1)\n\nggMarginal(plt, type=\"histogram\", fill = my_colors[3])","metadata":{"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig(21,22)\nggscatter(dat, x = \"protein\", y = \"RNA\", #color = \"day\",\n          shape = 20, alpha = 0.5, size = 5, color = my_colors[3],\n          add = \"reg.line\", conf.int = TRUE,\n          cor.coef = TRUE, cor.coef.size = 6,\n          cor.coeff.args = list(method = \"spearman\", label.sep = \"\\n\"),\n          add.params = list(color = \"black\", fill = \"lightgray\"),\n          xlab = prot_of_interest, ylab = RNA_name,\n          title = paste(\"Spearman correlation by cell type\"),\n          ggtheme = theme_bw(base_size = 20),\n          facet.by = \"sample_size_ct\", fill = my_colors[3]) + theme(aspect.ratio = 1)","metadata":{"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# 2. After dropouts removal","metadata":{}},{"cell_type":"markdown","source":"**Number of dropouts**","metadata":{}},{"cell_type":"code","source":"dat <- data.frame(\"RNA\" = my_RNA, \"protein\" = my_prot)\ndat <- merge(dat,select(metadata, c(day,donor,cell_type)) , by = 0)\n\ndropouts <- filter(dat, RNA == 0)\n\ncat(\"There are\", nrow(dropouts), \"rows with CD36 RNA dropout.\n\\nNumber and percent of dropouts:\\n\")\n\ncat(\"\\nBy cell type:\")\ndropouts_c <- dropouts %>%\n group_by(cell_type) %>% #day,donor, \n summarise(n = n(), .groups = \"keep\")\n\nall <- dat %>%\n group_by(cell_type) %>% #day,donor, \n summarise(n = n(), .groups = \"keep\")\n\ndropouts_c$all <- all$n\n\ndropouts_c %>%\n mutate(prop = paste0(n,\"/\",all, \" (\", round(n/all*100, 2), \"%)\") ) %>%\n select(-c(n, all)) %>%\n pivot_wider(names_from = \"cell_type\", values_from = \"prop\")\n\ncat(\"\\nBy cell type, day:\")\ndropouts_cdd <- dropouts %>%\n group_by(day, cell_type) %>% \n summarise(n = n(), .groups = \"keep\")\n\nall <- dat %>%\n group_by(day, cell_type) %>% # \n summarise(n = n(), .groups = \"keep\")\n\ndropouts_cdd$all <- all$n\n\ndropouts_cdd %>%\n mutate(prop = paste0(n,\"/\",all, \" (\", round(n/all*100, 2), \"%)\") ) %>%\n select(-c(n, all)) %>%\n pivot_wider(names_from = \"cell_type\", values_from = \"prop\")\n\n\ndat <- filter(dat, RNA != 0)","metadata":{"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Calculate correlations**","metadata":{}},{"cell_type":"code","source":"#calculate correlations by cell_type\nout <- NULL\nfor (ct in unique(metadata$cell_type)) {\n    meta <- filter(metadata, cell_type == ct)\n\n    prot <- my_prot[names(my_prot) %in% meta$Row.names]\n    RNA <- my_RNA[names(my_RNA) %in% meta$Row.names]\n    \n    RNA <- RNA[RNA !=0]\n    prot <- prot[names(prot) %in% names(RNA)[RNA != 0] ]\n    \n    corr <- cor(prot, RNA, method = \"spearman\") \n    \n    tmp <- data.frame(\"group\" = \"All data\",\n                      \"cell_type\" = ct,\n                     \"corr.coeff\" = corr)\n    out <- rbind(out, tmp)            \n}\nout_tab <- out %>%\n  pivot_wider(names_from = \"cell_type\", values_from = \"corr.coeff\")\n#out_tab\n\n#calculate correlations by cell_type for each day and donor\nfactor = c(\"donor\", \"day\")\nout <- NULL\nfor (ct in unique(metadata$cell_type)) {\n    for(fa in (factor)) {\n        metadata$filt <- metadata[, fa]\n        for(gr in unique(metadata$filt)) {\n            meta <- filter(metadata, cell_type == ct & filt == gr)\n\n            prot <- my_prot[names(my_prot) %in% meta$Row.names]\n            RNA <- my_RNA[names(my_RNA) %in% meta$Row.names]\n            \n            RNA <- RNA[RNA !=0]\n            prot <- prot[names(prot) %in% names(RNA)[RNA != 0] ]    \n    \n            corr <- cor(prot, RNA, method = \"spearman\") \n\n            tmp <- data.frame(\"group\" = gr,\n                              \"cell_type\" = ct,\n                             \"corr.coeff\" = corr)\n            out <- rbind(out, tmp)\n        }            \n    }  \n}\n\nout <- pivot_wider(out, names_from = \"cell_type\", values_from = \"corr.coeff\")\nout_tab <- rbind(out_tab, out)\n#out_tab\n\n#calculate correlations for all data, each day, donor\nfactor = c(\"donor\", \"day\")\nout <- NULL\nfor(fa in (factor)) {\n    metadata$filt <- metadata[, fa]\n    for(gr in unique(metadata$filt)) {\n        meta <- filter(metadata, filt == gr)\n\n        prot <- my_prot[names(my_prot) %in% meta$Row.names]\n        RNA <- my_RNA[names(my_RNA) %in% meta$Row.names]\n        \n        RNA <- RNA[RNA !=0]\n        prot <- prot[names(prot) %in% names(RNA)[RNA != 0] ]\n        \n        corr <- cor(prot, RNA, method = \"spearman\") \n\n        tmp <- data.frame(\"group\" = gr,\n                          \"cell_type\" = \"ALL\",\n                         \"corr.coeff\" = corr)\n        out <- rbind(out, tmp)    \n    }  \n}\nout <- pivot_wider(out, names_from = \"cell_type\", values_from = \"corr.coeff\")\nout <- rbind(data.frame(\n    \"group\" = \"All data\",\n    \"ALL\" = cor(my_prot[names(my_prot) %in% names(my_RNA)[my_RNA != 0] ],\n                my_RNA[my_RNA !=0], method = \"spearman\") ),\n             out)\n#out\n#combine all\nout_tab <- full_join(out, out_tab, by = \"group\")\ncat(\"Correlations between\",prot_of_interest,\"and\", RNA_name,\n    \": in all data and  by cell type + day, donor\")\nout_tab","metadata":{"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Calculate expression levels**","metadata":{}},{"cell_type":"code","source":"#calculate expression levels by cell_type\nout_prot <- NULL\nout_RNA <- 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    RNA <- my_RNA[names(my_RNA) %in% meta$Row.names]\n    \n    RNA <- RNA[RNA !=0]\n    prot <- prot[names(prot) %in% names(RNA)[RNA != 0] ] \n    \n    prot_expr <- mean(prot)\n    RNA_expr <- mean(RNA)\n    \n    tmp_prot <- data.frame(\"group\" = \"All data\",\n                           \"cell_type\" = ct,\n                           \"mean\" = prot_expr)\n    out_prot <- rbind(out_prot, tmp_prot)\n\n    tmp_RNA <- data.frame(\"group\" = \"All data\",\n                           \"cell_type\" = ct,\n                           \"mean\" = RNA_expr)\n    out_RNA <- rbind(out_RNA, tmp_RNA)\n}\n\nout_tab_prot <- out_prot %>%\n  pivot_wider(names_from = \"cell_type\", values_from = \"mean\")\nout_tab_RNA <- out_RNA %>%\n  pivot_wider(names_from = \"cell_type\", values_from = \"mean\")\n#out_tab_prot\n#out_tab_RNA\n\n#calculate expression levels by cell_type for each day and donor\nfactor = c(\"donor\", \"day\")\nout_prot <- NULL\nout_RNA <- NULL\nfor (ct in unique(metadata$cell_type)) {\n    for(fa in (factor)) {\n        metadata$filt <- metadata[, fa]\n        for(gr in unique(metadata$filt)) {\n            meta <- filter(metadata, cell_type == ct & filt == gr)\n\n            prot <- my_prot[names(my_prot) %in% meta$Row.names]\n            RNA <- my_RNA[names(my_RNA) %in% meta$Row.names]\n            \n            RNA <- RNA[RNA !=0]\n            prot <- prot[names(prot) %in% names(RNA)[RNA != 0] ] \n            \n            prot_expr <- mean(prot)\n            RNA_expr <- mean(RNA)\n\n            tmp_prot <- data.frame(\"group\" = gr,\n                                   \"cell_type\" = ct,\n                                   \"mean\" = prot_expr)\n            out_prot <- rbind(out_prot, tmp_prot)\n\n            tmp_RNA <- data.frame(\"group\" = gr,\n                                   \"cell_type\" = ct,\n                                   \"mean\" = RNA_expr)\n            out_RNA <- rbind(out_RNA, tmp_RNA)\n        }            \n    }  \n}\n\nout_prot <- pivot_wider(out_prot, names_from = \"cell_type\", values_from = \"mean\")\nout_tab_prot <- rbind(out_tab_prot, out_prot)\n\nout_RNA <- pivot_wider(out_RNA, names_from = \"cell_type\", values_from = \"mean\")\nout_tab_RNA <- rbind(out_tab_RNA, out_RNA)\n#out_tab_prot\n#out_tab_RNA\n\n#calculate expression levels for all data, each day, donor\nfactor = c(\"donor\", \"day\")\nout_prot <- NULL\nout_RNA <- NULL\nfor(fa in (factor)) {\n    metadata$filt <- metadata[, fa]\n    for(gr in unique(metadata$filt)) {\n        meta <- filter(metadata, filt == gr)\n\n        prot <- my_prot[names(my_prot) %in% meta$Row.names]\n        RNA <- my_RNA[names(my_RNA) %in% meta$Row.names]        \n        \n        RNA <- RNA[RNA !=0]\n        prot <- prot[names(prot) %in% names(RNA)[RNA != 0] ] \n        \n        prot_expr <- mean(prot)\n        RNA_expr <- mean(RNA)\n\n        tmp_prot <- data.frame(\"group\" = gr,\n                               \"cell_type\" = \"ALL\",\n                               \"mean\" = prot_expr)\n        out_prot <- rbind(out_prot, tmp_prot)\n\n        tmp_RNA <- data.frame(\"group\" = gr,\n                               \"cell_type\" = \"ALL\",\n                               \"mean\" = RNA_expr)\n        out_RNA <- rbind(out_RNA, tmp_RNA)    \n    }  \n}\nout_prot <- pivot_wider(out_prot, names_from = \"cell_type\", values_from = \"mean\")\nout_prot <- rbind(data.frame(\n    \"group\" = \"All data\",\n    \"ALL\" = mean(my_prot[names(my_prot) %in% names(my_RNA)[my_RNA != 0] ] ) ),\n                  out_prot)\n#out_prot\n#combine all\nout_tab_prot <- full_join(out_prot, out_tab_prot, by = \"group\")\ncat(\"Mean expression levels of\",prot_of_interest, \": in all data and  by cell type + day, donor\")\nout_tab_prot\n\nout_RNA <- pivot_wider(out_RNA, names_from = \"cell_type\", values_from = \"mean\")\nout_RNA <- rbind(data.frame(\n    \"group\" = \"All data\",\n    \"ALL\" = mean(my_RNA[my_RNA !=0])), out_RNA)\n#out_RNA\n#combine all\nout_tab_RNA <- full_join(out_RNA, out_tab_RNA, by = \"group\")\ncat(\"\\nMean expression levels of\",RNA_name, \": in all data and  by cell type + day, donor\")\nout_tab_RNA","metadata":{"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**PLot Correlations**","metadata":{}},{"cell_type":"code","source":"# add sample sizes\n# the number of cells with protein and RNA data is the same, so calculate by any column\nsample_size <- dat %>%\n  group_by(cell_type) %>%\n  summarize(num_prot=sum(!is.na(protein)),.groups = \"keep\" )\n\ndat <- dat %>%\n  left_join(sample_size, by = c(\"cell_type\")) %>%\n  mutate(sample_size_ct = paste0(cell_type, \"\\n\", \"n=\", num_prot)) %>%\n  select(-num_prot)\ncat(\"The head of the data frame made for plotting.\")\nhead(dat)","metadata":{"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig(12,12)\n\nplt <- ggscatter(dat, x = \"protein\", y = \"RNA\",\n          shape = 20, alpha = 0.3, size = 5, color = my_colors[3],\n          add = \"reg.line\", conf.int = TRUE,\n          cor.coef = TRUE, cor.coef.size = 6,\n          cor.coeff.args = list(method = \"spearman\", label.sep = \"\\n\"),\n          add.params = list(color = \"black\", fill = \"lightgray\"),\n          xlab = prot_of_interest, ylab = RNA_name,\n          title = paste0(\"Spearman correlation, dropouts are removed (n = \", n_distinct(dat$Row.names), \")\"),\n          ggtheme = theme_bw(base_size = 20)) + theme(aspect.ratio = 1)\n\nggMarginal(plt, type=\"histogram\", fill = my_colors[3])","metadata":{"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig(21,22)\nggscatter(dat, x = \"protein\", y = \"RNA\", #color = \"day\",\n          shape = 20, alpha = 0.3, size = 5, color = my_colors[3],\n          add = \"reg.line\", conf.int = TRUE,\n          cor.coef = TRUE, cor.coef.size = 6,\n          cor.coeff.args = list(method = \"spearman\", label.sep = \"\\n\"),\n          add.params = list(color = \"black\", fill = \"lightgray\"),\n          xlab = prot_of_interest, ylab = RNA_name,\n          title = paste(\"Spearman correlation by cell type, dropouts are removed\"),\n          ggtheme = theme_bw(base_size = 20),\n         facet.by = \"sample_size_ct\") + theme(aspect.ratio = 1)","metadata":{"_kg_hide-input":true,"_kg_hide-output":false,"trusted":true},"execution_count":null,"outputs":[]}]}