Last Updated: June 19, 2023
Part 2: QA/QC, filtering, and normalization
In this workshop, we are using the filtered feature barcode matrix. While this lowers the likelihood of encountering barcodes that are not cell-associated within our expression matrix, it is still good practice to perform quality assurance / control on the experiment.
Setup
First, we need to load the required libraries.
library(Seurat)
library(biomaRt)
library(ggplot2)
library(knitr)
library(kableExtra)
library(reshape2)
If you are continuing directly from part 1, the experiment.aggregate object is likely already in your workspace. In case you cleared your workspace at the end of the previous section, or are working on this project at a later date, after quitting and re-starting R, you can use the readRDS function to read your saved Seurat object from part 1.
experiment.aggregate <- readRDS("scRNA_workshop_1.rds")
experiment.aggregate
The seed is used to initialize pseudo-random functions. Some of the functions we will be using have pseudo-random elements. Setting a common seed ensures that all of us will get the same results, and that the results will remain stable when re-run.
set.seed(12345)
Display metadata by quantile
Using a few nested functions, we can produce prettier, more detailed, versions of the simple exploratory summary statistics we generated for the available metadata in the last section. In the code below, 5% quantile tables are produced for each metadata value, separated by sample identity.
Genes per cell
kable(do.call("cbind", tapply(experiment.aggregate$nFeature_RNA,
Idents(experiment.aggregate),quantile,probs=seq(0,1,0.05))),
caption = "5% Quantiles of Genes/Cell by Sample") %>% kable_styling()
| A001-C-007 | A001-C-104 | B001-A-301 | |
|---|---|---|---|
| 0% | 404.00 | 397.00 | 416.00 |
| 5% | 459.75 | 468.00 | 505.00 |
| 10% | 495.00 | 505.00 | 577.00 |
| 15% | 530.00 | 537.15 | 657.00 |
| 20% | 574.00 | 581.00 | 733.00 |
| 25% | 625.00 | 622.00 | 821.00 |
| 30% | 677.00 | 676.30 | 907.00 |
| 35% | 726.50 | 730.00 | 1006.55 |
| 40% | 783.00 | 791.00 | 1105.00 |
| 45% | 847.50 | 881.00 | 1228.85 |
| 50% | 927.00 | 959.00 | 1331.00 |
| 55% | 1020.25 | 1042.55 | 1447.00 |
| 60% | 1129.00 | 1127.00 | 1564.00 |
| 65% | 1241.50 | 1231.00 | 1680.00 |
| 70% | 1382.00 | 1358.00 | 1810.30 |
| 75% | 1570.00 | 1520.00 | 1955.75 |
| 80% | 1797.00 | 1718.80 | 2119.00 |
| 85% | 2101.75 | 1970.85 | 2325.00 |
| 90% | 2592.00 | 2344.70 | 2567.00 |
| 95% | 3735.00 | 3299.00 | 2959.35 |
| 100% | 12063.00 | 12064.00 | 8812.00 |
UMIs per cell
kable(do.call("cbind", tapply(experiment.aggregate$nCount_RNA,
Idents(experiment.aggregate),quantile,probs=seq(0,1,0.05))),
caption = "5% Quantiles of UMI/Cell by Sample") %>% kable_styling()
| A001-C-007 | A001-C-104 | B001-A-301 | |
|---|---|---|---|
| 0% | 500.00 | 500.00 | 500.00 |
| 5% | 536.00 | 542.00 | 595.00 |
| 10% | 576.50 | 588.00 | 698.00 |
| 15% | 632.25 | 631.00 | 803.00 |
| 20% | 686.00 | 692.20 | 909.00 |
| 25% | 756.00 | 746.00 | 1043.00 |
| 30% | 826.50 | 821.00 | 1179.90 |
| 35% | 894.25 | 899.00 | 1339.00 |
| 40% | 981.00 | 995.80 | 1514.60 |
| 45% | 1076.75 | 1106.45 | 1717.85 |
| 50% | 1203.50 | 1231.00 | 1913.00 |
| 55% | 1336.75 | 1371.00 | 2141.00 |
| 60% | 1493.00 | 1514.60 | 2399.00 |
| 65% | 1691.25 | 1686.30 | 2647.45 |
| 70% | 1920.00 | 1909.20 | 2960.10 |
| 75% | 2209.50 | 2188.25 | 3298.75 |
| 80% | 2690.00 | 2606.60 | 3704.00 |
| 85% | 3332.75 | 3124.00 | 4309.05 |
| 90% | 4419.00 | 4052.50 | 5087.10 |
| 95% | 7780.25 | 6615.75 | 6472.75 |
| 100% | 150805.00 | 149096.00 | 89743.00 |
Mitochondrial percentage per cell
kable(round(do.call("cbind", tapply(experiment.aggregate$percent_MT, Idents(experiment.aggregate),quantile,probs=seq(0,1,0.05))), digits = 3),
caption = "5% Quantiles of Percent Mitochondrial by Sample") %>% kable_styling()
| A001-C-007 | A001-C-104 | B001-A-301 | |
|---|---|---|---|
| 0% | 0.000 | 0.000 | 0.000 |
| 5% | 0.194 | 0.245 | 0.087 |
| 10% | 0.266 | 0.350 | 0.118 |
| 15% | 0.326 | 0.430 | 0.143 |
| 20% | 0.368 | 0.510 | 0.164 |
| 25% | 0.419 | 0.593 | 0.185 |
| 30% | 0.473 | 0.677 | 0.211 |
| 35% | 0.529 | 0.765 | 0.235 |
| 40% | 0.581 | 0.851 | 0.263 |
| 45% | 0.642 | 0.946 | 0.292 |
| 50% | 0.705 | 1.043 | 0.323 |
| 55% | 0.789 | 1.160 | 0.362 |
| 60% | 0.864 | 1.281 | 0.404 |
| 65% | 0.962 | 1.401 | 0.450 |
| 70% | 1.061 | 1.537 | 0.512 |
| 75% | 1.184 | 1.720 | 0.587 |
| 80% | 1.330 | 1.912 | 0.681 |
| 85% | 1.491 | 2.222 | 0.796 |
| 90% | 1.717 | 2.735 | 0.947 |
| 95% | 2.165 | 3.559 | 1.191 |
| 100% | 14.204 | 13.172 | 3.481 |
Visualize distribution of metadata values
Seurat has a number of convenient built-in functions for visualizing metadata. These functions produce ggplot objects, which can easily be modified using ggplot2. Of course, all of these visualizations can be reproduced with custom code as well, and we will include some examples of both modifying Seurat plots and generating plots from scratch as the analysis continues.
Violin plots
The VlnPlot function produces a composite plot with one panel for each element of the “features” vector. The data are grouped by the provided identity; by default, this is the active identity of the object, which can be accessed using the Idents() function, or in the “active.ident” slot.
VlnPlot(experiment.aggregate,
features = c("nFeature_RNA", "nCount_RNA","percent_MT"),
ncol = 1,
pt.size = 0.3)
Modifying Seurat plots
Modifying the ggplot objects produced by a Seurat plotting function works best on individual panels. Therefore, to recreate the function above with modifications, we can use lapply to create a list of plots. In some cases it may be more appropriate to create the plots individually so that different modifications can be applied to each plot.
lapply(c("nFeature_RNA", "nCount_RNA","percent_MT"), function(feature){
VlnPlot(experiment.aggregate,
features = feature,
pt.size = 0.01) +
scale_fill_viridis_d(option = "mako") # default colors are not colorblind-friendly
})
VlnPlot(experiment.aggregate, features = "nCount_RNA", pt.size = 0.01) +
scale_y_continuous(trans = "log10") +
scale_fill_viridis_d(option = "mako") +
ggtitle("log10(nCount_RNA)")
These can later be stitched together with another library, like patchwork, or cowplot.
Ridge plots
Ridge plots are very similar in appearance to violin plots turned on their sides. In some cases it may be more appropriate to create the plots individually so that appropriate transformations can be applied to each plot.
RidgePlot(experiment.aggregate, features="nFeature_RNA") +
scale_fill_viridis_d(option = "mako")
RidgePlot(experiment.aggregate, features="nCount_RNA") +
scale_x_continuous(trans = "log10") + # "un-squish" the distribution
scale_fill_viridis_d(option = "mako")
RidgePlot(experiment.aggregate, features="percent_MT") +
scale_fill_viridis_d(option = "mako") +
coord_cartesian(xlim = c(0, 10)) # zoom in on the lower end of the distribution
Custom plots
The Seurat built-in functions are useful and easy to interact with, but sometimes you may wish to visualize something for which a plotting function does not already exist. For example, we might want to see how many cells are expressing each gene over some UMI threshold.
The code below produces a ranked plot similar to the barcode inflection plots from the last section. On the x-axis are the genes arranged from most ubiquitously expressed to rarest. In a single cell dataset, many genes are expessed in a relatively small number of cells, or not at all. The y-axis displays the number of cells in which each gene is expressed.
Note: this function is SLOW You may want to skip this code block or run it while you take a break for a few minutes.
# retrieve count data
counts <- GetAssayData(experiment.aggregate)
# order genes from most to least ubiquitous
ranked.genes <- names(sort(Matrix::rowSums(counts >= 3), decreasing = TRUE))
# drop genes not expressed in any cell
ranked.genes <- ranked.genes[ranked.genes %in% names(which(Matrix::rowSums(counts >= 3) >= 1))]
# get number of cells for each sample
cell.counts <- sapply(ranked.genes, function(gene){
tapply(counts[gene,], experiment.aggregate$orig.ident, function(x){
sum(x >= 3)
})
})
cell.counts <- as.data.frame(t(cell.counts))
cell.counts$gene <- rownames(cell.counts)
cell.counts <- melt(cell.counts, variable.name = "sample", value.name = "count")
cell.counts$rank <- match(ranked.genes, cell.counts$gene)
# plot
ggplot(cell.counts, mapping = aes(x = rank, y = count, color = sample)) +
scale_x_continuous(trans = "log10") +
scale_y_continuous(trans = "log10") +
geom_smooth() +
scale_color_viridis_d(option = "mako") +
theme_classic() +
theme(legend.title = element_blank())
rm(counts, ranked.genes, cell.counts)
Scatter plots
Scatter plots allow us to visualize the relationships between the metadata variables.
# mitochondrial vs UMI
FeatureScatter(experiment.aggregate,
feature1 = "nCount_RNA",
feature2 = "percent_MT",
shuffle = TRUE) +
scale_color_viridis_d(option = "mako")
# mitochondrial vs genes
FeatureScatter(experiment.aggregate,
feature1 = "nFeature_RNA",
feature2 = "percent_MT",
shuffle = TRUE) +
scale_color_viridis_d(option = "mako")
# genes vs UMI
FeatureScatter(experiment.aggregate,
feature1 = "nCount_RNA",
feature2 = "nFeature_RNA",
shuffle = TRUE) +
scale_color_viridis_d(option = "mako")
Cell filtering
The goal of cell filtering is to remove cells with anomolous expression profiles, typically low UMI cells, which may correspond to low-quality cells or background barcodes that made it through the Cell Ranger filtration algorithm. It may also be appropriate to remove outlier cells with extremely high UMI counts.
In this case, the proposed cut-offs on the high end of the distributions are quite conservative, in part to reduce the size of the object and speed up analysis during the workshop.
The plots below display proposed filtering cut-offs.
FeatureScatter(experiment.aggregate,
feature1 = "nCount_RNA",
feature2 = "percent_MT",
shuffle = TRUE) +
geom_vline(xintercept = c(1000, 25000)) +
geom_hline(yintercept = 5) +
scale_color_viridis_d(option = "mako")
FeatureScatter(experiment.aggregate,
feature1 = "nFeature_RNA",
feature2 = "percent_MT",
shuffle = TRUE) +
geom_vline(xintercept = c(500, 7500)) +
geom_hline(yintercept = 5) +
scale_color_viridis_d(option = "mako")
FeatureScatter(experiment.aggregate,
feature1 = "nCount_RNA",
feature2 = "nFeature_RNA",
pt.size = 0.5,
shuffle = TRUE) +
geom_vline(xintercept = c(1000, 25000)) +
geom_hline(yintercept = c(500, 7500)) +
scale_color_viridis_d(option = "mako")
These filters can be put in place with the subset function.
table(experiment.aggregate$orig.ident)
# mitochondrial filter
experiment.filter <- subset(experiment.aggregate, percent_MT <= 5)
# UMI filter
experiment.filter <- subset(experiment.filter, nCount_RNA >= 1000 & nCount_RNA <= 25000)
# gene filter
experiment.filter <- subset(experiment.filter, nFeature_RNA >= 500 & nFeature_RNA <= 7500)
# filtering results
experiment.filter
table(experiment.filter$orig.ident)
# ridge plots
RidgePlot(experiment.filter, features="nFeature_RNA") +
scale_fill_viridis_d(option = "mako")
RidgePlot(experiment.filter, features="nCount_RNA") +
scale_x_continuous(trans = "log10") +
scale_fill_viridis_d(option = "mako")
RidgePlot(experiment.filter, features="percent_MT") +
scale_fill_viridis_d(option = "mako")
# use filtered results from now on
experiment.aggregate <- experiment.filter
rm(experiment.filter)
Play with the filtering parameters, and see how the results change. Is there a set of parameters you feel is more appropriate? Why?
Feature filtering
When creating the base Seurat object, we had the opportunity filter out some genes using the “min.cells” argument. At the time, we set that to 0. Since we didn’t filter our features then, we can apply a filter at this point, or if we did filter when the object was created, this would be an opportunity to be more aggressive with filtration. The custom code below provides a function that filters genes requiring a min.umi in at least min.cells, or takes a user-provided list of genes.
# define function
FilterGenes <- function(object, min.umi = NA, min.cells = NA, genes = NULL) {
genes.use = NA
if (!is.null(genes)) {
genes.use = intersect(rownames(object), genes)
} else if (min.cells & min.umi) {
num.cells = Matrix::rowSums(GetAssayData(object) >= min.umi)
genes.use = names(num.cells[which(num.cells >= min.cells)])
}
object = object[genes.use,]
object = LogSeuratCommand(object = object)
return(object)
}
# apply filter
experiment.filter <- FilterGenes(object = experiment.aggregate, min.umi = 2, min.cells = 10)
# filtering results
experiment.filter
experiment.aggregate <- experiment.filter
rm(experiment.filter)
Normalize the data
After filtering, the next step is to normalize the data. We employ a global-scaling normalization method LogNormalize that normalizes the gene expression measurements for each cell by the total expression, multiplies this by a scale factor (10,000 by default), and then log-transforms the data.
?NormalizeData
experiment.aggregate <- NormalizeData(
object = experiment.aggregate,
normalization.method = "LogNormalize",
scale.factor = 10000)
Cell cycle assignment
Cell cycle phase can be a significant source of variation in single cell and single nucleus experiments. There are a number of automated cell cycle stage detection methods available for single cell data. For this workshop, we will be using the built-in Seurat cell cycle function, CellCycleScoring. This tool compares gene expression in each cell to a list of cell cycle marker genes and scores each barcode based on marker expression. The phase with the highest score is selected for each barcode. Seurat includes a list of cell cycle genes in human single cell data.
s.genes <- cc.genes$s.genes
g2m.genes <- cc.genes$g2m.genes
For other species, a user-provided gene list may be substituted, or the orthologs of the human gene list used instead.
Do not run the code below for human experiments!
# mouse code DO NOT RUN for human data
convertHumanGeneList <- function(x){
require("biomaRt")
human = useEnsembl("ensembl",
dataset = "hsapiens_gene_ensembl",
mirror = "uswest")
mouse = useEnsembl("ensembl",
dataset = "mmusculus_gene_ensembl",
mirror = "uswest")
genes = getLDS(attributes = c("hgnc_symbol"),
filters = "hgnc_symbol",
values = x ,
mart = human,
attributesL = c("mgi_symbol"),
martL = mouse,
uniqueRows=T)
humanx = unique(genes[, 2])
print(head(humanx)) # print first 6 genes found to the screen
return(humanx)
}
# convert lists to mouse orthologs
s.genes <- convertHumanGeneList(cc.genes.updated.2019$s.genes)
g2m.genes <- convertHumanGeneList(cc.genes.updated.2019$g2m.genes)
Once an appropriate gene list has been identified, the CellCycleScoring function can be run.
experiment.aggregate <- CellCycleScoring(experiment.aggregate,
s.features = s.genes,
g2m.features = g2m.genes,
set.ident = TRUE)
table(experiment.aggregate@meta.data$Phase) %>%
kable(caption = "Number of Cells in each Cell Cycle Stage",
col.names = c("Stage", "Count"),
align = "c") %>%
kable_styling()
| Stage | Count |
|---|---|
| G1 | 3759 |
| G2M | 1136 |
| S | 1417 |
rm(s.genes, g2m.genes)
Because the “set.ident” argument was set to TRUE (this is also the default behavior), the active identity of the Seurat object was changed to the phase. To return the active identity to the sample identity, use the Idents function.
table(Idents(experiment.aggregate))
Idents(experiment.aggregate) <- "orig.ident"
table(Idents(experiment.aggregate))
Identify variable genes
The function FindVariableFeatures identifies the most highly variable genes (default 2000 genes) by fitting a line to the relationship of log(variance) and log(mean) using loess smoothing, uses this information to standardize the data, then calculates the variance of the standardized data. This helps avoid selecting genes that only appear variable due to their expression level.
?FindVariableFeatures
experiment.aggregate <- FindVariableFeatures(
object = experiment.aggregate,
selection.method = "vst")
length(VariableFeatures(experiment.aggregate))
top10 <- head(VariableFeatures(experiment.aggregate), 10)
top10
vfp1 <- VariableFeaturePlot(experiment.aggregate)
vfp1 <- LabelPoints(plot = vfp1, points = top10, repel = TRUE)
vfp1
How do the results change if you use selection.method = “dispersion” or selection.method = “mean.var.plot”?
FindVariableFeatures isn’t the only way to set the “variable features” of a Seurat object. Another reasonable approach is to select a set of “minimally expressed” genes.
min.value <- 2
min.cells <- 10
num.cells <- Matrix::rowSums(GetAssayData(experiment.aggregate, slot = "count") > min.value)
genes.use <- names(num.cells[which(num.cells >= min.cells)])
length(genes.use)
VariableFeatures(experiment.aggregate) <- genes.use
rm(min.value, min.cells, num.cells, genes.use)
Save the Seurat object and download the next Rmd file
saveRDS(experiment.aggregate, file="scRNA_workshop_2.rds")
download.file("https://raw.githubusercontent.com/ucdavis-bioinformatics-training/2023-June-Single-Cell-RNA-Seq-Analysis/main/data_analysis/scRNA_Workshop-PART3.Rmd", "scRNA_Workshop-PART3.Rmd")
Session Information
sessionInfo()