{"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":"This work is based on the package [tutorial](https://satijalab.org/seurat/archive/v3.2/pbmc3k_tutorial.html) and the [REMINHO's notebook](http://https://www.kaggle.com/code/reminho/basic-scrnaseq-tutorial) with another example of data processing with Seurat.\n\nData: 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)) \nsuppressMessages(library(ggrepel))\nsuppressMessages(library(ggplot2))\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,"execution":{"iopub.status.busy":"2023-01-18T07:24:02.471116Z","iopub.execute_input":"2023-01-18T07:24:02.509981Z","iopub.status.idle":"2023-01-18T07:24:02.532721Z"},"_kg_hide-input":true,"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":{"execution":{"iopub.status.busy":"2023-01-18T07:24:13.068171Z","iopub.execute_input":"2023-01-18T07:24:13.070574Z","iopub.status.idle":"2023-01-18T07:24:13.091389Z"},"_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, ] #uncomment touse 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\")\nmat_RNA[c(\"ENSG00000135218_CD36\", \"ENSG00000029534_ANK1\", \"ENSG00000160789_LMNA\"), 1:30]","metadata":{"execution":{"iopub.status.busy":"2023-01-18T07:24:17.086766Z","iopub.execute_input":"2023-01-18T07:24:17.088129Z","iopub.status.idle":"2023-01-18T07:25:17.740567Z"},"_kg_hide-input":true,"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,"execution":{"iopub.status.busy":"2023-01-18T07:25:21.851976Z","iopub.execute_input":"2023-01-18T07:25:21.853467Z","iopub.status.idle":"2023-01-18T07:25:22.644540Z"},"_kg_hide-input":true,"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":{"execution":{"iopub.status.busy":"2023-01-18T07:25:30.675908Z","iopub.execute_input":"2023-01-18T07:25:30.677250Z","iopub.status.idle":"2023-01-18T07:25:51.581460Z"},"_kg_hide-input":true,"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,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# 2. Standard pre-processing","metadata":{}},{"cell_type":"markdown","source":"The steps below encompass the standard pre-processing workflow for scRNA-seq data in Seurat. These represent the selection and filtration of cells based on QC metrics, data normalization and scaling, and the detection of highly variable features.","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-18T08:49:14.935476Z","iopub.execute_input":"2023-01-18T08:49:14.936918Z","iopub.status.idle":"2023-01-18T08:49:14.968161Z"},"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":{"execution":{"iopub.status.busy":"2023-01-18T07:26:10.523677Z","iopub.execute_input":"2023-01-18T07:26:10.525212Z","iopub.status.idle":"2023-01-18T07:26:11.022787Z"},"_kg_hide-input":true,"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":{"execution":{"iopub.status.busy":"2023-01-18T07:26:19.503744Z","iopub.execute_input":"2023-01-18T07:26:19.505272Z","iopub.status.idle":"2023-01-18T07:26:19.530298Z"},"_kg_hide-input":true,"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":{"execution":{"iopub.status.busy":"2023-01-18T07:26:25.125119Z","iopub.execute_input":"2023-01-18T07:26:25.126547Z","iopub.status.idle":"2023-01-18T07:26:25.629104Z"},"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## 2.2. Percentage of Largest Gene","metadata":{}},{"cell_type":"markdown","source":"From [here](http://www.bioinformatics.babraham.ac.uk/training/10XRNASeq/seurat_workflow.html#Setting_QC_Cutoffs): We run apply over the columns (cells) and calculate what percentage of the data comes from the single most observed gene. A high proportion of data dominated by a single gene is a metric which could either give biological context or indicate a technical problem, depending on what the gene is.\n\nThe calculation requires sparse to dense coercion which takes a lot of memory (impossible to do in the full dataset with ~70K cells).","metadata":{}},{"cell_type":"code","source":"largest <- seur_RNA\n\nlargest$largest_count <- apply(\n  largest@assays$RNA@counts, 2, max)\n\nlargest$largest_index <- apply(\n    largest@assays$RNA@counts, 2, which.max)\n\nlargest$largest_gene <- rownames(largest)[largest$largest_index]\nlargest$percent.Largest.Gene <- 100 * largest$largest_count / largest$nCount_RNA\n\ncat(\"Head of the data with largest gene calculated\")\nhead(largest)\nseur_RNA$largest_gene <- largest$largest_gene\nseur_RNA$percent.Largest.Gene <- largest$percent.Largest.Gene","metadata":{"execution":{"iopub.status.busy":"2023-01-18T07:26:30.247588Z","iopub.execute_input":"2023-01-18T07:26:30.249080Z","iopub.status.idle":"2023-01-18T07:27:06.022928Z"},"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"rm(largest, MT, rb)\ngc()","metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Let's see what the names of the largest genes are (largest_gene column) and in how many cells they were detected (n columnn). ","metadata":{}},{"cell_type":"code","source":"qc.metrics <- as_tibble(\n    seur_RNA[[]],\n    rownames = \"cell.id\"\n)\n\n#head(qc.metrics)\n\nqc.metrics %>%\n  group_by(largest_gene) %>%\n  count() %>%\n  arrange(desc(n))","metadata":{"execution":{"iopub.status.busy":"2023-01-18T07:27:15.870117Z","iopub.execute_input":"2023-01-18T07:27:15.871697Z","iopub.status.idle":"2023-01-18T07:27:15.930363Z"},"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## 2.3. Visualization","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":{"execution":{"iopub.status.busy":"2023-01-18T07:27:28.525966Z","iopub.execute_input":"2023-01-18T07:27:28.527515Z","iopub.status.idle":"2023-01-18T07:27:32.678672Z"},"_kg_hide-input":true,"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(18,12)\nplot1 <- FeatureScatter(seur_RNA, feature1 = \"nCount_RNA\", feature2 = \"nFeature_RNA\",pt.size = 0.03) + theme_classic(base_size = 24)\nplot2 <- FeatureScatter(seur_RNA, feature1 = \"nCount_RNA\", feature2 = \"percent.mt\",pt.size = 0.03) + theme_classic(base_size = 24)\nplot3 <- FeatureScatter(seur_RNA, feature1 = \"nCount_RNA\", feature2 = \"percent.rb\",pt.size = 0.03) + theme_classic(base_size = 24)\n\n# plot preselected cells^\n#non_na_cells <- rownames(seur_RNA@meta.data[!is.na(seur_RNA$percent.Largest.Gene), ]) + theme_classic(base_size = 24) \n# add cells = non_na_cells argument to FeatureScatter call\n\nplot4 <- FeatureScatter(seur_RNA,feature1 = \"nCount_RNA\", feature2 = \"percent.Largest.Gene\",pt.size = 0.03, ) + theme_classic(base_size = 24)\n\n# Same by color from metadata\n#FeatureScatter(seur_RNA, feature1 = \"nCount_RNA\", feature2 = \"nFeature_RNA\",\n#                        group.by = \"cell_type\",pt.size = 0.03) + theme_classic(base_size = 24)\n\nplot1 + plot2 + plot3 + plot4\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":{"execution":{"iopub.status.busy":"2023-01-18T07:27:37.031905Z","iopub.execute_input":"2023-01-18T07:27:37.034111Z","iopub.status.idle":"2023-01-18T07:27:41.191152Z"},"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"It seems also reasonable to check if the largest genes are correlated with mitochondrial or ribosomal.","metadata":{}},{"cell_type":"code","source":"fig(20, 8)\n\nplot1 <- FeatureScatter(seur_RNA,feature1 = \"percent.mt\", feature2 = \"percent.Largest.Gene\", pt.size = 0.03) + theme_classic(base_size = 24)\nplot2 <- FeatureScatter(seur_RNA,feature1 = \"percent.rb\", feature2 = \"percent.Largest.Gene\", pt.size = 0.03) + theme_classic(base_size = 24)\n\nplot1 + plot2","metadata":{"execution":{"iopub.status.busy":"2023-01-18T07:27:51.038511Z","iopub.execute_input":"2023-01-18T07:27:51.039799Z","iopub.status.idle":"2023-01-18T07:27:53.149333Z"},"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"From [this](http://holab-hku.github.io/Fundamental-scRNA/downstream.html) tutorial: Depend on the data we analyze, we can use different cutoff for nFeature_RNA and percent.mt. For example, we can filter away cells that have unique feature counts(genes) over 5,000 or less than 300, or cells that have > 10% mitochondrial counts, and see how QC metrics looks like. The filtration can also be made with SCTransform (see 3. Normalizing the data).","metadata":{}},{"cell_type":"code","source":"# Filtration\n#seur_RNA <- subset(seur_RNA, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mt < 15)","metadata":{"_kg_hide-input":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# 3. Normalizing the data","metadata":{}},{"cell_type":"markdown","source":"By default, Seurat package has a global-scaling normalization method \"LogNormalize\" in NormalizeData function 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. In 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 ))\n\nHowever, 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 time for the first time installation\nif (!requireNamespace(\"BiocManager\", quietly = TRUE)) install.packages(\"BiocManager\")\nBiocManager::install(\"glmGamPoi\")\ntoc()\ngc()","metadata":{"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"tic() # ~3.2 hours with 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,"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":{"execution":{"iopub.status.busy":"2023-01-18T07:37:15.129955Z","iopub.execute_input":"2023-01-18T07:37:15.131395Z","iopub.status.idle":"2023-01-18T07:37:15.884107Z"},"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"rm(plt, qc.metrics, plot1, plot2, plot3, plot4, 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,"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":"# Normalizing the data  with Seurat\n#seur_RNA <- NormalizeData(seur_RNA, normalization.method = \"LogNormalize\", scale.factor = 10000)\n# 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\", top10_genes)\n\n# plot variable features with and without labels of the top 10 genes\nfig(22,8)\nplot1 <- VariableFeaturePlot(seur_RNA) + theme_classic(base_size = 24)\nplot2 <- LabelPoints(plot = plot1, points = top10_genes, repel = TRUE, xnudge = 0, ynudge = 0) + theme_classic(base_size = 24)\nplot1 + plot2","metadata":{"execution":{"iopub.status.busy":"2023-01-18T07:39:30.221227Z","iopub.execute_input":"2023-01-18T07:39:30.222893Z","iopub.status.idle":"2023-01-18T07:39:31.996336Z"},"_kg_hide-input":true,"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. There are several (non)-linear methods available in Seurat. I will first perform a linear PCA and later demonstrate non-linear dimensionality reduction using t-distributed stochastic neighbor embedding (tSNE) and Uniform Manifold Approximation and Projection (UMAP).","metadata":{}},{"cell_type":"markdown","source":"## 5.1. Scaling the data","metadata":{}},{"cell_type":"markdown","source":"From [this](http://nbisweden.github.io/workshop-scRNAseq/labs/compiled/seurat/seurat_02_dim_reduction.html) tutorial: Since each gene has a different expression level, it means that genes with higher expression values will naturally have higher variation that will be captured by PCA. This means that we need to somehow give each gene a similar weight when performing PCA. The common practice is to center and scale each gene before performing PCA. This exact scaling is called Z-score normalization it is very useful for PCA, clustering and plotting heatmaps.\n\nFrom the [Seurat tutorial](https://satijalab.org/seurat/archive/v3.2/pbmc3k_tutorial.html): The linear transformation ('scaling') is a standard pre-processing step prior to dimensional reduction techniques like PCA and is an essential step in the Seurat workflow, but only on genes that will be used as input to PCA. \n\nIn Seurat the procedure is implemented in the ScaleData function described in detail [here](http://satijalab.org/seurat/). However, since we have used SCTransform function for data normalisation, we do not need additional scaling.","metadata":{}},{"cell_type":"code","source":"#tic()\n# Scaling all genes\n#all.genes <- rownames(var_features)\n#var_features <- ScaleData(var_features, features = all.genes)\n\n# Perform scaling on the previously identified variable features\n#var_features <- ScaleData(var_features)\n#toc()","metadata":{"execution":{"iopub.status.busy":"2023-01-14T23:07:07.175767Z","iopub.execute_input":"2023-01-14T23:07:07.178589Z","iopub.status.idle":"2023-01-14T23:07:22.543485Z"},"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## 5.2. PCA: Finding the most important features in principles components","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.","metadata":{}},{"cell_type":"code","source":"seur_RNA <- RunPCA(seur_RNA, npcs = 50, verbose = FALSE)","metadata":{"execution":{"iopub.status.busy":"2023-01-18T07:40:51.464577Z","iopub.execute_input":"2023-01-18T07:40:51.466264Z","iopub.status.idle":"2023-01-18T07:41:01.386567Z"},"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Seurat provides several useful ways of visualizing both cells and features that define the PCA, including VizDimReduction, DimPlot, and DimHeatmap.","metadata":{}},{"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\")\nprint(seur_RNA[[\"pca\"]], dims = 1:5, nfeatures = 5)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-01-18T08:56:13.285891Z","iopub.execute_input":"2023-01-18T08:56:13.287822Z","iopub.status.idle":"2023-01-18T08:56:13.357575Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig(22,15)\n# plot PCA, first 5 PCs\ncat(\"Top genes associated with reduction components\")\nVizDimLoadings(seur_RNA, dims = 1:6, reduction = \"pca\", ncol = 3)","metadata":{"execution":{"iopub.status.busy":"2023-01-18T08:57:38.735317Z","iopub.execute_input":"2023-01-18T08:57:38.738861Z","iopub.status.idle":"2023-01-18T08:57:40.023663Z"},"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## 5.3. Compute % variance explained by each PC","metadata":{}},{"cell_type":"markdown","source":"Seurat as of today does not provide an easy function for computing the percent variance explained by each principle component. Thanks to [this answer](https://github.com/satijalab/seurat/issues/982) in the Seurat Github repo, you can calculate it as follows if you wish to do so:","metadata":{}},{"cell_type":"code","source":"# get the percent variance explained by the individual principles components\nmat <- GetAssayData(seur_RNA, assay = \"SCT\", slot = \"scale.data\")\npca <- seur_RNA[[\"pca\"]]\n\n# Get the total variance:\ntotal_variance <- sum(matrixStats::rowVars(mat))\n\neigValues <- (pca@stdev)^2  ## EigenValues\nvarExplained <- eigValues / total_variance\n\n# calculate the % variance explained by the first 10 PCs\ncat(\"The first 50 principle components explain\", \n    scales::label_percent()(round(sum(varExplained[c(1:50)]), 4)), \"of the total variance in the dataset.\")","metadata":{"execution":{"iopub.status.busy":"2023-01-18T07:41:17.884568Z","iopub.execute_input":"2023-01-18T07:41:17.886223Z","iopub.status.idle":"2023-01-18T07:41:18.642913Z"},"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"###  \"Cumulative variance explained\" plot","metadata":{}},{"cell_type":"markdown","source":"If you want to visualize the cumulative variance explained by the first 50 PCs, I provide the following workaround, because there is no easy function in Seurat. If you want to visualize more PCs for some reason, you'd have to change the npcs = 50 argument in the RunPCA() call.","metadata":{}},{"cell_type":"code","source":"# custom plot size\nfig(14, 10)\n\n# store varEplained in a dataframe alongside the PC number\nvar_explained_df <- tibble(pc_number = 1:length(varExplained), var_explained = cumsum(varExplained))\n\n# plot this dataframe\nggplot(var_explained_df, aes(x = pc_number, y = var_explained)) +\ngeom_point(size = 2) +\ngeom_line(size = 0.8) +\nscale_y_continuous(labels = scales::percent_format(accuracy = 1)) +\nlabs(title = \"Cumulative Variance Explained by the first 50 PCs\",\n    x = \"PC\",\n    y = \"Cumulative Variance explained\") +\ntheme_classic(base_size = 24) +\ntheme(plot.title = element_text(size = 22),\n          plot.subtitle = element_text(size = 16),\n          axis.text.x= element_text(size = 15),\n          axis.text.y= element_text(size = 15), \n          axis.title=element_text(size = 18))","metadata":{"execution":{"iopub.status.busy":"2023-01-18T07:41:24.981830Z","iopub.execute_input":"2023-01-18T07:41:24.983356Z","iopub.status.idle":"2023-01-18T07:41:25.262076Z"},"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## 5.4. Plot the first two PCs","metadata":{}},{"cell_type":"code","source":"fig(25,8)\np1 <- DimPlot(seur_RNA, reduction = \"pca\", group.by = \"cell_type\") + theme_classic(base_size = 24)\np2 <- DimPlot(seur_RNA, reduction = \"pca\", group.by = \"day\") + theme_classic(base_size = 24)\np3 <- DimPlot(seur_RNA, reduction = \"pca\", group.by = \"donor\") + theme_classic(base_size = 24)\np1 + p2 + p3","metadata":{"execution":{"iopub.status.busy":"2023-01-18T07:41:32.178876Z","iopub.execute_input":"2023-01-18T07:41:32.180555Z","iopub.status.idle":"2023-01-18T07:41:34.862007Z"},"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"pca_color <- seur_RNA@reductions$pca@cell.embeddings %>%\n  as.data.frame %>% \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\nfig(22,12)\n\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='yellow', high='darkblue') +\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='yellow', high='darkblue') +\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='yellow', high='darkblue') +\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='yellow', high='darkblue') +\n      theme(aspect.ratio = 1)\n\n\nplt1 + plt2 + plt3 + plt4","metadata":{"execution":{"iopub.status.busy":"2023-01-18T07:41:40.878511Z","iopub.execute_input":"2023-01-18T07:41:40.880111Z","iopub.status.idle":"2023-01-18T07:41:45.307582Z"},"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## 5.5 Find top proteins correlated with PC_2 and their RNA","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, ] #uncomment touse xx% of cells for explaratory analysis\ngc()\n# be sure that there are the same cells in the same order in the two matrices\n#data.frame(\"PC_2_cells\" = names(PC_2), \"Prot_cells\" = rownames(mat_prot), RNA_cells = rownames(mat_RNA)) %>% head","metadata":{"_kg_hide-output":true,"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# calculate correlations\nprot_PC_2_corr <- \n   cor(mat_prot, PC_2, method = \"spearman\") %>% as.data.frame %>%\n   rename(corr_coeff = `V1`) %>%\n   mutate(Protein = rownames(.))\n   \nprot_PC_2_corr <- prot_PC_2_corr[order(prot_PC_2_corr$corr_coeff, decreasing = TRUE), ]\n\ntop20abs <- rbind(head(prot_PC_2_corr, n = 10),\n                  tail(prot_PC_2_corr, n = 10) )\n\ncat(\"Top 10 proteins positively correlated with PC_2\")\nhead(top20abs, n =10)\n\ncat(\"Top 10 proteins negatively correlated with PC_2\")\ntail(top20abs, n =10)","metadata":{"execution":{"iopub.status.busy":"2023-01-18T07:43:42.014084Z","iopub.execute_input":"2023-01-18T07:43:42.015651Z","iopub.status.idle":"2023-01-18T07:43:43.333666Z"},"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# extract Ensembl IDs to identify the RNA of the top 20 proteins and CD34 RNA\nensembl <- read.csv('/kaggle/input/research-project-01-around-multimodal-singlecell/CD_to_EnsemblSymbol_correspondence_NIPS2022.csv',row.names = 1)\nensembl <- ensembl[ensembl$CD_name %in% top20abs$Protein, ]\nensembl$CD_name <- factor(ensembl$CD_name, levels = top20abs$Protein)\nensembl <- ensembl[order(ensembl$CD_name), ]\nensembl","metadata":{"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# make long to plot with for loop all proteins\npca_color_long <- merge(pca_color, cbind(mat_prot[, rownames(top20abs)],\n                                        as.matrix(mat_RNA[, grepl(paste0(c(ensembl$ID.in.RNA.dataset, \"CD34\"), collapse = \"|\"), colnames(mat_RNA))])),\n                        by = 0)\n\npca_color_long <- pivot_longer(pca_color_long, cols = grep(\"CD|FceRIa|ENSG\", names(pca_color_long)),\n                          names_to = \"Protein\", values_to = \"Expression\")\n\npca_color_long <- merge(pca_color_long, top20abs, by = \"Protein\", all.x = TRUE)\npca_color_long <- pca_color_long %>% \n   mutate(Prot_cor = ifelse(!grepl(\"ENSG\", pca_color_long$Protein),\n                           paste0(Protein, \" (\", round(corr_coeff, 2), \")\"),\n                          pca_color_long$Protein) )\n#head(pca_color_long)\n#pca_color_long %>% group_by(Prot_cor) %>% summarise(n = n())","metadata":{"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"For proteins correlated with PC_2 we expect to see staining as with PC_2 itself (or reverse colors for negatively correlated proteins). Since the data were collected from mobilized peripheral CD34+\\nhematopoietic stem and progenitor cells (HSPCs), we also plot CD34 RNA.","metadata":{}},{"cell_type":"code","source":"fig(20,8)\nplt1 <- ggplot(pca_color, aes(x = PC_1, y = PC_2, col = PC_2) ) +\n      geom_point(size=0.8) +\n      theme_classic(base_size = 24) +\n      scale_color_gradient(low='yellow', high='darkblue') +\n      theme(aspect.ratio = 1) +\n      ggtitle(\"Reference staining\")\n\nplt2 <- ggplot(pca_color_long[grepl(\"CD34\", pca_color_long$Protein), ], aes(x = PC_1, y = PC_2, col = Expression) ) +\n      geom_point(size=0.8) +\n      theme_classic(base_size = 24) +\n      scale_color_gradient(low='yellow', high='darkblue') +\n      theme(aspect.ratio = 1) +\n      ggtitle(\"CD34 RNA\")\n\nplt1 + plt2","metadata":{"execution":{"iopub.status.busy":"2023-01-18T08:41:50.724871Z","iopub.execute_input":"2023-01-18T08:41:50.727642Z","iopub.status.idle":"2023-01-18T08:41:53.551284Z"},"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig(22,8)\n\nfor (p in ensembl$CD_name) {\n    \n    RNA <- ensembl[ensembl$CD_name == p, ]$ID.in.RNA.dataset    \n    \n    plt1 <- ggplot(pca_color_long[pca_color_long$Protein == p, ], aes(x = PC_1, y = PC_2, col = Expression) ) +\n      geom_point(size=0.8) +\n      theme_bw(base_size = 24) +\n      scale_color_gradient(low='yellow', high='darkblue') +\n      theme(aspect.ratio = 1) +\n      facet_wrap( ~ Prot_cor)\n    \n    plt2 <-  ggplot(pca_color_long[pca_color_long$Protein == RNA, ], aes(x = PC_1, y = PC_2, col = Expression) ) +\n      geom_point(size=0.8) +\n      theme_bw(base_size = 24) +\n      scale_color_gradient(low='yellow', high='darkblue') +\n      theme(aspect.ratio = 1) +\n      facet_wrap( ~ Prot_cor)\n    \n    print(plt1 + plt2)   \n}","metadata":{"execution":{"iopub.status.busy":"2023-01-18T08:31:12.850354Z","iopub.execute_input":"2023-01-18T08:31:12.851938Z","iopub.status.idle":"2023-01-18T08:32:01.569129Z"},"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# 6. What is the dimensionality of the dataset? (all parts below are unfinished)","metadata":{}},{"cell_type":"markdown","source":"As mentioned above, our dataset has a high dimensionality (each gene is a dimension). But this is way too complex and not meaningful for answering a question like: How many clusters of cells are there in the dataset?\nThis is why we try to find the true dimensionality of the dataset. There are several methods available and usually we have to test more than one and compare the results to determine the dimensionality.","metadata":{}},{"cell_type":"markdown","source":"## 6.1. Heatmap","metadata":{}},{"cell_type":"markdown","source":"In particular DimHeatmap allows for easy exploration of the primary sources of heterogeneity in a dataset, and can be useful when trying to decide which PCs to include for further downstream analyses. Both cells and features are ordered according to their PCA scores. Setting cells to a number plots the 'extreme' cells on both ends of the spectrum, which dramatically speeds plotting for large datasets. Though clearly a supervised analysis, we find this to be a valuable tool for exploring correlated feature sets.\n\nWe start with plotting heatmaps for the first 15 PCs. The DimHeatmap() function will group the cells by their gene erpression. That means that if the PC clearly divides the cells into distinct groups, the borders between the clusters are clearly defined by the colors. When the borders start to become blurry, this is one signal to \"cut off\" higher dimensions. The true dimensionality is decided by the user and it might have to be reconsidered where to draw the line. Don't be afraid to try out different \"cut off\" PC values.","metadata":{}},{"cell_type":"code","source":"fig(16,22)\nDimHeatmap(var_features, dims = 1:15, cells = 500, balanced = TRUE)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## 6.2. Jackstraw method","metadata":{}},{"cell_type":"markdown","source":"From the [tutorial](https://satijalab.org/seurat/archive/v3.2/pbmc3k_tutorial.html): To overcome the extensive technical noise in any single feature for scRNA-seq data, Seurat clusters cells based on their PCA scores, with each PC essentially representing a 'metafeature' that combines information across a correlated feature set. The top principal components therefore represent a robust compression of the dataset. However, how many componenets should we choose to include? 10? 20? 100?\n\nIn [Macosko et al](http://www.cell.com/abstract/S0092-8674(15)00549-8), we implemented a resampling test inspired by the JackStraw procedure. We randomly permute a subset of the data (1% by default) and rerun PCA, constructing a 'null distribution' of feature scores, and repeat this procedure. We identify 'significant' PCs as those who have a strong enrichment of low p-value features.\n\nTo perform the JackStraw method, we use the JackStraw() function. This function takes the Seurat object and two additional arguments: num.replicate (how many times you wish to permute the original PCA), and prop.freq (that proportion of the PCA scores you wish to permute randomly at each replicate step). We then score our permutes PCA with ScoreJackStraw() for the first 15 PCs and also plot the result with JackStrawplot().\nNOTE: This method is very computationally intensive and you should really consider twice if you want to use it on a large dataset!","metadata":{}},{"cell_type":"code","source":"# NOTE: This process can take a long time for big datasets, comment out for expediency. More\n# approximate techniques such as those implemented in ElbowPlot() can be used to reduce\n# computation time\ntic()\nvar_features <- JackStraw(var_features, num.replicate = 100, prop.freq = 0.05)\nvar_features <- ScoreJackStraw(var_features, dims = 1:20)\ntoc()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig(15,5)\nJackStrawPlot(var_features, dims = 1:20)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## 6.3. Elbow Plot method","metadata":{}},{"cell_type":"markdown","source":"An alternative heuristic method generates an 'Elbow plot': a ranking of principle components based on the percentage of variance explained by each one (ElbowPlot function).\n\nEssentially, it plots the % that each PC contributes to explaining the variation in the dataset. This means the higher the percentage, the more important the PC is. Typically, the first few PCs tend to explain quite a proportion of the variance which quickly drops and then is followed by the elbow of the plot as it turns into a more or less horizontal line toward the higher PCs.","metadata":{}},{"cell_type":"code","source":"fig(10,8)\nElbowPlot(var_features, ndims = 50)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# 7. Cluster the cells","metadata":{}},{"cell_type":"markdown","source":"Now that we determined the dimensionality of the dataset, we can use that knowledge to feed the clustering algorithm. The clustering of the cells is done with the functions FindNeighbors() and FindClusters() in Seurat. First, use FindNeighbors(): it takes a Seurat object and the proposed clusters (= dimensions) underlying the dataset. A K Nearest Neighbors graph is constructed based on the euclidean distance between points (=cells) in the PCA assuming there are n dimensions. The output from this function is the input of the FindClusters() function. This function applies modularity optimization to iteratively group cells together. You can choose between four different alogrithms to cluster the cells by. I will go with the Louvain algorithm which is the default. This is an iterative algorithm and it can therefore take some time to cluster the cells, espacially if you have large datasets.","metadata":{}},{"cell_type":"code","source":"fig(5,5)\nvar_features <- FindNeighbors(var_features, dims = 1:20)\nvar_features <- FindClusters(var_features, resolution = 0.5)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Look at cluster IDs of the first 5 cells\nhead(Idents(var_features), 10)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"If we now plot our PCA again, it should take up the proposed number of clusters that FindClusters() has identified.","metadata":{}},{"cell_type":"code","source":"fig(10,7)\n# PCA plot again with clusters encoded as colors\nDimPlot(var_features, reduction = \"pca\", dims = c(1,2))","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"","metadata":{}},{"cell_type":"markdown","source":"## 7.1. Finding differentially expressed features (cluster biomarkers)","metadata":{}},{"cell_type":"markdown","source":"Now that we found clusters, we can answer other research questions. These might be:\nWhat genes are differentially expressed in certain clusters compared to other clusters?\nOr: What genes define a cluster of cells? These genes are also called markers of a cluster.\n\nThe function we use has the surprising name FindMarkers(). We specify which cluster is of interest with the argument ident.1 and provide a precentage treshold min.pct. For example if we wanted to find the top 5 markers of the distant cluster 3:","metadata":{}},{"cell_type":"markdown","source":"I set the min.pct = 0.3 which means that the features above appear in at least 30% of all cells in all clusters and are therefore informative. You can set it to zero, but that might not be smart depending on your dataset and it increases the computing time. Also too high is not good as you might lose important information determining a cluster. The p_val is derived from a Wilcoxon Rank Sum test by default, but there are several tests available to determine differential expression.\nThe avg_logFC is the average of the log Fold change in gene expression of the cluster-of-interest (cluster 3 here) compared to all other clusters. This means that the gene MAG is on average...\n\ne^4.3=73 \n\n...fold higher expressed in cluster 3 compared to the other clusters.\n\n\nOf course, we can also compare one cluster to another like:","metadata":{}},{"cell_type":"code","source":"# find all markers of cluster 3\ncluster1.markers <- FindMarkers(var_features, ident.1 = 3, min.pct = 0.3)\nhead(cluster1.markers, n = 5)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# find all markers distinguishing cluster 5 from clusters 0 and 3\ncluster5.markers <- FindMarkers(var_features, ident.1 = 5, ident.2 = c(0, 3), min.pct = 0.25)\nhead(cluster5.markers, n = 5)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# find markers for every cluster compared to all remaining cells, report only the positive ones\nvar_features.markers <- FindAllMarkers(var_features, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.25)\nvar_features.markers %>% group_by(cluster) %>% top_n(n = 2, wt = avg_log2FC)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"cluster1.markers <- FindMarkers(var_features, ident.1 = 0, logfc.threshold = 0.25, test.use = \"roc\", only.pos = TRUE)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig(15, 7)\nVlnPlot(var_features, features = c(\"ENSG00000135218-CD36\", \"ENSG00000197956-S100A6\"))","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# you can plot raw counts as well\nfig(15, 7)\nVlnPlot(var_features, features = c(\"ENSG00000135218-CD36\", \"ENSG00000197956-S100A6\"), slot = \"counts\", log = TRUE)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig(15, 7)\nFeaturePlot(var_features, features = c(\"ENSG00000135218-CD36\", \"ENSG00000197956-S100A6\"))","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig(22, 10)\ntop10 <- var_features.markers %>% group_by(cluster) %>% top_n(n = 10, wt = avg_log2FC)\nDoHeatmap(var_features, features = top10$gene) + NoLegend()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## 7.2. Assigning cell type identity to clusters","metadata":{}},{"cell_type":"code","source":"fig(10,7)\nnew.cluster.ids <- c(\"Naive CD4 T\", \"Memory CD4 T\", \"CD14+ Mono\", \"B\", \"CD8 T\", \"FCGR3A+ Mono\", \n    \"NK\", \"DC\", \"Platelet\", \"X\")\nnames(new.cluster.ids) <- levels(var_features)\nvar_features <- RenameIdents(var_features, new.cluster.ids)\nDimPlot(var_features, reduction = \"umap\", label = TRUE, pt.size = 0.5) + NoLegend()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#saveRDS(var_features, file = \"../output/var_features3k_final.rds\")","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# 8. Run non-linear dimensional reduction (UMAP/tSNE)","metadata":{}},{"cell_type":"markdown","source":"Another non-linear technique to reducing the dimensionality of a dataset is called Uniform Manifold Approximation and Projection (UMAP). The alogrithm is fairly new.","metadata":{}},{"cell_type":"code","source":"# If you haven't installed UMAP, you can do so via reticulate::py_install(packages =\n# 'umap-learn')\nvar_features <- RunUMAP(var_features, dims = 1:10)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig(7,7)\n# note that you can set `label = TRUE` or use the LabelClusters function to help label\n# individual clusters\nDimPlot(var_features, reduction = \"umap\")","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#saveRDS(var_features, file = \"../output/var_features_tutorial.rds\")","metadata":{},"execution_count":null,"outputs":[]}]}