# Chapter 10 Feature selection

## 10.1 Motivation

We often use scRNA-seq data in exploratory analyses to characterize heterogeneity across cells. Procedures like clustering and dimensionality reduction compare cells based on their gene expression profiles, which involves aggregating per-gene differences into a single (dis)similarity metric between a pair of cells. The choice of genes to use 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. This aims to preserve interesting biological structure without the variance that obscures that structure. It also reduces the size of the dataset to improve computational efficiency of later steps.

The simplest approach to feature selection is to select the most variable genes 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). Several methods are available to quantify the variation per gene and to select the highly variable genes (HVGs), which we will discuss below. For demonstration, we will use the 10X PBMC dataset:

### loading ###
library(BiocFileCache)
bfc <- BiocFileCache("raw_data", ask = FALSE)
raw.path <- bfcrpath(bfc, file.path("http://cf.10xgenomics.com/samples",
"cell-exp/2.1.0/pbmc4k/pbmc4k_raw_gene_bc_matrices.tar.gz"))
untar(raw.path, exdir=file.path(tempdir(), "pbmc4k"))

library(DropletUtils)
fname <- file.path(tempdir(), "pbmc4k/raw_gene_bc_matrices/GRCh38")

### gene-annotation ###
library(scater)
rownames(sce.pbmc) <- uniquifyFeatureNames(rowData(sce.pbmc)$ID, rowData(sce.pbmc)$Symbol)

library(EnsDb.Hsapiens.v86)
location <- mapIds(EnsDb.Hsapiens.v86, keys=rowData(sce.pbmc)$ID, column="SEQNAME", keytype="GENEID") ### cell-detection ### set.seed(100) e.out <- emptyDrops(counts(sce.pbmc)) sce.pbmc <- sce.pbmc[,which(e.out$FDR <= 0.001)]

### quality-control ###
sce.pbmc <- calculateQCMetrics(sce.pbmc, feature_controls=list(Mito=which(location=="MT")))
high.mito <- isOutlier(sce.pbmc$pct_counts_Mito, nmads=3, type="higher") sce.pbmc <- sce.pbmc[,!high.mito] ### normalization ### library(scran) set.seed(1000) clusters <- quickCluster(sce.pbmc) sce.pbmc <- computeSumFactors(sce.pbmc, min.mean=0.1, cluster=clusters) sce.pbmc <- normalize(sce.pbmc) sce.pbmc ## class: SingleCellExperiment ## dim: 33694 3922 ## metadata(1): log.exprs.offset ## assays(2): counts logcounts ## rownames(33694): RP11-34P13.3 FAM138A ... AC213203.1 FAM231B ## rowData names(10): ID Symbol ... total_counts log10_total_counts ## colnames(3922): AAACCTGAGAAGGCCT-1 AAACCTGAGACAGACC-1 ... ## TTTGTCACAGGTCCAC-1 TTTGTCATCCCAAGAT-1 ## colData names(38): Sample Barcode ... ## pct_counts_in_top_200_features_Mito ## pct_counts_in_top_500_features_Mito ## reducedDimNames(0): ## spikeNames(0): As well as the 416B dataset: ### loading ### library(scRNAseq) sce.416b <- LunSpikeInData(which="416b") ### gene-annotation ### library(AnnotationHub) ens.mm.v97 <- AnnotationHub()[["AH73905"]] rowData(sce.416b)$ENSEMBL <- rownames(sce.416b)
rowData(sce.416b)$SYMBOL <- mapIds(ens.mm.v97, keys=rownames(sce.416b), keytype="GENEID", column="SYMBOL") rowData(sce.416b)$SEQNAME <- mapIds(ens.mm.v97, keys=rownames(sce.416b),
keytype="GENEID", column="SEQNAME")

library(scater)
rownames(sce.416b) <- uniquifyFeatureNames(rowData(sce.416b)$ENSEMBL, rowData(sce.416b)$SYMBOL)

