{"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":"# Exploratory analysis of several meaningful and meaningless associations","metadata":{}},{"cell_type":"markdown","source":"Continuing to study correlations  srarted [earlier](http://https://www.kaggle.com/code/antoninadolgorukova/mmscel-corr-analysis-with-r)\n\nHere I explore most pronounced correlations between\n- all proteins and RNA \n- several interesting proteins (found by [Alexander Chervov](http://https://www.kaggle.com/alexandervc)) and RNA\n\nI use RDS files with [sparce matrices of raw couts data](http://https://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) and randomly select 10% of cells.","metadata":{}},{"cell_type":"code","source":"# Libraries and some functions\nsuppressMessages(library(Matrix))\nsuppressMessages(library(dplyr)) \nsuppressMessages(library(tidyr))\nlibrary(ggstatsplot)\nsuppressMessages(library(rstatix))\nsuppressMessages(library(corrplot))\nlibrary(ggplot2)\nlibrary(tictoc)\nlibrary(visNetwork) # network visualisation\n\n# function for figure size adjusment\nfig <- function(width, heigth) {\n    options(repr.plot.width = width, repr.plot.height = heigth)\n}","metadata":{"execution":{"iopub.status.busy":"2022-11-26T19:53:00.185005Z","iopub.execute_input":"2022-11-26T19:53:00.188500Z","iopub.status.idle":"2022-11-26T19:53:02.210406Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Randomly select 10% of cells to save memory footprint\nset.seed(1)\nselected_prop <- 0.1 #select 10% of cells\nselected_cells <- sample(1:70988,\n                        floor(selected_prop * 70988),\n                        replace = FALSE)","metadata":{"execution":{"iopub.status.busy":"2022-11-26T19:53:06.798004Z","iopub.execute_input":"2022-11-26T19:53:06.838494Z","iopub.status.idle":"2022-11-26T19:53:06.860927Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Load data","metadata":{}},{"cell_type":"markdown","source":"For the surface protein levels, each row corresponds to a cell and each column to a protein (e.g. \"CD86\").","metadata":{}},{"cell_type":"code","source":"#load protein data\npath <- \"/kaggle/input/sparse-raw-counts-data-open-problems-multimodal/citeseq/sp_train_cite_targets_raw.rds\"\ndat_prot <- readRDS(path)\n\n#dgCMatrix to dataframe\ndat_prot <- as.matrix(dat_prot)\ndat_prot <- as.data.frame(dat_prot[selected_cells, ]) #select 10% of cells\n\n# take a look\ncat(\"\\nA data frame with\", ncol(dat_prot),\n    \"proteins in columns and\", nrow(dat_prot), \"cells in rows\")\nhead(dat_prot)","metadata":{"execution":{"iopub.status.busy":"2022-11-26T19:53:12.710628Z","iopub.execute_input":"2022-11-26T19:53:12.712563Z","iopub.status.idle":"2022-11-26T19:53:14.221307Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"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\").","metadata":{}},{"cell_type":"code","source":"#load RNA data\npath <- \"/kaggle/input/sparse-raw-counts-data-open-problems-multimodal/citeseq/sp_train_cite_inputs_raw.rds\"\ndat_RNA <- readRDS(path)\n\n#dgCMatrix to dataframe\ndat_RNA <- as.matrix(dat_RNA)\ndat_RNA <- as.data.frame(dat_RNA[selected_cells, ])  #select 10% of cells\n\n# select genes with at least one value in rows\ntic()\ndat_RNA <- \n  dat_RNA %>%\n  bind_rows(summarise_all(., ~sum(.))) # add summary row\ndat_RNA <- select(dat_RNA, which(dat_RNA[nrow(dat_RNA), ] != 0) ) # remove columns with all 0s\ndat_RNA <- dat_RNA[-nrow(dat_RNA), ] # remove summary row (overall ~1 min)\ntoc()\n\n# take a look\ncat(\"\\nA data frame with\", ncol(dat_RNA),\n    \"RNA in columns with at least one value and\", nrow(dat_RNA), \"cells in rows\")\nhead(dat_RNA)","metadata":{"execution":{"iopub.status.busy":"2022-11-26T19:53:20.183528Z","iopub.execute_input":"2022-11-26T19:53:20.185124Z","iopub.status.idle":"2022-11-26T19:55:01.896483Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"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\ncat(\"\\nA data frame with the indication of day,\ndonor, cell type, and technology in columns for each of the\",nrow(metadata), \"cells in rows\")\ncat(\"For citeseq technology, the number of cells per:\")\ncat(\"\\nDonor ID:\")\ntable(metadata[metadata$technology == \"citeseq\",]$donor)\ncat(\"\\nDay:\")\ntable(metadata[metadata$technology == \"citeseq\",]$day)\ncat(\"\\nCell type:\")\ntable(metadata[metadata$technology == \"citeseq\",]$cell_type)","metadata":{"execution":{"iopub.status.busy":"2022-11-26T19:55:16.776479Z","iopub.execute_input":"2022-11-26T19:55:16.778273Z","iopub.status.idle":"2022-11-26T19:55:17.611617Z"},"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":"markdown","source":"Now we have 2 data frames with raw counts of:\n1. protein levels counts - dat_prot\n2. RNA expression counts - dat_RNA  \n\nAnd a data frame with meta data - metadata","metadata":{}},{"cell_type":"markdown","source":"# Сorrelations between proteins and RNA","metadata":{}},{"cell_type":"code","source":"# Calculate correlations between all proteins and all RNA\ntic()\nprot_RNA_mcorr <- cor(dat_prot, dat_RNA, method = \"spearman\") # you can use \"pearson\" or \"kendall\"\n\nprot_RNA_mcorr <- \n  as.data.frame(prot_RNA_mcorr) %>%\n  mutate(protein = row.names(.)) %>% \n  pivot_longer(cols = 1:(ncol(.)-1), names_to = \"RNA\", values_to = \"corr.coeff\") %>%\n  separate(col = \"RNA\", into = c(\"EnsemblID\", \"RNA\"), sep = \"_\")\ntoc()","metadata":{"execution":{"iopub.status.busy":"2022-11-26T19:55:23.985555Z","iopub.execute_input":"2022-11-26T19:55:23.988009Z","iopub.status.idle":"2022-11-26T19:57:25.217194Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# prepare data for ifraph. Out is a data frame with\n# proteins, RNA and correlation coefficients\nprep_igraph <- function(out) {\n    \n    nodes <- \n     rbind(data.frame(label = unique(out$protein),\n                      group = \"protein\"),\n           data.frame(label = unique(out$RNA),\n                      group = \"RNA\") ) %>%\n     mutate(id = 1:nrow(.))\n\n    per_prot <- out %>%  \n     select(protein, RNA, corr.coeff) %>%\n     mutate(width = corr.coeff)\n\n    edges <- per_prot %>% \n     left_join(nodes, by = c(\"protein\" = \"label\")) %>% \n     rename(from = id) %>% \n     left_join(nodes, by = c(\"RNA\" = \"label\")) %>% \n     rename(to = id) %>%\n     filter(group.x == \"protein\" & group.y == \"RNA\") %>%\n     select(from, to, width)\n    \n    out <- list(nodes, edges)\n    names(out) <- c(\"nodes\", \"edges\")\n    return(out)\n}","metadata":{"execution":{"iopub.status.busy":"2022-11-26T21:05:23.472993Z","iopub.execute_input":"2022-11-26T21:05:23.474912Z","iopub.status.idle":"2022-11-26T21:05:23.499314Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# 1. All proteins with all RNA: most positively correlated combinations","metadata":{}},{"cell_type":"code","source":"# top 100\nt100_corr <-\n  prot_RNA_mcorr[order(prot_RNA_mcorr$corr.coeff, decreasing = TRUE), ]\nt100_corr <- t100_corr[1:100, ]\n\n# take a look at first 10 combinations \ncat(\"The range of top 100 positive correlations:\",\n    round(min(t100_corr$corr.coeff), digits = 2), \"-\",\n    round(max(t100_corr$corr.coeff), digits = 2))","metadata":{"execution":{"iopub.status.busy":"2022-11-26T21:17:57.845471Z","iopub.execute_input":"2022-11-26T21:17:57.847273Z","iopub.status.idle":"2022-11-26T21:17:59.033478Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Now let's visualise top most pronounced protein-RNA correlations (corr. coefficient > 0.5) with network analysis.","metadata":{}},{"cell_type":"code","source":"cut_off <- 0.5\n\nt100_corr %>%\n  filter(corr.coeff > cut_off) %>%\n  group_by(protein) %>%\n  summarise(top_corr = paste(RNA, collapse = \",\\n\"))","metadata":{"execution":{"iopub.status.busy":"2022-11-26T21:18:08.812305Z","iopub.execute_input":"2022-11-26T21:18:08.814077Z","iopub.status.idle":"2022-11-26T21:18:08.860704Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"dat <- filter(prot_RNA_mcorr, corr.coeff > cut_off)\ndat_graph <- prep_igraph(dat)\n\nvisNetwork(dat_graph$nodes, dat_graph$edges) %>% \n  visIgraphLayout(layout = \"layout_with_fr\") %>%\n  visLegend()","metadata":{"execution":{"iopub.status.busy":"2022-11-26T21:18:13.969300Z","iopub.execute_input":"2022-11-26T21:18:13.971852Z","iopub.status.idle":"2022-11-26T21:18:14.274532Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The same for most interesting proteins","metadata":{}},{"cell_type":"code","source":"cut_off <- 0.5\n\nm_int <-\n  prot_RNA_mcorr %>%\n  filter(protein %in% c(\n  \"CD41\", \"CD36\",\"CD115\", \"CD71\", \"CD88\")) %>%\n  filter(corr.coeff > cut_off)\n\nm_int %>%\n  group_by(protein) %>%\n  summarise(top_corr = paste(RNA, collapse = \",\\n\"))","metadata":{"execution":{"iopub.status.busy":"2022-11-26T21:18:34.827743Z","iopub.execute_input":"2022-11-26T21:18:34.829421Z","iopub.status.idle":"2022-11-26T21:18:34.968882Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"dat_graph <- prep_igraph(m_int)\n\nvisNetwork(dat_graph$nodes, dat_graph$edges) %>% \n  visIgraphLayout(layout = \"layout_with_fr\") %>%\n  visLegend()","metadata":{"execution":{"iopub.status.busy":"2022-11-26T21:18:38.163657Z","iopub.execute_input":"2022-11-26T21:18:38.165374Z","iopub.status.idle":"2022-11-26T21:18:38.369208Z"},"trusted":true},"execution_count":null,"outputs":[]}]}