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

RNA and proteins that may be involved in CD45 isoforms regulation in two CITEseq datasets

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

What we do here:

\n
\n The analysis workflow includes comparisons of correlations, dendrograms trees, expression, and scatter plots for RNAs and proteins vs CD45, CD45RA, and CD45RO in two different CITEseq datasets.\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 four healthy human donors.
    2. \n\n**Proteins (Targets dataset):** for the surface protein levels (n = 140), 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 (n = 22050), 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.\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. RDS files with sparse matrices of all RNA-RNA Spearman correlations (calculated here and stored in the dataset) \n
\n
","metadata":{}},{"cell_type":"code","source":"# Libraries\nsuppressPackageStartupMessages({\n \n library(Matrix)\n library(dplyr) \n library(tidyr)\n library(ggplot2)\n library(tictoc) #time measuring\n library(ggpubr) #ggscatter\n library(ggcorrplot) #ggcorrplot\n library(pheatmap) #heatmap\n library(ggvenn) #venn diagrams with labels\n library(psych)\n library(grid)\n library(patchwork)\n library(corrplot)\n library(\"gridExtra\") \n library(visNetwork) #interactive plots\n library(dendextend) #compare dendrograms https://cran.r-project.org/web/packages/dendextend/vignettes/dendextend.html\n\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,"execution":{"iopub.status.busy":"2023-04-23T17:15:51.270366Z","iopub.execute_input":"2023-04-23T17:15:51.303461Z","iopub.status.idle":"2023-04-23T17:15:51.319794Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Load the 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\n\n#All RNA- RNA correlations\n\nall_RNA_corr <- readRDS(\"/kaggle/input/proteinrna-vs-rna-spearman-correlation-data/allRNA_RNA_corr.RDS\")\n\ngc()","metadata":{"_uuid":"051d70d956493feee0c6d64651c6a088724dca2a","_execution_state":"idle","_kg_hide-output":true,"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-23T17:15:51.343906Z","iopub.execute_input":"2023-04-23T17:15:51.345423Z","iopub.status.idle":"2023-04-23T17:16:51.023149Z"},"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-23T17:16:51.025843Z","iopub.execute_input":"2023-04-23T17:16:51.027000Z","iopub.status.idle":"2023-04-23T17:16:51.100881Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Cell types in CITEseq2022 dataset: \nMasP = Mast Cell Progenitor \nMkP = Megakaryocyte Progenitor \nNeuP = Neutrophil Progenitor \nMoP = Monocyte Progenitor \nEryP = Erythrocyte Progenitor \nHSC = Hematopoietic 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":"prot_names <- c(\"CD45\",\"CD45RO\",\"CD45RA\")\n\n\n#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-23T17:16:51.102659Z","iopub.execute_input":"2023-04-23T17:16:51.103626Z","iopub.status.idle":"2023-04-23T17:16:51.116344Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#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\n#-------------ggscatter RNA vs protein--------------------------\n\n#Draw a ggscatter with the specified RNA and protein (prot argument) - data from all cells.\nmy_ggscatter <- function(dataset, prot, RNA, color = \"#033E8C\") {\n \n plt <- suppressMessages(ggscatter(dataset, 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(dataset), \")\"),\n ggtheme = theme_bw(base_size = 22) + theme(aspect.ratio = 1) ) )\n return(plt)\n} \n\n#-------------scatter plots for the two datasets with all cells' data--------------------------\n\n#Draw scatter plots for the two datasets (all cells)\n# for the specified RNAs or proteins (name.x and name.y arguments).\ncomb_paris <- function(name.x, name.y) {\n \n sbst1 <- set22 %>% select(all_of(c(name.x, name.y)))\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(name.x, name.y)))\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(name.x, \"vs\", name.y)\n \n plt <- ggscatter(plt_pair, x = name.x, y = name.y, 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#-------------scatter with facets with cell type in 2022 and all cells from 2023 data --------------------------\n\n#Draw scatter plots for the two datasets (2022 - by cell type and 2023 - all cells)\n# for the specified RNAs or proteins (name.x and name.y arguments).\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}\n\n#-------------correlograms by groups--------------------------\n\n#Plot correlograms for RNAs and proteins by groups specified with factor_column argument\nplot_corr_by_fct <- function(dat, factor_column, RNA_names, prot_names,\n lab_size = 4, tl.cex = 18, show_corr = TRUE) {\n \n fct_col = which(colnames(dat) == factor_column)\n \n all_plt <- list()\n for(fct in unique(dat[, fct_col])) {\n\n sbst <- filter(dat, dat[, fct_col] == fct )\n\n set_corr <- cor(select(sbst, c(all_of(c(RNA_names, prot_names)))), method = \"spearman\")\n \n plt <- ggcorrplot(set_corr, hc.order = TRUE, lab = show_corr,\n title = paste0(factor_column, \": \", fct, \", n = \", nrow(sbst)), lab_size = lab_size,\n ggtheme = theme_minimal(base_size = 24), tl.cex = tl.cex)\n all_plt[[fct]] <- plt\n }\n \n all_corr <- cor(select(dat, c(all_of(c(RNA_names, prot_names)))), method = \"spearman\")\n\n plt <- ggcorrplot(all_corr, hc.order = TRUE, lab = show_corr,\n title = paste0(\"All cells, n = \", nrow(dat)), lab_size = lab_size,\n ggtheme = theme_minimal(base_size = 24), tl.cex = tl.cex)\n all_plt[[\"all\"]] <- plt\n \n return(all_plt)\n}\n\n\n#-------------expression by groups--------------------------\n\n#Plot expression for the set of genes or proteins (gene_set argument) by groups specified with factor_column argument\nplot_expr_by_fct <- function(dat, gene_set, stats = median, factor_column, ncol = 3) {\n \n my_fun <- stats\n plt <- NULL\n for(g in gene_set) {\n \n plt1 <- dat %>% group_by(!! rlang::sym(factor_column)) %>%\n summarise(`Stats for expression` = my_fun(!! rlang::sym(g))) %>%\n mutate(gene_name = g)\n plt_all <- dat %>% summarise(`Stats for expression` = my_fun(!! rlang::sym(g))) %>%\n mutate(gene_name = g)\n plt_all[, factor_column] <- paste0(\"all cells\\nn = \", nrow(dat))\n \n plt <- rbind(plt, plt_all, plt1)\n }\n \n all_zeros <- plt %>% group_by(gene_name) %>% summarise(sum_exppr = sum(`Stats for expression`))\n all_zeros <- all_zeros[all_zeros$sum_exppr == 0, ]\n \n plt <- plt[!plt$gene_name %in% all_zeros$gene_name, ]\n \n out <- \n ggplot(plt, aes(x = !! rlang::sym(factor_column), y = `Stats for expression`) ) +\n geom_bar(stat = 'identity', col=\"grey\", fill = \"#BF9039\") + \n theme_bw(base_size = 18) +\n facet_wrap(~gene_name, scales = \"free\", ncol = ncol) +\n xlab(\"\") +\n ylab(paste0(deparse(substitute(stats)), \" expression\")) +\n labs(title = paste0(deparse(substitute(stats)), \" expression levels\")) \n if (nrow(all_zeros) >0) {\n \n print(paste0(\"Zero median expression in all goups:\") )\n cat(paste0('\"', all_zeros$gene_name, '\"', collapse = \", \"))\n } \n return(out)\n}\n\n#-------------graph network--------------------------\n\n#coloring of edges\ncolfunc_neg <- colorRampPalette(c(\"skyblue\",\"darkblue\"))\ncolfunc_pos <- colorRampPalette(c(\"red\", \"pink\"))\ncolfunc_zero <- colorRampPalette(c(\"lightgray\", \"gray\"))\n\n#common legend for igraphs\nledges <- data.frame(color = c(\"red\", \"blue\", \"gray\"),\n label = c(\"positive\",\"negative\", \"around zero\"))\n\n#Plot graph network for diff_dat dataframe\nprep_igraph <- function(diff_dat, width_coeff) {\n \n dat <- diff_dat %>%\n filter(RNA != \"PTPRC\") %>%\n select(CD45RO, CD45RA, PTPRC, cell_type, RNA) %>%\n pivot_longer(cols = c(\"CD45RO\", \"CD45RA\", \"PTPRC\"),\n names_to = \"Target\", values_to = \"corr.coef\")\n \n \n nodes <- \n rbind(data.frame(label = unique(dat$cell_type),\n group = \"Cell type\",\n shape = \"dot\",\n level = 1),\n data.frame(label = unique(dat$RNA),\n group = \"RNA\",\n shape = \"box\",\n level = 2),\n data.frame(label = unique(dat$Target),\n group = \"Target\",\n shape = \"box\",\n level = 3) ) %>%\n unique %>%\n mutate(id = 1:nrow(.))\n rownames(nodes) <- NULL\n \n edges1 <- dat %>% \n rename(group = cell_type,\n label = RNA) %>%\n mutate(width = 1, corr.coef = 1, Target = NA) %>%\n left_join(select(nodes, c(label, id)), by = c(\"group\" = \"label\")) %>%\n rename(from = id) %>%\n left_join(select(nodes, c(label, id)), by = c(\"label\")) %>% \n rename(to = id) %>% unique\n \n \n edges2 <- dat %>% \n rename(group = cell_type,\n label = RNA) %>%\n mutate(width = abs(corr.coef) * width_coeff) %>%\n left_join(select(nodes, c(label, id)), by = c(\"label\")) %>% \n rename(from = id) %>% \n left_join(select(nodes, c(label, id)), by = c(\"Target\" = \"label\")) %>% \n rename(to = id)\n \n edges <- rbind(edges1, edges2) %>%\n select(from, to, corr.coef, width)\n \n edges <- edges[order(edges$corr.coef, decreasing = TRUE), ]\n edges <- mutate(edges, association = \n ifelse(edges$corr.coef == 1, \"none\", \n ifelse(edges$corr.coef >= 0.1, \"positive\",\n ifelse(edges$corr.coef <= -0.1, \"negative\", \"zero\"))), \n color = NA)\n edges$association <- factor(edges$association, levels = c(\"none\",\"zero\", \"positive\",\"negative\"))\n\n edges[edges$association == \"zero\", ]$color <- \n colfunc_zero(nrow(edges[edges$association == \"zero\", ]))\n edges[edges$association == \"positive\", ]$color <- \n colfunc_pos(nrow(edges[edges$association == \"positive\", ]))\n edges[edges$association == \"negative\", ]$color <- \n colfunc_neg(nrow(edges[edges$association == \"negative\", ]))\n \n out <- list(nodes, edges)\n names(out) <- c(\"nodes\", \"edges\")\n return(out)\n}","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-23T17:16:51.119570Z","iopub.execute_input":"2023-04-23T17:16:51.120624Z","iopub.status.idle":"2023-04-23T17:16:51.156299Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#function to show top correlated proteins with the gene of interest and a heatmap\none_gene_vs_all_prots <- function(dat, gene_name, factor_column) {\n fct_col = which(colnames(dat) == factor_column)\n\n all_corr <- NULL\n for(fct in unique(dat[, fct_col])) {\n\n sbst <- filter(dat, dat[, fct_col] == fct )\n\n set_corr <- cor(sbst[, gene_name],\n select(sbst, c(colnames(mat_prot))), method = \"spearman\") %>%\n as.data.frame\n set_corr$cell_type <- fct\n\n all_corr <- rbind(all_corr, set_corr)\n }\n\n all_cells_corr <- cor(dat[, gene_name],\n select(dat, c(colnames(mat_prot))), method = \"spearman\") %>%\n as.data.frame\n all_cells_corr$cell_type <- \"All cells\"\n all_corr <- rbind(all_cells_corr, all_corr)\n\n fig(25, 7)\n mat_corr <- as.matrix(all_corr[, -141])\n rownames(mat_corr) <- all_corr$cell_type \n\n pheatmap(mat_corr,fontsize = 16,\n color = myColors,\n breaks = breaksList,\n angle_col = 90,\n main = paste0(\"Correlations of the \", gene_name, \" gene with all proteins in the dataset\"))\n\n all_corr <- all_corr %>%\n pivot_longer(id_cols = , \n cols= 1:140,\n names_to = \"Protein\",\n values_to = \"corr.coef\")\n\n top_corr <- NULL\n #gr = 'All cells'\n for(gr in unique(all_corr$cell_type)) {\n\n sbst <- all_corr[all_corr$cell_type == gr, ]\n sbst <- sbst %>%\n filter(abs(sbst$corr.coef) >= 0.1) %>%\n arrange(desc(corr.coef))\n top_corr <- rbind(top_corr, sbst)\n }\n\n cat(paste0(\"Number (n_cor) and names (Proteins) of the proteins with which the \", gene_name, \" RNA correlate with |R| > 0.1 in all cells and by cell types\"))\n top_corr %>% group_by(cell_type) %>%\n summarise(n_corr = n(),\n Proteins = list(paste0(Protein,\" (\", round(corr.coef, 2), \")\")))\n}","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-23T17:16:51.158990Z","iopub.execute_input":"2023-04-23T17:16:51.160324Z","iopub.status.idle":"2023-04-23T17:16:51.171812Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# dat <- set22\n# factor_column = \"cell_type\"\n# fct = \"MasP\"\n# RNA_names <- RNA_names22\n\n#correlograms by groups specified with factor_column argument\ndiff_corr_by_fct <- function(dat, factor_column, RNA_names, prot_names) {\n \n fct_col = which(colnames(dat) == factor_column)\n \n all_dat <- NULL\n for(fct in unique(dat[, fct_col])) {\n\n sbst <- filter(dat, dat[, fct_col] == fct )\n\n set_corr <- cor(select(sbst, c(all_of(c(RNA_names, prot_names)))), method = \"spearman\")\n \n dat_cor <- set_corr %>% as.data.frame %>%\n filter(!rownames(.) %in% prot_names) %>%\n mutate(RNA = rownames(.),\n \"RO\" = ifelse(CD45RO >= 0.1, \"+\",\n ifelse(CD45RO <= -0.1, \"-\", \"0\")),\n \"RA\" = ifelse(CD45RA >= 0.1, \"+\",\n ifelse(CD45RA <= -0.1, \"-\", \"0\")), \n \"CD45_RNA\" = ifelse(PTPRC == 1, NA,\n ifelse(PTPRC >= 0.1, \"+\",\n ifelse(PTPRC <= -0.1, \"-\", \"0\"))),\n \"cell_type\" = fct\n ) %>%\n filter(RO != RA)\n \n dat_cor <- dat_cor[order(dat_cor$RO, dat_cor$RA, dat_cor$CD45_RNA), ]\n dat_cor <- select(dat_cor, c(CD45RO, CD45RA, PTPRC, cell_type, RNA, RO, RA, CD45_RNA))\n rownames(dat_cor) <- NULL\n \n all_dat <- rbind(all_dat, dat_cor)\n\n }\n \n all_corr <- cor(select(dat, c(all_of(c(RNA_names, prot_names)))), method = \"spearman\")\n \n dat_cor <- all_corr %>% as.data.frame %>%\n filter(!rownames(.) %in% prot_names) %>%\n mutate(RNA = rownames(.),\n \"RO\" = ifelse(CD45RO >= 0.1, \"+\",\n ifelse(CD45RO <= -0.1, \"-\", \"0\")),\n \"RA\" = ifelse(CD45RA >= 0.1, \"+\",\n ifelse(CD45RA <= -0.1, \"-\", \"0\")), \n \"CD45_RNA\" = ifelse(PTPRC == 1, NA,\n ifelse(PTPRC >= 0.1, \"+\",\n ifelse(PTPRC <= -0.1, \"-\", \"0\"))),\n \"cell_type\" = \"All cells\"\n ) %>%\n filter(RO != RA)\n \n dat_cor <- dat_cor[order(dat_cor$RO, dat_cor$RA, dat_cor$CD45_RNA), ]\n dat_cor <- select(dat_cor, c(CD45RO, CD45RA, PTPRC,cell_type, RNA, RO, RA, CD45_RNA))\n rownames(dat_cor) <- NULL\n \n all_dat <- rbind(all_dat, dat_cor)\n \n return(all_dat)\n}","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-23T17:16:51.174664Z","iopub.execute_input":"2023-04-23T17:16:51.176011Z","iopub.status.idle":"2023-04-23T17:16:51.188463Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Prepare a subset of the CITEseq 2022 data with all protein-RNA pairs**","metadata":{}},{"cell_type":"code","source":"pairs <- read.csv(\"/kaggle/input/prot-rna-corr-analysis-for-2-citeseq-datasets/CITEseq22_112_prot_RNA_pairs.csv\")\npairs <- cbind(select(pairs[c(1:112), ], name) %>% rename(Protein = name),\n select(pairs[c(113:224), ], name, corr.coef) %>% rename(RNA = name))\npairs <- pairs %>%\n mutate(RNA = sub(\"\\\\).*\", \"\", sub(\".*\\\\(\", \"\", RNA)) )\n\ncat(\"The head of the data ferame with all protein-RNA pairs in the CITEseq22 dataset with Spearman correlation\")\npairs %>% head","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-23T17:16:51.191382Z","iopub.execute_input":"2023-04-23T17:16:51.192928Z","iopub.status.idle":"2023-04-23T17:16:51.262837Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# prepare RNA names for easy comparisons\nmy_RNA_names <- pairs %>% \n mutate(Protein = paste0(Protein, \"_RNA\")) %>%\n rename(My_Name = Protein)\nmy_RNA_names <- filter(my_RNA_names, !My_Name %in% c('CD45RO_RNA', 'CD45RA_RNA') )\ncat(\"The head of the data frame with My RNA Names\")\nmy_RNA_names %>% head","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-23T17:16:51.265479Z","iopub.execute_input":"2023-04-23T17:16:51.266725Z","iopub.status.idle":"2023-04-23T17:16:51.297945Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#make a set with all proteins and their RNA\nmy_RNA <- smat_RNA[, my_RNA_names$RNA] %>% as.matrix\n \n#remove Ensembl Gene IDs for easy comparison between datasets\nfor(i in colnames(my_RNA) ) {\n colnames(my_RNA)[colnames(my_RNA) == i] <- my_RNA_names[my_RNA_names$RNA == i, ]$My_Name\n}\n\nmy_RNA <- merge(metadata, my_RNA, by=0, all.y = TRUE)\nmy_RNA$day <- factor(my_RNA$day, levels = c(\"2\",\"3\",\"4\"))\n\nmy_prot <- mat_prot[, pairs$Protein] %>% as.data.frame %>% mutate(\"Row.names\" = rownames(.) )\n\npairs_set22 <- left_join(my_RNA, my_prot, by = \"Row.names\")\n\n# add sample sizes per cell type\nsample_size_ct <- pairs_set22 %>%\n group_by(cell_type) %>%\n summarize(ss_ct = n(),.groups = \"keep\" )\n\npairs_set22 <- pairs_set22 %>%\n left_join(sample_size_ct, by = c(\"cell_type\")) %>%\n mutate(sample_size_ct = paste0(cell_type, \"\\n\", \"n=\", ss_ct)) %>%\n select(-c(\"ss_ct\"))\n\ncat(\"CITEseq2022: Head of the data frame with\", ncol(pairs_set22)-6, \"RNA/proteins in columns X\", nrow(pairs_set22), \"cells in rows\")\nhead(pairs_set22)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-23T17:16:51.300617Z","iopub.execute_input":"2023-04-23T17:16:51.301847Z","iopub.status.idle":"2023-04-23T17:16:53.212727Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Check protein's RNA expression levels (plots are not shown)**","metadata":{}},{"cell_type":"code","source":"fig(25,35)\ngene_set <- c(my_RNA_names$My_Name)\n\nall_plots <- plot_expr_by_fct(pairs_set22, gene_set, stats = median, factor_column = \"sample_size_ct\", ncol = 3)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-23T17:16:53.215588Z","iopub.execute_input":"2023-04-23T17:16:53.216933Z","iopub.status.idle":"2023-04-23T17:16:55.254904Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# save names of the proteins and RNAs with low expression\nrem_RNA <- c(\"CD101_RNA\", \"CD112_RNA\", \"CD11a_RNA\", \"CD11b_RNA\", \"CD11c_RNA\",\n\"CD122_RNA\", \"CD123_RNA\", \"CD124_RNA\", \"CD127_RNA\",\n\"CD134_RNA\", \"CD137_RNA\", \"CD14_RNA\", \"CD141_RNA\",\n\"CD142_RNA\", \"CD146_RNA\", \"CD154_RNA\", \"CD155_RNA\",\n\"CD161_RNA\", \"CD163_RNA\", \"CD172a_RNA\", \"CD185_RNA\",\n\"CD19_RNA\", \"CD192_RNA\", \"CD194_RNA\", \"CD195_RNA\",\n\"CD1d_RNA\", \"CD2_RNA\", \"CD20_RNA\", \"CD21_RNA\", \"CD22_RNA\",\n\"CD223_RNA\", \"CD224_RNA\", \"CD226_RNA\", \"CD23_RNA\",\n\"CD24_RNA\", \"CD244_RNA\", \"CD25_RNA\", \"CD26_RNA\", \"CD268_RNA\",\n\"CD27_RNA\", \"CD270_RNA\", \"CD274_RNA\", \"CD278_RNA\", \"CD3_RNA\",\n\"CD303_RNA\", \"CD304_RNA\", \"CD319_RNA\", \"CD328_RNA\", \"CD335_RNA\",\n\"CD35_RNA\", \"CD352_RNA\", \"CD38_RNA\", \"CD39_RNA\", \"CD4_RNA\",\n\"CD40_RNA\", \"CD42b_RNA\", \"CD49a_RNA\", \"CD49b_RNA\", \"CD49f_RNA\",\n\"CD54_RNA\", \"CD56_RNA\", \"CD58_RNA\", \"CD62P_RNA\", \"CD64_RNA\",\n\"CD69_RNA\", \"CD7_RNA\", \"CD72_RNA\", \"CD83_RNA\", \"CD85j_RNA\",\n\"CD86_RNA\", \"CD88_RNA\", \"CD93_RNA\", \"CD94_RNA\", \"CD95_RNA\",\n\"CX3CR1_RNA\", \"IgD_RNA\", \"integrinB7_RNA\", \"KLRG1_RNA\")","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-23T17:16:55.256540Z","iopub.execute_input":"2023-04-23T17:16:55.257493Z","iopub.status.idle":"2023-04-23T17:16:55.267563Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Proteins and RNAs with maximal expression in MasP**","metadata":{}},{"cell_type":"code","source":"fig(25,25)\n\ngene_set <- colnames(pairs_set22)[6:227][!colnames(pairs_set22)[6:227] %in% rem_RNA]\ndat <- select(pairs_set22, -all_of(rem_RNA))\nfactor_column <- \"sample_size_ct\"\nplt <- NULL\nfor(g in gene_set) {\n\n plt1 <- dat %>% group_by(!! rlang::sym(factor_column)) %>%\n summarise(`Median expression` = median(!! rlang::sym(g))) %>%\n mutate(gene_name = g)\n plt_all <- dat %>% summarise(`Median expression` = median(!! rlang::sym(g))) %>%\n mutate(gene_name = g)\n plt_all[, factor_column] <- paste0(\"all cells\\nn = \", nrow(dat))\n\n plt <- rbind(plt, plt_all, plt1)\n}\n\nmax_masp_expr <- plt %>% group_by(gene_name) %>% \n filter(`Median expression` == max(`Median expression`))\nmax_masp_expr <- filter(max_masp_expr, sample_size_ct == \"MasP\\nn=8242\")\n\nggplot(plt[plt$gene_name %in% max_masp_expr$gene_name, ],\n aes(x = !! rlang::sym(factor_column), y = `Median expression`) ) +\n geom_bar(stat = 'identity', col=\"grey\", fill = \"#BF9039\") + \n theme_bw(base_size = 18) +\n xlab(\"\") +\n facet_wrap(~gene_name, scales = \"free\", ncol = 3) +\n labs(title = paste0(\"Median expression levels\"))","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-23T17:16:55.270324Z","iopub.execute_input":"2023-04-23T17:16:55.271579Z","iopub.status.idle":"2023-04-23T17:17:00.141222Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#

