In Ghosn-Lab/BatchNorm: Data Matrix Normalization and Merging Strategies Minimize Batch-specific Systemic Variation in scRNA-Seq Data
library(knitr)
library(dplyr)
library(Seurat)
library(ggplot2)
opts_chunk$set(
fig.pos = "!h", out.extra = "", warning = F,
fig.align = "center", fig.height = 7, fig.width = 7
)
Processing the Seurat Object:
Data filtering & normalization
library(BatchNorm)
# Import unfiltered Seurat PBMC4 object & extract "4A" (included with 'BatchNorm' package)
data(package = 'BatchNorm', PBMC4)
Idents(PBMC4) <- PBMC4[["orig.ident"]]
PBMC4A <- subset(PBMC4, idents = "Sample_4A")
rm(PBMC4)
# Run "standard" Seurat workflow
PBMC4A <- PBMC4A %>%
MitoFilter() %>%
NormalizeData(assay = "RNA", verbose = FALSE) %>%
FindVariableFeatures(selection.method = "vst", nfeatures = 2000, verbose = FALSE) %>%
ScaleData(verbose = FALSE) %>%
RunPCA(npcs = 30, verbose = FALSE)
Selecting the appropriate number of Principal Components for UMAP reduction
## Identify correct numbers of PCs
## (Takes up to 5 minutes. Not run while rendering vignette for time)
# PBMC4A.pca.test <- TestPCA(object = PBMC4A)
# PBMC4A.pca.test[, 1:20]
## 16 PCs with a z-score > 1
## Proceed with 16 PCs for dimensional reduction & clustering
## Visualize PCs plotted by standard deviation:
ElbowPlot(PBMC4A)
Selecting the appropriate number of Principal Components for UMAP reduction
PBMC4A <- PBMC4A %>%
RunUMAP(reduction = "pca", dims = 1:16, verbose = FALSE) %>%
FindNeighbors(reduction = "pca", dims = 1:16, verbose = FALSE) %>%
FindClusters(resolution = .8, verbose = FALSE)
# Visualize clusters on UMAP
DimPlot(PBMC4A, reduction = "umap", cols = colors.use) +
ggtitle("PBMC Clusters")
Biaxial Antigen Gating
B Cell Classification
# Produce scatter plot of CD20 vs CD19
d <- as.data.frame(FetchData(object = PBMC4A,
vars = c("adt_CD20", "adt_CD19", "seurat_clusters")))
## Format colnames to match FeatureScatter syntax
colnames(d) <- c("adt_CD20", "adt_CD19", "colors")
FeatureScatter(PBMC4A, "adt_CD20", "adt_CD19",
cols = colors.use, plot.cor = F) +
geom_density2d(data = d, color = "black", size = 1) +
theme_gray() +
labs(x = "CD20 (ADT)",
y = "CD19 (ADT)")
# clusters 6 & 9 contain B cells
# Subset B cells (CD20 vs CD19) and isolate Naive B cells (IGD+/CD27-)
bcells <- subset(PBMC4A, `adt_CD20` > 25 | `adt_CD19` > 10)
bcells[["IGD"]] <- colSums(bcells@assays$RNA@data[grep("^IGHD", rownames(bcells)), ])
d <- as.data.frame(FetchData(object = bcells,
vars = c("IGD", "adt_CD27", "seurat_clusters")))
## Format colnames to match FeatureScatter syntax
colnames(d) <- c("IGD", "adt_CD27", "colors")
FeatureScatter(bcells, "adt_CD27", "IGD", plot.cor = F,
cols = colors.use[c(2, 4, 7, 9, 10, 14, 15)]) +
geom_density2d(data = d, color = "black", size = 1) +
labs(x = "CD27 (ADT)",
y = "Total IGHD Expression (RNA)") + theme_gray() +
coord_cartesian(ylim = c(-0.3, 3.5), expand = T)
# cluster 6 contains Naive B cells
NK Cell Classification
# Remove B cells (already classified) and continue to classify non-B cell PBMCs
nonclassified <- subset(PBMC4A, idents = c(6, 9), invert = T)
d <- as.data.frame(FetchData(object = nonclassified,
vars = c("adt_CD56", "adt_CD14", "seurat_clusters")))
## Format colnames to match FeatureScatter syntax
colnames(d) <- c("adt_CD56", "adt_CD14", "colors")
FeatureScatter(nonclassified, "adt_CD14", "adt_CD56", plot.cor = F,
cols = colors.use[c(1:6, 8:9, 11:17)], pt.size = 3) +
geom_density2d(data = d, color = "black", size = 1) +
labs(x = "CD14 (ADT)",
y = "CD56 (ADT)") + theme_gray()
# clusters 3, 4, 5 & 12 appear to contain NK cells
# Cluster 12 appears to contain mostly CD56Hi NK cells
# We can further classify NK cells and confirm using CD3 ADT & NKG7 transcript
nkcells <- subset(PBMC4A, idents = c(3, 4, 5, 12))
d <- as.data.frame(FetchData(object = nkcells,
vars = c("NKG7", "adt_CD3", "seurat_clusters")))
## Format colnames to match FeatureScatter syntax
colnames(d) <- c("NKG7", "adt_CD3", "colors")
FeatureScatter(nkcells, "adt_CD3", "NKG7",
cols = colors.use[c(4, 5, 6, 14)], pt.size = 3, plot.cor = F) +
geom_density2d(data = d, color = "black", size = 1) +
labs(x = "CD3 (ADT)",
y = "NKG7 (RNA)") + theme_gray()
# cluster 12 = CD56Hi NK cells, Clusters 3 & 4 = NK Cells, Cluster 5 = NK-T Cells
Monocyte Classification
nonclassified <- subset(PBMC4A, idents = c(6, 9, 3, 4, 12, 5), invert = T)
d <- as.data.frame(FetchData(object = nonclassified,
vars = c("adt_CD16", "adt_CD14", "seurat_clusters")))
colnames(d) <- c("adt_CD16", "adt_CD14", "colors")
FeatureScatter(nonclassified, "adt_CD14", "adt_CD16",
cols = colors.use[c(1:3, 8:9, 11:12, 14:17)], pt.size = 3, plot.cor = F) +
geom_density2d(data = d, color = "black", size = 1) +
labs(x = "CD14 (ADT)",
y = "CD16 (ADT)") + theme_gray()
# Cluster 8 = classical monocytes, 13 = non-classical monocytes
HSPC Classification
nonclassified <- subset(PBMC4A, idents = c(6, 9, 3, 4, 12, 5, 8, 13), invert = T)
# Create a variable comparing only cluster 16 to all other unclassified clusters
nonclassified[["HSPC_gate"]] <- rep("Other Clusters", nrow(nonclassified[[]]))
HSPCs <- WhichCells(nonclassified, idents = 16)
nonclassified@meta.data[HSPCs, "HSPC_gate"] <- "16"
nonclassified@meta.data[, "HSPC_gate"] <- factor(nonclassified@meta.data[, "HSPC_gate"],
levels = c("Other Clusters", "16"))
VlnPlot(nonclassified, "CD34", cols = colors.use[c(1:3, 8, 11:12, 15:17)],
group.by = "HSPC_gate")
# Cluster 16 = HSPCs
T Cell Classification
nonclassified <- subset(PBMC4A, idents = c(6, 9, 3, 4, 12, 8, 5, 13, 16), invert = T)
d <- as.data.frame(FetchData(object = nonclassified,
vars = c("adt_CD8A", "adt_CD4", "seurat_clusters")))
colnames(d) <- c("adt_CD8A", "adt_CD4", "colors")
FeatureScatter(nonclassified, "adt_CD8A", "adt_CD4",
cols = colors.use[c(1:3, 8, 11:12, 15:16)], pt.size = 3, plot.cor = F) +
geom_density2d(data = d, color = "black", size = 1) +
labs(x = "CD8A (ADT)",
y = "CD4 (ADT)") + theme_gray()
# 0, 1, 2, 15 = CD4+ T-cell; 7, 11 = CD8+ T-cell
nonclassified <- subset(PBMC4A, idents = c(6, 9, 3, 4, 12, 8, 5, 13, 16, 7, 11), invert = T)
# Create a variable comparing only cluster 10 to all other unclassified clusters (other CD4+ T Cells)
nonclassified[["Treg_gate"]] <- rep("Other CD4+", nrow(nonclassified[[]]))
Tregs <- WhichCells(nonclassified, idents = 10)
nonclassified@meta.data[Tregs, "Treg_gate"] <- "10"
nonclassified@meta.data[, "Treg_gate"] <- factor(nonclassified@meta.data[, "Treg_gate"],
levels = c("Other CD4+", "10"))
VlnPlot(nonclassified, "FOXP3", cols = colors.use[c(11, 15)],
group.by = "Treg_gate")
# Cluster 10 = Treg cells
Plasmacytoid Dendritic Cell Classification
nonclassified <- subset(PBMC4A, idents = c(6, 9, 3, 4, 12, 8, 5, 13, 16, 7, 11), invert = T)
nonclassified[["pDC_gate"]] <- rep("Other CD4+", nrow(nonclassified[[]]))
PDCs <- WhichCells(nonclassified, idents = 14)
nonclassified@meta.data[PDCs, "pDC_gate"] <- "14"
nonclassified@meta.data[, "pDC_gate"] <- factor(nonclassified@meta.data[, "pDC_gate"],
levels = c("Other CD4+", "14"))
VlnPlot(nonclassified, 'LILRA4', cols = colors.use[c(11, 15)],
group.by = "pDC_gate")
# Cluster 14 = pDCs
Assigning Cell Classification Labels
# Assign cell classification labels to clusters, as follows:
# 6 = B Naive
# 9 = B memory
# 3, 4 = NK
# 5 = NK_T
# 12 = CD56Hi NK
# 8 = classical monocytes
# 13 = non-classical monocytes
# 16 = HSPCs
# 0, 1, 2, 15 = CD4+ T-cell
# 7, 11 = CD8+ T-cell
# 10 = TRegs
# 14 = DCs
Idents(PBMC4A) <- PBMC4A[["seurat_clusters"]]
Idents(PBMC4A) <- plyr::mapvalues(Idents(PBMC4A), from = c(0, 1, 2, 3, 4, 5,
6, 7, 8, 9, 10,
11, 12, 13, 14, 15,
16),
to = c('T_CD4', 'T_CD4', 'T_CD4', 'NK', 'NK', 'NK_T',
'B_Naive', 'T_CD8', 'Monocyte_Classical', 'B_Memory', 'TReg',
'T_CD8', 'NK_CD56Hi', 'Monocyte_NonClassical', 'Dendritic_Cells', 'T_CD4',
'HSPCs'))
# Order cell classifications
Idents(PBMC4A) <- factor(Idents(PBMC4A),
levels = c("B_Naive", "B_Memory", "T_CD4", "TReg",
"T_CD8", "NK_T", "NK", "NK_CD56Hi",
"Monocyte_Classical", "Monocyte_NonClassical",
"Dendritic_Cells", "HSPCs"))
PBMC4A[["Cell_Type"]] <- Idents(PBMC4A)
UMAPPlot(PBMC4A, cols = colors.use, pt.size = 2, label = F)
Ghosn-Lab/BatchNorm documentation built on Dec. 17, 2021, 9:32 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
library(knitr) library(dplyr) library(Seurat) library(ggplot2) opts_chunk$set( fig.pos = "!h", out.extra = "", warning = F, fig.align = "center", fig.height = 7, fig.width = 7 )
Processing the Seurat Object:
Data filtering & normalization
library(BatchNorm) # Import unfiltered Seurat PBMC4 object & extract "4A" (included with 'BatchNorm' package) data(package = 'BatchNorm', PBMC4) Idents(PBMC4) <- PBMC4[["orig.ident"]] PBMC4A <- subset(PBMC4, idents = "Sample_4A") rm(PBMC4) # Run "standard" Seurat workflow PBMC4A <- PBMC4A %>% MitoFilter() %>% NormalizeData(assay = "RNA", verbose = FALSE) %>% FindVariableFeatures(selection.method = "vst", nfeatures = 2000, verbose = FALSE) %>% ScaleData(verbose = FALSE) %>% RunPCA(npcs = 30, verbose = FALSE)
Selecting the appropriate number of Principal Components for UMAP reduction
## Identify correct numbers of PCs ## (Takes up to 5 minutes. Not run while rendering vignette for time) # PBMC4A.pca.test <- TestPCA(object = PBMC4A) # PBMC4A.pca.test[, 1:20] ## 16 PCs with a z-score > 1 ## Proceed with 16 PCs for dimensional reduction & clustering ## Visualize PCs plotted by standard deviation: ElbowPlot(PBMC4A)
Selecting the appropriate number of Principal Components for UMAP reduction
PBMC4A <- PBMC4A %>% RunUMAP(reduction = "pca", dims = 1:16, verbose = FALSE) %>% FindNeighbors(reduction = "pca", dims = 1:16, verbose = FALSE) %>% FindClusters(resolution = .8, verbose = FALSE) # Visualize clusters on UMAP DimPlot(PBMC4A, reduction = "umap", cols = colors.use) + ggtitle("PBMC Clusters")
Biaxial Antigen Gating
B Cell Classification
# Produce scatter plot of CD20 vs CD19 d <- as.data.frame(FetchData(object = PBMC4A, vars = c("adt_CD20", "adt_CD19", "seurat_clusters"))) ## Format colnames to match FeatureScatter syntax colnames(d) <- c("adt_CD20", "adt_CD19", "colors") FeatureScatter(PBMC4A, "adt_CD20", "adt_CD19", cols = colors.use, plot.cor = F) + geom_density2d(data = d, color = "black", size = 1) + theme_gray() + labs(x = "CD20 (ADT)", y = "CD19 (ADT)") # clusters 6 & 9 contain B cells # Subset B cells (CD20 vs CD19) and isolate Naive B cells (IGD+/CD27-) bcells <- subset(PBMC4A, `adt_CD20` > 25 | `adt_CD19` > 10) bcells[["IGD"]] <- colSums(bcells@assays$RNA@data[grep("^IGHD", rownames(bcells)), ]) d <- as.data.frame(FetchData(object = bcells, vars = c("IGD", "adt_CD27", "seurat_clusters"))) ## Format colnames to match FeatureScatter syntax colnames(d) <- c("IGD", "adt_CD27", "colors") FeatureScatter(bcells, "adt_CD27", "IGD", plot.cor = F, cols = colors.use[c(2, 4, 7, 9, 10, 14, 15)]) + geom_density2d(data = d, color = "black", size = 1) + labs(x = "CD27 (ADT)", y = "Total IGHD Expression (RNA)") + theme_gray() + coord_cartesian(ylim = c(-0.3, 3.5), expand = T) # cluster 6 contains Naive B cells
NK Cell Classification
# Remove B cells (already classified) and continue to classify non-B cell PBMCs nonclassified <- subset(PBMC4A, idents = c(6, 9), invert = T) d <- as.data.frame(FetchData(object = nonclassified, vars = c("adt_CD56", "adt_CD14", "seurat_clusters"))) ## Format colnames to match FeatureScatter syntax colnames(d) <- c("adt_CD56", "adt_CD14", "colors") FeatureScatter(nonclassified, "adt_CD14", "adt_CD56", plot.cor = F, cols = colors.use[c(1:6, 8:9, 11:17)], pt.size = 3) + geom_density2d(data = d, color = "black", size = 1) + labs(x = "CD14 (ADT)", y = "CD56 (ADT)") + theme_gray() # clusters 3, 4, 5 & 12 appear to contain NK cells # Cluster 12 appears to contain mostly CD56Hi NK cells # We can further classify NK cells and confirm using CD3 ADT & NKG7 transcript nkcells <- subset(PBMC4A, idents = c(3, 4, 5, 12)) d <- as.data.frame(FetchData(object = nkcells, vars = c("NKG7", "adt_CD3", "seurat_clusters"))) ## Format colnames to match FeatureScatter syntax colnames(d) <- c("NKG7", "adt_CD3", "colors") FeatureScatter(nkcells, "adt_CD3", "NKG7", cols = colors.use[c(4, 5, 6, 14)], pt.size = 3, plot.cor = F) + geom_density2d(data = d, color = "black", size = 1) + labs(x = "CD3 (ADT)", y = "NKG7 (RNA)") + theme_gray() # cluster 12 = CD56Hi NK cells, Clusters 3 & 4 = NK Cells, Cluster 5 = NK-T Cells
Monocyte Classification
nonclassified <- subset(PBMC4A, idents = c(6, 9, 3, 4, 12, 5), invert = T) d <- as.data.frame(FetchData(object = nonclassified, vars = c("adt_CD16", "adt_CD14", "seurat_clusters"))) colnames(d) <- c("adt_CD16", "adt_CD14", "colors") FeatureScatter(nonclassified, "adt_CD14", "adt_CD16", cols = colors.use[c(1:3, 8:9, 11:12, 14:17)], pt.size = 3, plot.cor = F) + geom_density2d(data = d, color = "black", size = 1) + labs(x = "CD14 (ADT)", y = "CD16 (ADT)") + theme_gray() # Cluster 8 = classical monocytes, 13 = non-classical monocytes
HSPC Classification
nonclassified <- subset(PBMC4A, idents = c(6, 9, 3, 4, 12, 5, 8, 13), invert = T) # Create a variable comparing only cluster 16 to all other unclassified clusters nonclassified[["HSPC_gate"]] <- rep("Other Clusters", nrow(nonclassified[[]])) HSPCs <- WhichCells(nonclassified, idents = 16) nonclassified@meta.data[HSPCs, "HSPC_gate"] <- "16" nonclassified@meta.data[, "HSPC_gate"] <- factor(nonclassified@meta.data[, "HSPC_gate"], levels = c("Other Clusters", "16")) VlnPlot(nonclassified, "CD34", cols = colors.use[c(1:3, 8, 11:12, 15:17)], group.by = "HSPC_gate") # Cluster 16 = HSPCs
T Cell Classification
nonclassified <- subset(PBMC4A, idents = c(6, 9, 3, 4, 12, 8, 5, 13, 16), invert = T) d <- as.data.frame(FetchData(object = nonclassified, vars = c("adt_CD8A", "adt_CD4", "seurat_clusters"))) colnames(d) <- c("adt_CD8A", "adt_CD4", "colors") FeatureScatter(nonclassified, "adt_CD8A", "adt_CD4", cols = colors.use[c(1:3, 8, 11:12, 15:16)], pt.size = 3, plot.cor = F) + geom_density2d(data = d, color = "black", size = 1) + labs(x = "CD8A (ADT)", y = "CD4 (ADT)") + theme_gray() # 0, 1, 2, 15 = CD4+ T-cell; 7, 11 = CD8+ T-cell nonclassified <- subset(PBMC4A, idents = c(6, 9, 3, 4, 12, 8, 5, 13, 16, 7, 11), invert = T) # Create a variable comparing only cluster 10 to all other unclassified clusters (other CD4+ T Cells) nonclassified[["Treg_gate"]] <- rep("Other CD4+", nrow(nonclassified[[]])) Tregs <- WhichCells(nonclassified, idents = 10) nonclassified@meta.data[Tregs, "Treg_gate"] <- "10" nonclassified@meta.data[, "Treg_gate"] <- factor(nonclassified@meta.data[, "Treg_gate"], levels = c("Other CD4+", "10")) VlnPlot(nonclassified, "FOXP3", cols = colors.use[c(11, 15)], group.by = "Treg_gate") # Cluster 10 = Treg cells
Plasmacytoid Dendritic Cell Classification
nonclassified <- subset(PBMC4A, idents = c(6, 9, 3, 4, 12, 8, 5, 13, 16, 7, 11), invert = T) nonclassified[["pDC_gate"]] <- rep("Other CD4+", nrow(nonclassified[[]])) PDCs <- WhichCells(nonclassified, idents = 14) nonclassified@meta.data[PDCs, "pDC_gate"] <- "14" nonclassified@meta.data[, "pDC_gate"] <- factor(nonclassified@meta.data[, "pDC_gate"], levels = c("Other CD4+", "14")) VlnPlot(nonclassified, 'LILRA4', cols = colors.use[c(11, 15)], group.by = "pDC_gate") # Cluster 14 = pDCs
Assigning Cell Classification Labels
# Assign cell classification labels to clusters, as follows: # 6 = B Naive # 9 = B memory # 3, 4 = NK # 5 = NK_T # 12 = CD56Hi NK # 8 = classical monocytes # 13 = non-classical monocytes # 16 = HSPCs # 0, 1, 2, 15 = CD4+ T-cell # 7, 11 = CD8+ T-cell # 10 = TRegs # 14 = DCs Idents(PBMC4A) <- PBMC4A[["seurat_clusters"]] Idents(PBMC4A) <- plyr::mapvalues(Idents(PBMC4A), from = c(0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16), to = c('T_CD4', 'T_CD4', 'T_CD4', 'NK', 'NK', 'NK_T', 'B_Naive', 'T_CD8', 'Monocyte_Classical', 'B_Memory', 'TReg', 'T_CD8', 'NK_CD56Hi', 'Monocyte_NonClassical', 'Dendritic_Cells', 'T_CD4', 'HSPCs')) # Order cell classifications Idents(PBMC4A) <- factor(Idents(PBMC4A), levels = c("B_Naive", "B_Memory", "T_CD4", "TReg", "T_CD8", "NK_T", "NK", "NK_CD56Hi", "Monocyte_Classical", "Monocyte_NonClassical", "Dendritic_Cells", "HSPCs")) PBMC4A[["Cell_Type"]] <- Idents(PBMC4A) UMAPPlot(PBMC4A, cols = colors.use, pt.size = 2, label = F)
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.