Adapted from https://www.broadinstitute.org/center-cell-circuits-computational-genomics-workshop by Inma Barrasa, BaRC, Whitehead Institute for Biomedical Research.

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")
 # dir.create("userName")   #Use your userName (eg. jdoe)
  1. Open the file with exercises and save it in your folder:
 #File -> Open File... -> "SingleCell_Seurat_2023.Rmd"
 
 #use "Save As..." to save the Rmd to your directory/folder
 #cd to your folder
 #setwd("/nfs/BaRC_training/userName")

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.

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.

library(Seurat)
library(dplyr)
library(Matrix)
library(gdata)
library(reshape)
library(ggplot2)
packageDescription("Seurat")$Version
## [1] "4.3.0"

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).

We will use Seurat to analyze this scRNA-seq sample:

http://satijalab.org/seurat/

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 contains various “slots” 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 slots 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")
## Warning: Feature names cannot have underscores ('_'), replacing with dashes
## ('-')

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

GetAssayData(object = seurat, slot = "counts")[1:10,1:10]
## 10 x 10 sparse Matrix of class "dgCMatrix"
##   [[ suppressing 10 column names 'AAACCTGAGCTAGTCT-1', 'AAACCTGAGGGCACTA-1', 'AAACCTGAGTACGTTC-1' ... ]]
##                                  
## 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 slot named meta.data, using the function “PercentageFeatureSet”. Then we plot the mit 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()
CombinePlots(plots = list(plot1, plot2))
## Warning: CombinePlots is being deprecated. Plots should now be combined using
## the patchwork system.

FeatureScatter is typically used to visualize gene-gene relationships, but can be used for anything calculated by the object, i.e. columns in , PC scores etc. 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 'Assay' [package "SeuratObject"] with 8 slots
##   .. .. .. ..@ 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
##   .. .. .. .. .. .. ..$ : chr [1:20213] "FO538757.2" "AP006222.2" "RP4-669L17.10" "RP11-206L10.9" ...
##   .. .. .. .. .. .. ..$ : chr [1:7109] "AAACCTGAGCTAGTCT-1" "AAACCTGAGGGCACTA-1" "AAACCTGAGTACGTTC-1" "AAACCTGAGTCCGGTC-1" ...
##   .. .. .. .. .. ..@ x       : num [1:13732088] 1 2 7 1 1 3 1 1 1 1 ...
##   .. .. .. .. .. ..@ factors : list()
##   .. .. .. ..@ data         :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
##   .. .. .. .. .. .. ..$ : chr [1:20213] "FO538757.2" "AP006222.2" "RP4-669L17.10" "RP11-206L10.9" ...
##   .. .. .. .. .. .. ..$ : chr [1:7109] "AAACCTGAGCTAGTCT-1" "AAACCTGAGGGCACTA-1" "AAACCTGAGTACGTTC-1" "AAACCTGAGTCCGGTC-1" ...
##   .. .. .. .. .. ..@ x       : num [1:13732088] 1 2 7 1 1 3 1 1 1 1 ...
##   .. .. .. .. .. ..@ factors : list()
##   .. .. .. ..@ scale.data   : num[0 , 0 ] 
##   .. .. .. ..@ key          : chr "rna_"
##   .. .. .. ..@ assay.orig   : NULL
##   .. .. .. ..@ var.features : logi(0) 
##   .. .. .. ..@ meta.features:'data.frame':   20213 obs. of  0 variables
##   .. .. .. ..@ misc         : list()
##   ..@ 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] 4 1 3
##   ..@ 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 throught 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)
## When using repel, set xnudge and ynudge to 0 for optimal results
#CombinePlots(plots = list(plot1, plot2))
plot1
## Warning: Transformation introduced infinite values in continuous x-axis

plot2
## Warning: Transformation introduced infinite values in continuous x-axis

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)
## Centering and scaling data matrix

Perform linear dimensional reduction

