{"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":"

Analysis of correlations of CD45 + isoforms and their RNA with all RNA (in progress).

","metadata":{}},{"cell_type":"markdown","source":"

What we do here:

\n
\n \n
","metadata":{}},{"cell_type":"markdown","source":"
\n\n

Data

\n\n
\n
    \n\n
  1. RDS files with sparse matrices of normalised counts data for Open Problems - Multimodal Single-Cell Integration (CITEseq2022). The dataset for this competition comprises single-cell multiomics data (n = 70988 cells) collected from mobilized peripheral CD34+ hematopoietic stem and progenitor cells (HSPCs) isolated from healthy human donors.\nIn total, the CITEseq2022 train dataset we use here, contains 140 surface proteins and 22050 RNA form 3 days, 3 donors, 7 cell types. The gene names are given by {EnsemblID}_{GeneName} where EnsemblID refers to the Ensembl Gene ID and GeneName to the gene name (e.g. \"ENSG00000159840_ZYX\").
    2. \n\n
  2. CSV files with log-normalised counts from the Multi-modal CITE-seq Prediction competition (CITEseq2023). The dataset consists of single-cell multiomics data: 25 proteins (Antibody-Derived Tags, ADT) and 639 highly variable genes (single-cell RNA-sequencing) in 4000 cells.
    3. \n \n
  3. CSV files with top300 lists of RNA that correlate positively or negatively with proteins and/or RNA of 112 protein-RNA pairs (lists were prepared in the Part 1 notebook)\n
\n
\n
","metadata":{}},{"cell_type":"code","source":"#Libraries\nsuppressPackageStartupMessages({\n \n library(Matrix)\n remotes::install_github(\"cysouw/qlcMatrix\")\n library(qlcMatrix) \n library(dplyr) \n library(tidyr)\n library(tibble) # add_column\n library(ggplot2)\n library(tictoc) #time measuring\n library(ggpubr) #ggscatter\n library(ggcorrplot) #ggcorrplot\n library(psych)\n library(grid)\n library(patchwork)\n library(corrplot)\n library(gridExtra) \n library(pheatmap)\n library(ggvenn) #venn diagrams with labels\n library(RColorBrewer)\n})\n\n# Function for figure size adjusment\nfig <- function(width, heigth) {\n options(repr.plot.width = width, repr.plot.height = heigth) }","metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"execution":{"iopub.status.busy":"2023-04-06T20:59:22.874815Z","iopub.execute_input":"2023-04-06T20:59:22.877196Z","iopub.status.idle":"2023-04-06T21:00:09.619471Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Load CITEseq data**","metadata":{}},{"cell_type":"code","source":"#CITEseq2022\n\n# 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\"\nsmat_RNA <- readRDS(path)\n\n# drop RNA with zero counts in all cells\nsmat_RNA <- smat_RNA[,unique(summary(smat_RNA)$j)]\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)\n\n\n#CITEseq2023\n\nmat_RNA23 <- t(as.matrix(read.csv(\"../input/machine-learning-challenge-2-prediction/training_set_rna.csv\", row.names = 1)))\nmat_prot23 <- t(as.matrix(read.csv(\"../input/machine-learning-challenge-2-prediction/training_set_adt.csv\", row.names = 1)))\n\ngc()","metadata":{"_uuid":"051d70d956493feee0c6d64651c6a088724dca2a","_execution_state":"idle","_kg_hide-output":true,"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-06T21:00:09.622474Z","iopub.execute_input":"2023-04-06T21:00:09.657660Z","iopub.status.idle":"2023-04-06T21:01:30.942043Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#describe data\npaste(\"**CITEseq2022 dataset**\")\ncat(\"\\nProteins:\", ncol(mat_prot),\n \"proteins in columns and\", nrow(mat_prot), \"cells in rows\")\ncat(\"\\nRNA:\", ncol(smat_RNA),\n \"RNA with non-zero expression in at least one cell in columns and\", nrow(smat_RNA), \"cells in rows\")\ncat(\"\\nMetadata:\")\nmetadata %>%\n group_by(technology) %>%\n summarise(donors = n_distinct(donor),\n days = n_distinct(day),\n cell_types = n_distinct(cell_type), \n .groups = \"keep\")\ncat(\"\\n\")\npaste(\"**CITEseq2023 dataset**\")\ncat(\"\\nProteins:\", ncol(mat_prot23),\n \"proteins in columns and\", nrow(mat_prot23), \"cells in rows\")\ncat(\"\\nRNA:\", ncol(mat_RNA23),\n \"RNA with non-zero expression in at least one cell in columns and\", nrow(mat_RNA23), \"cells in rows\")","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-06T21:01:30.944459Z","iopub.execute_input":"2023-04-06T21:01:30.945874Z","iopub.status.idle":"2023-04-06T21:01:31.054226Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Cell types in this dataset: \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":"markdown","source":"**Here are some functions for subsetting from the sparse and usual matrices and plotting**","metadata":{}},{"cell_type":"code","source":"#Function to subset needed columns, merge with metadata, and calculate sample sizes by cell_type\nmake_subset <- function(RNA_names, prot_names) {\n \n my_RNA <- smat_RNA[, RNA_names] %>% as.matrix\n \n #remove Ensembl Gene IDs for easy comparison between datasets\n colnames(my_RNA) <- gsub(\".*_\",\"\", colnames(my_RNA))\n \n my_RNA <- merge(metadata, my_RNA, by=0, all.y = TRUE)\n my_RNA$day <- factor(my_RNA$day, levels = c(\"2\",\"3\",\"4\"))\n \n mat_prot <- mat_prot[, prot_names] %>% as.data.frame %>% mutate(\"Row.names\" = rownames(.) )\n\n dat <- left_join(my_RNA, mat_prot, by = \"Row.names\")\n\n # add sample sizes per cell type\n sample_size_ct <- dat %>%\n group_by(cell_type) %>%\n summarize(ss_ct = n(),.groups = \"keep\" )\n \n # add sample sizes per day\n sample_size_da <- dat %>%\n group_by(day) %>%\n summarize(ss_da = n(),.groups = \"keep\" )\n \n # add sample sizes per donor\n sample_size_do <- dat %>%\n group_by(donor) %>%\n summarize(ss_do = n(),.groups = \"keep\" )\n\n dat <- dat %>%\n left_join(sample_size_ct, by = c(\"cell_type\")) %>%\n mutate(sample_size_ct = paste0(cell_type, \"\\n\", \"n=\", ss_ct)) %>%\n left_join(sample_size_da, by = c(\"day\")) %>%\n mutate(sample_size_day = paste0(day, \"\\n\", \"n=\", ss_da)) %>%\n left_join(sample_size_do, by = c(\"donor\")) %>%\n mutate(sample_size_donor = paste0(donor, \"\\n\", \"n=\", ss_do)) %>%\n select(-c(\"ss_ct\", \"ss_da\", \"ss_do\"))\n return(dat)\n}","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-06T21:01:31.057818Z","iopub.execute_input":"2023-04-06T21:01:31.059525Z","iopub.status.idle":"2023-04-06T21:01:31.076685Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# to draw a ggscatter of RNA vs protein use the function:\nmy_ggscatter <- function(set, prot, RNA, color = \"#033E8C\") {\n \n plt <- suppressMessages(ggscatter(set, x = prot, y = RNA, color = color,\n shape = 20, alpha = 0.2, 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\", color = \"black\",\n label.y.npc = \"bottom\",label.x.npc = \"right\", hjust = 1, vjust = -0.5),\n add.params = list(fill = \"lightgray\"),\n title = paste0(\"Spearman correlation,\\n(n = \", nrow(set), \")\"),\n ggtheme = theme_bw(base_size = 22) + theme(aspect.ratio = 1) ) )\n return(plt)\n} \n\n#to draw scatter plots for the two datasets we use the following functions\n# all data for both\ncomb_paris <- function(name1, name2) {\n \n sbst1 <- set22 %>% select(all_of(c(name1,name2)))\n names(sbst1) <- c(\"RNA\", \"Protein\")\n sbst1 <- sbst1 %>%\n mutate(data = \"CITEseq 2022\")\n\n sbst2<- as.data.frame(set23) %>% select(all_of(c(name1, name2)))\n names(sbst2) <- c(\"RNA\", \"Protein\")\n sbst2 <- sbst2 %>%\n mutate(data = \"CITEseq 2023\")\n\n plt_pair <- rbind(sbst1, sbst2)\n plt_pair$pair_name <- paste(name1, \"vs\", name2)\n \n plt <- ggscatter(plt_pair, x = \"Protein\", y = \"RNA\", color = \"#033E8C\",\n shape = 20, alpha = 0.2, 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\", color = \"black\",\n label.y.npc = \"bottom\",label.x.npc = \"right\", hjust = 1),\n add.params = list(fill = \"lightgray\"),\n facet.by = c(\"data\"), scales = \"free\", ncol = 4,\n title = paste0(\"Spearman correlation for \", unique(plt_pair$pair_name), \",\\n CITEseq22 (n = \", nrow(set22), \")\",\n \" vs CITEseq23 (n = \", nrow(set23), \")\"),\n ggtheme = theme_bw(base_size = 22) + theme(aspect.ratio = 1) ) \n \n g <- ggplot_gtable(ggplot_build(plt))\n strip_both <- which(grepl('strip-t', g$layout$name))\n fills <- c(\"#DFABA0\",\"#DFC79C\")\n k <- 1\n for (i in strip_both) {\n j <- which(grepl('rect', g$grobs[[i]]$grobs[[1]]$childrenOrder))\n g$grobs[[i]]$grobs[[1]]$children[[j]]$gp$fill <- fills[k]\n k <- k+1\n }\n grid.draw(g)\n}\n\n# by cell type for 2022 and 2023\nplot_2ds <- function(protein, RNA) {\n \n comp_pair <- \n rbind(set22 %>% \n select(all_of(c(RNA, protein)), \"sample_size_ct\"),\n set23 %>% as.data.frame %>% \n select(all_of(c(RNA, protein))) %>%\n mutate(sample_size_ct = paste0(\"CITEseq23 (n = \", nrow(set23), \")\") ) )\n comp_pair$sample_size_ct <- factor(comp_pair$sample_size_ct, levels = unique(comp_pair$sample_size_ct))\n\n fig(25,22)\n plt <- ggscatter(comp_pair, x = protein, y = RNA, color = \"#033E8C\",\n shape = 20, alpha = 0.2, 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\", color = \"black\",\n label.y.npc = \"bottom\",label.x.npc = \"right\", hjust = 1),\n add.params = list(fill = \"lightgray\"),\n facet.by = c(\"sample_size_ct\"), scales = \"free\",\n title = paste0(\"Spearman correlation, CITEseq22 data (n = \", nrow(set22),\n \") by cell type and CITEseq23 data (n = \", nrow(set23), \")\"),\n ggtheme = theme_bw(base_size = 26) + theme(aspect.ratio = 1) )\n\n #Change facet label text and background colour: \n #https://stackoverflow.com/questions/41631806/change-facet-label-text-and-background-colour\n g <- ggplot_gtable(ggplot_build(plt))\n strip_both <- which(grepl('strip-t', g$layout$name))\n strip_23 <- strip_both[strip_both %in% c(48)]\n strip_both <- strip_both[!strip_both %in% c(48, 49)]\n\n for (i in strip_both) {\n j <- which(grepl('rect', g$grobs[[i]]$grobs[[1]]$childrenOrder))\n g$grobs[[i]]$grobs[[1]]$children[[j]]$gp$fill <- \"#DFABA0\"\n }\n\n j <- which(grepl('rect', g$grobs[[strip_23]]$grobs[[1]]$childrenOrder))\n g$grobs[[strip_23]]$grobs[[1]]$children[[j]]$gp$fill <- \"#DFC79C\"\n\n grid.draw(g)\n}","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-06T21:01:31.081867Z","iopub.execute_input":"2023-04-06T21:01:31.084073Z","iopub.status.idle":"2023-04-06T21:01:31.105184Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#

Expression by cell type, day, donor

","metadata":{}},{"cell_type":"code","source":"RNA_names <- c(\"PTPRC\")\nRNA_names <- colnames(smat_RNA)[which(grepl(paste0(RNA_names, collapse = \"|\"),\n colnames(smat_RNA),\n ignore.case = T))] %>% as.character\nprot_names <- c(\"CD45\",\"CD45RO\",\"CD45RA\")\n\nset22 <- make_subset(RNA_names, prot_names)\n\ncat(\"The subset of the data with the RNAs and proteins of interest\")\nhead(set22)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-06T21:01:31.108902Z","iopub.execute_input":"2023-04-06T21:01:31.110793Z","iopub.status.idle":"2023-04-06T21:01:33.062022Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig(25,7)\n\ncat(\"Descriptive statistics of PTPRC, CD45, CD45RO, and CD45RA expression by cell type\")\n\nplt1 <- set22 %>% group_by(sample_size_ct) %>%\n get_summary_stats(c(\"PTPRC\"), show = c(\"min\", \"max\", \"median\"))\nplt2 <- set22 %>% group_by(sample_size_ct) %>%\n get_summary_stats(c(\"CD45\"), show = c(\"min\", \"max\", \"median\"))\nplt3 <- set22 %>% group_by(sample_size_ct) %>%\n get_summary_stats(c(\"CD45RO\"), show = c(\"min\", \"max\", \"median\"))\nplt4 <- set22 %>% group_by(sample_size_ct) %>%\n get_summary_stats(c(\"CD45RA\"), show = c(\"min\", \"max\", \"median\"))\n\ncbind(plt1, select(plt2, -c(\"sample_size_ct\")),\n select(plt3, -c(\"sample_size_ct\")),\n select(plt4, -c(\"sample_size_ct\")) )\n\nplt <- rbind(plt1, plt2, plt3, plt4)\n#plt\n\nggplot(plt, aes(x = sample_size_ct, y = median) ) +\n geom_bar(stat = 'identity', col=\"grey\", fill = \"#BF9039\") + \n theme_bw(base_size = 18) +\n xlab(\"\") +\n facet_wrap(~variable, scales = \"free\", ncol = 4) +\n ggtitle(paste0(\"Median expression levels of PTPRC, CD45RO, and CD45RA by cell type\"))","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-06T21:01:33.064980Z","iopub.execute_input":"2023-04-06T21:01:33.066517Z","iopub.status.idle":"2023-04-06T21:01:36.803691Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig(25,7)\ncat(\"Descriptive statistics of PTPRC, CD45, CD45RO, and CD45RA expression by day\")\n\nplt1 <- set22 %>% group_by(sample_size_day) %>%\n get_summary_stats(c(\"PTPRC\"), show = c(\"min\", \"max\", \"median\"))\nplt2 <- set22 %>% group_by(sample_size_day) %>%\n get_summary_stats(c(\"CD45\"), show = c(\"min\", \"max\", \"median\"))\nplt3 <- set22 %>% group_by(sample_size_day) %>%\n get_summary_stats(c(\"CD45RO\"), show = c(\"min\", \"max\", \"median\"))\nplt4 <- set22 %>% group_by(sample_size_day) %>%\n get_summary_stats(c(\"CD45RA\"), show = c(\"min\", \"max\", \"median\"))\n\ncbind(plt1, select(plt2, -c(\"sample_size_day\")),\n select(plt3, -c(\"sample_size_day\")),\n select(plt4, -c(\"sample_size_day\")) )\n\nplt <- rbind(plt1, plt2, plt3, plt4)\n#plt\n\nggplot(plt, aes(x = sample_size_day, y = median) ) +\n geom_bar(stat = 'identity', col=\"grey\", fill = \"#BF9039\") + \n theme_bw(base_size = 18) +\n xlab(\"\") +\n facet_wrap(~variable, scales = \"free\", ncol = 4) +\n ggtitle(paste0(\"Median expression levels of PTPRC, CD45RO, and CD45RA by day\")) +\n theme(aspect.ratio = 1)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-06T21:01:36.806207Z","iopub.execute_input":"2023-04-06T21:01:36.807624Z","iopub.status.idle":"2023-04-06T21:01:38.687910Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig(25,7)\ncat(\"Descriptive statistics of PTPRC, CD45RO, and CD45RA expression by donor\")\n\nplt1 <- set22 %>% group_by(sample_size_donor) %>%\n get_summary_stats(c(\"PTPRC\"), show = c(\"min\", \"max\", \"median\"))\nplt2 <- set22 %>% group_by(sample_size_donor) %>%\n get_summary_stats(c(\"CD45\"), show = c(\"min\", \"max\", \"median\"))\nplt3 <- set22 %>% group_by(sample_size_donor) %>%\n get_summary_stats(c(\"CD45RO\"), show = c(\"min\", \"max\", \"median\"))\nplt4 <- set22 %>% group_by(sample_size_donor) %>%\n get_summary_stats(c(\"CD45RA\"), show = c(\"min\", \"max\", \"median\"))\n\ncbind(plt1, select(plt2, -c(\"sample_size_donor\")),\n select(plt3, -c(\"sample_size_donor\")),\n select(plt4, -c(\"sample_size_donor\")) )\n\nplt <- rbind(plt1, plt2, plt3, plt4)\n#plt\n\nggplot(plt, aes(x = sample_size_donor, y = median) ) +\n geom_bar(stat = 'identity', col=\"grey\", fill = \"#BF9039\") + \n theme_bw(base_size = 18) +\n xlab(\"\") +\n facet_wrap(~variable, scales = \"free\", ncol = 4) +\n ggtitle(paste0(\"Median expression levels of PTPRC, CD45RO, and CD45RA by donor\")) +\n theme(aspect.ratio = 1)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-06T21:01:38.690596Z","iopub.execute_input":"2023-04-06T21:01:38.692146Z","iopub.status.idle":"2023-04-06T21:01:40.642773Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"
\n 💡 Note: \n \n
","metadata":{}},{"cell_type":"markdown","source":"#

Top RNAs correlated (|R| > 0.1) with CD45 isoforms (RO and RA) and/or their RNA (PTPRC) in all cells

","metadata":{}},{"cell_type":"markdown","source":"
\n 📌 Task: \n
\n In this notebook top 300 lists of RNA that correlate positively and negatively with proteins and/or RNA of protein-RNA pairs were matched and separated by:
\n- Common RNA - present in both top300 lists (for the protein and the RNA of the pair)
\n- Unique for protein - present only in the top300 list for the protein
\n- Unique for RNA - present only in the top300 list for the RNA

\nHere we take the top by the absolute correlation coefficient cutoff = 0.1 and visualise intersections between RNA lists\n
\n
\n\n\n\n","metadata":{}},{"cell_type":"markdown","source":"**Load top300 Spearman correlation lists with common and unique RNA for CD45 isoforms and their RNA**","metadata":{}},{"cell_type":"code","source":"pos_corr <- read.csv(paste0(\"/kaggle/input/mmscel-difference-in-prot-rna-correlations-p1/top300_pos_norib_mit_EIF_corr.csv\"))\npos_corr <- pos_corr %>% filter(grepl(\"CD45R\", Pair),\n corr.coef > 0.1)\ncat(\"Range of positive correlations, n - number of RNA in the list\")\npos_corr %>% pivot_wider(names_from = \"RNA_list\", id_cols = -RNA,\n values_from = \"corr.coef\",\n values_fn = function(x) paste(round(min(x), 3), \"-\", round(max(x), 3), \", n = \", length(x)) )\n \n \nneg_corr <- read.csv(paste0(\"/kaggle/input/mmscel-difference-in-prot-rna-correlations-p1/top300_neg_norib_mit_EIF_corr.csv\"))\nneg_corr <- neg_corr %>% filter(grepl(\"CD45R\", Pair),\n corr.coef < -0.1)\n#neg_corr <- arrange(neg_corr, corr.coef)\ncat(\"\\n\\nRange of negative correlations, n - number of RNA in the list\")\nneg_corr %>% pivot_wider(names_from = \"RNA_list\", id_cols = -RNA,\n values_from = \"corr.coef\",\n values_fn = function(x) paste(round(min(x), 3), \"-\", round(max(x), 3), \", n = \", length(x)) )","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-06T21:01:40.664928Z","iopub.execute_input":"2023-04-06T21:01:40.666383Z","iopub.status.idle":"2023-04-06T21:01:41.202998Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Venn diagrams to visualise intersections between RNA lists**","metadata":{}},{"cell_type":"code","source":"#Extract RNA names\npos_ROrna <- pos_corr %>% filter(grepl(\"CD45RO\", Pair) & RNA_list == \"Unique for RNA\") %>% select(RNA, corr.coef)\npos_RArna <- pos_corr %>% filter(grepl(\"CD45RA\", Pair) & RNA_list == \"Unique for RNA\") %>% select(RNA, corr.coef)\npos_RO <- pos_corr %>% filter(grepl(\"CD45RO\", Pair) &RNA_list == \"Unique for protein\") %>% select(RNA, corr.coef) \npos_RA <- pos_corr %>% filter(grepl(\"CD45RA\", Pair) &RNA_list == \"Unique for protein\") %>% select(RNA, corr.coef) \npos_com <- pos_corr %>% filter(grepl(\"CD45RA\", Pair) &RNA_list == \"Common for prot and RNA\") %>% select(RNA, corr.coef) \n\nneg_ROrna <- neg_corr %>% filter(grepl(\"CD45RO\", Pair) & RNA_list == \"Unique for RNA\") %>% select(RNA, corr.coef) \nneg_RArna <- neg_corr %>% filter(grepl(\"CD45RA\", Pair) & RNA_list == \"Unique for RNA\") %>% select(RNA, corr.coef)\nneg_RO <- neg_corr %>% filter(grepl(\"CD45RO\", Pair) &RNA_list == \"Unique for protein\") %>% select(RNA, corr.coef)\nneg_RA <- neg_corr %>% filter(grepl(\"CD45RA\", Pair) &RNA_list == \"Unique for protein\") %>% select(RNA, corr.coef)\nneg_com <- neg_corr %>% filter(grepl(\"CD45RA\", Pair) &RNA_list == \"Common for prot and RNA\") %>% select(RNA, corr.coef) ","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-06T21:01:41.205356Z","iopub.execute_input":"2023-04-06T21:01:41.206718Z","iopub.status.idle":"2023-04-06T21:01:41.470525Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig(25, 10)\nplt <- list(\"RO (+)\" = pos_RO$RNA, \n \"RA (+)\" = pos_RA$RNA,\n \"RO (-)\" = neg_RO$RNA, \n \"RA (-)\" = neg_RA$RNA)\n\np1 <- ggvenn(plt, show_elements = FALSE, label_sep = \"\\n\",\n text_size = 6,\n set_name_size = 8,\n fill_color = c(\"#BF9039\", \"#2D735F\", \"#BF9039\", \"#2D735F\")) +\nggtitle(paste(\"CD45RO vs CD45RA\")) +\ntheme(plot.title = element_text(size = 20, face = \"bold\"))\n\nplt <- list(\"RO (+)\" = pos_RO$RNA, \n \"RNA (+)\" = pos_ROrna$RNA,\n \"RO (-)\" = neg_RO$RNA, \n \"RNA (-)\" = neg_ROrna$RNA)\n\np2 <- ggvenn(plt, show_elements = FALSE, label_sep = \"\\n\",\n text_size = 6,\n set_name_size = 8,\n fill_color = c(\"#BF9039\", \"#2D735F\", \"#BF9039\", \"#2D735F\")) +\nggtitle(paste(\"CD45RO vs PTPRC\")) +\ntheme(plot.title = element_text(size = 20, face = \"bold\"))\n\nplt <- list(\"RA (+)\" = pos_RA$RNA, \n \"RNA (+)\" = pos_RArna$RNA,\n \"RA (-)\" = neg_RA$RNA, \n \"RNA (-)\" = neg_RArna$RNA)\n\np3 <- ggvenn(plt, show_elements = FALSE, label_sep = \"\\n\",\n text_size = 6,\n set_name_size = 8,\n fill_color = c(\"#BF9039\", \"#2D735F\", \"#BF9039\", \"#2D735F\")) +\nggtitle(paste(\"CD45RA vs PTPRC\")) +\ntheme(plot.title = element_text(size = 20, face = \"bold\"))\n\np1 + p2 + p3 + plot_annotation(\n title = \"Intersections between top RNA correlatied (|R| > 0.1) with CD45 isoforms (RO and RA) or their RNA (PTPRC)\",\n subtitle = \"+ positive correlation\\n- negative correlation\\n\") & \ntheme(text = element_text(size = 22) )","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-06T21:01:41.472960Z","iopub.execute_input":"2023-04-06T21:01:41.474312Z","iopub.status.idle":"2023-04-06T21:01:42.743028Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"
\n 💡 Note: \n \n
","metadata":{}},{"cell_type":"markdown","source":"#

All correlations

","metadata":{}},{"cell_type":"markdown","source":"#

- All cells

","metadata":{}},{"cell_type":"markdown","source":"
\n 📌 Task: \n
\n Here we look at how the RNAs that correlate best with the protein correlate with its RNA and vice versa.\n
\n
","metadata":{}},{"cell_type":"markdown","source":"**Functions for memory efficient Spearman correlation on sparse matrices**","metadata":{}},{"cell_type":"code","source":"# Memory efficient spearman correlation on sparse matrices\n\nSparsifiedRanks2 <- function(X) {\n if (class(X)[1] != \"dgCMatrix\") {\n #X <- as(object = X, Class = \"dgCMatrix\")\n X <- Matrix(X, sparse = TRUE)\n }\n non_zeros_per_col <- diff(x = X@p)\n n_zeros_per_col <- nrow(x = X) - non_zeros_per_col\n offsets <- (n_zeros_per_col - 1) / 2\n #in case there is no zeros (protein matrix)\n if (any(offsets < 0)) {\n offsets[offsets < 0] <- 0\n }\n \n x <- X@x\n ## split entries to columns\n col_lst <- split(x = x, f = rep.int(1:ncol(X), non_zeros_per_col))\n ## calculate sparsified ranks and do shifting\n sparsified_ranks <- unlist(x = lapply(X = seq_along(col_lst), \n FUN = function(i) rank(x = col_lst[[i]]) + offsets[i]))\n ## Create template rank matrix\n X.ranks <- X\n X.ranks@x <- sparsified_ranks\n return(X.ranks)\n}\n\nSparseSpearmanCor2 <- function(X, Y = NULL, cov = FALSE) {\n\n # Get sparsified ranks\n rankX <- SparsifiedRanks2(X)\n if (is.null(Y)){\n # Calculate pearson correlation on rank matrices\n return (qlcMatrix::corSparse(X=rankX, cov=cov))\n }\n rankY <- SparsifiedRanks2(Y)\n return(qlcMatrix::corSparse( X = rankX, Y = rankY, cov = cov))\n}","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-06T21:01:42.745399Z","iopub.execute_input":"2023-04-06T21:01:42.746847Z","iopub.status.idle":"2023-04-06T21:01:42.762265Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Calculate correlations between CD45RA and its RNA vs all RNA in all cells**","metadata":{}},{"cell_type":"code","source":"# Calculate RNA-RNA correlations in all cells\nmy_genes <- cbind(mat_prot[, \"CD45RA\"],\n as.matrix(smat_RNA[, \"ENSG00000081237_PTPRC\"]))\ngene_names <- c(\"CD45RA\", \"PTPRC\")\n\nRNA_cor_all <- SparseSpearmanCor2(my_genes, smat_RNA)\ncolnames(RNA_cor_all) <- colnames(smat_RNA) \nRNA_cor_all <- cbind(Name = gene_names,\n group = \"all\", RNA_cor_all,\n sample_size = nrow(my_genes))","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-06T21:01:42.764591Z","iopub.execute_input":"2023-04-06T21:01:42.765973Z","iopub.status.idle":"2023-04-06T21:04:13.163747Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#reshape to long\nall_corr_long <- RNA_cor_all %>% as.data.frame %>%\n pivot_longer(cols = 3:(ncol(.)-1), names_to = \"RNA\", values_to = \"corr.coef\") %>%\n pivot_wider(names_from = \"Name\", values_from = \"corr.coef\")\nall_corr_long[, 4:5] <- apply(all_corr_long[, 4:5], 2, as.numeric)\n\n#remove correlations with ENSG00000081237_PTPRC\nall_corr_long <- all_corr_long[all_corr_long$RNA != \"ENSG00000081237_PTPRC\", ]\n\nall_corr_long$PTPRC_rank <- \n ifelse(all_corr_long$PTPRC > 0, rank(all_corr_long$PTPRC*-1), rank(all_corr_long$PTPRC))\nall_corr_long$CD45RA_rank <- \n ifelse(all_corr_long$CD45RA > 0, rank(all_corr_long$CD45RA*-1), rank(all_corr_long$CD45RA))\n\n# all_corr_long$PTPRC_rank <- \n# ifelse(all_corr_long$PTPRC_rank > quantile(all_corr_long$CD45RA_rank, 0.8))\n\nhead(all_corr_long)\n#options(scipen = 1000)\n#cor.test(as.matrix(smat_RNA[, \"ENSG00000162378_ZYG11B\"]), mat_prot[, \"CD45RA\"], method = \"spearman\")","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-06T21:04:13.166208Z","iopub.execute_input":"2023-04-06T21:04:13.167536Z","iopub.status.idle":"2023-04-06T21:04:13.519085Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Descriptive statistics for correlation coefficients (CD45RA protein and RNA vs all RNA) and their ranks**","metadata":{}},{"cell_type":"code","source":"fig(25, 6)\n#Correlation coefficients in all cells\n\ntbl1 <- all_corr_long %>% \n get_summary_stats(c(\"CD45RA\", \"PTPRC\"),\n show = c(\"n\", \"min\", \"max\", \"median\", \"mean\", \"sd\"))\ntbl1 <- tableGrob(tbl1, theme = ttheme_minimal(base_size = 20), rows = NULL)\np1 <- ggplot(all_corr_long, aes(PTPRC)) +\n geom_histogram(bins = 40, fill=\"#BF5841\", col=\"grey\") + \n theme_bw(base_size = 26)\n\np2 <- ggplot(all_corr_long, aes(CD45RA)) +\n geom_histogram(bins = 40, fill=\"#BF5841\", col=\"grey\") + \n theme_bw(base_size = 20)\n\ngrid.arrange(\n p1, p2, tbl1, ncol=3, \n top=textGrob(\"Descriptive statistics for correlation coefficients - vs all RNA, in all cells\", gp=gpar(fontsize = 24), x = 0.25) )","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-06T21:04:13.521347Z","iopub.execute_input":"2023-04-06T21:04:13.522637Z","iopub.status.idle":"2023-04-06T21:04:14.296722Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Subset all RNA positively and negatively correlated (|R| > 0.1) with the CD45RA protein and/or its RNA**","metadata":{}},{"cell_type":"code","source":"RA_pos_corr <- pos_corr %>% filter(grepl(\"CD45RA\", Pair)) %>% select(\"RNA\", \"RNA_list\")\nset22_pos <- all_corr_long %>% filter(RNA %in% RA_pos_corr$RNA)\nset22_pos <- set22_pos %>% left_join(RA_pos_corr, by = \"RNA\")\n#set22_pos$RNA <- gsub(\".*_\", \"\", set22_pos$RNA)\ncat(\"The data frame with \", nrow(set22_pos), \"RNA, positively correlated with the protein and/or its RNA\")\nhead(set22_pos)\n\nRA_neg_corr <- neg_corr %>% filter(grepl(\"CD45RA\", Pair)) %>% select(\"RNA\", \"RNA_list\")\nset22_neg <- all_corr_long %>% filter(RNA %in% RA_neg_corr$RNA)\nset22_neg <- set22_neg %>% left_join(RA_neg_corr, by = \"RNA\")\n#set22_neg$RNA <- gsub(\".*_\", \"\", set22_neg$RNA)\ncat(\"\\nThe data frame with \", nrow(set22_neg), \"RNA, negatively correlated with the protein and/or its RNA\")\nhead(set22_neg)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-06T21:04:14.299156Z","iopub.execute_input":"2023-04-06T21:04:14.300535Z","iopub.status.idle":"2023-04-06T21:04:14.412941Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig(25, 12)\n\n#--------Positive------------\n\n#Correlation coefficients\ntbl1 <- set22_pos %>% \n get_summary_stats(c(\"CD45RA\", \"PTPRC\"),\n show = c(\"n\", \"min\", \"max\", \"median\", \"mean\", \"sd\"))\ntbl1 <- tableGrob(tbl1, theme = ttheme_minimal(base_size = 20), rows = NULL)\np1 <- ggplot(set22_pos, aes(PTPRC)) +\n geom_histogram(bins = 40, fill=\"#BF5841\", col=\"grey\") + \n theme_bw(base_size = 20) +\n ggtitle(\"Correlation coefficients\")\n\np2 <- ggplot(set22_pos, aes(CD45RA)) +\n geom_histogram(bins = 40, fill=\"#BF5841\", col=\"grey\") + \n theme_bw(base_size = 20) +\n ggtitle(\"\")\n\n#Ranks of correlation coefficients\ntbl2 <- set22_pos %>% \n get_summary_stats(c(\"CD45RA_rank\", \"PTPRC_rank\"),\n show = c(\"n\", \"min\", \"max\", \"median\", \"mean\", \"sd\"))\ntbl2 <- tableGrob(tbl2, theme = ttheme_minimal(base_size = 20), rows = NULL)\np3 <- ggplot(set22_pos, aes(PTPRC_rank)) +\n geom_histogram(bins = 40, fill=\"#BF5841\", col=\"grey\") + \n theme_bw(base_size = 20) +\n ggtitle(\"Correlation coefficient ranks\")\n\np4 <- ggplot(set22_pos, aes(CD45RA_rank)) +\n geom_histogram(bins = 40, fill=\"#BF5841\", col=\"grey\") + \n theme_bw(base_size = 20) +\n ggtitle(\"\")\n\ngrid.arrange(\n p1, p2, tbl1, \n p3, p4, tbl2, ncol=3,\n top=textGrob(paste0(nrow(set22_pos),\n \" RNA from the top 300 positively correlated with the CD45RA protein and/or its RNA\"),\n gp=gpar(fontsize = 26), x = 0.35) )\n\n#--------Negative------------\n\n#Correlation coefficients\ntbl1 <- set22_neg %>% \n get_summary_stats(c(\"CD45RA\", \"PTPRC\"),\n show = c(\"n\", \"min\", \"max\", \"median\", \"mean\", \"sd\"))\ntbl1 <- tableGrob(tbl1, theme = ttheme_minimal(base_size = 20), rows = NULL)\np1 <- ggplot(set22_neg, aes(PTPRC)) +\n geom_histogram(bins = 40, fill=\"#BF5841\", col=\"grey\") + \n theme_bw(base_size = 20) +\n ggtitle(\"Correlation coefficients\")\n\np2 <- ggplot(set22_neg, aes(CD45RA)) +\n geom_histogram(bins = 40, fill=\"#BF5841\", col=\"grey\") + \n theme_bw(base_size = 20) +\n ggtitle(\"\")\n\n#Ranks of correlation coefficients\ntbl2 <- set22_neg %>% \n get_summary_stats(c(\"CD45RA_rank\", \"PTPRC_rank\"),\n show = c(\"n\", \"min\", \"max\", \"median\", \"mean\", \"sd\"))\ntbl2 <- tableGrob(tbl2, theme = ttheme_minimal(base_size = 20), rows = NULL)\np3 <- ggplot(set22_neg, aes(PTPRC_rank)) +\n geom_histogram(bins = 40, fill=\"#BF5841\", col=\"grey\") + \n theme_bw(base_size = 20) +\n ggtitle(\"Correlation coefficient ranks\")\n\np4 <- ggplot(set22_neg, aes(CD45RA_rank)) +\n geom_histogram(bins = 40, fill=\"#BF5841\", col=\"grey\") + \n theme_bw(base_size = 20) +\n ggtitle(\"\")\n\ngrid.arrange(\n p1, p2, tbl1, \n p3, p4, tbl2, ncol=3,\n top=textGrob(paste0(nrow(set22_neg),\n \" RNA from the top 300 negatively correlated with the CD45RA protein and/or its RNA\"),\n gp=gpar(fontsize = 26), x = 0.35) )\n","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-06T21:04:14.415328Z","iopub.execute_input":"2023-04-06T21:04:14.416633Z","iopub.status.idle":"2023-04-06T21:04:17.451254Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig(25,8)\n\n#--------Positive---------\np1 <- ggplot(set22_pos, aes(x = `CD45RA`, y = `PTPRC`, color = `RNA_list`)) +\n geom_point() +\n theme_classic(base_size = 22) +\n geom_hline(yintercept = 0, linetype = \"dashed\", color = \"black\") +\n geom_vline(xintercept = 0,linetype = \"dashed\", color = \"black\") +\n scale_color_manual(values=c(\"#2D735F\", \"#BF9039\", \"#BF5841\")) +\n theme(aspect.ratio = 1) +\n ggtitle(\"Correlation coeffecients\")\n\np2 <- ggplot(set22_pos, aes(x = `CD45RA_rank`, y = `PTPRC_rank`, color = `RNA_list`)) +\n geom_point() +\n theme_classic(base_size = 22) +\n geom_hline(yintercept = 0, linetype = \"dashed\", color = \"black\") +\n geom_vline(xintercept = 0,linetype = \"dashed\", color = \"black\") +\n scale_color_manual(values=c(\"#2D735F\", \"#BF9039\", \"#BF5841\")) +\n scale_x_log10() +\n scale_y_log10() + \n theme(aspect.ratio = 1) +\n ggtitle(\"Correlation coeffecient ranks\")\n\np1 + p2 + plot_annotation(\n title = paste0(nrow(set22_pos), \" RNA from the top 300 positively correlated with the CD45RA protein and/or its RNA\")) & \n theme(text = element_text(size = 22) )\n\n#--------Negative---------\np1 <- ggplot(set22_neg, aes(x = `CD45RA`, y = `PTPRC`, color = `RNA_list`)) +\n geom_point() +\n theme_classic(base_size = 22) +\n geom_hline(yintercept = 0, linetype = \"dashed\", color = \"black\") +\n geom_vline(xintercept = 0,linetype = \"dashed\", color = \"black\") +\n scale_color_manual(values=c(\"#2D735F\", \"#BF9039\", \"#BF5841\")) +\n scale_x_reverse() +\n scale_y_reverse() +\n theme(aspect.ratio = 1) +\n ggtitle(\"Correlation coeffecients\")\n\np2 <- ggplot(set22_neg, aes(x = `CD45RA_rank`, y = `PTPRC_rank`, color = `RNA_list`)) +\n geom_point() +\n theme_classic(base_size = 22) +\n scale_color_manual(values=c(\"#2D735F\", \"#BF9039\", \"#BF5841\")) +\n scale_x_log10() +\n scale_y_log10() + \n theme(aspect.ratio = 1) +\n ggtitle(\"Correlation coeffecient ranks\")\n\np1 + p2 + plot_annotation(\n title = paste0(nrow(set22_neg), \" RNA from the top 300 negatively correlated with the CD45RA protein and/or its RNA\")) & \n theme(text = element_text(size = 22) )","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-06T21:04:17.453995Z","iopub.execute_input":"2023-04-06T21:04:17.455595Z","iopub.status.idle":"2023-04-06T21:04:19.733369Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# sbst <- filter(set22_pos, CD45RA > 0.1 & PTPRC < 0.05 & PTPRC > -0.05)\n# cat(\"Correlate positively with CD45RA (> 0.1) and around zero with its RNA (PTPRC):\", nrow(sbst) )\n#sbst %>% select(RNA) %>% paste0(collapse = \", \")\n\n# sbst <- filter(set22_pos, CD45RA < -0.1 & PTPRC < 0.05 & PTPRC > -0.05)\n# cat(\"\\n\\nCorrelate negatively with CD45RA (< -0.1) and around zero with its RNA (PTPRC):\", nrow(sbst) )\n#sbst %>% select(RNA) %>% paste0(collapse = \", \")","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-06T21:04:19.735946Z","iopub.execute_input":"2023-04-06T21:04:19.737370Z","iopub.status.idle":"2023-04-06T21:04:19.803202Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#**Clustermaps for each of the 6 RNA lists**\n#fix colors range to have same colors on each heatmap\nbreaksList = seq(-1, 1, by = 0.01)\nmyColors <- c(colorRampPalette(c(\"darkblue\",\"white\"))(100), colorRampPalette(c(\"white\",\"red\"))(100))\n\nprot_names <- c(\"CD45\",\"CD45RO\",\"CD45RA\")\n\nmy_corr_plot <- function(RNA_names, prot_names, title) {\n \n set22 <- make_subset(RNA_names, prot_names)\n RNA_names <- gsub(\".*_\",\"\", RNA_names)\n set22_corr <- cor(select(set22, c(all_of(c(RNA_names, prot_names)))), method = \"spearman\")\n \n plt <- pheatmap(set22_corr,fontsize = 14,\n color = myColors,\n breaks = breaksList, #col = colorRampPalette(c(\"darkblue\",\"white\",\"red\"))(200),\n main = title)\n #return(plt)\n}","metadata":{"execution":{"iopub.status.busy":"2023-04-06T21:04:19.805771Z","iopub.execute_input":"2023-04-06T21:04:19.807222Z","iopub.status.idle":"2023-04-06T21:04:19.831211Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sbst <- filter(set22_pos, CD45RA > 0.1 & PTPRC < -0.1)\ncat(\"RNA in the clastermap:\", nrow(sbst) )\nsbst %>% select(RNA) %>% paste0(collapse = \", \")\n\nfig(10,10)\nRNA_names <- c(sbst$RNA, \"ENSG00000081237_PTPRC\")\n\nmy_corr_plot(RNA_names, prot_names, title = \n paste(length(RNA_names)-1,\n \"RNA, positively correlated\\nwith CD45RA (> 0.1) and negatively (< -0.1) with its RNA (PTPRC)\"))","metadata":{"execution":{"iopub.status.busy":"2023-04-06T21:58:07.063897Z","iopub.execute_input":"2023-04-06T21:58:07.065921Z","iopub.status.idle":"2023-04-06T21:58:09.249256Z"},"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sbst <- filter(set22_pos, CD45RA < -0.1 & PTPRC > 0.1)\ncat(\"RNA in the clastermap:\", nrow(sbst) )\nsbst %>% select(RNA) %>% paste0(collapse = \", \")\n\nfig(14,14)\nRNA_names <- c(sbst$RNA, \"ENSG00000081237_PTPRC\")\n\nmy_corr_plot(RNA_names, prot_names, title = \n paste(length(RNA_names)-1,\n \"RNA, negatively correlated\\nwith CD45RA (< -0.1) and positively (> 0.1) with its RNA (PTPRC)\"))\n","metadata":{"execution":{"iopub.status.busy":"2023-04-06T21:59:42.563014Z","iopub.execute_input":"2023-04-06T21:59:42.564800Z","iopub.status.idle":"2023-04-06T21:59:45.827544Z"},"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"
\n 💡 Note: \n \n
","metadata":{}},{"cell_type":"markdown","source":"\n#

- By cell type, day and donor

","metadata":{}},{"cell_type":"markdown","source":"
\n 📌 Task: \n
\n Examine how different are the assosoations of CD45RA and its RNA with other RNAs in different cell types, cells, donors.\n
\n
","metadata":{}},{"cell_type":"markdown","source":"**Calculate correlations between CD45RA and its RNA vs all RNA by cell type day, and donor**","metadata":{}},{"cell_type":"code","source":"gc() ; \n# function fo calculate RNA-RNA correlations by cell type\n#test_genes <- colnames(smat_RNA) #c(\"ENSG00000081237_PTPRC\", \"ENSG00000121410_A1BG\")\n\ncorr_by_fct <- function(factor_col) {\n \n RNA_cor_by_ct <- NULL\n chunk_size <- 1000\n # Loop over factor levels\n for(fct in unique(metadata[, factor_col])) {\n\n #fct = \"BP\"\n cells <- rownames(metadata %>% filter(.data[[factor_col]] == fct))\n\n # Preallocate correlation matrix for current cell type\n cor_by_fct <- matrix(nrow = length(gene_names), \n ncol = ncol(smat_RNA)+2)\n \n # Loop over columns of the sparse matrix with all RNAs\n for(i in seq(1, ncol(smat_RNA), chunk_size)) {\n\n # Select a chunk of columns\n chunk <- i:min(i+chunk_size-1, ncol(smat_RNA))\n\n # Select only the columns with non-zero values\n has_nonzero <- apply(smat_RNA[cells, chunk], 2, function(x) any(x != 0))\n chunk_has_nonzero <- chunk[has_nonzero]\n\n if (length(chunk_has_nonzero) > 0) {\n cor_chunk <- SparseSpearmanCor2(my_genes[cells, ], smat_RNA[cells, chunk_has_nonzero])\n colnames(cor_chunk) <- colnames(smat_RNA)[chunk_has_nonzero]\n cor_by_fct[, chunk_has_nonzero+1] <- cor_chunk\n }\n\n # Fill in NAs for the columns without non-zero values\n chunk_no_nonzero <- chunk[!has_nonzero]\n cor_by_fct[, chunk_no_nonzero+1] <- NA\n }\n # Add cell type column and store in RNA_cor_by_ct matrix\n cor_by_fct[, 1] <- fct\n cor_by_fct[, ncol(cor_by_fct)] <- length(cells)\n colnames(cor_by_fct) <- c(\"group\", colnames(smat_RNA), \"sample size\")\n \n\n RNA_cor_by_ct <- rbind(RNA_cor_by_ct, cor_by_fct)\n\n } \n RNA_cor_by_ct <- cbind(\"Name\" = rep(gene_names, times = n_distinct(metadata[, factor_col])),\n RNA_cor_by_ct)\n return(RNA_cor_by_ct)\n}","metadata":{"_kg_hide-output":true,"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-06T21:04:23.603237Z","iopub.execute_input":"2023-04-06T21:04:23.604716Z","iopub.status.idle":"2023-04-06T21:04:24.338539Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"metadata$day_donor <- paste0(\"Day: \", metadata$day,\", Donor: \", metadata$donor)\ntic()\nRNA_cor_by_ct <- corr_by_fct(factor_col = \"cell_type\")\nRNA_cor_by_da <- corr_by_fct(factor_col = \"day\")\nRNA_cor_by_do <- corr_by_fct(factor_col = \"donor\")\nRNA_cor_by_dado <- corr_by_fct(factor_col = \"day_donor\")\ntoc()\nRNA_cor_by_fct <- rbind(RNA_cor_all, RNA_cor_by_ct, RNA_cor_by_da, RNA_cor_by_do, RNA_cor_by_dado)\n\nrm(RNA_cor_by_ct, RNA_cor_by_da, RNA_cor_by_do, RNA_cor_all, RNA_cor_by_dado)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-06T21:04:24.341095Z","iopub.execute_input":"2023-04-06T21:04:24.342536Z","iopub.status.idle":"2023-04-06T21:23:11.169026Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# #check calculations - compare with the results of cor function\n# cells <- rownames(metadata[metadata$day == \"2\", ])\n# #test_gene <- \"ENSG00000081237_PTPRC\"\n# test_prot <- c(\"CD45RA\", \"CD45RO\")\n# #test_genes <- c(\"ENSG00000162378_ZYG11B\", \"ENSG00000159840_ZYX\", \"ENSG00000074755_ZZEF1\")\n# test_genes <- c(\"ENSG00000081237_PTPRC\", \"ENSG00000121410_A1BG\", \"ENSG00000175899_A2M\")\n\n# tst <- cor(mat_prot[cells ,test_prot], as.matrix(smat_RNA[cells, test_genes]), method = \"spearman\")\n# tst <- t(tst)\n# colnames(tst) <- paste0(colnames(tst), \"_dense_corr\")\n\n# #sparse_corr <- SparseSpearmanCor2(Matrix(mat_prot[ ,test_prot], sparse = TRUE), smat_RNA[, test_genes])\n# #sparse_corr <- SparseSpearmanCor2(mat_prot[cells ,test_prot], smat_RNA[cells, 1:1000])\n# #rownames(sparse_corr) <- paste0(test_prot, \"_sparse_corr\")\n# #colnames(sparse_corr) <- colnames(smat_RNA[, test_genes]) \n# #tst <- cbind(tst, t(sparse_corr))\n\n# tst <- cbind(tst, \"CD45RA_sparse_loop_corr\" = RNA_cor_by_da[1,test_genes])\n\n# tst <- tst %>% as.data.frame %>% mutate_if(is.character, as.numeric)\n# #tst$sparse_corr <- all_corr[1:1000, test_gene] # to easy for me :)\n# cat(\"Comparison of Spearman correlation coefficients for one RNA\",test_prot,\"\n# calculated in the usual way on a dense matrix (dense_corr)\n# and those calculated on a sparse matrix using the SparseSpearmanCor2 function (sparse_corr).\")\n\n# tst %>% head\n\n# ggplot(as.data.frame(tst), aes(x = CD45RA_dense_corr, y = CD45RA_sparse_loop_corr)) +\n# geom_point() +\n# theme_bw(base_size = 22) +\n# theme(aspect.ratio = 1)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-06T21:23:11.173157Z","iopub.execute_input":"2023-04-06T21:23:11.174826Z","iopub.status.idle":"2023-04-06T21:23:11.191344Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#reshape to long\nall_corr_long <- RNA_cor_by_fct %>% as.data.frame %>%\n pivot_longer(cols = 3:(ncol(.)-1), names_to = \"RNA\", values_to = \"corr.coef\")\nall_corr_long$corr.coef <- as.numeric(all_corr_long$corr.coef)\n\ncat(\"The resulting data frame with correlations between CD45RA and its RNA vs all RNA in all cells and by cell type day, and donor\")\nall_corr_long %>% head","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-06T21:23:11.195157Z","iopub.execute_input":"2023-04-06T21:23:11.197006Z","iopub.status.idle":"2023-04-06T21:23:12.031122Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Calculate p-values**","metadata":{}},{"cell_type":"code","source":"tmp <- NULL\nfor (gr in unique(all_corr_long$group)) {\n for(nm in unique(all_corr_long$Name)) {\n \n sbst <- all_corr_long[all_corr_long$group == gr & all_corr_long$Name == nm, ]\n ss <- unique(sbst$sample_size) %>% as.numeric\n\n #for protein\n sbst$t_stats <- unlist(abs(sbst[ , \"corr.coef\"])* sqrt((ss-2)/(1-abs(sbst[ , \"corr.coef\"])^2)), use.names = FALSE)\n sbst$p.value <- 2*pt(sbst$t_stats, ss-2, lower.tail=FALSE)\n sbst$Bonf_p.value <- p.adjust(sbst$p.value, method = \"bonferroni\")\n\n sbst <- select(sbst, -c(\"t_stats\", \"p.value\"))\n tmp <- rbind(tmp, sbst)\n \n }\n}\nall_corr_long <- tmp\nall_corr_long <- all_corr_long[all_corr_long$RNA != \"ENSG00000081237_PTPRC\", ]\nrm(tmp)\n\n#combine group and sample size in one column\nall_corr_long <- all_corr_long\n\nall_corr_long <- cbind(\n \"Group\" = paste0(all_corr_long$group,\" (n = \", all_corr_long$sample_size, \")\"),\n select(all_corr_long, -c(\"group\", \"sample_size\"))\n)\n\ncat(\"Data frame with correlations and statistics for CD45RA protein and its RNA vs all RNA in all cells and by cell type, day, donor\")\nhead(all_corr_long)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-06T21:23:12.034848Z","iopub.execute_input":"2023-04-06T21:23:12.036495Z","iopub.status.idle":"2023-04-06T21:23:17.262483Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Descriptive Statistics for Top Correlations in Subgroups**","metadata":{}},{"cell_type":"code","source":"all_res <- NULL\nfor (gr in unique(all_corr_long$Group)) {\n for(nm in unique(all_corr_long$Name)) {\n \n sbst <- all_corr_long[all_corr_long$Group == gr & all_corr_long$Name == nm, ] %>% na.omit\n sbst <- sbst[order(-sbst$corr.coef), ]\n \n res <- rbind(\n sbst %>% head(n = 300) %>% mutate(\"Subset\" = \"Top 300 (+)\"),\n sbst %>% tail(n = 300) %>% mutate(\"Subset\" = \"Top 300 (-)\"),\n# sbst %>% filter(corr.coef < quantile(sbst$corr.coef, 0.8) &\n# corr.coef > quantile(sbst$corr.coef, 0.2)) %>%\n# mutate(\"Subset\" = \"Between 20 and 80 persentiles\"),\n sbst %>% filter(corr.coef > 0.1) %>% mutate(\"Subset\" = \"R > 0.1\"),\n sbst %>% filter(corr.coef < -0.1) %>% mutate(\"Subset\" = \"R < -0.1\")\n )\n \n all_res <- rbind(all_res, res)\n }\n}\n\n#d_stats <- all_res %>% group_by(Name, Group, Subset) %>% get_summary_stats(show = c(\"n\", \"min\", \"max\", \"median\", \"mean\", \"sd\"))","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-06T21:23:17.264829Z","iopub.execute_input":"2023-04-06T21:23:17.266164Z","iopub.status.idle":"2023-04-06T21:23:22.167037Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"##

By cell type

","metadata":{}},{"cell_type":"code","source":"show_stats <- all_res %>% \n filter(grepl(paste0(paste(unique(metadata$cell_type), collapse = \"|\"), \"|all\"), Group),\n Name == \"CD45RA\")\ncat(\"Descriptive statistics for correlation coefficients for the top 300 RNAs correlated with CD45RA by cell type.\")\nshow_stats %>%\n group_by(Name, Group, Subset) %>% \n get_summary_stats(corr.coef, show = c(\"n\", \"min\", \"max\", \"median\")) %>%\n mutate(stats = ifelse(grepl(\"Top\", Subset), \n paste0(median, \" (\", min, \" to \", max, \"), \", \"n=\", n ), n)) %>%\n pivot_wider(id_cols = -c(variable, n, min, max, median),\n names_from = \"Subset\", values_from = \"stats\") %>%\n mutate(\"+/- ratio\" = round(as.numeric(`R > 0.1`)/as.numeric(`R < -0.1`), 1)) %>%\n select(c(\"Group\", \"Name\", \"R > 0.1\", \"R < -0.1\", \"+/- ratio\", \"Top 300 (-)\", \"Top 300 (+)\"))","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-06T21:23:22.169551Z","iopub.execute_input":"2023-04-06T21:23:22.170986Z","iopub.status.idle":"2023-04-06T21:23:24.611534Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig(20, 12)\nshow_stats[show_stats$Subset %in% c(\"Top 300 (+)\", \"Top 300 (-)\"), ] %>%\n mutate(yintercept = ifelse(Subset == \"Top 300 (+)\", 0.1, -0.1)) %>%\n ggplot(aes(x = Group, y = corr.coef, fill = Name)) +\n geom_violin() +\n scale_fill_manual(values=c(\"#BF5841\")) +\n geom_boxplot(width=0.1, color=\"black\", alpha=0.2,\n outlier.shape = NA) +\n xlab(\"\") +\n ggtitle(\"Distribution of Correlation Coefficients for the Top 300 RNAs Correlated with CD45RA by Cell Type\") +\n geom_hline(aes(yintercept = yintercept), linetype=\"dashed\", color = \"black\") +\n facet_wrap(~Subset, scales = \"free\", ncol = 1) + \n theme_bw(base_size = 22) +\n theme(legend.position = \"none\")","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-06T21:23:24.613976Z","iopub.execute_input":"2023-04-06T21:23:24.615387Z","iopub.status.idle":"2023-04-06T21:23:26.042629Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"show_stats <- all_res %>% \n filter(grepl(paste0(paste(unique(metadata$cell_type), collapse = \"|\"), \"|all\"), Group),\n Name == \"PTPRC\")\ncat(\"Descriptive statistics for correlation coefficients for the top 300 RNAs correlated with PTPRC by cell type.\")\nshow_stats %>%\n group_by(Name, Group, Subset) %>% \n get_summary_stats(corr.coef, show = c(\"n\", \"min\", \"max\", \"median\")) %>%\n mutate(stats = ifelse(grepl(\"Top\", Subset), \n paste0(median, \" (\", min, \" to \", max, \"), \", \"n=\", n ), n)) %>%\n pivot_wider(id_cols = -c(variable, n, min, max, median),\n names_from = \"Subset\", values_from = \"stats\") %>%\n mutate(\"+/- ratio\" = round(as.numeric(`R > 0.1`)/as.numeric(`R < -0.1`), 1)) %>%\n select(c(\"Group\", \"Name\", \"R > 0.1\", \"R < -0.1\", \"+/- ratio\", \"Top 300 (-)\", \"Top 300 (+)\"))","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-06T21:23:26.045559Z","iopub.execute_input":"2023-04-06T21:23:26.047236Z","iopub.status.idle":"2023-04-06T21:23:28.642143Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig(20, 12)\nshow_stats[show_stats$Subset %in% c(\"Top 300 (+)\", \"Top 300 (-)\"), ] %>%\n mutate(yintercept = ifelse(Subset == \"Top 300 (+)\", 0.1, -0.1)) %>%\n ggplot(aes(x = Group, y = corr.coef, fill = Name)) +\n geom_violin() +\n scale_fill_manual(values=c(\"#BF5841\")) +\n geom_boxplot(width=0.1, color=\"black\", alpha=0.2,\n outlier.shape = NA) +\n xlab(\"\") +\n ggtitle(\"Distribution of Correlation Coefficients for the Top 300 RNAs Correlated with PTPRC by Cell Type\") +\n geom_hline(aes(yintercept = yintercept), linetype=\"dashed\", color = \"black\") +\n facet_wrap(~Subset, scales = \"free\", ncol = 1) + \n theme_bw(base_size = 22) +\n theme(legend.position = \"none\")","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-06T21:23:28.644916Z","iopub.execute_input":"2023-04-06T21:23:28.646265Z","iopub.status.idle":"2023-04-06T21:23:30.033926Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"##

By day and donor

","metadata":{}},{"cell_type":"code","source":"#all_res$Group %>% unique\nshow_groups <- c('all (n = 70988)', '2 (n = 21942)','3 (n = 20901)','4 (n = 28145)',\n '32606 (n = 23986)','13176 (n = 22199)','31800 (n = 24803)')\nshow_stats <- all_res %>% \n filter(Group %in% show_groups, Name == \"CD45RA\") %>%\n mutate(Group = factor(Group, levels = show_groups))\n\ncat(\"Descriptive statistics for correlation coefficients for the top 300 RNAs correlated with CD45RA by day and donor.\")\nshow_stats %>%\n group_by(Name, Group, Subset) %>% \n get_summary_stats(corr.coef, show = c(\"n\", \"min\", \"max\", \"median\")) %>%\n mutate(stats = ifelse(grepl(\"Top\", Subset), \n paste0(median, \" (\", min, \" to \", max, \"), \", \"n=\", n ), n)) %>%\n pivot_wider(id_cols = -c(variable, n, min, max, median),\n names_from = \"Subset\", values_from = \"stats\") %>%\n mutate(\"+/- ratio\" = round(as.numeric(`R > 0.1`)/as.numeric(`R < -0.1`), 1)) %>%\n select(c(\"Group\", \"Name\", \"R > 0.1\", \"R < -0.1\", \"+/- ratio\", \"Top 300 (-)\", \"Top 300 (+)\"))","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-06T21:23:30.037232Z","iopub.execute_input":"2023-04-06T21:23:30.038724Z","iopub.status.idle":"2023-04-06T21:23:31.970377Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig(20, 12)\nshow_stats[show_stats$Subset %in% c(\"Top 300 (+)\", \"Top 300 (-)\"), ] %>%\n mutate(yintercept = ifelse(Subset == \"Top 300 (+)\", 0.1, -0.1)) %>%\n ggplot(aes(x = Group, y = corr.coef, fill = Name)) +\n geom_violin() +\n scale_fill_manual(values=c(\"#BF5841\")) +\n geom_boxplot(width=0.1, color=\"black\", alpha=0.2,\n outlier.shape = NA) +\n xlab(\"\") +\n ggtitle(\"Distribution of Correlation Coefficients for the Top 300 RNAs Correlated with CD45RA by Day and Donor\") +\n geom_hline(aes(yintercept = yintercept), linetype=\"dashed\", color = \"black\") +\n facet_wrap(~Subset, scales = \"free\", ncol = 1) + \n theme_bw(base_size = 22) +\n theme(legend.position = \"none\")","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-06T21:23:31.973064Z","iopub.execute_input":"2023-04-06T21:23:31.974444Z","iopub.status.idle":"2023-04-06T21:23:33.229426Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"show_groups <- c('all (n = 70988)', '2 (n = 21942)','3 (n = 20901)','4 (n = 28145)',\n '32606 (n = 23986)','13176 (n = 22199)','31800 (n = 24803)')\nshow_stats <- all_res %>% \n filter(Group %in% show_groups, Name == \"PTPRC\") %>%\n mutate(Group = factor(Group, levels = show_groups))\n\ncat(\"Descriptive statistics for correlation coefficients for the top 300 RNAs correlated with PTPRC by day and donor.\")\nshow_stats %>%\n group_by(Name, Group, Subset) %>% \n get_summary_stats(corr.coef, show = c(\"n\", \"min\", \"max\", \"median\")) %>%\n mutate(stats = ifelse(grepl(\"Top\", Subset), \n paste0(median, \" (\", min, \" to \", max, \"), \", \"n=\", n ), n)) %>%\n pivot_wider(id_cols = -c(variable, n, min, max, median),\n names_from = \"Subset\", values_from = \"stats\") %>%\n mutate(\"+/- ratio\" = round(as.numeric(`R > 0.1`)/as.numeric(`R < -0.1`), 1)) %>%\n select(c(\"Group\", \"Name\", \"R > 0.1\", \"R < -0.1\", \"+/- ratio\", \"Top 300 (-)\", \"Top 300 (+)\"))","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-06T21:23:33.231812Z","iopub.execute_input":"2023-04-06T21:23:33.233134Z","iopub.status.idle":"2023-04-06T21:23:35.137899Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig(20, 12)\nshow_stats[show_stats$Subset %in% c(\"Top 300 (+)\", \"Top 300 (-)\"), ] %>%\n mutate(yintercept = ifelse(Subset == \"Top 300 (+)\", 0.1, -0.1)) %>%\n ggplot(aes(x = Group, y = corr.coef, fill = Name)) +\n geom_violin() +\n scale_fill_manual(values=c(\"#BF5841\")) +\n geom_boxplot(width=0.1, color=\"black\", alpha=0.2,\n outlier.shape = NA) +\n xlab(\"\") +\n ggtitle(\"Distribution of Correlation Coefficients for the Top 300 RNAs Correlated with PTPRC by Day and Donor\") +\n geom_hline(aes(yintercept = yintercept), linetype=\"dashed\", color = \"black\") +\n facet_wrap(~Subset, scales = \"free\", ncol = 1) + \n theme_bw(base_size = 22) +\n theme(legend.position = \"none\")","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-06T21:23:35.140255Z","iopub.execute_input":"2023-04-06T21:23:35.141617Z","iopub.status.idle":"2023-04-06T21:23:36.413478Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"##

By each day*donor combination

","metadata":{}},{"cell_type":"code","source":"#all_res$Group %>% unique\nshow_stats <- all_res %>% \n filter(grepl(\"Day: \", Group), Name == \"CD45RA\") %>%\n mutate(Group = gsub(\", \", \"\\n\", Group)) %>%\n mutate(Group = gsub(\" \\\\(n\", \"\\n(n\", Group))\n\ncat(\"Descriptive statistics for correlation coefficients for the top 300 RNAs correlated with CD45RA by each day*donor combination.\")\nshow_stats %>%\n group_by(Name, Group, Subset) %>% \n get_summary_stats(corr.coef, show = c(\"n\", \"min\", \"max\", \"median\")) %>%\n mutate(stats = ifelse(grepl(\"Top\", Subset), \n paste0(median, \" (\", min, \" to \", max, \"), \", \"n=\", n ), n)) %>%\n pivot_wider(id_cols = -c(variable, n, min, max, median),\n names_from = \"Subset\", values_from = \"stats\") %>%\n mutate(\"+/- ratio\" = round(as.numeric(`R > 0.1`)/as.numeric(`R < -0.1`), 1)) %>%\n select(c(\"Group\", \"Name\", \"R > 0.1\", \"R < -0.1\", \"+/- ratio\", \"Top 300 (-)\", \"Top 300 (+)\"))","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-06T21:23:36.415951Z","iopub.execute_input":"2023-04-06T21:23:36.417367Z","iopub.status.idle":"2023-04-06T21:23:39.018821Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig(20, 12)\nshow_stats[show_stats$Subset %in% c(\"Top 300 (+)\", \"Top 300 (-)\"), ] %>%\n mutate(yintercept = ifelse(Subset == \"Top 300 (+)\", 0.1, -0.1)) %>%\n ggplot(aes(x = Group, y = corr.coef, fill = Name)) +\n geom_violin() +\n scale_fill_manual(values=c(\"#BF5841\")) +\n geom_boxplot(width=0.1, color=\"black\", alpha=0.2,\n outlier.shape = NA) +\n xlab(\"\") +\n ggtitle(\"Distribution of Correlation Coefficients for the Top 300 RNAs Correlated with CD45RA by each Day*Donor Combination\") +\n geom_hline(aes(yintercept = yintercept), linetype=\"dashed\", color = \"black\") +\n facet_wrap(~Subset, scales = \"free\", ncol = 1) + \n theme_bw(base_size = 22) +\n theme(legend.position = \"none\")","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-06T21:23:39.021435Z","iopub.execute_input":"2023-04-06T21:23:39.022937Z","iopub.status.idle":"2023-04-06T21:23:41.280444Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"show_stats <- all_res %>% \n filter(grepl(\"Day: \", Group), Name == \"PTPRC\") %>%\n mutate(Group = gsub(\", \", \"\\n\", Group)) %>%\n mutate(Group = gsub(\" \\\\(n\", \"\\n(n\", Group))\n\ncat(\"Descriptive statistics for correlation coefficients for the top 300 RNAs correlated with PTPRC by each day*donor combination.\")\nshow_stats %>%\n group_by(Name, Group, Subset) %>% \n get_summary_stats(corr.coef, show = c(\"n\", \"min\", \"max\", \"median\")) %>%\n mutate(stats = ifelse(grepl(\"Top\", Subset), \n paste0(median, \" (\", min, \" to \", max, \"), \", \"n=\", n ), n)) %>%\n pivot_wider(id_cols = -c(variable, n, min, max, median),\n names_from = \"Subset\", values_from = \"stats\") %>%\n mutate(\"+/- ratio\" = round(as.numeric(`R > 0.1`)/as.numeric(`R < -0.1`), 1)) %>%\n select(c(\"Group\", \"Name\", \"R > 0.1\", \"R < -0.1\", \"+/- ratio\", \"Top 300 (-)\", \"Top 300 (+)\")) ","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-06T21:23:41.283041Z","iopub.execute_input":"2023-04-06T21:23:41.284492Z","iopub.status.idle":"2023-04-06T21:23:43.776484Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig(20, 12)\nshow_stats[show_stats$Subset %in% c(\"Top 300 (+)\", \"Top 300 (-)\"), ] %>%\n mutate(yintercept = ifelse(Subset == \"Top 300 (+)\", 0.1, -0.1)) %>%\n ggplot(aes(x = Group, y = corr.coef, fill = Name)) +\n geom_violin() +\n scale_fill_manual(values=c(\"#BF5841\")) +\n geom_boxplot(width=0.1, color=\"black\", alpha=0.2,\n outlier.shape = NA) +\n xlab(\"\") +\n ggtitle(\"Distribution of Correlation Coefficients for the Top 300 RNAs Correlated with PTPRC by each Day*Donor Combination\") +\n geom_hline(aes(yintercept = yintercept), linetype=\"dashed\", color = \"black\") +\n facet_wrap(~Subset, scales = \"free\", ncol = 1) + \n theme_bw(base_size = 22) +\n theme(legend.position = \"none\")","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-06T21:23:43.778845Z","iopub.execute_input":"2023-04-06T21:23:43.780182Z","iopub.status.idle":"2023-04-06T21:23:45.205583Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"
\n 💡 Note: \n \n
","metadata":{}},{"cell_type":"markdown","source":"#

Intersections in RNA lists

","metadata":{}},{"cell_type":"markdown","source":"#

- For CD45RA protein and for CD45 RNA (separately) by cell type, day and donor

","metadata":{}},{"cell_type":"markdown","source":"**Subset RNA positively and negatively correlated with the CD45RA protein and/or its RNA by cell type**","metadata":{}},{"cell_type":"code","source":"all_corr_wide <- NULL\nfor (nm in unique(all_corr_long$Name)) {\n \n nm_sbst <- all_corr_long[all_corr_long$Name == nm, ]\n \n name_corr_wide <- data.frame(\n Name = rep(nm, each = n_distinct(all_corr_long$RNA)))\n for(gr in unique(nm_sbst$Group)) {\n \n \n sbst <- nm_sbst[nm_sbst$Group == gr, ] #%>% na.omit\n sbst <- sbst[order(-abs(sbst$corr.coef)), ] #%>% head(3)\n sbst <- sbst %>% select(\"RNA\") \n names(sbst) <- gr\n name_corr_wide <- cbind(name_corr_wide, sbst)\n }\n all_corr_wide <- rbind(all_corr_wide, name_corr_wide)\n}\n\nrownames(all_corr_wide) <- NULL\nhead(all_corr_wide)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-06T21:48:08.356503Z","iopub.execute_input":"2023-04-06T21:48:08.358197Z","iopub.status.idle":"2023-04-06T21:48:11.552614Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Calculate number of intersections**","metadata":{}},{"cell_type":"code","source":"list_l_seq <- function(max_lenght) {\n \n s <- c(10, 25, 50, 100, 200, 500, 1000, 2500, 5000, 10000,\n ceiling(seq(1, max_lenght,\n length.out = max_lenght/100)))\n s <- s %>% sort %>% unique\n return(s)\n}","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-06T21:48:31.528638Z","iopub.execute_input":"2023-04-06T21:48:31.531370Z","iopub.status.idle":"2023-04-06T21:48:31.556509Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"res_num <- NULL\nfor(nm in unique(all_corr_wide$Name)) {\n \n dat <- all_corr_wide[all_corr_wide$Name == nm, ]\n len <- list_l_seq(n_distinct(all_corr_long$RNA))\n\n #resulting data frame to fill in\n res <- data.frame(data = colnames(dat)[3:ncol(dat)], #1st column - protein 2nd - reference data\n Name = nm)\n new_cols <- paste(rep(\"top_\", length(len)), len, sep = \"\")\n res[ , new_cols] <- NA\n\n #fill the data frame\n for(i in len){\n\n sbst <- dat[1:i, ] \n col_num <- which(grepl(paste0(\"_\", i, \"$\"), colnames(res)))\n\n #for each column after the first with the reference data\n for (d in res$data) {\n\n int <- intersect(sbst[ , 2], sbst[ , d]) %>% length\n res[res$data == d, col_num] <- int\n } \n } \n res_num <- rbind(res_num, res)\n}\n\nhead(res_num, n = 5)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-06T21:48:31.667722Z","iopub.execute_input":"2023-04-06T21:48:31.669280Z","iopub.status.idle":"2023-04-06T21:48:49.121872Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Calculate proportion of intersections**","metadata":{}},{"cell_type":"code","source":"res_prop <- NULL\nfor(nm in unique(all_corr_wide$Name)) {\n \n dat <- all_corr_wide[all_corr_wide$Name == nm, ]\n len <- list_l_seq(n_distinct(all_corr_long$RNA))\n\n #resulting data frame to fill in\n res <- data.frame(data = colnames(dat)[3:ncol(dat)], #1st column - protein 2nd - reference data\n Name = nm)\n new_cols <- paste(rep(\"top_\", length(len)), len, sep = \"\")\n res[ , new_cols] <- NA\n\n #fill the data frame\n for(i in len){\n\n sbst <- dat[1:i, ] \n col_num <- which(grepl(paste0(\"_\", i, \"$\"), colnames(res)))\n\n #for each column after the first with the reference data\n for (d in res$data) {\n\n int <- intersect(sbst[ , 2], sbst[ , d]) %>% length\n res[res$data == d, col_num] <- int/i\n } \n } \n res_prop <- rbind(res_prop, res)\n}\n\nhead(res_prop, n = 5)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-06T21:48:49.257431Z","iopub.execute_input":"2023-04-06T21:48:49.259121Z","iopub.status.idle":"2023-04-06T21:49:06.718566Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt_prop <- res_prop %>% pivot_longer(cols = 3:ncol(res_prop), \n names_to = \"Top_length\",\n values_to = \"Prop_intersections\")\nplt_prop <- mutate(plt_prop, Top_length = as.numeric(gsub(\"top_\", \"\", Top_length)))\n\n#head(plt_prop)\n\nplt_num <- res_num %>% pivot_longer(cols = 3:ncol(res_num), \n names_to = \"Top_length\",\n values_to = \"Number_intersections\")\nplt_num <- mutate(plt_num, Top_length = as.numeric(gsub(\"top_\", \"\", Top_length)))\n#head(plt_num)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-06T21:49:06.722272Z","iopub.execute_input":"2023-04-06T21:49:06.723972Z","iopub.status.idle":"2023-04-06T21:49:06.785295Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"all_res <- plt_num %>% left_join(plt_prop, by = c(\"data\",\"Name\", \"Top_length\"))\nall_res <- filter(all_res, Top_length %in% c(10, 25, 50, 100, 150, 200, 500, max(Top_length)))\n#all_res <- filter(all_res, Top_length %in% c(seq(0, max(Top_length), by = 50), max(Top_length)))\nall_res <- pivot_wider(all_res, id_cols = -c(\"Number_intersections\"), names_from = \"data\", values_from = c( \"Prop_intersections\"))\n\nall_res","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-06T21:49:06.787823Z","iopub.execute_input":"2023-04-06T21:49:06.789223Z","iopub.status.idle":"2023-04-06T21:49:07.037182Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"##

Intersections by Cell Type

","metadata":{}},{"cell_type":"code","source":"fig(25,8)\n# by cell type\n\ndat_num <- plt_num[grepl(paste0(unique(metadata$cell_type), collapse = \"|\"), plt_num$data), ]\ndat_prop <- plt_prop[grepl(paste0(unique(metadata$cell_type), collapse = \"|\"), plt_prop$data), ]\n\nfor(nm in unique(dat_num$Name)) {\n \n sbst_num <- dat_num[dat_num$Name == nm, ]\n sbst_prop <- dat_prop[dat_prop$Name == nm, ]\n \n p1 <- ggplot(sbst_num, aes(x = Top_length, y=Number_intersections, color = data)) +\n geom_line()+\n geom_point(size = 1) +\n ylab(\"Number of intersections\") +\n xlab(\"RNA list length (N)\") +\n theme_bw(base_size = 20) +\n theme(aspect.ratio = 1)\n\n p2 <- ggplot(sbst_prop, aes(x = Top_length, y=Prop_intersections, color = data)) +\n geom_line()+\n geom_point(size = 1) +\n ylab(\"Proportion of intersections\") +\n xlab(\"RNA list length (N)\") +\n theme_bw(base_size = 20) +\n theme(aspect.ratio = 1)\n\n# p3 <- ggplot(filter(sbst_num, Top_length < 101), aes(x = Top_length, y=Number_intersections, color = data)) +\n# geom_line()+\n# geom_point(size = 1) +\n# ylab(\"Number of intersections\") +\n# xlab(\"RNA list length (N)\") +\n# ylim(0,100) +\n# theme_bw(base_size = 20) +\n# theme(aspect.ratio = 1)\n\n# p4 <- ggplot(filter(sbst_prop, Top_length < 101), aes(x = Top_length, y=Prop_intersections, color = data)) +\n# geom_line()+\n# geom_point(size = 1) +\n# ylab(\"Prop of intersections\") +\n# xlab(\"RNA list length (N)\") +\n# theme_bw(base_size = 20) +\n# theme(aspect.ratio = 1)\n\n print(p1 + p2 + \n #p3 + p4 + \n plot_annotation(\n title = paste0(\"Intersection between RNA lists for each cell type and RNA list for all cells: \", nm),\n subtitle = \"RNA are sorted by absolute value of Spearman correlation coefficient\") & \n theme(text = element_text(size = 22) ) )\n \n cat(\"\\n\\n\")\n print(paste0(\"------------\", nm, \"------------\"))\n cat(\"\\n\\n\")\n \n}","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-06T21:49:07.039569Z","iopub.execute_input":"2023-04-06T21:49:07.040927Z","iopub.status.idle":"2023-04-06T21:49:10.905659Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"hsc_corr_wide <- all_corr_wide[, -2]\n\nres_num <- NULL\nfor(nm in unique(hsc_corr_wide$Name)) {\n \n dat <- hsc_corr_wide[hsc_corr_wide$Name == nm, ]\n len <- list_l_seq(n_distinct(all_corr_long$RNA))\n\n #resulting data frame to fill in\n res <- data.frame(data = colnames(dat)[3:ncol(dat)], #1st column - protein 2nd - reference data\n Name = nm)\n new_cols <- paste(rep(\"top_\", length(len)), len, sep = \"\")\n res[ , new_cols] <- NA\n\n #fill the data frame\n for(i in len){\n\n sbst <- dat[1:i, ] \n col_num <- which(grepl(paste0(\"_\", i, \"$\"), colnames(res)))\n\n #for each column after the first with the reference data\n for (d in res$data) {\n\n int <- intersect(sbst[ , 2], sbst[ , d]) %>% length\n res[res$data == d, col_num] <- int\n } \n } \n res_num <- rbind(res_num, res)\n}\n\n\n\nres_prop <- NULL\nfor(nm in unique(hsc_corr_wide$Name)) {\n \n dat <- hsc_corr_wide[hsc_corr_wide$Name == nm, ]\n len <- list_l_seq(n_distinct(all_corr_long$RNA))\n\n #resulting data frame to fill in\n res <- data.frame(data = colnames(dat)[3:ncol(dat)], #1st column - protein 2nd - reference data\n Name = nm)\n new_cols <- paste(rep(\"top_\", length(len)), len, sep = \"\")\n res[ , new_cols] <- NA\n\n #fill the data frame\n for(i in len){\n\n sbst <- dat[1:i, ] \n col_num <- which(grepl(paste0(\"_\", i, \"$\"), colnames(res)))\n\n #for each column after the first with the reference data\n for (d in res$data) {\n\n int <- intersect(sbst[ , 2], sbst[ , d]) %>% length\n res[res$data == d, col_num] <- int/i\n } \n } \n res_prop <- rbind(res_prop, res)\n}\n\nplt_prop_hsc <- res_prop %>% pivot_longer(cols = 3:ncol(res_prop), \n names_to = \"Top_length\",\n values_to = \"Prop_intersections\")\nplt_prop_hsc <- mutate(plt_prop_hsc, Top_length = as.numeric(gsub(\"top_\", \"\", Top_length)))\n\nplt_num_hsc <- res_num %>% pivot_longer(cols = 3:ncol(res_num), \n names_to = \"Top_length\",\n values_to = \"Number_intersections\")\nplt_num_hsc <- mutate(plt_num_hsc, Top_length = as.numeric(gsub(\"top_\", \"\", Top_length)))\n\n\n\n\nfig(25,8)\n# by cell type\n\ndat_num <- plt_num_hsc[grepl(paste0(unique(metadata$cell_type), collapse = \"|\"), plt_num_hsc$data), ]\ndat_prop <- plt_prop_hsc[grepl(paste0(unique(metadata$cell_type), collapse = \"|\"), plt_prop_hsc$data), ]\n\nfor(nm in unique(dat_num$Name)) {\n \n sbst_num <- dat_num[dat_num$Name == nm, ]\n sbst_prop <- dat_prop[dat_prop$Name == nm, ]\n \n p1 <- ggplot(sbst_num, aes(x = Top_length, y=Number_intersections, color = data)) +\n geom_line()+\n geom_point(size = 1) +\n ylab(\"Number of intersections\") +\n xlab(\"RNA list length (N)\") +\n theme_bw(base_size = 20) +\n theme(aspect.ratio = 1)\n\n p2 <- ggplot(sbst_prop, aes(x = Top_length, y=Prop_intersections, color = data)) +\n geom_line()+\n geom_point(size = 1) +\n ylab(\"Proportion of intersections\") +\n xlab(\"RNA list length (N)\") +\n theme_bw(base_size = 20) +\n theme(aspect.ratio = 1)\n\n# p3 <- ggplot(filter(sbst_num, Top_length < 101), aes(x = Top_length, y=Number_intersections, color = data)) +\n# geom_line()+\n# geom_point(size = 1) +\n# ylab(\"Number of intersections\") +\n# xlab(\"RNA list length (N)\") +\n# ylim(0,100) +\n# theme_bw(base_size = 20) +\n# theme(aspect.ratio = 1)\n\n# p4 <- ggplot(filter(sbst_prop, Top_length < 101), aes(x = Top_length, y=Prop_intersections, color = data)) +\n# geom_line()+\n# geom_point(size = 1) +\n# ylab(\"Prop of intersections\") +\n# xlab(\"RNA list length (N)\") +\n# theme_bw(base_size = 20) +\n# theme(aspect.ratio = 1)\n\n print(p1 + p2 + \n #p3 + p4 + \n plot_annotation(\n title = paste0(\"Intersection between RNA lists for each cell type and RNA list for HSC cells: \", nm),\n subtitle = \"RNA are sorted by absolute value of Spearman correlation coefficient\") & \n theme(text = element_text(size = 22) ) )\n \n cat(\"\\n\\n\")\n print(paste0(\"------------\", nm, \"------------\"))\n cat(\"\\n\\n\")\n \n}","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-06T21:49:10.908520Z","iopub.execute_input":"2023-04-06T21:49:10.910115Z","iopub.status.idle":"2023-04-06T21:49:48.305086Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"##

Intersections by Day

","metadata":{}},{"cell_type":"code","source":"fig(25,8)\n# by cell type\n\ndat_num <- plt_num[plt_num$data %in% c('2 (n = 21942)','3 (n = 20901)','4 (n = 28145)'), ]\ndat_prop <- plt_prop[plt_prop$data %in% c('2 (n = 21942)','3 (n = 20901)','4 (n = 28145)'), ]\n\nfor(nm in unique(dat_num$Name)) {\n \n sbst_num <- dat_num[dat_num$Name == nm, ]\n sbst_prop <- dat_prop[dat_prop$Name == nm, ]\n \n p1 <- ggplot(sbst_num, aes(x = Top_length, y=Number_intersections, color = data)) +\n geom_line()+\n geom_point(size = 1) +\n ylab(\"Number of intersections\") +\n xlab(\"RNA list length (N)\") +\n theme_bw(base_size = 20) +\n theme(aspect.ratio = 1)\n\n p2 <- ggplot(sbst_prop, aes(x = Top_length, y=Prop_intersections, color = data)) +\n geom_line()+\n geom_point(size = 1) +\n ylab(\"Proportion of intersections\") +\n xlab(\"RNA list length (N)\") +\n theme_bw(base_size = 20) +\n theme(aspect.ratio = 1)\n\n print(p1 + p2 + \n #p3 + p4 + \n plot_annotation(\n title = paste0(\"Intersection between RNA lists for each day and RNA list for all cells: \", nm),\n subtitle = \"RNA are sorted by absolute value of Spearman correlation coefficient\") & \n theme(text = element_text(size = 22) ) )\n \n cat(\"\\n\\n\")\n print(paste0(\"------------\", nm, \"------------\"))\n cat(\"\\n\\n\")\n \n}","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-06T21:49:48.308000Z","iopub.execute_input":"2023-04-06T21:49:48.310066Z","iopub.status.idle":"2023-04-06T21:49:51.443352Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"##

Intersections By Donor

","metadata":{}},{"cell_type":"code","source":"fig(25,8)\n# by donor\n\ndat_num <- plt_num[plt_num$data %in% c('32606 (n = 23986)','13176 (n = 22199)','31800 (n = 24803)'), ]\ndat_prop <- plt_prop[plt_prop$data %in% c('32606 (n = 23986)','13176 (n = 22199)','31800 (n = 24803)'), ]\n\nfor(nm in unique(dat_num$Name)) {\n \n sbst_num <- dat_num[dat_num$Name == nm, ]\n sbst_prop <- dat_prop[dat_prop$Name == nm, ]\n \n p1 <- ggplot(sbst_num, aes(x = Top_length, y=Number_intersections, color = data)) +\n geom_line()+\n geom_point(size = 1) +\n ylab(\"Number of intersections\") +\n xlab(\"RNA list length (N)\") +\n theme_bw(base_size = 20) +\n theme(aspect.ratio = 1)\n\n p2 <- ggplot(sbst_prop, aes(x = Top_length, y=Prop_intersections, color = data)) +\n geom_line()+\n geom_point(size = 1) +\n ylab(\"Proportion of intersections\") +\n xlab(\"RNA list length (N)\") +\n theme_bw(base_size = 20) +\n theme(aspect.ratio = 1)\n\n print(p1 + p2 + \n #p3 + p4 + \n plot_annotation(\n title = paste0(\"Intersection between RNA lists for each donor and RNA list for all cells: \", nm),\n subtitle = \"RNA are sorted by absolute value of Spearman correlation coefficient\") & \n theme(text = element_text(size = 22) ) )\n \n cat(\"\\n\\n\")\n print(paste0(\"------------\", nm, \"------------\"))\n cat(\"\\n\\n\")\n \n}","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-06T21:49:51.445948Z","iopub.execute_input":"2023-04-06T21:49:51.447417Z","iopub.status.idle":"2023-04-06T21:49:53.816579Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"##

Intersections by all Day*Donor combinations

","metadata":{}},{"cell_type":"code","source":"fig(25,8)\n# by donor\n\ndat_num <- plt_num[grepl(\"Day\", plt_num$data), ]\ndat_prop <- plt_prop[grepl(\"Day\", plt_prop$data), ]\n\nfor(nm in unique(dat_num$Name)) {\n \n sbst_num <- dat_num[dat_num$Name == nm, ]\n sbst_prop <- dat_prop[dat_prop$Name == nm, ]\n \n p1 <- ggplot(sbst_num, aes(x = Top_length, y=Number_intersections, color = data)) +\n geom_line()+\n geom_point(size = 1) +\n ylab(\"Number of intersections\") +\n xlab(\"RNA list length (N)\") +\n theme_bw(base_size = 20) +\n theme(aspect.ratio = 1)\n\n p2 <- ggplot(sbst_prop, aes(x = Top_length, y=Prop_intersections, color = data)) +\n geom_line()+\n geom_point(size = 1) +\n ylab(\"Proportion of intersections\") +\n xlab(\"RNA list length (N)\") +\n theme_bw(base_size = 20) +\n theme(aspect.ratio = 1)\n\n print(p1 + p2 + \n #p3 + p4 + \n plot_annotation(\n title = paste0(\"Intersection between RNA lists for each day and donor and RNA list for all cells: \", nm),\n subtitle = \"RNA are sorted by absolute value of Spearman correlation coefficient\") & \n theme(text = element_text(size = 22) ) )\n \n cat(\"\\n\\n\")\n print(paste0(\"------------\", nm, \"------------\"))\n cat(\"\\n\\n\")\n \n}","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-06T21:49:53.819679Z","iopub.execute_input":"2023-04-06T21:49:53.821268Z","iopub.status.idle":"2023-04-06T21:49:57.167634Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#

- For CD45RA protein and its RNA in different cell types, days and donors

","metadata":{}},{"cell_type":"code","source":"#gr = \"all (n = 70988)\"\n#n = 300\n\ntop_corr_analysis <- function(data, n,\n corr.list = c(\"all\", \"pos\", \"neg\"),\n include_vect = NULL) {\n\n comp_by_fct <- data.frame(\n \"Group\" = unique(data$Group),\n \"Common RNA\" = NA,\n \"N of Common RNA\" = NA,\n \"Unique for protein\" = NA,\n \"Unique for RNA\" = NA,\n check.names = FALSE\n )\n \n for (gr in unique(all_corr_long$Group)) {\n \n top_prot <- filter(data, Group == gr & Name == \"CD45RA\")\n top_RNA <- filter(data, Group == gr & Name == \"PTPRC\")\n \n\n if(!is.null(include_vect)) {\n top_prot <- top_prot %>% filter(RNA %in% include_vect)\n top_RNA <- top_RNA %>% filter(RNA %in% include_vect)\n } \n \n if (corr.list == \"all\") {\n \n top_prot <- top_prot[order(-abs(top_prot$corr.coef)), ] %>% head(n)\n top_RNA <- top_RNA[order(-abs(top_RNA$corr.coef)), ] %>% head(n)\n \n } else if(corr.list == \"pos\") {\n \n top_prot <- top_prot[order(-top_prot$corr.coef), ] %>% head(n)\n top_RNA <- top_RNA[order(-top_RNA$corr.coef), ] %>% head(n)\n \n } else if (corr.list == \"neg\"){\n \n top_prot <- top_prot[order(top_prot$corr.coef), ] %>% head(n)\n top_RNA <- top_RNA[order(top_RNA$corr.coef), ] %>% head(n)\n \n } \n \n com <- top_prot[top_prot$RNA %in% intersect(top_prot$RNA, top_RNA$RNA), ]\n un_prot <- top_prot[top_prot$RNA %in% setdiff(top_prot$RNA, top_RNA$RNA), ]\n un_rna <- top_RNA[top_RNA$RNA %in% setdiff(top_RNA$RNA, top_prot$RNA), ]\n \n if(nrow(com)!=0) com$name_cor <- paste0(com[,\"RNA\"], \": \", com[, \"corr.coef\"])\n if(nrow(un_prot)!=0) un_prot$name_cor <- paste0(un_prot[,\"RNA\"], \": \", un_prot[, \"corr.coef\"])\n if(nrow(un_rna)!=0) un_rna$name_cor <- paste0(un_rna[, \"RNA\"], \": \", un_rna[, \"corr.coef\"])\n \n comp_by_fct[comp_by_fct$Group == gr, ]$`Common RNA` <-\n ifelse(is.null(com$name_cor), NA, list(com$name_cor) )\n comp_by_fct[comp_by_fct$Group == gr, ]$`N of Common RNA` <-\n nrow(com)\n comp_by_fct[comp_by_fct$Group == gr, ]$`Unique for protein` <-\n ifelse(is.null(un_prot$name_cor), NA, list(un_prot$name_cor) )\n comp_by_fct[comp_by_fct$Group == gr, ]$`Unique for RNA` <-\n ifelse(is.null(un_rna$name_cor), NA, list(un_rna$name_cor) )\n }\n return(comp_by_fct)\n}\n\n#reshape to long format\n","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-06T21:49:57.170186Z","iopub.execute_input":"2023-04-06T21:49:57.171670Z","iopub.status.idle":"2023-04-06T21:49:57.188222Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"comp_to_long <- function(comp_by_fct) {\n \n dat <- rbind(\n data.frame(\n Group = rep(comp_by_fct$Group,\n times = ifelse(comp_by_fct$`N of Common RNA` == 0, \n 1, comp_by_fct$`N of Common RNA`)),\n RNA = unlist(comp_by_fct$`Common RNA`),\n RNA_list = \"Common for prot and RNA\"),\n data.frame(\n Group = rep(comp_by_fct$Group,\n times = n-comp_by_fct$`N of Common RNA`),\n RNA = unlist(comp_by_fct$`Unique for protein`),\n RNA_list = \"Unique for protein\"),\n data.frame(\n Group = rep(comp_by_fct$Group,\n times = n-comp_by_fct$`N of Common RNA`),\n RNA = unlist(comp_by_fct$`Unique for RNA`),\n RNA_list = \"Unique for RNA\")\n ) #%>% na.omit\n\n dat <- separate(dat, col = RNA, into = c(\"RNA\", \"corr.coef\"), sep = \": \")\n dat$corr.coef <- as.numeric(dat$corr.coef)\n\n# dat %>%\n# pivot_wider(id_cols = !corr.coef, names_from = \"RNA_list\",\n# values_from = \"RNA\", values_fn = list)\n return(dat)\n \n}","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-06T21:49:57.190750Z","iopub.execute_input":"2023-04-06T21:49:57.192186Z","iopub.status.idle":"2023-04-06T21:49:57.205973Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#cat(\"Common and unique RNA in the top 300 correlated RNA lists\")\nn = 300\ncomp_by_fct <- top_corr_analysis(all_corr_long, corr.list = \"pos\", n = n)\n\ncat(\"The number and proportion of common RNAs in two RNA lists: one for the protein and one for the RNA, each with 300 RNAs, sorted by absolute correlation coefficient.\")\ncomp_by_fct %>% select(c(\"Group\", \"N of Common RNA\")) %>% #filter(grepl(paste(unique(metadata$day_donor), collapse = \"|\"), Group))\n mutate(`Prop of Common RNA` = round(`N of Common RNA`/n, 2))\n\ncomp_by_fct <- comp_to_long(comp_by_fct) %>% na.omit","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-06T21:49:57.208281Z","iopub.execute_input":"2023-04-06T21:49:57.209647Z","iopub.status.idle":"2023-04-06T21:49:59.895525Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-06T21:49:59.897890Z","iopub.execute_input":"2023-04-06T21:49:59.899245Z","iopub.status.idle":"2023-04-06T21:49:59.910770Z"},"trusted":true},"execution_count":null,"outputs":[]}]}