{"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":"# Top 100 correlated RNA & proteins for CD36 and their intersections","metadata":{}},{"cell_type":"markdown","source":"What we do here:\n- make 4 data sets with top correlated (by absolute value of r): \n1. RNAs to CD36 RNA (n = 100 pairs)\n1. RNAs to CD36 protein (n = 100 pairs) \n1. Proteins to CD36 RNA (n = 20 pairs)\n1. Proteins to CD36 protein (n = 20 pairs)\n- visualize intersections in these data sets (Interactive d3, Heatmaps, UpSet plot) ","metadata":{}},{"cell_type":"markdown","source":"**Data:** RDS files with [sparse matrices of normalised counts](http://https//www.kaggle.com/datasets/stautxie/sparse-measurement-data-open-problems-multimodal) data for [Open Problems - Multimodal Single-Cell Integration](http://http//www.kaggle.com/competitions/open-problems-multimodal). The dataset for this competition comprises single-cell multiomics data collected from mobilized peripheral CD34+ hematopoietic stem and progenitor cells (HSPCs) isolated from four healthy human donors.","metadata":{}},{"cell_type":"code","source":"# Libraries\nsuppressMessages(library(Matrix))\nsuppressMessages(library(dplyr)) \nsuppressMessages(library(tidyr)) #pivot_longer/wider\nsuppressMessages(library(ggplot2))\nsuppressMessages(library(tictoc))\n\n# Function for figure size adjusment\nfig <- function(width, heigth) {\n    options(repr.plot.width = width, repr.plot.height = heigth) }","metadata":{"_uuid":"051d70d956493feee0c6d64651c6a088724dca2a","_execution_state":"idle","_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-12-11T11:33:16.347123Z","iopub.execute_input":"2022-12-11T11:33:16.349586Z","iopub.status.idle":"2022-12-11T11:33:17.833842Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Load data**","metadata":{}},{"cell_type":"markdown","source":"Targets: for the surface protein levels, each row corresponds to a cell (e.g. \"45006fe3e4c8\") and each column to a protein (e.g. \"CD86\").\n\nInputs: For the RNA counts, each row corresponds to a cell (e.g. \"45006fe3e4c8\") and each column to a gene. The column format for a gene is given by {EnsemblID}_{GeneName} where EnsemblID refers to the Ensembl Gene ID and GeneName to the gene name (e.g. \"ENSG00000159840_ZYX\").\n\nMetadata: Donor and cell types. The train data consists of both gene expression (RNA) and surface protein data for days 2,3,4 for donors 1-3 (donor IDs: 32606,13176, and 31800), the public test data consists of RNA for days 2,3,4 for donor 4 (donor ID: 27678) and the private test data consists data from day 7 from all donors.","metadata":{}},{"cell_type":"code","source":"# Load Metadata (day, donor, cell type, and technology)\nmetadata <- read.csv('../input/open-problems-multimodal/metadata.csv',row.names=1)\n\nmetadata <- \n  metadata %>% \n  filter(technology == \"citeseq\" &\n         donor %in% c(\"13176\", \"31800\", \"32606\") &\n         day %in% c(\"2\",\"3\",\"4\")) %>%\n  mutate_all(as.character) %>%\n  mutate(\"Row.names\" = row.names(.))\n\n# Load RNA normalized data\n#path <- \"/kaggle/input/sparse-raw-counts-data-open-problems-multimodal/citeseq/sp_train_cite_inputs_raw.rds\"\npath <- \"/kaggle/input/sparse-measurement-data-open-problems-multimodal/sp_train_cite_inputs.rds\"\nmat_RNA <- readRDS(path)\n\n# dgCMatrix to matrix\nmat_RNA <- as.matrix(mat_RNA)\n\n# Load protein normalized data\n#path <- \"/kaggle/input/sparse-raw-counts-data-open-problems-multimodal/citeseq/sp_train_cite_targets_raw.rds\"\npath <- \"/kaggle/input/sparse-measurement-data-open-problems-multimodal/sp_train_cite_targets.rds\"\nmat_prot <- readRDS(path)\n\n#dgCMatrix to matrix\nmat_prot <- as.matrix(mat_prot)\ngc()","metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Select the protein/RNA of interest**","metadata":{}},{"cell_type":"code","source":"# Find the protein and RNA of interest\nprot_name <- c(\"CD36\")\nRNA_name <- colnames(mat_RNA)[which(grepl(paste0(\"_\", prot_name,'$', collapse = \"|\"),\n                                     colnames(mat_RNA),\n                                     ignore.case = T))]","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-12-11T11:34:24.576593Z","iopub.execute_input":"2022-12-11T11:34:24.577928Z","iopub.status.idle":"2022-12-11T11:34:24.604467Z"},"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') }\ngc()","metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# 1. CD36 RNA - RNAs","metadata":{}},{"cell_type":"markdown","source":"Find the top 100 RNAs correlated with CD36 RNA","metadata":{}},{"cell_type":"markdown","source":"**Correlations in all cells**","metadata":{}},{"cell_type":"code","source":"# Function to calculate correlations in a loop (less RAM consuming, but slow)\nloop_corr_calc <- function(one_mol, all_mol, mol_of_interest) {\n    all_corr <- NULL\n    for (i in 1:ncol(all_mol)) {\n          \n        corr <- cor(one_mol, all_mol[, i], method = \"spearman\") %>%\n          suppressWarnings # “the standard deviation is zero” - means zero counts in all cells\n        #reshape to long format\n        corr <- \n          data.frame(\"var1\" = mol_of_interest,\n                    \"var2\" = colnames(all_mol)[i],\n                    \"corr.coeff\" = corr)\n        all_corr <- rbind(all_corr, corr)\n    }\n    return(all_corr)\n}","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-12-11T11:34:24.958969Z","iopub.execute_input":"2022-12-11T11:34:24.960488Z","iopub.status.idle":"2022-12-11T11:34:24.973480Z"},"jupyter":{"source_hidden":true},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"tic()\nrna_rna_corr <- loop_corr_calc(one_mol = mat_RNA[, RNA_name],\n                              all_mol = mat_RNA,\n                              mol_of_interest = RNA_name)\ntoc()\n#sort by absolute value of correlation coefficient\nrna_rna_corr <- \n  rna_rna_corr[order(abs(rna_rna_corr$corr.coeff), decreasing = TRUE), ] \n#rna_rna_corr <- filter(rna_rna_corr, var1 != var2) # to remove CD36-CD36 pair\nnames(rna_rna_corr) <- c(\"RNA1\", \"RNA2\", \"corr.coeff\")\n#print some info\ncat(nrow(rna_rna_corr), \"RNA-RNA pairs with correlations from\",\n   round(min(rna_rna_corr$corr.coeff, na.rm = TRUE),2), \"to\",\n    round(max(rna_rna_corr[rna_rna_corr$RNA1 != rna_rna_corr$RNA2,]$corr.coeff,\n              na.rm = TRUE), 2), \"\\nmean =\",\n   round(mean(rna_rna_corr$corr.coeff, na.rm = TRUE), 4),\"\\nfor\",\n   nrow(rna_rna_corr[is.na(rna_rna_corr$corr.coeff),]),\n    \"the standard deviation of couns is zero, r = NA\"   )\n\nrna_rna_corr <- na.omit(rna_rna_corr)\n\n# take a look at top 10 positive correlations\ncat(\"\\nTop 10 positively correlated pairs.\")\nhead(rna_rna_corr[order(rna_rna_corr$corr.coeff, decreasing = TRUE), ], n = 10)\n\n# take a look at top 10 negative correlations\ncat(\"\\nTop 10 negatively correlated pairs.\")\nhead(rna_rna_corr[order(rna_rna_corr$corr.coeff, decreasing = FALSE), ], n = 10)","metadata":{"_kg_hide-input":true,"jupyter":{"source_hidden":true},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#save top 100 pairs\ntop_rna_rna_corr <- rna_rna_corr[1:100, ]\nwrite.csv(top_rna_rna_corr, 'top100corr_СD36RNA-RNAs.csv', row.names = F)\ngc()","metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# 2. CD36 RNA - proteins","metadata":{}},{"cell_type":"markdown","source":"Find the top 100 proteins correlated with CD36 RNA","metadata":{}},{"cell_type":"markdown","source":"**Correlations in all cells**","metadata":{}},{"cell_type":"code","source":"tic()\nrna_prot_corr <- loop_corr_calc(one_mol = mat_RNA[, RNA_name],\n                              all_mol = mat_prot,\n                              mol_of_interest = RNA_name)\ntoc()\n#sort by absolute value of correlation coefficient\nrna_prot_corr <- rna_prot_corr[order(abs(rna_prot_corr$corr.coeff), decreasing = TRUE), ]\nnames(rna_prot_corr) <- c(\"RNA\", \"Protein\", \"corr.coeff\")\n#print some info\ncat(nrow(rna_prot_corr), \"RNA-protein pairs with correlations from\",\n   round(min(rna_prot_corr$corr.coeff, na.rm = TRUE),2), \"to\",\n    round(max(rna_prot_corr$corr.coeff, na.rm = TRUE), 2), \"\\nmean =\",\n   round(mean(rna_prot_corr$corr.coeff, na.rm = TRUE), 4),\"\\nfor\",\n   nrow(rna_prot_corr[is.na(rna_prot_corr$corr.coeff),]),\n    \"the standard deviation of couns is zero, r = NA\"   )\n\nrna_prot_corr <- na.omit(rna_prot_corr)\n\n# take a look at top 10 positive correlations\ncat(\"\\nTop 10 positively correlated pairs.\")\nhead(rna_prot_corr[order(rna_prot_corr$corr.coeff, decreasing = TRUE), ], n = 10)\n\n# take a look at top 10 negative correlations\ncat(\"\\nTop 10 negatively correlated pairs.\")\nhead(rna_prot_corr[order(rna_prot_corr$corr.coeff, decreasing = FALSE), ], n = 10)","metadata":{"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#save top 20 pairs\ntop_rna_prot_corr <- rna_prot_corr[1:20, ]\nwrite.csv(top_rna_prot_corr, 'top20corr_СD36RNA-prots.csv', row.names = F)\ngc()","metadata":{"_kg_hide-output":true,"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# 3. CD36 protein - RNAs","metadata":{}},{"cell_type":"markdown","source":"Find the top 100 RNAs correlated with CD36 protein","metadata":{}},{"cell_type":"markdown","source":"**Correlations in all cells**","metadata":{}},{"cell_type":"code","source":"tic()\nprot_rna_corr <- loop_corr_calc(one_mol = mat_prot[, prot_name],\n                              all_mol = mat_RNA,\n                              mol_of_interest = prot_name)\ntoc()\n#sort by absolute value of correlation coefficient\nprot_rna_corr <- prot_rna_corr[order(abs(prot_rna_corr$corr.coeff), decreasing = TRUE), ]\nnames(prot_rna_corr) <- c(\"Protein\", \"RNA\", \"corr.coeff\")\n#print some info\ncat(nrow(prot_rna_corr), \"Protein-RNA pairs with correlations from\",\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\nprot_rna_corr <- na.omit(prot_rna_corr)\n\n# take a look at top 10 positive correlations\ncat(\"\\nTop 10 positively correlated pairs.\")\nhead(prot_rna_corr[order(prot_rna_corr$corr.coeff, decreasing = TRUE), ], n = 10)\n\n# take a look at top 10 negative correlations\ncat(\"\\nTop 10 negatively correlated pairs.\")\nhead(prot_rna_corr[order(prot_rna_corr$corr.coeff, decreasing = FALSE), ], n = 10)","metadata":{"_kg_hide-input":true,"jupyter":{"source_hidden":true},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#save top 100 pairs\ntop_prot_rna_corr <- prot_rna_corr[1:100, ]\nwrite.csv(top_prot_rna_corr, 'top100corr_СD36prot-RNAs.csv', row.names = F)\ngc()","metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# 4. CD36 protein - proteins","metadata":{}},{"cell_type":"markdown","source":"Find the top 100 proteins correlated with CD36 protein","metadata":{}},{"cell_type":"markdown","source":"**Correlations in all cells**","metadata":{}},{"cell_type":"code","source":"tic()\nprot_prot_corr <- loop_corr_calc(one_mol = mat_prot[, prot_name],\n                              all_mol = mat_prot,\n                              mol_of_interest = prot_name)\ntoc()\n#sort by absolute value of correlation coefficient\nprot_prot_corr <- prot_prot_corr[order(abs(prot_prot_corr$corr.coeff), decreasing = TRUE), ]\n#prot_prot_corr <- filter(prot_prot_corr, var1 != var2)  # to remove CD36-CD36 pair\nnames(prot_prot_corr) <- c(\"Protein1\", \"Protein2\", \"corr.coeff\")\n#print some info\ncat(nrow(prot_prot_corr), \"Protein-Protein pairs with correlations from\",\n   round(min(prot_prot_corr$corr.coeff, na.rm = TRUE),2), \"to\",\n    round(max(prot_prot_corr[prot_prot_corr$Protein1 != prot_prot_corr$Protein2, ]$corr.coeff, na.rm = TRUE), 2), \"\\nmean =\",\n   round(mean(prot_prot_corr$corr.coeff, na.rm = TRUE), 4),\"\\nfor\",\n   nrow(prot_prot_corr[is.na(prot_prot_corr$corr.coeff),]),\n    \"the standard deviation of couns is zero, r = NA\"   )\n\nprot_prot_corr <- na.omit(prot_prot_corr)\n\n# take a look at top 10 positive correlations\ncat(\"\\nTop 10 positively correlated pairs.\")\nhead(prot_prot_corr[order(prot_prot_corr$corr.coeff, decreasing = TRUE), ], n = 10)\n\n# take a look at top 10 negative correlations\ncat(\"\\nTop 10 negatively correlated pairs.\")\nhead(prot_prot_corr[order(prot_prot_corr$corr.coeff, decreasing = FALSE), ], n = 10)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-12-11T12:50:21.858670Z","iopub.execute_input":"2022-12-11T12:50:21.860306Z","iopub.status.idle":"2022-12-11T12:50:27.078112Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#save top 30 pairs\ntop_prot_prot_corr <- prot_prot_corr[1:20, ]\nwrite.csv(top_prot_prot_corr, 'top20corr_СD36prot-prots.csv', row.names = F)\ngc()","metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# 5. Add labels to RNA (TF)","metadata":{}},{"cell_type":"markdown","source":"Load data table with Transcription Factors","metadata":{}},{"cell_type":"code","source":"tf_tab <- read.csv('/kaggle/input/genes-information/TranscriptionFactorsHuman_LambertDatabaseExtract_v_101.csv')\nhead(tf_tab, n = 3)","metadata":{"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#RNA\nrna_rna_corr <- \n  mutate(rna_rna_corr,\n         Ensembl.ID = gsub(\"_.*\", \"\", rna_rna_corr$RNA2),\n         is_in_tf_tab = ifelse(Ensembl.ID %in% tf_tab$Ensembl.ID, \"yes\", \"no\"))\nrna_rna_corr$is.TF <- NA\n\nfor(rna in (rna_rna_corr[rna_rna_corr$is_in_tf_tab == \"yes\", ]$Ensembl.ID)) {\n    \n    rna_rna_corr[rna_rna_corr$is_in_tf_tab == \"yes\" &\n                 rna_rna_corr$Ensembl.ID == rna, ]$is.TF <- \n    tolower(tf_tab[tf_tab$Ensembl.ID == rna, ]$Is.TF.)\n}\n\nrna_rna_corr[is.na(rna_rna_corr$is.TF), ]$is.TF <- \"unknown\"\nrna_rna_corr$RNA2 <- ifelse(rna_rna_corr$is.TF == \"yes\", paste0(rna_rna_corr$RNA2, \"<--- TF\"), rna_rna_corr$RNA2)\nrna_rna_corr <- select(rna_rna_corr, -Ensembl.ID)\ncat(rna_rna_corr %>% filter(is.TF == \"yes\") %>% nrow,\n    \"Transcription Factors found in the RNA data set (is TF? = yes)\" )\n      \n#Protein\nprot_rna_corr <- \n  mutate(prot_rna_corr,\n         Ensembl.ID = gsub(\"_.*\", \"\", prot_rna_corr$RNA),\n         is_in_tf_tab = ifelse(Ensembl.ID %in% tf_tab$Ensembl.ID, \"yes\", \"no\"))\nprot_rna_corr$is.TF <- NA\n\nfor(rna in (prot_rna_corr[prot_rna_corr$is_in_tf_tab == \"yes\", ]$Ensembl.ID)) {\n    \n    prot_rna_corr[prot_rna_corr$is_in_tf_tab == \"yes\" &\n                 prot_rna_corr$Ensembl.ID == rna, ]$is.TF <- \n    tolower(tf_tab[tf_tab$Ensembl.ID == rna, ]$Is.TF.)\n}\n\nprot_rna_corr[is.na(prot_rna_corr$is.TF), ]$is.TF <- \"unknown\"\nprot_rna_corr$RNA <- ifelse(prot_rna_corr$is.TF == \"yes\", paste0(prot_rna_corr$RNA, \"<--- TF\"), prot_rna_corr$RNA)\n\nprot_rna_corr <- select(prot_rna_corr, -Ensembl.ID)\ncat(prot_rna_corr %>% filter(is.TF == \"yes\") %>% nrow,\n    \"Transcription Factors found in the protein data set (is TF? = yes)\" )","metadata":{"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"cat(\"Transcription factors in the protein CD36 data set\")\nprot_rna_corr[1:20,]  %>% filter(is.TF == \"yes\")","metadata":{"execution":{"iopub.status.busy":"2022-12-11T20:35:05.003551Z","iopub.execute_input":"2022-12-11T20:35:05.005335Z","iopub.status.idle":"2022-12-11T20:35:05.047833Z"},"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"cat(\"Transcription factors in the RNA CD36 data set\")\nrna_rna_corr[1:100,]  %>% filter(is.TF == \"yes\")","metadata":{"execution":{"iopub.status.busy":"2022-12-11T20:35:00.491959Z","iopub.execute_input":"2022-12-11T20:35:00.495287Z","iopub.status.idle":"2022-12-11T20:35:00.553037Z"},"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# 5. Intersection visualization","metadata":{}},{"cell_type":"markdown","source":"An article about tools for intersection and visualization of multiple gene or genomic region sets: https://europepmc.org/article/PMC/5452382","metadata":{}},{"cell_type":"markdown","source":"## 5.1. Interactive d3 visualization","metadata":{}},{"cell_type":"markdown","source":"- Interactive d3 visualization library \"[visNetwork](https://cran.r-project.org/web/packages/visNetwork/vignettes/Introduction-to-visNetwork.html)\" \nThe plot shows top n RNA and proteins correlated with CD36 RNA and protein (separately for top positive and top negative correlations). In the plot you can select a protein/RNA (label) and what associations it has. With these plots it's easy to visualise inique associations for CD36 RNA and protein (cirqles that are connected to only one of them).","metadata":{}},{"cell_type":"code","source":"library(visNetwork) # network visualisation","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-12-11T12:50:51.689776Z","iopub.execute_input":"2022-12-11T12:50:51.691418Z","iopub.status.idle":"2022-12-11T12:50:51.706681Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Combine all positive correlations (for visualisation we use top = num)\nnum = 20\ndat_igraph <- list(rna_rna_corr, rna_prot_corr,\n                  prot_rna_corr, prot_prot_corr)\n\ndat_igraph <- \n  lapply(dat_igraph, function(x){\n    x %>% filter(corr.coeff > 0) %>% head(n = num) %>%\n    rename_with(.cols = 1, ~\"label\") %>%\n    rename_with(.cols = 2, ~\"var2\")\n})\n\n# Prepare data for igraph\nnodes <- \n  lapply(dat_igraph, function(x){\n      rbind(x[1], rename(x[2], label = var2) ) }) %>% \n  bind_rows %>%\n  unique %>% \n  mutate(id = 1:nrow(.),\n         group = ifelse(grepl(\"ENSG00\", label), \"RNA\", \"Protein\"))\n \nedges <- bind_rows(dat_igraph) %>%  \n  mutate(width = corr.coeff) %>%\n  left_join(nodes, by = c(\"label\")) %>% \n  rename(from = id) %>% \n  left_join(nodes, by = c(\"var2\" = \"label\")) %>% \n  rename(to = id) %>%\n  select(from, to, corr.coeff, width)\n\n#igraph\ncat(\"Top\",num,\"positive correlations between\", RNA_name, \"and\", prot_name,\n    \"vs all RNA and proteins.\\n\", nrow(nodes), \"unique entries\\n\",\n   nrow(edges), \"pairs\")\nvisNetwork(nodes, edges) %>%\nvisOptions(highlightNearest = TRUE,\n             selectedBy = \"label\") %>%\nvisLegend() #%>%\n#visIgraphLayout(layout = \"layout_with_fr\")","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-12-11T20:12:35.452546Z","iopub.execute_input":"2022-12-11T20:12:35.455467Z","iopub.status.idle":"2022-12-11T20:12:35.720321Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Combine all negative correlations (for visualisation we use top = num)\nnum = 20\ndat_igraph <- list(rna_rna_corr, rna_prot_corr,\n                  prot_rna_corr, prot_prot_corr)\n\ndat_igraph <- \n  lapply(dat_igraph, function(x){\n    x %>% filter(corr.coeff < 0) %>% head(n = num) %>%\n    rename_with(.cols = 1, ~\"label\") %>%\n    rename_with(.cols = 2, ~\"var2\")\n})\n\n# Prepare data for igraph\nnodes <- \n  lapply(dat_igraph, function(x){\n      rbind(x[1], rename(x[2], label = var2) ) }) %>% \n  bind_rows %>%\n  unique %>% \n  mutate(id = 1:nrow(.),\n         group = ifelse(grepl(\"ENSG00\", label), \"RNA\", \"Protein\"))\n \nedges <- bind_rows(dat_igraph) %>%  \n  mutate(width = corr.coeff) %>%\n  left_join(nodes, by = c(\"label\")) %>% \n  rename(from = id) %>% \n  left_join(nodes, by = c(\"var2\" = \"label\")) %>% \n  rename(to = id) %>%\n  select(from, to, corr.coeff, width)\n\n#igraph\ncat(\"Top\",num,\"negative correlations between\", RNA_name, \"and\", prot_name,\n    \"vs all RNA and proteins.\\n\", nrow(nodes), \"unique entries\\n\",\n   nrow(edges), \"pairs\")\nvisNetwork(nodes, edges) %>%\nvisOptions(highlightNearest = TRUE,\n             selectedBy = \"label\") %>%\nvisLegend() #%>%\n#visIgraphLayout(layout = \"layout_with_fr\")","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-12-11T20:14:15.657998Z","iopub.execute_input":"2022-12-11T20:14:15.660888Z","iopub.status.idle":"2022-12-11T20:14:16.001008Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## 5.2. ComplexHeatmap","metadata":{}},{"cell_type":"markdown","source":"- [ComplexHeatmap package](https://jokergoo.github.io/ComplexHeatmap-reference/book/upset-plot.html) - to visualize associations between different sources of data sets and reveal potential patterns (look for intersections in full top 100 lists).","metadata":{}},{"cell_type":"code","source":"suppressMessages(remotes::install_github(\"jokergoo/ComplexHeatmap\")) #takes time to install\nlibrary(ComplexHeatmap)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-12-11T11:57:44.741580Z","iopub.execute_input":"2022-12-11T11:57:44.742952Z","iopub.status.idle":"2022-12-11T12:01:06.860722Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#RNA\n# Find intersections\nnum = 100\nRNA <- rbind(rename(rna_rna_corr[1:num,], label = RNA1, RNA = RNA2),\n            rename(prot_rna_corr[1:num,], label = Protein)) %>%\n       pivot_wider(names_from = \"label\", values_from = \"corr.coeff\") #%>% na.omit\n \nRNA_res <- as.matrix(select(RNA, c(ENSG00000135218_CD36, CD36)) )\nrownames(RNA_res) <- RNA$RNA\n\n#print some stats\ndstats <- pivot_longer(RNA, cols = c(\"ENSG00000135218_CD36\", \"CD36\")) %>% na.omit\ndstats$association = ifelse(dstats$value > 0, \"positive\", \"negative\")\ndstats %>% group_by(name, association) %>% summarise(n = n(),\n                                             median = median(value),\n                                             min = min(value),\n                                             max = max(value), .groups = \"keep\")\n#plot heatmap\nfig(12,25)\nHeatmap(RNA_res, name = \"corr.coeff\",na_col = \"black\",\n       row_title = \"Top 100 RNAs correlated with CD36 protein and RNA\",\n       row_title_gp = gpar(fontsize = 20), cluster_rows = FALSE,\n       heatmap_legend_param = list(legend_height = unit(6, \"cm\"),\n                                   grid_width = unit(0.5, \"cm\"),\n                                   at = c(-1,-0.5,-0.2, 0, 0.2,0.5, 1),\n                                  labels_gp = gpar(fontsize = 16))\n)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-12-11T20:37:43.876977Z","iopub.execute_input":"2022-12-11T20:37:43.878889Z","iopub.status.idle":"2022-12-11T20:37:44.609679Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Contrast assosiations can also be checked with simple filtration of the data","metadata":{}},{"cell_type":"code","source":"cat(\"RNAs with opposite associations with RNA and protein\")\nselect(RNA, -c(is_in_tf_tab, is.TF)) %>% filter(ENSG00000135218_CD36 > 0 & CD36 < 0 |\n               ENSG00000135218_CD36 < 0 & CD36 > 0)\n\ncat(\"\\nRNAs with r for RNA and protein different by > o.2\")\nselect(RNA, -c(is_in_tf_tab, is.TF)) %>% filter(abs(ENSG00000135218_CD36 - CD36) > 0.2)","metadata":{"execution":{"iopub.status.busy":"2022-12-11T20:39:32.162519Z","iopub.execute_input":"2022-12-11T20:39:32.165423Z","iopub.status.idle":"2022-12-11T20:39:32.298054Z"},"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#Protein\n# Find intersections\nnum = 20\nprot <- rbind(rename(rna_prot_corr[1:num,], label = RNA),\n            rename(prot_prot_corr[1:num,], label = Protein1, Protein = Protein2)) %>%\n       pivot_wider(names_from = \"label\", values_from = \"corr.coeff\") #%>% na.omit\nprot_res <- as.matrix(select(prot, c(ENSG00000135218_CD36, CD36)) )\nrownames(prot_res) <- prot$Protein\n\n#print some stats\ndstats <- pivot_longer(prot, cols = c(\"ENSG00000135218_CD36\", \"CD36\")) %>% na.omit\ndstats$association = ifelse(dstats$value > 0, \"positive\", \"negative\")\ndstats %>% group_by(name, association) %>% summarise(n = n(),\n                                             median = median(value),\n                                             min = min(value),\n                                             max = max(value), .groups = \"keep\")\n#plot heatmap\nfig(12,8)\n#fig(12,25)\nHeatmap(prot_res, name = \"corr.coeff\",na_col = \"black\",\n       row_title = \"Top 20 proteins correlated with CD36 protein and RNA\",\n       row_title_gp = gpar(fontsize = 20), cluster_rows = FALSE,\n       heatmap_legend_param = list(legend_height = unit(6, \"cm\"),\n                                   grid_width = unit(0.5, \"cm\"),\n                                   at = c(-1,-0.5,-0.2, 0, 0.2,0.5, 1),\n                                  labels_gp = gpar(fontsize = 16))\n)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-12-11T20:45:53.956889Z","iopub.execute_input":"2022-12-11T20:45:53.961365Z","iopub.status.idle":"2022-12-11T20:45:54.378045Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Contrast assosiations can also be checked with simple filtration of the data","metadata":{}},{"cell_type":"code","source":"cat(\"Proteins with opposite associations with RNA and protein\")\nprot %>% filter(ENSG00000135218_CD36 > 0 & CD36 < 0 |\n               ENSG00000135218_CD36 < 0 & CD36 > 0)\n\ncat(\"Proteins with r for RNA and protein different by > o.2\")\nprot %>% filter(abs(ENSG00000135218_CD36 - CD36) > 0.2)","metadata":{"execution":{"iopub.status.busy":"2022-12-11T20:41:08.197053Z","iopub.execute_input":"2022-12-11T20:41:08.199814Z","iopub.status.idle":"2022-12-11T20:41:08.364723Z"},"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## 5.3. RVenn (setmap)","metadata":{}},{"cell_type":"markdown","source":"- [RVenn](https://cran.r-project.org/web/packages/RVenn/index.html) - package for set operations on multiple sets; setmap function shows the presence/absence of the elements among all the sets and cluster both the sets and the elements based on Jaccard distances.","metadata":{}},{"cell_type":"code","source":"suppressMessages(library(RVenn))","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-12-11T12:01:11.672368Z","iopub.execute_input":"2022-12-11T12:01:11.674611Z","iopub.status.idle":"2022-12-11T12:01:11.718389Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"num = 100\nintr <- list(\"RNA CD36\" = rna_rna_corr[1:num,]$RNA2,\n             \"Prot CD36\" = prot_rna_corr[1:num,]$RNA)  \n#str(intr)\nintrVenn <- Venn(intr)\n\ncat(\"The union of the top 100 RNAs:\", unite(intrVenn) %>% length)\nfig(12,20)\nsetmap(intrVenn, title = \"A clustered heatmap showing presence/absence of the RNAs in the sets for CD36 protein and its RNA\",\n      set_fontsize = 16, element_fontsize = 12)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-12-11T20:42:59.871798Z","iopub.execute_input":"2022-12-11T20:42:59.873653Z","iopub.status.idle":"2022-12-11T20:43:00.321449Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"num = 20\nintr <- list(\"RNA CD36\" = rna_prot_corr[1:num,]$Protein,\n             \"Prot CD36\" = prot_prot_corr[1:num,]$Protein2)  \n#str(intr)\nintrVenn <- Venn(intr)\n\ncat(\"The union of the top 20 Proteins:\", unite(intrVenn) %>% length)\nfig(12,8)\nsetmap(intrVenn, title = \"A clustered heatmap showing presence/absence of the proteins in the sets for CD36 protein and its RNA\",\n      set_fontsize = 16, element_fontsize = 12)","metadata":{"execution":{"iopub.status.busy":"2022-12-11T20:46:09.500886Z","iopub.execute_input":"2022-12-11T20:46:09.503506Z","iopub.status.idle":"2022-12-11T20:46:09.705611Z"},"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## 5.4. UpSet","metadata":{}},{"cell_type":"markdown","source":"UpSet plot provides an efficient way to visualize intersections of multiple sets compared to the traditional approaches, i.e. the Venn Diagram. It is implemented in the [UpSetR package](http://cran.r-project.org/web/packages/UpSetR/index.html) in R and re-implemented with the [ComplexHeatmap package](http://jokergoo.github.io/ComplexHeatmap-reference/book/upset-plot.html) with some improvements. Nice explanation is [here](http://upset.app/)","metadata":{}},{"cell_type":"code","source":"num = 100\nupst_rna <- list(\"RNA CD36 - RNAs\" = rna_rna_corr[1:num,]$RNA2,\n                 \"Prot CD36 - RNAs\" = prot_rna_corr[1:num,]$RNA)\n#str(upst_rna)\nupst_rna <- make_comb_mat(upst_rna)\n\ncat(\"Intersecions of correlated RNAs between the sets for CD36 protein and its RNA\")\nfig(7,5)\nUpSet(upst_rna, top_annotation = upset_top_annotation(upst_rna, add_numbers = TRUE),\n     pt_size = unit(7, \"mm\"), lwd = 2)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-12-11T20:47:38.238108Z","iopub.execute_input":"2022-12-11T20:47:38.240948Z","iopub.status.idle":"2022-12-11T20:47:38.522877Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"num = 20\nupst_prot <- list(\"RNA CD36 - prot\" = rna_prot_corr[1:num,]$Protein,\n                  \"Prot CD36 - prot\" = prot_prot_corr[1:num,]$Protein2)\n#str(upst_prot)\nupst_prot <- make_comb_mat(upst_prot)\n\ncat(\"Intersecions of correlated Proteins between the sets for CD36 protein and its RNA\")\nfig(7,5)\nUpSet(upst_prot, top_annotation = upset_top_annotation(upst_prot, add_numbers = TRUE),\n     pt_size = unit(7, \"mm\"), lwd = 2)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-12-11T20:48:10.170833Z","iopub.execute_input":"2022-12-11T20:48:10.172411Z","iopub.status.idle":"2022-12-11T20:48:10.416643Z"},"trusted":true},"execution_count":null,"outputs":[]}]}