seurat <- RunPCA(seurat, features = VariableFeatures(object = seurat))
## PC_ 1 
## Positive:  KRT7, TACSTD2, SLPI, FXYD3, CLDN7, S100P, MPZL2, PTGES, KRT19, SERINC2 
##     ELF3, DSP, C19orf33, KRT17, LCN2, DDR1, KLF5, GPRC5A, KRT18, MDK 
##     S100A14, S100A16, ID1, TRIM29, IFITM3, TMPRSS4, CLDN4, SPINT1, SERPINB3, PERP 
## Negative:  CD74, HLA-DPA1, HLA-DPB1, HLA-DQA1, HLA-DRB1, HLA-DRA, HERPUD1, IGHA1, TNF, PIM2 
##     ITM2A, CCL5, CD7, KIAA0125, PLAC8, GZMA, STAG3, KLRB1, AC079767.4, IGKC 
##     DERL3, CTSW, C12orf75, BATF, JCHAIN, GZMK, CST7, HOPX, HLA-DRB5, TRBC1 
## PC_ 2 
## Positive:  AIF1, LYZ, TYROBP, FCER1G, SERPINA1, FCGR2A, MS4A6A, CD14, CD68, CSF1R 
##     CPVL, FCGRT, FPR1, MAFB, SPI1, CLEC7A, CD163, RNF130, PILRA, C1orf162 
##     IGSF6, PLXDC2, CST3, LILRB2, SLC11A1, FTL, NPC2, LST1, MNDA, NCF2 
## Negative:  RPS18, PERP, AQP3, KRT7, KRT17, CD24, FXYD3, DSP, C19orf33, SLPI 
##     S100A2, KRT18, PTGES, KLF5, S100P, S100A14, LCN2, S100A16, CDKN2A, ELF3 
##     TRIM29, KRT19, DDR1, SERPINB3, TACSTD2, GPRC5A, KRT6A, RAB25, MDK, TFAP2A 
## PC_ 3 
## Positive:  TPSD1, CPA3, TPSAB1, TPSB2, MS4A2, GATA2, HPGDS, SLC18A2, KIT, HDC 
##     RP11-354E11.2, IL1RL1, CLU, VWA5A, MAOB, CNRIP1, RP11-32B5.7, HPGD, LTC4S, RGS13 
##     SLC45A3, GMPR, NTRK1, GALC, TMEM176B, ENPP3, SDPR, CALB2, CD9, TIMP3 
## Negative:  S100A8, S100A12, S100A9, SPRR1B, HLA-DRA, PDZK1IP1, NAMPT, CD74, SPRR2A, TMPRSS11D 
##     HLA-DRB1, RDH10, ADGRF1, ERO1A, SCEL, FCN1, CSTA, PRSS22, F3, HLA-DPA1 
##     APOBEC3A, MUC4, RAB31, RNF39, VCAN, HLA-DPB1, TMPRSS4, BCL2A1, CLEC4E, ELF3 
## PC_ 4 
## Positive:  DHRS9, CPA3, TPSD1, TPSAB1, SPRR1B, TPSB2, MS4A2, GATA2, HDC, SPRR2A 
##     HPGDS, SLC18A2, IL1RL1, PDZK1IP1, RDH10, RP11-354E11.2, KIT, ADGRF1, HPGD, SCEL 
##     CNRIP1, RNF39, VWA5A, MUC4, ERO1A, TMPRSS11D, PRSS22, RGS13, RP11-32B5.7, MAOB 
## Negative:  SFRP1, MSLN, GAS6, FBLN1, PCSK1, LINC00473, IGFBP4, IGFBP2, FOLR1, IGFBP7 
##     BCAM, TNS4, SOD3, TUSC3, CPE, RAMP1, ADIRF, ANXA8, RAET1G, EPCAM 
##     HTRA1, INHA, PON3, RP11-49I11.1, SMOC1, FEZ1, ARG2, NNMT, CAV2, ANXA8L1 
## PC_ 5 
## Positive:  SLC11A1, VCAN, TREM1, MCEMP1, LINC01272, AQP9, CD7, PLIN2, SPP1, C5AR1 
##     SDC2, ANGPTL4, CXCL8, ANPEP, GZMA, CTSL, UPP1, MMP19, NKG7, CCL5 
##     THBS1, KLRB1, PLAUR, RGCC, PRF1, CTSW, FCGR3A, CD300E, GPNMB, FCN1 
## Negative:  HLA-DQA1, HLA-DPB1, HLA-DPA1, HLA-DQB2, CD1E, S100B, HLA-DRB1, HLA-DRB5, HLA-DRA, FCER1A 
##     PKIB, CD1A, CD74, PLD4, GSN, NDRG2, SERPINF1, CLEC10A, CD1C, APOL4 
##     KCNMB1, CD207, ENPP2, HLA-DQA2, PAK1, AXL, SERPINF2, CEACAM4, SIGLEC10, P2RY13
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, KRT7, KRT17 
## PC_ 3 
## Positive:  TPSD1, CPA3, TPSAB1, TPSB2, MS4A2 
## Negative:  S100A8, S100A12, S100A9, SPRR1B, HLA-DRA 
## PC_ 4 
## Positive:  DHRS9, CPA3, TPSD1, TPSAB1, SPRR1B 
## Negative:  SFRP1, MSLN, GAS6, FBLN1, PCSK1 
## PC_ 5 
## Positive:  SLC11A1, VCAN, TREM1, MCEMP1, LINC01272 
## Negative:  HLA-DQA1, HLA-DPB1, HLA-DPA1, HLA-DQB2, CD1E
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.