CD53 vs CD45, CD45RA, CD45RO

","metadata":{}},{"cell_type":"markdown","source":"CD53 protein controls CD45 function and T cell activation and is required for CD45RO expression and mobility (link).\nIn addition, CD53 is shown to stabilize CD45 on the membrane and is required for optimal phosphatase activity and subsequent LCK activation.","metadata":{}},{"cell_type":"markdown","source":"**Subset data from the two datasets**","metadata":{}},{"cell_type":"code","source":"#CITEseq2022\nRNA_pattern <- c(\"PTPRC$\", \"CD53\", \"LCK\")\nRNA_names22 <- colnames(smat_RNA)[which(grepl(paste0(RNA_pattern, collapse = \"|\"),\n colnames(smat_RNA),\n ignore.case = T))] %>% as.character\n\ncat(length(RNA_names22)-1, \"Found genes:\")\nRNA_names22[RNA_names22 != \"ENSG00000081237_PTPRC\"]\n\nset22 <- make_subset(RNA_names22, prot_names)\nRNA_names22 <- gsub(\".*_\",\"\", RNA_names22)\n\ncat(\"CITEseq2022: Head of the data frame with\", ncol(set22)-8, \"RNA/proteins in columns X\", nrow(set22), \"cells in rows\")\nhead(set22)\n\n#CITEseq2023\nRNA_names23 <- colnames(mat_RNA23)[which(grepl(paste0(RNA_pattern, collapse = \"|\"),\n colnames(mat_RNA23),\n ignore.case = T))]\nprot_names23 <- colnames(mat_prot23)[grep(paste(prot_names, collapse = \"|\"), colnames(mat_prot23))]\n\nset23 <- cbind(CD53 = mat_RNA23[, RNA_names23],\n mat_prot23[, prot_names23]) %>% as.data.frame\ncat(\"\\nCITEseq2023: Head of the data frame with\", ncol(set23), \"RNA/proteins in columns X\", nrow(set23), \"cells in rows\")\nhead(set23)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-23T17:17:00.162211Z","iopub.execute_input":"2023-04-23T17:17:00.163819Z","iopub.status.idle":"2023-04-23T17:17:01.172509Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Check expression levels**","metadata":{}},{"cell_type":"code","source":"fig(25,10)\ngene_set <- c(RNA_names22, 'CD45', 'CD45RO', 'CD45RA')\n\nplot_expr_by_fct(set22, gene_set, stats = median, factor_column = \"sample_size_ct\")","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-23T17:17:01.175318Z","iopub.execute_input":"2023-04-23T17:17:01.176628Z","iopub.status.idle":"2023-04-23T17:17:01.934035Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig(7,7)\n\nsum_expr_RNA23 <- apply(select(set23, -grep(\"CD45\", colnames(set23))), 2, median) %>% sort(decreasing = TRUE)\nsum_expr_RNA23 <- data.frame(\"Median_expression\" = sum_expr_RNA23,\n \"Name\" = names(sum_expr_RNA23))\nsum_expr_RNA23$Name <- factor(sum_expr_RNA23$Name, levels = unique(sum_expr_RNA23$Name ))\n\nggplot(sum_expr_RNA23, aes(x = Name, y = Median_expression)) +\n geom_bar(stat = 'identity', fill=\"#BF9039\", col=\"grey\") + \n theme_bw(base_size = 22) +\n xlab(\"\") +\n ggtitle(paste0(\"CITEseq 2023\")) +\n theme(axis.text.x = element_text(angle = 45, hjust=1))\n\ncat(\"Zero median expression:\\n\")\ncat(paste0('\"', sum_expr_RNA23[sum_expr_RNA23$Median_expression == 0,]$Name, '\"', collapse = \", \"))","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-23T17:17:01.935849Z","iopub.execute_input":"2023-04-23T17:17:01.936954Z","iopub.status.idle":"2023-04-23T17:17:02.155292Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Select only well expressed genes**","metadata":{}},{"cell_type":"code","source":"RNA_names22 <- RNA_names22[!RNA_names22 %in% c(\n \"LCK\")]\n\nRNA_names23 <- colnames(set23)[!colnames(set23) %in% c(\n \"LCK\", \"CD53\")]","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-23T17:17:02.157032Z","iopub.execute_input":"2023-04-23T17:17:02.158004Z","iopub.status.idle":"2023-04-23T17:17:02.168423Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### CITEseq22","metadata":{}},{"cell_type":"code","source":"cat(length(prot_names), \"Proteins:\\n\")\ncat(paste0('\"', prot_names, '\"', collapse = \", \"))\ncat(\"\\n\\n\", length(RNA_names22), \"RNAs:\\n\")\ncat(paste0('\"', RNA_names22, '\"', collapse = \", \"))","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-23T17:17:02.171257Z","iopub.execute_input":"2023-04-23T17:17:02.172624Z","iopub.status.idle":"2023-04-23T17:17:02.195631Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig(25, 10)\nall_plt <- plot_corr_by_fct(set22, \"cell_type\", RNA_names22, prot_names)\n\ncat(\"CITEseq 2022, Spearman correlation\")\ndo.call(\"grid.arrange\", c(all_plt, ncol = 4))","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-23T17:17:02.197364Z","iopub.execute_input":"2023-04-23T17:17:02.198418Z","iopub.status.idle":"2023-04-23T17:17:03.963410Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Compare dendrogram trees**","metadata":{}},{"cell_type":"markdown","source":"The function cor.dendlist() is used to compute “Cophenetic” correlation matrix between a list of trees. The value can range between -1 to 1. With near 0 values meaning that the two trees are not statistically similar ([tutorial](http://www.sthda.com/english/articles/28-hierarchical-clustering-essentials/91-comparing-dendrograms-essentials/?utm_content=buffer05e1c&utm_medium=social&utm_source=twitter.com&utm_campaign=buffer)). The all.equal.dendrogram function makes a global comparison of two or more dendrograms trees.","metadata":{}},{"cell_type":"code","source":"dendlst <- dendlist()\nfor(ct in unique(set22$cell_type)) {\n \n sbst <- set22 %>% filter(cell_type == ct) %>%\n select(all_of(RNA_names22), all_of(prot_names)) %>% \n cor(method = \"spearman\") %>% \n dist %>% hclust(\"complete\") %>% as.dendrogram\n dendlst[[ct]] <- sbst\n}\n\ndendlst[[\"All_cells\"]] <- set22 %>%\n select(all_of(RNA_names22), all_of(prot_names)) %>% \n cor(method = \"spearman\") %>% \n dist %>% hclust(\"complete\") %>% as.dendrogram\n\ncors <- cor.dendlist(dendlst)\n# Print correlation matrix\n#round(cors, 2)\nfig(25, 7)\nggcorrplot(cors, hc.order = TRUE, lab = TRUE,\n title = paste0(\"Correlation matrix between dendrograms by cell type\"), lab_size = 6,\n ggtheme = theme_minimal(base_size = 18), tl.cex = 18)\n\n#all.equal(dendlst$MasP, dendlst$All_cells)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-23T17:17:03.965320Z","iopub.execute_input":"2023-04-23T17:17:03.966362Z","iopub.status.idle":"2023-04-23T17:17:04.593755Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Differential correlations with RO and RA isoforms**","metadata":{}},{"cell_type":"code","source":"cat(\"Positive is for R >= 0.1, negative for R <= -0.1 and zero is |R| < 0.1\")\n\ndiff_dat <- diff_corr_by_fct(set22, factor_column = \"cell_type\", RNA_names = RNA_names22, prot_names)\ndiff_dat","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-23T17:17:04.595839Z","iopub.execute_input":"2023-04-23T17:17:04.597035Z","iopub.status.idle":"2023-04-23T17:17:05.034533Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"graph <- prep_igraph(diff_dat, width_coeff = 5)\n\nvisNetwork(graph$nodes, graph$edges, width=\"100%\", height = 700,\n main=paste0(\"Network of RNA differentially correlated with RO and RA isoforms (targets)\")) %>% \n visOptions(highlightNearest = list(enabled = TRUE, algorithm = \"hierarchical\"), selectedBy = \"label\") %>%\n visHierarchicalLayout() %>%\n visLegend(addEdges = ledges)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-23T17:17:05.036438Z","iopub.execute_input":"2023-04-23T17:17:05.037564Z","iopub.status.idle":"2023-04-23T17:17:05.459465Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"
\n 💡 Note:
\n LCK has zero expression in all cell subgroups. There is no CD53 protein in the dataset, but its RNA is well expressed in all cell types except MkP and NeuP.
\n CD53 RNA shows differential correlation with CD45 isoforms in: \n \n The strongest correlations of CD53 are with CD45RA in BP and MoP, but this may be due to small sample sizes.\n
","metadata":{}},{"cell_type":"markdown","source":"### Comparison of the two datasets","metadata":{}},{"cell_type":"code","source":"fig(12,8)\ncomb_paris(name.x = \"CD45RO\", name.y = \"CD53\")","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-23T17:17:05.461536Z","iopub.execute_input":"2023-04-23T17:17:05.462888Z","iopub.status.idle":"2023-04-23T17:17:09.233474Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"protein = \"CD45RO\"\nRNA = \"CD53\"\ncat(RNA, \"vs\", protein)\n\nplot_2ds(protein = protein, RNA = RNA)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-23T17:17:09.235791Z","iopub.execute_input":"2023-04-23T17:17:09.237254Z","iopub.status.idle":"2023-04-23T17:17:13.941368Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig(12,8)\ncomb_paris(name.y = \"CD53\", name.x = \"CD45RA\")","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-23T17:17:13.943242Z","iopub.execute_input":"2023-04-23T17:17:13.944571Z","iopub.status.idle":"2023-04-23T17:17:17.529212Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"protein = \"CD45RA\"\nRNA = \"CD53\"\ncat(RNA, \"vs\", protein)\n\nplot_2ds(protein = protein, RNA = RNA)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-23T17:17:17.531069Z","iopub.execute_input":"2023-04-23T17:17:17.532096Z","iopub.status.idle":"2023-04-23T17:17:22.147837Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Correlation with other proteins","metadata":{}},{"cell_type":"code","source":"fig(25, 12)\n\nRNA_pattern <- c(\"PTPRC$\", \"CD53\")\nRNA_names22_set <- colnames(smat_RNA)[which(grepl(paste0(RNA_pattern, collapse = \"|\"),\n colnames(smat_RNA),\n ignore.case = T))] %>% as.character\n\none_vs_all <- make_subset(RNA_names22_set, colnames(mat_prot))\n\none_gene_vs_all_prots(one_vs_all, gene_name = \"CD53\", factor_column <- \"sample_size_ct\")","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-23T17:17:22.149477Z","iopub.execute_input":"2023-04-23T17:17:22.150463Z","iopub.status.idle":"2023-04-23T17:17:28.772706Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#

Transcription factors vs CD45, CD45RA, CD45RO

