Adapted from the material from ‘Center for Cell Circuits Computational Genomics Workshop’ by Inma Barrasa, BaRC, Whitehead Institute for Biomedical Research. Further updated to compatible with Seurat version 5 by Xinlei Gao.

Follow these steps:

  1. Go to this address in your browser and login to https://rstudio.wi.mit.edu
  2. Type the following commands:
 #setwd("/nfs/BaRC_training/SingleCellRNA-seq")
 #userName = 'userName'  #Use your own userName (eg. jdoe)
 #dir.create(userName)
  1. Open the file with exercises and save it in your folder:
 #File -> Open File... -> "SingleCell_Seurat_2024.Rmd"
 
 #use "Save As..." to save the Rmd to your directory/folder

Load necessary packages

When loading libraries, we are asking R to load code for us written by someone else. It is a convenient way to leverage and reproduce methodology developed by others.

We will use Seurat to analyze this scRNA-seq sample: http://satijalab.org/seurat/

In 2023, Seurat has been updated from version 4 to version 5. There are some changes on data structure. In this tutorial, we will use Seurat 5 to demonstrate its usage in single-cell RNA-seq analysis. Seurat5 will become the default version on the slurm cluster and RStudio server at Whitehead Institute soon.

Notably, for those who are familiar with early versions of Seurat and have been using Seurat4 for their projects, there is also a way to get back to Seurat4.

#set the R library path to the new one: 2023q4
set_lib_paths <- function(lib_vec) {
  lib_vec <- normalizePath(lib_vec, mustWork = TRUE)
  shim_fun <- .libPaths
  shim_env <- new.env(parent = environment(shim_fun))
  shim_env$.Library <- character()
  shim_env$.Library.site <- character()
  environment(shim_fun) <- shim_env
  shim_fun(lib_vec)
}

set_lib_paths(c("/nfs/apps/lib/R/4.2-focal/site-library.2023q4", "/opt/R/4.2.1/lib/R/library"))
#load Seurat 5 and other required packages from new R library '2023q4'
library(Seurat)
library(dplyr)
library(Matrix)
library(gdata)
library(reshape)
library(ggplot2)
packageDescription("Seurat")$Version
## [1] "5.0.1"
# an alternative way of loading Seurat 5 and other required packages
# require(Matrix, lib.loc="/nfs/apps/lib/R/4.2-focal/site-library.2023q4/")
# require(SeuratObject, lib.loc="/nfs/apps/lib/R/4.2-focal/site-library.2023q4/")
# require(Seurat, lib.loc="/nfs/apps/lib/R/4.2-focal/site-library.2023q4/")
# library(dplyr)
# library(Matrix)
# library(gdata)
# library(reshape)
# library(ggplot2)
# packageDescription("Seurat")$Version

Goal

Introduction

Data produced in a single cell RNA-seq experiment has several interesting characteristics that make it distinct from data produced in a bulk population RNA-seq experiment. Two characteristics that are important to keep in mind when working with scRNA-Seq are drop-out (the excessive amount of zeros due to limiting mRNA) and the potential for quality control (QC) metrics to be confounded with biology. This combined with the ability to measure heterogeniety from cells in samples has shifted the field away from the typical analysis in population-based RNA-Seq. Here we demonstrate some approaches to quality control, followed by identifying and analyzing cell subsets.

