{"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 (random 40% of cells):\n\n- Normalisation using Pearson residuals from negative binomial regression with [sctransform packadge](http://satijalab.org/seurat/articles/sctransform_v2_vignette.html#identify-differential-expressed-genes-across-conditions-1) (reasons are discussed [here](http://www.kaggle.com/datasets/alexandervc/research-project-01-around-multimodal-singlecell/discussion/378460))\n- PCA with [Seurat](http://satijalab.org/seurat/)\n- Plot the protein of interest vs two first principal components (3D interactive scatterplot with [plotly](http://plotly.com/r/))\n- Plot PCA_2 vs top 100 RNA correlated with the protein of interest (csv with all correlations are generated [in this notebook](http://www.kaggle.com/code/antoninadolgorukova/mmscel-protein-analysis))\n\n[Referense notebook](http://www.kaggle.com/code/reminho/basic-scrnaseq-tutorial), [Seurat package tutorial](https://satijalab.org/seurat/archive/v3.2/pbmc3k_tutorial.html)\n\n**Data:** RDS files with sparse matrices of [raw counts](https://www.kaggle.com/datasets/antoninadolgorukova/sparse-raw-counts-data-open-problems-multimodal) for [Open Problems - Multimodal Single-Cell Integration](http://http//www.kaggle.com/competitions/open-problems-multimodal) (MmSCel). 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":"suppressMessages(library(Matrix))\n#devtools::install_github(\"satijalab/seurat\", ref = \"develop\")\nsuppressMessages(library(Seurat))\nsuppressMessages(library(tictoc)) \nsuppressMessages(library(tidyverse))\n\n# visualisation packages\nsuppressMessages(library(ggrepel))\nsuppressMessages(library(ggplot2))\nsuppressMessages(library(patchwork))\nsuppressMessages(library(ggpubr)) #ggscatter \nsuppressMessages(library(ggExtra)) #ggMarginal\nsuppressMessages(library(gridExtra)) # side by side plots\nsuppressMessages(library(plotly)) #3D plot\n\nfig <- function(width, heigth){\n  options(repr.plot.width = width, repr.plot.height = heigth)\n}","metadata":{"_uuid":"051d70d956493feee0c6d64651c6a088724dca2a","_execution_state":"idle","_kg_hide-output":true,"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-01-24T06:59:34.426965Z","iopub.execute_input":"2023-01-24T06:59:34.429484Z","iopub.status.idle":"2023-01-24T06:59:43.377412Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# 1. 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":"markdown","source":"To create a Seurat object we need either a matrix-like object with unnormalized data or an Assay-derived object with number of molecules for each feature (i.e. gene) in rows that are detected in each cell (in columns). Thus, we load and transpose the MmSCel data.","metadata":{}},{"cell_type":"code","source":"# To randomly select x% of cells to save memory footprint\nset.seed(1)\nselected_prop <- 0.4\nselected_cells <- sample(1:70988,\n                        floor(selected_prop * 70988),\n                        replace = FALSE)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-01-24T06:59:43.381139Z","iopub.execute_input":"2023-01-24T06:59:43.432888Z","iopub.status.idle":"2023-01-24T06:59:43.461939Z"},"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, ]\nmat_RNA <- t(mat_RNA)\n\ncat(nrow(mat_RNA), \"rows with RNA names and\", ncol(mat_RNA), \"columns with cell IDs\")\nmat_RNA[c(\"ENSG00000135218_CD36\", \"ENSG00000029534_ANK1\", \"ENSG00000160789_LMNA\"), 1:30]","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-01-24T06:59:43.465766Z","iopub.execute_input":"2023-01-24T06:59:43.468008Z","iopub.status.idle":"2023-01-24T07:00:47.672934Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Additional cell-level metadata are stored in a data.frame where the rows are cell names and the columns are additional metadata fields. Row names in the metadata have to match the column names of the counts matrix.","metadata":{}},{"cell_type":"code","source":"# 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":{"_kg_hide-output":false,"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-01-24T07:00:47.676205Z","iopub.execute_input":"2023-01-24T07:00:47.678382Z","iopub.status.idle":"2023-01-24T07:00:49.067233Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Now we create the Seurat object with metadata. As in [this notebook](https://www.kaggle.com/code/reminho/basic-scrnaseq-tutorial) we filtered for genes that are found in at least 100 cells (min.cells = 100) and a single cell should express at least 100 of our target genes (min.features = 100).","metadata":{}},{"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)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-01-24T07:00:49.070683Z","iopub.execute_input":"2023-01-24T07:00:49.072973Z","iopub.status.idle":"2023-01-24T07:01:18.082195Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"rm(mat_RNA, metadata,path,selected_prop)\ngc()\nfor (thing in ls()) { message(thing); print(object.size(get(thing)), units='auto') }","metadata":{"_kg_hide-output":true,"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-01-24T07:01:18.084280Z","iopub.execute_input":"2023-01-24T07:01:18.085513Z","iopub.status.idle":"2023-01-24T07:01:18.868110Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# 2. Standard pre-processing","metadata":{}},{"cell_type":"markdown","source":"## 2.1. Mitochondrial/ribosomal contamination","metadata":{}},{"cell_type":"markdown","source":"From [here](http://www.bioinformatics.babraham.ac.uk/training/10XRNASeq/seurat_workflow.html): Single cell datasets can be filled with large numbers of reads coming from mitochondria. These often indicate a sick cell undergoing apoptosis. We want to check for this. We can find the gene names as the rownames of the @assays$RNA@counts slot of the Seurat object and we identify the mitochondrial genes by their names starting with “MT-”. Be aware that in other species the naming of the mitochondrial genes may not be the same so you’d need to adjust the pattern.","metadata":{}},{"cell_type":"code","source":"#check pattern for mitochondrial genes\nMT <- grep(\"MT-\",rownames(seur_RNA@assays$RNA@counts),value = TRUE)\n\ncat(\"There are\", length(MT), \"mitochondrial genes in the dataset:\")\nMT","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-01-24T07:01:18.870205Z","iopub.execute_input":"2023-01-24T07:01:18.871437Z","iopub.status.idle":"2023-01-24T07:01:18.902307Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# calculate mitochondrial percentage and asign this to the \"meta.data\" dataframe\nseur_RNA[[\"percent.mt\"]] <- PercentageFeatureSet(seur_RNA, pattern = \"MT-\")\n\n# check the dataframe with the new column, cells with the most mitochondrial contamination first\ncat(\"Mitochondrial contamination varies between\", round(min(seur_RNA@meta.data$percent.mt),2), \"%\",\n   \"and\", round(max(seur_RNA@meta.data$percent.mt),2), \"%\",\n   \"\\n\\n Number of cells with mitochondrial contamination > 5%:\",\n    nrow(seur_RNA@meta.data[seur_RNA@meta.data$percent.mt > 5, ]), \"of\",nrow(seur_RNA@meta.data),\n   \"\\n\\n Number of cells with mitochondrial contamination > 10%:\",\n    nrow(seur_RNA@meta.data[seur_RNA@meta.data$percent.mt > 10, ]), \"of\",nrow(seur_RNA@meta.data))\n\ncat(\"Head of the data\")\nseur_RNA@meta.data %>%\n  arrange(desc(percent.mt)) %>%\n  head(n = 5)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-01-24T07:01:18.904458Z","iopub.execute_input":"2023-01-24T07:01:18.905660Z","iopub.status.idle":"2023-01-24T07:01:19.528836Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"From [here](http://www.bioinformatics.babraham.ac.uk/training/10XRNASeq/seurat_workflow.html): Ribosomal genes also tend to be very highly represented, and can vary between cell types, so it can be instructive to see how prevalent they are in the data. These are ribosomal protein genes rather than the actual rRNA, so they’re more a measure of the translational activity of the cell rather than the cleanliness of the polyA selection.","metadata":{}},{"cell_type":"code","source":"#check pattern for ribosomal genes\nrb <- grep(\"-RP\",rownames(seur_RNA@assays$RNA@counts),value = TRUE)\n#rb\ncat(\"There are\", length(rb), \"ribosomal genes in the dataset.\")","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-01-24T07:01:19.530860Z","iopub.execute_input":"2023-01-24T07:01:19.532060Z","iopub.status.idle":"2023-01-24T07:01:19.555465Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# calculate percentage of ribosomal RNA and asign this to the \"meta.data\" dataframe\nseur_RNA[[\"percent.rb\"]] <- PercentageFeatureSet(seur_RNA, pattern = \"-RP\")\n\n# check the dataframe with the new column, cells with the most mitochondrial contamination first\ncat(\"Ribosomal contamination varies between\", round(min(seur_RNA@meta.data$percent.rb),2), \"%\",\n   \"and\", round(max(seur_RNA@meta.data$percent.rb),2), \"%\",\n   \"\\n\\n Number of cells with ribosomal contamination > 30%:\",\n    nrow(seur_RNA@meta.data[seur_RNA@meta.data$percent.rb > 30, ]), \"of\",nrow(seur_RNA@meta.data),\n   \"\\n\\n Number of cells with ribosomal contamination > 40%:\",\n    nrow(seur_RNA@meta.data[seur_RNA@meta.data$percent.rb > 40, ]), \"of\",nrow(seur_RNA@meta.data))\n\ncat(\"Head of the data\")\nseur_RNA@meta.data %>%\n  arrange(desc(percent.rb)) %>%\n  head(n = 5)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-01-24T07:01:19.557503Z","iopub.execute_input":"2023-01-24T07:01:19.558746Z","iopub.status.idle":"2023-01-24T07:01:20.247651Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"rm(MT, rb)\ngc()","metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"execution":{"iopub.status.busy":"2023-01-24T07:01:20.249670Z","iopub.execute_input":"2023-01-24T07:01:20.250948Z","iopub.status.idle":"2023-01-24T07:01:20.633946Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## 2.2. Visualization QC","metadata":{}},{"cell_type":"markdown","source":"Next, we visualize QC metrics.","metadata":{}},{"cell_type":"code","source":"# visualize the quality of our dataset with Vlnplot().\nfig(25,9)\nVlnPlot(seur_RNA, features = c(\"nFeature_RNA\", \"nCount_RNA\", \"percent.mt\", \"percent.rb\"), ncol = 4)\n\ncat(\"The dataset contains cells with\", min(seur_RNA@meta.data$nFeature_RNA), \"up to\", max(seur_RNA@meta.data$nFeature_RNA), \n    \"features (genes) and\", min(seur_RNA@meta.data$nCount_RNA), \"up to\", max(seur_RNA@meta.data$nCount_RNA), \"gene copies.\")","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-01-24T07:01:20.635980Z","iopub.execute_input":"2023-01-24T07:01:20.637186Z","iopub.status.idle":"2023-01-24T07:01:26.347221Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# FeatureScatter is typically used to visualize feature-feature relationships, but can be used\n# for anything calculated by the object, i.e. columns in object metadata, PC scores etc.\nfig(30,8)\nplt1 <- FeatureScatter(seur_RNA, feature1 = \"nCount_RNA\", feature2 = \"nFeature_RNA\",pt.size = 0.03) + theme_classic(base_size = 24)\nplt2 <- FeatureScatter(seur_RNA, feature1 = \"nCount_RNA\", feature2 = \"percent.mt\",pt.size = 0.03) + theme_classic(base_size = 24)\nplt3 <- FeatureScatter(seur_RNA, feature1 = \"nCount_RNA\", feature2 = \"percent.rb\",pt.size = 0.03) + theme_classic(base_size = 24)\n\nplt1 + plt2 + plt3\ncat(\"nFeature_RNA is the number of genes detected in each cell.\nnCount_RNA is the total number of molecules detected within a cell.\nAnd each dot in the following plots represents a cell.\nThe number above each plot denotes the correlations between x-axis and y-axis.\")","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-01-24T07:01:26.349471Z","iopub.execute_input":"2023-01-24T07:01:26.350780Z","iopub.status.idle":"2023-01-24T07:01:30.918104Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# 3. Normalizing the data","metadata":{}},{"cell_type":"markdown","source":"As discussed [here](http://www.kaggle.com/datasets/alexandervc/research-project-01-around-multimodal-singlecell/discussion/378460), 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)) outperform other methods in application for common downstream analytical tasks such as variable gene selection, dimensional reduction, and differential expression. We also compared different normalization techniques in [this notebook](http://www.kaggle.com/code/antoninadolgorukova/mmscel-raw-data-normalization-with-r).\n\nFor this we use SCTransform function from [sctransform packadge](http://github.com/saketkc/sctransform) (the last update - v2 regularization, vignette is [here](http://satijalab.org/seurat/articles/sctransform_v2_vignette.html#identify-differential-expressed-genes-across-conditions-1), accompanying article - [here](http://www.biorxiv.org/content/10.1101/2021.07.07.451498v1). \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() # takes ~10 min for the first time installation\nif (!requireNamespace(\"BiocManager\", quietly = TRUE)) install.packages(\"BiocManager\")\nBiocManager::install(\"glmGamPoi\")\ntoc()\ngc()","metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"execution":{"iopub.status.busy":"2023-01-24T07:01:30.921479Z","iopub.execute_input":"2023-01-24T07:01:30.923795Z","iopub.status.idle":"2023-01-24T07:10:11.281390Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"tic() # ~4 min or 3.2 hours if conserve.memory = TRUE on 40% of cells\nseur_RNA <- SCTransform(seur_RNA, vst.flavor = \"v2\",\n                           method=\"glmGamPoi\", verbose = FALSE,\n                           #conserve.memory = TRUE # also note that return.only.var.genes = TRUE is the default\n                       ) \ntoc()\ngc()","metadata":{"_kg_hide-output":true,"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-01-24T07:10:11.283432Z","iopub.execute_input":"2023-01-24T07:10:11.284639Z","iopub.status.idle":"2023-01-24T07:14:27.206146Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Let's plot variance of Pearson residuals vs mean expresseion of genes to see how this normalization method handle technical variability and stabilize the variance across genes with different expression levels.","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(pr.var = seur_RNA@assays$SCT@SCTModel.list$model1@feature.attributes$residual_variance,\n                  mean = seur_RNA@assays$SCT@SCTModel.list$model1@feature.attributes$gmean,\n                  gene = seur_RNA@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":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-01-24T07:14:27.208516Z","iopub.execute_input":"2023-01-24T07:14:27.209837Z","iopub.status.idle":"2023-01-24T07:14:28.146841Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"rm(plt, top_var_genes)\ngc()\nfor (thing in ls()) { message(thing); print(object.size(get(thing)), units='auto') }","metadata":{"_kg_hide-output":true,"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-01-24T07:14:28.149024Z","iopub.execute_input":"2023-01-24T07:14:28.150317Z","iopub.status.idle":"2023-01-24T07:14:28.956307Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# 4. Identification of highly variable features (feature selection)","metadata":{}},{"cell_type":"markdown","source":"We next calculate a subset of features that exhibit high cell-to-cell variation in the dataset (i.e, they are highly expressed in some cells, and lowly expressed in others). We and others have found that focusing on these genes in downstream analysis helps to highlight biological signal in single-cell datasets.\n\nIn Seurat the procedure is implemented in the FindVariableFeatures function described in detail [here](http://satijalab.org/seurat/). However, since we have used SCTransform function from sctransform packadge, we already have fariable features (3000 by default) stored in","metadata":{}},{"cell_type":"code","source":"# Identify the 10 most variable genes with Seurat\ntop10_genes <- head(VariableFeatures(seur_RNA), 10)\ncat(\"The top 10 most variable genes among cells are:\\n\\n\", top10_genes)\n\n# plot variable features with and without labels of the top 10 genes\nfig(25,10)\nplt1 <- VariableFeaturePlot(seur_RNA) + theme_classic(base_size = 24)\nplt2 <- LabelPoints(plot = plt1, points = top10_genes, repel = TRUE, xnudge = 0, ynudge = 0) + theme_classic(base_size = 24)\nplt1 + plt2","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-01-24T07:14:28.958517Z","iopub.execute_input":"2023-01-24T07:14:28.959789Z","iopub.status.idle":"2023-01-24T07:14:31.343885Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# 5. Dimensionality Reduction: PCA","metadata":{}},{"cell_type":"markdown","source":"Single-cell gene expression holds great potential to identify which features are important, but the problem is that these are high dimensional datasets. Essentially, every gene represents one dimension, one feature. This makes them quite challenging to understand and to plot. In most cases this is a far too complicated representation of reality anyway. The goal of performing dimensionality reduction is to find the approximate underlying dimensionality of the dataset.","metadata":{}},{"cell_type":"markdown","source":"Next we perform PCA. By default, only the previously determined variable features are used as input, but can be defined using features argument if you wish to choose a different subset. Note that the features must be present in the scaled data. Any requested features that are not scaled or have 0 variance will be dropped, and the PCA will be run using the remaining features. Since we have used SCTransform function for data normalisation, we do not need additional scaling.","metadata":{}},{"cell_type":"code","source":"seur_RNA <- RunPCA(seur_RNA, npcs = 50, verbose = FALSE)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-01-24T07:14:31.346638Z","iopub.execute_input":"2023-01-24T07:14:31.348224Z","iopub.status.idle":"2023-01-24T07:14:47.908764Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Examine and visualize PCA results a few different ways\ncat(\"5 most important positive and negative expression features in the first 5 principle components\\n\")\nprint(seur_RNA[[\"pca\"]], dims = 1:5, nfeatures = 5)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-01-24T07:14:47.913442Z","iopub.execute_input":"2023-01-24T07:14:47.917237Z","iopub.status.idle":"2023-01-24T07:14:47.955872Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# 6. Plot the first two PCs","metadata":{}},{"cell_type":"code","source":"fig(25,8)\n#Standart plots with seurat DimPlot function\nplt1 <- DimPlot(seur_RNA, reduction = \"pca\", group.by = \"cell_type\") + theme_classic(base_size = 24)\nplt2 <- DimPlot(seur_RNA, reduction = \"pca\", group.by = \"day\") + theme_classic(base_size = 24)\nplt3 <- DimPlot(seur_RNA, reduction = \"pca\", group.by = \"donor\") + theme_classic(base_size = 24)\nplt1 + plt2 + plt3","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-01-24T07:14:47.958662Z","iopub.execute_input":"2023-01-24T07:14:47.960269Z","iopub.status.idle":"2023-01-24T07:14:51.713037Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"pca_color <- seur_RNA@reductions$pca@cell.embeddings %>%\n  as.data.frame %>% \n  select(PC_1, PC_2) %>%\n  mutate(percent_mt = seur_RNA$percent.mt,\n         percent_rb = seur_RNA$percent.rb,\n        nCount_RNA = seur_RNA$nCount_RNA,\n        nFeature_RNA = seur_RNA$nFeature_RNA,\n        day = as.factor(seur_RNA$day),\n        donor = as.factor(seur_RNA$donor),\n        cell_type = seur_RNA$cell_type)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-01-24T07:14:51.716316Z","iopub.execute_input":"2023-01-24T07:14:51.718275Z","iopub.status.idle":"2023-01-24T07:14:51.849167Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig(25,15)\n#Ploting with ggplot\nplt1 <- ggplot(pca_color, aes(x = PC_1, y = PC_2, col = percent_mt)) +\n      geom_point(size=0.8) +\n      theme_classic(base_size = 24) +\n      scale_color_gradient(low='darkblue', high='yellow') +\n      theme(aspect.ratio = 1)\n\nplt2 <- ggplot(pca_color, aes(x = PC_1, y = PC_2, col = percent_rb)) +\n      geom_point(size=0.8) +\n      theme_classic(base_size = 24) +\n      scale_color_gradient(low='darkblue', high='yellow') +\n      theme(aspect.ratio = 1)\n\nplt3 <- ggplot(pca_color, aes(x = PC_1, y = PC_2, col = nCount_RNA)) +\n      geom_point(size=0.8) +\n      theme_classic(base_size = 24) +\n      scale_color_gradient(low='darkblue', high='yellow') +\n      theme(aspect.ratio = 1)\n\nplt4 <- ggplot(pca_color, aes(x = PC_1, y = PC_2, col = nFeature_RNA)) +\n      geom_point(size=0.8) +\n      theme_classic(base_size = 24) +\n      scale_color_gradient(low='darkblue', high='yellow') +\n      theme(aspect.ratio = 1)\n\n\nplt1 + plt2 + plt3 + plt4","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-01-24T07:14:51.851754Z","iopub.execute_input":"2023-01-24T07:14:51.853271Z","iopub.status.idle":"2023-01-24T07:14:57.929690Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# 7. PCA_2 vs top 100 RNA correlated with the protein of interest","metadata":{}},{"cell_type":"markdown","source":"- Load normalized protein and RNA data and sample the same 40% of cells","metadata":{}},{"cell_type":"code","source":"PC_2 <- Embeddings(object = seur_RNA, reduction = \"pca\")[, \"PC_2\"]\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)\nmat_prot <- mat_prot[selected_cells, ]\n\n# Load RNA normalized data\npath <- \"/kaggle/input/sparse-measurement-data-open-problems-multimodal/sp_train_cite_inputs.rds\"\nmat_RNA <- readRDS(path)\nmat_RNA <- mat_RNA[selected_cells, ]\n\ngc()","metadata":{"_kg_hide-output":true,"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-01-24T07:14:57.932121Z","iopub.execute_input":"2023-01-24T07:14:57.933468Z","iopub.status.idle":"2023-01-24T07:16:10.258198Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"- Select the protein of interest","metadata":{}},{"cell_type":"code","source":"prot_of_interest <- \"CD36\"","metadata":{"execution":{"iopub.status.busy":"2023-01-24T08:16:24.769029Z","iopub.execute_input":"2023-01-24T08:16:24.770765Z","iopub.status.idle":"2023-01-24T08:16:24.783502Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"- Load protein-RNA Spearman correlations data (generated [in this notebook](http://www.kaggle.com/code/antoninadolgorukova/mmscel-protein-analysis)) and select correlation coefficients for the protein of interest","metadata":{}},{"cell_type":"code","source":"# Load correlations data (all proteins vs all RNA, Spearman correlation)\npath <- \"/kaggle/input/mmscel-protein-analysis/Prot-RNA_corr_all.csv\"\nprot_RNA_corr <- read.csv(path)\n\nprot_RNA_corr <- select(prot_RNA_corr, c(\"RNA\", all_of(prot_of_interest)) ) %>% na.omit\nnames(prot_RNA_corr)[names(prot_RNA_corr) == prot_of_interest] <- \"corr_coeff\"\nprot_RNA_corr <- prot_RNA_corr[order(prot_RNA_corr$corr_coeff, decreasing = TRUE), ]\n\ncat(\"Top 10 RNA positively correlated with\", prot_of_interest)\nhead(prot_RNA_corr, n = 10)\n\ncat(\"Top 10 RNA negatively correlated with\", prot_of_interest)\ntail(prot_RNA_corr, n = 10)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-01-24T08:16:27.293447Z","iopub.execute_input":"2023-01-24T08:16:27.295166Z","iopub.status.idle":"2023-01-24T08:16:31.587057Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"- Combine all the data together (top 100 RNA, the protein of interest) with PC_1 and PC_2 data","metadata":{}},{"cell_type":"code","source":"prot_RNA_corr <- prot_RNA_corr[order(abs(prot_RNA_corr$corr_coeff), decreasing = TRUE), ]\n\n# merge PC_2 with top 100 RNA (norm data)\nmy_RNA <- mat_RNA[, prot_RNA_corr[1:100, \"RNA\"]]\nmy_RNA <- as.matrix(my_RNA)\n\npca_color_long <- merge(pca_color, my_RNA, by = 0)\n\n#merge with CD protein\nmy_prot <- mat_prot[, prot_of_interest] %>% as.data.frame\ncolnames(my_prot) <- prot_of_interest\nmy_prot$Row.names <- rownames(my_prot)\n\npca_color_long <- merge(pca_color_long, my_prot, by = \"Row.names\", all.x = TRUE)\n\n# make long to plot with for loop all RNA\npca_color_long <- pivot_longer(pca_color_long, cols = (ncol(pca_color)+2):ncol(pca_color_long),\n                          names_to = \"RNA\", values_to = \"Expression\")\n\n# add correlation coefficients\npca_color_long <- merge(pca_color_long, prot_RNA_corr, by = \"RNA\", all.x = TRUE)\n\n# add column with RNA name and correlation coefficient in braskets\npca_color_long <- pca_color_long %>% \n   mutate(RNA_cor = paste0(RNA, \" (\", round(corr_coeff, 2), \")\") ) %>%\n   rename(mol = RNA)\npca_color_long$RNA_cor <- gsub(\"\\\\(NA)\", \"\", pca_color_long$RNA_cor)\n#head(pca_color_long)\n#pca_color_long %>% group_by(RNA_cor) %>% summarise(n = n())","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-01-24T08:16:38.459263Z","iopub.execute_input":"2023-01-24T08:16:38.460922Z","iopub.status.idle":"2023-01-24T08:17:20.896716Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"rm(my_RNA,PC_2, prot_RNA_corr)\ngc()","metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"execution":{"iopub.status.busy":"2023-01-24T08:17:36.139311Z","iopub.execute_input":"2023-01-24T08:17:36.140931Z","iopub.status.idle":"2023-01-24T08:17:38.772435Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"- Visualise the CD vs PC_2/PC_1 relationships in all and EryP cells","metadata":{}},{"cell_type":"code","source":"#all cells\nall_cells <- filter(pca_color_long, grepl(prot_of_interest, mol))\ncat(\"The head of all data\")\nhead(all_cells)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-01-24T08:17:40.920634Z","iopub.execute_input":"2023-01-24T08:17:40.922182Z","iopub.status.idle":"2023-01-24T08:17:42.193555Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig(24,10)\nplt <- ggscatter(all_cells[all_cells$mol == prot_of_interest,], x = \"Expression\", y = \"PC_2\", color = \"PC_1\",\n          shape = 20, size = 3, #, alpha = 0.3\n          add = \"reg.line\", conf.int = TRUE,\n          cor.coef = TRUE, cor.coef.size = 6,\n          cor.coeff.args = list(method = \"spearman\", label.sep = \"\\n\"),\n          add.params = list(fill = \"lightgray\"),\n          xlab = prot_of_interest,\n          title = paste0(\"Spearman correlation, all cells (n = \", nrow(all_cells[all_cells$mol == prot_of_interest,]), \")\"),\n          ggtheme = theme_bw(base_size = 20))\n\nplt1 <- ggMarginal(plt, type=\"histogram\")\n\nplt <- ggscatter(all_cells[grepl(\"ENSG\", all_cells$mol),], x = \"Expression\", y = \"PC_2\", color = \"PC_1\",\n          shape = 20, size = 3, #, alpha = 0.3\n          add = \"reg.line\", conf.int = TRUE,\n          cor.coef = TRUE, cor.coef.size = 6,\n          cor.coeff.args = list(method = \"spearman\", label.sep = \"\\n\"),\n          add.params = list(fill = \"lightgray\"),\n          xlab = paste(prot_of_interest, \"RNA\"),\n          title = paste0(\"Spearman correlation, all cells (n = \", nrow(all_cells[all_cells$mol == prot_of_interest,]), \")\"),\n          ggtheme = theme_bw(base_size = 20))\n\nplt2 <- ggMarginal(plt, type=\"histogram\")\n\ngrid.arrange(plt1, plt2, ncol=2)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-01-24T08:17:47.573753Z","iopub.execute_input":"2023-01-24T08:17:47.575252Z","iopub.status.idle":"2023-01-24T08:17:53.527521Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig(15,20) \n\nplt <- plot_ly(sample_n(all_cells[all_cells$mol == prot_of_interest,], 5000), x = ~PC_1, y = ~PC_2, z = ~Expression, color = ~cell_type,\n       size = 0.5) %>% \n   add_markers() %>%\n   layout(title = paste(prot_of_interest, \"expression vs first two principal components (all cells, sample = 5000)\"))\n\nhtmlwidgets::saveWidget(partial_bundle(plt), file = paste0(prot_of_interest, \"_scatter3d.HTML\"), selfcontained = TRUE)\nutils::browseURL(paste0(prot_of_interest, \"_scatter3d.HTML\"))","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-01-24T08:18:07.344548Z","iopub.execute_input":"2023-01-24T08:18:07.346067Z","iopub.status.idle":"2023-01-24T08:18:08.433974Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig(15,20) \n\nplt <- plot_ly(sample_n(all_cells[grepl(\"ENSG\", all_cells$mol),], 5000), x = ~PC_1, y = ~PC_2, z = ~Expression, color = ~cell_type,\n       size = 0.5) %>% \n   add_markers() %>%\n   layout(title = paste(prot_of_interest, \"RNA expression vs first two principal components (all cells, sample = 5000)\"))\n\nhtmlwidgets::saveWidget(partial_bundle(plt), file = paste0(prot_of_interest, \"RNA_scatter3d.HTML\"), selfcontained = TRUE)\nutils::browseURL(paste0(prot_of_interest, \"RNA_scatter3d.HTML\"))","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-01-24T08:18:19.880141Z","iopub.execute_input":"2023-01-24T08:18:19.881665Z","iopub.status.idle":"2023-01-24T08:18:20.868194Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# EryP cells\neryp <- filter(pca_color_long, grepl(prot_of_interest, mol) & cell_type == \"EryP\")\neryp <- eryp %>% select(\"Row.names\", \"PC_1\", \"PC_2\", \"day\", \"donor\",\"cell_type\", \"mol\",\"Expression\") %>%\n  pivot_wider(names_from = \"mol\", values_from = \"Expression\")\ncat(\"The head the subset with EryP cells, the protein of interest and its RNA\")\nhead(eryp)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-01-24T08:07:27.409664Z","iopub.execute_input":"2023-01-24T08:07:27.411541Z","iopub.status.idle":"2023-01-24T08:07:28.944466Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig(24,8)\nplt <- ggscatter(eryp, x = prot_of_interest, y = \"PC_2\", color = \"PC_1\",\n          shape = 20, size = 3, #alpha = 0.3, \n          add = \"reg.line\", conf.int = TRUE,\n          cor.coef = TRUE, cor.coef.size = 6,\n          cor.coeff.args = list(method = \"spearman\", label.sep = \"\\n\"),\n          title = paste0(\"Spearman correlation,\\nEryP cells (n = \", nrow(eryp), \")\"),\n          ggtheme = theme_bw(base_size = 20)) +\n          theme(aspect.ratio=1)\n\nplt1 <- ggMarginal(plt, type=\"histogram\")\n\nplt <- ggscatter(eryp, x = prot_of_interest, y = \"PC_1\", color = \"PC_2\",\n          shape = 20, size = 3, #alpha = 0.3,\n          add = \"reg.line\", conf.int = TRUE,\n          cor.coef = TRUE, cor.coef.size = 6,\n          cor.coeff.args = list(method = \"spearman\", label.sep = \"\\n\"),\n          title = paste0(\"Spearman correlation,\\nEryP cells (n = \", nrow(eryp), \")\"),\n          ggtheme = theme_bw(base_size = 20)) +\n          theme(aspect.ratio=1)\n\nplt2 <- ggMarginal(plt, type=\"histogram\")\n\ngrid.arrange(plt1, plt2, ncol=2) %>% suppressMessages","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-01-24T08:07:32.215418Z","iopub.execute_input":"2023-01-24T08:07:32.217711Z","iopub.status.idle":"2023-01-24T08:07:34.685475Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"- Fit linear regression model for CD vs PC_2, PC_1 in EryP cells","metadata":{}},{"cell_type":"code","source":"reg_eryp = lm(get(prot_of_interest) ~ PC_2 + PC_1, data = eryp) \nsummary(reg_eryp)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-01-24T08:07:42.530945Z","iopub.execute_input":"2023-01-24T08:07:42.532651Z","iopub.status.idle":"2023-01-24T08:07:42.556147Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"- Predict the protein of interest by PC_2 and PC_1 in EryP cells","metadata":{}},{"cell_type":"code","source":"eryp$Predicted_CD_expression <- predict(reg_eryp, select(eryp, PC_2, PC_1))\nhead(eryp)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-01-24T08:07:46.448186Z","iopub.execute_input":"2023-01-24T08:07:46.449736Z","iopub.status.idle":"2023-01-24T08:07:46.510767Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"cat(\"The Spearman correlation between\", prot_of_interest, \"expression and its predicted expression (EryP cells) = \",\n    cor(eryp$Predicted_CD_expression, eryp[, prot_of_interest]))\ncat(\"\\nThe Spearman correlation between\", prot_of_interest, \"expression and its RNA (EryP cells)= \",\n    cor(eryp[, names(eryp)[grepl(\"ENSG\", names(eryp))]], eryp[, prot_of_interest]))","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-01-24T08:07:49.818041Z","iopub.execute_input":"2023-01-24T08:07:49.819599Z","iopub.status.idle":"2023-01-24T08:07:49.837415Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"- Calculate median expression of each RNA for binarisation","metadata":{}},{"cell_type":"code","source":"# calculate median expression of each RNA for binarisation\ntic()\npca_color_long$`Expression\\nto median` <- NA\nfor(mol in unique(pca_color_long$mol) ) {\n    \n    tmp <- pca_color_long[pca_color_long$mol == mol, ]\n    med <- median(tmp$Expression)\n    pca_color_long[pca_color_long$mol == mol, ]$`Expression\\nto median` <-\n      ifelse(pca_color_long[pca_color_long$mol == mol, ]$Expression > med, \"higher\", \"lower\")\n}\n\n# exclude \"outliers\" - extrim values\n\npca_color_long$`Expression\\n(corrected)` <- NA\nfor(mol in unique(pca_color_long$mol) ) {\n    \n    tmp <- pca_color_long[pca_color_long$mol == mol, ]\n    q5 <- quantile(tmp$Expression, 0.05)\n    q95 <- quantile(tmp$Expression, 0.95)\n    pca_color_long[pca_color_long$mol == mol, ]$`Expression\\n(corrected)` <-\n      ifelse(pca_color_long[pca_color_long$mol == mol, ]$Expression > q95, q95,\n      ifelse(pca_color_long[pca_color_long$mol == mol, ]$Expression < q5, q5,\n            pca_color_long[pca_color_long$mol == mol, ]$Expression) )\n}\ntoc()\n\nrm(med, mol, q5, q95,tmp)\ngc()","metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"execution":{"iopub.status.busy":"2023-01-24T08:08:07.559348Z","iopub.execute_input":"2023-01-24T08:08:07.560915Z","iopub.status.idle":"2023-01-24T08:09:43.416125Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"- Plot top 100 RNA correlated with the protein of interest","metadata":{}},{"cell_type":"code","source":"fig(22,6)\n\ncat(\"RNA expression was corrected using 5 and 95 percentiles as a cut-off\\n\\n\")\npca_color_long <- pca_color_long[order(pca_color_long$corr_coeff, decreasing = TRUE), ]\n\nfor (p in unique(pca_color_long$mol) ) {\n    \n    plt1 <- ggplot(pca_color_long[pca_color_long$mol == p, ],\n                   aes(x = PC_1, y = PC_2, col = `Expression\\n(corrected)`) ) +\n      geom_point(size=0.8) +\n      theme_classic(base_size = 18) +\n      scale_color_gradient(low='darkblue', high='yellow') +\n      theme(aspect.ratio = 1) +\n      facet_wrap( ~ RNA_cor)\n    \n    plt2 <- ggplot(pca_color_long[pca_color_long$mol == p, ],\n                   aes(x = PC_1, y = PC_2, col = `Expression\\nto median`) ) +\n      geom_point(size=0.8) +\n      theme_classic(base_size = 18) +\n      scale_colour_manual(values = c('yellow', 'darkblue')) +\n      theme(aspect.ratio = 1) +\n      facet_wrap( ~ RNA_cor)\n    \n   print(plt1 + plt2 + plot_layout(ncol = 2))\n   print(paste0(\"-----------------------------\", p, \"----------------------------------\"))\n}","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-01-24T08:12:51.342659Z","iopub.execute_input":"2023-01-24T08:12:51.345183Z","iopub.status.idle":"2023-01-24T08:13:05.635240Z"},"trusted":true},"execution_count":null,"outputs":[]}]}