1 Visualizing Seurat objects
library(Seurat)
dir = "data/filtered_gene_bc_matrices/hg19"
pbmc.data <- Read10X(data.dir = dir)
pbmc <- CreateSeuratObject(counts = pbmc.data, project = "pbmc3k",
min.cells=3, min.features=200)
pbmc## An object of class Seurat
## 13714 features across 2700 samples within 1 assay
## Active assay: RNA (13714 features, 0 variable features)
## 1 layer present: countspbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, pattern = "^MT-")
pbmc <- subset(pbmc,
subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mt < 5
)pbmc <- NormalizeData(pbmc, normalization.method = "LogNormalize",
scale.factor = 10000)
pbmc <- ScaleData(pbmc)pbmc <- FindVariableFeatures(pbmc, selection.method = "vst",
nfeatures = 2000)
pbmc <- RunPCA(pbmc, features = VariableFeatures(object = pbmc))
pbmc <- RunUMAP(pbmc, dims = 1:10)## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
##
## Number of nodes: 2638
## Number of edges: 95927
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.8728
## Number of communities: 9
## Elapsed time: 0 seconds## Assigning cell type identity to clusters
cluster.ids <- c("Naive CD4 T", "CD14+ Mono", "Memory CD4 T",
"B", "CD8 T", "FCGR3A+ Mono", "NK", "DC", "Platelet")
names(cluster.ids) <- levels(pbmc)
pbmc <- RenameIdents(pbmc, cluster.ids)1.1 Dimensional reduction plot
# DimPlot(pbmc, reduction = "umap",
# label = TRUE, pt.size = 0.5)
library(ggplot2)
library(ggsc)
sc_dim(pbmc) + sc_dim_geom_label()
# using vector point layer when geom = geom_bgpoint
sc_dim(pbmc, geom = geom_bgpoint, pointsize=1) + sc_dim_geom_label()
p <- sc_dim(pbmc) +
sc_dim_geom_label(geom = shadowtext::geom_shadowtext,
color='black', bg.color='white')
p
f <- sc_dim(pbmc, bg_colour = 'black') +
sc_dim_geom_label(geom = shadowtext::geom_shadowtext,
color = 'black', bg.color='white')
f
f2 <- sc_dim(pbmc,
mapping = aes(bg_colour = ident),
geom = geom_bgpoint,
pointsize=1,
gap_line_width=.2
) +
sc_dim_geom_label(geom = shadowtext::geom_shadowtext,
color = 'black', bg.color='white')
f2
The bg_colour aesthetic has its own scales: scale_bg_colour_discrete() is the default, scale_bg_colour_manual() sets the colours by hand, and scale_bg_color_manual() is the same function under the American spelling.
By default sc_dim() and sc_feature() draw the cells with sc_geom_point(), a rasterised point layer that stays fast with tens of thousands of cells. The geom_bgpoint() layer used above is the vector version that adds the background colour, and geom_scattermore2() is the rasterised version of that same layer, for a bg_colour plot with more cells than the vector layer can draw comfortably:
sc_dim(pbmc, geom = geom_scattermore2, pointsize = 1,
mapping = aes(bg_colour = ident), gap_line_width = .2) +
sc_dim_geom_label(geom = shadowtext::geom_shadowtext,
color = 'black', bg.color = 'white')
The number of cells in each clusters can be easily Visualized using sc_dim_count(). The colors of the bar plot is consistent with the dimensional reduction plot.
library(dplyr)
top20 <- pbmc.markers |>
dplyr::group_by(cluster) |>
dplyr::filter(avg_log2FC > 1) |>
dplyr::slice_head(n = 20) |>
dplyr::ungroup()
library(clusterProfiler)
library(enrichplot)
gg <- bitr(top20$gene, 'SYMBOL', 'ENTREZID', 'org.Hs.eg.db')
top20 <- merge(top20, gg, by.x='gene', by.y = 'SYMBOL')
kk <- compareCluster(ENTREZID~cluster, data = top20, fun=enrichKEGG)g <- dotplot(kk, label_format=100) +
aes(x=sub("\n.*", "", Cluster)) +
xlab("Cell Clusters") +
ggtitle(NULL) +
theme(axis.text.x = element_text(angle=30, hjust=1))
p3 <- p2 + coord_cartesian() +
ggfun::theme_noxaxis() +
xlab(NULL)
insert_top(g, p3, height=.2)
1.2 Visualize ‘features’ on a dimensional reduction plot
features = c("MS4A1", "GNLY", "CD3E",
"CD14", "FCER1A", "FCGR3A",
"LYZ", "PPBP", "CD8A")
# FeaturePlot(pbmc,'CD4')
sc_feature(pbmc, 'CD4')
Here is the real ‘features’ on dimensional plot
sc_dim(pbmc, alpha=.3) +
ggnewscale::new_scale_color() +
sc_dim_geom_feature(pbmc, features, mapping=aes(color=features)) +
scale_color_viridis_d()
sc_dim(pbmc, alpha=.3) +
ggnewscale::new_scale_color() +
sc_dim_geom_feature(pbmc, features, mapping=aes(color=features)) +
scale_color_viridis_d()
sc_feature(pbmc,
mapping = aes(bg_colour = RNA_snn_res.0.5),
features = features,
reduction = 'umap',
geom = geom_bgpoint,
pointsize = 1
) +
ggnewscale::new_scale_color() +
sc_dim_geom_ellipse(
mapping = aes(group = RNA_snn_res.0.5, colour=RNA_snn_res.0.5),
) +
sc_dim_geom_label(
geom = shadowtext::geom_shadowtext,
mapping = aes(label = RNA_snn_res.0.5, colour = RNA_snn_res.0.5),
bg.color = 'white',
size = 5
)
p1 <- sc_feature(pbmc, features='CD3E', density = FALSE) + scale_color_viridis_c()
p2 <- sc_feature(pbmc, features='CD3E', density = TRUE) + scale_color_viridis_c()
p1 | p2
f1 <- sc_feature(pbmc, features = 'CD3E',
mapping = aes(bg_colour = ident),
geom = geom_bgpoint,
density = FALSE,
bg_line_width = .4,
gap_line_width = .2,
size = .5) +
scale_color_viridis_c()
f2 <- sc_feature(pbmc, features = 'CD3E',
mapping = aes(bg_colour = ident),
geom = geom_bgpoint,
density = TRUE,
bg_line_width = .4,
gap_line_width = .22,
size = .8) +
scale_color_viridis_c()
f1 | f2
1.3 Visualize ellipse on a dimensional reduction plot
sc_dim_geom_ellipse() adds one ellipse per cluster to a dimensional reduction plot. The default geom is stat_ellipse(), and the level argument sets its confidence level:
The geom argument also takes any other layer function that draws a group. With ggforce::geom_mark_hull(), the cells of each cluster are outlined with a hull instead of an ellipse:
sel <- c("B", "NK", "FCGR3A+ Mono", "DC")
sc_dim(subset(pbmc, idents = sel)) +
sc_dim_geom_ellipse(geom = ggforce::geom_mark_hull, concavity = 2) +
sc_dim_geom_label(geom = shadowtext::geom_shadowtext,
color = 'black', bg.color = 'white', size = 5)
level is specific to stat_ellipse() and is only forwarded to a layer that accepts it, so geom_mark_hull() is not passed an unknown parameter. Arguments of the chosen layer, such as concavity, are supplied through ....
A hull follows the cells rather than a covariance ellipse, so it is best used for clusters that are compact and well separated; a single cell far from the rest of its cluster is drawn as a long thin spike. The clusters drawn here are a subset for that reason. sc_dim_geom_label() combines with either layer and places each label at a robust centre estimated from the coordinates of its cluster.
1.4 Visualize selected clusters
sc_dim(pbmc, color='grey', geom = geom_bgpoint, pointsize=1) +
sc_dim_geom_sub(subset=selected, mapping = aes(bg_colour = ident)) +
sc_dim_geom_label(geom = shadowtext::geom_shadowtext,
mapping = aes(subset = ident %in% selected),
color='black', bg.color='white')
1.5 Dot plot for selected features
features = c("MS4A1", "GNLY", "CD3E",
"CD14", "FCER1A", "FCGR3A",
"LYZ", "PPBP", "CD8A")
# DotPlot(pbmc, features = features,
# group.by = 'seurat_clusters')
sc_dot(pbmc, features=features)
1.6 Violin plot of gene expression
## allows applying an user-defined function to transform/filter the data
sc_violin(pbmc, 'CD4', .fun=function(d) dplyr::filter(d, value > 0)) +
ggforce::geom_sina(size=.1) +
guides(x = guide_axis(angle = -45))
#VlnPlot(pbmc, features)
sc_violin(pbmc, features) +
theme(axis.text.x = element_text(angle=45, hjust=1))
1.7 Spatial features
library(SeuratData)
if (!requireNamespace("SVP", quietly=TRUE)){
remotes::install_github("YuLab-SMU/SVP")
}
# InstallData("stxBrain")
brain <- LoadData("stxBrain", type = "anterior1")
# Normalization
brain <- SCTransform(brain, assay = "Spatial", verbose = FALSE)## SpatialFeaturePlot(brain, features = c("Hpca", "Ttr"))
f1 <- sc_spatial(brain, features = c("Hpca", "Ttr"), image.mirror.axis = 'v', geom = geom_bgpoint, pointsize=1.5)
f2 <- sc_spatial(brain, features=c("Hpca", "Ttr"), image.mirror.axis='v', density=TRUE, geom=geom_bgpoint, pointsize=1.5)
plot_list(f1, f2, ncol=1)
To mark the spatially variable genes on such a plot, SVP::runLISA() computes a local indicator of spatial association for each spot. Its result is drawn by plot_lisa_feature(), which takes the SpatialExperiment and the runLISA() result and returns the spatial plot with the significant spots outlined. The same annotation can be built layer by layer with sc_geom_annot(), which is what a density plot needs.
{library(SingleCellExperiment)
library(SpatialExperiment)} |> suppressPackageStartupMessages()
spe.brain <- brain |> as.SingleCellExperiment()
spe.brain <- as(spe.brain, "SpatialExperiment")
spatialCoords(spe.brain) <- GetTissueCoordinates(brain)[,seq(2)] |> as.matrix()
spe.brain## class: SpatialExperiment
## dim: 17668 2696
## metadata(0):
## assays(2): counts logcounts
## rownames(17668): Xkr4 Sox17 ... CAAA01118383.1
## CAAA01147332.1
## rowData names(0):
## colnames(2696): AAACAAGTATCTCCCA-1 AAACACCAATAACTGC-1
## ... TTGTTTCACATCCAGG-1 TTGTTTCCATACAACT-1
## colData names(9): orig.ident nCount_Spatial ... ident
## sample_id
## reducedDimNames(0):
## mainExpName: SCT
## altExpNames(1): Spatial
## spatialCoords names(2) : x y
## imgData names(0):library(SVP)
res.lisa1 <- runLISA(spe.brain, features = c("Hpca", "Ttr"), assay.type = 'logcounts', weight.method='knn', k=20)
da.lisa1 <- res.lisa1 |> lapply(function(x)x |> tibble::rownames_to_column(var='.BarcodeID')) |> dplyr::bind_rows(.id='features')
head(da.lisa1)## features .BarcodeID Gi E.Gi
## 1 Hpca AAACAAGTATCTCCCA-1 0.0002755891 0.0003710575
## 2 Hpca AAACACCAATAACTGC-1 0.0001295426 0.0003710575
## 3 Hpca AAACAGAGCGACTCCT-1 0.0004076779 0.0003710575
## 4 Hpca AAACAGCTTTCAGAAG-1 0.0001373876 0.0003710575
## 5 Hpca AAACAGGGTCTATATT-1 0.0001281104 0.0003710575
## 6 Hpca AAACATGGTGAGAGGA-1 0.0001670363 0.0003710575
## Var.Gi Z.Gi Pr (z != E(Gi)) padj
## 1 1.257273e-09 -2.692433 7.093284e-03 1.537258e-02
## 2 1.254907e-09 -6.817704 9.250707e-12 1.458474e-10
## 3 1.258419e-09 1.032309 3.019273e-01 4.011808e-01
## 4 1.256547e-09 -6.591942 4.341096e-11 4.938225e-10
## 5 1.254907e-09 -6.858134 6.976566e-12 1.210529e-10
## 6 1.254907e-09 -5.759298 8.446424e-09 4.773912e-08
## cluster.no.test cluster.test
## 1 Low Low
## 2 Low Low
## 3 High NoSign
## 4 Low Low
## 5 Low Low
## 6 Low Low# plot_lisa_feature() draws the annotation in one call
plot_lisa_feature(spe.brain, lisa.res = res.lisa1, features = c("Hpca", "Ttr"),
geom = geom_bgpoint, pointsize = 1.5, hlpointsize = 1.5,
gap_line_width = .3, bg_line_width = .3)
# the same annotation is added layer by layer with sc_geom_annot(),
# which is what a density plot needs
f2 + sc_geom_annot(
data = da.lisa1,
mapping = aes(bg_colour = cluster.test, subset = cluster.test == 'High'),
pointsize = 1.5,
bg_line_width = .3,
gap_line_width = .3) +
scale_bg_colour_manual(values = 'black')
f3 <- sc_spatial(brain, features=c("Hpca", "Ttr"), image.mirror.axis='v',
density=TRUE, joint=TRUE)
f4 <- sc_spatial(brain, features=c("Hpca", "Ttr"), image.mirror.axis='v',
density=TRUE, joint=TRUE, common.legend=FALSE)
plot_list(f3, f4, ncol=1)
genes <- rownames(brain) |> sample(10)
sc_spatial(brain, features=genes, plot.pie=TRUE,
pie.radius.scale=.35, image.plot=FALSE,
color = NA
)