{"metadata":{"kernelspec":{"name":"ir","display_name":"R","language":"R"},"language_info":{"name":"R","codemirror_mode":"r","pygments_lexer":"r","mimetype":"text/x-r-source","file_extension":".r","version":"4.0.5"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"

RNA-RNA correlation analysis

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

What we do here:

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

Data

\n\n
\n
    \n\n
  1. RDS files with sparse matrices of normalised counts data for Open Problems - Multimodal Single-Cell Integration (CITEseq2022). The dataset for this competition comprises single-cell multiomics data (n = 70988 cells) collected from mobilized peripheral CD34+ hematopoietic stem and progenitor cells (HSPCs) isolated from four healthy human donors (normalised counts data).\nIn total, the CITEseq2022 train dataset contains 140 surface proteins and 22050 RNA form 3 days, 3 donors, 7 cell types. The gene names are given by {EnsemblID}_{GeneName} where EnsemblID refers to the Ensembl Gene ID and GeneName to the gene name (e.g. \"ENSG00000159840_ZYX\").\n
  2. CSV files with 112 protein-RNA pairs form the competition precalculated in this notebook.\n
  3. RDS files with sparse matrices of all RNA-RNA Spearman correlations (calculated here and stored in the dataset)\n
\n
\n

Output

\n \n
\n
    \n \n
  1. CSV file with 112 RNA vs all RNA Spearman correlations (112 RNA with corresponding protein product in the data set)
    2. \n
  2. CSV file with 112 RNA vs all RNA Spearman correlations for 5 different sample sizes
    3. \n
\n
\n
","metadata":{}},{"cell_type":"code","source":"# Libraries\nsuppressPackageStartupMessages({\n library(Matrix)\n library(dplyr) \n library(tidyr)\n library(ggplot2)\n library(tictoc)\n library(rhdf5)\n library(rstatix) #get_summary_stats\n library(patchwork) # arrange plots\n\n remotes::install_github(\"cysouw/qlcMatrix\")\n library(qlcMatrix)\n \n})\n\n# Function for figure size adjusment\nfig <- function(width, heigth) {\n options(repr.plot.width = width, repr.plot.height = heigth) }","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-03-27T16:55:33.129470Z","iopub.execute_input":"2023-03-27T16:55:33.131457Z","iopub.status.idle":"2023-03-27T16:55:33.415159Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Load data**","metadata":{}},{"cell_type":"code","source":"# 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)\ntoc()\ngc(full = TRUE)\n\nprot_RNA_corr <- read.csv(\"/kaggle/input/prot-rna-corr-analysis-for-2-citeseq-datasets/CITEseq22_112_prot_RNA_pairs.csv\",\n check.names = FALSE)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-03-27T11:06:35.719572Z","iopub.execute_input":"2023-03-27T11:06:35.757385Z","iopub.status.idle":"2023-03-27T11:07:21.577933Z"},"_kg_hide-output":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Functions**","metadata":{}},{"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,"execution":{"iopub.status.busy":"2023-03-27T11:07:21.580311Z","iopub.execute_input":"2023-03-27T11:07:21.581696Z","iopub.status.idle":"2023-03-27T11:07:21.599776Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Drop RNA with zero counts in all cells**","metadata":{}},{"cell_type":"code","source":"tic()\nmat_RNA <- mat_RNA[,unique(summary(mat_RNA)$j)]\ncat(ncol(mat_RNA), \"RNA have non-zero expression in at least one cell\\n\")\ncat(\"The resulting data consists of\", ncol(mat_RNA), \"RNA in\", nrow(mat_RNA), \"cells\\n\")\ntoc()","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-03-27T11:07:21.602075Z","iopub.execute_input":"2023-03-27T11:07:21.603413Z","iopub.status.idle":"2023-03-27T11:07:43.278648Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#

112 RNA vs all RNA correlations

","metadata":{}},{"cell_type":"markdown","source":"Here we calculate RNA-RNA Spearman correlations for 5 differend samples of cells. \nFirst, we subset 112 RNAs, that have corresponding protein product in the data set.","metadata":{}},{"cell_type":"code","source":"RNA_set <- filter(prot_RNA_corr, grepl(\"RNA\", name)) %>% select(name)\nRNA_set <- sub(\"\\\\).*\", \"\", sub(\".*\\\\(\", \"\", RNA_set$name)) %>% unique\nRNA_set <- mat_RNA[ , RNA_set]\n\ncat(\"The resulting data consists of\", ncol(RNA_set), \"RNA in\", nrow(RNA_set), \"cells\\n\")","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-03-27T11:07:43.281056Z","iopub.execute_input":"2023-03-27T11:07:43.282500Z","iopub.status.idle":"2023-03-27T11:07:43.967032Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"list_l_seq <- function(max_lenght) {\n \n s <- ceiling(seq(4000, max_lenght, length.out = max_lenght/10000))\n s <- s %>% sort %>% unique\n return(s)\n}","metadata":{"execution":{"iopub.status.busy":"2023-03-27T11:07:43.970768Z","iopub.execute_input":"2023-03-27T11:07:43.972519Z","iopub.status.idle":"2023-03-27T11:07:43.986955Z"},"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Next, we calculate correlations for 5 different samples of 4000, 10000, 20000, 50000, and 70988 cells. ","metadata":{}},{"cell_type":"code","source":"gc()\n\nRNA_RNA_corr_by_sample <- NULL\ncorr_per_ss <- function(N) {\n \n sampled_cells <- sample(1:nrow(mat_RNA), N)\n RNA_RNA_corr <- SparseSpearmanCor2(mat_RNA[sampled_cells, ], RNA_set[sampled_cells, ])\n colnames(RNA_RNA_corr) <- colnames(RNA_set) #%>% length\n rownames(RNA_RNA_corr) <- colnames(mat_RNA) #%>% length\n RNA_RNA_corr <- cbind(\"sample_size\" = N,\n RNA_RNA_corr)\n return(RNA_RNA_corr)\n}\n\nRNA_RNA_corr <- corr_per_ss(4000)\nRNA_RNA_corr_by_sample <- rbind(RNA_RNA_corr_by_sample, RNA_RNA_corr)\ngc()\nRNA_RNA_corr <- corr_per_ss(10000)\nRNA_RNA_corr_by_sample <- rbind(RNA_RNA_corr_by_sample, RNA_RNA_corr)\ngc()\nRNA_RNA_corr <- corr_per_ss(20000)\nRNA_RNA_corr_by_sample <- rbind(RNA_RNA_corr_by_sample, RNA_RNA_corr)\ngc()\nRNA_RNA_corr <- corr_per_ss(50000)\nRNA_RNA_corr_by_sample <- rbind(RNA_RNA_corr_by_sample, RNA_RNA_corr)\ngc()\nRNA_RNA_corr <- corr_per_ss(70988)\nRNA_RNA_corr_by_sample <- rbind(RNA_RNA_corr_by_sample, RNA_RNA_corr)\ngc()\n\n#RNA_RNA_corr_by_sample[, \"sample_size\"] %>% table","metadata":{"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Write 112RNA - allRNA Spearman correrelation matrix into csv file**","metadata":{}},{"cell_type":"code","source":"cat(\"The correlation matrix for\", ncol(RNA_RNA_corr), \"RNA vs all (\", nrow(RNA_RNA_corr), \") RNA with at least one non-zero entry\")\nhead(RNA_RNA_corr)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-03-27T12:00:05.806327Z","iopub.execute_input":"2023-03-27T12:00:05.808009Z","iopub.status.idle":"2023-03-27T12:00:05.846033Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"write.csv(RNA_RNA_corr[, -1],\"allRNA_112RNA_corr.csv\", row.names=TRUE)\nrm(RNA_RNA_corr) ; gc()","metadata":{"execution":{"iopub.status.busy":"2023-03-27T12:00:39.062411Z","iopub.execute_input":"2023-03-27T12:00:39.064174Z","iopub.status.idle":"2023-03-27T12:00:42.520464Z"},"_kg_hide-input":false,"_kg_hide-output":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#

Significant correlations

","metadata":{}},{"cell_type":"markdown","source":"We calculate p-values as in the [cor.test](http://www.rdocumentation.org/packages/stats/versions/3.6.2/topics/cor.test) function. For ***n >= 1290*** p-values are computed using the asymptotic t approximation. Therefore, we first estimate t statistic as:\n\n$$t = rho * \\sqrt{\\frac{(n-2)}{1 − rho^2}}$$\n\nwhere ***rho*** is a Spearman correlation coefficientm and ***n*** is the number of observations (70988 in this data set), and then use [pt function](http://www.rdocumentation.org/packages/stats/versions/3.6.2/topics/TDist) to get p-values. We use the Bonferroni correction and the \"BH\" method (Benjamini, Hochberg) for p-value correction. According to Stats::p.adjust [documentation](http://www.rdocumentation.org/packages/stats/versions/3.6.2/topics/p.adjust), BH method controls the false discovery rate, the expected proportion of false discoveries amongst the rejected hypotheses. The false discovery rate is a less stringent condition than the family-wise error rate, so this method is more powerful than the others.","metadata":{}},{"cell_type":"code","source":"tic()\ncor_stats <- NULL\nfor(ss in unique(RNA_RNA_corr_by_sample[, \"sample_size\"])) {\n \n tmp <- subset(RNA_RNA_corr_by_sample, \n RNA_RNA_corr_by_sample[, \"sample_size\"] == ss) %>%\n as.data.frame %>%\n mutate(RNA2 = rownames(.))\n \n for (g in colnames(tmp)[-c(1, ncol(tmp))]) {\n\n #subset correlations with all RNA\n cor <- select(tmp, c(\"RNA2\", all_of(g), \"sample_size\") )\n names(cor)[2] <- \"corr.coeff\"\n cor$RNA1 <- g\n\n # remove RNA correlations with itself\n cor <- filter(cor, RNA1 != RNA2)\n \n #Estimate T Statistic and calculate p-values\n cor$t_stats <- abs(cor[,2])* sqrt((ss-2)/(1-abs(cor[,2])^2))\n cor$p.value <- 2*pt(cor$t_stats, ss-2, lower.tail=FALSE)\n cor$BH_p.value <- p.adjust(cor$p.value, method = \"BH\") \n cor$Bonf_p.value <- p.adjust(cor$p.value, method = \"bonferroni\")\n cor_stats <- rbind(cor_stats, cor)\n }\n}\ntoc()\ncat(\"The resulting data frame with Spearman correlations between\", \n n_distinct(cor_stats$RNA1), \"RNA vs all other\", n_distinct(cor_stats$RNA2), \"RNA\" , \"and statistics for\",\n n_distinct(cor_stats$sample_size),\"different samples, including the whole dataset\")\ncor_stats <- cor_stats %>% select(\"sample_size\", \"RNA1\", \"RNA2\", \"corr.coeff\", everything()) \nhead(cor_stats) \n\ncor_stats %>% group_by(sample_size) %>% summarise(n = n())","metadata":{"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Write 112RNA - allRNA Spearman correrelations with stats for 6 samples into csv and rds files**","metadata":{}},{"cell_type":"code","source":"write.csv(cor_stats,\"allRNA_112RNA_corr_by6ss.csv\", row.names=TRUE)\nsaveRDS(cor_stats, \"allRNA_112RNA_corr_by6ss.RDS\")","metadata":{"_kg_hide-input":false,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#cor_stats <- read.csv(\"/kaggle/working/allRNA_112RNA_corr_by6ss.csv\")\ncat(\"Descriptive statistics of p-values with and without adjustment\")\ncor_stats %>% as.data.frame %>% \n group_by(sample_size) %>%\n #get_summary_stats(c(\"p.value\", \"BH_p.value\", \"Bonf_p.value\"), type = \"common\", digits = 5) %>%\n summarise(\"n_corr\" = sum(!is.na(corr.coeff)),\n \"Mean_p.value\" = mean(p.value, na.rm = TRUE),\n \"Mean_BH_p.value\" = mean(BH_p.value, na.rm = TRUE),\n \"Mean_Bonf_p.value\" = mean(Bonf_p.value, na.rm = TRUE))\ncor_stats <- na.omit(cor_stats)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-03-27T17:16:20.587702Z","iopub.execute_input":"2023-03-27T17:16:20.589631Z","iopub.status.idle":"2023-03-27T17:16:47.984138Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"cor_sign_p <- NULL\n\nfor(s in unique(cor_stats$sample_size)) {\n \n for(rna in unique(cor_stats$RNA1)) {\n \n sbst <- cor_stats[cor_stats$sample_size == s & cor_stats$RNA1 == rna, ]\n\n if(nrow(sbst) > 0) {\n\n p <- filter(sbst, p.value < 0.05)\n p <- filter(p, p.value == max(p.value))\n p$Adjustment <- \"none\"\n\n cor_sign_p <- rbind(cor_sign_p, p)\n\n p <- filter(sbst, BH_p.value < 0.05)\n p <- filter(p, BH_p.value == max(BH_p.value))\n p$Adjustment <- \"BH\"\n\n cor_sign_p <- rbind(cor_sign_p, p)\n\n p <- filter(sbst, Bonf_p.value < 0.05)\n p <- filter(p, Bonf_p.value == max(Bonf_p.value))\n p$Adjustment <- \"Bonferroni\"\n\n cor_sign_p <- rbind(cor_sign_p, p) \n } \n } \n}\n\n\ncat(\"Minimal value of the absolute correlation coefficient\\nor significant correlations with and without p-value adjustment \")\ncor_sign_p %>% group_by(sample_size, Adjustment) %>%\n summarise(corr.coeff = min(abs(corr.coeff)), .groups = \"keep\") %>%\n pivot_wider(names_from = \"Adjustment\", values_from = \"corr.coeff\")","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-03-27T17:16:47.986688Z","iopub.execute_input":"2023-03-27T17:16:47.988070Z","iopub.status.idle":"2023-03-27T17:16:48.806984Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Unsignificant cutoff with Benjamini & Hochberg p-value correction**","metadata":{}},{"cell_type":"code","source":"fig(25,15)\nn = 20000\n\nsbst <- mutate(cor_stats, sample_size = factor(sample_size))\nsbst <- filter(sbst, corr.coeff < 0.1 & corr.coeff > -0.1) \ncutoff <- filter(cor_stats, BH_p.value < 0.05)\ncutoff <- min(abs(cutoff$corr.coeff))\n\nsbst <- sbst[sample(1:nrow(sbst), n), ]\n\np1 <- ggplot(sbst, aes(x = BH_p.value, y = corr.coeff, color = sample_size)) +\n geom_point() +\n geom_hline(yintercept = c(cutoff, -1*cutoff),\n linetype=\"dashed\") +\n annotate(\"text\", x = 0.3, y = cutoff,\n label = paste0(\"Cutoff with\\nBenjamini & Hochberg p-value correction:\\nfor p < 0.05, min(R) =\", round(cutoff, 4)),\n size = 8, vjust = - 3 ) +\n annotate(\"rect\", xmin = 0.5, xmax = Inf, ymin = -1*cutoff, ymax = cutoff, alpha = .2) +\n geom_vline(xintercept = 0.5, linetype=\"dashed\") +\n theme_bw(base_size = 20)\n\np2 <- ggplot(sbst, aes(x = BH_p.value, y = corr.coeff, color = sample_size)) +\n geom_point() +\n geom_hline(yintercept = c(cutoff, -1*cutoff),\n linetype=\"dashed\") +\n annotate(\"rect\", xmin = 0.05, xmax = Inf, ymin = -1*cutoff, ymax = cutoff, alpha = .2) +\n geom_vline(xintercept = 0.05, linetype=\"dashed\") +\n xlim(c(0, 0.055)) +\n ylim(c(-0.1, 0.1)) +\n ggtitle(\"Zoom in\") +\n theme_bw(base_size = 20)\n\np1 + p2 +\n plot_layout(widths = c(3, 1)) +\n plot_annotation(\n title = paste0(\"112 RNA vs all RNA in the CITEseq22 dataset, sample = \", n)) & \ntheme(text = element_text(size = 22) )","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-03-27T17:36:20.151613Z","iopub.execute_input":"2023-03-27T17:36:20.153260Z","iopub.status.idle":"2023-03-27T17:36:31.253137Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Unsignificant cutoff with Bonferroni p-value correction**","metadata":{}},{"cell_type":"code","source":"fig(25,10)\nn = 20000\n\n#sbst <- cor_stats\n#sbst <- filter(cor_stats, corr.coeff < 0.1 & corr.coeff > -0.1) \ncutoff <- filter(cor_stats, Bonf_p.value < 0.05)\ncutoff <- min(abs(cutoff$corr.coeff))\n\nsbst <- sbst[sample(1:nrow(sbst), n), ]\n\np1 <- ggplot(sbst, aes(x = Bonf_p.value, y = corr.coeff, color = sample_size)) +\n geom_point() +\n geom_hline(yintercept = c(cutoff, -1*cutoff),\n linetype=\"dashed\") +\n annotate(\"text\", x = 0.3, y = cutoff,\n label = paste0(\"Cutoff with Bonferroni p-value correction:\\nfor p < 0.05, min(R) =\", round(cutoff, 4)),\n size = 8, vjust = -3 ) +\n annotate(\"rect\", xmin = 0.05, xmax = Inf, ymin = -1*cutoff, ymax = cutoff, alpha = .2) +\n geom_vline(xintercept = 0.05, linetype=\"dashed\") +\n theme_bw(base_size = 20)\n\np2 <- ggplot(sbst, aes(x = Bonf_p.value, y = corr.coeff, color = sample_size)) +\n geom_point() +\n geom_hline(yintercept = c(cutoff, -1*cutoff),\n linetype=\"dashed\") +\n annotate(\"rect\", xmin = 0.05, xmax = Inf, ymin = -1*cutoff, ymax = cutoff, alpha = .2) +\n geom_vline(xintercept = 0.05, linetype=\"dashed\") +\n xlim(c(0, 0.055)) +\n ggtitle(\"Zoom in\") +\n theme_bw(base_size = 20)\n\np1 + p2 +\n plot_layout(widths = c(3, 1)) +\n plot_annotation(\n title = paste0(\"112 RNA vs all RNA in the CITEseq22 dataset, sample = \", n)) & \ntheme(text = element_text(size = 22) )","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-03-27T17:36:42.169157Z","iopub.execute_input":"2023-03-27T17:36:42.171927Z","iopub.status.idle":"2023-03-27T17:36:44.896082Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#ls()\nrm('cor','cor_sign_p','cor_stats','cutoff','g','n','p','p1','p2',\n 'prot_RNA_corr','rna','RNA_RNA_corr','RNA_set','sbst', 'tmp') ; gc()","metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"execution":{"iopub.status.busy":"2023-03-27T17:37:25.876686Z","iopub.execute_input":"2023-03-27T17:37:25.878436Z","iopub.status.idle":"2023-03-27T17:37:32.742046Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#

Validate all RNA - all RNA correlations calculation

","metadata":{}},{"cell_type":"code","source":"tic() #\n# Load Metadata (day, donor, cell type, and technology)\nRNA_RNA_corr <- readRDS('/kaggle/input/proteinrna-vs-rna-spearman-correlation-data/allRNA_RNA_corr.RDS')\nall_corr <- as.matrix(RNA_RNA_corr)\n\ntoc()\n\nhead(all_corr, n = 10)\n\n# #save \n# file_name <- \"allRNA_RNA_corr\"\n# saveRDS(RNA_RNA_corr, file = paste0(file_name, \".RDS\"))\n\n# h5write(colnames(RNA_RNA_corr), paste0(file_name, \".h5\"), name=\"col_and_row_names\")\n# h5write(RNA_RNA_corr@i, paste0(file_name, \".h5\"), name=\"row_ind_i\")\n# h5write(RNA_RNA_corr@p, paste0(file_name, \".h5\"), name=\"col_ind_p\")\n# h5write(RNA_RNA_corr@x, paste0(file_name, \".h5\"), name=\"corr_coeff_x\")\n\n# cat(\"The resuling data structure\")\n# h5ls(paste0(file_name, \".h5\"))\n# cat(\"\\n\\n\\nThe tables are stored in the output\")","metadata":{"_kg_hide-output":true,"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Comparison of Spearman correlation coefficients calculated on dense and sparse matrices**","metadata":{}},{"cell_type":"code","source":"test_gene <- \"ENSG00000135218_CD36\"\n\ntst <- cor(mat_RNA[,test_gene], as.matrix(mat_RNA[, 1:1000]), method = \"spearman\")\ntst <- pivot_longer(as_tibble(tst),cols = 1:ncol(tst), names_to = \"RNA\", values_to = \"dense_corr\")\n\ntst$sparse_corr <- NA\nfor (rna in tst$RNA) {\n tst[tst$RNA == rna, ]$sparse_corr <-\n all_corr[rna, test_gene] \n}\n\ntst$abs_difference <- abs(tst$dense_corr - tst$dense_corr)\n\n#tst$sparse_corr <- all_corr[1:1000, test_gene] # to easy for me :)\ncat(\"Comparison of Spearman correlation coefficients for one RNA\",test_gene,\"\ncalculated in the usual way on a dense matrix (dense_corr)\nand those calculated on a sparse matrix using the SparseSpearmanCor2 function (sparse_corr).\")\n\ntst[tst$sparse_corr != 0, ] %>% head","metadata":{"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig(25,10)\ncat(\"The Spearman correlation coefficients for one RNA (ENSG00000135218_CD36)\ncalculated on a dense (x axis) and sparse matrices (y axis) are the same.\nExcept that we have replaced all values between -0.05 and 0.05 with zeros.\")\ntst[tst$RNA == test_gene, ]$dense_corr <- 0\n\nggplot(tst, aes(x = dense_corr, y = sparse_corr)) +\n geom_point() +\n theme_bw(base_size = 22) +\n theme(aspect.ratio = 1)","metadata":{"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# w <- ncol(all_corr)\n# all_corr <- all_corr[, colSums(all_corr) != 0]\n# cat(\"There are\",ncol(all_corr),\"/\", w, \"RNA with at least one non-zero value across\",\n# nrow(all_corr), \"cells.\")\n# rm(w)","metadata":{"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# # reshape to long format\n# all_corr <- pivot_longer(all_corr[], cols = 1:(ncol(all_corr)-1),\n# names_to = \"RNA2\", values_to = \"corr.coeff\") \n\n# #remove lower part of a correlation matrix\n# all_corr <- filter(all_corr, corr.coeff != 0)\n\n# cat(\"The head of the data frame with\", nrow(all_corr), \"protein-protein pairs\",\n# \"\\nr (correlation coefficient) range:\",\n# round(min(all_corr$corr.coeff, na.rm = TRUE),2), \"to\",\n# round(max(all_corr$corr.coeff, na.rm = TRUE), 2), \"\\nmean =\",\n# round(mean(all_corr$corr.coeff, na.rm = TRUE), 4) )\n\n# head(all_corr)","metadata":{"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# zeros <- apply(all_corr_RNA[2:ncol(all_corr_RNA)],2, function(x) length(x[x == 0]))\n# cat(\"Of\", length(zeros), \"RNA\", length(zeros[zeros == max(zeros)]),\n# \"have zero correlarion with all RNA in the dataset\\nthe remaining RNA have from\",\n# min(zeros), \"to\", max(zeros[zeros != max(zeros)]), \"zero correlarions with all RNA.\" )","metadata":{"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The tables are stored in the output.","metadata":{}}]}