Chapter 3 Feature selection
3.1 Motivation
We often use scRNA-seq data in exploratory analyses to characterize heterogeneity across the cell population. Procedures like clustering (Chapter 6) and dimensionality reduction (Chapter 4) compare cells based on their gene expression profiles, typically by aggregating per-gene differences into a single (dis)similarity metric between a pair of cells. The choice of genes in this calculation has a major impact on the behavior of the metric and the performance of downstream methods. We want to select genes that contain useful information about the biology of the system while removing genes that contain random noise. In addition to improving signal, this reduces the size of the data to improve computational efficiency of later steps.
3.2 Selecting highly variable genes
3.2.1 Modelling the mean-variance trend
Our aim is to select the most highly variable genes (HVGs) based on their expression across the population. This assumes that genuine biological differences will manifest as increased variation in the affected genes, compared to other genes that are only affected by technical noise or a baseline level of “uninteresting” biological variation (e.g., from transcriptional bursting). To demonstrate, we’ll load the classic PBMC dataset from 10X Genomics (Zheng et al. 2017):
# Loading in raw data from the 10X output files.
library(DropletTestFiles)
raw.path.10x <- getTestFile("tenx-2.1.0-pbmc4k/1.0.0/filtered.tar.gz")
dir.path.10x <- file.path(tempdir(), "pbmc4k")
untar(raw.path.10x, exdir=dir.path.10x)
library(DropletUtils)
fname.10x <- file.path(dir.path.10x, "filtered_gene_bc_matrices/GRCh38")
sce.10x <- read10xCounts(fname.10x, col.names=TRUE)
# Applying our default QC with outlier-based thresholds.
library(scrapper)
is.mito.10x <- grepl("^MT-", rowData(sce.10x)$Symbol)
sce.qc.10x <- quickRnaQc.se(sce.10x, subsets=list(MT=is.mito.10x))
sce.qc.10x <- sce.qc.10x[,sce.qc.10x$keep]
# Computing log-normalized expression values.
sce.norm.10x <- normalizeRnaCounts.se(sce.qc.10x, size.factors=sce.qc.10x$sum)
sce.norm.10x## class: SingleCellExperiment
## dim: 33694 4147
## metadata(2): Samples qc
## assays(2): counts logcounts
## rownames(33694): ENSG00000243485 ENSG00000237613 ... ENSG00000277475
## ENSG00000268674
## rowData names(2): ID Symbol
## colnames(4147): AAACCTGAGACAGACC-1 AAACCTGAGCGCCTCA-1 ...
## TTTGTCAGTTAAGACA-1 TTTGTCATCCCAAGAT-1
## colData names(7): Sample Barcode ... keep sizeFactor
## reducedDimNames(0):
## mainExpName: NULL
## altExpNames(0):We quantify per-gene variation by computing the variance of the log-normalized expression values for each gene across all cells (Lun, McCarthy, and Marioni 2016). This is motivated by the use of log-transformed values for quantifying differences between cells in downstream steps like clustering and dimensionality reduction (Section 2.3.1). Genes with the largest variances in log-space will contribute most to these differences and thus be most informative for distinguishing subpopulations. Alternatively, other metrics based on untransformed data (e.g., the coefficient of variation) could be used here but may focus on genes with smaller average log-fold changes.
To account the mean-variance relationship in the log-transformed expression values, we fit a trend to the variances with respect to the mean (Figure 3.1). The hump-like shape of this trend is typical of log-transformed RNA-seq counts. Specifically, the variance increases linearly as the mean increases from zero, as larger variances are obviously possible when the counts are not all equal to zero. On the other hand, the relative contribution of sampling noise decreases at high abundances, resulting in a downward trend to a plateau (Law et al. 2014). The peak represents the point at which these two competing effects cancel each other out.
sce.var.10x <- chooseRnaHvgs.se(sce.norm.10x)
# Let's have a peek at the statistics for the top HVGs.
rd.10x <- rowData(sce.var.10x)
ordered.residual.10x <- order(rd.10x$residuals, decreasing=TRUE)
rd.10x[head(ordered.residual.10x),c("Symbol", "means", "variances", "fitted", "residuals", "hvg")]## DataFrame with 6 rows and 6 columns
## Symbol means variances fitted residuals hvg
## <character> <numeric> <numeric> <numeric> <numeric> <logical>
## ENSG00000090382 LYZ 1.86235 4.92849 0.757074 4.17141 TRUE
## ENSG00000163220 S100A9 1.82739 4.34678 0.759089 3.58769 TRUE
## ENSG00000143546 S100A8 1.60171 4.18943 0.762348 3.42708 TRUE
## ENSG00000204287 HLA-DRA 2.08390 3.74121 0.735784 3.00542 TRUE
## ENSG00000019582 CD74 2.89021 3.41935 0.608014 2.81134 TRUE
## ENSG00000101439 CST3 1.40312 2.89004 0.753198 2.13684 TRUE## [1] 4000plot(rd.10x$means, rd.10x$variances,
col=ifelse(is.hvg.10x, "black", "grey"), pch=16, cex=1,
xlab="Mean of log-expression", ylab="Variance of log-expression")
legend("topright", col=c("black", "grey"), pch=16, legend=c("HVG", "not HVG"))
# Just using approxfun() to make a nice-looking curve for us.
trend.10x <- approxfun(rd.10x$means, rd.10x$fitted)
curve(trend.10x, add=TRUE, col="dodgerblue", lwd=2)
Figure 3.1: Variance of the log-normalized expression values across all genes in the PBMC data set, as a function of the mean. Each point represents a gene, colored according to whether it was chosen as a HVG. The blue line represents the trend fitted to all genes.
We define HVGs as the top \(H\) genes with the largest residuals above the trend, where \(H\) is typically between 1000 and 5000 (defaulting to 4000). Our assumption is that, at any given mean, the variation in expression for most genes is driven by uninteresting processes like sampling noise. The fitted value of the trend at any given gene’s mean represents a mean-dependent estimate of its uninteresting variation, while the residuals represent the “interesting” variation for each gene and can be used as the metric for HVG selection. By comparison, if we just used each gene’s variance directly, the choice of HVGs would be driven more by the gene’s abundance than its biological heterogeneity. (In other words, the log-transformation is not a variance-stabilizing transformation7.) This would bias us against lower-abundance genes that exhibit increased biological variation.
Once we have our top HVGs, we can use them in downstream steps like principal components analysis. We’ll discuss this more in Chapter 4, but it’s as simple as using only the subset of HVGs in the analysis:
3.2.2 Choosing the number of HVGs
How many HVGs should we use in our downstream analyses, i.e., what is the “best” value of \(H\)? A larger set of HVGs will reduce the risk of discarding interesting biological signal by retaining more potentially relevant genes, at the cost of adding noise from irrelevant genes that might obscure that signal. It’s difficult to determine the optimal trade-off for any given application as the distinction between noise and signal is context-dependent. For example, variation in the activation status of certain immune cells may not be interesting when we only want to identify the cell types; the former can even interfere with the latter by encouraging the formation of clusters based on activation strength instead.
Our recommendation is to simply pick a “reasonable” \(H\) - usually somewhere between 1000 and 5000 - and proceed with the rest of the analysis. If we can answer our scientific question, then our choice is good enough; if not, we can just try another value. There’s nothing wrong with trying different parameters during data exploration8. In fact, different choices of \(H\) can provide new perspectives of the same dataset by changing the balance between signal and noise, so we might discover new population structure that would not be apparent with other parameters. Don’t spend too much time worrying about obtaining the “optimal” value.
If we really want to ensure that all biological structure is preserved, we could define the set of HVGs as all genes with variances above the trend. This avoids any judgement calls about the definition of “interesting” variation, giving an opportunity for weaker population structure to manifest. It is most useful for rare and/or weakly-separated subpopulations where the relevant marker genes are not variable enough to sneak into the top \(H\) genes. The obvious cost is that more noise is also captured, which can reduce the resolution of subpopulations; and we need to perform more computational work in each downstream step, as more genes are involved.
# Setting top=Inf to select all genes with positive residuals.
hvgs.all.10x <- chooseHighlyVariableGenes(rd.10x$residuals, top=Inf)
length(hvgs.all.10x)## [1] 47383.3 Blocking on uninteresting factors
Larger datasets may contain multiple blocks of cells that exhibit uninteresting differences in gene expression, e.g., batch effects, variability between donors. We are not interested in HVGs that are driven by these differences; instead, we want to focus on genes that are highly variable within each block. Let’s demonstrate using some trophoblast scRNA-seq data generated across two plates (Lun et al. 2017):
library(scRNAseq)
sce.tropho <- LunSpikeInData("tropho")
table(sce.tropho$block) # i.e., the plate of origin. ##
## 20160906 20170201
## 96 96# Computing the QC metrics.
library(scrapper)
is.mito.tropho <- which(any(seqnames(rowRanges(sce.tropho))=="MT"))
sce.qc.tropho <- quickRnaQc.se(
sce.tropho,
subsets=list(MT=is.mito.tropho),
altexp.proportions="ERCC",
block=sce.tropho$block
)
sce.qc.tropho <- sce.qc.tropho[,sce.qc.tropho$keep]
# Computing log-normalized expression values.
sce.norm.tropho <- normalizeRnaCounts.se(
sce.qc.tropho,
size.factors=sce.qc.tropho$sum,
block=sce.qc.tropho$block
)Setting block= instructs chooseRnaHvgs.se() to compute the mean and variance for each gene within each plate.
This ensures that any systematic technical differences between plates (e.g., in sequencing depth) will not inflate the variance estimates.
It will also fit a separate trend for each plate, which accommodates differences in the mean-variance relationships between plates.
In this case, we observe only minor differences between the trends in Figure 3.2, which indicates that the experiment was tightly replicated across plates.
sce.var.tropho <- chooseRnaHvgs.se(
sce.norm.tropho,
block=sce.qc.tropho$block,
include.per.block=TRUE # only needed for plotting.
)
per.block <- rowData(sce.var.tropho)$per.block
par(mfrow=c(1,2))
for (block in colnames(per.block)) {
current <- per.block[[block]]
plot(current$means, current$variances,
xlab="Mean of log-expression",
ylab="Variance of log-expression",
main=block, pch=16, cex=0.5)
trend <- approxfun(current$means, current$fitted)
curve(trend, add=TRUE, col="dodgerblue", lwd=2)
}
Figure 3.2: Variance of the log-normalized expression values across all genes in the trophoblast data set, as a function of the mean after blocking on the plate of origin. Each plot represents the results for a single plate. Each point represents a gene and the fitted trend is shown in blue.
chooseRnaHVgs.se() also combines information across blocks by reporting the weighted mean of each statistic,
where the weight is determined by the size of number of cells in each plate.
We use the mean residuals to select our top HVGs as described previously.
This ensures that each block contributes some information about the variability of each gene.
High variability in any block can increase the residual for a gene, giving it an opportunity to be selected as a HVG.
## DataFrame with 46603 rows and 4 columns
## means variances fitted residuals
## <numeric> <numeric> <numeric> <numeric>
## ENSMUSG00000102693 0.0000000 0.000000 0.000000 0.000000
## ENSMUSG00000064842 0.0000000 0.000000 0.000000 0.000000
## ENSMUSG00000051951 0.0351964 0.190773 0.141455 0.049318
## ENSMUSG00000102851 0.0000000 0.000000 0.000000 0.000000
## ENSMUSG00000103377 0.0993010 0.536793 0.380165 0.156627
## ... ... ... ... ...
## ENSMUSG00000094431 0 0 0 0
## ENSMUSG00000094621 0 0 0 0
## ENSMUSG00000098647 0 0 0 0
## ENSMUSG00000096730 0 0 0 0
## ENSMUSG00000095742 0 0 0 0## [1] 4000Alternatively, we could focus on genes are consistently variable within each block by asking chooseRnaHvgs.se() to compute a quantile instead of a weighted mean.
For example, we could report the minimum residual across blocks, which means that genes will only be considered as HVGs if they have large positive residuals in each block.
This tends to scale poorly as it becomes too stringent with a large number of blocks.
sce.var.min.tropho <- chooseRnaHvgs.se(
sce.norm.tropho,
block=sce.qc.tropho$block,
more.var.args=list(
block.average.policy="quantile",
block.quantile=0 # i.e., minimum.
)
)
rowData(sce.var.min.tropho)[,c("means", "variances", "fitted", "residuals")] # minimum across blocks.## DataFrame with 46603 rows and 4 columns
## means variances fitted residuals
## <numeric> <numeric> <numeric> <numeric>
## ENSMUSG00000102693 0.0000000 0.000000 0.000000 0.0000
## ENSMUSG00000064842 0.0000000 0.000000 0.000000 0.0000
## ENSMUSG00000051951 0.0000000 0.000000 0.000000 0.0000
## ENSMUSG00000102851 0.0000000 0.000000 0.000000 0.0000
## ENSMUSG00000103377 0.0780524 0.469098 0.311943 0.1561
## ... ... ... ... ...
## ENSMUSG00000094431 0 0 0 0
## ENSMUSG00000094621 0 0 0 0
## ENSMUSG00000098647 0 0 0 0
## ENSMUSG00000096730 0 0 0 0
## ENSMUSG00000095742 0 0 0 0It is generally expected that block= will be used for uninteresting factors of variation.
In this case, the plate of origin is a technical factor that should be ignored.
However, imagine instead that each plate corresponds to a different treatment condition.
In such cases, we might not use block= to ensure that our variance estimates can capture the differences between treatments.
This decision is discussed in more detail in Chapter 8.
3.4 Customizing the trend fit
The trend fit in chooseRnaHvgs.se() is based on the LOWESS non-parametric smoother (Cleveland 1979) with some modifications.
LOWESS slides a window across the x-coordinates and performs a linear regression within each window to obtain the fitted value for the point at the window’s center.
The size of the window varies between points, expanding or contracting until it contains a specified number/proportion of all points in the dataset.
Standard LOWESS mostly works well but is suboptimal in x-axis intervals that contain very few points - hence, the modifications.
To demonstrate, let’s have a look at a human pancreas dataset from Segerstolpe et al. (2016):
library(scRNAseq)
sce.seger <- SegerstolpePancreasData()
# For simplicity, we'll focus on one of the donors.
sce.seger <- sce.seger[,sce.seger$individual=="H2"]
# For reasons unknown to us, the data supplied by the authors contain
# duplicated row names, so we'll just get rid of those to avoid confusion.
sce.seger <- sce.seger[!duplicated(rownames(sce.seger)),]
# Running QC. Seems like they don't have any data for the mitochondrial genes,
# unfortunately, but they do have spike-ins so we'll just use those instead.
library(scrapper)
sce.qc.seger <- quickRnaQc.se(sce.seger, subsets=list(), altexp.proportions="ERCC")
# Computing log-normalized expression values.
sce.norm.seger <- normalizeRnaCounts.se(sce.qc.seger, size.factors=sce.qc.seger$sum)The drawbacks of standard LOWESS manifest in typical scRNA-seq datasets where there are very few genes at high abundances.
This forces the LOWESS window to expand to contain more points, reducing sensitivity of the fitted trend to the behavior of the high-abundance genes (Figure 3.3).
In chooseRnaHvgs.se(), we modify the definition of each window to (i) contain fewer points and (ii) have a minimum width along the x-axis.
The former improves sensitivity in sparse intervals while the latter reduces the risk of overfitting in dense intervals.
This results in a fitted trend that is more faithful to the high-abundance genes.
sce.var.seger <- chooseRnaHvgs.se(sce.norm.seger)
rd.seger <- rowData(sce.var.seger)
plot(rd.seger$means, rd.seger$variances,
xlab="Mean of log-expression", ylab="Variance of log-expression", pch=16, cex=0.5)
trend.seger <- approxfun(rd.seger$means, rd.seger$fitted)
curve(trend.seger, add=TRUE, col="dodgerblue", lwd=2)
# fitVarianceTrend() is the underlying function used by chooseRnaHvgs.se() to
# fit the trend; we set use.min.width=FALSE to forcibly use standard LOWESS.
fit.standard.seger <- fitVarianceTrend(rd.seger$means, rd.seger$variances, use.min.width=FALSE)
trend.standard.seger <- approxfun(rd.seger$means, fit.standard.seger$fitted)
curve(trend.standard.seger, add=TRUE, col="salmon", lwd=2)
legend("topright", lwd=2, col=c("dodgerblue", "salmon"), legend=c("modified", "standard"))
Figure 3.3: Variance of the log-normalized expression values across all genes in one donor of the Segerstople pancreas data set, as a function of the mean. Each point represents a gene while the lines represent trends fitted with standard (red) or modified LOWESS (blue).
It’s worth noting that many of the default settings in chooseRnaHvgs.se() are specifically tuned for scRNA-seq data.
In particular, we assume that the input data is on the same scale as the log-transformed counts, and that the number of available genes is on the order of several thousands or more.
These assumptions may not be valid in other contexts where standard LOWESS may be more appropriate.
For example, we use standard LOWESS to fit a mean-variance trend to spike-in data in Section 3.6,
where the chooseRnaHvgs.se() defaults are not suitable due to the lower number of spike-in transcripts.
All that said, the precise parametrization of the trend fitting doesn’t actually matter all that much. There are so few genes in these sparse intervals that their (lack of) selection as HVGs won’t have a major effect on downstream analyses. But sometimes it’s just nice to look at some well-fitted curves.
3.5 Selecting a priori genes of interest
A blunt yet effective feature selection strategy is to use pre-defined sets of interesting genes. The aim is to focus on specific aspects of biological heterogeneity that may be masked by other factors when using unsupervised methods for HVG selection. For example, to study transcriptional changes during the earliest stages of cell fate commitment (Messmer et al. 2019), we might focus only on lineage markers to avoid interference from variability in other pathways (e.g., cell cycle, metabolism). Using scRNA-seq data in this manner is conceptually equivalent to a fluorescence activated cell sorting (FACS) experiment, with the convenience of being able to (re)define the features of interest at any time. We provide some examples of a priori selection based on MSigDB gene sets (Liberzon et al. 2015) below:
library(msigdbr)
c7.sets <- msigdbr(species = "Homo sapiens", category = "C7")
head(unique(c7.sets$gs_name))## [1] "ANDERSON_BLOOD_CN54GP140_ADJUVANTED_WITH_GLA_AF_AGE_18_45YO_1DY_DN"
## [2] "ANDERSON_BLOOD_CN54GP140_ADJUVANTED_WITH_GLA_AF_AGE_18_45YO_1DY_UP"
## [3] "ANDERSON_BLOOD_CN54GP140_ADJUVANTED_WITH_GLA_AF_AGE_18_45YO_3DY_DN"
## [4] "ANDERSON_BLOOD_CN54GP140_ADJUVANTED_WITH_GLA_AF_AGE_18_45YO_3DY_UP"
## [5] "ANDERSON_BLOOD_CN54GP140_ADJUVANTED_WITH_GLA_AF_AGE_18_45YO_6HR_DN"
## [6] "ANDERSON_BLOOD_CN54GP140_ADJUVANTED_WITH_GLA_AF_AGE_18_45YO_6HR_UP"# Using the Goldrath sets to distinguish CD8 subtypes
cd8.sets <- c7.sets[grep("GOLDRATH", c7.sets$gs_name),]
cd8.genes <- rownames(sce.10x) %in% cd8.sets$ensembl_gene
summary(cd8.genes)## Mode FALSE TRUE
## logical 32851 843# Using GSE11924 to distinguish between T helper subtypes
th.sets <- c7.sets[grep("GSE11924", c7.sets$gs_name),]
th.genes <- rownames(sce.10x) %in% th.sets$ensembl_gene
summary(th.genes)## Mode FALSE TRUE
## logical 31722 1972# Using GSE11961 to distinguish between B cell subtypes
b.sets <- c7.sets[grep("GSE11961", c7.sets$gs_name),]
b.genes <- rownames(sce.10x) %in% b.sets$ensembl_gene
summary(b.genes)## Mode FALSE TRUE
## logical 27995 5699Don’t be ashamed to take advantage of prior biological knowledge during feature selection to address specific hypotheses! We say this because a common refrain in genomics is that the data analysis should be “unbiased”, i.e., free from any biological preconceptions. Which is fine and all, but such “biases” are already present at every stage, starting with experimental design and ending with the interpretation of the data. So if we already know what we’re looking for, why not make life simpler and just go for it? Of course, the downside of focusing on pre-defined genes is that it will limit our capacity to detect novel or unexpected aspects of variation. Thus, this kind of focused analysis should be complementary to (rather than a replacement for) the unsupervised feature selection strategies discussed above.
We can also invert this reasoning to remove genes that are unlikely to be of interest prior to downstream analyses. This eliminates unwanted variation that could mask relevant biology and interfere with interpretation of the results. Ribosomal protein genes or mitochondrial genes are common candidates for removal, especially in situations with varying levels of cell damage within a population. For immune cell subsets, we might also be inclined to remove immunoglobulin genes and T cell receptor genes for which clonal expression introduces irrelevant population structure.
# Identifying ribosomal proteins:
ribo.discard <- grepl("^RP[SL]\\d+", rowData(sce.10x)$Symbol)
sum(ribo.discard)## [1] 99# A more curated approach for identifying ribosomal protein genes:
c2.sets <- msigdbr(species = "Homo sapiens", category = "C2")
ribo.set <- c2.sets[c2.sets$gs_name=="KEGG_RIBOSOME",]$ensembl_gene
ribo.discard <- rownames(sce.10x) %in% ribo.set
sum(ribo.discard)## [1] 87library(AnnotationHub)
edb <- AnnotationHub()[["AH73881"]]
anno <- select(edb, keys=rowData(sce.10x)$ID, keytype="GENEID",
columns="TXBIOTYPE")
# Removing immunoglobulin variable chains:
igv.set <- anno$GENEID[anno$TXBIOTYPE %in% c("IG_V_gene", "IG_V_pseudogene")]
igv.discard <- rownames(sce.10x) %in% igv.set
sum(igv.discard)## [1] 326# Removing TCR variable chains:
tcr.set <- anno$GENEID[anno$TXBIOTYPE %in% c("TR_V_gene", "TR_V_pseudogene")]
tcr.discard <- rownames(sce.10x) %in% tcr.set
sum(tcr.discard)## [1] 138We tend to err on the side of caution and abstain from preemptive filtering on biological function until these genes are demonstrably problematic in downstream analyses.
3.6 Quantifying technical noise
Back in the old days, everyone was obsessed with modelling the gene-wise variability in scRNA-seq data (Brennecke et al. 2013; Vallejos, Marioni, and Richardson 2015; Kim et al. 2015). Spike-in transcripts were critical to this effort as they allowed us to decompose each gene’s variance into technical and biological components. As spike-ins should not be subject to biological effects, the variance in spike-in expression could be used as an estimate of the technical component. Subtracting the spike-in variance from the variance of an endogenous gene at a similar abundance would yield an estimate of the biological component. Sadly, those days are gone and people don’t care about variance decomposition anymore. But for old times’ sake, we’ll demonstrate how to do this with the Zeisel et al. (2015) dataset:
library(scRNAseq)
sce.zeisel <- ZeiselBrainData()
is.mito.zeisel <- rowData(sce.zeisel)$featureType=="mito"
# Performing some QC to set up the dataset prior to normalization.
library(scrapper)
sce.qc.zeisel <- quickRnaQc.se(sce.zeisel, subsets=list(MT=is.mito.zeisel), altexp.proportions="ERCC")
sce.qc.zeisel <- sce.qc.zeisel[,sce.qc.zeisel$keep]We compute log-normalized expression values for endogenous genes and spike-in transcripts with their respective size factors.
Unlike Section 2.5, we still use the library size factors for the endogenous genes as we are not currently interested in changes in total RNA content.
Both sets of size factors are centered to preserve the scale of the original counts, ensuring that normalized abundances are comparable between genes and spike-ins.
(This is a bit more complicated with blocking, as the mean of the spike-in factors within each block must be scaled to the mean of the library size factors in that block;
this is handled by setting block= in normalizeRnaCountsWithSpikeIns.se(), but for brevity, we won’t show that here.)
sce.norm.zeisel <- normalizeRnaCountsWithSpikeIns.se(
sce.qc.zeisel,
spike.altexps="ERCC",
use.spike.ins.for.endogenous=FALSE # spike-in factors only used to normalize spike-in counts.
)We fit a mean-dependent trend to the variances of the spike-in transcripts (Figure 3.4). At any given mean, the fitted value of the spike-in trend represents an estimate of the techical component of the variance. This assumes that an endogenous gene is subject to the same technical noise as a spike-in transcript of the same abundance. The residual from the trend represents each gene’s biological component of variation, which can be used to select HVGs as previously described.
sce.var.zeisel <- chooseRnaHvgsWithSpikeIns.se(sce.norm.zeisel, spike.altexp="ERCC")
var.gene.zeisel <- rowData(sce.var.zeisel)
summary(var.gene.zeisel$residuals)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -0.0439129 -0.0003047 0.0155817 0.1348219 0.1285275 15.0815244## Mode FALSE TRUE
## logical 16006 4000plot(var.gene.zeisel$means, var.gene.zeisel$variances,
col=ifelse(var.gene.zeisel$hvg, "black", "grey"),
xlab="Mean of log-expression", ylab="Variance of log-expression", pch=16, cex=0.5)
legend("topright", col=c("black", "grey"), pch=16, legend=c("HVG", "not HVG"))
var.spike.zeisel <- rowData(altExp(sce.var.zeisel, "ERCC"))
points(var.spike.zeisel$means, var.spike.zeisel$variances, col="dodgerblue", pch=4)
trend.spike.zeisel <- approxfun(var.spike.zeisel$means, var.spike.zeisel$fitted)
curve(trend.spike.zeisel, add=TRUE, col="dodgerblue", lwd=2)
# Fit a trend to the endogenous genes for comparison.
fit.endogenous.zeisel <- fitVarianceTrend(var.gene.zeisel$means, var.gene.zeisel$variances)
trend.endogenous.zeisel <- approxfun(var.gene.zeisel$means, fit.endogenous.zeisel$fitted)
curve(trend.endogenous.zeisel, add=TRUE, col="salmon", lwd=2, lty=2)
Figure 3.4: Variance of endogenous genes (black for HVGs, grey otherwise) and spike-in transcripts (blue crosses) in the Zeisel brain dataset, as a function of the mean. Trends are fitted to the spike-in transcripts (blue full) and endogenous genes (red dashed).
In practice, using the biological component from the spike-in trend doesn’t provide much benefit over the residuals from the trend fitted to the endogenous genes. Oh well.
Session information
## R version 4.6.1 (2026-06-24)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 24.04.4 LTS
##
## Matrix products: default
## BLAS: /home/biocbuild/bbs-3.24-bioc/R/lib/libRblas.so
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.12.0 LAPACK version 3.12.0
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_GB LC_COLLATE=C
## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## time zone: America/New_York
## tzcode source: system (glibc)
##
## attached base packages:
## [1] stats4 stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] AnnotationHub_4.3.2 BiocFileCache_3.3.0
## [3] dbplyr_2.6.0 msigdbr_26.1.0
## [5] ensembldb_2.37.3 AnnotationFilter_1.37.0
## [7] GenomicFeatures_1.65.0 AnnotationDbi_1.75.2
## [9] scRNAseq_2.27.0 scrapper_1.7.3
## [11] DropletUtils_1.33.0 SingleCellExperiment_1.35.2
## [13] SummarizedExperiment_1.43.0 Biobase_2.73.2
## [15] GenomicRanges_1.65.1 Seqinfo_1.3.0
## [17] IRanges_2.47.2 S4Vectors_0.51.6
## [19] BiocGenerics_0.59.10 generics_0.1.4
## [21] MatrixGenerics_1.25.0 matrixStats_1.5.0
## [23] DropletTestFiles_1.23.0
##
## loaded via a namespace (and not attached):
## [1] DBI_1.3.0 bitops_1.1-0
## [3] httr2_1.3.0 rlang_1.3.0
## [5] magrittr_2.0.5 otel_0.2.0
## [7] gypsum_1.9.0 compiler_4.6.1
## [9] RSQLite_3.53.3 DelayedMatrixStats_1.35.0
## [11] png_0.1-9 vctrs_0.7.3
## [13] ProtGenerics_1.45.0 pkgconfig_2.0.3
## [15] crayon_1.5.3 fastmap_1.2.0
## [17] XVector_0.53.0 scuttle_1.23.1
## [19] Rsamtools_2.29.0 rmarkdown_2.31
## [21] UCSC.utils_1.9.0 purrr_1.2.2
## [23] bit_4.6.0 xfun_0.60
## [25] cachem_1.1.0 beachmat_2.29.0
## [27] cigarillo_1.3.1 GenomeInfoDb_1.49.1
## [29] jsonlite_2.0.0 blob_1.3.0
## [31] rhdf5filters_1.25.2 DelayedArray_0.39.3
## [33] Rhdf5lib_2.1.0 BiocParallel_1.47.0
## [35] parallel_4.6.1 R6_2.6.1
## [37] bslib_0.11.0 limma_3.69.2
## [39] rtracklayer_1.73.0 jquerylib_0.1.4
## [41] assertthat_0.2.1 Rcpp_1.1.2
## [43] bookdown_0.47 knitr_1.51
## [45] R.utils_2.13.0 BiocBaseUtils_1.15.1
## [47] Matrix_1.7-6 tidyselect_1.2.1
## [49] abind_1.4-8 yaml_2.3.12
## [51] codetools_0.2-20 curl_7.1.0
## [53] alabaster.sce_1.13.0 lattice_0.22-9
## [55] tibble_3.3.1 withr_3.0.3
## [57] KEGGREST_1.53.6 evaluate_1.0.5
## [59] alabaster.schemas_1.13.0 ExperimentHub_3.3.1
## [61] Biostrings_2.81.6 pillar_1.11.1
## [63] BiocManager_1.30.27 filelock_1.0.3
## [65] RCurl_1.98-1.19 BiocVersion_3.24.0
## [67] alabaster.base_1.13.1 sparseMatrixStats_1.25.0
## [69] alabaster.ranges_1.13.0 glue_1.8.1
## [71] lazyeval_0.2.3 alabaster.matrix_1.13.0
## [73] tools_4.6.1 BiocIO_1.23.3
## [75] BiocNeighbors_2.7.2 GenomicAlignments_1.49.1
## [77] locfit_1.5-9.12 babelgene_22.9
## [79] XML_3.99-0.23 rhdf5_2.57.3
## [81] grid_4.6.1 edgeR_4.11.4
## [83] HDF5Array_1.41.0 restfulr_0.0.17
## [85] cli_3.6.6 rappdirs_0.3.4
## [87] S4Arrays_1.13.0 dplyr_1.2.1
## [89] alabaster.se_1.13.0 R.methodsS3_1.8.2
## [91] sass_0.4.10 digest_0.6.39
## [93] SparseArray_1.13.2 dqrng_0.4.1
## [95] rjson_0.2.23 memoise_2.0.1
## [97] htmltools_0.5.9 R.oo_1.27.1
## [99] lifecycle_1.0.5 h5mread_1.5.0
## [101] httr_1.4.8 statmod_1.5.2
## [103] bit64_4.8.2References
Brennecke, P., S. Anders, J. K. Kim, A. A. Kołodziejczyk, X. Zhang, V. Proserpio, B. Baying, et al. 2013. “Accounting for technical noise in single-cell RNA-seq experiments.” Nat. Methods 10 (11): 1093–5.
Cleveland, W. S. 1979. “Robust Locally Weighted Regression and Smoothing Scatterplots.” J. Am. Stat. Assoc. 74 (368): 829–36.
Kim, J. K., A. A. Kołodziejczyk, T. Illicic, S. A. Teichmann, and J. C. Marioni. 2015. “Characterizing noise structure in single-cell RNA-seq distinguishes genuine from technical stochastic allelic expression.” Nat. Commun. 6: 8687.
Law, C. W., Y. Chen, W. Shi, and G. K. Smyth. 2014. “voom: Precision weights unlock linear model analysis tools for RNA-seq read counts.” Genome Biol. 15 (2): R29.
Liberzon, A., C. Birger, H. Thorvaldsdóttir, M. Ghandi, J. P. Mesirov, and P. Tamayo. 2015. “The Molecular Signatures Database (MSigDB) hallmark gene set collection.” Cell. Syst. 1 (6): 417–25.
Lun, A. T. L., F. J. Calero-Nieto, L. Haim-Vilmovsky, B. Gottgens, and J. C. Marioni. 2017. “Assessing the reliability of spike-in normalization for analyses of single-cell RNA sequencing data.” Genome Res. 27 (11): 1795–1806.
Lun, A. T. L., D. J. McCarthy, and J. C. Marioni. 2016. “A step-by-step workflow for low-level analysis of single-cell RNA-seq data.” F1000Res. 5 (August).
Messmer, T., F. von Meyenn, A. Savino, F. Santos, H. Mohammed, A. T. L. Lun, J. C. Marioni, and W. Reik. 2019. “Transcriptional heterogeneity in naive and primed human pluripotent stem cells at single-cell resolution.” Cell. Rep. 26 (4): 815–24.
Segerstolpe, A., A. Palasantza, P. Eliasson, E. M. Andersson, A. C. Andreasson, X. Sun, S. Picelli, et al. 2016. “Single-cell transcriptome profiling of human pancreatic islets in health and type 2 diabetes.” Cell Metab. 24 (4): 593–607.
Vallejos, C. A., J. C. Marioni, and S. Richardson. 2015. “BASiCS: Bayesian analysis of single-cell sequencing data.” PLoS Comput. Biol. 11 (6): e1004333.
Zeisel, A., A. B. Munoz-Manchado, S. Codeluppi, P. Lonnerberg, G. La Manno, A. Jureus, S. Marques, et al. 2015. “Brain structure. Cell types in the mouse cortex and hippocampus revealed by single-cell RNA-seq.” Science 347 (6226): 1138–42.
Zheng, G. X., J. M. Terry, P. Belgrader, P. Ryvkin, Z. W. Bent, R. Wilson, S. B. Ziraldo, et al. 2017. “Massively parallel digital transcriptional profiling of single cells.” Nat. Commun. 8 (January): 14049.
And even VSTs need to deal with a trend at small counts. A gene with a mean of zero will also have a variance of zero, and the variance has to become positive with increasing mean, so at some point, a trend will manifest.↩︎
It would be another story if we were doing some confirmatory analysis with rigorous hypothesis testing. In such cases, it would be improper to shop around for the best parameters that gives us the result that we want. Fortunately, single-cell analyses have looser standards as we often don’t know what we’re looking for.↩︎