","metadata":{}},{"cell_type":"markdown","source":"PU.1, GATA1, c-Myb, Mafk, GATA2, Tal1 and Mxi1 are all well-known regulators of multilineage blood cell specification ([link](http://www.sciencedirect.com/science/article/pii/S0301472X14001969#bib34)). ","metadata":{}},{"cell_type":"markdown","source":"**Subset data from the two datasets**","metadata":{}},{"cell_type":"code","source":"#CITEseq2022\nRNA_pattern <- c(\"PTPRC$\", \"GATA\",\"Tal1\",\"STAT\", \"Myb$\", \"Mafk\", \"Mxi1\",\n \"T-bet\",\n \"BATF$\", \"IRF4\", \"AHR$\",\"c-Maf\",\"IRF4\")\nRNA_names22 <- colnames(smat_RNA)[which(grepl(paste0(RNA_pattern, collapse = \"|\"),\n colnames(smat_RNA),\n ignore.case = T))] %>% as.character\n\ncat(length(RNA_names22)-1, \"Found genes:\")\nRNA_names22[RNA_names22 != \"ENSG00000081237_PTPRC\"]\n\nset22 <- make_subset(RNA_names22, prot_names)\nRNA_names22 <- gsub(\".*_\",\"\", RNA_names22)\n\ncat(\"CITEseq2022: Head of the data frame with\", ncol(set22)-8, \"RNA/proteins in columns X\", nrow(set22), \"cells in rows\")\nhead(set22)\n\n#CITEseq2023\nRNA_names23 <- colnames(mat_RNA23)[which(grepl(paste0(RNA_pattern, collapse = \"|\"),\n colnames(mat_RNA23),\n ignore.case = T))]\nprot_names23 <- colnames(mat_prot23)[grep(paste(prot_names, collapse = \"|\"), colnames(mat_prot23))]\n\nset23 <- cbind(PTPRC = mat_RNA23[, RNA_names23],\n mat_prot23[, prot_names23]) %>% as.data.frame\ncat(\"\\nCITEseq2023: Head of the data frame with\", ncol(set23), \"RNA/proteins in columns X\", nrow(set23), \"cells in rows\")\nhead(set23)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-23T17:17:28.775506Z","iopub.execute_input":"2023-04-23T17:17:28.776868Z","iopub.status.idle":"2023-04-23T17:17:29.801998Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### CITEseq2022","metadata":{}},{"cell_type":"markdown","source":"**Check expression levels**","metadata":{}},{"cell_type":"code","source":"fig(25,20)\ngene_set <- c(RNA_names22[RNA_names22 != \"PTPRC\"])\n\nplot_expr_by_fct(set22, gene_set, factor_column = \"sample_size_ct\")","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-23T17:17:29.803932Z","iopub.execute_input":"2023-04-23T17:17:29.805076Z","iopub.status.idle":"2023-04-23T17:17:31.590129Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Select only well expressed genes**","metadata":{}},{"cell_type":"code","source":"RNA_names22 <- RNA_names22[!RNA_names22 %in% c(\"AHR\", \"GATA2-AS1\", \"GATA3\", \"GATA5\", \"GATAD2B\", \"IRF4\", \"MAFK\", \"MXI1\", \"STAT2\", \"STAT4\", \"STAT5B\")] ","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-23T17:17:31.592170Z","iopub.execute_input":"2023-04-23T17:17:31.593420Z","iopub.status.idle":"2023-04-23T17:17:31.605313Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Correlograms**","metadata":{}},{"cell_type":"code","source":"cat(length(prot_names), \"Proteins:\\n\")\ncat(paste0('\"', prot_names, '\"', collapse = \", \"))\ncat(\"\\n\\n\", length(RNA_names22), \"RNAs:\\n\")\ncat(paste0('\"', RNA_names22, '\"', collapse = \", \"))","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-23T17:17:31.607331Z","iopub.execute_input":"2023-04-23T17:17:31.636442Z","iopub.status.idle":"2023-04-23T17:17:31.662396Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig(25, 35)\nall_plt <- plot_corr_by_fct(set22, \"cell_type\", RNA_names22, prot_names)\n\ncat(\"CITEseq 2022, Spearman correlation\")\ndo.call(\"grid.arrange\", c(all_plt, ncol = 2))","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-23T17:17:31.665908Z","iopub.execute_input":"2023-04-23T17:17:31.667907Z","iopub.status.idle":"2023-04-23T17:17:34.504453Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Compare dendrogram trees**","metadata":{}},{"cell_type":"code","source":"dendlst <- dendlist()\nfor(ct in unique(set22$cell_type)) {\n \n sbst <- set22 %>% filter(cell_type == ct) %>%\n select(all_of(RNA_names22), all_of(prot_names)) %>% \n cor(method = \"spearman\") %>% \n dist %>% hclust(\"complete\") %>% as.dendrogram\n dendlst[[ct]] <- sbst\n}\n\ndendlst[[\"All_cells\"]] <- set22 %>%\n select(all_of(RNA_names22), all_of(prot_names)) %>% \n cor(method = \"spearman\") %>% \n dist %>% hclust(\"complete\") %>% as.dendrogram\n\ncors <- cor.dendlist(dendlst)\n# Print correlation matrix\n#round(cors, 2)\nfig(25, 7)\nggcorrplot(cors, hc.order = TRUE, lab = TRUE,\n title = paste0(\"Correlation matrix between dendrograms by cell type\"), lab_size = 6,\n ggtheme = theme_minimal(base_size = 18), tl.cex = 18)\n\n#all.equal(dendlst$MasP, dendlst$All_cells)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-23T17:17:34.507130Z","iopub.execute_input":"2023-04-23T17:17:34.508699Z","iopub.status.idle":"2023-04-23T17:17:35.543845Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Differential correlations with RO and RA isoforms**","metadata":{}},{"cell_type":"code","source":"diff_dat <- diff_corr_by_fct(set22, factor_column = \"cell_type\", RNA_names = RNA_names22, prot_names)\n\ncat(\"Positive is for R >= 0.1, negative for R <= -0.1 and zero is |R| < 0.1\")\ndiff_dat #%>% filter(!RA == \"-\" )#& cell_type == \"All cells\")","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-23T17:17:35.545643Z","iopub.execute_input":"2023-04-23T17:17:35.546609Z","iopub.status.idle":"2023-04-23T17:17:36.283217Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"graph <- prep_igraph(diff_dat, width_coeff = 3)\n\nvisNetwork(graph$nodes, graph$edges, width=\"100%\", height = 700,\n main=paste0(\"Network of RNA differentially correlated with RO and RA isoforms (targets)\")) %>% \n visOptions(highlightNearest = list(enabled = TRUE, algorithm = \"hierarchical\"), selectedBy = \"label\") %>%\n visHierarchicalLayout() %>%\n visLegend(addEdges = ledges)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-23T17:17:36.285901Z","iopub.execute_input":"2023-04-23T17:17:36.286988Z","iopub.status.idle":"2023-04-23T17:17:36.486271Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"
\n 💡 Note: \n The most common differential correlation is negative with CD45RA and positive or around zero with CD45RO:

\n By genes:\n \n By cell type\n\n Positive with CD45RA and around zero with CD45RO:\n \nThe strongest correlations of CD53 are with CD45RA in BP and MoP, but this may be due to small sample sizes.\n \n
","metadata":{}},{"cell_type":"markdown","source":"### Correlation with other proteins","metadata":{}},{"cell_type":"code","source":"fig(25, 12)\n\nRNA_pattern <- c(\"PTPRC$\", \"_GATA2$\")\nRNA_names22_set <- colnames(smat_RNA)[which(grepl(paste0(RNA_pattern, collapse = \"|\"),\n colnames(smat_RNA),\n ignore.case = T))] %>% as.character\n\none_vs_all <- make_subset(RNA_names22_set, colnames(mat_prot))\n\none_gene_vs_all_prots(one_vs_all, gene_name = \"GATA2\", factor_column <- \"sample_size_ct\")","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-23T17:17:36.489322Z","iopub.execute_input":"2023-04-23T17:17:36.491239Z","iopub.status.idle":"2023-04-23T17:17:42.463028Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#

HNRNP**: Alternative Splicing of CD45

","metadata":{}},{"cell_type":"markdown","source":"CD45 transcripts undergo extensive alternative splicing in which exons 4, 5, and 6 are variably excluded. Primary naive T cells and B cells express the larger isoforms and are referred to as RA+. In contrast, activated and memory T cells express the shortest isoform, CD45RO. Unlike T cells, which readily undergo CD45 alternative splicing upon activation, B cells are refractory to the CD45RA > RO switch. The data confirm an important role for hnRNPLL in CD45 alternative splicing in T cells and also in B cells. Reference: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC2791692/","metadata":{}},{"cell_type":"markdown","source":"**Subset data from the two datasets**","metadata":{}},{"cell_type":"code","source":"#CITEseq2022\nRNA_pattern <- c(\"PTPRC$\", \"hnRNP\") #,\" ENSG00000143889\"\nRNA_names22 <- colnames(smat_RNA)[which(grepl(paste0(RNA_pattern, collapse = \"|\"),\n colnames(smat_RNA),\n ignore.case = T))] %>% as.character\n\ncat(length(RNA_names22)-1, \"Found genes:\")\nRNA_names22[RNA_names22 != \"ENSG00000081237_PTPRC\"]\n\nset22 <- make_subset(RNA_names22, prot_names)\nRNA_names22 <- gsub(\".*_\",\"\", RNA_names22)\n\ncat(\"CITEseq2022: Head of the data frame with\", ncol(set22)-8, \"RNA/proteins in columns X\", nrow(set22), \"cells in rows\")\nhead(set22)\n\n#CITEseq2023\nRNA_names23 <- colnames(mat_RNA23)[which(grepl(paste0(RNA_pattern, collapse = \"|\"),\n colnames(mat_RNA23),\n ignore.case = T))]\nprot_names23 <- colnames(mat_prot23)[grep(paste(prot_names, collapse = \"|\"), colnames(mat_prot23))]\n\nset23 <- cbind(mat_RNA23[, RNA_names23],\n mat_prot23[, prot_names23]) %>% as.data.frame\ncat(\"CITEseq2023: Head of the data frame with\", ncol(set23), \"RNA/proteins in columns X\", nrow(set23), \"cells in rows\")\nhead(set23)\n","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-23T17:17:42.465895Z","iopub.execute_input":"2023-04-23T17:17:42.467792Z","iopub.status.idle":"2023-04-23T17:17:43.514742Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Check expression levels**","metadata":{}},{"cell_type":"code","source":"fig(25,30)\ngene_set <- c(RNA_names22[RNA_names22 != \"PTPRC\"])\n\nplot_expr_by_fct(set22, gene_set, factor_column = \"sample_size_ct\")","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-23T17:17:43.516384Z","iopub.execute_input":"2023-04-23T17:17:43.517353Z","iopub.status.idle":"2023-04-23T17:17:46.695420Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig(7,7)\n\nsum_expr_RNA23 <- apply(select(set23, grep(\"HNRNP|PTPRC\", colnames(set23))), 2, median) %>% sort(decreasing = TRUE)\nsum_expr_RNA23 <- data.frame(\"Median_expression\" = sum_expr_RNA23,\n \"Name\" = names(sum_expr_RNA23))\nsum_expr_RNA23$Name <- factor(sum_expr_RNA23$Name, levels = unique(sum_expr_RNA23$Name ))\n#hline <- sum_expr_RNA23[sum_expr_RNA23$Name == \"PTPRC\",]$Median_expression\n\nggplot(sum_expr_RNA23, aes(x = Name, y = Median_expression)) +\n geom_bar(stat = 'identity', fill=\"#BF9039\", col=\"grey\") + \n# geom_hline(yintercept=sum_expr_RNA23[sum_expr_RNA23$Name == \"PTPRC\",]$Median_expression,\n# linetype=\"dashed\")+\n theme_bw(base_size = 22) +\n xlab(\"\") +\n ggtitle(paste0(\"CITEseq 2023\")) +\n theme(axis.text.x = element_text(angle = 45, hjust=1))\n\ncat(\"Zero median expression:\\n\")\ncat(paste0('\"', sum_expr_RNA23[sum_expr_RNA23$Median_expression == 0,]$Name, '\"', collapse = \", \"))","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-23T17:17:46.698062Z","iopub.execute_input":"2023-04-23T17:17:46.699506Z","iopub.status.idle":"2023-04-23T17:17:46.925910Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Select only well expressed genes**","metadata":{}},{"cell_type":"code","source":"RNA_names22 <- RNA_names22[!RNA_names22 %in% c(\n \"HNRNPA1L2\", \"HNRNPA1P10\", \"HNRNPA1P14\", \"HNRNPA1P16\", \"HNRNPA1P34\", \"HNRNPA1P40\",\n \"HNRNPA1P5\", \"HNRNPA1P50\", \"HNRNPA1P54\", \"HNRNPA1P68\", \"HNRNPA1P76\", \"HNRNPA3P11\",\n \"HNRNPA3P2\", \"HNRNPA3P6\", \"HNRNPCP1\", \"HNRNPCP2\", \"HNRNPCP4\", \"HNRNPCP6\", \"HNRNPH1P1\",\n \"HNRNPKP2\", \"HNRNPKP4\", \"HNRNPMP1\", \"HNRNPRP1\")]\ncat(length(RNA_names22), \"RNA remained in CITEseq22\")\n\n#save RNA names to use after\nHNRNP_names22 <- RNA_names22[RNA_names22 != \"PTPRC\"]\n\nRNA_names23 <- colnames(set23)[!colnames(set23) %in% c(\n \"HNRNPU\", \"HNRNPA3\", \"HNRNPA0\", \"HNRNPH1\", \"HNRNPC\")]\ncat(\"\\n\", length(RNA_names23), \"RNA remained in CITEseq23\")","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-23T17:17:46.928312Z","iopub.execute_input":"2023-04-23T17:17:46.929726Z","iopub.status.idle":"2023-04-23T17:17:46.952753Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### CITEseq22","metadata":{}},{"cell_type":"code","source":"cat(length(prot_names), \"Proteins:\\n\")\ncat(paste0('\"', prot_names, '\"', collapse = \", \"))\ncat(\"\\n\\n\", length(RNA_names22), \"RNAs:\\n\")\ncat(paste0('\"', RNA_names22, '\"', collapse = \", \"))","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-23T17:17:46.954830Z","iopub.execute_input":"2023-04-23T17:17:46.956459Z","iopub.status.idle":"2023-04-23T17:17:46.976863Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig(25, 45)\nall_plt <- plot_corr_by_fct(set22, \"cell_type\", RNA_names22, prot_names, show_cor = TRUE)\n\ncat(\"CITEseq 2022, Spearman correlation\")\ndo.call(\"grid.arrange\", c(all_plt, ncol = 2))","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-23T17:17:46.978514Z","iopub.execute_input":"2023-04-23T17:17:46.979570Z","iopub.status.idle":"2023-04-23T17:17:50.865748Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Compare dendrogram trees**","metadata":{}},{"cell_type":"code","source":"dendlst <- dendlist()\nfor(ct in unique(set22$cell_type)) {\n \n sbst <- set22 %>% filter(cell_type == ct) %>%\n select(all_of(RNA_names22), all_of(prot_names)) %>% \n cor(method = \"spearman\") %>% \n dist %>% hclust(\"complete\") %>% as.dendrogram\n dendlst[[ct]] <- sbst\n}\n\ndendlst[[\"All_cells\"]] <- set22 %>%\n select(all_of(RNA_names22), all_of(prot_names)) %>% \n cor(method = \"spearman\") %>% \n dist %>% hclust(\"complete\") %>% as.dendrogram\n\ncors <- cor.dendlist(dendlst)\n# Print correlation matrix\n#round(cors, 2)\nfig(25, 7)\nggcorrplot(cors, hc.order = TRUE, lab = TRUE,\n title = paste0(\"Correlation matrix between dendrograms by cell type\"), lab_size = 6,\n ggtheme = theme_minimal(base_size = 18), tl.cex = 18)\n\n#all.equal(dendlst$MasP, dendlst$All_cells)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-23T17:17:50.867936Z","iopub.execute_input":"2023-04-23T17:17:50.869374Z","iopub.status.idle":"2023-04-23T17:17:52.420493Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"HNRNP_clust <- RNA_names22[!RNA_names22 %in% c('HNRNPA1P48','HNRNPA1','HNRNPUL1','HNRNPH1','HNRNPUL2','HNRNPLL','HNRNPL','PTPRC')]\ncat(\"14 genes, forming a cluster in all cells:\\n\")\ncat(paste0('\"', HNRNP_clust, '\"', collapse = \", \"))","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-23T17:17:52.422327Z","iopub.execute_input":"2023-04-23T17:17:52.423461Z","iopub.status.idle":"2023-04-23T17:17:52.438760Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The gene function is taken from Table 1 in https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4947485/ unless otherwise indicated: \n* 'HNRNPA0' - Involved in post-transcriptional regulation of cytokines mRNAs (https://www.genecards.org/cgi-bin/carddisp.pl?gene=HNRNPA0) \n* 'HNRNPA2B1' - Splicing \n* 'HNRNPA3' - paralog of HNRNPA0, Plays a role in cytoplasmic trafficking of RNA. Binds to the cis-acting response element, A2RE. May be involved in pre-mRNA splicing (https://www.genecards.org/cgi-bin/carddisp.pl?gene=HNRNPA3&keywords=HNRNPA3) \n* 'HNRNPAB' - it is not a member of the HNRNP A/B subfamily of HNRNPs, but groups together closely with HNRNPD/AUF1 and HNRNPDL (https://en.wikipedia.org/wiki/HNRNPAB)\n* 'HNRNPC' - May play a role in the early steps of spliceosome assembly and pre-mRNA splicing (https://www.genecards.org/cgi-bin/carddisp.pl?gene=HNRNPC&keywords=HNRNPC)\n* 'HNRNPD' - mRNA decay, Telomere maintenance\n* 'HNRNPDL' - Acts as a transcriptional regulator. Promotes transcription repression (https://www.genecards.org/cgi-bin/carddisp.pl?gene=HNRNPDL&keywords=HNRNPDL)\n* 'HNRNPF' - Splicing\n* 'HNRNPH2' - ?\n* 'HNRNPH3' - Involved in the splicing process and participates in early heat shock-induced splicing arrest (https://www.genecards.org/cgi-bin/carddisp.pl?gene=HNRNPH3&keywords=HNRNPH3)\n* 'HNRNPK' - Translational regulation, Transcriptional regulation, mRNA stability, Splicing\n* 'HNRNPM' - Splicing\n* 'HNRNPR' - Transcriptional regulation\n* 'HNRNPU' - Splicing, Transcriptional regulation","metadata":{}},{"cell_type":"markdown","source":"**Differential correlations with RO and RA isoforms**","metadata":{}},{"cell_type":"code","source":"cat(\"Positive is for R >= 0.1, negative for R <= -0.1 and zero is |R| < 0.1\")\ndiff_dat <- diff_corr_by_fct(set22, factor_column = \"cell_type\", RNA_names = RNA_names22, prot_names)\ndiff_dat # %>% filter(cell_type == \"NeuP\")","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-23T17:17:52.440474Z","iopub.execute_input":"2023-04-23T17:17:52.441535Z","iopub.status.idle":"2023-04-23T17:17:53.608455Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"graph <- prep_igraph(diff_dat, width_coeff = 5)\n\nvisNetwork(graph$nodes, graph$edges, width=\"100%\", height = 700,\n main=paste0(\"Network of RNA differentially correlated with RO and RA isoforms (targets)\")) %>% \n visOptions(highlightNearest = list(enabled = TRUE, algorithm = \"hierarchical\"), selectedBy = \"label\") %>%\n visHierarchicalLayout() %>%\n visLegend(addEdges = ledges)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-23T17:17:53.610286Z","iopub.execute_input":"2023-04-23T17:17:53.611364Z","iopub.status.idle":"2023-04-23T17:17:53.840765Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"
\n 💡 Note:
\n The most pronounced opposite correlations:
\n \n Differential correlation: positive with CD45RA and around zero with CD45RO:
\n\n Differential correlation: negative with CD45RO and around zero with CD45RA:
\n\n Differential correlation: negative with CD45RA and around zero with CD45RO:
\n\nThe strongest correlations of CD53 are with CD45RA in BP and MoP, but this may be due to small sample sizes.\n \n
","metadata":{}},{"cell_type":"markdown","source":"### Comparison of the two datasets","metadata":{}},{"cell_type":"code","source":"#calculate correlations\nset23 <- select(set23, all_of(RNA_names23))\n\nsbst22 <- set22 %>% select(names(set23))\nset22_corr <- cor(sbst22, method = \"spearman\")\ncat(\"The correlation matrix CITEseq2022\")\nset22_corr\n\nset23_corr <- cor(set23, method = \"spearman\")\ncat(\"\\nThe correlation matrix CITEseq2023\")\nset23_corr","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-23T17:17:53.843007Z","iopub.execute_input":"2023-04-23T17:17:53.844344Z","iopub.status.idle":"2023-04-23T17:17:54.025707Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig(25,7)\nplt1 <- ggcorrplot(set22_corr, hc.order = FALSE, lab = TRUE,\n title = \"CITEseq 2022\", lab_size = 6, p.mat = NULL, #as.matrix(set22_p), insig = \"pch\", pch = 0,\n ggtheme = theme_minimal(base_size = 24), tl.cex = 18)\nplt2 <- ggcorrplot(set23_corr, hc.order = FALSE, lab = TRUE,\n title = \"CITEseq 2023\", lab_size = 6, p.mat = NULL, #set23_p, insig = \"pch\",\n ggtheme = theme_minimal(base_size = 24), tl.cex = 18)\nplt1 + plt2","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-23T17:17:54.027379Z","iopub.execute_input":"2023-04-23T17:17:54.028347Z","iopub.status.idle":"2023-04-23T17:17:54.985771Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Comparison of two dendrograms**","metadata":{}},{"cell_type":"code","source":"fig(25, 7)\n# Compute distance matrix, hierarchical clusterings and create two dendrograms in a list\ndend1 <- set22_corr %>% dist %>% hclust(\"complete\") %>% as.dendrogram\ndend2 <- set23_corr %>% dist %>% hclust(\"complete\") %>% as.dendrogram\ndend_list <- dendlist(dend1, dend2)\n# Align and plot two dendrograms side by side\ndendlist(dend1, dend2) %>%\n untangle(method = \"step1side\") %>% # Find the best alignment layout\n tanglegram(cex_main = 3, lab.cex = 2,\n margin_inner = 15,\n common_subtrees_color_branches = TRUE,\n main_left = \"CITEseq 2022\",\n main_right = \"CITEseq 2023\")\ncor_res <- cor.dendlist(dend_list, method = \"cophenetic\")\ncat(\"Cophenetic correlation between the two dendrograms\", round(cor_res[2,1], 2))\nall.equal(dend1, dend2)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-23T17:17:54.987509Z","iopub.execute_input":"2023-04-23T17:17:54.988513Z","iopub.status.idle":"2023-04-23T17:17:55.259673Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Scatter plots for HNRNPLL vs RA and RO by cell type in CITEseq 2022","metadata":{}},{"cell_type":"code","source":"fig(25,8)\n\np1 <- ggscatter(set22, x = \"CD45RO\", y = \"HNRNPLL\", 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 ggtheme = theme_bw(base_size = 22) + theme(aspect.ratio = 1) ) \n\np2 <- ggscatter(set22, x = \"CD45RA\", y = \"HNRNPLL\", 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 ggtheme = theme_bw(base_size = 22) + theme(aspect.ratio = 1) ) \n\np1 + p2 + plot_annotation(\n title = \"CITEseq2022, Spearman correlation, all cells\") & \ntheme(text = element_text(size = 22) )","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-23T17:17:55.261328Z","iopub.execute_input":"2023-04-23T17:17:55.262317Z","iopub.status.idle":"2023-04-23T17:18:01.339110Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig(25,25)\nggscatter(set22, x = \"CD45RO\", y = \"HNRNPLL\", 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(\"CITEseq2022, Spearman correlation by cell type: CD45RO\"),\n ggtheme = theme_bw(base_size = 26) + theme(aspect.ratio = 1) )","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-23T17:18:01.340793Z","iopub.execute_input":"2023-04-23T17:18:01.341803Z","iopub.status.idle":"2023-04-23T17:18:05.522947Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig(25,25)\nggscatter(set22, x = \"CD45RA\", y = \"HNRNPLL\", 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(\"CITEseq2022, Spearman correlation by cell type: CD45RA\"),\n ggtheme = theme_bw(base_size = 26) + theme(aspect.ratio = 1) )","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-23T17:18:05.524929Z","iopub.execute_input":"2023-04-23T17:18:05.526456Z","iopub.status.idle":"2023-04-23T17:18:09.653275Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Corretation of HNRNPLL with all proteins","metadata":{}},{"cell_type":"code","source":"fig(25, 12)\n\nRNA_pattern <- c(\"PTPRC$\", \"_HNRNPLL$\")\nRNA_names22_set <- colnames(smat_RNA)[which(grepl(paste0(RNA_pattern, collapse = \"|\"),\n colnames(smat_RNA),\n ignore.case = T))] %>% as.character\n\none_vs_all <- make_subset(RNA_names22_set, colnames(mat_prot))\n\none_gene_vs_all_prots(one_vs_all, gene_name = \"HNRNPLL\", factor_column <- \"sample_size_ct\")","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-23T17:18:09.655037Z","iopub.execute_input":"2023-04-23T17:18:09.656006Z","iopub.status.idle":"2023-04-23T17:18:15.906869Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Corretation of the mean of 14 HNRNP** genes with other proteins in CITEseq 2022","metadata":{}},{"cell_type":"code","source":"mean_of_genes <- select(set22, c(\"Row.names\", all_of(HNRNP_clust)))\nmean_of_genes$Mean_14HNRNP <- apply(mean_of_genes[, -1], 1, mean)\n\nall_prots <- merge(metadata, mat_prot, by=0, all.y = TRUE)\nall_prots$day <- factor(my_RNA$day, levels = c(\"2\",\"3\",\"4\"))\n\nmean_of_genes <- left_join(all_prots, select(mean_of_genes, c(\"Row.names\", \"Mean_14HNRNP\")), by = \"Row.names\")\n\n# add sample sizes per cell type\nsample_size_ct <- dat %>%\n group_by(cell_type) %>%\n summarize(ss_ct = n(),.groups = \"keep\" )\n\nmean_of_genes <- mean_of_genes %>%\n left_join(sample_size_ct, by = c(\"cell_type\")) %>%\n mutate(sample_size_ct = paste0(cell_type, \"\\n\", \"n=\", ss_ct)) %>%\n select(-c(\"ss_ct\"))\n\nfig(25, 12)\n\none_gene_vs_all_prots(mean_of_genes, gene_name = \"Mean_14HNRNP\", factor_column <- \"sample_size_ct\")","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-23T17:18:15.908617Z","iopub.execute_input":"2023-04-23T17:18:15.909589Z","iopub.status.idle":"2023-04-23T17:18:22.003189Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Correlation of 14 HNRNP** genes forming the cluster with other RNA","metadata":{}},{"cell_type":"code","source":"set22_corr <- all_RNA_corr\nset22_corr <- set22_corr[ , grepl(paste0(\"_\", paste0(HNRNP_clust, collapse = \"$|\"), \"$\"), colnames(set22_corr))] %>% as.matrix\n\ncat(\"The correlation matrix of the\", length(colnames(set22_corr)), \"HNRNP** genes with \", nrow(set22_corr),\"RNA in the CITEseq2022\")\nset22_corr %>% head","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-23T17:18:22.005093Z","iopub.execute_input":"2023-04-23T17:18:22.006168Z","iopub.status.idle":"2023-04-23T17:18:22.307118Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#calculate sum and mean correlation of each row (each of 21k genes)\nmat_rowsums <- rowSums(set22_corr)\nmat_rowmeans <- rowMeans(set22_corr)\n\n#mat_rowsums %>% head\n#mat_rowsums[mat_rowsums!=0] %>% length\n#mat_rowmeans[abs(mat_rowmeans) > 0.2] %>% length\n\nfig(25, 20)\ntop_corr <- set22_corr[rownames(set22_corr) %in% names(mat_rowmeans[abs(mat_rowmeans) > 0.3]), ]\ntop_corr <- top_corr[!grepl(paste0(RNA_names22[RNA_names22 != 'PTPRC'], collapse = \"$|\"), rownames(top_corr)), ]\n\ncat(\"Heatmap for all\", nrow(top_corr), \"RNA with absolute mean correlation with 14 HNRNP** genes > 0.3\")\n\npheatmap(top_corr, fontsize = 14, fontsize_row = 14,\n color = myColors,\n breaks = breaksList,\n angle_col = 90,\n main = \"\")","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-23T17:18:22.308994Z","iopub.execute_input":"2023-04-23T17:18:22.310030Z","iopub.status.idle":"2023-04-23T17:18:22.806094Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"cutoff <- 0.4\ntop_corr <- set22_corr[rownames(set22_corr) %in% names(mat_rowmeans[abs(mat_rowmeans) > cutoff]), ]\ntop_corr <- top_corr[!grepl(paste0(RNA_names22[RNA_names22 != 'PTPRC'], collapse = \"$|\"), rownames(top_corr)), ]\ncat(\"Mean abs correlation with the 14 HNRNP** genes >\", cutoff, \":\", nrow(top_corr), \"\\n\")\ncat(paste0('\"', gsub(\".*_\", \"\", rownames(top_corr)), '\"', collapse = \", \"))\n\ncutoff <- 0.3\ntop_corr <- set22_corr[rownames(set22_corr) %in% names(mat_rowmeans[mat_rowmeans > cutoff]), ]\ntop_corr <- top_corr[!grepl(paste0(RNA_names22[RNA_names22 != 'PTPRC'], collapse = \"$|\"), rownames(top_corr)), ]\ncat(\"\\n\\nMean correlation with the 14 HNRNP** genes >\", cutoff, \":\", nrow(top_corr), \"\\n\")\ncat(paste0('\"', gsub(\".*_\", \"\", rownames(top_corr)), '\"', collapse = \", \"))\n\ncutoff <- -0.3\ntop_corr <- set22_corr[rownames(set22_corr) %in% names(mat_rowmeans[mat_rowmeans < cutoff]), ]\ntop_corr <- top_corr[!grepl(paste0(RNA_names22[RNA_names22 != 'PTPRC'], collapse = \"$|\"), rownames(top_corr)), ]\ncat(\"\\n\\nMean correlation with the 14 HNRNP** genes <\", cutoff, \":\", nrow(top_corr), \"\\n\")\ncat(paste0('\"', gsub(\".*_\", \"\", rownames(top_corr)), '\"', collapse = \", \"))","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-23T17:18:22.807908Z","iopub.execute_input":"2023-04-23T17:18:22.808991Z","iopub.status.idle":"2023-04-23T17:18:22.846046Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"set22_corr <- as.data.frame(set22_corr)\nset22_corr$RNA <- rownames(set22_corr)\n\n#cat(\"Descriptive statistics of correlations of the\", ncol(set22_corr)-1, \"HNRNP** genes vs All RNA\")\n#set22_corr %>% get_summary_stats(type = \"common\")","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-23T17:18:22.847627Z","iopub.execute_input":"2023-04-23T17:18:22.848574Z","iopub.status.idle":"2023-04-23T17:18:22.864247Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**RNAs with mean correlation with the 14 HNRNPs genes > 0.3 vs CD45 proteins and PTPRC**","metadata":{}},{"cell_type":"code","source":"#CITEseq2022\nRNA_pattern <- c(\"PTPRC\", \"ANP32A\", \"ANP32E\", \"CCT2\", \"CCT5\", \"CCT6A\", \"EIF5A\", \"GSPT1\",\n \"HMGB1\", \"HSP90AA1\", \"HSP90AB1\", \"HSPA4\", \"HSPD1\", \"HSPH1\", \"KPNB1\", \"LMNB1\",\n \"NASP\", \"NCL\", \"NUCKS1\", \"NUDC\", \"PA2G4\", \"PTGES3\", \"SERBP1\", \"SET\", \"SFPQ\",\n \"SMARCA5\", \"SUPT16H\", \"SYNCRIP\", \"TMPO\", \"TPM3\", \"XRCC5\", \"YWHAE$\")\nRNA_names22 <- colnames(smat_RNA)[which(grepl(paste0(\"_\", RNA_pattern, collapse = \"$|\"),\n colnames(smat_RNA),\n ignore.case = T))] %>% as.character\n\nset22 <- make_subset(RNA_names22, prot_names)\nRNA_names22 <- gsub(\".*_\",\"\", RNA_names22)\nRNA_names22 <- RNA_names22[RNA_names22 != \"PTPRC\"]\n\ncat(length(prot_names), \"Proteins:\\n\")\ncat(paste0('\"', prot_names, '\"', collapse = \", \"))\ncat(\"\\n\\n\", length(RNA_names22), \"RNAs:\\n\")\ncat(paste0('\"', RNA_names22, '\"', collapse = \", \"))\n\nset22_corr <- cor(select(set22, all_of(RNA_names22)),\n select(set22, all_of(prot_names), \"PTPRC\"),\n method = \"spearman\")\nset22_corr[abs(set22_corr) < 0.1] <- 0\nset22_corr %>% as.data.frame %>% mutate(RNA = rownames(.)) %>%\n pivot_longer(cols = c(all_of(prot_names), \"PTPRC\")) %>%\n filter(value != 0) %>%\n pivot_wider(names_from = \"RNA\", values_from = \"value\")","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-23T17:18:22.865968Z","iopub.execute_input":"2023-04-23T17:18:22.866927Z","iopub.status.idle":"2023-04-23T17:18:24.876619Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**RNAs with mean correlation with the 14 HNRNPs genes < -0.3 vs CD45 proteins and PTPRC**","metadata":{}},{"cell_type":"code","source":"#CITEseq2022\nRNA_pattern <- c(\"PTPRC\", \"EEF1A1\", \"LRRC75A-AS1\", \"RACK1\", \"RPL10\", \"RPL11\", \"RPL12\", \"RPL13\",\n \"RPL18\", \"RPL18A\", \"RPL19\", \"RPL26\", \"RPL28\", \"RPL3\", \"RPL30\", \"RPL32\",\n \"RPL34\", \"RPL35A\", \"RPL37\", \"RPL39\", \"RPL41\", \"RPL7A\", \"RPLP1\", \"RPS12\",\n \"RPS14\", \"RPS15\", \"RPS15A\", \"RPS18\", \"RPS19\", \"RPS2\", \"RPS24\", \"RPS27\",\n \"RPS27A\", \"RPS28\", \"RPS3\", \"RPS3A\", \"RPS5\", \"RPS6\", \"RPS8\", \"RPS9\")\nRNA_names22 <- colnames(smat_RNA)[which(grepl(paste0(\"_\", RNA_pattern, collapse = \"$|\"),\n colnames(smat_RNA),\n ignore.case = T))] %>% as.character\n\nset22 <- make_subset(RNA_names22, prot_names)\nRNA_names22 <- gsub(\".*_\",\"\", RNA_names22)\nRNA_names22 <- RNA_names22[RNA_names22 != \"PTPRC\"]\n\ncat(length(prot_names), \"Proteins:\\n\")\ncat(paste0('\"', prot_names, '\"', collapse = \", \"))\ncat(\"\\n\\n\", length(RNA_names22), \"RNAs:\\n\")\ncat(paste0('\"', RNA_names22, '\"', collapse = \", \"))\n\nset22_corr <- cor(select(set22, all_of(RNA_names22)),\n select(set22, all_of(prot_names), \"PTPRC\"),\n method = \"spearman\")\nset22_corr[abs(set22_corr) < 0.1] <- 0\n\ncat(\"All RNA correlated with at least one of the CD45 proteins or PTPRC with |R| > 0.1\")\nset22_corr %>% as.data.frame %>% mutate(RNA = rownames(.)) %>%\n pivot_longer(cols = c(all_of(prot_names), \"PTPRC\")) %>%\n filter(value != 0) %>%\n pivot_wider(names_from = \"RNA\", values_from = \"value\")","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-23T17:18:24.878321Z","iopub.execute_input":"2023-04-23T17:18:24.879305Z","iopub.status.idle":"2023-04-23T17:18:26.895187Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Corretation of HNRNP** genes with protein-RNA pairs in 2022 dataset","metadata":{}},{"cell_type":"markdown","source":"**Make a set with genes of interest and merge with the pairs_set22 (all proteins with an RNA pair form the CITEseq22 dataset)**","metadata":{}},{"cell_type":"code","source":"HNRNP_names22 <- colnames(smat_RNA)[which(grepl(paste0(HNRNP_names22, collapse = \"$|\"),\n colnames(smat_RNA),\n ignore.case = T))] %>% as.character #%>% length\nHNRNP_set <- smat_RNA[, HNRNP_names22] %>% as.matrix\ncolnames(HNRNP_set) <- gsub(\".*_\", \"\", colnames(HNRNP_set))\n\nHNRNP_set <- HNRNP_set %>% as.data.frame %>% mutate(\"Row.names\" = rownames(.) )\n\nHNRNP_set <- left_join(pairs_set22, HNRNP_set, by = \"Row.names\")\ncat(\"CITEseq2022: Head of the data frame with\", ncol(HNRNP_set)-6, \"RNA/proteins in columns X\", nrow(HNRNP_set), \"cells in rows\")\nhead(HNRNP_set)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-23T17:18:26.896821Z","iopub.execute_input":"2023-04-23T17:18:26.897753Z","iopub.status.idle":"2023-04-23T17:18:27.396412Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Select pairs with only well expressed RNA**","metadata":{}},{"cell_type":"code","source":"RNA_names22 <- my_RNA_names$My_Name[!my_RNA_names$My_Name %in% rem_RNA]\nRNA_names22 <- append(RNA_names22, gsub(\".*_\", \"\", HNRNP_names22))\n\ncat(\"RNA remained:\", length(RNA_names22))\n\nrem_prot <- sub(\"_.*\",\"\", rem_RNA)\nmy_prot_names <- pairs$Protein[!pairs$Protein %in% rem_prot]\ncat(\"\\nProteins remained:\", length(my_prot_names))","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-23T17:18:27.398029Z","iopub.execute_input":"2023-04-23T17:18:27.398997Z","iopub.status.idle":"2023-04-23T17:18:27.417035Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Correlogram**","metadata":{}},{"cell_type":"code","source":"cat(length(my_prot_names), \"Proteins:\\n\")\nmy_prot_names\n\ncat(\"\\n\", length(RNA_names22[!RNA_names22 %in% gsub(\".*_\", \"\", HNRNP_names22)]), \"RNAs:\\n\")\nRNA_names22[!RNA_names22 %in% gsub(\".*_\", \"\", HNRNP_names22)]\n\ncat(\"\\n\", length(HNRNP_names22), \"HNRNPs:\")\ngsub(\".*_\", \"\", HNRNP_names22)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-23T17:18:27.418669Z","iopub.execute_input":"2023-04-23T17:18:27.419601Z","iopub.status.idle":"2023-04-23T17:18:27.464101Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig(40, 55)\nall_plt <- plot_corr_by_fct(HNRNP_set, \"cell_type\", RNA_names22, my_prot_names, show_corr = FALSE, tl.cex = 12)\n\ncat(\"CITEseq 2022, Spearman correlation\")\ndo.call(\"grid.arrange\", c(all_plt, ncol = 2)) #[length(all_plt)]","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-23T17:18:27.466001Z","iopub.execute_input":"2023-04-23T17:18:27.467113Z","iopub.status.idle":"2023-04-23T17:18:36.004449Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Differential correlations with all well expressed proteins and their RNA**","metadata":{}},{"cell_type":"code","source":"dat <- select(HNRNP_set, \"cell_type\", all_of(RNA_names22), all_of(my_prot_names))\nall_dat <- NULL\n#fct = \"MasP\"\nfor(fct in unique(dat[, \"cell_type\"])) {\n\n sbst <- filter(dat, dat[, \"cell_type\"] == fct )\n\n set_corr <- cor(select(sbst, all_of(gsub(\".*_\", \"\", HNRNP_names22))),\n select(sbst, all_of(c(RNA_names22[grepl(\"_RNA\", RNA_names22)]))),\n method = \"spearman\")\n \n set_corr <- cbind(set_corr, \n \"CD45RO_RNA\" = set_corr[, \"CD45_RNA\"],\n \"CD45RA_RNA\" = set_corr[, \"CD45_RNA\"])\n\n dat_cor <- set_corr %>% as.data.frame %>%\n mutate(HNRNP_name = rownames(.),\n cell_type = fct) %>%\n pivot_longer(cols = -c(\"HNRNP_name\", \"cell_type\"), names_to = \"Pair\", values_to = \"cor_with_RNA\") %>%\n mutate(Pair = gsub(\"_.*\", \"\", Pair))\n \n set_corr <- cor(select(sbst, all_of(gsub(\".*_\", \"\", HNRNP_names22))),\n select(sbst, all_of(my_prot_names)),\n method = \"spearman\")\n \n dat_prot_cor <- set_corr %>% as.data.frame %>%\n mutate(HNRNP_name = rownames(.),\n cell_type = fct) %>%\n pivot_longer(cols = -c(\"HNRNP_name\", \"cell_type\"), names_to = \"Pair\", values_to = \"cor_with_prot\") \n \n dat_cor <- merge(dat_cor, dat_prot_cor, by = c(\"Pair\", \"HNRNP_name\", \"cell_type\"))\n\n all_dat <- rbind(all_dat, dat_cor)\n\n}\n\nset_corr <- cor(select(dat, all_of(gsub(\".*_\", \"\", HNRNP_names22))),\n select(dat, all_of(c(RNA_names22[grepl(\"_RNA\", RNA_names22)]))),\n method = \"spearman\")\n\nset_corr <- cbind(set_corr, \n \"CD45RO_RNA\" = set_corr[, \"CD45_RNA\"],\n \"CD45RA_RNA\" = set_corr[, \"CD45_RNA\"])\n\ndat_cor <- set_corr %>% as.data.frame %>%\n mutate(HNRNP_name = rownames(.),\n cell_type = \"All cells\") %>%\n pivot_longer(cols = -c(\"HNRNP_name\", \"cell_type\"), names_to = \"Pair\", values_to = \"cor_with_RNA\") %>%\nmutate(Pair = gsub(\"_.*\", \"\", Pair))\n\nset_corr <- cor(select(dat, all_of(gsub(\".*_\", \"\", HNRNP_names22))),\n select(dat, all_of(my_prot_names)),\n method = \"spearman\")\n\ndat_prot_cor <- set_corr %>% as.data.frame %>%\n mutate(HNRNP_name = rownames(.),\n cell_type = \"All cells\") %>%\n pivot_longer(cols = -c(\"HNRNP_name\", \"cell_type\"), names_to = \"Pair\", values_to = \"cor_with_prot\") \n\ndat_cor <- merge(dat_cor, dat_prot_cor, by = c(\"Pair\", \"HNRNP_name\", \"cell_type\"))\nall_dat <- rbind(all_dat, dat_cor)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-23T17:18:36.007137Z","iopub.execute_input":"2023-04-23T17:18:36.008206Z","iopub.status.idle":"2023-04-23T17:18:40.046175Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"- HNRNP, which correlate positively (>0.1) with RNA and negatively (< -0.1) with protein","metadata":{}},{"cell_type":"code","source":"all_dat %>% \n filter(abs(cor_with_RNA) > 0.1 & abs(cor_with_prot) > 0.1) %>%\n filter(cor_with_RNA > 0 & cor_with_prot < 0) %>%\n mutate(cor = paste0(\"RNA: \", round(cor_with_RNA, 2),\"; Prot: \", round(cor_with_prot, 2))) %>%\n pivot_wider(id_cols = -c(\"cor_with_RNA\", \"cor_with_prot\"),\n names_from = \"HNRNP_name\", values_from = \"cor\")","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-23T17:18:40.048456Z","iopub.execute_input":"2023-04-23T17:18:40.049613Z","iopub.status.idle":"2023-04-23T17:18:40.109743Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"- HNRNP, which correlate negatively (< -0.1) with RNA and positively (>0.1) with protein","metadata":{}},{"cell_type":"code","source":"all_dat %>% \n filter(abs(cor_with_RNA) > 0.1 & abs(cor_with_prot) > 0.1) %>%\n filter(cor_with_RNA < 0 & cor_with_prot > 0) %>%\n mutate(cor = paste0(\"RNA: \", round(cor_with_RNA, 2),\"; Prot: \", round(cor_with_prot, 2))) %>%\n pivot_wider(id_cols = -c(\"cor_with_RNA\", \"cor_with_prot\"),\n names_from = \"HNRNP_name\", values_from = \"cor\")","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-23T17:18:40.111437Z","iopub.execute_input":"2023-04-23T17:18:40.112468Z","iopub.status.idle":"2023-04-23T17:18:40.190202Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#

Splicing genes from GO terms (CD45 isoforms)

","metadata":{}},{"cell_type":"markdown","source":"Found here: https://www.kaggle.com/code/antoninadolgorukova/mmscel-cd45-iso-go-enrichment","metadata":{}},{"cell_type":"markdown","source":"**Subset data from the two datasets**","metadata":{}},{"cell_type":"code","source":"#CITEseq2022\nRNA_pattern <- c(\"PTPRC\", 'SNRPD2', 'POLR2E', 'SNRNP25', 'TXNL4A', 'SNRNP40', 'SNRPE',\n 'POLR2I', 'SNRPD3', 'POLR2F', 'SNRPG', 'SF3B5', 'POLR2L', 'SNU13',\n 'SNRPB', 'SRSF7', 'SRSF2', 'SNRPD1', 'SNRPF', 'YBX1', 'SRRM1',\n 'SNRNP70', 'PPIH', 'LSM3', 'HNRNPF', 'LSM7', 'SNRPA1',\n 'ALYREF', 'STRAP', 'SRSF3', 'LSM5', 'SNRPC', 'DDX39A',\n 'UBL5', 'SRPK1', 'HSPA8', 'C1QBP')\nRNA_names22 <- colnames(smat_RNA)[which(grepl(paste0(RNA_pattern, collapse = \"$|\"),\n colnames(smat_RNA),\n ignore.case = T))] %>% as.character\n\ncat(length(RNA_names22)-1, \"splicing genes and PTPRC:\")\nRNA_names22\n\nset22 <- make_subset(RNA_names22, prot_names)\nRNA_names22 <- gsub(\".*_\",\"\", RNA_names22)\n\ncat(\"CITEseq2022: Head of the data frame with\", ncol(set22)-8, \"RNA/proteins in columns X\", nrow(set22), \"cells in rows\")\nhead(set22)\n\n#CITEseq2023\nRNA_names23 <- colnames(mat_RNA23)[which(grepl(paste0(RNA_pattern, collapse = \"|\"),\n colnames(mat_RNA23),\n ignore.case = T))]\nprot_names23 <- colnames(mat_prot23)[grep(paste(prot_names, collapse = \"|\"), colnames(mat_prot23))]\n\nset23 <- cbind(mat_RNA23[, RNA_names23],\n mat_prot23[, prot_names23]) %>% as.data.frame\ncat(\"CITEseq2023: Head of the data frame with\", ncol(set23), \"RNA/proteins in columns X\", nrow(set23), \"cells in rows\")\nhead(set23)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-23T17:18:40.191957Z","iopub.execute_input":"2023-04-23T17:18:40.192963Z","iopub.status.idle":"2023-04-23T17:18:41.369979Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Check expression levels**","metadata":{}},{"cell_type":"code","source":"fig(25,50)\ngene_set <- c(RNA_names22[RNA_names22 != \"PTPRC\"])\n\nplot_expr_by_fct(set22, gene_set, factor_column = \"sample_size_ct\", ncol = 3)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-23T17:18:41.371718Z","iopub.execute_input":"2023-04-23T17:18:41.372801Z","iopub.status.idle":"2023-04-23T17:18:46.236275Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig(7,7)\n\nsum_expr_RNA23 <- apply(select(set23, -grep(\"CD45\", colnames(set23))), 2, median) %>% sort(decreasing = TRUE)\nsum_expr_RNA23 <- data.frame(\"Median_expression\" = sum_expr_RNA23,\n \"Name\" = names(sum_expr_RNA23))\nsum_expr_RNA23$Name <- factor(sum_expr_RNA23$Name, levels = unique(sum_expr_RNA23$Name ))\n#hline <- sum_expr_RNA23[sum_expr_RNA23$Name == \"PTPRC\",]$Median_expression\n\nggplot(sum_expr_RNA23, aes(x = Name, y = Median_expression)) +\n geom_bar(stat = 'identity', fill=\"#BF9039\", col=\"grey\") + \n# geom_hline(yintercept=sum_expr_RNA23[sum_expr_RNA23$Name == \"PTPRC\",]$Median_expression,\n# linetype=\"dashed\")+\n theme_bw(base_size = 22) +\n xlab(\"\") +\n ggtitle(paste0(\"CITEseq 2023\")) +\n theme(axis.text.x = element_text(angle = 45, hjust=1))\n\ncat(\"Zero median expression:\\n\")\ncat(paste0('\"', sum_expr_RNA23[sum_expr_RNA23$Median_expression == 0,]$Name, '\"', collapse = \", \"))","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-23T17:18:46.238167Z","iopub.execute_input":"2023-04-23T17:18:46.239171Z","iopub.status.idle":"2023-04-23T17:18:46.448093Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Select only well expressed genes**","metadata":{}},{"cell_type":"code","source":"RNA_names23 <- colnames(set23)[!colnames(set23) %in% c(\n \"SRRM1\", \"SNRPE\", \"SRSF3\", \"SF3B5\", \"POLR2L\", \"HSPA8\", \"SRSF2\", \"LSM7\")]\ncat(\"Remained in CITEseq23:\")\nlength(RNA_names23)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-23T17:18:46.449912Z","iopub.execute_input":"2023-04-23T17:18:46.450925Z","iopub.status.idle":"2023-04-23T17:18:46.466989Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### CITEseq22","metadata":{}},{"cell_type":"code","source":"cat(length(prot_names), \"Proteins:\\n\")\ncat(paste0('\"', prot_names, '\"', collapse = \", \"))\ncat(\"\\n\\n\", length(RNA_names22), \"RNAs:\\n\")\ncat(paste0('\"', RNA_names22, '\"', collapse = \", \"))","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-23T17:18:46.468702Z","iopub.execute_input":"2023-04-23T17:18:46.469757Z","iopub.status.idle":"2023-04-23T17:18:46.488371Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig(30, 50)\nall_plt <- plot_corr_by_fct(set22, \"cell_type\", RNA_names22, prot_names, show_corr = FALSE)\n\ncat(\"CITEseq 2022, Spearman correlation\")\ndo.call(\"grid.arrange\", c(all_plt, ncol = 2)) ","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-23T17:18:46.490325Z","iopub.execute_input":"2023-04-23T17:18:46.491336Z","iopub.status.idle":"2023-04-23T17:18:50.701099Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Compare dendrogram trees**","metadata":{}},{"cell_type":"code","source":"dendlst <- dendlist()\nfor(ct in unique(set22$cell_type)) {\n \n sbst <- set22 %>% filter(cell_type == ct) %>%\n select(all_of(RNA_names22), all_of(prot_names)) %>% \n cor(method = \"spearman\") %>% \n dist %>% hclust(\"complete\") %>% as.dendrogram\n dendlst[[ct]] <- sbst\n}\n\ndendlst[[\"All_cells\"]] <- set22 %>%\n select(all_of(RNA_names22), all_of(prot_names)) %>% \n cor(method = \"spearman\") %>% \n dist %>% hclust(\"complete\") %>% as.dendrogram\n\ncors <- cor.dendlist(dendlst)\n# Print correlation matrix\n#round(cors, 2)\nfig(25, 7)\nggcorrplot(cors, hc.order = TRUE, lab = TRUE,\n title = paste0(\"Correlation matrix between dendrograms by cell type\"), lab_size = 6,\n ggtheme = theme_minimal(base_size = 18), tl.cex = 18)\n\n#all.equal(dendlst$MasP, dendlst$All_cells)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-23T17:18:50.702940Z","iopub.execute_input":"2023-04-23T17:18:50.704012Z","iopub.status.idle":"2023-04-23T17:18:52.930844Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Differential correlations with RO and RA isoforms**","metadata":{}},{"cell_type":"code","source":"diff_dat <- diff_corr_by_fct(set22, factor_column = \"cell_type\", RNA_names = RNA_names22, prot_names)\n\ncat(\"Only positive and negative are shown\")\ndiff_dat %>% filter(RA %in% c(\"-\", \"+\") & RO %in% c(\"-\", \"+\")) #filter(cell_type == \"All cells\")","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-23T17:18:52.933802Z","iopub.execute_input":"2023-04-23T17:18:52.935671Z","iopub.status.idle":"2023-04-23T17:18:54.665800Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"graph <- prep_igraph(diff_dat, width_coeff = 3)\n\nvisNetwork(graph$nodes, graph$edges, width=\"100%\", height = 700,\n main=paste0(\"Network of RNA differentially correlated with RO and RA isoforms (targets)\")) %>% \n visOptions(highlightNearest = list(enabled = TRUE, algorithm = \"hierarchical\"), selectedBy = \"label\") %>%\n visHierarchicalLayout() %>%\n visLegend(addEdges = ledges)","metadata":{"execution":{"iopub.status.busy":"2023-04-23T17:18:54.667707Z","iopub.execute_input":"2023-04-23T17:18:54.668720Z","iopub.status.idle":"2023-04-23T17:18:54.843038Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"
\n 💡 Note:
\n The most pronounced opposite correlations - positive with CD45RA and negative with CD45RO in MasP:
\n \n Differential correlation: negative with CD45RO and around zero with CD45RA:
\n\n Differential correlation: negative with CD45RA and around zero with CD45RO:
\n\n Differential correlation: positive with CD45RA and around zero with CD45RO:
\n\nThe strongest correlations of CD53 are with CD45RA in BP and MoP, but this may be due to small sample sizes.\n \n
","metadata":{}},{"cell_type":"markdown","source":"### Compare two CITEseq datasets","metadata":{}},{"cell_type":"code","source":"sbst22 <- set22 %>% select(names(set23))\nset22_corr <- cor(sbst22, method = \"spearman\")\ncat(\"The correlation matrix CITEseq2022\")\nset22_corr\n\nset23_corr <- cor(set23, method = \"spearman\")\ncat(\"\\nThe correlation matrix CITEseq2023\")\nset23_corr","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-23T17:18:54.845010Z","iopub.execute_input":"2023-04-23T17:18:54.846079Z","iopub.status.idle":"2023-04-23T17:18:55.172636Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig(25,12)\nplt1 <- ggcorrplot(set22_corr, hc.order = FALSE, lab = TRUE,\n title = \"CITEseq 2022\", lab_size = 6, p.mat = NULL, #as.matrix(set22_p), insig = \"pch\", pch = 0,\n ggtheme = theme_minimal(base_size = 24), tl.cex = 18)\nplt2 <- ggcorrplot(set23_corr, hc.order = FALSE, lab = TRUE,\n title = \"CITEseq 2023\", lab_size = 6, p.mat = NULL, #set23_p, insig = \"pch\",\n ggtheme = theme_minimal(base_size = 24), tl.cex = 18)\nplt1 + plt2","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-23T17:18:55.174603Z","iopub.execute_input":"2023-04-23T17:18:55.175713Z","iopub.status.idle":"2023-04-23T17:18:55.893066Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Comparison of two dendrograms**","metadata":{}},{"cell_type":"code","source":"fig(25, 7)\n# Compute distance matrix, hierarchical clusterings and create two dendrograms in a list\ndend1 <- set22_corr %>% dist %>% hclust(\"complete\") %>% as.dendrogram\ndend2 <- set23_corr %>% dist %>% hclust(\"complete\") %>% as.dendrogram\ndend_list <- dendlist(dend1, dend2)\n# Align and plot two dendrograms side by side\ndendlist(dend1, dend2) %>%\n untangle(method = \"step1side\") %>% # Find the best alignment layout\n tanglegram(cex_main = 3, lab.cex = 2,\n margin_inner = 15,\n common_subtrees_color_branches = TRUE,\n main_left = \"CITEseq 2022\",\n main_right = \"CITEseq 2023\")\ncor_res <- cor.dendlist(dend_list, method = \"cophenetic\")\ncat(\"Cophenetic correlation between the two dendrograms\", round(cor_res[2,1], 2))\nall.equal(dend1, dend2)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-23T17:18:55.894983Z","iopub.execute_input":"2023-04-23T17:18:55.896048Z","iopub.status.idle":"2023-04-23T17:18:56.237108Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Corretation of the mean of all splicing genes with other proteins in CITEseq 2022","metadata":{"_kg_hide-input":true}},{"cell_type":"code","source":"mean_of_genes <- select(set22, c(\"Row.names\", RNA_names22[RNA_names22 != 'PTPRC']))\nmean_of_genes$Mean_splicing <- apply(mean_of_genes[, -1], 1, mean)\n\nall_prots <- merge(metadata, mat_prot, by=0, all.y = TRUE)\nall_prots$day <- factor(my_RNA$day, levels = c(\"2\",\"3\",\"4\"))\n\nmean_of_genes <- left_join(all_prots, select(mean_of_genes, c(\"Row.names\", \"Mean_splicing\")), by = \"Row.names\")\n\n# add sample sizes per cell type\nsample_size_ct <- dat %>%\n group_by(cell_type) %>%\n summarize(ss_ct = n(),.groups = \"keep\" )\n\nmean_of_genes <- mean_of_genes %>%\n left_join(sample_size_ct, by = c(\"cell_type\")) %>%\n mutate(sample_size_ct = paste0(cell_type, \"\\n\", \"n=\", ss_ct)) %>%\n select(-c(\"ss_ct\"))\n\nfig(25, 12)\n\none_gene_vs_all_prots(mean_of_genes, gene_name = \"Mean_splicing\", factor_column <- \"sample_size_ct\")","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-23T17:18:56.238878Z","iopub.execute_input":"2023-04-23T17:18:56.239879Z","iopub.status.idle":"2023-04-23T17:19:02.551198Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Corretation of the splicing genes with all RNA","metadata":{}},{"cell_type":"code","source":"set22_corr <- all_RNA_corr\nset22_corr <- set22_corr[ , grepl(paste0(\"_\", paste0(RNA_names22[RNA_names22 != 'PTPRC'], collapse = \"$|\"), \"$\"), colnames(set22_corr))] %>% as.matrix\n\ncat(\"The correlation matrix of the\", length(colnames(set22_corr)), \"splicing genes with \", nrow(set22_corr),\"RNA in the CITEseq2022\")\nset22_corr %>% head","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-23T17:19:02.553854Z","iopub.execute_input":"2023-04-23T17:19:02.555522Z","iopub.status.idle":"2023-04-23T17:19:02.908811Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#calculate sum and mean correlation of each row (each of 21k genes)\nmat_rowsums <- rowSums(set22_corr)\nmat_rowmeans <- rowMeans(set22_corr)\n\n#mat_rowsums %>% head\n#mat_rowsums[mat_rowsums!=0] %>% length\n#mat_rowmeans[abs(mat_rowmeans) > 0.2] %>% length\n\ncat(\"Heatmap for all other RNA with absolute mean correlation with the splicing genes > 0.3\")\n\nfig(25, 15)\ntop_corr <- set22_corr[rownames(set22_corr) %in% names(mat_rowmeans[abs(mat_rowmeans) > 0.3]), ]\ntop_corr <- top_corr[!grepl(paste0(RNA_names22[RNA_names22 != 'PTPRC'], collapse = \"$|\"), rownames(top_corr)), ]\n\npheatmap(top_corr, fontsize = 14, fontsize_row = 14,\n color = myColors,\n breaks = breaksList,\n angle_col = 90,\n main = \"\")","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-23T17:19:02.910411Z","iopub.execute_input":"2023-04-23T17:19:02.911375Z","iopub.status.idle":"2023-04-23T17:19:03.355820Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"cutoff <- 0.4\ntop_corr <- set22_corr[rownames(set22_corr) %in% names(mat_rowmeans[abs(mat_rowmeans) > cutoff]), ]\ntop_corr <- top_corr[!grepl(paste0(RNA_names22[RNA_names22 != 'PTPRC'], collapse = \"$|\"), rownames(top_corr)), ]\ncat(\"Mean abs correlation with the 36 splicing genes >\", cutoff, \":\", nrow(top_corr), \"\\n\")\ncat(paste0('\"', gsub(\".*_\", \"\", rownames(top_corr)), '\"', collapse = \", \"))\n\ncutoff <- 0.3\ntop_corr <- set22_corr[rownames(set22_corr) %in% names(mat_rowmeans[mat_rowmeans > cutoff]), ]\ntop_corr <- top_corr[!grepl(paste0(RNA_names22[RNA_names22 != 'PTPRC'], collapse = \"$|\"), rownames(top_corr)), ]\ncat(\"\\n\\nMean correlation with the 36 splicing genes >\", cutoff, \":\", nrow(top_corr), \"\\n\")\ncat(paste0('\"', gsub(\".*_\", \"\", rownames(top_corr)), '\"', collapse = \", \"))\n\ncutoff <- -0.3\ntop_corr <- set22_corr[rownames(set22_corr) %in% names(mat_rowmeans[mat_rowmeans < cutoff]), ]\ntop_corr <- top_corr[!grepl(paste0(RNA_names22[RNA_names22 != 'PTPRC'], collapse = \"$|\"), rownames(top_corr)), ]\ncat(\"\\n\\nMean correlation with the 36 splicing genes <\", cutoff, \":\", nrow(top_corr), \"\\n\")\ncat(paste0('\"', gsub(\".*_\", \"\", rownames(top_corr)), '\"', collapse = \", \"))","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-23T17:19:03.357915Z","iopub.execute_input":"2023-04-23T17:19:03.359125Z","iopub.status.idle":"2023-04-23T17:19:03.400737Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"set22_corr <- as.data.frame(set22_corr)\nset22_corr$RNA <- rownames(set22_corr)\n\n#cat(\"Descriptive statistics of correlations of all\", ncol(set22_corr)-1, \"splicing genes vs All RNA\")\n#set22_corr %>% get_summary_stats(type = \"common\")","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-23T17:19:03.402622Z","iopub.execute_input":"2023-04-23T17:19:03.403609Z","iopub.status.idle":"2023-04-23T17:19:03.426903Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Corretation of splicing genes with protein-RNA pairs in 2022 dataset","metadata":{}},{"cell_type":"markdown","source":"**Make a set with genes of interest and merge with the pairs_set22 (all proteins with an RNA pair form the CITEseq22 dataset)**","metadata":{}},{"cell_type":"code","source":"spl_names22 <- RNA_names22[RNA_names22 != \"PTPRC\"]\nspl_names22 <- colnames(smat_RNA)[which(grepl(paste0(\"_\", spl_names22,\"$\", collapse = \"|\"),\n colnames(smat_RNA),\n ignore.case = T))] %>% as.character #%>% length\nspl_set <- smat_RNA[, spl_names22] %>% as.matrix\ncolnames(spl_set) <- gsub(\".*_\", \"\", colnames(spl_set))\n\nspl_set <- spl_set %>% as.data.frame %>% mutate(\"Row.names\" = rownames(.) )\n\nspl_set <- left_join(pairs_set22, spl_set, by = \"Row.names\")\ncat(\"CITEseq2022: Head of the data frame with\", ncol(spl_set)-6, \"RNA/proteins in columns X\", nrow(spl_set), \"cells in rows\")\nhead(spl_set)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-23T17:19:03.428833Z","iopub.execute_input":"2023-04-23T17:19:03.429882Z","iopub.status.idle":"2023-04-23T17:19:04.012856Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Select pairs with only well expressed RNA**","metadata":{}},{"cell_type":"code","source":"RNA_names22 <- my_RNA_names$My_Name[!my_RNA_names$My_Name %in% rem_RNA]\nRNA_names22 <- append(RNA_names22, gsub(\".*_\", \"\", spl_names22))\n\ncat(\"RNA remained:\", length(RNA_names22))\n\nrem_prot <- sub(\"_.*\",\"\", rem_RNA)\nmy_prot_names <- pairs$Protein[!pairs$Protein %in% rem_prot]\ncat(\"\\nProteins remained:\", length(my_prot_names))","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-23T17:19:04.015697Z","iopub.execute_input":"2023-04-23T17:19:04.016981Z","iopub.status.idle":"2023-04-23T17:19:04.039956Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Correlogram**","metadata":{}},{"cell_type":"code","source":"cat(length(my_prot_names), \"Proteins:\\n\")\nmy_prot_names\n\ncat(\"\\n\", length(RNA_names22[!RNA_names22 %in% gsub(\".*_\", \"\", spl_names22)]), \"RNAs:\\n\")\nRNA_names22[!RNA_names22 %in% gsub(\".*_\", \"\", spl_names22)]\n\ncat(\"\\n\", length(spl_names22), \"Splicing genes:\")\ngsub(\".*_\", \"\", spl_names22)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-23T17:19:04.041769Z","iopub.execute_input":"2023-04-23T17:19:04.042827Z","iopub.status.idle":"2023-04-23T17:19:04.073468Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig(40, 65)\nall_plt <- plot_corr_by_fct(spl_set, \"cell_type\", RNA_names22, my_prot_names, show_corr = FALSE, tl.cex = 12)\n\ncat(\"CITEseq 2022, Spearman correlation\")\ndo.call(\"grid.arrange\", c(all_plt, ncol = 2)) #[length(all_plt)]","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-23T17:19:04.075324Z","iopub.execute_input":"2023-04-23T17:19:04.076425Z","iopub.status.idle":"2023-04-23T17:19:13.630350Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Differential correlations with all well expressed proteins and their RNA**","metadata":{}},{"cell_type":"code","source":"dat <- select(spl_set, \"cell_type\", all_of(RNA_names22), all_of(my_prot_names))\nall_dat <- NULL\n#fct = \"MasP\"\nfor(fct in unique(dat[, \"cell_type\"])) {\n\n sbst <- filter(dat, dat[, \"cell_type\"] == fct )\n\n set_corr <- cor(select(sbst, all_of(gsub(\".*_\", \"\", spl_names22))),\n select(sbst, all_of(c(RNA_names22[grepl(\"_RNA\", RNA_names22)]))),\n method = \"spearman\")\n \n set_corr <- cbind(set_corr, \n \"CD45RO_RNA\" = set_corr[, \"CD45_RNA\"],\n \"CD45RA_RNA\" = set_corr[, \"CD45_RNA\"])\n\n dat_cor <- set_corr %>% as.data.frame %>%\n mutate(spl_name = rownames(.),\n cell_type = fct) %>%\n pivot_longer(cols = -c(\"spl_name\", \"cell_type\"), names_to = \"Pair\", values_to = \"cor_with_RNA\") %>%\n mutate(Pair = gsub(\"_.*\", \"\", Pair))\n \n set_corr <- cor(select(sbst, all_of(gsub(\".*_\", \"\", spl_names22))),\n select(sbst, all_of(my_prot_names)),\n method = \"spearman\")\n \n dat_prot_cor <- set_corr %>% as.data.frame %>%\n mutate(spl_name = rownames(.),\n cell_type = fct) %>%\n pivot_longer(cols = -c(\"spl_name\", \"cell_type\"), names_to = \"Pair\", values_to = \"cor_with_prot\") \n \n dat_cor <- merge(dat_cor, dat_prot_cor, by = c(\"Pair\", \"spl_name\", \"cell_type\"))\n\n all_dat <- rbind(all_dat, dat_cor)\n\n}\n\nset_corr <- cor(select(dat, all_of(gsub(\".*_\", \"\", spl_names22))),\n select(dat, all_of(c(RNA_names22[grepl(\"_RNA\", RNA_names22)]))),\n method = \"spearman\")\n\nset_corr <- cbind(set_corr, \n \"CD45RO_RNA\" = set_corr[, \"CD45_RNA\"],\n \"CD45RA_RNA\" = set_corr[, \"CD45_RNA\"])\n\ndat_cor <- set_corr %>% as.data.frame %>%\n mutate(spl_name = rownames(.),\n cell_type = \"All cells\") %>%\n pivot_longer(cols = -c(\"spl_name\", \"cell_type\"), names_to = \"Pair\", values_to = \"cor_with_RNA\") %>%\nmutate(Pair = gsub(\"_.*\", \"\", Pair))\n\nset_corr <- cor(select(dat, all_of(gsub(\".*_\", \"\", spl_names22))),\n select(dat, all_of(my_prot_names)),\n method = \"spearman\")\n\ndat_prot_cor <- set_corr %>% as.data.frame %>%\n mutate(spl_name = rownames(.),\n cell_type = \"All cells\") %>%\n pivot_longer(cols = -c(\"spl_name\", \"cell_type\"), names_to = \"Pair\", values_to = \"cor_with_prot\") \n\ndat_cor <- merge(dat_cor, dat_prot_cor, by = c(\"Pair\", \"spl_name\", \"cell_type\"))\nall_dat <- rbind(all_dat, dat_cor)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-23T17:19:13.632688Z","iopub.execute_input":"2023-04-23T17:19:13.634308Z","iopub.status.idle":"2023-04-23T17:19:18.905026Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"- Splicing genes, which correlate positively (>0.1) with RNA and negatively (< -0.1) with protein","metadata":{}},{"cell_type":"code","source":"all_dat %>% \n filter(abs(cor_with_RNA) > 0.1 & abs(cor_with_prot) > 0.1) %>%\n filter(cor_with_RNA > 0 & cor_with_prot < 0) %>%\n mutate(cor = paste0(\"RNA: \", round(cor_with_RNA, 2),\"; Prot: \", round(cor_with_prot, 2))) %>%\n pivot_wider(id_cols = -c(\"cor_with_RNA\", \"cor_with_prot\"),\n names_from = \"spl_name\", values_from = \"cor\")","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-23T17:19:18.907216Z","iopub.execute_input":"2023-04-23T17:19:18.908630Z","iopub.status.idle":"2023-04-23T17:19:18.969010Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"- Splicing genes, which correlate negatively (< -0.1) with RNA and positively (>0.1) with protein (CD45RA, BP and MoP are not shown)","metadata":{}},{"cell_type":"code","source":"all_dat %>% \n filter(abs(cor_with_RNA) > 0.1 & abs(cor_with_prot) > 0.1) %>%\n filter(cor_with_RNA < 0 & cor_with_prot > 0) %>%\n mutate(cor = paste0(\"RNA: \", round(cor_with_RNA, 2),\"; Prot: \", round(cor_with_prot, 2))) %>%\n pivot_wider(id_cols = -c(\"cor_with_RNA\", \"cor_with_prot\"),\n names_from = \"spl_name\", values_from = \"cor\") %>%\n filter(!cell_type %in% c(\"BP\", \"MoP\"))","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-23T17:19:18.971041Z","iopub.execute_input":"2023-04-23T17:19:18.992732Z","iopub.status.idle":"2023-04-23T17:19:19.063893Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#

Genes of the ubiquitin-proteasome pathway (UPS) vs CD45, CD45RA, CD45RO

","metadata":{}},{"cell_type":"markdown","source":"**Subset data from the two datasets**","metadata":{}},{"cell_type":"code","source":"RNA_pattern <- c(\"PTPRC$\", 'RPS27A','DUB','OTU','CUE',\n 'UBB', 'UBC', 'UBA52', 'RING','UBA','UCH', 'USP', \n 'JAMM', 'MINDI', 'ULP', 'UBD', 'CUE', 'UMI','Beta-Prp',\n 'Jab_MPN','ULD'\n )\n\nRNA_names22 <- colnames(smat_RNA)[which(grepl(paste0(\"_\", RNA_pattern, collapse = \"|\"),\n colnames(smat_RNA),\n ignore.case = T))] %>% as.character\ncat(length(RNA_names22)-1, \"Found genes:\")\nRNA_names22[RNA_names22 != \"ENSG00000081237_PTPRC\"]\n\nset22 <- make_subset(RNA_names22, prot_names)\nRNA_names22 <- gsub(\".*_\",\"\", RNA_names22)\n\ncat(\"CITEseq2022: Head of the data frame with\", ncol(set22)-8, \"RNA/proteins in columns X\", nrow(set22), \"cells in rows\")\nhead(set22)\n\n#CITEseq2023\nRNA_names23 <- colnames(mat_RNA23)[which(grepl(paste0(RNA_pattern, collapse = \"|\"),\n colnames(mat_RNA23),\n ignore.case = T))]\nprot_names23 <- colnames(mat_prot23)[grep(paste(prot_names23, collapse = \"|\"), colnames(mat_prot23))]\n\nset23 <- cbind(mat_RNA23[, RNA_names23],\n mat_prot23[, prot_names23]) %>% as.data.frame\ncat(\"CITEseq2023: Head of the data frame with\", ncol(set23), \"RNA/proteins in columns X\", nrow(set23), \"cells in rows\")\nhead(set23)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-23T17:19:19.065913Z","iopub.execute_input":"2023-04-23T17:19:19.066956Z","iopub.status.idle":"2023-04-23T17:19:20.335808Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Check expression levels**","metadata":{}},{"cell_type":"code","source":"fig(25,50)\ngene_set <- c(RNA_names22[RNA_names22 != \"PTPRC\"])\n\nplot_expr_by_fct(set22, gene_set, factor_column = \"sample_size_ct\", ncol = 3)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-23T17:19:20.338165Z","iopub.execute_input":"2023-04-23T17:19:20.339234Z","iopub.status.idle":"2023-04-23T17:19:26.110865Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig(7,7)\n\nsum_expr_RNA23 <- apply(select(set23, -grep(\"CD45\", colnames(set23))), 2, median) %>% sort(decreasing = TRUE)\nsum_expr_RNA23 <- data.frame(\"Median_expression\" = sum_expr_RNA23,\n \"Name\" = names(sum_expr_RNA23))\nsum_expr_RNA23$Name <- factor(sum_expr_RNA23$Name, levels = unique(sum_expr_RNA23$Name ))\n#hline <- sum_expr_RNA23[sum_expr_RNA23$Name == \"PTPRC\",]$Median_expression\n\nggplot(sum_expr_RNA23, aes(x = Name, y = Median_expression)) +\n geom_bar(stat = 'identity', fill=\"#BF9039\", col=\"grey\") + \n# geom_hline(yintercept=sum_expr_RNA23[sum_expr_RNA23$Name == \"PTPRC\",]$Median_expression,\n# linetype=\"dashed\")+\n theme_bw(base_size = 22) +\n xlab(\"\") +\n ggtitle(paste0(\"CITEseq 2023\")) +\n theme(axis.text.x = element_text(angle = 45, hjust=1))\n\ncat(\"Zero median expression:\\n\")\ncat(paste0('\"', sum_expr_RNA23[sum_expr_RNA23$Median_expression == 0,]$Name, '\"', collapse = \", \"))","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-23T17:19:26.113311Z","iopub.execute_input":"2023-04-23T17:19:26.114501Z","iopub.status.idle":"2023-04-23T17:19:26.341432Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Select only well expressed genes**","metadata":{}},{"cell_type":"code","source":"RNA_names22 <- RNA_names22[!RNA_names22 %in% c(\n \"CUEDC1\", \"DUBR\", \"OTUB2\", \"OTUD1\", \"OTUD3\", \"OTUD4\", \"OTUD4P1\", \"OTUD5\", \"OTUD6A\",\n \"OTUD6B\", \"OTUD7A\", \"OTUD7B\", \"OTULIN\", \"OTULINL\", \"RING1\", \"RPS27AP1\", \"RPS27AP11\",\n \"RPS27AP12\", \"RPS27AP2\", \"UBA52P5\", \"UBA52P6\", \"UBA6-AS1\", \"UBA7\", \"UBAC2-AS1\", \"UBALD1\",\n \"UBAP1\", \"UBAP1L\", \"UBASH3A\", \"UBASH3B\", \"UBBP4\", \"UCHL1\", \"USP12\", \"USP12-AS2\", \"USP13\",\n \"USP18\", \"USP19\", \"USP2\", \"USP2-AS1\", \"USP20\", \"USP21\", \"USP24\", \"USP25\", \"USP27X\", \"USP27X-AS1\",\n \"USP28\", \"USP3-AS1\", \"USP30\", \"USP30-AS1\", \"USP31\", \"USP32\", \"USP32P1\", \"USP32P2\", \"USP32P3\",\n \"USP35\", \"USP36\", \"USP37\", \"USP38\", \"USP40\", \"USP41\", \"USP42\", \"USP43\", \"USP44\", \"USP45\", \"USP46\",\n \"USP46-AS1\", \"USP49\", \"USP51\", \"USP53\", \"USP54\", \"USP6\", \"USP6NL\", \"USP9Y\", \"USPL1\")]\ncat(length(RNA_names22), \"RNA remained in CITEseq22\")\n\nRNA_names23 <- colnames(set23)[!colnames(set23) %in% c(\n \"DUSP1\", \"TUBA1B\")]\ncat(\"\\n\", length(RNA_names23), \"RNA remained in CITEseq23\")","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-23T17:19:26.343364Z","iopub.execute_input":"2023-04-23T17:19:26.344415Z","iopub.status.idle":"2023-04-23T17:19:26.360438Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### CITEseq22","metadata":{}},{"cell_type":"markdown","source":"**Correlogram**","metadata":{}},{"cell_type":"code","source":"cat(length(prot_names), \"Proteins:\\n\")\ncat(paste0('\"', prot_names, '\"', collapse = \", \"))\ncat(\"\\n\\n\", length(RNA_names22), \"RNAs:\\n\")\ncat(paste0('\"', RNA_names22, '\"', collapse = \", \"))","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-23T17:19:26.362103Z","iopub.execute_input":"2023-04-23T17:19:26.363217Z","iopub.status.idle":"2023-04-23T17:19:26.380273Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig(25, 45)\nall_plt <- plot_corr_by_fct(set22, \"cell_type\", RNA_names22, prot_names, show_corr = FALSE)\n\ncat(\"CITEseq 2022, Spearman correlation\")\ndo.call(\"grid.arrange\", c(all_plt, ncol = 2)) ","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-23T17:19:26.381999Z","iopub.execute_input":"2023-04-23T17:19:26.383135Z","iopub.status.idle":"2023-04-23T17:19:30.402014Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Compare dendrogram trees**","metadata":{}},{"cell_type":"code","source":"dendlst <- dendlist()\nfor(ct in unique(set22$cell_type)) {\n \n sbst <- set22 %>% filter(cell_type == ct) %>%\n select(all_of(RNA_names22), all_of(prot_names)) %>% \n cor(method = \"spearman\") %>% \n dist %>% hclust(\"complete\") %>% as.dendrogram\n dendlst[[ct]] <- sbst\n}\n\ndendlst[[\"All_cells\"]] <- set22 %>%\n select(all_of(RNA_names22), all_of(prot_names)) %>% \n cor(method = \"spearman\") %>% \n dist %>% hclust(\"complete\") %>% as.dendrogram\n\ncors <- cor.dendlist(dendlst)\n# Print correlation matrix\n#round(cors, 2)\nfig(25, 7)\nggcorrplot(cors, hc.order = TRUE, lab = TRUE,\n title = paste0(\"Correlation matrix between dendrograms by cell type\"), lab_size = 6,\n ggtheme = theme_minimal(base_size = 18), tl.cex = 18)\n\n#all.equal(dendlst$MasP, dendlst$All_cells)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-23T17:19:30.404425Z","iopub.execute_input":"2023-04-23T17:19:30.405857Z","iopub.status.idle":"2023-04-23T17:19:32.756969Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Differential correlations with RO and RA isoforms**","metadata":{}},{"cell_type":"code","source":"diff_dat <- diff_corr_by_fct(set22, factor_column = \"cell_type\", RNA_names = RNA_names22, prot_names)\ncat(\"Positive is for R >= 0.1, negative for R <= -0.1 and zero is |R| < 0.1\")\n# cat(\"\\nPositive with RA or RO:\")\n# diff_cor %>% filter(RA %in% c(\"+\"))\n# diff_cor %>% filter(RO %in% c(\"+\"))\ndiff_dat %>% filter(cell_type %in% c(\"MasP\", \"NeuP\"))","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-23T17:19:32.758725Z","iopub.execute_input":"2023-04-23T17:19:32.759820Z","iopub.status.idle":"2023-04-23T17:19:34.519117Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"graph <- prep_igraph(diff_dat, width_coeff = 3)\n\nvisNetwork(graph$nodes, graph$edges, width=\"100%\", height = 700,\n main=paste0(\"Network of RNA differentially correlated with RO and RA isoforms (targets)\")) %>% \n visOptions(highlightNearest = list(enabled = TRUE, algorithm = \"hierarchical\"), selectedBy = \"label\") %>%\n visHierarchicalLayout() %>%\n visLegend(addEdges = ledges)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-23T17:19:34.520921Z","iopub.execute_input":"2023-04-23T17:19:34.521959Z","iopub.status.idle":"2023-04-23T17:19:34.693059Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Comparison of the two datasets","metadata":{}},{"cell_type":"code","source":"#calculate correlations\nsbst22 <- set22 %>% select(intersect(names(set23),names(set22)))\nsbst23 <- set23 %>% select(intersect(names(set23),names(set22)))\n\nset22_corr <- cor(sbst22, method = \"spearman\")\ncat(\"The correlation matrix CITEseq2022\")\nset22_corr\n\nset23_corr <- cor(sbst23, method = \"spearman\")\ncat(\"\\nThe correlation matrix CITEseq2023\")\nset23_corr","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-23T17:19:34.711481Z","iopub.execute_input":"2023-04-23T17:19:34.712872Z","iopub.status.idle":"2023-04-23T17:19:34.891225Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig(25,7)\nplt1 <- ggcorrplot(set22_corr, hc.order = FALSE, lab = TRUE,\n title = \"CITEseq 2022\", lab_size = 6, p.mat = NULL, #as.matrix(set22_p), insig = \"pch\", pch = 0,\n ggtheme = theme_minimal(base_size = 24), tl.cex = 18)\nplt2 <- ggcorrplot(set23_corr, hc.order = FALSE, lab = TRUE,\n title = \"CITEseq 2023\", lab_size = 6, p.mat = NULL, #set23_p, insig = \"pch\",\n ggtheme = theme_minimal(base_size = 24), tl.cex = 18)\nplt1 + plt2","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-23T17:19:34.893050Z","iopub.execute_input":"2023-04-23T17:19:34.894083Z","iopub.status.idle":"2023-04-23T17:19:35.454993Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Comparison of two dendrograms**","metadata":{}},{"cell_type":"code","source":"fig(25, 7)\n# Compute distance matrix, hierarchical clusterings and create two dendrograms in a list\ndend1 <- set22_corr %>% dist %>% hclust(\"complete\") %>% as.dendrogram\ndend2 <- set23_corr %>% dist %>% hclust(\"complete\") %>% as.dendrogram\ndend_list <- dendlist(dend1, dend2)\n# Align and plot two dendrograms side by side\ndendlist(dend1, dend2) %>%\n untangle(method = \"step1side\") %>% # Find the best alignment layout\n tanglegram(cex_main = 3, lab.cex = 2,\n margin_inner = 15,\n common_subtrees_color_branches = TRUE,\n main_left = \"CITEseq 2022\",\n main_right = \"CITEseq 2023\")\ncor_res <- cor.dendlist(dend_list, method = \"cophenetic\")\ncat(\"Cophenetic correlation between the two dendrograms\", round(cor_res[2,1], 2))\nall.equal(dend1, dend2)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-23T17:19:35.457960Z","iopub.execute_input":"2023-04-23T17:19:35.459774Z","iopub.status.idle":"2023-04-23T17:19:35.704433Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Corretation of the RPS27A gene with other proteins in CITEseq 2022","metadata":{}},{"cell_type":"code","source":"RNA_pattern <- c(\"PTPRC$\", \"RPS27A\")\nRNA_names22_set <- colnames(smat_RNA)[which(grepl(paste0(RNA_pattern, collapse = \"|\"),\n colnames(smat_RNA),\n ignore.case = T))] %>% as.character\n\none_vs_all <- make_subset(RNA_names22_set, colnames(mat_prot))\n\nfig(25, 12)\n\none_gene_vs_all_prots(one_vs_all, gene_name = \"RPS27A\", factor_column <- \"sample_size_ct\")","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-23T17:19:35.707397Z","iopub.execute_input":"2023-04-23T17:19:35.709188Z","iopub.status.idle":"2023-04-23T17:19:41.740534Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Corretation of the well-expressed genes of the ubiquitin-proteasome pathway with all other RNAs in CITEseq 2022","metadata":{}},{"cell_type":"code","source":"set22_corr <- all_RNA_corr\nset22_corr <- set22_corr[ , grepl(paste0(\"_\", paste0(RNA_names22[RNA_names22 != 'PTPRC'], collapse = \"$|_\"), \"$\"), colnames(set22_corr))] %>% as.matrix\n\ncat(\"The correlation matrix of the\", length(colnames(set22_corr)), \"UPS genes with \", nrow(set22_corr),\"RNA in the CITEseq2022\")\nset22_corr %>% head","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-23T17:19:41.742348Z","iopub.execute_input":"2023-04-23T17:19:41.743482Z","iopub.status.idle":"2023-04-23T17:19:42.100044Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#calculate sum and mean correlation of each row (each of 21k genes)\nmat_rowsums <- rowSums(set22_corr)\nmat_rowmeans <- rowMeans(set22_corr)\n\n#mat_rowsums %>% head\n#mat_rowsums[mat_rowsums!=0] %>% length\n#mat_rowmeans[abs(mat_rowmeans) > 0.1] %>% length\n\ncat(\"Heatmap for all other RNA with absolute mean correlation with UPS genes > 0.1\")\nfig(25, 15)\ntop_corr <- set22_corr[rownames(set22_corr) %in% names(mat_rowmeans[abs(mat_rowmeans) > 0.1]), ]\ntop_corr <- top_corr[!grepl(paste0(RNA_names22[RNA_names22 != 'PTPRC'], collapse = \"$|\"), rownames(top_corr)), ]\n\npheatmap(top_corr, fontsize = 14,\n color = myColors,\n breaks = breaksList,\n angle_col = 90,\n main = \"\")","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-23T17:19:42.102319Z","iopub.execute_input":"2023-04-23T17:19:42.103627Z","iopub.status.idle":"2023-04-23T17:19:42.523761Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"cutoff <- 0.2\ntop_corr <- set22_corr[rownames(set22_corr) %in% names(mat_rowmeans[abs(mat_rowmeans) > cutoff]), ]\ntop_corr <- top_corr[!grepl(paste0(RNA_names22[RNA_names22 != 'PTPRC'], collapse = \"$|\"), rownames(top_corr)), ]\ncat(\"Mean abs correlation with the 38 UPS genes >\", cutoff, \":\", nrow(top_corr), \"\\n\")\ncat(paste0('\"', gsub(\".*_\", \"\", rownames(top_corr)), '\"', collapse = \", \"))\n\ncutoff <- 0.1\ntop_corr <- set22_corr[rownames(set22_corr) %in% names(mat_rowmeans[mat_rowmeans > cutoff]), ]\ntop_corr <- top_corr[!grepl(paste0(RNA_names22[RNA_names22 != 'PTPRC'], collapse = \"$|\"), rownames(top_corr)), ]\ncat(\"\\n\\nMean correlation with the 38 UPS genes >\", cutoff, \":\", nrow(top_corr), \"\\n\")\ncat(paste0('\"', gsub(\".*_\", \"\", rownames(top_corr)), '\"', collapse = \", \"))\n\ncutoff <- -0.1\ntop_corr <- set22_corr[rownames(set22_corr) %in% names(mat_rowmeans[mat_rowmeans < cutoff]), ]\ntop_corr <- top_corr[!grepl(paste0(RNA_names22[RNA_names22 != 'PTPRC'], collapse = \"$|\"), rownames(top_corr)), ]\ncat(\"\\n\\nMean correlation with the 38 UPS genes <\", cutoff, \":\", nrow(top_corr), \"\\n\")\ncat(paste0('\"', gsub(\".*_\", \"\", rownames(top_corr)), '\"', collapse = \", \"))","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-23T17:19:42.525673Z","iopub.execute_input":"2023-04-23T17:19:42.526738Z","iopub.status.idle":"2023-04-23T17:19:42.568407Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Corretation of the UPS genes with protein-RNA pairs in 2022 dataset","metadata":{}},{"cell_type":"markdown","source":"**Make a set with genes of interest and merge with the pairs_set22 (all proteins with an RNA pair form the CITEseq22 dataset)**","metadata":{}},{"cell_type":"code","source":"UPS_names22 <- RNA_names22[RNA_names22 != 'PTPRC']\nUPS_names22 <- colnames(smat_RNA)[which(grepl(paste0(\"_\", UPS_names22,\"$\", collapse = \"|\"),\n colnames(smat_RNA),\n ignore.case = T))] %>% as.character #%>% length\nUPS_set <- smat_RNA[, UPS_names22] %>% as.matrix\ncolnames(UPS_set) <- gsub(\".*_\", \"\", colnames(UPS_set))\n\nUPS_set <- UPS_set %>% as.data.frame %>% mutate(\"Row.names\" = rownames(.) )\n\nUPS_set <- left_join(pairs_set22, UPS_set, by = \"Row.names\")\ncat(\"CITEseq2022: Head of the data frame with\", ncol(UPS_set)-6, \"RNA/proteins in columns X\", nrow(UPS_set), \"cells in rows\")\nhead(UPS_set)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-23T17:19:42.571315Z","iopub.execute_input":"2023-04-23T17:19:42.572661Z","iopub.status.idle":"2023-04-23T17:19:43.155431Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Select pairs with only well expressed RNA**","metadata":{}},{"cell_type":"code","source":"RNA_names22 <- my_RNA_names$My_Name[!my_RNA_names$My_Name %in% rem_RNA]\nRNA_names22 <- append(RNA_names22, gsub(\".*_\", \"\", UPS_names22))\n\ncat(\"RNA remained:\", length(RNA_names22))\n\nrem_prot <- sub(\"_.*\",\"\", rem_RNA)\nmy_prot_names <- pairs$Protein[!pairs$Protein %in% rem_prot]\ncat(\"\\nProteins remained:\", length(my_prot_names))","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-23T17:19:43.158130Z","iopub.execute_input":"2023-04-23T17:19:43.159306Z","iopub.status.idle":"2023-04-23T17:19:43.199960Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Correlogram**","metadata":{}},{"cell_type":"code","source":"cat(length(my_prot_names), \"Proteins:\\n\")\nmy_prot_names\n\ncat(\"\\n\", length(RNA_names22[!RNA_names22 %in% gsub(\".*_\", \"\", UPS_names22)]), \"RNAs:\\n\")\nRNA_names22[!RNA_names22 %in% gsub(\".*_\", \"\", UPS_names22)]\n\ncat(\"\\n\", length(UPS_names22), \"UPS genes:\")\ngsub(\".*_\", \"\", UPS_names22)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-23T17:19:43.202322Z","iopub.execute_input":"2023-04-23T17:19:43.203790Z","iopub.status.idle":"2023-04-23T17:19:43.235851Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig(40, 65)\nall_plt <- plot_corr_by_fct(UPS_set, \"cell_type\", RNA_names22, my_prot_names, show_corr = FALSE, tl.cex = 12)\n\ncat(\"CITEseq 2022, Spearman correlation\")\ndo.call(\"grid.arrange\", c(all_plt, ncol = 2)) #[length(all_plt)]","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-23T17:19:43.238793Z","iopub.execute_input":"2023-04-23T17:19:43.240147Z","iopub.status.idle":"2023-04-23T17:19:53.492363Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Differential correlations with all well expressed proteins and their RNA**","metadata":{}},{"cell_type":"code","source":"dat <- select(UPS_set, \"cell_type\", all_of(RNA_names22), all_of(my_prot_names))\nall_dat <- NULL\n#fct = \"MasP\"\nfor(fct in unique(dat[, \"cell_type\"])) {\n\n sbst <- filter(dat, dat[, \"cell_type\"] == fct )\n\n set_corr <- cor(select(sbst, all_of(gsub(\".*_\", \"\", UPS_names22))),\n select(sbst, all_of(c(RNA_names22[grepl(\"_RNA\", RNA_names22)]))),\n method = \"spearman\")\n \n set_corr <- cbind(set_corr, \n \"CD45RO_RNA\" = set_corr[, \"CD45_RNA\"],\n \"CD45RA_RNA\" = set_corr[, \"CD45_RNA\"])\n\n dat_cor <- set_corr %>% as.data.frame %>%\n mutate(UPS_name = rownames(.),\n cell_type = fct) %>%\n pivot_longer(cols = -c(\"UPS_name\", \"cell_type\"), names_to = \"Pair\", values_to = \"cor_with_RNA\") %>%\n mutate(Pair = gsub(\"_.*\", \"\", Pair))\n \n set_corr <- cor(select(sbst, all_of(gsub(\".*_\", \"\", UPS_names22))),\n select(sbst, all_of(my_prot_names)),\n method = \"spearman\")\n \n dat_prot_cor <- set_corr %>% as.data.frame %>%\n mutate(UPS_name = rownames(.),\n cell_type = fct) %>%\n pivot_longer(cols = -c(\"UPS_name\", \"cell_type\"), names_to = \"Pair\", values_to = \"cor_with_prot\") \n \n dat_cor <- merge(dat_cor, dat_prot_cor, by = c(\"Pair\", \"UPS_name\", \"cell_type\"))\n\n all_dat <- rbind(all_dat, dat_cor)\n\n}\n\nset_corr <- cor(select(dat, all_of(gsub(\".*_\", \"\", UPS_names22))),\n select(dat, all_of(c(RNA_names22[grepl(\"_RNA\", RNA_names22)]))),\n method = \"spearman\")\n\nset_corr <- cbind(set_corr, \n \"CD45RO_RNA\" = set_corr[, \"CD45_RNA\"],\n \"CD45RA_RNA\" = set_corr[, \"CD45_RNA\"])\n\ndat_cor <- set_corr %>% as.data.frame %>%\n mutate(UPS_name = rownames(.),\n cell_type = \"All cells\") %>%\n pivot_longer(cols = -c(\"UPS_name\", \"cell_type\"), names_to = \"Pair\", values_to = \"cor_with_RNA\") %>%\nmutate(Pair = gsub(\"_.*\", \"\", Pair))\n\nset_corr <- cor(select(dat, all_of(gsub(\".*_\", \"\", UPS_names22))),\n select(dat, all_of(my_prot_names)),\n method = \"spearman\")\n\ndat_prot_cor <- set_corr %>% as.data.frame %>%\n mutate(UPS_name = rownames(.),\n cell_type = \"All cells\") %>%\n pivot_longer(cols = -c(\"UPS_name\", \"cell_type\"), names_to = \"Pair\", values_to = \"cor_with_prot\") \n\ndat_cor <- merge(dat_cor, dat_prot_cor, by = c(\"Pair\", \"UPS_name\", \"cell_type\"))\nall_dat <- rbind(all_dat, dat_cor)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-23T17:19:53.494147Z","iopub.execute_input":"2023-04-23T17:19:53.495178Z","iopub.status.idle":"2023-04-23T17:19:58.797399Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"- UPS genes, which correlate positively (>0.1) with RNA and negatively (< -0.1) with protein","metadata":{}},{"cell_type":"code","source":"all_dat %>% \n filter(abs(cor_with_RNA) > 0.1 & abs(cor_with_prot) > 0.1) %>%\n filter(cor_with_RNA > 0 & cor_with_prot < 0) %>%\n mutate(cor = paste0(\"RNA: \", round(cor_with_RNA, 2),\"; Prot: \", round(cor_with_prot, 2))) %>%\n pivot_wider(id_cols = -c(\"cor_with_RNA\", \"cor_with_prot\"),\n names_from = \"UPS_name\", values_from = \"cor\")","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-23T17:19:58.799248Z","iopub.execute_input":"2023-04-23T17:19:58.800325Z","iopub.status.idle":"2023-04-23T17:19:58.877346Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"- UPS genes, which correlate negatively (< -0.1) with RNA and positively (>0.1) with protein","metadata":{}},{"cell_type":"code","source":"all_dat %>% \n filter(abs(cor_with_RNA) > 0.1 & abs(cor_with_prot) > 0.1) %>%\n filter(cor_with_RNA < 0 & cor_with_prot > 0) %>%\n mutate(cor = paste0(\"RNA: \", round(cor_with_RNA, 2),\"; Prot: \", round(cor_with_prot, 2))) %>%\n pivot_wider(id_cols = -c(\"cor_with_RNA\", \"cor_with_prot\"),\n names_from = \"UPS_name\", values_from = \"cor\")","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-23T17:19:58.879611Z","iopub.execute_input":"2023-04-23T17:19:58.880668Z","iopub.status.idle":"2023-04-23T17:19:58.939009Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#

Ribosomal RNAs

","metadata":{}},{"cell_type":"markdown","source":"Here we examine 4 clusters of ribosomal genes identified [here](http://https://www.kaggle.com/code/antoninadolgorukova/mmscel-cd45related-pca-3d-gif)","metadata":{}},{"cell_type":"markdown","source":"### Cluster 1","metadata":{}},{"cell_type":"markdown","source":"**Subset data from the two datasets**","metadata":{}},{"cell_type":"code","source":"RNA_names22 <- c(\"ENSG00000081237_PTPRC\",\n\"ENSG00000132383_RPA1\", \"ENSG00000165502_RPL36AL\", \n\"ENSG00000163902_RPN1\", \"ENSG00000118705_RPN2\", \n\"ENSG00000141425_RPRD1A\", \"ENSG00000187051_RPS19BP1\",\n\"ENSG00000177189_RPS6KA3\", \"ENSG00000108443_RPS6KB1\") %>% as.character\ncat(length(RNA_names22)-1, \"Found genes:\")\nRNA_names22[RNA_names22 != \"ENSG00000081237_PTPRC\"]\n\nset22 <- make_subset(RNA_names22, prot_names)\nRNA_names22 <- gsub(\".*_\",\"\", RNA_names22)\n\ncat(\"CITEseq2022: Head of the data frame with\", ncol(set22)-8, \"RNA/proteins in columns X\", nrow(set22), \"cells in rows\")\nhead(set22)\n\n#CITEseq2023\nRNA_names23 <- colnames(mat_RNA23)[which(grepl(paste0(gsub(\"_.*\", \"\", RNA_names22), collapse = \"|\"),\n colnames(mat_RNA23),\n ignore.case = T))]\nprot_names23 <- colnames(mat_prot23)[grep(paste(prot_names, collapse = \"|\"), colnames(mat_prot23))]\n\nset23 <- cbind(mat_RNA23[, RNA_names23],\n mat_prot23[, prot_names23]) %>% as.data.frame\ncat(\"\\nCITEseq2023: Head of the data frame with\", ncol(set23), \"RNA/proteins in columns X\", nrow(set23), \"cells in rows\")\nhead(set23)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-23T17:19:58.940849Z","iopub.execute_input":"2023-04-23T17:19:58.941870Z","iopub.status.idle":"2023-04-23T17:19:59.892046Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Check expression levels**","metadata":{}},{"cell_type":"code","source":"fig(25,12)\ngene_set <- c(RNA_names22[!RNA_names22 == \"PTPRC\"])\n\nplot_expr_by_fct(set22, gene_set, factor_column = \"sample_size_ct\", ncol = 3)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-23T17:19:59.893967Z","iopub.execute_input":"2023-04-23T17:19:59.895052Z","iopub.status.idle":"2023-04-23T17:20:00.900144Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig(7,7)\n\nsum_expr_RNA23 <- apply(select(set23, -grep(\"CD45\", colnames(set23))), 2, median) %>% sort(decreasing = TRUE)\nsum_expr_RNA23 <- data.frame(\"Median_expression\" = sum_expr_RNA23,\n \"Name\" = names(sum_expr_RNA23))\nsum_expr_RNA23$Name <- factor(sum_expr_RNA23$Name, levels = unique(sum_expr_RNA23$Name ))\n#hline <- sum_expr_RNA23[sum_expr_RNA23$Name == \"PTPRC\",]$Median_expression\n\nggplot(sum_expr_RNA23, aes(x = Name, y = Median_expression)) +\n geom_bar(stat = 'identity', fill=\"#BF9039\", col=\"grey\") + \n# geom_hline(yintercept=sum_expr_RNA23[sum_expr_RNA23$Name == \"PTPRC\",]$Median_expression,\n# linetype=\"dashed\")+\n theme_bw(base_size = 22) +\n xlab(\"\") +\n ggtitle(paste0(\"CITEseq 2023\")) +\n theme(axis.text.x = element_text(angle = 45, hjust=1))\n\ncat(\"Zero median expression:\\n\")\ncat(paste0('\"', sum_expr_RNA23[sum_expr_RNA23$Median_expression == 0,]$Name, '\"', collapse = \", \"))","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-23T17:20:00.902071Z","iopub.execute_input":"2023-04-23T17:20:00.903146Z","iopub.status.idle":"2023-04-23T17:20:01.120629Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### CITEseq22","metadata":{}},{"cell_type":"code","source":"cat(length(prot_names), \"Proteins:\\n\")\ncat(paste0('\"', prot_names, '\"', collapse = \", \"))\ncat(\"\\n\\n\", length(RNA_names22), \"RNAs:\\n\")\ncat(paste0('\"', RNA_names22, '\"', collapse = \", \"))","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-23T17:20:01.122444Z","iopub.execute_input":"2023-04-23T17:20:01.123475Z","iopub.status.idle":"2023-04-23T17:20:01.143569Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig(25, 12)\nall_plt <- plot_corr_by_fct(set22, \"cell_type\", RNA_names22, prot_names, show_corr = FALSE)\n\ncat(\"CITEseq 2022, Spearman correlation\")\ndo.call(\"grid.arrange\", c(all_plt, ncol = 4)) ","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-23T17:20:01.146365Z","iopub.execute_input":"2023-04-23T17:20:01.147492Z","iopub.status.idle":"2023-04-23T17:20:03.017391Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Compare dendrogram trees**","metadata":{}},{"cell_type":"code","source":"dendlst <- dendlist()\nfor(ct in unique(set22$cell_type)) {\n \n sbst <- set22 %>% filter(cell_type == ct) %>%\n select(all_of(RNA_names22), all_of(prot_names)) %>% \n cor(method = \"spearman\") %>% \n dist %>% hclust(\"complete\") %>% as.dendrogram\n dendlst[[ct]] <- sbst\n}\n\ndendlst[[\"All_cells\"]] <- set22 %>%\n select(all_of(RNA_names22), all_of(prot_names)) %>% \n cor(method = \"spearman\") %>% \n dist %>% hclust(\"complete\") %>% as.dendrogram\n\ncors <- cor.dendlist(dendlst)\n# Print correlation matrix\n#round(cors, 2)\nfig(25, 7)\nggcorrplot(cors, hc.order = TRUE, lab = TRUE,\n title = paste0(\"Correlation matrix between dendrograms by cell type\"), lab_size = 6,\n ggtheme = theme_minimal(base_size = 18), tl.cex = 18)\n\n#all.equal(dendlst$MasP, dendlst$All_cells)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-23T17:20:03.019770Z","iopub.execute_input":"2023-04-23T17:20:03.021199Z","iopub.status.idle":"2023-04-23T17:20:03.906120Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Differential correlations with RO and RA isoforms**","metadata":{}},{"cell_type":"code","source":"diff_cor <- diff_corr_by_fct(set22, factor_column = \"cell_type\", RNA_names = RNA_names22, prot_names)\ncat(\"Positive is for R >= 0.1, negative for R <= -0.1 and zero is |R| < 0.1\")\ndiff_cor# %>% filter(RA == \"-\")","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-23T17:20:03.907744Z","iopub.execute_input":"2023-04-23T17:20:03.908665Z","iopub.status.idle":"2023-04-23T17:20:04.509088Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Comparison of the two datasets","metadata":{}},{"cell_type":"code","source":"#calculate correlations\nsbst22 <- set22 %>% select(intersect(names(set23),names(set22)))\nsbst23 <- set23 %>% select(intersect(names(set23),names(set22)))\n\nset22_corr <- cor(sbst22, method = \"spearman\")\ncat(\"The correlation matrix CITEseq2022\")\nset22_corr\n\nset23_corr <- cor(sbst23, method = \"spearman\")\ncat(\"\\nThe correlation matrix CITEseq2023\")\nset23_corr","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-23T17:20:04.510765Z","iopub.execute_input":"2023-04-23T17:20:04.511739Z","iopub.status.idle":"2023-04-23T17:20:04.624864Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig(25,5)\nplt1 <- ggcorrplot(set22_corr, hc.order = FALSE, lab = TRUE,\n title = \"CITEseq 2022\", lab_size = 6, p.mat = NULL, #as.matrix(set22_p), insig = \"pch\", pch = 0,\n ggtheme = theme_minimal(base_size = 24), tl.cex = 18)\nplt2 <- ggcorrplot(set23_corr, hc.order = FALSE, lab = TRUE,\n title = \"CITEseq 2023\", lab_size = 6, p.mat = NULL, #set23_p, insig = \"pch\",\n ggtheme = theme_minimal(base_size = 24), tl.cex = 18)\nplt1 + plt2","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-23T17:20:04.626681Z","iopub.execute_input":"2023-04-23T17:20:04.627651Z","iopub.status.idle":"2023-04-23T17:20:05.125329Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Comparison of two dendrograms**","metadata":{}},{"cell_type":"code","source":"fig(25, 7)\n# Compute distance matrix, hierarchical clusterings and create two dendrograms in a list\ndend1 <- set22_corr %>% dist %>% hclust(\"complete\") %>% as.dendrogram\ndend2 <- set23_corr %>% dist %>% hclust(\"complete\") %>% as.dendrogram\ndend_list <- dendlist(dend1, dend2)\n# Align and plot two dendrograms side by side\ndendlist(dend1, dend2) %>%\n untangle(method = \"step1side\") %>% # Find the best alignment layout\n tanglegram(cex_main = 3, lab.cex = 2,\n margin_inner = 15,\n common_subtrees_color_branches = TRUE,\n main_left = \"CITEseq 2022\",\n main_right = \"CITEseq 2023\")\ncor_res <- cor.dendlist(dend_list, method = \"cophenetic\")\ncat(\"Cophenetic correlation between the two dendrograms\", round(cor_res[2,1], 2))\nall.equal(dend1, dend2)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-23T17:20:05.127029Z","iopub.execute_input":"2023-04-23T17:20:05.127965Z","iopub.status.idle":"2023-04-23T17:20:05.332794Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Cluster 2","metadata":{}},{"cell_type":"markdown","source":"**Subset data from the two datasets**\n","metadata":{}},{"cell_type":"code","source":"RNA_names22 <- c(\"ENSG00000081237_PTPRC\",\n\"ENSG00000117748_RPA2\", \"ENSG00000106399_RPA3\" , \n\"ENSG00000197498_RPF2\" , \"ENSG00000185834_RPL12P4\", \n\"ENSG00000265681_RPL17\", \"ENSG00000213442_RPL18AP3\",\n\"ENSG00000163584_RPL22L1\", \"ENSG00000037241_RPL26L1\", \n\"ENSG00000236264_RPL26P30\", \"ENSG00000131469_RPL27\", \n\"ENSG00000163923_RPL39L\", \"ENSG00000232573_RPL3P4\", \n\"ENSG00000256338_RPL41P2\", \"ENSG00000237550_RPL9P9\", \n\"ENSG00000148688_RPP30\", \"ENSG00000197728_RPS26\" , \n\"ENSG00000224631_RPS27AP16\", \"ENSG00000185088_RPS27L\", \n\"ENSG00000129824_RPS4Y1\", \"ENSG00000162302_RPS6KA4\", \n\"ENSG00000175634_RPS6KB2\", \"ENSG00000007376_RPUSD1\" , \n\"ENSG00000166133_RPUSD2\", \"ENSG00000156990_RPUSD3\") %>% as.character\n\ncat(length(RNA_names22)-1, \"Found genes:\")\nRNA_names22[RNA_names22 != \"ENSG00000081237_PTPRC\"]\n\nset22 <- make_subset(RNA_names22, prot_names)\nRNA_names22 <- gsub(\".*_\",\"\", RNA_names22)\n\ncat(\"CITEseq2022: Head of the data frame with\", ncol(set22)-8, \"RNA/proteins in columns X\", nrow(set22), \"cells in rows\")\nhead(set22)\n\n#CITEseq2023\nRNA_names23 <- colnames(mat_RNA23)[which(grepl(paste0(gsub(\"_.*\", \"\", RNA_names22), collapse = \"$|\"),\n colnames(mat_RNA23),\n ignore.case = T))]\nprot_names23 <- colnames(mat_prot23)[grep(paste(prot_names, collapse = \"|\"), colnames(mat_prot23))]\n\nset23 <- cbind(mat_RNA23[, RNA_names23],\n mat_prot23[, prot_names23]) %>% as.data.frame\ncat(\"\\nCITEseq2023: Head of the data frame with\", ncol(set23), \"RNA/proteins in columns X\", nrow(set23), \"cells in rows\")\nhead(set23)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-23T17:20:05.335648Z","iopub.execute_input":"2023-04-23T17:20:05.336909Z","iopub.status.idle":"2023-04-23T17:20:06.345780Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Check expression levels**","metadata":{}},{"cell_type":"code","source":"fig(25,25)\ngene_set <- c(RNA_names22[!RNA_names22 == \"PTPRC\"])\n\nplot_expr_by_fct(set22, gene_set, factor_column = \"sample_size_ct\", ncol = 3)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-23T17:20:06.347319Z","iopub.execute_input":"2023-04-23T17:20:06.348234Z","iopub.status.idle":"2023-04-23T17:20:08.652999Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig(7,7)\n\nsum_expr_RNA23 <- apply(select(set23, -grep(\"CD45\", colnames(set23))), 2, median) %>% sort(decreasing = TRUE)\nsum_expr_RNA23 <- data.frame(\"Median_expression\" = sum_expr_RNA23,\n \"Name\" = names(sum_expr_RNA23))\nsum_expr_RNA23$Name <- factor(sum_expr_RNA23$Name, levels = unique(sum_expr_RNA23$Name ))\n#hline <- sum_expr_RNA23[sum_expr_RNA23$Name == \"PTPRC\",]$Median_expression\n\nggplot(sum_expr_RNA23, aes(x = Name, y = Median_expression)) +\n geom_bar(stat = 'identity', fill=\"#BF9039\", col=\"grey\") + \n# geom_hline(yintercept=sum_expr_RNA23[sum_expr_RNA23$Name == \"PTPRC\",]$Median_expression,\n# linetype=\"dashed\")+\n theme_bw(base_size = 22) +\n xlab(\"\") +\n ggtitle(paste0(\"CITEseq 2023\")) +\n theme(axis.text.x = element_text(angle = 45, hjust=1))\n\ncat(\"Zero median expression:\\n\")\ncat(paste0('\"', sum_expr_RNA23[sum_expr_RNA23$Median_expression == 0,]$Name, '\"', collapse = \", \"))","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-23T17:20:08.654807Z","iopub.execute_input":"2023-04-23T17:20:08.655910Z","iopub.status.idle":"2023-04-23T17:20:08.866813Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### CITEseq22","metadata":{}},{"cell_type":"code","source":"cat(length(prot_names), \"Proteins:\\n\")\ncat(paste0('\"', prot_names, '\"', collapse = \", \"))\ncat(\"\\n\\n\", length(RNA_names22), \"RNAs:\\n\")\ncat(paste0('\"', RNA_names22, '\"', collapse = \", \"))","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-23T17:20:08.868427Z","iopub.execute_input":"2023-04-23T17:20:08.869360Z","iopub.status.idle":"2023-04-23T17:20:08.886069Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig(25, 35)\nall_plt <- plot_corr_by_fct(set22, \"cell_type\", RNA_names22, prot_names, show_corr = FALSE, tl.cex = 16)\n\ncat(\"CITEseq 2022, Spearman correlation\")\ndo.call(\"grid.arrange\", c(all_plt, ncol = 2)) ","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-23T17:20:08.887675Z","iopub.execute_input":"2023-04-23T17:20:08.888579Z","iopub.status.idle":"2023-04-23T17:20:11.910462Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Compare dendrogram trees**","metadata":{}},{"cell_type":"code","source":"dendlst <- dendlist()\nfor(ct in unique(set22$cell_type)) {\n \n sbst <- set22 %>% filter(cell_type == ct) %>%\n select(all_of(RNA_names22), all_of(prot_names)) %>% \n cor(method = \"spearman\") %>% \n dist %>% hclust(\"complete\") %>% as.dendrogram\n dendlst[[ct]] <- sbst\n}\n\ndendlst[[\"All_cells\"]] <- set22 %>%\n select(all_of(RNA_names22), all_of(prot_names)) %>% \n cor(method = \"spearman\") %>% \n dist %>% hclust(\"complete\") %>% as.dendrogram\n\ncors <- cor.dendlist(dendlst)\n# Print correlation matrix\n#round(cors, 2)\nfig(25, 7)\nggcorrplot(cors, hc.order = TRUE, lab = TRUE,\n title = paste0(\"Correlation matrix between dendrograms by cell type\"), lab_size = 6,\n ggtheme = theme_minimal(base_size = 18), tl.cex = 18)\n\n#all.equal(dendlst$MasP, dendlst$All_cells)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-23T17:20:11.912862Z","iopub.execute_input":"2023-04-23T17:20:11.914175Z","iopub.status.idle":"2023-04-23T17:20:13.483407Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Differential correlations with RO and RA isoforms**","metadata":{}},{"cell_type":"code","source":"diff_cor <- diff_corr_by_fct(set22, factor_column = \"cell_type\", RNA_names = RNA_names22, prot_names)\ncat(\"Positive is for R >= 0.1, negative for R <= -0.1 and zero is |R| < 0.1\")\ndiff_cor %>% filter(cell_type %in% c(\"MasP\", \"NeuP\", \"All cells\"))","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-23T17:20:13.485921Z","iopub.execute_input":"2023-04-23T17:20:13.487429Z","iopub.status.idle":"2023-04-23T17:20:14.668312Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Comparison of the two datasets","metadata":{}},{"cell_type":"code","source":"#calculate correlations\nsbst22 <- set22 %>% select(intersect(names(set23),names(set22)))\nsbst23 <- set23 %>% select(intersect(names(set23),names(set22)))\n\nset22_corr <- cor(sbst22, method = \"spearman\")\ncat(\"The correlation matrix CITEseq2022\")\nset22_corr\n\nset23_corr <- cor(sbst23, method = \"spearman\")\ncat(\"\\nThe correlation matrix CITEseq2023\")\nset23_corr","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-23T17:20:14.670007Z","iopub.execute_input":"2023-04-23T17:20:14.670978Z","iopub.status.idle":"2023-04-23T17:20:14.847002Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig(25,7)\nplt1 <- ggcorrplot(set22_corr, hc.order = FALSE, lab = TRUE,\n title = \"CITEseq 2022\", lab_size = 6, p.mat = NULL, #as.matrix(set22_p), insig = \"pch\", pch = 0,\n ggtheme = theme_minimal(base_size = 24), tl.cex = 18)\nplt2 <- ggcorrplot(set23_corr, hc.order = FALSE, lab = TRUE,\n title = \"CITEseq 2023\", lab_size = 6, p.mat = NULL, #set23_p, insig = \"pch\",\n ggtheme = theme_minimal(base_size = 24), tl.cex = 18)\nplt1 + plt2","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-23T17:20:14.848593Z","iopub.execute_input":"2023-04-23T17:20:14.849503Z","iopub.status.idle":"2023-04-23T17:20:15.399633Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Comparison of two dendrograms**","metadata":{}},{"cell_type":"code","source":"fig(25, 7)\n# Compute distance matrix, hierarchical clusterings and create two dendrograms in a list\ndend1 <- set22_corr %>% dist %>% hclust(\"complete\") %>% as.dendrogram\ndend2 <- set23_corr %>% dist %>% hclust(\"complete\") %>% as.dendrogram\ndend_list <- dendlist(dend1, dend2)\n# Align and plot two dendrograms side by side\ndendlist(dend1, dend2) %>%\n untangle(method = \"step1side\") %>% # Find the best alignment layout\n tanglegram(cex_main = 3, lab.cex = 2,\n margin_inner = 15,\n common_subtrees_color_branches = TRUE,\n main_left = \"CITEseq 2022\",\n main_right = \"CITEseq 2023\")\ncor_res <- cor.dendlist(dend_list, method = \"cophenetic\")\ncat(\"Cophenetic correlation between the two dendrograms\", round(cor_res[2,1], 2))\nall.equal(dend1, dend2)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-23T17:20:15.401856Z","iopub.execute_input":"2023-04-23T17:20:15.403026Z","iopub.status.idle":"2023-04-23T17:20:15.642413Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Cluster 3","metadata":{}},{"cell_type":"markdown","source":"**Subset data from the two datasets**","metadata":{}},{"cell_type":"code","source":"RNA_names22 <- c(\"ENSG00000081237_PTPRC\",\n\"ENSG00000147403_RPL10\", \"ENSG00000198755_RPL10A\", \"ENSG00000142676_RPL11\" ,\n\"ENSG00000197958_RPL12\" , \"ENSG00000167526_RPL13\" ,\"ENSG00000188846_RPL14\" ,\n\"ENSG00000174748_RPL15\", \"ENSG00000063177_RPL18\" ,\"ENSG00000105640_RPL18A\",\n\"ENSG00000108298_RPL19\" , \"ENSG00000122026_RPL21\" ,\"ENSG00000116251_RPL22\" ,\n\"ENSG00000198242_RPL23A\", \"ENSG00000114391_RPL24\" ,\"ENSG00000161970_RPL26\" ,\n\"ENSG00000108107_RPL28\", \"ENSG00000162244_RPL29\" , \"ENSG00000100316_RPL3\" , \n\"ENSG00000156482_RPL30\", \"ENSG00000144713_RPL32\" , \"ENSG00000109475_RPL34\" ,\n\"ENSG00000136942_RPL35\" , \"ENSG00000182899_RPL35A\" ,\"ENSG00000130255_RPL36\" ,\n\"ENSG00000241343_RPL36A\", \"ENSG00000145592_RPL37\" ,\"ENSG00000197756_RPL37A\",\n\"ENSG00000198918_RPL39\", \"ENSG00000229117_RPL41\" , \"ENSG00000122406_RPL5\" , \n\"ENSG00000089009_RPL6\" , \"ENSG00000148303_RPL7A\", \"ENSG00000161016_RPL8\" , \n\"ENSG00000163682_RPL9\" , \"ENSG00000089157_RPLP0\" , \"ENSG00000137818_RPLP1\" ,\n\"ENSG00000177600_RPLP2\" , \"ENSG00000124614_RPS10\" ,\"ENSG00000112306_RPS12\",\n\"ENSG00000110700_RPS13\", \"ENSG00000164587_RPS14\" , \"ENSG00000115268_RPS15\",\n\"ENSG00000134419_RPS15A\", \"ENSG00000105193_RPS16\" ,\"ENSG00000231500_RPS18\", \n\"ENSG00000105372_RPS19\", \"ENSG00000140988_RPS2\", \"ENSG00000171858_RPS21\", \n\"ENSG00000186468_RPS23\", \"ENSG00000138326_RPS24\" , \"ENSG00000118181_RPS25\", \n\"ENSG00000177954_RPS27\", \"ENSG00000143947_RPS27A\", \"ENSG00000233927_RPS28\", \n\"ENSG00000213741_RPS29\", \"ENSG00000149273_RPS3\", \"ENSG00000145425_RPS3A\",\n\"ENSG00000198034_RPS4X\", \"ENSG00000083845_RPS5\" , \"ENSG00000137154_RPS6\", \n\"ENSG00000171863_RPS7\", \"ENSG00000142937_RPS8\", \"ENSG00000170889_RPS9\", \n\"ENSG00000168028_RPSA\" ) %>% as.character\n\ncat(length(RNA_names22)-1, \"Found genes:\")\nRNA_names22[RNA_names22 != \"ENSG00000081237_PTPRC\"]\n\nset22 <- make_subset(RNA_names22, prot_names)\nRNA_names22 <- gsub(\".*_\",\"\", RNA_names22)\n\ncat(\"CITEseq2022: Head of the data frame with\", ncol(set22)-8, \"RNA/proteins in columns X\", nrow(set22), \"cells in rows\")\nhead(set22)\n\n#CITEseq2023\nRNA_names23 <- colnames(mat_RNA23)[which(grepl(paste0(gsub(\"_.*\", \"\", RNA_names22), collapse = \"$|\"),\n colnames(mat_RNA23),\n ignore.case = T))]\nprot_names23 <- colnames(mat_prot23)[grep(paste(prot_names, collapse = \"|\"), colnames(mat_prot23))]\n\nset23 <- cbind(mat_RNA23[, RNA_names23],\n mat_prot23[, prot_names23]) %>% as.data.frame\ncat(\"\\nCITEseq2023: Head of the data frame with\", ncol(set23), \"RNA/proteins in columns X\", nrow(set23), \"cells in rows\")\nhead(set23)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-23T17:20:15.645307Z","iopub.execute_input":"2023-04-23T17:20:15.647214Z","iopub.status.idle":"2023-04-23T17:20:16.997824Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Check expression levels**","metadata":{}},{"cell_type":"code","source":"fig(25,50)\ngene_set <- c(RNA_names22[!RNA_names22 == \"PTPRC\"])\n\nplot_expr_by_fct(set22, gene_set, factor_column = \"sample_size_ct\", ncol = 3)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-23T17:20:16.999429Z","iopub.execute_input":"2023-04-23T17:20:17.000435Z","iopub.status.idle":"2023-04-23T17:20:23.832899Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig(25,7)\n\nsum_expr_RNA23 <- apply(select(set23, -grep(\"CD45\", colnames(set23))), 2, median) %>% sort(decreasing = TRUE)\nsum_expr_RNA23 <- data.frame(\"Median_expression\" = sum_expr_RNA23,\n \"Name\" = names(sum_expr_RNA23))\nsum_expr_RNA23$Name <- factor(sum_expr_RNA23$Name, levels = unique(sum_expr_RNA23$Name ))\n#hline <- sum_expr_RNA23[sum_expr_RNA23$Name == \"PTPRC\",]$Median_expression\n\nggplot(sum_expr_RNA23, aes(x = Name, y = Median_expression)) +\n geom_bar(stat = 'identity', fill=\"#BF9039\", col=\"grey\") + \n# geom_hline(yintercept=sum_expr_RNA23[sum_expr_RNA23$Name == \"PTPRC\",]$Median_expression,\n# linetype=\"dashed\")+\n theme_bw(base_size = 22) +\n xlab(\"\") +\n ggtitle(paste0(\"CITEseq 2023\")) +\n theme(axis.text.x = element_text(angle = 45, hjust=1))\n\ncat(\"Zero median expression:\\n\")\ncat(paste0('\"', sum_expr_RNA23[sum_expr_RNA23$Median_expression == 0,]$Name, '\"', collapse = \", \"))","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-23T17:20:23.835505Z","iopub.execute_input":"2023-04-23T17:20:23.837015Z","iopub.status.idle":"2023-04-23T17:20:24.202219Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### CITEseq22","metadata":{}},{"cell_type":"code","source":"cat(length(prot_names), \"Proteins:\\n\")\ncat(paste0('\"', prot_names, '\"', collapse = \", \"))\ncat(\"\\n\\n\", length(RNA_names22), \"RNAs:\\n\")\ncat(paste0('\"', RNA_names22, '\"', collapse = \", \"))","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-23T17:20:24.203993Z","iopub.execute_input":"2023-04-23T17:20:24.205067Z","iopub.status.idle":"2023-04-23T17:20:24.223021Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig(30, 50)\nall_plt <- plot_corr_by_fct(set22, \"cell_type\", RNA_names22, prot_names, show_corr = FALSE, tl.cex = 14)\n\ncat(\"CITEseq 2022, Spearman correlation\")\ndo.call(\"grid.arrange\", c(all_plt, ncol = 2)) ","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-23T17:20:24.224683Z","iopub.execute_input":"2023-04-23T17:20:24.225659Z","iopub.status.idle":"2023-04-23T17:20:30.479255Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Compare dendrogram trees**","metadata":{}},{"cell_type":"code","source":"dendlst <- dendlist()\nfor(ct in unique(set22$cell_type)) {\n \n sbst <- set22 %>% filter(cell_type == ct) %>%\n select(all_of(RNA_names22), all_of(prot_names)) %>% \n cor(method = \"spearman\") %>% \n dist %>% hclust(\"complete\") %>% as.dendrogram\n dendlst[[ct]] <- sbst\n}\n\ndendlst[[\"All_cells\"]] <- set22 %>%\n select(all_of(RNA_names22), all_of(prot_names)) %>% \n cor(method = \"spearman\") %>% \n dist %>% hclust(\"complete\") %>% as.dendrogram\n\ncors <- cor.dendlist(dendlst)\n# Print correlation matrix\n#round(cors, 2)\nfig(25, 7)\nggcorrplot(cors, hc.order = TRUE, lab = TRUE,\n title = paste0(\"Correlation matrix between dendrograms by cell type\"), lab_size = 6,\n ggtheme = theme_minimal(base_size = 18), tl.cex = 18)\n\n#all.equal(dendlst$MasP, dendlst$All_cells)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-23T17:20:30.481938Z","iopub.execute_input":"2023-04-23T17:20:30.483498Z","iopub.status.idle":"2023-04-23T17:20:34.408477Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Differential correlations with RO and RA isoforms**","metadata":{}},{"cell_type":"code","source":"diff_cor <- diff_corr_by_fct(set22, factor_column = \"cell_type\", RNA_names = RNA_names22, prot_names)\ncat(\"Positive is for R >= 0.1, negative for R <= -0.1 and zero is |R| < 0.1\")\ndiff_cor %>% filter(cell_type %in% c(\"All cells\")) #\"MasP\", \"NeuP\", ","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-23T17:20:34.410332Z","iopub.execute_input":"2023-04-23T17:20:34.411337Z","iopub.status.idle":"2023-04-23T17:20:37.558797Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Comparison of the two datasets","metadata":{}},{"cell_type":"code","source":"#calculate correlations\nsbst22 <- set22 %>% select(intersect(names(set23),names(set22)))\nsbst23 <- set23 %>% select(intersect(names(set23),names(set22)))\n\nset22_corr <- cor(sbst22, method = \"spearman\")\ncat(\"The correlation matrix CITEseq2022\")\nset22_corr %>% head\n\nset23_corr <- cor(sbst23, method = \"spearman\")\ncat(\"\\nThe correlation matrix CITEseq2023\")\nset23_corr %>% head","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-23T17:20:37.561283Z","iopub.execute_input":"2023-04-23T17:20:37.562709Z","iopub.status.idle":"2023-04-23T17:20:39.247395Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig(25,12)\nplt1 <- ggcorrplot(set22_corr, hc.order = FALSE, lab = FALSE,\n title = \"CITEseq 2022\", lab_size = 3, p.mat = NULL, #as.matrix(set22_p), insig = \"pch\", pch = 0,\n ggtheme = theme_minimal(base_size = 24), tl.cex = 14)\nplt2 <- ggcorrplot(set23_corr, hc.order = FALSE, lab = FALSE,\n title = \"CITEseq 2023\", lab_size = 3, p.mat = NULL, #set23_p, insig = \"pch\",\n ggtheme = theme_minimal(base_size = 24), tl.cex = 14)\nplt1 + plt2","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-23T17:20:39.249292Z","iopub.execute_input":"2023-04-23T17:20:39.250304Z","iopub.status.idle":"2023-04-23T17:20:40.255075Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Comparison of two dendrograms**","metadata":{}},{"cell_type":"code","source":"fig(25, 15)\n# Compute distance matrix, hierarchical clusterings and create two dendrograms in a list\ndend1 <- set22_corr %>% dist %>% hclust(\"complete\") %>% as.dendrogram\ndend2 <- set23_corr %>% dist %>% hclust(\"complete\") %>% as.dendrogram\ndend_list <- dendlist(dend1, dend2)\n# Align and plot two dendrograms side by side\ndendlist(dend1, dend2) %>%\n untangle(method = \"step1side\") %>% # Find the best alignment layout\n tanglegram(cex_main = 3, lab.cex = 2,\n margin_inner = 15,\n common_subtrees_color_branches = TRUE,\n main_left = \"CITEseq 2022\",\n main_right = \"CITEseq 2023\")\ncor_res <- cor.dendlist(dend_list, method = \"cophenetic\")\ncat(\"Cophenetic correlation between the two dendrograms\", round(cor_res[2,1], 2))\nall.equal(dend1, dend2)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-23T17:20:40.256935Z","iopub.execute_input":"2023-04-23T17:20:40.257958Z","iopub.status.idle":"2023-04-23T17:20:43.485473Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Cluster 4","metadata":{}},{"cell_type":"markdown","source":"**Subset data from the two datasets**","metadata":{}},{"cell_type":"code","source":"RNA_names22 <- c(\"ENSG00000081237_PTPRC\",\n \"ENSG00000142541_RPL13A\",\n \"ENSG00000125691_RPL23\" ,\n \"ENSG00000166441_RPL27A\",\n\"ENSG00000071082_RPL31\" ,\n \"ENSG00000172809_RPL38\" ,\n \"ENSG00000174444_RPL4\", \n\"ENSG00000147604_RPL7\" ,\n \"ENSG00000142534_RPS11\",\n \"ENSG00000008988_RPS20\" ) %>% as.character\n\ncat(length(RNA_names22)-1, \"Found genes:\")\nRNA_names22[RNA_names22 != \"ENSG00000081237_PTPRC\"]\n\nset22 <- make_subset(RNA_names22, prot_names)\nRNA_names22 <- gsub(\".*_\",\"\", RNA_names22)\n\ncat(\"CITEseq2022: Head of the data frame with\", ncol(set22)-8, \"RNA/proteins in columns X\", nrow(set22), \"cells in rows\")\nhead(set22)\n\n#CITEseq2023\nRNA_names23 <- colnames(mat_RNA23)[which(grepl(paste0(gsub(\"_.*\", \"\", RNA_names22), collapse = \"$|\"),\n colnames(mat_RNA23),\n ignore.case = T))]\nprot_names23 <- colnames(mat_prot23)[grep(paste(prot_names, collapse = \"|\"), colnames(mat_prot23))]\n\nset23 <- cbind(mat_RNA23[, RNA_names23],\n mat_prot23[, prot_names23]) %>% as.data.frame\ncat(\"\\nCITEseq2023: Head of the data frame with\", ncol(set23), \"RNA/proteins in columns X\", nrow(set23), \"cells in rows\")\nhead(set23)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-23T17:20:43.487368Z","iopub.execute_input":"2023-04-23T17:20:43.488513Z","iopub.status.idle":"2023-04-23T17:20:44.449611Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Check expression levels**","metadata":{}},{"cell_type":"code","source":"fig(25,15)\ngene_set <- c(RNA_names22[!RNA_names22 == \"PTPRC\"])\n\nplot_expr_by_fct(set22, gene_set, factor_column = \"sample_size_ct\", ncol = 3)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-23T17:20:44.451332Z","iopub.execute_input":"2023-04-23T17:20:44.452400Z","iopub.status.idle":"2023-04-23T17:20:45.698675Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig(7,7)\n\nsum_expr_RNA23 <- apply(select(set23, -grep(\"CD45\", colnames(set23))), 2, median) %>% sort(decreasing = TRUE)\nsum_expr_RNA23 <- data.frame(\"Median_expression\" = sum_expr_RNA23,\n \"Name\" = names(sum_expr_RNA23))\nsum_expr_RNA23$Name <- factor(sum_expr_RNA23$Name, levels = unique(sum_expr_RNA23$Name ))\n#hline <- sum_expr_RNA23[sum_expr_RNA23$Name == \"PTPRC\",]$Median_expression\n\nggplot(sum_expr_RNA23, aes(x = Name, y = Median_expression)) +\n geom_bar(stat = 'identity', fill=\"#BF9039\", col=\"grey\") + \n# geom_hline(yintercept=sum_expr_RNA23[sum_expr_RNA23$Name == \"PTPRC\",]$Median_expression,\n# linetype=\"dashed\")+\n theme_bw(base_size = 22) +\n xlab(\"\") +\n ggtitle(paste0(\"CITEseq 2023\")) +\n theme(axis.text.x = element_text(angle = 45, hjust=1))\n\ncat(\"Zero median expression:\\n\")\ncat(paste0('\"', sum_expr_RNA23[sum_expr_RNA23$Median_expression == 0,]$Name, '\"', collapse = \", \"))","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-23T17:20:45.700760Z","iopub.execute_input":"2023-04-23T17:20:45.701938Z","iopub.status.idle":"2023-04-23T17:20:45.938213Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### CITEseq22","metadata":{}},{"cell_type":"code","source":"cat(length(prot_names), \"Proteins:\\n\")\ncat(paste0('\"', prot_names, '\"', collapse = \", \"))\ncat(\"\\n\\n\", length(RNA_names22), \"RNAs:\\n\")\ncat(paste0('\"', RNA_names22, '\"', collapse = \", \"))","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-23T17:20:45.941352Z","iopub.execute_input":"2023-04-23T17:20:45.943406Z","iopub.status.idle":"2023-04-23T17:20:45.966220Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig(25, 12)\nall_plt <- plot_corr_by_fct(set22, \"cell_type\", RNA_names22, prot_names, show_corr = FALSE, tl.cex = 16)\n\ncat(\"CITEseq 2022, Spearman correlation\")\ndo.call(\"grid.arrange\", c(all_plt, ncol = 4)) ","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-23T17:20:45.968431Z","iopub.execute_input":"2023-04-23T17:20:45.969840Z","iopub.status.idle":"2023-04-23T17:20:48.016713Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Compare dendrogram trees**","metadata":{}},{"cell_type":"code","source":"dendlst <- dendlist()\nfor(ct in unique(set22$cell_type)) {\n \n sbst <- set22 %>% filter(cell_type == ct) %>%\n select(all_of(RNA_names22), all_of(prot_names)) %>% \n cor(method = \"spearman\") %>% \n dist %>% hclust(\"complete\") %>% as.dendrogram\n dendlst[[ct]] <- sbst\n}\n\ndendlst[[\"All_cells\"]] <- set22 %>%\n select(all_of(RNA_names22), all_of(prot_names)) %>% \n cor(method = \"spearman\") %>% \n dist %>% hclust(\"complete\") %>% as.dendrogram\n\ncors <- cor.dendlist(dendlst)\n# Print correlation matrix\n#round(cors, 2)\nfig(25, 7)\nggcorrplot(cors, hc.order = TRUE, lab = TRUE,\n title = paste0(\"Correlation matrix between dendrograms by cell type\"), lab_size = 6,\n ggtheme = theme_minimal(base_size = 18), tl.cex = 18)\n\n#all.equal(dendlst$MasP, dendlst$All_cells)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-23T17:20:48.019070Z","iopub.execute_input":"2023-04-23T17:20:48.020544Z","iopub.status.idle":"2023-04-23T17:20:49.005973Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Differential correlations with RO and RA isoforms**","metadata":{}},{"cell_type":"code","source":"diff_cor <- diff_corr_by_fct(set22, factor_column = \"cell_type\", RNA_names = RNA_names22, prot_names)\ncat(\"Positive is for R >= 0.1, negative for R <= -0.1 and zero is |R| < 0.1\")\ndiff_cor %>% filter(cell_type %in% c(\"MasP\", \"NeuP\", \"All cells\")) #","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-23T17:20:49.008711Z","iopub.execute_input":"2023-04-23T17:20:49.010508Z","iopub.status.idle":"2023-04-23T17:20:49.697173Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Comparison of the two datasets","metadata":{}},{"cell_type":"code","source":"#calculate correlations\nsbst22 <- set22 %>% select(intersect(names(set23),names(set22)))\nsbst23 <- set23 %>% select(intersect(names(set23),names(set22)))\n\nset22_corr <- cor(sbst22, method = \"spearman\")\ncat(\"The correlation matrix CITEseq2022\")\nset22_corr\n\nset23_corr <- cor(sbst23, method = \"spearman\")\ncat(\"\\nThe correlation matrix CITEseq2023\")\nset23_corr","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-23T17:20:49.699926Z","iopub.execute_input":"2023-04-23T17:20:49.701142Z","iopub.status.idle":"2023-04-23T17:20:50.012826Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig(25,10)\nplt1 <- ggcorrplot(set22_corr, hc.order = FALSE, lab = TRUE,\n title = \"CITEseq 2022\", lab_size = 6, p.mat = NULL, #as.matrix(set22_p), insig = \"pch\", pch = 0,\n ggtheme = theme_minimal(base_size = 24), tl.cex = 18)\nplt2 <- ggcorrplot(set23_corr, hc.order = FALSE, lab = TRUE,\n title = \"CITEseq 2023\", lab_size = 6, p.mat = NULL, #set23_p, insig = \"pch\",\n ggtheme = theme_minimal(base_size = 24), tl.cex = 18)\nplt1 + plt2","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-23T17:20:50.015075Z","iopub.execute_input":"2023-04-23T17:20:50.016625Z","iopub.status.idle":"2023-04-23T17:20:50.753151Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Comparison of two dendrograms**","metadata":{}},{"cell_type":"code","source":"fig(25, 7)\n# Compute distance matrix, hierarchical clusterings and create two dendrograms in a list\ndend1 <- set22_corr %>% dist %>% hclust(\"complete\") %>% as.dendrogram\ndend2 <- set23_corr %>% dist %>% hclust(\"complete\") %>% as.dendrogram\ndend_list <- dendlist(dend1, dend2)\n# Align and plot two dendrograms side by side\ndendlist(dend1, dend2) %>%\n untangle(method = \"step1side\") %>% # Find the best alignment layout\n tanglegram(cex_main = 3, lab.cex = 2,\n margin_inner = 15,\n common_subtrees_color_branches = TRUE,\n main_left = \"CITEseq 2022\",\n main_right = \"CITEseq 2023\")\ncor_res <- cor.dendlist(dend_list, method = \"cophenetic\")\ncat(\"Cophenetic correlation between the two dendrograms\", round(cor_res[2,1], 2))\nall.equal(dend1, dend2)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-23T17:20:50.755632Z","iopub.execute_input":"2023-04-23T17:20:50.756847Z","iopub.status.idle":"2023-04-23T17:20:51.081787Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#","metadata":{"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]}]}