{"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- subset 10% cells for fast explaratory analysis\n- use methods: log normalization, centered log ratio transformation, log1p, Pearson residuals\n- comparison of gene expression\n\n**Data:** RDS files with sparse matrices of [normalised counts data](http://www.kaggle.com/datasets/stautxie/sparse-measurement-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. The train data we use here 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).****","metadata":{}},{"cell_type":"code","source":"suppressMessages(library(Matrix))\nsuppressMessages(library(Seurat))\nsuppressMessages(library(tictoc)) \nsuppressMessages(library(tidyverse)) \nsuppressMessages(library(ggrepel)) \n\n\nfig <- function(width, heigth){\n  options(repr.plot.width = width, repr.plot.height = heigth)\n}","metadata":{"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Load data and create the Seurat object**","metadata":{}},{"cell_type":"markdown","source":"**Proteins (Targets dataset):** for the surface protein levels, each row corresponds to a cell (e.g. \"45006fe3e4c8\") and each column to a protein (e.g. \"CD86\").\n\n**RNA (Inputs dataset):** For the RNA counts, each row corresponds to a cell (e.g. \"45006fe3e4c8\") and each column to a gene. The column format for a gene is given by {EnsemblID}_{GeneName} where EnsemblID refers to the Ensembl Gene ID and GeneName to the gene name (e.g. \"ENSG00000159840_ZYX\").\n\n**Metadata:** Donor and cell types. The train data consists of both gene expression (RNA) and surface protein data for days 2,3,4 for donors 1-3 (donor IDs: 32606,13176, and 31800), the public test data consists of RNA for days 2,3,4 for donor 4 (donor ID: 27678) and the private test data consists data from day 7 from all donors.","metadata":{}},{"cell_type":"code","source":"# To randomly select x% of cells to save memory footprint\nset.seed(1)\nselected_prop <- 0.1\nselected_cells <- sample(1:70988,\n                        floor(selected_prop * 70988),\n                        replace = FALSE)","metadata":{"_uuid":"051d70d956493feee0c6d64651c6a088724dca2a","_execution_state":"idle","execution":{"iopub.status.busy":"2023-01-16T18:24:37.509909Z","iopub.execute_input":"2023-01-16T18:24:37.511393Z","iopub.status.idle":"2023-01-16T18:24:37.528040Z"},"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Load raw RNA data\npath <- \"/kaggle/input/sparse-raw-counts-data-open-problems-multimodal/citeseq/sp_train_cite_inputs_raw.rds\"\nmat_RNA <- readRDS(path)\nmat_RNA <- mat_RNA[selected_cells, ] #use xx% of cells for explaratory analysis\nmat_RNA <- t(mat_RNA)\n\ncat(nrow(mat_RNA), \"rows with RNA names and\", ncol(mat_RNA), \"columns with cell IDs\")\n\n# Load and filter metadata\nmetadata <- read.csv('../input/open-problems-multimodal/metadata.csv',row.names = 1)\n#summary(mutate_all(metadata, as.factor))\n\n# Filter citeseq, donor, day\nmetadata <- filter(metadata, technology == \"citeseq\" &\n                   donor %in% c(\"13176\", \"31800\", \"32606\") &\n                  day %in% c(\"2\",\"3\",\"4\"))\nmetadata <- filter(metadata, row.names(metadata) %in% colnames(mat_RNA)) #use xx% of cells for explaratory analysis\ncat(nrow(metadata), \"rows with cell IDs\")\nsummary(mutate_all(metadata, as.factor))\n#head(metadata)","metadata":{"execution":{"iopub.status.busy":"2023-01-16T18:24:40.398395Z","iopub.execute_input":"2023-01-16T18:24:40.400152Z","iopub.status.idle":"2023-01-16T18:25:18.898226Z"},"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Initialize the Seurat object with the raw (non-normalized data)\ntic()\nseur_RNA <- CreateSeuratObject(counts = mat_RNA, project = \"MmSCel\",\n                               min.cells = 100, min.features = 100, \n                               meta.data = metadata)\ntoc()\nseur_RNA\ncat(\"meta data:\\n\")\nstr(seur_RNA@meta.data)\n\ncat(\"\\n\\nThe top 6 rows of the data:\")\nhead(as.data.frame(seur_RNA@assays$RNA@data))","metadata":{"execution":{"iopub.status.busy":"2023-01-16T18:25:26.027894Z","iopub.execute_input":"2023-01-16T18:25:26.029378Z","iopub.status.idle":"2023-01-16T18:25:36.313549Z"},"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#clean up\nfor (thing in ls()) { message(thing); print(object.size(get(thing)), units='auto') }\nrm(mat_RNA, metadata,path,selected_prop)\ngc()","metadata":{"_kg_hide-output":true,"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Log normalization","metadata":{}},{"cell_type":"markdown","source":"By default, we employ a global-scaling normalization method \"LogNormalize\" that normalizes the feature expression measurements for each cell by the total expression, multiplies this by a scale factor (10,000 by default), and log-transforms the result. Normalized values are stored in seur_RNA[[\"RNA\"]]@data. \n\nIn the MmSCel competition, RNA counts were [library-size normalized](http://scanpy.readthedocs.io/en/stable/generated/scanpy.pp.normalize_per_cell.html) (each cell normalized by total counts over all genes, so that every cell has the same total count after normalization) and [log1p transformed](http://scanpy.readthedocs.io/en/stable/generated/scanpy.pp.log1p.html) (X = log ⁡ ( X + 1 ))","metadata":{}},{"cell_type":"code","source":"# this method divides the feature expression measures for each cell by the total expression of the cell\n# then then result is multiplied by 10'000 (scale.factor) and log transformed (normalization.method)\n\nlognorm <- NormalizeData(seur_RNA, normalization.method = \"LogNormalize\", scale.factor = 10000) # default settings are shown for clarity\n\ncat(\"\\n\\nThe top 6 rows of the data after Log normalization:\")\nhead(as.data.frame(lognorm@assays$RNA@data))","metadata":{"execution":{"iopub.status.busy":"2023-01-16T18:25:49.277633Z","iopub.execute_input":"2023-01-16T18:25:49.279033Z","iopub.status.idle":"2023-01-16T18:25:57.675292Z"},"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Centered log ratio transformation","metadata":{}},{"cell_type":"markdown","source":"As [here](http://www.bioinformatics.babraham.ac.uk/training/10XRNASeq/seurat_workflow.html#Normalisation,_Selection_and_Scaling), simplistic normalisation of raw RNA counts doesn’t actually normalise the quantitative data very well because it’s so biased by the proportion of zero values in the dataset. \n\nWe can try the normalisation again, this time using a centered log ratio transformation - more similar to the sort of size factor based normalisation which is used for many RNA-Seq experiments. The margin=2 option means that it normalises per cell instead of per gene.","metadata":{}},{"cell_type":"code","source":"# The centered log ratio transformation\nclr <- NormalizeData(seur_RNA, normalization.method = \"CLR\", margin = 2)\n\ncat(\"\\n\\nThe top 6 rows of the data after centered log ratio transformation:\")\nhead(as.data.frame(clr@assays$RNA@data))","metadata":{"execution":{"iopub.status.busy":"2023-01-16T18:26:03.977368Z","iopub.execute_input":"2023-01-16T18:26:03.978913Z","iopub.status.idle":"2023-01-16T18:26:25.305904Z"},"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Pearson residuals","metadata":{}},{"cell_type":"markdown","source":"Another option is to compute Pearson residuals from negative binomial regression, recommend e.g. by Hafemeister and Satija (Genome Biol 20:296, 2019, [vignette](http://satijalab.org/seurat/articles/sctransform_vignette.html)). This procedure omits the need for heuristic steps including pseudocount addition or log-transformation and in application for common downstream analytical tasks such as variable gene selection, dimensional reduction, and differential expression outperform other methods (discussed [here](http://www.kaggle.com/datasets/alexandervc/research-project-01-around-multimodal-singlecell/discussion/378460)).\n\nFor this we use SCTransform function from [sctransform packadge](http://github.com/saketkc/sctransform) ([v2 regularization vignette](http://satijalab.org/seurat/articles/sctransform_v2_vignette.html#identify-differential-expressed-genes-across-conditions-1)). \n- Note that this single command replaces NormalizeData(), ScaleData(), and FindVariableFeatures().\n- Transformed data will be available in the SCT assay, which is set as the default after running sctransform\n- During normalization, we can also remove confounding sources of variation, for example, mitochondrial mapping percentage","metadata":{}},{"cell_type":"code","source":"tic()\ndevtools::install_github(\"satijalab/seurat\", ref = \"develop\")\ninstall.packages(\"BiocManager\")\nBiocManager::install(\"glmGamPoi\")\ntoc()\ngc()","metadata":{"_kg_hide-output":true,"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"tic()\npears_resid <- SCTransform(seur_RNA, vst.flavor = \"v2\",\n                           method=\"glmGamPoi\", verbose = FALSE,\n                          return.only.var.genes = FALSE)\ntoc()\n# SCTransform returns a Seurat object with a new assay (named SCT by default) with\n# counts being (corrected) counts, data being log1p(counts), scale.data being pearson residuals; \n# sctransform::vst intermediate results are saved in misc slot of the new assay.\n\ncat(\"\\n\\nThe top 6 rows of the data after log1p transformation:\")\nhead(as.data.frame(pears_resid@assays$SCT@data))\ncat(\"\\n\\nThe top 6 rows of the data after normalisation using Pearson residuals:\")\nhead(as.data.frame(pears_resid@assays$SCT@scale.data))","metadata":{"execution":{"iopub.status.busy":"2023-01-16T21:55:14.019439Z","iopub.execute_input":"2023-01-16T21:55:14.021257Z","iopub.status.idle":"2023-01-16T21:56:32.639349Z"},"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Compare results","metadata":{}},{"cell_type":"markdown","source":"## Most highly expressed genes","metadata":{}},{"cell_type":"code","source":"cat(\"Top 10 most highly expressed genes (maximal mean expressions per cell)\")\ndata.frame(\n    \"RawCounts\" = names(apply(seur_RNA@assays$RNA@data,1,mean) %>% sort(decreasing = TRUE) %>% head(n=10)),\n    \"LogNorm\" = names(apply(lognorm@assays$RNA@data,1,mean) %>% sort(decreasing = TRUE) %>% head(n=10)),\n    \"CenteredLogRatio\" = names(apply(clr@assays$RNA@data,1,mean) %>% sort(decreasing = TRUE) %>% head(n=10)),\n    \"Log1p\" = names(apply(pears_resid@assays$SCT@data,1,mean) %>% sort(decreasing = TRUE) %>% head(n=10)),\n    \"PearsonResiduals\" = names(apply(pears_resid@assays$SCT@scale.data,1,mean) %>% sort(decreasing = TRUE) %>% head(n=10)))","metadata":{"execution":{"iopub.status.busy":"2023-01-16T21:56:46.054830Z","iopub.execute_input":"2023-01-16T21:56:46.056583Z","iopub.status.idle":"2023-01-16T21:57:15.665840Z"},"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"With Pearson residuals we see MALAT1 - a transcript highly expressed in the nucleus everywhere apart from platelets. The high ammounts of MALAT1 may indicate dead/dying cells ([link](http://kb.10xgenomics.com/hc/en-us/articles/360004729092-Why-do-I-see-high-levels-of-Malat1-in-my-gene-expression-data-)). This transcript is usualy one of the most highly expressed genes in scRNA-seq data.","metadata":{}},{"cell_type":"markdown","source":"## Gene expression histograms","metadata":{}},{"cell_type":"code","source":"fig(22,7)\n# set seed and put two plots in one figure\n#set.seed(123) - if you need to sample data\npar(mfrow=c(1,5),\n    mar=c(8,8,4,2),\n    mgp=c(5,1,0), cex.lab=2, cex.axis=2, cex.main=2, cex.sub=2\n   )\ncat(\"Gene expression by different normalization methods (zeros are not included)\")\n\n# original expression distribution\nRaw = as.vector(seur_RNA[['RNA']]@counts)# %>% sample(100)\nRaw = Raw[Raw != 0]\nhist(Raw, \n    xlab = \"Raw data\\n- RNA counts\")\n\n# expression distribution after log normalization\nlogNorm = as.vector(lognorm[['RNA']]@data)# %>% sample(100)\nlogNorm = logNorm[logNorm != 0]\nhist(logNorm, \n    xlab = \"NormalizeData fun,\\nnormalization.method = LogNormalize\")\n\n# expression distribution after centered log ratio transformation\nCLR = as.vector(clr[['RNA']]@data)# %>% sample(100)\nCLR = CLR[CLR != 0]\nhist(CLR, \n    xlab = \"NormalizeData fun,\\nnormalization.method = CLR\")\n\n# Log1p by SCTransform\nSCTLog1p = as.vector(pears_resid[['SCT']]@data)# %>% sample(100)\nSCTLog1p = SCTLog1p[SCTLog1p != 0]\nhist(SCTLog1p, \n    xlab = \"SCTransform fun,\\nlog1p\")\n\n# pearson residuals by SCTransform\nSCT_PR = as.vector(pears_resid[['SCT']]@scale.data)# %>% sample(100)\nSCT_PR = SCT_PR[SCT_PR != 0]\nhist(SCT_PR, \n    xlab = \"SCTransform fun,\\nPearson residuals\")","metadata":{"execution":{"iopub.status.busy":"2023-01-16T21:57:42.825820Z","iopub.execute_input":"2023-01-16T21:57:42.827262Z","iopub.status.idle":"2023-01-16T21:57:57.724681Z"},"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Distributions of expression values","metadata":{}},{"cell_type":"code","source":"fig(20,7)\ncat(\"Distributions of the first 100 cells' expression values\")\nplot_dens <- function(dat) {    \n    \n    plt <- as_tibble(dat[,1:100]) %>%\n          pivot_longer(\n            cols = everything(),\n            names_to = \"cell\",\n            values_to = \"Expression\" ) %>%\n          ggplot(aes(x = Expression, group = cell)) +\n          geom_density() +\n          coord_cartesian(ylim = c(0, 5), xlim = c(-1, 1)) +\n          theme_minimal(base_size = 24) \n\n    return(plt)    \n}\n\nplot_dens(seur_RNA[['RNA']]@counts) + ggtitle(\"Raw data - RNA counts\")\nplot_dens(lognorm[['RNA']]@data) + ggtitle(\"NormalizeData fun, normalization.method = LogNormalize\")\nplot_dens(clr[['RNA']]@data) + ggtitle(\"NormalizeData fun, normalization.method = CLR\")\nplot_dens(pears_resid[['SCT']]@data) + ggtitle(\"SCTransform fun, log1p\")\nplot_dens(pears_resid[['SCT']]@scale.data) + ggtitle(\"SCTransform fun, Pearson residuals\")","metadata":{"execution":{"iopub.status.busy":"2023-01-16T22:23:28.763896Z","iopub.execute_input":"2023-01-16T22:23:28.766061Z","iopub.status.idle":"2023-01-16T22:23:39.070279Z"},"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Variance of Pearson residuals vs mean expresseion of genes","metadata":{}},{"cell_type":"code","source":"# the plot was constructed using the code for fig.4 in Hafemeister & Satija 2019, https://doi.org/10.1186%2Fs13059-019-1874-1\nfig(20,10)\nplt <- data.frame(var = apply(seur_RNA@assays$RNA@data,1,var),\n                  mean = apply(seur_RNA@assays$RNA@data,1,mean),\n                  gene = unlist(seur_RNA@assays$RNA@data@Dimnames[1]))\n\ntop_var_genes <- plt %>% filter(rank(-var) <= 15) %>%\n  mutate(gene = gsub(\".*\\\\-\",\"\", gene))\n\nggplot(plt, aes(mean, var)) + \n  geom_point(shape=16, alpha=0.5, size=1) +\n  geom_point(data = top_var_genes, shape=16, size=0.9, color='red') +\n  scale_y_continuous(trans='sqrt') +\n  scale_x_continuous(trans='log10', breaks = c(0.001, 0.01, 0.1, 1, 10, 100), labels = MASS::rational) +\n  xlab('Gene mean') + ylab('Variance') +\n  theme_minimal(base_size = 24) +\n  geom_text_repel(data = top_var_genes, aes(label = gene), size = 7, color='red', max.overlaps = 15) +\n  ggtitle(\"Raw data - RNA counts\")","metadata":{"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# the plot was constructed using the code for fig.4 in Hafemeister & Satija 2019, https://doi.org/10.1186%2Fs13059-019-1874-1\nfig(20,10)\nplt <- data.frame(var = apply(clr@assays$RNA@data,1,var),\n                  mean = apply(clr@assays$RNA@data,1,mean),\n                  gene = rownames(clr@assays$RNA@data))\n\ntop_var_genes <- plt %>% filter(rank(-var) <= 15) %>%\n  mutate(gene = gsub(\".*\\\\-\",\"\", gene))\n\nggplot(plt, aes(mean, var)) + \n  geom_point(shape=16, alpha=0.5, size=1) +\n  geom_point(data = top_var_genes, shape=16, size=0.9, color='red') +\n  scale_y_continuous(trans='sqrt') +\n  scale_x_continuous(trans='log10', breaks = c(0.001, 0.01, 0.1, 1, 10, 100), labels = MASS::rational) +\n  xlab('Gene mean') + ylab('Variance') +\n  theme_minimal(base_size = 24) +\n  geom_text_repel(data = top_var_genes, aes(label = gene), size = 7, color='red', max.overlaps = 15) +\n  ggtitle(\"NormalizeData fun, normalization.method = CLR\")","metadata":{"execution":{"iopub.status.busy":"2023-01-16T22:58:42.202956Z","iopub.execute_input":"2023-01-16T22:58:42.204580Z","iopub.status.idle":"2023-01-16T22:58:55.269402Z"},"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# the plot was constructed using the code for fig.4 in Hafemeister & Satija 2019, https://doi.org/10.1186%2Fs13059-019-1874-1\nfig(20,10)\nplt <- data.frame(pr.var = pears_resid@assays$SCT@SCTModel.list$model1@feature.attributes$residual_variance,\n                  mean = pears_resid@assays$SCT@SCTModel.list$model1@feature.attributes$gmean,\n                  gene = pears_resid@assays$SCT@data@Dimnames[1] %>% unlist)\n\ntop_var_genes <- plt %>% filter(rank(-pr.var) <= 15) %>%\n  mutate(gene = gsub(\".*\\\\-\",\"\", gene))\n\nggplot(plt, aes(mean, pr.var)) + \n  geom_point(shape=16, alpha=0.5, size=0.8) +\n  geom_point(data = top_var_genes, shape=16, size=0.9, color='red') +\n  scale_y_continuous(trans='sqrt') +\n  scale_x_continuous(trans='log10', breaks = c(0.001, 0.01, 0.1, 1, 10, 100), labels = MASS::rational) +\n  xlab('Gene mean') + ylab('Residual variance') +\n  theme_minimal(base_size = 24) +\n  geom_text_repel(data = top_var_genes, aes(label = gene), size = 7, color='red') +\n  ggtitle(\"SCTransform fun, Pearson residuals\")","metadata":{"execution":{"iopub.status.busy":"2023-01-16T22:58:29.831405Z","iopub.execute_input":"2023-01-16T22:58:29.833182Z","iopub.status.idle":"2023-01-16T22:58:30.720598Z"},"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"As in the article [Hafemeister & Satija 2019](http://doi.org/10.1186%2Fs13059-019-1874-1), who suggested usind regularized negative binomial regression, the variance of Pearson residuals is independent of gene abundance, demonstrating that the GLM has successfully captured the mean-variance relationship inherent in the data.","metadata":{}}]}