{"metadata":{"kernelspec":{"name":"ir","display_name":"R","language":"R"},"language_info":{"name":"R","codemirror_mode":"r","pygments_lexer":"r","mimetype":"text/x-r-source","file_extension":".r","version":"4.0.5"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"What we do here:\n\n- Memory efficient calculation of spearman correlations on sparse matrices for all cells, by day, donor, cell type, and by cell type for each donor*day and save in the output data\n- make interactive visualisation of most pronounced protein-protein/RNA correlations (**all cells**)\n- levels of the protein and RNA of interest (moved to [this notebook](http://www.kaggle.com/code/antoninadolgorukova/cd36-protein-rna-corr-expression-dropouts))\n- threshold of random variability in correlation coefficients (moved to [this notebook](http://www.kaggle.com/code/antoninadolgorukova/mmscel-corr-rand-threshold))\n\nData: RDS files with sparse matrices of [normalised counts data](http://www.kaggle.com/datasets/stautxie/sparse-measurement-data-open-problems-multimodal) or [raw counts](http://www.kaggle.com/datasets/antoninadolgorukova/sparse-raw-counts-data-open-problems-multimodal) for [Open Problems - Multimodal Single-Cell Integration]http://www.kaggle.com/competitions/open-problems-multimodal). The dataset for this competition comprises single-cell multiomics data collected from mobilized peripheral CD34+ hematopoietic stem and progenitor cells (HSPCs) isolated from four healthy human donors.\n\nFor interactive plots: d3 visualization library \"[visNetwork](http://datastorm-open.github.io/visNetwork/)\".\n\n**Versions:**  \n17: prot_of_interest = CD36  \n18: prot_of_interest = CD115  \n19: prot_of_interest = CD71  \n20: prot_of_interest = CD88  \n21: prot_of_interest = CD44  \n22: prot_of_interest = CD45  \n23: prot_of_interest = CD45RA","metadata":{}},{"cell_type":"code","source":"# Libraries\nsuppressMessages(library(Matrix))\nsuppressMessages(library(dplyr)) \nsuppressMessages(library(tidyr))\nsuppressMessages(library(rstatix))\n#suppressMessages(library(corrplot))\nsuppressMessages(library(ggplot2))\nsuppressMessages(library(tictoc))\nlibrary(visNetwork) # network visualisation\n\nremotes::install_github(\"cysouw/qlcMatrix\")\nsuppressMessages(library(qlcMatrix))\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,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Load data**","metadata":{}},{"cell_type":"markdown","source":"Targets: for the surface protein levels, each row corresponds to a cell (e.g. \"45006fe3e4c8\") and each column to a protein (e.g. \"CD86\").\n\nInputs: For the RNA counts, each row corresponds to a cell (e.g. \"45006fe3e4c8\") and each column to a gene. The column format for a gene is given by {EnsemblID}_{GeneName} where EnsemblID refers to the Ensembl Gene ID and GeneName to the gene name (e.g. \"ENSG00000159840_ZYX\").\n\nMetadata: Donor and cell types. The train data consists of both gene expression (RNA) and surface protein data for days 2,3,4 for donors 1-3 (donor IDs: 32606,13176, and 31800), the public test data consists of RNA for days 2,3,4 for donor 4 (donor ID: 27678) and the private test data consists data from day 7 from all donors.\n","metadata":{}},{"cell_type":"code","source":"tic() # ~ 1 min\n# Load Metadata (day, donor, cell type, and technology)\nmetadata <- read.csv('../input/open-problems-multimodal/metadata.csv',row.names=1)\nmetadata <- filter(metadata, technology == \"citeseq\" &\n                   donor %in% c(\"13176\", \"31800\", \"32606\") &\n                   day %in% c(\"2\",\"3\",\"4\")) %>%\n  mutate(\"Row.names\" = row.names(.))\n# Load protein data\n# raw data\n#path <- \"/kaggle/input/sparse-raw-counts-data-open-problems-multimodal/citeseq/sp_train_cite_targets_raw.rds\"\n# normalized data\npath <- \"/kaggle/input/sparse-measurement-data-open-problems-multimodal/sp_train_cite_targets.rds\"\nmat_prot <- readRDS(path)\n\n#dgCMatrix to matrix\n#mat_prot <- as.matrix(mat_prot)\n#mat_prot <- mat_prot[, prot_of_interest]\n\n# Load RNA data\n# raw data\n#path <- \"/kaggle/input/sparse-raw-counts-data-open-problems-multimodal/citeseq/sp_train_cite_inputs_raw.rds\"\n# normalized data\npath <- \"/kaggle/input/sparse-measurement-data-open-problems-multimodal/sp_train_cite_inputs.rds\"\nmat_RNA <- readRDS(path)\n\n# dgCMatrix to matrix\n#mat_RNA <- as.matrix(mat_RNA)\ntoc()\ngc()","metadata":{"_kg_hide-output":true,"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"For Memory efficient calculation of spearman correlations on sparse matrices we use Saket Choudhary's functions ([associated notebook](http://github.com/saketkc/blog/blob/main/2022-03-10/SparseSpearmanCorrelation2.ipynb), [explanation](http://saket-choudhary.me/blog/2022/03/09/sparsespearman/)) (THANK YOU).","metadata":{}},{"cell_type":"code","source":"# Memory efficient spearman correlation on sparse matrices\n\nSparsifiedRanks2 <- function(X) {\n  if (class(X)[1] != \"dgCMatrix\") {\n    X <- as(object = X, Class = \"dgCMatrix\")\n  }\n  non_zeros_per_col <- diff(x = X@p)\n  n_zeros_per_col <- nrow(x = X) - non_zeros_per_col\n  offsets <- (n_zeros_per_col - 1) / 2\n  x <- X@x\n  ## split entries to columns\n  col_lst <- split(x = x, f = rep.int(1:ncol(X), non_zeros_per_col))\n  ## calculate sparsified ranks and do shifting\n  sparsified_ranks <- unlist(x = lapply(X = seq_along(col_lst), \n                                        FUN = function(i) rank(x = col_lst[[i]]) + offsets[i]))\n  ## Create template rank matrix\n  X.ranks <- X\n  X.ranks@x <- sparsified_ranks\n  return(X.ranks)\n}\n\nSparseSpearmanCor2 <- function(X, Y = NULL, cov = FALSE) {\n\n  # Get sparsified ranks\n  rankX <- SparsifiedRanks2(X)\n  if (is.null(Y)){\n    # Calculate pearson correlation on rank matrices\n    return (corSparse(X=rankX, cov=cov))\n    }\n  rankY <- SparsifiedRanks2(Y)\n  return(corSparse( X = rankX, Y = rankY, cov = cov))\n}","metadata":{"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Sample sizes**","metadata":{}},{"cell_type":"markdown","source":"Let's see now many cells in the data and how many cells of different types was sampled at each day and in each donor.","metadata":{}},{"cell_type":"code","source":"#calculate sample size by cell_type\nout <- NULL\nfor (ct in unique(metadata$cell_type)) {\n    meta <- filter(metadata, cell_type == ct)\n\n    tmp <- data.frame(\"group\" = \"All data\",\n                      \"cell_type\" = ct,\n                      \"sample_size\" = n_distinct(meta$Row.names) )\n    out <- rbind(out, tmp)            \n}\nout_tab <- out %>%\n  pivot_wider(names_from = \"cell_type\", values_from = \"sample_size\")\n#out_tab\n\n#calculate sample size by cell_type for each day and donor\nfactor = c(\"donor\", \"day\")\nout <- NULL\nfor (ct in unique(metadata$cell_type)) {\n    for(fa in (factor)) {\n        metadata$filt <- metadata[, fa]\n        for(gr in unique(metadata$filt)) {\n            meta <- filter(metadata, cell_type == ct & filt == gr)\n\n            tmp <- data.frame(\"group\" = gr,\n                              \"cell_type\" = ct,\n                             \"sample_size\" = n_distinct(meta$Row.names) )\n            out <- rbind(out, tmp)\n        }            \n    }  \n}\n\nout <- pivot_wider(out, names_from = \"cell_type\", values_from = \"sample_size\")\nout_tab <- rbind(out_tab, out)\n#out_tab\n\n#calculate sample size for all data, each day, donor\nfactor = c(\"donor\", \"day\")\nout <- NULL\nfor(fa in (factor)) {\n    metadata$filt <- metadata[, fa]\n    for(gr in unique(metadata$filt)) {\n        meta <- filter(metadata, filt == gr)\n\n        tmp <- data.frame(\"group\" = gr,\n                          \"cell_type\" = \"ALL\",\n                         \"sample_size\" = n_distinct(meta$Row.names))\n        out <- rbind(out, tmp)    \n    }  \n}\nout <- pivot_wider(out, names_from = \"cell_type\", values_from = \"sample_size\")\nout <- rbind(data.frame(\n    \"group\" = \"All data\",\n    \"ALL\" = n_distinct(metadata$Row.names)),\n             out )\n#out\n#combine all\nout_tab <- full_join(out, out_tab, by = \"group\")\ncat(\"Sample_size: all data and  by cell type + day, donor\")\nout_tab","metadata":{"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Cell types:  \nMasP = Mast Cell Progenitor  \nMkP = Megakaryocyte Progenitor  \nNeuP = Neutrophil Progenitor  \nMoP = Monocyte Progenitor  \nEryP = Erythrocyte Progenitor  \nHSC = Hematoploetic Stem Cell  \nBP = B-Cell Progenitor  ","metadata":{"execution":{"iopub.status.busy":"2022-11-19T10:42:38.565155Z","iopub.execute_input":"2022-11-19T10:42:38.567236Z","iopub.status.idle":"2022-11-19T10:42:38.601462Z"}}},{"cell_type":"code","source":"rm(out,meta, tmp,fa,gr,ct, out_tab) ; gc()","metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Protein - protein correlations","metadata":{}},{"cell_type":"markdown","source":"# 1. All cells","metadata":{}},{"cell_type":"code","source":"tic()\nprot_prot_corr <- SparseSpearmanCor2(mat_prot)\nprot_prot_corr[lower.tri(prot_prot_corr, diag = TRUE)] <- 0\nprot_prot_corr <- as.data.frame(prot_prot_corr) \nnames(prot_prot_corr) <- colnames(mat_prot)\nprot_prot_corr$Protein1 <- colnames(mat_prot)\ntoc()\n\n# reshape to long format\nprot_prot_corr <-  pivot_longer(prot_prot_corr, cols = 1:(ncol(prot_prot_corr)-1),\n                                names_to = \"Protein2\", values_to = \"corr.coeff\") \n\n#remove lower part of a correlation matrix\nprot_prot_corr <- filter(prot_prot_corr, corr.coeff != 0)\n\ncat(\"The head of the data frame with\", nrow(prot_prot_corr), \"protein-protein pairs\",\n  \"\\nr (correlation coefficient) range:\",\n   round(min(prot_prot_corr$corr.coeff, na.rm = TRUE),2), \"to\",\n    round(max(prot_prot_corr$corr.coeff, na.rm = TRUE), 2), \"\\nmean =\",\n   round(mean(prot_prot_corr$corr.coeff, na.rm = TRUE), 4)  )\n\nhead(prot_prot_corr)","metadata":{"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# take a look at top 10 positive correlations between proteins\ncat(\"\\nTop 10 positively correlated proteins.\")\nhead(prot_prot_corr[order(prot_prot_corr$corr.coeff, decreasing = TRUE), ], n = 10)","metadata":{"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# take a look at top 10 negative correlations between proteins\ncat(\"\\nTop 10 negatively correlated proteins.\")\nhead(prot_prot_corr[order(prot_prot_corr$corr.coeff, decreasing = FALSE), ], n = 10)","metadata":{"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig(5,5)\nhist(prot_prot_corr$corr.coeff,\n    main = paste(\"Histogram of correlation coefficients\\nfor\", nrow(prot_prot_corr),\n                 \"protein-protein pairs\"),\n    xlab = \"Correlation coefficients\")","metadata":{"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Interactive visualisation of top correlations","metadata":{}},{"cell_type":"code","source":"# prepare data for igraph. Out is the data frame with\n# proteins, RNA and correlation coefficients\nprep_igraph <- function(out, width_coeff) {\n    \n    nodes <- \n      rbind(select(out, c(Protein1)) %>% rename(label = Protein1),\n            select(out, c(Protein2)) %>% rename(label = Protein2) ) %>%\n      unique %>%\n      mutate(id = 1:nrow(.))\n            \n    edges <- out %>%  \n      mutate(width = abs(corr.coeff)*width_coeff) %>%\n      left_join(nodes, by = c(\"Protein1\" = \"label\")) %>% \n      rename(from = id) %>% \n      left_join(nodes, by = c(\"Protein2\" = \"label\")) %>% \n      rename(to = id) %>%\n      select(from, to, corr.coeff, width)\n\n    out <- list(nodes, edges)\n    names(out) <- c(\"nodes\", \"edges\")\n    return(out)\n}\n\n#coloring of edges\ncolfunc_neg <- colorRampPalette(c(\"skyblue\",\"darkblue\"))\ncolfunc_pos <- colorRampPalette(c(\"red\", \"pink\"))","metadata":{"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# to show most pronounced correlations, set cut_off\ncut_off = 0.45\ndat_graph <- filter(prot_prot_corr, abs(corr.coeff) > cut_off)\ndat_graph <- dat_graph[order(abs(dat_graph$corr.coeff), decreasing = TRUE), ]\n\n# to show top correlations filter:\n#dat_graph <- dat_all[order(dat_all$corr.coeff, decreasing = TRUE), ]\n#dat_graph <- dat_all[1:10,] # top 10\n\n#dat_graph$group <- paste0(prot_of_interest)\ngraph <- prep_igraph(dat_graph, width_coeff = 10)\n\nedges <- graph$edges\nedges <- mutate(edges, association = \n                ifelse(edges$corr.coeff > 0, \"positive\",\"negative\"), \n                color = NA)\nedges$association <- factor(edges$association, levels = c( \"positive\",\"negative\"))\nedges[edges$association == \"positive\", ]$color <- \n  colfunc_pos(nrow(edges[edges$association == \"positive\", ]))\nedges[edges$association == \"negative\", ]$color <- \n  colfunc_neg(nrow(edges[edges$association == \"negative\", ]))\n\nledges <- data.frame(color = c(\"red\", \"blue\"),\n                     label = levels(edges$association))\n\ncat(nrow(edges), \"pairs of\",nrow(graph$nodes),\n    \"proteins\",\n    \"\\nThe hue and thickness of the lines correspond to the strength of the association (the darker and thicker, the stronger)\",\n    \"\\nr cutoff = \", cut_off)\nvisNetwork(graph$nodes, edges, width = \"100%\") %>%\nvisOptions(highlightNearest = TRUE, selectedBy = \"label\") %>%\nvisLegend(addEdges = ledges) \n  #visIgraphLayout(layout = \"layout_with_fr\")","metadata":{"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#save protein-RNA correlations\nwrite.csv(prot_prot_corr, 'Prot-prot_corr_all.csv', row.names = F)\n\n#remove unused variables\nfor (thing in ls()) { message(thing); print(object.size(get(thing)), units='auto') }\nrm(prot_prot_corr, ledges, edges, graph, dat_graph) ; gc()","metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# 2. By day, donor, cell type","metadata":{}},{"cell_type":"code","source":"tic() # ~30 sec\ncorr_fct_all <- NULL\nfor (fa in c(\"donor\", \"day\", \"cell_type\")) {\n    \n    metadata$filt <- metadata[, fa]\n    corr_by_gr <- NULL\n    for (gr in unique(metadata$filt)) {\n        meta <- filter(metadata, filt == gr)\n\n        out <- SparseSpearmanCor2(mat_prot[rownames(mat_prot) %in% meta$Row.names, ])\n        \n        out[lower.tri(out, diag = TRUE)] <- 0\n        out <- as.data.frame(out) \n        names(out) <- colnames(mat_prot)\n        out$Protein1 <- colnames(mat_prot)\n        out$factor = fa\n        out$group <- gr\n        out$sample_size <- nrow(mat_prot[rownames(mat_prot) %in% meta$Row.names, ])\n        \n        corr_by_gr <- rbind(corr_by_gr, out)\n    }\n    corr_fct_all <- rbind(corr_fct_all, corr_by_gr)\n}\n\ntoc()\n\n# order factor levels\ncorr_fct_all$group <- factor(corr_fct_all$group,\n                         levels = c(\"13176\", \"31800\",\"32606\",\n                                    \"2\",\"3\",\"4\",\n                                    \"EryP\", \"HSC\", \"MasP\",\"MkP\", \"MoP\", \"BP\",\"NeuP\"))\n#head(corr_fct_all)\n\n# reshape to long format\ncorr_fct_all <-  pivot_longer(corr_fct_all, cols = 1:(ncol(corr_fct_all)-4),\n                                names_to = \"Protein2\", values_to = \"corr.coeff\") \n\n#remove pairs made of the same protein\ncorr_fct_all <- filter(corr_fct_all, corr.coeff != 0)\n\ncat(\"The head of the data frame with\", nrow(corr_fct_all), \"protein-protein pairs\",\n  \"\\nr (correlation coefficient) range:\",\n   round(min(corr_fct_all$corr.coeff, na.rm = TRUE),2), \"to\",\n    round(max(corr_fct_all$corr.coeff, na.rm = TRUE), 2), \"\\nmean =\",\n   round(mean(corr_fct_all$corr.coeff, na.rm = TRUE), 4)  )\n\nhead(corr_fct_all)","metadata":{"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"First, let's see how many cell there are in each subgroup (sample_size) and for how many pairs (n.of.corrs) correlation coefficients are calculated.","metadata":{}},{"cell_type":"code","source":"corr_fct_all %>%\n group_by(factor, group, sample_size) %>%\n summarise(n.of.corrs = sum(!is.na(corr.coeff)), .groups = \"keep\")","metadata":{"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Interactive visualisation of top correlations","metadata":{}},{"cell_type":"markdown","source":"Interactive visualisation was made using d3 visualization library \"[visNetwork](http://datastorm-open.github.io/visNetwork/)\". The code can be easily adjusted with \"cut_off\" and \"n\" variables to show either most strong or top positive correlations.  \n​\nYou can select a protein (label) and see all the connections in each donor, day, cell type.","metadata":{}},{"cell_type":"code","source":"# several functions to make ighaphs and tables\n\n# data - a df with correlation coefficients for protein-protein|RNA pairs,\n# calculated for subgroups by donor, day, cell type (names in group).\n# factor - \"donor\", \"day\" or \"cell type\",\n# cut_off if a cutoff for correlation coefficients\n# n - max number of pairs to show\n\nprep_df <- function(data, fact, cut_off = 0, n) {\n    dat_graph <- filter(data, factor %in% fact & abs(corr.coeff) > cut_off)\n\n    all_tmp <- NULL\n    for (f in unique(dat_graph$group)) {\n\n        tmp <- dat_graph[dat_graph$group == f, ]\n        tmp <- tmp[order(abs(tmp$corr.coeff), decreasing = TRUE), ]\n        tmp <- head(tmp, n = n)\n        all_tmp <- rbind(all_tmp, tmp)\n    }\n    dat_graph <- all_tmp\n    #dat_graph <- rename(dat_graph, group = group)\n    return(dat_graph)\n}\n\n# print a table with the range and the mean of all shown correlation coefficients \ntab_for_igraph <- function(df) {\n    res <- dat_graph %>%\n      group_by(group, sample_size) %>%\n      summarize(n = n(),\n               corr_range = paste(round(min(corr.coeff), 2), \"-\",\n                                  round(max(corr.coeff), 2)),\n                mean_corr = round(mean(corr.coeff), 2),\n               .groups = \"keep\")\n    names(res)[names(res) == \"group\"] <- fact\n    return(res)\n}","metadata":{"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# prepare data for igraph. Out is the data frame with\n# proteins, and correlation coefficients\nprep_igraph <- function(out, width_coeff) {\n    \n    nodes <- \n      rbind(select(out, c(Protein1)) %>% rename(label = Protein1),\n            select(out, c(Protein2)) %>% rename(label = Protein2) ) %>%\n      unique %>%\n      mutate(id = 1:nrow(.))\n            \n    edges <- out %>%  \n      select(Protein1, Protein2, corr.coeff, group) %>%\n      mutate(width = abs(corr.coeff)*width_coeff) %>%\n      left_join(nodes, by = c(\"Protein1\" = \"label\")) %>% \n      rename(from = id) %>% \n      left_join(nodes, by = c(\"Protein2\" = \"label\")) %>% \n      rename(to = id) %>%\n      select(from, to, corr.coeff, width, group)\n\n    out <- list(nodes, edges)\n    names(out) <- c(\"nodes\", \"edges\")\n    return(out)\n}\n\n#coloring of edges\ncolfunc_neg <- colorRampPalette(c(\"skyblue\",\"darkblue\"))\ncolfunc_pos <- colorRampPalette(c(\"red\", \"pink\"))","metadata":{"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fact = \"cell_type\"\ncut_off = 0.2\nn = 20\n# to show most pronounced correlations, set cut_off, to show top, set n\ndat_graph <- prep_df(data = corr_fct_all, fact, cut_off, n)\n\n# make data for igraph\ngraph <- prep_igraph(dat_graph, width_coeff = 10)\n\n# color edges\nedges <- graph$edges\ncolors <- data.frame(\"labels\" = unique(edges$group),\n           color = 1:length(unique(edges$group)) )\n\nedges <- edges %>% left_join(colors, by = c(\"group\" = \"labels\"))\nedges$color <- palette()[edges$color]\n\nledges <- data.frame(color = unique(edges$color),\n                     label = unique(edges$group))\n\n# type some info about the graph\ncat(nrow(graph$edges), \"pairs of proteins with\",n_distinct(graph$nodes$label)-n_distinct(dat_graph$group),\"other proteins\",\n    \"\\nr cutoff = \", cut_off, \"\\nshown up to\", n, \"pairs per\", fact) \n\n# plot the graph\nvisNetwork(graph$nodes, edges, width = \"100%\")  %>%\n  visLegend(addEdges = ledges, main = fact) %>%\n  visOptions(highlightNearest = TRUE,\n             selectedBy = \"label\")\n\ntab_for_igraph(dat_graph)","metadata":{"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#save all protein-protein correlations by day, donor, cell type in csv\nwrite.csv(corr_fct_all, 'Prot-prot_corr_13gr.csv', row.names = F)\n#remove unused variables\nfor (thing in ls()) { message(thing); print(object.size(get(thing)), units='auto') }\nrm(dat_graph, edges, graph,corr_fct_all, corr_by_gr, fa, gr, meta)\ngc()","metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Protein - RNA correlations","metadata":{}},{"cell_type":"markdown","source":"Choose the protein of interest and store in prot_of_interest variable. The script below will show correlation analysis for this protein.","metadata":{}},{"cell_type":"code","source":"prot_of_interest <- \"CD45RA\"","metadata":{"_kg_hide-input":false,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# 1.All cells","metadata":{}},{"cell_type":"code","source":"tic() #~7 min\nall_corr <- SparseSpearmanCor2(mat_RNA, mat_prot)\nall_corr <- as.data.frame(all_corr) \nnames(all_corr) <- colnames(mat_prot)\nall_corr$RNA <- colnames(mat_RNA)\ntoc()\n\n# select pairs with the protein of interest\nprot_RNA_corr <- select(all_corr, all_of(prot_of_interest), RNA)\nprot_RNA_corr$protein <- prot_of_interest\nnames(prot_RNA_corr)[names(prot_RNA_corr) == prot_of_interest] <- \"corr.coeff\"\n\ncat(\"The head of the data frame with\", nrow(prot_RNA_corr), \"protein-RNA pairs\",\n  \"\\nr (correlation coefficient) range:\",\n   round(min(prot_RNA_corr$corr.coeff, na.rm = TRUE),2), \"to\",\n    round(max(prot_RNA_corr$corr.coeff, na.rm = TRUE), 2), \"\\nmean =\",\n   round(mean(prot_RNA_corr$corr.coeff, na.rm = TRUE), 4),\"\\nfor\",\n   nrow(prot_RNA_corr[is.na(prot_RNA_corr$corr.coeff),]),\n    \"the standard deviation of couns is zero, r = NA\"   )\n\nhead(prot_RNA_corr)\nprot_RNA_corr <- na.omit(prot_RNA_corr)","metadata":{"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#save protein-RNA correlations\nwrite.csv(all_corr, 'Prot-RNA_corr_all.csv', row.names = F)\nrm(all_corr) ; gc()","metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# take a look at top 10 positive correlations between the protein and RNA\ncat(\"\\nTop 10 positively correlated protein-RNA pairs.\")\nhead(prot_RNA_corr[order(prot_RNA_corr$corr.coeff, decreasing = TRUE), ], n = 10)","metadata":{"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# take a look at top 10 negative correlations between the protein and RNA\ncat(\"\\nTop 10 negatively correlated protein-RNA pairs.\")\nhead(prot_RNA_corr[order(prot_RNA_corr$corr.coeff, decreasing = FALSE), ], n = 10)","metadata":{"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig(5,5)\nhist(prot_RNA_corr$corr.coeff,\n    main = paste(\"Histogram of correlation coefficients\\nfor\", nrow(prot_RNA_corr),\n                 \"protein-RNA pairs\"),\n    xlab = \"Correlation coefficients\")","metadata":{"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Interactive visualisation of top correlations","metadata":{}},{"cell_type":"markdown","source":"According to calculations in [this notebook](https://www.kaggle.com/code/antoninadolgorukova/mmscel-corr-rand-threshold), let's assume that correlations with r > 0.2 are rather non-random, and see if they are stable across different donors, days, cell types.","metadata":{"_kg_hide-input":true}},{"cell_type":"markdown","source":"Interactive visualisation was made using d3 visualization library \"[visNetwork](http://datastorm-open.github.io/visNetwork/)\". The code can be easily adjusted with \"cut_off\" and \"n\" variables to show either most strong or top positive correlations.  \n\nYou can select a protein/RNA (label) and see all the connections in each donor, day, cell type.","metadata":{}},{"cell_type":"code","source":"# prepare data for igraph. Out is the data frame with\n# proteins, RNA and correlation coefficients\nprep_igraph <- function(out, width_coeff) {\n    \n    out <- mutate(out, protein = paste0(protein,\"(\",group,\")\" ))\n    nodes <- \n      rbind(select(out, c(RNA)) %>% rename(label = RNA) %>% mutate(group = \"RNA\"),\n            select(out, c(protein, group)) %>% rename(label = protein) ) %>%\n      unique %>%\n      mutate(id = 1:nrow(.))\n            \n    edges <- out %>%  \n      select(RNA, protein, corr.coeff, group) %>%\n      mutate(width = abs(corr.coeff)*width_coeff) %>%\n      left_join(nodes, by = c(\"RNA\" = \"label\")) %>% \n      rename(from = id) %>% \n      left_join(nodes, by = c(\"protein\" = \"label\")) %>% \n      rename(to = id) %>%\n      select(from, to, corr.coeff, width, group)\n\n    out <- list(nodes, edges)\n    names(out) <- c(\"nodes\", \"edges\")\n    return(out)\n}\n\n#coloring of edges\ncolfunc_neg <- colorRampPalette(c(\"skyblue\",\"darkblue\"))\ncolfunc_pos <- colorRampPalette(c(\"red\", \"pink\"))","metadata":{"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# to show most pronounced correlations, set cut_off\ncut_off = 0.2\ndat_graph <- filter(prot_RNA_corr, abs(corr.coeff) > cut_off)\n\n# to show top correlations filter:\n#dat_graph <- dat_all[order(dat_all$corr.coeff, decreasing = TRUE), ]\n#dat_graph <- dat_all[1:10,] # top 10\n\ndat_graph$group <- paste0(prot_of_interest)\ngraph <- prep_igraph(dat_graph, width_coeff = 10)\n\nedges <- graph$edges\nedges <- mutate(edges, association = \n                ifelse(edges$corr.coeff > 0, \"positive\",\"negative\"), \n                color = NA)\nedges$association <- factor(edges$association, levels = c( \"positive\",\"negative\"))\nedges[edges$association == \"positive\", ]$color <- \n  colfunc_pos(nrow(edges[edges$association == \"positive\", ]))\nedges[edges$association == \"negative\", ]$color <- \n  colfunc_neg(nrow(edges[edges$association == \"negative\", ]))\n\nledges <- data.frame(color = c(\"red\", \"blue\"),\n                     label = levels(edges$association))\n\ncat(\"Correlations between\", prot_of_interest,\n    \"protein and\",nrow(graph$nodes)-1, \"RNA\",\n    \"\\nThe hue and thickness of the lines correspond to the strength of the association (the darker and thicker, the stronger)\",\n    \"\\nr cutoff = \", cut_off)\nvisNetwork(graph$nodes, edges, width = \"100%\") %>%\nvisLegend(addEdges = ledges) \n  #visIgraphLayout(layout = \"layout_with_fr\")","metadata":{"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#remove unused variables\nfor (thing in ls()) { message(thing); print(object.size(get(thing)), units='auto') }\nrm(prot_RNA_corr,dat_graph, graph, edges) ; gc()","metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# 2.By day, donor, cell type","metadata":{}},{"cell_type":"code","source":"#calculate expression levels by cell_type\n\nmy_prot <- as.matrix(mat_prot)\nmy_prot <- my_prot[, prot_of_interest]\n\nout_prot <- NULL\nfor (ct in unique(metadata$cell_type)) {\n    \n    meta <- filter(metadata, cell_type == ct)\n\n    prot <- my_prot[names(my_prot) %in% meta$Row.names]\n    prot_expr <- mean(prot)\n    \n    tmp_prot <- data.frame(\"group\" = \"All data\",\n                           \"cell_type\" = ct,\n                           \"mean\" = prot_expr)\n    out_prot <- rbind(out_prot, tmp_prot)\n}\n\nout_tab_prot <- out_prot %>%\n  pivot_wider(names_from = \"cell_type\", values_from = \"mean\")\n#out_tab_prot\n\n#calculate expression levels by cell_type for each day and donor\nfactor = c(\"donor\", \"day\")\nout_prot <- NULL\nfor (ct in unique(metadata$cell_type)) {\n    for(fa in (factor)) {\n        metadata$filt <- metadata[, fa]\n        for(gr in unique(metadata$filt)) {\n            meta <- filter(metadata, cell_type == ct & filt == gr)\n\n            prot <- my_prot[names(my_prot) %in% meta$Row.names]\n            prot_expr <- mean(prot)\n            \n            tmp_prot <- data.frame(\"group\" = gr,\n                                   \"cell_type\" = ct,\n                                   \"mean\" = prot_expr)\n            out_prot <- rbind(out_prot, tmp_prot)\n        }            \n    }  \n}\n\nout_prot <- pivot_wider(out_prot, names_from = \"cell_type\", values_from = \"mean\")\nout_tab_prot <- rbind(out_tab_prot, out_prot)\n\n#out_tab_prot\n\n#calculate expression levels for all data, each day, donor\nfactor = c(\"donor\", \"day\")\nout_prot <- NULL\nfor(fa in (factor)) {\n    metadata$filt <- metadata[, fa]\n    for(gr in unique(metadata$filt)) {\n        meta <- filter(metadata, filt == gr)\n\n        prot <- my_prot[names(my_prot) %in% meta$Row.names]\n        prot_expr <- mean(prot)\n        \n        tmp_prot <- data.frame(\"group\" = gr,\n                               \"cell_type\" = \"ALL\",\n                               \"mean\" = prot_expr)\n        out_prot <- rbind(out_prot, tmp_prot) \n    }  \n}\nout_prot <- pivot_wider(out_prot, names_from = \"cell_type\", values_from = \"mean\")\nout_prot <- rbind(data.frame(\n    \"group\" = \"All data\",\n    \"ALL\" = mean(my_prot)), out_prot)\n#out_prot\n#combine all\nout_tab_prot <- full_join(out_prot, out_tab_prot, by = \"group\")\ncat(\"Mean expression levels of\",prot_of_interest, \": in all data and  by cell type + day, donor\")\nout_tab_prot","metadata":{"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"tic() # ~20 min\ncorr_fct_all <- NULL\nfor (fa in c(\"donor\", \"day\", \"cell_type\")) {\n    \n    metadata$filt <- metadata[, fa]\n    corr_by_gr <- NULL\n    for (gr in unique(metadata$filt)) {\n        meta <- filter(metadata, filt == gr)\n\n        out <- SparseSpearmanCor2(mat_RNA[rownames(mat_RNA) %in% meta$Row.names, ],\n                                  mat_prot[rownames(mat_prot) %in% meta$Row.names, ])\n        out <- as.data.frame(out) \n        names(out) <- colnames(mat_prot)\n        out$RNA <- colnames(mat_RNA)\n        out$factor = fa\n        out$group <- gr\n        out$sample_size <- nrow(mat_RNA[rownames(mat_RNA) %in% meta$Row.names, ])\n        corr_by_gr <- rbind(corr_by_gr, out)\n    }\n    corr_fct_all <- rbind(corr_fct_all, corr_by_gr)\n}\n\ntoc()\n\n# order factor levels\ncorr_fct_all$group <- factor(corr_fct_all$group,\n                         levels = c(\"13176\", \"31800\",\"32606\",\n                                    \"2\",\"3\",\"4\",\n                                    \"EryP\", \"HSC\", \"MasP\",\"MkP\", \"MoP\", \"BP\",\"NeuP\"))\n#head(corr_fct_all)","metadata":{"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"First, let's see how many cell there are in each subgroup (sample_size) and for how many pairs (n.of.corrs) correlation coefficients are calculated.","metadata":{}},{"cell_type":"code","source":"corr_fct_all %>%\n group_by(factor, group, sample_size) %>%\n summarise(n.of.corrs = sum(!is.na(CD36)), .groups = \"keep\")","metadata":{"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# select pairs with the protein of interest\ncorr_fct <- select(corr_fct_all, c(all_of(prot_of_interest), RNA, factor, group, sample_size) )\ncorr_fct$protein <- prot_of_interest\nnames(corr_fct)[names(corr_fct) == prot_of_interest] <- \"corr.coeff\"","metadata":{"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Interactive visualisation of top correlations","metadata":{}},{"cell_type":"code","source":"# several functions to make ighaphs and tables\n\n# data - a df with correlation coefficients for protein-protein|RNA pairs,\n# calculated for subgroups by donor, day, cell type (names in group).\n# factor - \"donor\", \"day\" or \"cell type\",\n# cut_off if a cutoff for correlation coefficients\n# n - max number of pairs to show\n\nprep_df <- function(data, fact, cut_off = 0, n) {\n    dat_graph <- filter(data, factor %in% fact & abs(corr.coeff) > cut_off)\n\n    all_tmp <- NULL\n    for (f in unique(dat_graph$group)) {\n\n        tmp <- dat_graph[dat_graph$group == f, ]\n        tmp <- tmp[order(abs(tmp$corr.coeff), decreasing = TRUE), ]\n        tmp <- head(tmp, n = n)\n        all_tmp <- rbind(all_tmp, tmp)\n    }\n    dat_graph <- all_tmp\n    #dat_graph <- rename(dat_graph, group = group)\n    return(dat_graph)\n}\n\n# print a table with the range and the mean of all shown correlation coefficients \ntab_for_igraph <- function(df) {\n    res <- dat_graph %>%\n      group_by(group, sample_size) %>%\n      summarize(n = n(),\n               corr_range = paste(round(min(corr.coeff), 2), \"-\",\n                                  round(max(corr.coeff), 2)),\n                mean_corr = round(mean(corr.coeff), 2),\n               .groups = \"keep\")\n    names(res)[names(res) == \"group\"] <- fact\n    return(res)\n}","metadata":{"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**By donor**","metadata":{}},{"cell_type":"code","source":"fact = \"donor\"\ncut_off = 0.2\nn = 50\n# to show most pronounced correlations, set cut_off, to show top, set n\ndat_graph <- prep_df(data = corr_fct, fact, cut_off, n)\n\n# make data for igraph\ngraph <- prep_igraph(dat_graph, width_coeff = 15)\n\n# color edges\nedges <- graph$edges\nedges <- edges[order(edges$corr.coeff, decreasing = TRUE), ]\nedges <- mutate(edges, association = \n                ifelse(edges$corr.coeff > 0, \"positive\",\"negative\"), \n                color = NA)\nedges$association <- factor(edges$association, levels = c( \"positive\",\"negative\"))\nedges[edges$association == \"positive\", ]$color <- \n  colfunc_pos(nrow(edges[edges$association == \"positive\", ]))\nedges[edges$association == \"negative\", ]$color <- \n  colfunc_neg(nrow(edges[edges$association == \"negative\", ]))\n\nledges <- data.frame(color = c(\"red\", \"blue\"),\n                     label = levels(edges$association))\n\n# type some info about the graph\ncat(nrow(graph$edges), \"pairs of\", prot_of_interest,\n    \"protein with\",n_distinct(graph$nodes$label)-n_distinct(dat_graph$group),\"unique RNA\",\n    \"\\nr cutoff = \", cut_off, \"\\nshown up to\", n, \"pairs per\", fact) \n\n# plot the graph\nvisNetwork(graph$nodes, edges, width = \"100%\")  %>%\n  visLegend(addEdges = ledges, main = fact) %>%\n  visOptions(highlightNearest = TRUE,\n             selectedBy = \"label\")\n\ntab_for_igraph(dat_graph)","metadata":{"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Top 100 correlations by donor with correlation coefficients (no cutoff) in brackets and their intersections.","metadata":{}},{"cell_type":"code","source":"fact = \"donor\"\ncut_off = -1\nn = 100\n# to show most pronounced correlations, set cut_off, to show top, set n\nres <- prep_df(data = corr_fct, fact, cut_off, n)\n\nres_don <- \n  group_by(res, group) %>%\n  summarize(RNA = paste0(gsub(\".*_\",\"\",RNA),\" (\", round(corr.coeff,2),\")\"),\n           .groups = \"keep\") %>%\n  pivot_wider(names_from = \"group\", values_from = RNA, values_fn = list)\n\nt1 <- gsub(\"\\\\s*\\\\([^\\\\)]+\\\\)\",\"\",unlist(res_don$`13176`))\nt2 <- gsub(\"\\\\s*\\\\([^\\\\)]+\\\\)\",\"\",unlist(res_don$`31800`))\nt3 <- gsub(\"\\\\s*\\\\([^\\\\)]+\\\\)\",\"\",unlist(res_don$`32606`))\n\nres_don$`In all donors` <- t1[t1 %in% t2 & t1 %in% t3] %>% list\nres_don\ncat(length(unlist(res_don$`In all donors`)), \"RNA are in the top 100 in each donor\")","metadata":{"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**By day**","metadata":{}},{"cell_type":"code","source":"fact = \"day\"\ncut_off = 0.2\nn = 50\n# to show most pronounced correlations, set cut_off, to show top, set n\ndat_graph <- prep_df(data = corr_fct, fact, cut_off, n)\n\n# make data for igraph\ngraph <- prep_igraph(dat_graph, width_coeff = 15)\n\n# color edges\nedges <- graph$edges\nedges <- edges[order(edges$corr.coeff, decreasing = TRUE), ]\nedges <- mutate(edges, association = \n                ifelse(edges$corr.coeff > 0, \"positive\",\"negative\"), \n                color = NA)\nedges$association <- factor(edges$association, levels = c( \"positive\",\"negative\"))\nedges[edges$association == \"positive\", ]$color <- \n  colfunc_pos(nrow(edges[edges$association == \"positive\", ]))\nedges[edges$association == \"negative\", ]$color <- \n  colfunc_neg(nrow(edges[edges$association == \"negative\", ]))\n\nledges <- data.frame(color = c(\"red\", \"blue\"),\n                     label = levels(edges$association))\n\n# type some info about the graph\ncat(nrow(graph$edges), \"pairs of\", prot_of_interest,\n    \"protein with\",n_distinct(graph$nodes$label)-n_distinct(dat_graph$group),\"unique RNA\",\n    \"\\nr cutoff = \", cut_off, \"\\nshown up to\", n, \"pairs per\", fact) \n\n# plot the graph\nvisNetwork(graph$nodes, edges, width = \"100%\")  %>%\n  visLegend(addEdges = ledges, main = fact) %>%\n  visOptions(highlightNearest = TRUE,\n             selectedBy = \"label\")\n\ntab_for_igraph(dat_graph)","metadata":{"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Top 100 correlations by day with correlation coefficients (no cutoff) in brackets and their intersections.","metadata":{}},{"cell_type":"code","source":"fact = \"day\"\ncut_off = -1\nn = 100\n# to show most pronounced correlations, set cut_off, to show top, set n\nres <- prep_df(data = corr_fct, fact, cut_off, n)\n\nres_day <- \n  group_by(res, group) %>%\n  summarize(RNA = paste0(gsub(\".*_\",\"\",RNA),\" (\", round(corr.coeff,2),\")\"),\n           .groups = \"keep\") %>%\n  pivot_wider(names_from = \"group\", values_from = RNA, values_fn = list)\n\nt1 <- gsub(\"\\\\s*\\\\([^\\\\)]+\\\\)\",\"\",unlist(res_day$`2`))\nt2 <- gsub(\"\\\\s*\\\\([^\\\\)]+\\\\)\",\"\",unlist(res_day$`3`))\nt3 <- gsub(\"\\\\s*\\\\([^\\\\)]+\\\\)\",\"\",unlist(res_day$`4`))\n\nres_day$`In all days` <- t1[t1 %in% t2 & t1 %in% t3] %>% list\nres_day\ncat(length(unlist(res_day$`In all days`)), \"RNA are in the top 100 in each day\")","metadata":{"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**By cell type**","metadata":{}},{"cell_type":"markdown","source":"1. With cutoff 0.2","metadata":{}},{"cell_type":"code","source":"fact = \"cell_type\"\ncut_off = 0.2\nn = 50\n# to show most pronounced correlations, set cut_off, to show top, set n\ndat_graph <- prep_df(data = corr_fct, fact, cut_off, n)\n\n# make data for igraph\ngraph <- prep_igraph(dat_graph, width_coeff = 15)\n\n# color edges\nedges <- graph$edges\nedges <- edges[order(edges$corr.coeff, decreasing = TRUE), ]\nedges <- mutate(edges, association = \n                ifelse(edges$corr.coeff > 0, \"positive\",\"negative\"), \n                color = NA)\nedges$association <- factor(edges$association, levels = c( \"positive\",\"negative\"))\nedges[edges$association == \"positive\", ]$color <- \n  colfunc_pos(nrow(edges[edges$association == \"positive\", ]))\nedges[edges$association == \"negative\", ]$color <- \n  colfunc_neg(nrow(edges[edges$association == \"negative\", ]))\n\nledges <- data.frame(color = c(\"red\", \"blue\"),\n                     label = levels(edges$association))\n\n# type some info about the graph\ncat(nrow(graph$edges), \"pairs of\", prot_of_interest,\n    \"protein with\",n_distinct(graph$nodes$label)-n_distinct(dat_graph$group),\"unique RNA\",\n    \"\\nr cutoff = \", cut_off, \"\\nshown up to\", n, \"pairs per\", fact) \n\n# plot the graph\nvisNetwork(graph$nodes, edges, width = \"100%\")  %>%\n  visLegend(addEdges = ledges, main = fact) %>%\n  visOptions(highlightNearest = TRUE,\n             selectedBy = \"label\")\n\ntab_for_igraph(dat_graph)","metadata":{"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"2. Plot top RNA 100 without cutoff.","metadata":{}},{"cell_type":"code","source":"fact = \"cell_type\"\ncut_off = -1\nn = 100\n# to show most pronounced correlations, set cut_off, to show top, set n\ndat_graph <- prep_df(data = corr_fct, fact, cut_off, n)\n\n# make data for igraph\ngraph <- prep_igraph(dat_graph, width_coeff = 15)\n\n# color edges\nedges <- graph$edges\nedges <- edges[order(edges$corr.coeff, decreasing = TRUE), ]\nedges <- mutate(edges, association = \n                ifelse(edges$corr.coeff > 0, \"positive\",\"negative\"), \n                color = NA)\nedges$association <- factor(edges$association, levels = c( \"positive\",\"negative\"))\nedges[edges$association == \"positive\", ]$color <- \n  colfunc_pos(nrow(edges[edges$association == \"positive\", ]))\nedges[edges$association == \"negative\", ]$color <- \n  colfunc_neg(nrow(edges[edges$association == \"negative\", ]))\n\nledges <- data.frame(color = c(\"red\", \"blue\"),\n                     label = levels(edges$association))\n\n# type some info about the graph\ncat(nrow(graph$edges), \"pairs of\", prot_of_interest,\n    \"protein with\",n_distinct(graph$nodes$label)-n_distinct(dat_graph$group),\"unique RNA\",\n    \"\\nr cutoff = \", cut_off, \"\\nshown up to\", n, \"pairs per\", fact) \n\n# plot the graph\nvisNetwork(graph$nodes, edges, width = \"100%\")  %>%\n  visLegend(addEdges = ledges, main = fact) %>%\n  visOptions(highlightNearest = TRUE,\n             selectedBy = \"label\")\n\ntab_for_igraph(dat_graph)","metadata":{"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Top 100 correlations by cell type with correlation coefficients (no cutoff) in brackets and their intersections.","metadata":{}},{"cell_type":"code","source":"fact = \"cell_type\"\ncut_off = -1\nn = 100\n# to show most pronounced correlations, set cut_off, to show top, set n\nres <- prep_df(data = corr_fct, fact, cut_off, n)\n\nres <- \n  group_by(res, group) %>%\n  summarize(RNA = paste0(gsub(\".*_\",\"\",RNA),\" (\", round(corr.coeff,2),\")\"),\n           .groups = \"keep\") %>%\n  pivot_wider(names_from = \"group\", values_from = RNA, values_fn = list)\n\nt1 <- gsub(\"\\\\s*\\\\([^\\\\)]+\\\\)\",\"\",unlist(res$`BP`))\nt2 <- gsub(\"\\\\s*\\\\([^\\\\)]+\\\\)\",\"\",unlist(res$`EryP`))\nt3 <- gsub(\"\\\\s*\\\\([^\\\\)]+\\\\)\",\"\",unlist(res$`HSC`))\nt4 <- gsub(\"\\\\s*\\\\([^\\\\)]+\\\\)\",\"\",unlist(res$`MasP`))\nt5 <- gsub(\"\\\\s*\\\\([^\\\\)]+\\\\)\",\"\",unlist(res$`MkP`))\nt6 <- gsub(\"\\\\s*\\\\([^\\\\)]+\\\\)\",\"\",unlist(res$`MoP`))\nt7 <- gsub(\"\\\\s*\\\\([^\\\\)]+\\\\)\",\"\",unlist(res$`NeuP`))\n\nres$`In all cells` <- t1[t1 %in% t2 & t1 %in% t3\n                       & t1 %in% t4\n                       & t1 %in% t5\n                       & t1 %in% t6\n                       & t1 %in% t7] %>% list\nres\ncat(length(unlist(res$`In all cells`)), \"RNA are in the top 100 in each cell type\")","metadata":{"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Let's see if there are intersections between top 100 RNA in different cell types: with which RNA the protein of interest correlates in most types of cells?**","metadata":{}},{"cell_type":"code","source":"fact = c(\"day\",\"donor\",\"cell_type\")\ncut_off = -1\nn = 100\n# to show most pronounced correlations, set cut_off, to show top, set n\nw_csv <- prep_df(data = corr_fct, fact, cut_off, n)\n\nintersections <- filter(w_csv, factor == \"cell_type\") %>%\n      select(-c(sample_size, factor)) %>%\n      pivot_wider(names_from = \"group\", values_from = \"corr.coeff\") %>%\n      as.data.frame %>%\n      mutate(intersections = rowSums(!is.na(.[-(1:2)])) )\n\ncat(\"There are:\", nrow(intersections), \"unique RNA\\n\",\n    \"RNA with more than 1 intersection:\", nrow(filter(intersections, intersections > 1)),\"\\n\",\n   \"RNA with at least 4 intersections:\", nrow(filter(intersections, intersections > 3)) )\n \nintersections[order(intersections$intersections, decreasing = T), ] %>%\n filter(intersections > 3)","metadata":{"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Let's see if there are intersections between RNA by day, donor, and most cell types**","metadata":{}},{"cell_type":"code","source":"res_cell <- intersections[intersections$intersections > 3, ]$RNA\nres_cell <- gsub(\".*_\",\"\",res_cell)\n\ncat(\"Stable\",length(unlist(res_don$`In all donors`)),\n    \"RNA\", \"by donor:\\n\\n\", unlist(res_don$`In all donors`))\ncat(\"\\n\\nStable\",length(unlist(res_day$`In all days`)),\n    \"RNA\", \"by day\\n\", unlist(res_day$`In all days`))\ncat(\"\\n\\nStable\",length(res_cell),\n    \"RNA\", \"in 4 cell types:\\n\", res_cell)\ncat(\"\\n\\nIntersections in day and donor and 4 cell types\")\nres_cell[res_cell %in% unlist(res_day$`In all days`) & res_cell %in% unlist(res_don$`In all donors`)]","metadata":{"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Write tables with correlations","metadata":{}},{"cell_type":"markdown","source":"**By donor, day, cell type - 13 groups**","metadata":{}},{"cell_type":"code","source":"#correlations by donor, day, cell type\ncorr_fct %>%\n group_by(factor,group, sample_size) %>%\n summarise(n.of.corrs = sum(!is.na(corr.coeff)), .groups = \"keep\")","metadata":{"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"write.csv(corr_fct_all, 'Prot-RNA_corr_13gr.csv', row.names = F)\n#remove unused variables\nfor (thing in ls()) { message(thing); print(object.size(get(thing)), units='auto') }\nrm(dat_graph,graph,intersections,w_csv)\ngc()","metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**For each cell type at each day - 21 groups**","metadata":{}},{"cell_type":"code","source":"tic() #~ 6 min\n\ncorr_fct_all <- NULL\nfor (ct in unique(metadata$cell_type)) {\n    corr_by_gr <- NULL\n    for (da in c(\"2\", \"3\", \"4\")) {\n\n        meta <- filter(metadata, cell_type == ct & day == da)\n\n        out <- SparseSpearmanCor2(mat_RNA[rownames(mat_RNA) %in% meta$Row.names, ],\n                                  mat_prot[rownames(mat_prot) %in% meta$Row.names, ])\n        \n        out <- as.data.frame(out) \n        names(out) <- colnames(mat_prot)\n        out$RNA <- colnames(mat_RNA)\n        out$day = da\n        out$cell_type <- ct\n        out$sample_size <- nrow(mat_RNA[rownames(mat_RNA) %in% meta$Row.names, ])\n        corr_by_gr <- rbind(corr_by_gr, out)\n       \n    }\n    corr_fct_all <- rbind(corr_fct_all, corr_by_gr)\n}\ntoc()\n\ncorr_fct_all %>%\n  group_by(cell_type, day,sample_size) %>%\n  summarise(n.of.corrs = sum(!is.na(CD36)), .groups = \"keep\") %>%\n  pivot_wider(names_from = \"cell_type\", values_from = c(\"n.of.corrs\"))","metadata":{"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"write.csv(corr_fct_all, 'Prot-RNA_corr_21gr.csv', row.names = F)","metadata":{"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**For each cell type in each donor at each day - 63 groups**","metadata":{}},{"cell_type":"code","source":"tic() #~ 7 min\n\ncorr_fct_all <- NULL\nfor (ct in unique(metadata$cell_type)) {\n    \n    corr_by_gr <- NULL\n    for (da in unique(metadata$day)) {\n        for(do in unique(metadata$donor)) {\n\n            meta <- filter(metadata, cell_type == ct & day == da & donor == do)\n            if(nrow(meta) >0) {\n                \n                out <- SparseSpearmanCor2(mat_RNA[rownames(mat_RNA) %in% meta$Row.names, ],\n                                      mat_prot[rownames(mat_prot) %in% meta$Row.names, ])\n\n                out <- as.data.frame(out) \n                names(out) <- colnames(mat_prot)\n                out$RNA <- colnames(mat_RNA)\n                out$day = da\n                out$donor = do\n                out$cell_type <- ct\n                out$sample_size <- nrow(mat_RNA[rownames(mat_RNA) %in% meta$Row.names, ])\n                corr_by_gr <- rbind(corr_by_gr, out)\n                \n            }            \n        }       \n    }\n    corr_fct_all <- rbind(corr_fct_all, corr_by_gr)   \n}\ntoc()\n\ncorr_fct_all %>%\n  group_by(cell_type, day, donor, sample_size) %>%\n  summarise(n.of.corrs = sum(!is.na(CD36)), .groups = \"keep\") %>%\n  pivot_wider(names_from = \"cell_type\", values_from = c(\"n.of.corrs\"))","metadata":{"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"write.csv(corr_fct_all, 'Prot-RNA_corr_63gr.csv', row.names = F)","metadata":{"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The tables are stored in the output.","metadata":{}}]}