Why HVG Methods Disagree
Highly variable gene (HVG) selection is not a direct measurement of a single, method-independent property. A dispersion-based method compares observed variability with a mean-variance trend; a negative-binomial (NB) deviance method scores how poorly a fitted count model represents a gene; and a residual-based method summarizes variability after scaling counts relative to a fitted expectation. Because these approaches define and adjust for expected variation differently, they need not rank genes alike.
Disagreement is useful diagnostic information, not by itself evidence that one method is wrong. It can indicate that a gene lies in a region where one method is especially sensitive, that the fitted model or preprocessing differs, or that a few cells or a technical pattern drive the score. The goal is to explain the difference in context, not to choose a winner by majority vote.
A comparison is interpretable only when methods use the same cells, the same eligible genes, and appropriately comparable input data. Earlier material on filtering low-expression genes and choosing an expression scale is relevant here: if one method ranks a different set of genes, or uses counts while another uses transformed values, apparent disagreement can be caused by the comparison setup rather than by the ranking criterion itself.
What the Three Rankings Emphasize
Dispersion-based approaches start with variability relative to a gene’s mean. A raw variance-to-mean ratio is not generally comparable across the full expression range, because count variance changes with mean. Binned or regularized approaches address this by comparing a gene with genes of similar mean or with a smooth mean-variance trend. Their rankings therefore depend on the estimated trend and on how scores are standardized or regularized.
An NB deviance score measures lack of fit to a specified NB model. For a gene, let \(Y_i\) be its count in cell \(i\), \(\mu_i\) its fitted mean (including any modeled offset), and \(\theta\) its NB size parameter, with variance \(\mu_i+\mu_i^2/\theta\). The contribution to the NB deviance is
where terms involving \(Y_i\log(Y_i/\mu_i)\) are defined as zero when \(Y_i=0\). Summing contributions gives a measure of discrepancy between observed counts and fitted means. Larger deviance can reflect cell-to-cell structure not captured by the model, but it can also reflect outliers, misspecification, or technical effects. The ranking is meaningful only with a stated model and comparable score definition across genes.
A Pearson residual for a cell-gene pair is the observed-minus-fitted count divided by its model-based standard deviation:
A residual-based ranking often summarizes the spread of these residuals across cells, such as their variance. Under a correctly specified model, residuals should not show systematic structure left unexplained by the fitted mean and variance. However, a high residual spread is not proof of biological variation: model misspecification, extreme observations, or unmodeled technical covariates can also produce it.
These interpretations rely on comparable model fitting and preprocessing. For NB scores and residuals, the count model should be appropriate enough for the data, the offset should represent the intended exposure or library-size adjustment, and fitted parameters should not be dominated by a few observations. For dispersion scores, the mean-variance adjustment should be estimated from a sufficiently informative set of genes. None of the rankings alone tests whether a gene is biologically relevant.
Comparing Ranks and Locating Disagreement
First compare ranks on the shared gene universe, rather than comparing raw scores with different scales. A rank correlation summarizes agreement in ordering; a list overlap summarizes agreement near a chosen cutoff. Neither is sufficient alone. Two methods can have moderate overall correlation but disagree substantially among the top genes, or have a high correlation while differing at a cutoff that matters for the analysis.
For rankings without ties, Spearman’s rank correlation can be computed from rank differences. If there are \(m\) genes and \(d_g\) is the difference between the two ranks for gene \(g\), then
This expression assumes no tied ranks. With ties, compute the ordinary correlation of the rank values using a consistent tie-handling rule. The coefficient ranges from -1 to 1: 1 means identical rank order, 0 means no linear association between the rank values, and -1 means reversed order. It does not indicate why rankings differ.
Worked Example: Quantifying Rank Disagreement
Suppose four hypothetical genes have these ranks in two methods, with rank 1 highest:
| Gene | Dispersion rank | Residual rank | Rank difference \(d_g\) | \(d_g^2\) |
|---|---|---|---|---|
| A | 1 | 2 | -1 | 1 |
| B | 2 | 4 | -2 | 4 |
| C | 3 | 1 | 2 | 4 |
| D | 4 | 3 | 1 | 1 |
There are \(m=4\) genes and \(\sum d_g^2=1+4+4+1=10\). Therefore,
The rank correlation is zero for these four rank pairs. In particular, gene C is ranked third by dispersion but first by residual variability, while gene B moves from second to fourth. This small example does not establish that either method is generally unreliable; it identifies genes to investigate.
Next locate discordant genes on a mean-variance plot. Plot each gene’s mean expression against an appropriate variability summary, then distinguish the observed quantity from the method’s expected trend or model-based reference. Label or highlight genes with large rank changes. Ask whether the changes cluster at low means, high means, a particular range of fitted values, or among genes with a few extreme counts. Such patterns can point to different score behavior or a failure of assumptions.
Use the same cells, gene universe, count or expression scale, and documented filtering rules. Verify that each method’s score has the same direction, so that a larger value consistently means a higher rank.
Calculate a rank correlation and inspect top-list overlap at the planned selection size. Identify individual genes whose ranks shift sharply.
Plot discordant genes against mean expression and the relevant variability or model-based score. Check whether disagreements concentrate in a range where the trend or fit is weak.
Inspect expression across cells and technical groups, including library size and batch, as well as plausible cell populations. Look for outliers, zero-heavy patterns, and variability that is confined to a small group.
Describe which genes and patterns drive disagreement. Treat agreement as a robustness check, not as proof that the selected genes are biologically meaningful.
Worked Diagnostics
Worked Example: A Mean-Variance Pattern Behind a Rank Shift
Consider hypothetical gene G. Across 100 cells its mean normalized count is 2.0 and its observed variance is 8.0. Its variance-to-mean ratio is \(8.0/2.0=4.0\). A fitted mean-variance trend predicts variance 3.0 for genes at this mean, so the variance excess is \(8.0-3.0=5.0\). Suppose its dispersion-based rank is 8, while its NB deviance rank is 76.
The ranks can differ because the scores are not interchangeable. The dispersion score emphasizes variability relative to the trend at the gene’s mean; the NB deviance score reflects discrepancy from a particular fitted count model, including its mean structure and variance assumptions. The numbers alone do not tell which explanation applies.
Suppose inspection shows that 90 cells have counts near 1 and 10 cells have counts near 11. The mean is \((90(1)+10(11))/100=2.0\), as stated, and the sample variance is approximately 9.09 (the sum of squared deviations is \(90(1-2)^2+10(11-2)^2=900\), divided by 99). That is not the stated variance of 8.0, so this proposed cell-level pattern is inconsistent with the summary and cannot be used to explain it. A valid diagnosis must use the actual cell-level counts or a consistent summary. This check prevents a plausible-sounding but numerically impossible story.
If verified data instead showed a small high-count subgroup, one would ask whether it corresponds to a coherent cell type, a batch, or a technical artifact. A mean-variance plot points to where a discrepancy occurs; examining cells and metadata is needed to interpret it.
Worked Example: Distinguishing Technical Structure From a Candidate Cell-State Signal
In hypothetical data, gene H is ranked 12th by NB deviance, 35th by regularized dispersion, and 18th by residual variability. Its residuals are mostly near zero, but 15 of 120 cells have large positive residuals. Those 15 cells all belong to one annotated cell population and are distributed across three sequencing batches. The calculation \(15/120=0.125\) shows that 12.5% of cells carry most of the positive residual pattern.
The concentration in one population and presence across multiple batches make a cell-state-associated pattern plausible, but do not prove a biological cause. Check whether those cells also differ in library size, whether the pattern remains after the modeled normalization or offsets, and whether the gene’s counts are supported by more than a handful of extreme observations. If instead all 15 cells came from one batch, technical confounding would be a more immediate concern.
The rank ordering is consistent with a gene whose cell-specific departures are prominent under the NB model, while its mean-adjusted dispersion is less exceptional. The appropriate conclusion is that H merits contextual inspection, not that NB deviance has identified a causal marker.
A practical alternative to selecting genes by only one ranking is to compare rankings and investigate their disagreement before choosing a feature set. A consensus approach can favor genes ranked highly by several methods, but it may suppress real signals that only one score captures. Conversely, choosing the union of top lists can increase sensitivity while admitting genes driven by method-specific artifacts. These are trade-offs, not universally superior solutions; the choice should match the analysis goal and be reported transparently.
Failure Modes and Interpretation
Common diagnostics include plotting counts or residuals by library size and batch, checking whether top genes are dominated by a few cells, and comparing their expression across plausible cell groups. A suspicious mean-variance trend, unstable fitted NB parameters, extreme residuals, or a rank shift concentrated among barely expressed genes warrants closer review. Such checks identify possible problems; they do not automatically establish their cause.
A particularly important distinction is between technical structure and biological heterogeneity. A gene varying across cells can reflect a genuine cell population, but can also reflect sequencing depth, batch, ambient RNA, or model misspecification. A pattern that tracks a technical variable is a warning, not conclusive proof of artifact; a pattern aligned with an annotation is evidence to investigate, not proof of function. Expression and metadata should be considered together.
- Comparing raw scores across methods. Scores can have different scales and meanings. Compare ranks or carefully defined standardized scores instead.
- Using rank correlation as a complete diagnosis. A single summary hides which genes disagree and whether disagreement is concentrated near the selection cutoff.
- Calling a high deviance “biological variability.” Deviance measures mismatch to a fitted model; technical effects and misspecification can also increase it.
- Ignoring the comparison universe. Different filters or expression scales can create apparent rank changes. Restrict comparison to a common, documented set.
- Assuming consensus guarantees validity. Multiple methods can share the same problematic input or miss the same rare-cell signal.
Use disagreement as a prompt for targeted checks: verify data and model alignment, identify where in the mean-variance range ranks diverge, and inspect the cells and technical covariates driving the score. The ultimate limitation remains that an unsupervised variability ranking cannot determine biological relevance on its own.
Check Your Understanding
Answer each question, then compare with the brief solutions.
- Two methods rank 200 genes on the same scale, but their top-30 overlap is small. Name two additional checks that help diagnose this disagreement.
- A gene has high NB deviance and its residuals are extreme only in cells from one sequencing batch. What can and cannot be concluded?
- For five untied genes, two rank lists have squared rank differences summing to 10. Calculate Spearman’s rank correlation using the no-ties formula.
- Why should a rank correlation not be interpreted as proof that two methods select the same useful genes?
Solutions. (1) Calculate overall rank correlation and inspect genes with large rank shifts; plot discordant genes against mean expression and check cell-level patterns by batch, library size, and cell population. (2) The score indicates mismatch to the fitted model, and its batch concentration raises a technical concern; it does not prove that the gene is nonbiological or establish a causal batch effect. (3) With \(m=5\), \(\rho_s=1-6(10)/(5(25-1))=1-60/120=0.500\). (4) Correlation summarizes ordering, can obscure disagreement near a selection cutoff, and says nothing by itself about biological interpretation or technical artifacts.