Load in the data

This vignette demonstrates some useful features for interacting with the Seurat object. For demonstration purposes, we will be using the 2,700 PBMC object that is created in the first guided tutorial. You can download the pre-computed object here. To simulate the scenario where we have two replicates, we will randomly assign half the cells in each cluster to be from “rep1” and other half from “rep2”.

library(Seurat)
pbmc <- readRDS(file = "../data/pbmc3k_final.rds")

# pretend that cells were originally assigned to one of two replicates (we assign randomly here)
# if your cells do belong to multiple replicates, and you want to add this info to the Seurat
# object create a data frame with this information (similar to replicate.info below)
set.seed(42)
pbmc$replicate <- sample(c("rep1", "rep2"), size = ncol(pbmc), replace = TRUE)

Switch identity class between cluster ID and replicate

# Plot UMAP, coloring cells by cell type (currently stored in object@ident)
DimPlot(pbmc, reduction = "umap")

# How do I create a UMAP plot where cells are colored by replicate?  First, store the current
# identities in a new column of meta.data called CellType
pbmc$CellType <- Idents(pbmc)
# Next, switch the identity class of all cells to reflect replicate ID
Idents(pbmc) <- "replicate"
DimPlot(pbmc, reduction = "umap")

# alternately : DimPlot(pbmc, reduction = 'umap', group.by = 'replicate') you can pass the
# shape.by to label points by both replicate and cell type

# Switch back to cell type labels
Idents(pbmc) <- "CellType"

Tabulate cells by cluster ID, replicate, or both

# How many cells are in each cluster
table(Idents(pbmc))
## 
##  Naive CD4 T Memory CD4 T   CD14+ Mono            B        CD8 T 
##          697          483          480          344          271 
## FCGR3A+ Mono           NK           DC     Platelet 
##          162          155           32           14
# How many cells are in each replicate?
table(pbmc$replicate)
## 
## rep1 rep2 
## 1327 1311
# What proportion of cells are in each cluster?
prop.table(table(Idents(pbmc)))
## 
##  Naive CD4 T Memory CD4 T   CD14+ Mono            B        CD8 T 
##  0.264215315  0.183093252  0.181956027  0.130401820  0.102729340 
## FCGR3A+ Mono           NK           DC     Platelet 
##  0.061410159  0.058756634  0.012130402  0.005307051
# How does cluster membership vary by replicate?
table(Idents(pbmc), pbmc$replicate)
##               
##                rep1 rep2
##   Naive CD4 T   337  360
##   Memory CD4 T  243  240
##   CD14+ Mono    240  240
##   B             166  178
##   CD8 T         152  119
##   FCGR3A+ Mono   84   78
##   NK             83   72
##   DC             18   14
##   Platelet        4   10
prop.table(table(Idents(pbmc), pbmc$replicate), margin = 2)
##               
##                       rep1        rep2
##   Naive CD4 T  0.253956292 0.274599542
##   Memory CD4 T 0.183119819 0.183066362
##   CD14+ Mono   0.180859081 0.183066362
##   B            0.125094197 0.135774218
##   CD8 T        0.114544084 0.090770404
##   FCGR3A+ Mono 0.063300678 0.059496568
##   NK           0.062547099 0.054919908
##   DC           0.013564431 0.010678871
##   Platelet     0.003014318 0.007627765

Selecting particular cells and subsetting the Seurat object