### quality-control ###
mito <- which(rowData(sce.416b)$SEQNAME=="MT") sce.416b <- calculateQCMetrics(sce.416b, feature_controls=list(Mt=mito)) combined <- paste0(sce.416b$block, ":", sce.416b$phenotype) libsize.drop <- isOutlier(sce.416b$total_counts, nmads=3, type="lower",
log=TRUE, batch=combined)
feature.drop <- isOutlier(sce.416b$total_features_by_counts, nmads=3, type="lower", log=TRUE, batch=combined) spike.drop <- isOutlier(sce.416b$pct_counts_ERCC, nmads=3, type="higher",
batch=combined)

keep <- !(libsize.drop | feature.drop | spike.drop)
sce.416b <- sce.416b[,keep]

### normalization ###
library(scran)
sce.416b <- computeSumFactors(sce.416b)
sce.416b <- computeSpikeFactors(sce.416b, general.use=FALSE)
sce.416b <- normalize(sce.416b)
sce.416b
## class: SingleCellExperiment
## dim: 46703 188
## assays(2): counts logcounts
## rownames(46703): ENSMUSG00000102693 ENSMUSG00000064842 ... SIRV7
##   CBFB-MYH11-mcherry
## rowData names(14): Length ENSEMBL ... total_counts
##   log10_total_counts
## colnames(188): SLX-9555.N701_S502.C89V9ANXX.s_1.r_1
##   SLX-9555.N701_S503.C89V9ANXX.s_1.r_1 ...
##   SLX-11312.N712_S507.H5H5YBBXX.s_8.r_1
##   SLX-11312.N712_S517.H5H5YBBXX.s_8.r_1
## colData names(63): Source Name cell line ...
##   pct_counts_in_top_200_features_SIRV
##   pct_counts_in_top_500_features_SIRV
## reducedDimNames(0):
## spikeNames(2): ERCC SIRV

## 10.2 Quantifying per-gene variation

### 10.2.1 Variance of the log-counts

The simplest approach to quantifying per-gene variation is to simply compute the variance of the log-normalized expression values (referred to as “log-counts” for simplicity) for each gene across all cells in the population. This has an advantage in that the feature selection is based on the same log-values that are used for later downstream steps. Genes with the largest variances in log-values will contribute the most to the Euclidean distances between cells. By using log-values here, we ensure that our quantitative definition of heterogeneity is consistent throughout the entire analysis.

Calculation of the per-gene variance is simple but feature selection requires modelling of the mean-variance relationship. As discussed briefly in the last chapter, the log-transformation does not achieve perfect variance stabilization. This means that the variance of a gene is more affected by its abundance than the underlying biological heterogeneity. To account for this effect, we use the trendVar() function to fit a trend to the variance with respect to abundance across all genes (Figure 10.1).

library(scran)

# No spike-ins, so setting 'use.spikes=FALSE'.
fit.pbmc <- trendVar(sce.pbmc, use.spikes=FALSE)

