Why Measure Stability?
A highly variable gene (HVG) ranking is calculated from a particular set of cells. If a modest change in the sampled cells substantially changes which genes are selected, a single ranked list can give a misleading impression of precision. Repeating the selection after subsampling asks a focused question: how consistently would this procedure select each gene if the observed cells were perturbed by the same sampling design?
This is a sensitivity analysis of a selection procedure, not a test of whether a gene is truly variable in a wider biological population. Its result depends on which cells are eligible for resampling, the subsample size, the ranking method, and the selection rule. The earlier material on choosing a feature count and selection threshold motivates fixing the number selected; the earlier material on preserving rare-cell signals is relevant when deciding how cells should be sampled.
Definition and Resampling Design
For example, \(f_g=0.80\) means that gene \(g\) was selected in 80% of the resampling replicates under this particular design. It does not mean there is an 80% probability that the gene is biologically important. Nor does it show that the gene will be selected in an independent experiment.
Before calculating frequencies, specify the complete procedure. This includes the eligible cells, subsample size, whether sampling is with or without replacement, any strata used in sampling, the expression scale and ranking method, and the selection rule. If exactly \(K\) genes are selected in every replicate, use that same \(K\) throughout. If the rule is instead a fixed threshold, state it and allow the number selected to vary. Resolve ties consistently.
Subsampling without replacement within a replicate avoids duplicating cells in that replicate. Independent replicates can nevertheless reuse cells from the full data set. If cells come from distinct samples or batches, a design may sample within each stratum to retain the intended representation. The choice should reflect the question: a design that changes cell-type composition is testing both ranking sensitivity and composition sensitivity.
What the Frequency Estimates
Imagine repeatedly drawing a subsample from the specified design and rerunning the fixed selection procedure on the observed data. Let \(p_g\) be the probability that gene \(g\) is selected in one such draw, conditional on the observed data and design. Then \(f_g\) estimates \(p_g\). It summarizes sensitivity to the resampling perturbation, not uncertainty about an unobserved biological truth.
When the \(B\) subsample draws are independently randomized from the same design, the indicators \(I_{1g},\ldots,I_{Bg}\) are independent Bernoulli variables with success probability \(p_g\). Therefore,
The first equality follows from linearity of expectation: each indicator has expectation \(p_g\), so \(E(f_g\mid D)=B^{-1}Bp_g=p_g\). For the variance, independence makes the variance of the sum equal to the sum of the variances. Each indicator has variance \(p_g(1-p_g)\), giving \(\operatorname{Var}(f_g\mid D)=B^{-2}Bp_g(1-p_g)\). Replacing the unknown \(p_g\) by \(f_g\) gives the Monte Carlo standard error estimate \(\sqrt{f_g(1-f_g)/B}\).
This standard error describes the numerical precision of the estimated frequency due to using a finite number of random subsamples. It is not a confidence interval for biological reproducibility. Increasing \(B\) reduces Monte Carlo error but does not undo limitations of the observed cells, the selection rule, or the resampling design.
If exactly \(K\) genes are selected per replicate, then \(\sum_g I_{bg}=K\) for every \(b\), and thus \(\sum_g f_g=K\). Consequently, a low frequency for one gene is partly a result of competition among genes for a limited number of selection slots. Frequencies in a fixed-\(K\) list should be interpreted as stability of selection under that constraint.
Worked Example: Computing Selection Frequencies
Suppose a hypothetical analysis draws \(B=10\) independent subsamples, reranks the genes in each, and selects the top \(K=3\) each time. The selection counts are: Gene A, 10; Gene B, 8; Gene C, 6; Gene D, 4; and Gene E, 2.
For Gene C, the frequency is \(f_C=6/10=0.60\). Gene A has \(f_A=10/10=1.00\), while Gene E has \(f_E=2/10=0.20\). Thus A was selected in every replicate, C in six, and E in two. Under this subsampling design, A is more stable than C, and C more stable than E.
The counts sum to \(10+8+6+4+2=30\), which equals \(BK=10\times3\). This is a useful check: selecting exactly three genes in each of ten replicates must produce 30 total gene selections. The conclusion is conditional on the fixed ranking method, top-three rule, eligible cells, and subsample size; it does not establish that A is more biologically relevant than the other genes.
Comparing Whole Selected Sets
Per-gene frequencies show which genes recur, but a summary of overlap between pairs of selected sets can show how much the lists change as a whole. For sets \(A\) and \(C\), their Jaccard similarity is the number of shared genes divided by the number of distinct genes appearing in either set:
The value ranges from 0 (no shared genes) to 1 (identical sets). With equal-sized top-\(K\) lists, it gives a compact description of set overlap, but it does not reveal which genes are consistently selected across all replicates. Pairwise Jaccard values should therefore complement, not replace, gene-level frequencies.
Worked Example: Interpreting Set Overlap
In two hypothetical replicates, the selected sets are \(A=\{G1,G2,G3,G4\}\) and \(C=\{G2,G3,G5,G6\}\). The shared genes are \(G2\) and \(G3\), so \(|A\cap C|=2\). The union contains \(G1,G2,G3,G4,G5,G6\), so \(|A\cup C|=6\).
Therefore \(J(A,C)=2/6=0.333\). The two lists share one third of the genes appearing in either list. This indicates substantial set variation between these two replicates, but the value alone cannot tell whether the same two genes recur in most other replicates. For that question, inspect the selection frequencies over all replicates.
Monte Carlo Precision and Choice of Replicate Count
A frequency near 0.5 is intrinsically more variable across independent resampling runs than a frequency near 0 or 1, for the same \(B\), because \(p_g(1-p_g)\) is largest at \(p_g=0.5\). The number of replicates should be large enough that the estimated frequencies are not dominated by Monte Carlo noise. It is often useful to calculate the Monte Carlo standard error for genes near a practical stability boundary, rather than choosing \(B\) without checking precision.
Worked Example: Monte Carlo Standard Error
Suppose a hypothetical gene is selected in 15 of \(B=20\) independently randomized subsamples. Its observed frequency is \(f_g=15/20=0.75\). The estimated Monte Carlo standard error is
The estimated frequency is 0.75, with Monte Carlo standard error about 0.097, or 9.68 percentage points. This indicates that 20 runs provide limited numerical precision for the frequency. It is not a statement that the gene has a 75% probability of being a true HVG, and the standard error does not account for uncertainty from collecting a different biological sample.
The approximation is appropriate here because each independently randomized replicate applies the same rule and design to the same observed data. If the subsampling scheme changes between replicates, the indicators need not share one common selection probability, and this simple standard error no longer has its stated interpretation.
Diagnostics, Failure Modes, and Alternatives
Review the resampling design and the resulting selections before interpreting a stability summary. Useful checks include verifying the number of cells in every replicate, confirming that the same preprocessing and ranking rule were applied, and checking that the selected-gene count matches the stated rule. For a fixed-\(K\) selection, confirm that every replicate selects \(K\) genes and that the frequencies sum to \(K\).
Cell composition can be a major source of instability. In a hypothetical data set of 100 cells containing 5 cells of a rare type, a simple random subsample of 40 cells has probability
Thus, about 7.3% of such subsamples would contain no rare-type cells at all. A gene that varies mainly within that cell type may then be selected inconsistently—not necessarily because its within-type variability is weak, but because the relevant cells are sometimes absent. Depending on the aim, stratified subsampling or a larger subsample may better preserve the population of interest. Report the design so readers can distinguish these explanations.
Cells from the same biological sample are not independent biological replicates. Resampling cells from one observed collection measures sensitivity to cell-level perturbations conditional on that collection; it does not quantify variation across donors or experiments. When the scientific target is reproducibility across biological samples, the resampling unit and analysis must reflect those samples.
Subsampling without replacement is a common choice for perturbing the observed cell set without duplicating cells within a run. Bootstrap resampling draws cells with replacement and can be used as a different perturbation, but its duplicate-cell structure changes the ranking input. Repeated split-half analyses are another option, though each half uses fewer cells and so may be more sensitive to sample size. These approaches answer related but not identical sensitivity questions; their results should not be treated as interchangeable.
- Unstable membership: Many genes have intermediate frequencies, or pairwise set overlap is low. This indicates sensitivity under the chosen perturbation, not proof that the ranking method is invalid.
- Rare-cell loss: Replicates differ in whether a rare population is represented. Consider whether the design should preserve cell-type or sample strata.
- Changing procedures: If normalization, ranking, or the number selected changes between runs, observed frequency mixes sampling sensitivity with procedure changes. Hold these choices fixed when measuring the stability of one procedure.
- Too few runs: Frequencies jump noticeably when additional replicates are added. Increase \(B\) and monitor Monte Carlo precision.
Common Mistake
The central practical value of stability analysis is prioritization and diagnosis. Genes selected nearly every time are robust to the tested cell-level perturbation; genes selected intermittently deserve closer inspection of their rank, expression, and representation across relevant cell groups. A stable list can still reflect systematic bias, while an unstable list may reflect a deliberately challenging sampling design. Stability is evidence about sensitivity, not a substitute for biological interpretation.
Check Your Understanding
Answer each question, then compare your reasoning with the solutions.
- A gene is selected in 18 of 30 subsamples. What is its selection frequency, and what does it mean?
- In every replicate, exactly 50 genes are selected. What should the sum of gene-selection frequencies equal?
- A gene has frequency 0.90. Does this mean there is a 90% probability it is a true HVG? Explain.
- Why might a gene expressed mainly in a rare cell type have a low or unstable selection frequency under unstratified subsampling?
- For two selected sets with 6 shared genes and 14 genes in their union, calculate the Jaccard similarity.
Solutions. (1) \(f_g=18/30=0.60\); the gene was selected in 60% of replicates under the stated design and procedure. (2) The frequencies sum to 50, since each replicate contributes 50 selections. (3) No. The frequency describes conditional stability under resampling, not a probability that the gene is biologically true. (4) Some subsamples may contain few or no cells of that type, changing the information available to the ranking procedure. (5) \(J=6/14\approx0.429\), so the sets share about 42.9% of the genes in their union.