For this tutorial, we will be analyzing the a dataset of Non-Small Cell Lung Cancer Cells (NSCLC) freely available from 10X Genomics (https://support.10xgenomics.com/single-cell-vdj/datasets/2.2.0/vdj_v1_hs_nsclc_5gex), using the Seurat R package (http://satijalab.org/seurat/), a popular and powerful set of tools to conduct scRNA-seq analysis in R. In this dataset, there are 7802 single cells that were sequenced on the Illumina NovaSeq 6000. Please note this tutorial borrows heavily from Seurat’s tutorials (https://satijalab.org/seurat/vignettes.html), so feel free to go through them in more detail.

Read in NSCLC counts matrix.

The data for Non-Small Cell Lung Cancer Cells (NSCLC) is freely available from 10X Genomics (https://support.10xgenomics.com/single-cell-vdj/datasets/2.2.0/vdj_v1_hs_nsclc_5gex).

One very popular product for 3’ single cell transcriptomics sequencing is the 10X assay. Sequencing data from 10X can be transformed to gene count matrices using the Cell Ranger pipeline. Here we go through many of the steps we saw earlier but using the Seurat R library. The Seurat R library makes available many useful functionalities. We will be using the same data as before.

We start by reading in the counts matrix generated by the Cell Ranger count program.

counts_matrix_filename = "/nfs/BaRC_training/SingleCellRNA-seq/filtered_gene_bc_matrices/GRCh38/"
counts <- Read10X(data.dir = counts_matrix_filename)

Seurat works with the count matrix in the sparse format, which is far more memory efficient. Almost all our analysis will be on the single object, of class Seurat, created above. This object can contain 1 or several assays and each assay contains various “layers” (which is called “slots” in earlier versions). that will store not only the raw count data, but also the results from various computations below. This has the advantage that we do not need to keep track of inidividual variables of interest - they can all be collapsed into a single object as long as these “layers” are pre-defined. This object was made with the CreateSeuratObject(). Here we are filtering out genes that are are expressed in 2 or fewer cells and are filtering cells with complexity less than 350 genes.

seurat <- CreateSeuratObject(counts, min.cells = 3, min.features = 350, project = "10X_NSCLC")

The “LayerData” function allows us to get the counts/raw data stored on the “counts” layer. As we proceed on our analysis other layers like “data” (normalized data) and “scale.data” (data scaled by gene) will be populated.

LayerData(seurat, assay="RNA", layer='counts')[1:10,1:10]
## 10 x 10 sparse Matrix of class "dgCMatrix"
##                                  
## FO538757.2    . . . . 1 2 . . . .
## AP006222.2    . . . . . . . . . .
## RP4-669L17.10 . . . . . . . . . .
## RP11-206L10.9 . . . . . . . . . .
## LINC00115     . . . . . . . . . .
## FAM41C        . . . . . . . . . .
## SAMD11        . . . . . 1 . . . .
## NOC2L         . . . . 2 4 . 1 . .
## KLHL17        . . . . . . . . . .
## PLEKHN1       . . . . . . . . . .

Preprocessing step 1 : Filter out low-quality cells

The Seurat object initialization step above only considered cells that expressed at least 350 genes. Additionally, we would like to exclude cells that are damaged. A common metric to judge this (although by no means the only one) is the relative expression of mitochondrially derived genes.

Here we calculate the percent mitochondrial reads and add it to the Seurat object in the layer named meta.data, using the function “PercentageFeatureSet”. Then we plot the mitochodrial percentage using the violin plot function provided by Seurat.

# The number of genes and UMIs (nFeature_RNA and nCount_RNA) are automatically calculated for every object by Seurat.

# this is step adds the % mit content to the metadata
seurat[["percent.mt"]] <- PercentageFeatureSet(seurat, pattern = "^MT-")

head(seurat@meta.data)
##                    orig.ident nCount_RNA nFeature_RNA percent.mt
## AAACCTGAGCTAGTCT-1  10X_NSCLC       3605         1184  20.832178
## AAACCTGAGGGCACTA-1  10X_NSCLC       3828         1387   3.448276
## AAACCTGAGTACGTTC-1  10X_NSCLC       6457         1784   2.865108
## AAACCTGAGTCCGGTC-1  10X_NSCLC       3075         1092   6.276423
## AAACCTGCACCAGGTC-1  10X_NSCLC       9399         2625   3.202468
## AAACCTGCACCTCGTT-1  10X_NSCLC      48838         5602   2.321962
VlnPlot(seurat, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3, pt.size = 0.1) + NoLegend()

# FeatureScatter is typically used to visualize feature-feature relationships, but can be used
# for anything calculated by the object, i.e. columns in object metadata, PC scores etc.

plot1 <- FeatureScatter(seurat, feature1 = "nCount_RNA", feature2 = "percent.mt") + NoLegend()
plot2 <- FeatureScatter(seurat, feature1 = "nCount_RNA", feature2 = "nFeature_RNA") + NoLegend()
plot1 + plot2

Since there is a rare subset of cells with an outlier level of high mitochondrial percentage and also low UMI content, we filter these as well. The title of these plots is the Pearson correlation between the two variables: percent.mito vs nCount_RNA or nFeature_RNA (number of genes) vs nCount_RNA.

Examine contents of Seurat object

str(seurat)
## Formal class 'Seurat' [package "SeuratObject"] with 13 slots
##   ..@ assays      :List of 1
##   .. ..$ RNA:Formal class 'Assay5' [package "SeuratObject"] with 8 slots
##   .. .. .. ..@ layers    :List of 1
##   .. .. .. .. ..$ counts:Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
##   .. .. .. .. .. .. ..@ i       : int [1:13732088] 22 37 100 102 105 132 146 171 177 185 ...
##   .. .. .. .. .. .. ..@ p       : int [1:7110] 0 1184 2571 4355 5447 8072 13674 15861 17504 18605 ...
##   .. .. .. .. .. .. ..@ Dim     : int [1:2] 20213 7109
##   .. .. .. .. .. .. ..@ Dimnames:List of 2
##   .. .. .. .. .. .. .. ..$ : NULL
##   .. .. .. .. .. .. .. ..$ : NULL
##   .. .. .. .. .. .. ..@ x       : num [1:13732088] 1 2 7 1 1 3 1 1 1 1 ...
##   .. .. .. .. .. .. ..@ factors : list()
##   .. .. .. ..@ cells     :Formal class 'LogMap' [package "SeuratObject"] with 1 slot
##   .. .. .. .. .. ..@ .Data: logi [1:7109, 1] TRUE TRUE TRUE TRUE TRUE TRUE ...
##   .. .. .. .. .. .. ..- attr(*, "dimnames")=List of 2
##   .. .. .. .. .. .. .. ..$ : chr [1:7109] "AAACCTGAGCTAGTCT-1" "AAACCTGAGGGCACTA-1" "AAACCTGAGTACGTTC-1" "AAACCTGAGTCCGGTC-1" ...
##   .. .. .. .. .. .. .. ..$ : chr "counts"
##   .. .. .. .. .. ..$ dim     : int [1:2] 7109 1
##   .. .. .. .. .. ..$ dimnames:List of 2
##   .. .. .. .. .. .. ..$ : chr [1:7109] "AAACCTGAGCTAGTCT-1" "AAACCTGAGGGCACTA-1" "AAACCTGAGTACGTTC-1" "AAACCTGAGTCCGGTC-1" ...
##   .. .. .. .. .. .. ..$ : chr "counts"
##   .. .. .. ..@ features  :Formal class 'LogMap' [package "SeuratObject"] with 1 slot
##   .. .. .. .. .. ..@ .Data: logi [1:20213, 1] TRUE TRUE TRUE TRUE TRUE TRUE ...
##   .. .. .. .. .. .. ..- attr(*, "dimnames")=List of 2
##   .. .. .. .. .. .. .. ..$ : chr [1:20213] "FO538757.2" "AP006222.2" "RP4-669L17.10" "RP11-206L10.9" ...
##   .. .. .. .. .. .. .. ..$ : chr "counts"
##   .. .. .. .. .. ..$ dim     : int [1:2] 20213 1
##   .. .. .. .. .. ..$ dimnames:List of 2
##   .. .. .. .. .. .. ..$ : chr [1:20213] "FO538757.2" "AP006222.2" "RP4-669L17.10" "RP11-206L10.9" ...
##   .. .. .. .. .. .. ..$ : chr "counts"
##   .. .. .. ..@ default   : int 1
##   .. .. .. ..@ assay.orig: chr(0) 
##   .. .. .. ..@ meta.data :'data.frame':  20213 obs. of  0 variables
##   .. .. .. ..@ misc      :List of 1
##   .. .. .. .. ..$ calcN: logi TRUE
##   .. .. .. ..@ key       : chr "rna_"
##   ..@ meta.data   :'data.frame': 7109 obs. of  4 variables:
##   .. ..$ orig.ident  : Factor w/ 1 level "10X_NSCLC": 1 1 1 1 1 1 1 1 1 1 ...
##   .. ..$ nCount_RNA  : num [1:7109] 3605 3828 6457 3075 9399 ...
##   .. ..$ nFeature_RNA: int [1:7109] 1184 1387 1784 1092 2625 5602 2187 1643 1101 2717 ...
##   .. ..$ percent.mt  : num [1:7109] 20.83 3.45 2.87 6.28 3.2 ...
##   ..@ active.assay: chr "RNA"
##   ..@ active.ident: Factor w/ 1 level "10X_NSCLC": 1 1 1 1 1 1 1 1 1 1 ...
##   .. ..- attr(*, "names")= chr [1:7109] "AAACCTGAGCTAGTCT-1" "AAACCTGAGGGCACTA-1" "AAACCTGAGTACGTTC-1" "AAACCTGAGTCCGGTC-1" ...
##   ..@ graphs      : list()
##   ..@ neighbors   : list()
##   ..@ reductions  : list()
##   ..@ images      : list()
##   ..@ project.name: chr "10X_NSCLC"
##   ..@ misc        : list()
##   ..@ version     :Classes 'package_version', 'numeric_version'  hidden list of 1
##   .. ..$ : int [1:3] 5 0 1
##   ..@ commands    : list()
##   ..@ tools       : list()

These are the slots in the Seurat object. Some of the slots are automatically updated by Seurat as you move through analysis. Take a moment to look through the information, knowing the slots allow you to leverage work Seurat has already done for you.

VlnPlot(object = seurat, features = c("nFeature_RNA"), group.by = c('orig.ident'), pt.size = 0.1) +NoLegend()

Here we plot the number of genes per cell by what Seurat calls orig.ident. Identity is a concept that is used in the Seurat object to refer to the cell identity. In this case, the cell identity is 10X_NSCLC, but after we cluster the cells, the cell identity will be whatever cluster the cell belongs to. We will see how identity updates as we go through the analysis.

Next, let’s filter the cells based on the quality control metrics.

seurat <- subset(seurat, subset = nFeature_RNA > 350 & nFeature_RNA < 5000 & percent.mt < 10)

Preprocessing step 2 : Expression normalization

After removing unwanted genes cells from the dataset, the next step is to normalize the data. By default, we employ a global-scaling normalization method “LogNormalize” that normalizes the gene expression measurements for each cell by the total expression, multiplies this by a scale factor (10,000 by default), and log-transforms the result. The log-transformation is a commonly used transformation that has many desirable properties, such as variance stabilization. Seurat uses natural log.

seurat <- NormalizeData(seurat, normalization.method = "LogNormalize", scale.factor = 10000)

Detection of variable genes across the single cells

Seurat calculates highly variable genes and focuses on these for downstream analysis. FindVariableFeatures, by default, it returns 2,000 features per dataset. More detailed can be found here: https://www.ncbi.nlm.nih.gov/pubmed/31178118

seurat <- FindVariableFeatures(seurat, selection.method = "vst", nfeatures = 2000)
length(VariableFeatures(seurat))
## [1] 2000
# Identify the 10 most highly variable genes
top10 <- head(VariableFeatures(seurat), 10)

# plot variable features with and without labels
plot1 <- VariableFeaturePlot(seurat)
plot2 <- LabelPoints(plot = plot1, points = top10, repel = TRUE)
plot1

plot2

Scaling the data

Next, we apply a linear transformation (‘scaling’) that is a standard pre-processing step prior to dimensional reduction techniques like PCA. The ScaleData function:

Shifts the expression of each gene, so that the mean expression across cells is 0
Scales the expression of each gene, so that the variance across cells is 1
    This step gives equal weight to all the genes in downstream analyses, so that highly-expressed genes do not dominate
The results of this are stored in seurat[["RNA"]]@scale.data
all.genes <- rownames(seurat)
seurat <- ScaleData(seurat, features = all.genes)

Perform linear dimensionality reduction

seurat <- RunPCA(seurat, features = VariableFeatures(object = seurat))
print(seurat[["pca"]], dims = 1:5, nfeatures = 5)
## PC_ 1 
## Positive:  KRT7, TACSTD2, SLPI, FXYD3, CLDN7 
## Negative:  CD74, HLA-DPA1, HLA-DPB1, HLA-DQA1, HLA-DRB1 
## PC_ 2 
## Positive:  AIF1, LYZ, TYROBP, FCER1G, SERPINA1 
## Negative:  RPS18, PERP, AQP3, CD24, KRT7 
## PC_ 3 
## Positive:  TPSD1, CPA3, TPSAB1, TPSB2, MS4A2 
## Negative:  S100A8, S100A12, S100A9, HLA-DRA, SPRR1B 
## PC_ 4 
## Positive:  SFRP1, MSLN, GAS6, FBLN1, PCSK1 
## Negative:  DHRS9, SPRR1B, CPA3, TPSD1, TPSAB1 
## PC_ 5 
## Positive:  SLC11A1, VCAN, TREM1, MCEMP1, LINC01272 
## Negative:  HLA-DQA1, HLA-DQB2, HLA-DPB1, CD1E, HLA-DPA1
VizDimLoadings(seurat, dims = 1:2, reduction = "pca")

Visualize the PCA plot

DimPlot(seurat, reduction = "pca") + NoLegend()

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.

DimHeatmap(seurat, dims = 1:3, cells = 500, balanced = TRUE)

Now that we have performed PCA, we can plot metadata and genes on the cell groups. Here we use the FeaturePlot function to plot on the PCA. We specify plotting on the principal components using reduction.use, specifically the first two PCA components (specified in dim.1 and dim.2).

FeaturePlot(seurat,  reduction='pca', features=c("nFeature_RNA",  "percent.mt"))

How would you plot the same plot using component 2 and component 3? If we wanted to save the Seurat object we could use the save function (and later load the object with load()). We are not going to do this now, but this is one of the ways to share analysis, by saving and sharing Seurat R objects.

# save.image("./Seurat.RData")

Determine the ‘dimensionality’ of the dataset

PCA will create many components of variation but not all are useful. One way to select principal components that explain the variation in our data is to use an elbow or scree plot.

ElbowPlot(seurat)

From this plot we are going to select the first 20 principal components to use in the t-SNE visualization.

Visualization by tSNE and UMAP

tSNE is a way to reduce data dimensionality and visualize data, but it is not a clustering approach. Please read How to use tSNE effectively (https://distill.pub/2016/misread-tsne/)

UMAP is a much faster visualization than tSNE and may preserve more of the global structure in your data. The preprint for the method is available at (https://arxiv.org/abs/1802.03426).

seurat <- RunTSNE(seurat,dims = 1:20) # runs Barnes-hut t-SNE
TSNEPlot(object = seurat)  + NoLegend()

seurat <- RunUMAP(seurat, dims = 1:10)
DimPlot(seurat, reduction = "umap")  + NoLegend()

Cluster the cells

Seurat applies a graph-based clustering approach, building upon initial strategies in Macosko et al (https://www.cell.com/fulltext/S0092-8674(15)00549-8). Importantly, the distance metric which drives the clustering analysis (based on previously identified PCs) remains the same. However, its approach to partioning the cellular distance matrix into clusters has dramatically improved. See https://satijalab.org/seurat/articles/pbmc3k_tutorial.

## Find clusters

# save.SNN = T saves the SNN so that the clustering algorithm can be rerun using the same graph
# but with a different resolution value.
#Find clusters 
seurat <- FindNeighbors(object = seurat, dims = 1:20)
seurat <- FindClusters(object = seurat, resolution = 0.5)
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
## 
## Number of nodes: 6122
## Number of edges: 223642
## 
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.9036
## Number of communities: 14
## Elapsed time: 0 seconds
# note that you can set label=T to help label individual clusters
TSNEPlot(object = seurat, label=T)+ NoLegend()

Compare tSNE and UMAP visualizations

Here we compare how the clusters look when mapped onto either the tSNE or UMAP visualizations.

p1 <- TSNEPlot(object = seurat, label=T) + NoLegend()
p2 <- UMAPPlot(object = seurat, label=T) + NoLegend()

p1 + p2

#save the plots to a pdf file
pdf("./TSNE_UMAPByCluster.pdf", w=11, h=8.5)
p1 + p2
dev.off()
## png 
##   2
pdf("./TSNE_Cluster.pdf", w=11, h=8.5)
p1
dev.off()
## png 
##   2

Mapping gene expression onto cell clusters and assigning identity to clusters

We can now visualize how an individual gene’s expression changes across the different cell subsets. We use FeaturePlot in these visualizations. Plotting known marker genes, helps identifying some of the clusters.

# Plotting genes of interest
# Monocyte/Macrophage markers
FeaturePlot(object = seurat, features = c("CD14", "LYZ","FCGR3A"), cols = c("grey", "blue"), reduction = "tsne")

pdf("./TSNE_MonocyteMacroph.pdf", w=11, h=8.5)
FeaturePlot(object = seurat, features = c("CD14", "LYZ","FCGR3A"), cols = c("grey", "blue"), reduction = "tsne")
dev.off()
## png 
##   2
pdf("./UMAP_MonocyteMacroph.pdf", w=11, h=8.5)
FeaturePlot(object = seurat, features = c("CD14", "LYZ","FCGR3A"), cols = c("grey", "blue"), reduction = "umap")
dev.off()
## png 
##   2
# T cell markers 
genes1 <- c("CD3D", "CD3E", "CD4","CD8A")
#treg markers
genes2 <- c( "FOXP3", "IL7R") 
FeaturePlot(object = seurat, features = genes1, cols = c("grey", "blue"), reduction = "tsne")

FeaturePlot(object = seurat, features = genes1, cols = c("grey", "blue"), reduction = "umap")

FeaturePlot(object = seurat, features = genes2, cols = c("grey", "blue"), reduction = "tsne")

FeaturePlot(object = seurat, features = genes2, cols = c("grey", "blue"), reduction = "umap")

pdf("./T_cell_markers1.pdf", w=11, h=8.5)
FeaturePlot(object = seurat, features = genes1, cols = c("grey", "blue"), reduction = "tsne")
dev.off()
## png 
##   2
pdf("./T_cell_markers2.pdf", w=11, h=8.5)
FeaturePlot(object = seurat, features = genes2, cols = c("grey", "blue"), reduction = "tsne")
dev.off()
## png 
##   2
pdf("./T_cell_markers3.pdf", w=11, h=8.5)
FeaturePlot(object = seurat, features = genes1, cols = c("grey", "blue"), reduction = "umap")
dev.off()
## png 
##   2
pdf("./T_cell_markers4.pdf", w=11, h=8.5)
FeaturePlot(object = seurat, features = genes2, cols = c("grey", "blue"), reduction = "umap")
dev.off()
## png 
##   2
# NK cell markers
genes <- c("GZMA", "GZMB", "PRF1", "NKG7", "GNLY")
FeaturePlot(object = seurat, features = genes, cols = c("grey", "blue"), reduction = "tsne")

FeaturePlot(object = seurat, features = genes, cols = c("grey", "blue"), reduction = "umap")

pdf("./NK_cell_markers.pdf", w=11, h=8.5)
FeaturePlot(object = seurat, features = genes, cols = c("grey", "blue"), reduction = "tsne")
FeaturePlot(object = seurat, features = genes, cols = c("grey", "blue"), reduction = "umap")
dev.off()
## png 
##   2
# B cell markers
genes <- genes <- c("CD19", "MS4A1")
FeaturePlot(object = seurat, features = genes, cols = c("grey", "blue"), reduction = "tsne")

FeaturePlot(object = seurat, features = genes, cols = c("grey", "blue"), reduction = "umap")

pdf("./B_cell_markers.pdf", w=11, h=8.5)
FeaturePlot(object = seurat, features = genes, cols = c("grey", "blue"), reduction = "tsne")
FeaturePlot(object = seurat, features = genes, cols = c("grey", "blue"), reduction = "umap")
dev.off()
## png 
##   2

Labeling cell subsets

We want to label each cell subset according to what cell type it is. We can update the identity slot to these new identities.

#current.cluster.ids <- c(0:14)
new.cluster.ids <- c("B cells", "CD4+ T cells", "2", "CD8+ T cells", "Monocytes/Macrophages", "Monocytes/Macrophages", "NK cells", "Monocytes/Macrophages", "Treg cells" , "9", "10", "11", "12", "13", "14") 

names(new.cluster.ids) <- levels(seurat)
seurat <- RenameIdents(seurat, new.cluster.ids)
DimPlot(seurat, reduction = "tsne", label = TRUE, pt.size = 0.1) + NoLegend()

pdf("./NewLabels.pdf", w=11, h=8.5)
DimPlot(seurat, reduction = "tsne", label = TRUE, pt.size = 0.1) + NoLegend()
dev.off()
## png 
##   2

Using differentially expressed genes to help identify cell subsets

In the above code, we have identified several clusters of immune cells. However we also expect to see cell clusters consisting of stromal cells and/or malignant cells in this NSCLC sample. For some of these cell subsets, we may not have pre-existing marker genes that suggest the cell type. To help identify these cells, we look for genes that are differentially expressed in these clusters relative to all other clusters.

MIN_LOG2FOLD_CHANGE = 0.6 # set to minimum required average log fold change in gene expression.
MIN_PCT_CELLS_EXPR_GENE = .25  # minimum percent of cells that must express gene in either clstr.
#these values above are a little more stringent than the defaults
#you can set them to defaults, logfc.threshold = 0.25, min.pct = 0.1,  

# Here we find all markers for 
all.markers = FindAllMarkers(seurat,
                             min.pct = MIN_PCT_CELLS_EXPR_GENE,
                             logfc.threshold = MIN_LOG2FOLD_CHANGE,
                             only.pos = TRUE)

Here we found markers for all clusters using a DE test. There are other ways of selecting markers, feel free to read the original Seurat tutorial for more details. Here we use the Wilcoxon Rank Sum test, “wilcox”. There are many other options, including MAST.

Let’s look at the top markers for each cluster sorted by p-value.

# sort all the markers by p-value
all.markers.sortedByPval = all.markers[order(all.markers$p_val),]

# take a look at the top most significant markers
head(all.markers.sortedByPval)
##          p_val avg_log2FC pct.1 pct.2 p_val_adj cluster     gene
## MS4A1        0   5.221317 0.963 0.052         0 B cells    MS4A1
## CD79A        0   4.340079 0.964 0.093         0 B cells    CD79A
## BANK1        0   5.073961 0.840 0.059         0 B cells    BANK1
## HLA-DQA1     0   1.549903 0.979 0.282         0 B cells HLA-DQA1
## CD19         0   5.181863 0.672 0.042         0 B cells     CD19
## HLA-DQB1     0   1.337430 0.965 0.343         0 B cells HLA-DQB1

Make a heatmap showing the top 10 marker genes for each cluster

We can also make heat maps of the top markers for each cluster.

top10 <- all.markers.sortedByPval %>%  group_by(cluster)  %>% do(head(., n=10))
DoHeatmap(object = seurat, features = top10$gene)+ NoLegend() +
    theme(axis.text.y = element_text(color = "black", size = 5))

pdf("./HeatMaptop10Markers.pdf", w=11, h=8.5)
DoHeatmap(object = seurat, features = top10$gene)+ NoLegend() +
    theme(axis.text.y = element_text(color = "black", size = 5))
dev.off()
## png 
##   2

Based on these differentially expressed genes, can you identify additional cell subsets?

You can save the object at this point so that it can easily be loaded back into an R session without having to rerun the computationally intensive steps performed above, or easily shared with collaborators.

save.image("./Seurat.RData")
# examine the top 4 markers in the context of the tSNE plots:
FeaturePlot(seurat, features = all.markers.sortedByPval$gene[1:4], reduction = "tsne")

Get genes uniquely differentially expressed in each cluster.

In the tests above for differentially expressed genes, some differentially expressed genes may be common across clusters. For example, the 3 clusters of T cells all have IL7R (T cell development gene) as a differentially expressed gene. To find genes that are differentially expressed and unique to the cell cluster, we can use the below code.

# Find differentially expressed genes that are only differentially expressed in one cluster of cells.
genes_uniquely_DE = all.markers.sortedByPval %>% dplyr::filter(avg_log2FC >= MIN_LOG2FOLD_CHANGE) %>% group_by(gene) %>%  summarize(n=n()) %>%  dplyr::filter(n==1)

# Sort DE genes by p value.
genes_uniquely_DE.markers.sortedByPval =
  all.markers.sortedByPval[all.markers.sortedByPval$gene %in% genes_uniquely_DE$gene & 
                             all.markers.sortedByPval$avg_log2FC >= MIN_LOG2FOLD_CHANGE,]

# Get top 3 unique genes for each cell cluster.
top_marker_each = genes_uniquely_DE.markers.sortedByPval %>% dplyr::group_by(cluster) %>% do(head(., n=3))  # number of top markers for each cluster.

print(top_marker_each)
## # A tibble: 36 × 7
## # Groups:   cluster [12]
##        p_val avg_log2FC pct.1 pct.2 p_val_adj cluster      gene     
##        <dbl>      <dbl> <dbl> <dbl>     <dbl> <fct>        <chr>    
##  1 0              5.18  0.672 0.042 0         B cells      CD19     
##  2 0              5.41  0.636 0.029 0         B cells      LINC00926
##  3 0              5.03  0.543 0.034 0         B cells      TNFRSF13B
##  4 2.89e-204      3.99  0.297 0.025 5.84e-200 CD4+ T cells MAL      
##  5 7.04e-200      3.65  0.343 0.041 1.42e-195 CD4+ T cells LEF1     
##  6 1.99e-182      0.709 0.995 0.995 4.03e-178 CD4+ T cells TPT1     
##  7 0              8.04  0.86  0.02  0         2            KRT7     
##  8 0              7.57  0.895 0.083 0         2            S100A2   
##  9 0              7.58  0.825 0.015 0         2            C19orf33 
## 10 7.01e- 79      2.94  0.253 0.046 1.42e- 74 CD8+ T cells NELL2    
## # ℹ 26 more rows
DoHeatmap(object = seurat, features = top_marker_each$gene)+ NoLegend() +
    theme(axis.text.y = element_text(color = "black", size = 5))

pdf("./HeatMaptop3Markers.pdf", w=11, h=8.5)
DoHeatmap(object = seurat, features = top_marker_each$gene)+ NoLegend() +
    theme(axis.text.y = element_text(color = "black", size = 5))
dev.off()
## png 
##   2

Based on these differentially expressed genes, can you identify additional cell subsets?

Visualizing marker genes: violin plots

We can look at marker gene expression using violin plots.

#uncomment the loop below if you want to make a plot for each of the marker genes.
#for (i in 1:length(top_marker_each$gene)) {
#  print(VlnPlot(seurat, features = top_marker_each$gene[i]), pt.size = 0.1) +NoLegend()
#}
#look at the firt 5 marker genes:
for (i in 1:5) {
  print(VlnPlot(seurat, features = top_marker_each$gene[i]), pt.size = 0.1) +NoLegend()
}

Other ways to visualize marker genes: dot plots

Dotplots are a concise way of exploring percent of cells expressing a gene and the gene expression intensity. This gives much of the same information a heatmap gives us without the resolution of looking at every cell. This can shield us from some of the noise caused by drop out and scales better with larger datasets.

DotPlot(seurat, features=unique(top_marker_each$gene),  dot.scale = 6)  + RotatedAxis() +theme(axis.text.x = element_text(color = "black", size = 7))

DotPlot(seurat, features=unique(top_marker_each$gene),  dot.scale = 6)  + RotatedAxis() +theme(axis.text.x = element_text(color = "black", size = 7)) + NoLegend()

pdf("./DotPlot.pdf", w=11, h=8.5)
DotPlot(seurat, features=unique(top_marker_each$gene),  dot.scale = 6)  + RotatedAxis() +theme(axis.text.x = element_text(color = "black", size = 7))
dev.off()
## png 
##   2