plot(fit.pbmc$mean, fit.pbmc$var, xlab="Mean of log-expression",
ylab="Variance of log-expression")
curve(fit.pbmc$trend(x), col="dodgerblue", add=TRUE, lwd=2) At any given abundance, we assume that the expression profiles of most genes are dominated by random technical noise (see 10.2.3 for details). Under this assumption, our trend represents an estimate of the technical noise as a function of abundance. We then use the decomposeVar() function to break down the total variance of each gene into the technical component, i.e., the fitted value of the trend at that gene’s abundance; and the biological component, defined as the difference between the total variance and the technical component. This biological component represents the “interesting” variation for each gene and can be used as the metric for HVG selection. dec.pbmc <- decomposeVar(fit=fit.pbmc) # Ordering by most interesting genes for inspection. dec.pbmc[order(dec.pbmc$bio, decreasing=TRUE),] 
## DataFrame with 33694 rows and 6 columns
##                     mean             total                bio
##                <numeric>         <numeric>          <numeric>
## LYZ      1.9806721473164  5.12675017749729   4.23343856586781
## S100A9  1.95576910151395  4.62353584269683   3.73059071698752
## S100A8  1.72447657712873  4.49362852270639   3.60800589885247
## HLA-DRA 2.10051754708271  3.73885530331428   2.84304344288055
## CD74    2.90746230772648   3.3292370802074   2.30788747256288
## ...                  ...               ...                ...
## MT-CO2  4.19188502683674 0.544465230103077 -0.635044758465113
## PTMA    3.84384732827719 0.498416795149973 -0.650614922570656
## HLA-B   4.51931577517642 0.499881887649168 -0.701096896119891
## TMSB4X   6.1055295386377 0.470802972021666 -0.776306567179052
## B2M     5.97628515673904  0.34502768951129 -0.899609719380582
##                      tech   p.value       FDR
##                 <numeric> <numeric> <numeric>
## LYZ     0.893311611629482         0         0
## S100A9  0.892945125709309         0         0
## S100A8  0.885622623853927         0         0
## HLA-DRA 0.895811860433723         0         0
## CD74     1.02134960764452         0         0
## ...                   ...       ...       ...
## MT-CO2   1.17950998856819         1         1
## PTMA     1.14903171772063         1         1
## HLA-B    1.20097878376906         1         1
## TMSB4X   1.24710953920072         1         1
## B2M      1.24463740889187         1         1

(Careful readers will notice that some genes have negative biological components. Negative components of variation have no obvious interpretation and can be ignored for most applications. They are inevitable when fitting a trend to the per-gene variances, as approximately half of the genes will lie below the trend.)

### 10.2.2 Coefficient of variation

An alternative approach to quantification uses the squared coefficient of variation (CV2) of the normalized expression values prior to log-transformation. The CV2 is a widely used metric for describing variation in non-negative data and is closely related to the dispersion parameter of the negative binomial distribution in packages like edgeR and DESeq2. We compute the CV2 for each gene in the PBMC dataset using the improvedCV2() function, which provides a more robust implementation of the approach described by Brennecke et al. (2013).

# TODO: upgrade scran to use the latest version of this.
# No spike-ins, so setting 'spike.type=NA'.
cv2.pbmc <- improvedCV2(sce.pbmc, spike.type=NA) 

This allows us to model the mean-variance relationship when considering the relevance of each gene (Figure 10.2). Again, our assumption is that most genes contain random noise and that the trend captures mostly technical variation. Large CV$$^2$$ values that deviate strongly from the trend are likely to represent genes affected by biological structure.