# What are the cell names of all NK cells?
WhichCells(pbmc, idents = "NK")
##   [1] "AAACCGTGTATGCG" "AAATTCGATTCTCA" "AACCTTACGCGAGA" "AACGCCCTCGTACA"
##   [5] "AACGTCGAGTATCG" "AAGATTACCTCAAG" "AAGCAAGAGCTTAG" "AAGCAAGAGGTGTT"
##   [9] "AAGTAGGATACAGC" "AATACTGAATTGGC" "AATCCTTGGTGAGG" "AATCTCTGCTTTAC"
##  [13] "ACAAATTGTTGCGA" "ACAACCGAGGGATG" "ACAATTGATGACTG" "ACACCCTGGTGTTG"
##  [17] "ACAGGTACTGGTGT" "ACCTGGCTAAGTAG" "ACGAACACCTTGTT" "ACGATCGAGGACTT"
##  [21] "ACGCAATGGTTCAG" "ACGCTGCTGTTCTT" "ACGGAACTCAGATC" "ACGTGATGTGACAC"
##  [25] "ACGTTGGAGCCAAT" "ACTGCCACTCCGTC" "ACTGGCCTTCAGTG" "ACTTCAACGTAGGG"
##  [29] "AGAACAGAAATGCC" "AGATATACCCGTAA" "AGATTCCTGTTCAG" "AGCCTCTGCCAATG"
##  [33] "AGCGATTGAGATCC" "AGGATGCTTTAGGC" "AGGGACGAGTCAAC" "AGTAATACATCACG"
##  [37] "AGTCACGATGAGCT" "AGTTTGCTACTGGT" "ATACCACTGCCAAT" "ATACTCTGGTATGC"
##  [41] "ATCCCGTGCAGTCA" "ATCTTTCTTGTCCC" "ATGAAGGACTTGCC" "ATGATAACTTCACT"
##  [45] "ATGATATGGTGCTA" "ATGGACACGCATCA" "ATGGGTACATCGGT" "ATTAACGATGAGAA"
##  [49] "ATTCCAACTTAGGC" "CAAGGTTGTCTGGA" "CAATCTACTGACTG" "CACCACTGGCGAAG"
##  [53] "CACGGGTGGAGGAC" "CAGATGACATTCTC" "CAGCAATGGAGGGT" "CAGCGGACCTTTAC"
##  [57] "CAGCTCTGTGTGGT" "CAGTTTACACACGT" "CATCAGGACTTCCG" "CATCAGGATAGCCA"
##  [61] "CATGAGACGTTGAC" "CATTACACCAACTG" "CATTTCGAGATACC" "CCTCGAACACTTTC"
##  [65] "CGACCACTAAAGTG" "CGACCACTGCCAAT" "CGAGGCTGACGCTA" "CGCCGAGAGCTTAG"
##  [69] "CGGCGAACGACAAA" "CGGCGAACTACTTC" "CGGGCATGTCTCTA" "CGTACCTGGCATCA"
##  [73] "CGTGTAGACGATAC" "CGTGTAGAGTTACG" "CGTGTAGATTCGGA" "CTAAACCTCTGACA"
##  [77] "CTAACGGAACCGAT" "CTACGCACTGGTCA" "CTACTCCTATGTCG" "CTAGTTACGAAACA"
##  [81] "CTATACTGCTACGA" "CTATACTGTCTCAT" "CTCGACTGGTTGAC" "CTGAGAACGTAAAG"
##  [85] "CTTTAGTGACGGGA" "GAACCAACTTCCGC" "GAAGTGCTAAACGA" "GAATGCACCTTCGC"
##  [89] "GAATTAACGTCGTA" "GACGGCACACGGGA" "GAGCGCTGAAGATG" "GAGGTACTGACACT"
##  [93] "GAGGTGGATCCTCG" "GATAGAGAAGGGTG" "GATCCCTGACCTTT" "GCACACCTGTGCTA"
##  [97] "GCACCACTTCCTTA" "GCACTAGAGTCGTA" "GCAGGGCTATCGAC" "GCCGGAACGTTCTT"
## [101] "GCCTACACAGTTCG" "GCGCATCTTGCTCC" "GCGCGATGGTGCAT" "GGAAGGTGGCGAGA"
## [105] "GGACGCTGTCCTCG" "GGAGGCCTCGTTGA" "GGCAAGGAAAAAGC" "GGCATATGCTTATC"
## [109] "GGCCGAACTCTAGG" "GGCTAAACACCTGA" "GGGTTAACGTGCAT" "GGTGGAGAAACGGG"
## [113] "GTAGTGTGAGCGGA" "GTCGACCTGAATGA" "GTGATTCTGGTTCA" "GTGTATCTAGTAGA"
## [117] "GTTAAAACCGAGAG" "GTTCAACTGGGACA" "GTTGACGATATCGG" "TAACTCACTCTACT"
## [121] "TAAGAGGACTTGTT" "TAATGCCTCGTCTC" "TACGGCCTGGGACA" "TACTACTGATGTCG"
## [125] "TACTCTGAATCGAC" "TACTGTTGAGGCGA" "TAGCATCTCAGCTA" "TAGCCCACAGCTAC"
## [129] "TAGGGACTGAACTC" "TAGTGGTGAAGTGA" "TAGTTAGAACCACA" "TATGAATGGAGGAC"
## [133] "TATGGGTGCATCAG" "TATTTCCTGGAGGT" "TCAACACTGTTTGG" "TCAGACGACGTTAG"
## [137] "TCCCGAACACAGTC" "TCCTAAACCGCATA" "TCGATTTGCAGCTA" "TCTAACACCAGTTG"
## [141] "TGATAAACTCCGTC" "TGCACAGACGACAT" "TGCCACTGCGATAC" "TGCTGAGAGAGCAG"
## [145] "TGGAACACAAACAG" "TGGTAGACCCTCAC" "TGTAATGACACAAC" "TGTAATGAGGTAAA"
## [149] "TTACTCGATCTACT" "TTAGTCTGCCAACA" "TTCCAAACTCCCAC" "TTCCCACTTGAGGG"
## [153] "TTCTAGTGGAGAGC" "TTCTGATGGAGACG" "TTGTCATGGACGGA"
# How can I extract expression matrix for all NK cells (perhaps, to load into another package)
nk.raw.data <- as.matrix(GetAssayData(pbmc, slot = "counts")[, WhichCells(pbmc, ident = "NK")])

