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)
# 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"
# 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
# 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()
# 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)