plot(cv2.pbmc$mean, cv2.pbmc$cv2, log="xy")
o <- order(cv2.pbmc$mean) lines(cv2.pbmc$mean[o], cv2.pbmc$trend[o], col="dodgerblue", lwd=2) For each gene, we quantify the deviation from the trend in terms of the ratio of its CV2 to the fitted value of trend at its abundance. This is more appropriate than the directly subtracting the trend from the CV2, as the magnitude of the ratio is not affected by the mean. diff <- cv2.pbmc$cv2/cv2.pbmc$trend cv2.pbmc[order(diff, decreasing=TRUE),] ## DataFrame with 33694 rows and 6 columns ## mean cv2 trend ## <numeric> <numeric> <numeric> ## PRTFDC1 0.349193768088249 3781.11597077882 3.7476790639203 ## GNG11 2.0189430911897 415.713443913274 0.870369442833297 ## RCAN1 0.194958683724188 3033.04562507858 6.50000581102497 ## SDPR 1.60637825330209 467.847955624938 1.02491132642899 ## BEND2 0.210970404091042 2644.94583077305 6.02707265545054 ## ... ... ... ... ## AC023491.2 0 NaN Inf ## AC233755.2 0 NaN Inf ## AC233755.1 0 NaN Inf ## AC213203.1 0 NaN Inf ## FAM231B 0 NaN Inf ## ratio p.value FDR ## <numeric> <numeric> <numeric> ## PRTFDC1 1008.92203048559 0 0 ## GNG11 477.628721155479 0 0 ## RCAN1 466.621986696395 0 0 ## SDPR 456.476520027366 0 0 ## BEND2 438.844192193554 0 0 ## ... ... ... ... ## AC023491.2 NaN 1 1 ## AC233755.2 NaN 1 1 ## AC233755.1 NaN 1 1 ## AC213203.1 NaN 1 1 ## FAM231B NaN 1 1 Both the CV2 and the variance of log-counts are effective metrics for quantifying variation in gene expression. The CV2 tends to give higher rank to low-abundance HVGs driven by upregulation in rare subpopulations, for which the increase in variance on the raw scale is stronger than that on the log-scale. However, the variation described by the CV2 is less directly relevant to downstream procedures operating on the log-counts, and the reliance on the ratio can assign high rank to uninteresting genes with low absolute variance. We generally prefer the use of the variance of log-counts and will use it in the following sections, though the many of the same principles apply to procedures based on the CV2. ### 10.2.3 Quantifying technical noise The use of a trend fitted to endogenous genes assumes that the expression profiles of most genes are dominated by random technical noise. In practice, all expressed genes will exhibit some non-zero level of biological variability due to events like transcriptional bursting. This suggests that our estimates of the technical component are likely to be inflated. It would be more appropriate to consider these estimates as technical noise plus “uninteresting” biological variation, under the assumption that most genes are unaffected by the relevant heterogeneity in the population. This assumption is generally reasonable but may be problematic in some scenarios where many genes at a particular abundance are affected by a biological process. For example, strong upregulation of cell type-specific genes may result in an enrichment of HVGs at high abundances. This would inflate the fitted trend at high abundances and compromise the detection of the affected genes. We can avoid this problem by fitting a mean-dependent trend to the variance of the spike-in transcripts, if they are available (Figure 10.3). The premise here is that spike-ins should not be affected by biological variation, so the fitted value of the spike-in trend should represent a better estimate of the technical component for each gene. fit.416b <- trendVar(sce.416b) dec.416b <- decomposeVar(sce.416b, fit.416b) dec.416b[order(dec.416b$bio, decreasing=TRUE),]
## DataFrame with 46703 rows and 6 columns
##                                mean             total               bio
##                           <numeric>         <numeric>         <numeric>
## Lyz2               6.57237822734111  13.9096268809543  12.0532209419185
## Ccl9               6.72849586510843  13.2061002387558   11.530457079581
## Top2a              5.77874153915785  13.9853581173099  10.9827877439828
## Cd200r3            4.84945724546366  15.5251666159147  10.9183993789504
## ENSMUSG00000096842 4.21058861964987  15.9330829391269  10.5640743777755
## ...                             ...               ...               ...
## ENSMUSG00000097554 3.21315842136385 0.869510710256039 -5.25440177296372
## ENSMUSG00000083011 2.91045997741136 0.828657162029132 -5.29256475478502
## ENSMUSG00000082064 3.60090056018739 0.607899262223165 -5.32548110628434
## ENSMUSG00000083007 2.71208464889622 0.762293733047904  -5.3432512239974
## ENSMUSG00000078087 3.14961465503528 0.738996070206387 -5.39306841637665
##                                tech               p.value
##                           <numeric>             <numeric>
## Lyz2               1.85640593903582 8.54995227701758e-185
## Ccl9               1.67564315917488 1.53437515087334e-198
## Top2a               3.0025703733271  9.88028287012378e-89
## Cd200r3            4.60676723696428  2.14422927380176e-49
## ENSMUSG00000096842 5.36900856135144  3.89126290980577e-38
## ...                             ...                   ...
## ENSMUSG00000097554 6.12391248321976                     1
## ENSMUSG00000083011 6.12122191681415                     1
## ENSMUSG00000082064  5.9333803685075                     1
## ENSMUSG00000083007 6.10554495704531                     1
## ENSMUSG00000078087 6.13206448658304                     1
##                                      FDR
##                                <numeric>
## Lyz2               2.65641317278752e-181
## Ccl9               7.94533550347791e-195
## Top2a               5.90334234460575e-86
## Cd200r3             3.85828807244237e-47
## ENSMUSG00000096842  3.98567948678215e-36
## ...                                  ...
## ENSMUSG00000097554                     1
## ENSMUSG00000083011                     1
## ENSMUSG00000082064                     1
## ENSMUSG00000083007                     1
## ENSMUSG00000078087                     1
plot(dec.416b$mean, dec.416b$total, xlab="Mean of log-expression",
ylab="Variance of log-expression")
is.spike <- isSpike(sce.416b)
points(dec.416b$mean[is.spike], dec.416b$total[is.spike], col="red", pch=16)
curve(fit.416b$trend(x), col="dodgerblue", add=TRUE, lwd=2) In the absence of spike-in data, one can attempt to create a trend by making some distributional assumptions about the noise. For example, UMI counts typically exhibit near-Poisson noise when only technical effects are considered. This can be used to construct a mean-variance trend in the log-counts (Figure 10.4) with the makeTechTrend() function. Note the increased residuals of the high-abundance genes, which can be interpreted as the amount of biological variation that was assumed to be “uninteresting” when fitting the gene-based trend in Figure 10.1. tech.trend <- makeTechTrend(x=sce.pbmc) fit.pbmc2 <- fit.pbmc fit.pbmc2$trend <- tech.trend # overwrite trend.
dec.pbmc2 <- decomposeVar(fit=fit.pbmc2)
dec.pbmc2 <- dec.pbmc2[order(dec.pbmc2$bio, decreasing=TRUE),] head(dec.pbmc2) ## DataFrame with 6 rows and 6 columns ## mean total bio ## <numeric> <numeric> <numeric> ## LYZ 1.9806721473164 5.12675017749729 4.4957050381929 ## S100A9 1.95576910151395 4.62353584269683 3.98858788456279 ## S100A8 1.72447657712873 4.49362852270639 3.82909592851813 ## HLA-DRA 2.10051754708271 3.73885530331428 3.12807677300474 ## CD74 2.90746230772648 3.3292370802074 2.87761651579676 ## CST3 1.49489535215485 2.97218277928419 2.29406052725448 ## tech p.value FDR ## <numeric> <numeric> <numeric> ## LYZ 0.631045139304399 0 0 ## S100A9 0.634947958134047 0 0 ## S100A8 0.664532594188262 0 0 ## HLA-DRA 0.610778530309538 0 0 ## CD74 0.451620564410635 0 0 ## CST3 0.678122252029714 0 0 plot(dec.pbmc2$mean, dec.pbmc2$total, pch=16, xlab="Mean of log-expression", ylab="Variance of log-expression") curve(fit.pbmc2$trend(x), col="dodgerblue", add=TRUE)

