Seurat - Interaction Tips

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 FCGR3A+ Mono 
##          697          483          480          344          271          162 
##           NK           DC     Platelet 
##          155           32           14

# How many cells are in each replicate?
table(pbmc$replicate)

## 
## rep1 rep2 
## 1348 1290

# What proportion of cells are in each cluster?
prop.table(table(Idents(pbmc)))

## 
##  Naive CD4 T Memory CD4 T   CD14+ Mono            B        CD8 T FCGR3A+ Mono 
##  0.264215315  0.183093252  0.181956027  0.130401820  0.102729340  0.061410159 
##           NK           DC     Platelet 
##  0.058756634  0.012130402  0.005307051

# How does cluster membership vary by replicate?
table(Idents(pbmc), pbmc$replicate)

##               
##                rep1 rep2
##   Naive CD4 T   354  343
##   Memory CD4 T  249  234
##   CD14+ Mono    232  248
##   B             173  171
##   CD8 T         154  117
##   FCGR3A+ Mono   81   81
##   NK             81   74
##   DC             18   14
##   Platelet        6    8

prop.table(table(Idents(pbmc), pbmc$replicate), margin = 2)

##               
##                       rep1        rep2
##   Naive CD4 T  0.262611276 0.265891473
##   Memory CD4 T 0.184718101 0.181395349
##   CD14+ Mono   0.172106825 0.192248062
##   B            0.128338279 0.132558140
##   CD8 T        0.114243323 0.090697674
##   FCGR3A+ Mono 0.060089021 0.062790698
##   NK           0.060089021 0.057364341
##   DC           0.013353116 0.010852713
##   Platelet     0.004451039 0.006201550

Selecting particular cells and subsetting the Seurat object

# What are the cell names of all NK cells?
WhichCells(pbmc, idents = "NK")

