A Workflow Is a Set of Decisions, Not a Single Ranking
Selecting informative genes requires more than sorting a variability score. A defensible workflow specifies which genes are eligible, what expression scale is analyzed, how expected variability depends on mean expression, how genes are ranked, and whether the resulting choices are reproducible under resampling. These choices matter because a gene can appear variable simply because its mean is high, or can rank highly in the full data while depending on only a few cells.
The example below is hypothetical and deliberately small enough to calculate by hand. Its purpose is to show how the steps fit together, not to recommend its particular thresholds or number of selected genes for a real study. The earlier material on filtering low-expression genes, mean-variance modeling, ranking methods, and subsampling stability supplies the ideas used here; the emphasis is on applying them coherently.
Worked Example: Filter a Small Count Matrix
Suppose eight hypothetical cells have equal library totals of 100 counts. The table shows five candidate genes and an aggregate “Other genes” remainder that makes each total 100. The remainder is not a single gene and is not ranked. Equal library totals mean that no library-size adjustment is needed for this illustration; unequal totals would require a suitable normalization step before comparing expression variability.
| Cell | Gene A | Gene B | Gene C | Gene D | Gene E | Other genes | Total |
|---|---|---|---|---|---|---|---|
| 1 | 10 | 2 | 0 | 1 | 0 | 87 | 100 |
| 2 | 10 | 3 | 0 | 1 | 2 | 84 | 100 |
| 3 | 10 | 4 | 0 | 1 | 0 | 85 | 100 |
| 4 | 10 | 5 | 0 | 1 | 2 | 82 | 100 |
| 5 | 10 | 6 | 0 | 1 | 0 | 83 | 100 |
| 6 | 10 | 7 | 0 | 1 | 2 | 80 | 100 |
| 7 | 10 | 8 | 12 | 1 | 0 | 69 | 100 |
| 8 | 10 | 9 | 12 | 1 | 2 | 66 | 100 |
For this example, a gene is eligible if its count is at least 2 in at least two cells and its total count is at least 6. This is a deliberately explicit filter, not a universal standard. Gene A passes: it has count 10 in all eight cells and a total of 80. Gene B passes: it has at least 2 in all eight cells and a total of 44. Gene C passes in the full data: it has count 12 in two cells and a total of 24. Gene E passes: it has count 2 in four cells and a total of 8. Gene D fails because it never reaches count 2, even though it is detected at count 1 in every cell. Thus the full-data eligible set is A, B, C, and E.
Model Variability Relative to Mean Expression
Raw variance is not a fair ranking by itself when expected variability changes with mean expression. Here we use each gene’s sample mean \(\bar{x}\) and sample variance \(s^2\) across the eight cells, and compare the observed variance with a specified mean-variance trend. In a real analysis, the trend would generally be estimated from a suitable collection of genes using the chosen expression scale and method. For transparent arithmetic in this toy example, take the illustrative trend to be \(\hat{v}(\bar{x})=\bar{x}+0.25\bar{x}^2\).
The denominator puts the difference on a common scale for this defined score. A positive value means the observed variance exceeds the trend prediction; a negative value means it falls below. The score does not establish that a gene is biologically important, and its meaning depends on the trend and expression scale.
Worked Example: Calculate and Rank Mean-Adjusted Variability
For Gene B, the eight counts are 2, 3, 4, 5, 6, 7, 8, and 9. Their mean is \(\bar{x}=44/8=5.5\). The sum of squared deviations from 5.5 is 42, so the sample variance is \(s^2=42/(8-1)=6\). The trend prediction and score are
Gene B therefore falls below the illustrative trend. For Gene C, six zeroes and two values of 12 give \(\bar{x}=24/8=3\). The sum of squared deviations is \(6(0-3)^2+2(12-3)^2=216\), hence \(s^2=216/7\approx30.857\). Its trend prediction is \(3+0.25(3)^2=5.25\), and \(R_C=(30.857-5.25)/\sqrt{5.25}\approx11.18\).
The same calculation for Gene A gives \(\bar{x}=10\), \(s^2=0\), \(\hat{v}(10)=35\), and \(R_A=-35/\sqrt{35}\approx-5.916\). For Gene E, four zeroes and four twos give \(\bar{x}=1\), \(s^2=8/7\approx1.143\), \(\hat{v}(1)=1.25\), and \(R_E=(1.143-1.25)/\sqrt{1.25}\approx-0.096\). The ranking is therefore C, E, B, A. If the planned feature count is three, the full-data candidate set is C, E, and B.
Why the Ranking Behaves This Way
This example also gives a short justification for mean adjustment. If genes with different means are ranked only by observed variance, larger expected variance at larger means can dominate the list. The trend supplies a gene-specific reference level. Subtracting that level, \(s^2-\hat{v}(\bar{x})\), asks whether a gene varies more than expected for its mean; dividing by \(\sqrt{\hat{v}(\bar{x})}\) scales that excess. Consequently, a large score indicates excess variability relative to this particular trend, not large variability on every possible scale.
The calculation assumes that the cells being summarized are the intended analysis population, that the expression values are comparable across cells, and that the trend is an adequate description of expected variability over the relevant mean range. The toy data have equal library totals, which supports direct comparison of their counts. In real data, library-size differences, the chosen transformation or count model, and the trend fit must be addressed before interpreting the ranking. The score itself makes no independence claim, supplies no uncertainty interval, and does not adjust for batch or biological covariates.
Check Stability Without Freezing Eligibility
A full-data ranking can be sensitive to a small number of cells. To assess that sensitivity, draw a specified number of cells repeatedly, reapply the expression filter, recompute the ranking, and record selected genes. Eligibility must be recomputed in each draw: a gene that passes in the complete data may fail when some cells are omitted. Conversely, the selection frequencies have a denominator of all draws, while the number of draws in which a gene was eligible can be reported separately.
Worked Example: Subsample and Report Selection Stability
Continue the hypothetical workflow with 100 subsamples, each containing 6 of the 8 cells drawn without replacement. Six is an integer subsample size, or 75% of the cells. In each subsample, apply the same filter and ranking rule, then record the top three eligible genes. At least A, B, and E are eligible in every draw: A and B always have counts meeting the rule, while any six of the eight cells contain at least two of E’s four positive cells. Gene C is eligible only when both cells with count 12 are included.
The probability that a random six-cell subsample contains both C-positive cells is the number of ways to choose the remaining four cells from the other six, divided by the number of six-cell subsets:
Suppose the 100-draw hypothetical output has C eligible in 54 draws. The top-three selections are C, B, and E in those 54 draws, and A, B, and E in the other 46. Thus B and E are each selected 100 times, C is selected 54 times, and A is selected 46 times. There are exactly \(100\times3=300\) recorded selections, and the gene-specific counts sum to \(100+100+54+46=300\). This total reflects three recorded selections per draw; it does not mean that the same genes were eligible in every draw.
A stability summary could report the selection frequency as \(100/100=100\%\) for B and E, \(54/100=54\%\) for C, and \(46/100=46\%\) for A. For C, it is also useful to report eligibility in \(54/100\) draws and selection in \(54/54\) eligible draws. That conditional summary answers a different question from the overall selection frequency: C ranks highly whenever it can be considered, but it is not available in every subsample.
For this hypothetical analysis, suppose the selection plan was to take the top three genes in the full data and retain those candidates if their overall selection frequency is at least 50%. Then C, E, and B all meet that rule: their frequencies are 54%, 100%, and 100%. A is not added merely because it appears in 46% of draws; the feature-count rule restricted the candidate list to the top three full-data genes. The final set is therefore C, E, and B, with the explicit qualification that C’s stability reflects both its high rank when eligible and its limited eligibility under subsampling.
Diagnostics, Alternatives, and Failure Modes
A workflow is only useful if its choices can be checked. Before accepting a selected set, inspect the following:
- Filter behavior: report how many genes pass in the full data and how often candidate genes pass in subsamples. A feature that frequently disappears at filtering should not be described as consistently selected.
- Trend fit: inspect whether the trend follows the broad mean-variance pattern, especially near the mean range of selected genes. A poor fit can make scores reflect trend error rather than unusual variation.
- Ranking sensitivity: compare the selected set under reasonable alternative trend or expression-scale choices. Large changes are a reason to qualify the list, not to hide the alternatives.
- Cell influence: identify whether a high score depends on one or two cells. Subsampling helps reveal sensitivity, but a small number of cells also makes any stability estimate coarse.
- Biological and redundant features: consider whether a rare-cell signal is scientifically relevant and whether several selected genes carry nearly duplicate information. A purely numerical rank does not resolve either question.
The standardized excess-variance score used here is only one ranking approach. Binned dispersion compares variability among genes in similar expression bins; regularized dispersion stabilizes such comparisons. Negative-binomial deviance and Pearson residual approaches use model-based departures or residual variation. Each method has its own assumptions, scale, and tuning choices. Agreement across methods can support robustness, while disagreement can indicate sensitivity to modeling decisions; it is not a reason to average incomparable scores without justification.
Common mistakes include ranking raw variance without accounting for mean, applying the full-data eligible set unchanged in every subsample, and calling a selection frequency a probability that a gene is truly informative. A frequency is the proportion of specified resampling runs in which the gene was selected under the stated pipeline. It is conditional on the data, filter, ranking method, feature-count rule, and resampling design. It is not a biological probability or a p-value.
The final feature set is also not a list of genes proven to be markers, causal drivers, or differentially expressed. Selection is a preprocessing decision for a downstream analysis. It can miss rare signals, favor correlated or redundant genes, and change with the cell population or normalization. Report enough detail for another analyst to reproduce the decision: the expression scale and normalization, filter thresholds, trend and ranking rule, number of genes selected, resampling size and repetitions, whether filtering was reapplied, stability threshold, and any deliberate exceptions.
Check Your Understanding
Answer each question, then compare your reasoning with the brief solutions.
- Under the example’s filter, why is Gene D excluded even though it has a nonzero count in every cell?
- For Gene B, explain why its variance of 6 does not make it the most variable gene under the defined ranking.
- In six-of-eight subsampling, why can Gene C be eligible in some draws and ineligible in others?
- What does a 54% overall selection frequency for Gene C mean in this workflow, and what additional denominator clarifies its behavior?
- If only two features could be selected, which genes would the full-data ranking choose, and what limitation should be reported?
Solutions. (1) D never reaches count 2, so it fails the requirement of count at least 2 in at least two cells. (2) The ranking compares variance with its mean-dependent trend prediction; B has predicted variance about 13.06 and score about -1.954, while C’s score is about 11.18. (3) C has nonzero counts in exactly two cells and meets the detection rule only when both are included. (4) C is selected in 54 of all 100 draws; report that it was eligible in 54 draws and selected in 54 of those 54 eligible draws. (5) The top two full-data genes are C and E. C’s eligibility is sensitive to which cells are sampled, so the report should distinguish its full-data rank from its subsampling frequency and eligibility frequency.