### 10.2.4 Accounting for blocking factors

#### 10.2.4.2 Using a design matrix

The use of block-specific trends is the recommended approach for experiments with a single blocking factor. However, this is not practical for studies involving a large number of blocking factors and/or covariates. In such cases, we can use the design= argument to specify a design matrix with uninteresting factors of variation. decomposeVar() will then focus on genes with large residual variances, i.e., beyond that explained by the design matrix. We illustrate again with the 416B data set, blocking on the plate of origin and oncogene induction.

design <- model.matrix(~factor(block) + phenotype, colData(sce.416b))
fit.416b.2 <- trendVar(sce.416b, design=design)
dec.416b.2 <- decomposeVar(sce.416b, fit.416b.2)
dec.416b.2[order(dec.416b.2$bio, decreasing=TRUE),] ## DataFrame with 46703 rows and 6 columns ## mean total bio ## <numeric> <numeric> <numeric> ## Ccnb2 5.89821323258258 9.79093386488465 7.0893697170353 ## Lyz2 6.57237822734111 8.84765039552134 7.07843003454142 ## Gem 5.86432958703213 9.6816178335956 6.92675481460231 ## ENSMUSG00000076617 6.37255828888088 7.87882496529147 5.86185528899622 ## Idh1 5.98634090489272 8.40011571572379 5.83482283821305 ## ... ... ... ... ## Gm7429 3.4731788564726 0.253468796782048 -5.74461281024374 ## ENSMUSG00000083061 3.57276296115205 0.169207180855588 -5.75979426776692 ## ENSMUSG00000081740 2.54482110453333 0.271685903926373 -5.81003429108432 ## ENSMUSG00000103903 3.0916185774635 0.198365087332871 -5.92335585818509 ## ENSMUSG00000095762 2.84338170706645 0.205600789787391 -5.93183550024739 ## tech p.value ## <numeric> <numeric> ## Ccnb2 2.70156414784935 3.19116390139292e-56 ## Lyz2 1.76922036097993 8.9427450119608e-99 ## Gem 2.75486301899329 4.98139083384899e-53 ## ENSMUSG00000076617 2.01696967629525 1.37607077983069e-64 ## Idh1 2.56529287751074 3.43618243897962e-46 ## ... ... ... ## Gm7429 5.99808160702579 1 ## ENSMUSG00000083061 5.92900144862251 1 ## ENSMUSG00000081740 6.08172019501069 1 ## ENSMUSG00000103903 6.12172094551797 1 ## ENSMUSG00000095762 6.13743629003478 1 ## FDR ## <numeric> ## Ccnb2 1.32786609339746e-53 ## Lyz2 1.81203342842357e-95 ## Gem 1.7196499142274e-50 ## ENSMUSG00000076617 8.43821087147755e-62 ## Idh1 8.08787102960638e-44 ## ... ... ## Gm7429 1 ## ENSMUSG00000083061 1 ## ENSMUSG00000081740 1 ## ENSMUSG00000103903 1 ## ENSMUSG00000095762 1 This strategy is simple but somewhat inaccurate as it does not consider the mean expression in each blocking level. Briefly, the technical component is estimated as the fitted value of the trend at the average abundance for each gene. However, the true technical component is the average of the fitted values at the per-block means, which may be quite different for strong batch effects and non-linear mean-variance relationships. The multiBlockVar() approach is safer and should be preferred in all situations where it is applicable. ## 10.3 Selecting highly variable genes ### 10.3.1 Based on the largest metrics Once we have quantified the per-gene variation, the next step is to select the subset of HVGs to use in downstream analyses. The simplest approach is to simply take the top $$X$$ genes with the largest values for the relevant measure of variation. For decomposeVar(), this would be the genes with the largest biological components: hvg.pbmc.var <- head(order(dec.pbmc$bio, decreasing=TRUE), 1000)
hvg.pbmc.var <- rownames(dec.pbmc)[hvg.pbmc.var]