# Can I create a Seurat object based on expression of a feature or value in object metadata?
subset(pbmc, subset = MS4A1 > 1)
## An object of class Seurat 
## 13714 features across 414 samples within 1 assay 
## Active assay: RNA (13714 features)
##  2 dimensional reductions calculated: pca, umap
subset(pbmc, subset = replicate == "rep2")
## An object of class Seurat 
## 13714 features across 1311 samples within 1 assay 
## Active assay: RNA (13714 features)
##  2 dimensional reductions calculated: pca, umap
# Can I create a Seurat object of just the NK cells and B cells?
subset(pbmc, idents = c("NK", "B"))
## An object of class Seurat 
## 13714 features across 499 samples within 1 assay 
## Active assay: RNA (13714 features)
##  2 dimensional reductions calculated: pca, umap
# Can I create a Seurat object of all cells except the NK cells and B cells?
subset(pbmc, idents = c("NK", "B"), invert = TRUE)
## An object of class Seurat 
## 13714 features across 2139 samples within 1 assay 
## Active assay: RNA (13714 features)
##  2 dimensional reductions calculated: pca, umap
# note that if you wish to perform additional rounds of clustering after subsetting we recommend
# re-running FindVariableFeatures() and ScaleData()

Calculating the average gene expression within a cluster

# How can I calculate the average expression of all cells within a cluster?
cluster.averages <- AverageExpression(pbmc)
head(cluster.averages[["RNA"]][, 1:5])
Naive CD4 T Memory CD4 T CD14+ Mono B CD8 T
AL627309.1 0.0061287 0.0059273 0.0485434 0.0000000 0.0205459
AP006222.2 0.0000000 0.0082061 0.0108847 0.0000000 0.0119149
RP11-206L10.2 0.0074531 0.0000000 0.0000000 0.0206503 0.0000000
RP11-206L10.9 0.0000000 0.0000000 0.0105012 0.0000000 0.0000000
LINC00115 0.0191189 0.0246905 0.0375374 0.0388854 0.0194828
NOC2L 0.4974632 0.3598115 0.2725375 0.5865349 0.5570490
# Return this information as a Seurat object (enables downstream plotting and analysis) First,
# replace spaces with underscores '_' so ggplot2 doesn't fail
orig.levels <- levels(pbmc)
Idents(pbmc) <- gsub(pattern = " ", replacement = "_", x = Idents(pbmc))
orig.levels <- gsub(pattern = " ", replacement = "_", x = orig.levels)
levels(pbmc) <- orig.levels
cluster.averages <- AverageExpression(pbmc, return.seurat = TRUE)
cluster.averages
## An object of class Seurat 
## 13714 features across 9 samples within 1 assay 
## Active assay: RNA (13714 features)
# How can I plot the average expression of NK cells vs. CD8 T cells?  Pass do.hover = T for an
# interactive plot to identify gene outliers
CellScatter(cluster.averages, cell1 = "NK", cell2 = "CD8_T")

# How can I calculate expression averages separately for each replicate?
cluster.averages <- AverageExpression(pbmc, return.seurat = TRUE, add.ident = "replicate")
CellScatter(cluster.averages, cell1 = "CD8_T_rep1", cell2 = "CD8_T_rep2")

# You can also plot heatmaps of these 'in silico' bulk datasets to visualize agreement between
# replicates
DoHeatmap(cluster.averages, features = unlist(TopFeatures(pbmc[["pca"]], balanced = TRUE)), size = 3, 
    draw.lines = FALSE)