##   [1] "AAACCGTGTATGCG-1" "AAATTCGATTCTCA-1" "AACCTTACGCGAGA-1"
##   [4] "AACGCCCTCGTACA-1" "AACGTCGAGTATCG-1" "AAGATTACCTCAAG-1"
##   [7] "AAGCAAGAGCTTAG-1" "AAGCAAGAGGTGTT-1" "AAGTAGGATACAGC-1"
##  [10] "AATACTGAATTGGC-1" "AATCCTTGGTGAGG-1" "AATCTCTGCTTTAC-1"
##  [13] "ACAAATTGTTGCGA-1" "ACAACCGAGGGATG-1" "ACAATTGATGACTG-1"
##  [16] "ACACCCTGGTGTTG-1" "ACAGGTACTGGTGT-1" "ACCTGGCTAAGTAG-1"
##  [19] "ACGAACACCTTGTT-1" "ACGATCGAGGACTT-1" "ACGCAATGGTTCAG-1"
##  [22] "ACGCTGCTGTTCTT-1" "ACGGAACTCAGATC-1" "ACGTGATGTGACAC-1"
##  [25] "ACGTTGGAGCCAAT-1" "ACTGCCACTCCGTC-1" "ACTGGCCTTCAGTG-1"
##  [28] "ACTTCAACGTAGGG-1" "AGAACAGAAATGCC-1" "AGATATACCCGTAA-1"
##  [31] "AGATTCCTGTTCAG-1" "AGCCTCTGCCAATG-1" "AGCGATTGAGATCC-1"
##  [34] "AGGATGCTTTAGGC-1" "AGGGACGAGTCAAC-1" "AGTAATACATCACG-1"
##  [37] "AGTCACGATGAGCT-1" "AGTTTGCTACTGGT-1" "ATACCACTGCCAAT-1"
##  [40] "ATACTCTGGTATGC-1" "ATCCCGTGCAGTCA-1" "ATCTTTCTTGTCCC-1"
##  [43] "ATGAAGGACTTGCC-1" "ATGATAACTTCACT-1" "ATGATATGGTGCTA-1"
##  [46] "ATGGACACGCATCA-1" "ATGGGTACATCGGT-1" "ATTAACGATGAGAA-1"
##  [49] "ATTCCAACTTAGGC-1" "CAAGGTTGTCTGGA-1" "CAATCTACTGACTG-1"
##  [52] "CACCACTGGCGAAG-1" "CACGGGTGGAGGAC-1" "CAGATGACATTCTC-1"
##  [55] "CAGCAATGGAGGGT-1" "CAGCGGACCTTTAC-1" "CAGCTCTGTGTGGT-1"
##  [58] "CAGTTTACACACGT-1" "CATCAGGACTTCCG-1" "CATCAGGATAGCCA-1"
##  [61] "CATGAGACGTTGAC-1" "CATTACACCAACTG-1" "CATTTCGAGATACC-1"
##  [64] "CCTCGAACACTTTC-1" "CGACCACTAAAGTG-1" "CGACCACTGCCAAT-1"
##  [67] "CGAGGCTGACGCTA-1" "CGCCGAGAGCTTAG-1" "CGGCGAACGACAAA-1"
##  [70] "CGGCGAACTACTTC-1" "CGGGCATGTCTCTA-1" "CGTACCTGGCATCA-1"
##  [73] "CGTGTAGACGATAC-1" "CGTGTAGAGTTACG-1" "CGTGTAGATTCGGA-1"
##  [76] "CTAAACCTCTGACA-1" "CTAACGGAACCGAT-1" "CTACGCACTGGTCA-1"
##  [79] "CTACTCCTATGTCG-1" "CTAGTTACGAAACA-1" "CTATACTGCTACGA-1"
##  [82] "CTATACTGTCTCAT-1" "CTCGACTGGTTGAC-1" "CTGAGAACGTAAAG-1"
##  [85] "CTTTAGTGACGGGA-1" "GAACCAACTTCCGC-1" "GAAGTGCTAAACGA-1"
##  [88] "GAATGCACCTTCGC-1" "GAATTAACGTCGTA-1" "GACGGCACACGGGA-1"
##  [91] "GAGCGCTGAAGATG-1" "GAGGTACTGACACT-1" "GAGGTGGATCCTCG-1"
##  [94] "GATAGAGAAGGGTG-1" "GATCCCTGACCTTT-1" "GCACACCTGTGCTA-1"
##  [97] "GCACCACTTCCTTA-1" "GCACTAGAGTCGTA-1" "GCAGGGCTATCGAC-1"
## [100] "GCCGGAACGTTCTT-1" "GCCTACACAGTTCG-1" "GCGCATCTTGCTCC-1"
## [103] "GCGCGATGGTGCAT-1" "GGAAGGTGGCGAGA-1" "GGACGCTGTCCTCG-1"
## [106] "GGAGGCCTCGTTGA-1" "GGCAAGGAAAAAGC-1" "GGCATATGCTTATC-1"
## [109] "GGCCGAACTCTAGG-1" "GGCTAAACACCTGA-1" "GGGTTAACGTGCAT-1"
## [112] "GGTGGAGAAACGGG-1" "GTAGTGTGAGCGGA-1" "GTCGACCTGAATGA-1"
## [115] "GTGATTCTGGTTCA-1" "GTGTATCTAGTAGA-1" "GTTAAAACCGAGAG-1"
## [118] "GTTCAACTGGGACA-1" "GTTGACGATATCGG-1" "TAACTCACTCTACT-1"
## [121] "TAAGAGGACTTGTT-1" "TAATGCCTCGTCTC-1" "TACGGCCTGGGACA-1"
## [124] "TACTACTGATGTCG-1" "TACTCTGAATCGAC-1" "TACTGTTGAGGCGA-1"
## [127] "TAGCATCTCAGCTA-1" "TAGCCCACAGCTAC-1" "TAGGGACTGAACTC-1"
## [130] "TAGTGGTGAAGTGA-1" "TAGTTAGAACCACA-1" "TATGAATGGAGGAC-1"
## [133] "TATGGGTGCATCAG-1" "TATTTCCTGGAGGT-1" "TCAACACTGTTTGG-1"
## [136] "TCAGACGACGTTAG-1" "TCCCGAACACAGTC-1" "TCCTAAACCGCATA-1"
## [139] "TCGATTTGCAGCTA-1" "TCTAACACCAGTTG-1" "TGATAAACTCCGTC-1"
## [142] "TGCACAGACGACAT-1" "TGCCACTGCGATAC-1" "TGCTGAGAGAGCAG-1"
## [145] "TGGAACACAAACAG-1" "TGGTAGACCCTCAC-1" "TGTAATGACACAAC-1"
## [148] "TGTAATGAGGTAAA-1" "TTACTCGATCTACT-1" "TTAGTCTGCCAACA-1"
## [151] "TTCCAAACTCCCAC-1" "TTCCCACTTGAGGG-1" "TTCTAGTGGAGAGC-1"
## [154] "TTCTGATGGAGACG-1" "TTGTCATGGACGGA-1"

# 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, 2000 variable features)
##  2 dimensional reductions calculated: pca, umap

subset(pbmc, subset = replicate == "rep2")

## An object of class Seurat 
## 13714 features across 1290 samples within 1 assay 
## Active assay: RNA (13714 features, 2000 variable 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, 2000 variable 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, 2000 variable 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, 0 variable 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")