For improvedCV2(), this would instead be the genes with the largest ratios:

hvg.pbmc.cv2 <- head(order(diff, decreasing=TRUE), 1000)
hvg.pbmc.cv2 <- rownames(dec.pbmc)[hvg.pbmc.cv2]

The main advantage of this approach is that the user can easily control the number of genes retained. This ensures that the computational complexity of downstream calculations is easily predicted. It is also fairly easy to translate the choice of $$X$$ into a biological statement. Recall our trend-fitting assumption that most genes are not differentially expressed between cell types or states in our population. If we quantify this assumption into a statement that, e.g., no more than 5% of genes are differentially expressed, we can simply take the top 5% of genes with the largest variance.

The main disadvantage of this approach that it turns HVG selection into a competition between genes, whereby a subset of very highly variable genes can push other informative genes out of the top set. This can be problematic if a single subpopulation is very different from the others. In such cases, the top set will be dominated by differentially expressed genes involving the outlier subpopulation, compromising resolution of heterogeneity between the other populations. Similar problems are encountered when the magnitude of the chosen variance measure varies with abundance.

The choice of $$X$$ is also fairly arbitrary, with any value from 500 to 5000 being considered “reasonable”. We have chosen $$X=1000$$ in the code above though there is no particular a priori reason for doing so. A larger $$X$$ will reduce the risk of discarding signal at the cost of increasing noise that obscures signal. Our recommendation is to simply pick an arbitrary $$X$$ and proceed with the rest of the analysis, with the intention of testing other choices later, rather than spending much time worrying about obtaining the “optimal” value9.