The JackStrawPlot function provides a visualization tool for comparing the distribution of p-values for each PC with a uniform distribution (dashed line). ‘Significant’ PCs will show a strong enrichment of features with low p-values (solid curve above the dashed line). This process takes a while, so instead of doing 100 replicates, we do 10 replicates.

#Uncoment below if you want to run 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.
#seurat <- JackStraw(seurat, num.replicate = 100)
#seurat <- ScoreJackStraw(seurat, dims = 1:20)
#JackStrawPlot(seurat, dims = 1:15)

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)
## Warning: The default method for RunUMAP has changed from calling Python UMAP via reticulate to the R-native UWOT using the cosine metric
## To use Python UMAP via reticulate, set umap.method to 'umap-learn' and metric to 'correlation'
## This message will be shown once per session
## 11:49:58 UMAP embedding parameters a = 0.9922 b = 1.112
## 11:49:59 Read 6122 rows and found 10 numeric columns
## 11:49:59 Using Annoy for neighbor search, n_neighbors = 30
## 11:49:59 Building Annoy index with metric = cosine, n_trees = 50
## 0%   10   20   30   40   50   60   70   80   90   100%
## [----|----|----|----|----|----|----|----|----|----|
## **************************************************|
## 11:49:59 Writing NN index file to temp file /lab/solexa_scratch/ibarrasa/tmpRstudio//RtmpQrkL8V/file224be236becde8
## 11:49:59 Searching Annoy index using 1 thread, search_k = 3000
## 11:50:01 Annoy recall = 100%
## 11:50:02 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 30
## 11:50:03 Initializing from normalized Laplacian + noise (using irlba)
## 11:50:03 Commencing optimization for 500 epochs, with 252596 positive edges
## 11:50:11 Optimization finished
DimPlot(seurat, reduction = "umap")  + NoLegend()

Cluster the cells

We use a graph-based community detection algorithm to cluster cells into cell subsets. Seurat v3 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, Seurat v3 approach to partioning the cellular distance matrix into clusters has dramatically improved. See https://satijalab.org/seurat/v3.1/pbmc3k_tutorial.html

## 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)
## Computing nearest neighbor graph
## Computing SNN
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: 224626
## 
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.9025
## 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()

CombinePlots(plots = list(p1, p2))
## Warning: CombinePlots is being deprecated. Plots should now be combined using
## the patchwork system.

#save the plots to a pdf file
pdf("./TSNE_UMAPByCluster.pdf", w=11, h=8.5)
CombinePlots(plots = list(p1, p2))
## Warning: CombinePlots is being deprecated. Plots should now be combined using
## the patchwork system.
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)
## Warning: Cannot find identity NA
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

How did this work? We took the original cluster identities (originally by default numbered) and mapped them to new values using the mapvalues command.

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)
## Calculating cluster B cells
## Calculating cluster CD4+ T cells
## Calculating cluster 2
## Calculating cluster CD8+ T cells
## Calculating cluster Monocytes/Macrophages
## Calculating cluster NK cells
## Calculating cluster Treg cells
## Calculating cluster 9
## Calculating cluster 10
## Calculating cluster 11
## Calculating cluster 12
## Calculating cluster 13

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   3.458122 0.962 0.052         0 B cells     MS4A1
## CD79A         0   3.030606 0.963 0.093         0 B cells     CD79A
## BANK1         0   2.397567 0.839 0.059         0 B cells     BANK1
## CD19          0   1.937300 0.672 0.042         0 B cells      CD19
## LINC00926     0   1.873539 0.635 0.029         0 B cells LINC00926
## CD79B         0   1.842993 0.711 0.114         0 B cells     CD79B

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: 33 × 7
## # Groups:   cluster [11]
##        p_val avg_log2FC pct.1 pct.2 p_val_adj cluster      gene     
##        <dbl>      <dbl> <dbl> <dbl>     <dbl> <fct>        <chr>    
##  1 0              2.40  0.839 0.059 0         B cells      BANK1    
##  2 0              1.94  0.672 0.042 0         B cells      CD19     
##  3 0              1.87  0.635 0.029 0         B cells      LINC00926
##  4 0              1.49  0.438 0.029 0         CD4+ T cells CD40LG   
##  5 2.40e-186      1.13  0.362 0.056 4.86e-182 CD4+ T cells TCF7     
##  6 8.64e-180      0.662 0.996 0.995 1.75e-175 CD4+ T cells TPT1     
##  7 0              5.72  0.912 0.107 0         2            KRT17    
##  8 0              5.42  0.892 0.083 0         2            S100A2   
##  9 0              5.14  0.724 0.02  0         2            LCN2     
## 10 4.97e- 69      0.844 0.263 0.057 1.01e- 64 CD8+ T cells THEMIS   
## # … with 23 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