### 10.3.2 Based on a fixed threshold

Another approach to feature selection is to set a fixed threshold of one of the metrics. This is most commonly done with the (adjusted) $$p$$-value reported by each of the above methods. The $$p$$-value for each gene is generated by testing against the null hypothesis that the variance is equal to the trend. For example, we might define our HVGs as all genes that have adjusted $$p$$-values below 0.05.

hvg.pbmc.var.2 <- rownames(dec.pbmc)[dec.pbmc$FDR <= 0.05] length(hvg.pbmc.var.2) ## [1] 4282 This approach is simple to implement and - if the test holds its size - it controls the false discovery rate (FDR). That is, it returns a subset of genes where the proportion of false positives is expected to be below the specified threshold. This can occasionally be useful in applications where the HVGs themselves are of interest. For example, if a collaborator were to take a list of HVGs back to the bench to verify the existence of heterogeneous expression for some of the genes, we would want to control the FDR in that list. The downside of this approach is that it is less predictable than the top $$X$$ strategy. The number of genes returned depends on the type II error rate of the test and the severity of the multiple testing correction. One might obtain no genes or every gene at a given FDR threshold, depending on the circumstances. Moreover, control of the FDR is usually not helpful at this stage of the analysis. We are not interpreting the individual HVGs themselves but are only using them for feature selection prior to downstream steps. There is no reason to think that a 5% threshold on the FDR yields a more suitable compromise between bias and noise. ### 10.3.3 Keeping all genes above the trend As the title suggests, this involves keeping all genes above the trend. The rationale is to only remove the obviously uninteresting genes with variances below the trend. By doing so, we avoid the need to make any judgement calls regarding what level of variation is interesting enough to retain. This approach represents one extreme of the bias-variance trade-off where bias is minimized at the cost of maximizing noise. For decomposeVar(), it equates to keeping all positive biological components: hvg.pbmc.var.3 <- rownames(dec.pbmc)[dec.pbmc$bio > 0]
length(hvg.pbmc.var.3)
## [1] 8110

For improvedCV2(), this involves keeping all ratios above 1:

hvg.pbmc.cv2.3 <- rownames(cv2.pbmc)[diff > 0]
length(hvg.pbmc.cv2.3)
## [1] 33694

This strategy is the most conservative as it does not discard any potential biological signal. Weak or secondary population structure is given the chance to manifest as the affected genes are retained. This makes it useful for reliable automated processing of diverse data sets where the primary factor of variation in one data set is a secondary factor in another data set (and thus overlooked by the top $$X$$ approach). The obvious cost is that more noise is also captured, which can reduce the resolution of otherwise well-separated populations. From a practical perspective, the use of more genes involves more computational work in each downstream step.

## 10.4 Putting it all together

The few lines of code below will select the top 10% of genes with the highest biological components.

fit.pbmc <- trendVar(sce.pbmc, use.spikes=FALSE)
dec.pbmc <- decomposeVar(fit=fit.pbmc)
chosen <- rownames(dec.pbmc)[order(dec.pbmc\$bio, decreasing=TRUE)]
chosen <- head(chosen, nrow(dec.pbmc) * 0.1)
length(chosen)
## [1] 3369

We can then subset the SingleCellExperiment to only retain our selection of HVGs. This ensures that downstream methods will only use these genes for their calculations.

sce.pbmc <- sce.pbmc[chosen,]
dim(sce.pbmc)
## [1] 3369 3922

Alternatively, some methods may allow users to pass in the full SingleCellExperiment object and specify the genes to use via an extra argument like subset.row=. This may be more convenient in the context of the overall analysis, where genes outside of this subset may still be of interest during DE analyses or for visualization.

## 10.5 Session Info

R version 3.6.0 (2019-04-26)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 14.04.6 LTS

Matrix products: default
BLAS/LAPACK: /app/easybuild/software/OpenBLAS/0.2.18-GCC-5.4.0-2.26-LAPACK-3.6.1/lib/libopenblas_prescottp-r0.2.18.so

locale:
[1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C
[3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8
[5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8
[7] LC_PAPER=en_US.UTF-8       LC_NAME=C
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C

attached base packages:
[1] stats4    parallel  stats     graphics  grDevices utils     datasets
[8] methods   base

other attached packages:
[1] scran_1.13.9                SingleCellExperiment_1.7.0
[3] SummarizedExperiment_1.15.5 DelayedArray_0.11.4
[5] BiocParallel_1.19.0         matrixStats_0.54.0
[7] Biobase_2.45.0              GenomicRanges_1.37.14
[9] GenomeInfoDb_1.21.1         IRanges_2.19.10
[11] S4Vectors_0.23.17           BiocGenerics_0.31.5
[13] BiocStyle_2.13.2            Cairo_1.5-10

loaded via a namespace (and not attached):
[1] Rcpp_1.0.2               rsvd_1.0.2
[3] locfit_1.5-9.1           lattice_0.20-38
[5] assertthat_0.2.1         digest_0.6.20
[7] R6_2.4.0                 dynamicTreeCut_1.63-1
[9] evaluate_0.14            ggplot2_3.2.0
[11] pillar_1.4.2             zlibbioc_1.31.0
[13] rlang_0.4.0              lazyeval_0.2.2
[15] irlba_2.3.3              Matrix_1.2-17
[17] rmarkdown_1.14           BiocNeighbors_1.3.3
[19] statmod_1.4.32           stringr_1.4.0
[21] igraph_1.2.4.1           RCurl_1.95-4.12
[23] munsell_0.5.0            vipor_0.4.5
[25] compiler_3.6.0           BiocSingular_1.1.5
[27] xfun_0.8                 pkgconfig_2.0.2
[29] ggbeeswarm_0.6.0         htmltools_0.3.6
[31] tidyselect_0.2.5         gridExtra_2.3
[33] tibble_2.1.3             GenomeInfoDbData_1.2.1
[35] bookdown_0.12            edgeR_3.27.9
[37] viridisLite_0.3.0        crayon_1.3.4
[39] dplyr_0.8.3              bitops_1.0-6
[41] grid_3.6.0               gtable_0.3.0
[43] magrittr_1.5             scales_1.0.0
[45] dqrng_0.2.1              stringi_1.4.3
[47] XVector_0.25.0           viridis_0.5.1
[49] limma_3.41.15            scater_1.13.9
[51] DelayedMatrixStats_1.7.1 tools_3.6.0
[53] beeswarm_0.2.3           glue_1.3.1
[55] purrr_0.3.2              yaml_2.2.0
[57] colorspace_1.4-1         BiocManager_1.30.4
[59] knitr_1.23              

### Bibliography

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.

1. When everything’s said and done, the “optimal” value is one that gets you a figure for your manuscript.