Highly Variable Gene Selection for scRNA-seq PCA: A Foundational Guide for Robust Single-Cell Analysis

Abstract

This article provides a comprehensive guide for researchers and bioinformaticians on the critical practice of filtering highly variable genes (HVGs) prior to principal component analysis (PCA) in single-cell RNA sequencing (scRNA-seq). Covering foundational principles to advanced applications, we detail the necessity of HVG selection to combat the 'curse of dimensionality' and enhance biological signal capture. We explore and compare established methods from Seurat and Scanpy, alongside innovative approaches like GLP and model-based dimensionality reduction. The content further addresses common pitfalls, optimization strategies for challenging datasets, and quantitative frameworks for method validation. This guide is indispensable for ensuring accurate clustering, trajectory inference, and differential expression analysis in biomedical and clinical research.

Why HVG Selection is Your First Defense Against the Curse of Dimensionality

Understanding the 'Curse of Dimensionality' in Single-Cell Data

Single-cell RNA sequencing (scRNA-seq) represents a transformative technology in biomedical research, enabling the quantification of gene expression across thousands of individual cells simultaneously [1]. While this generates unprecedented opportunities for discovering cellular heterogeneity, it also introduces significant computational challenges, primarily the curse of dimensionality (COD). In the context of scRNA-seq analysis, the curse of dimensionality describes statistical problems that arise when analyzing data with extremely high dimensions—typically tens of thousands of genes measured per cell [1] [2]. This phenomenon differs from the computational complexity curse of dimensionality referenced in computer science, instead focusing on how noise accumulation in high dimensions obscures true biological signals [1].

The fundamental issue stems from the nature of scRNA-seq data itself, which is characterized by substantial technical noise including detection failures (dropouts) and substantial variations in detection levels across transcripts [1]. As scRNA-seq data measures over 10,000 genes per cell, with each feature bearing technical noise, the accumulation of these noise components creates severe problems in downstream analyses [1]. The high dimensionality means that even small noise components in each feature accumulate to dominate the distance metrics and statistical analyses essential for interpreting single-cell data. Understanding and addressing COD is therefore not merely an optional optimization but a fundamental prerequisite for extracting biologically meaningful insights from scRNA-seq experiments.

Theoretical Framework: How the Curse of Dimensionality Manifests in scRNA-seq Data

Types of Curse of Dimensionality in scRNA-seq

In scRNA-seq analysis, the curse of dimensionality manifests through three distinct but interconnected statistical problems that adversely affect data interpretation [1]:

• COD1: Loss of Closeness - This phenomenon occurs when noise accumulation in high-dimensional space obscures the true distances between cells. In high-dimensional data with noise, the errors in Euclidean distance grow with the dimensionality, eventually overwhelming the actual biological differences between cells. This impairment makes detailed classification of high-dimensional data impossible, as evidenced by the increasingly elongated "legs" in dendrograms of unsupervised hierarchical clustering as dimensionality increases [1]. Similar problems affect other distance metrics, including correlation distance, causing conventional data analysis methods based on distances to fail in identifying true data structures [1].

• COD2: Inconsistency of Statistics - Noise accumulation creates adverse effects on data statistics, particularly contribution rates in principal component analysis (PCA) and clustering metrics like the Silhouette score [1]. High-dimensional statistics theory demonstrates that data variances of PCA-transformed data (eigenvalues of the data covariance matrix) may not converge to true variances for high-dimensional data with noise and low sample sizes [1]. This statistical inconsistency leads to false inferences during data analysis and interpretation.

• COD3: Inconsistency of Principal Components - When substantial variation exists in noise scale across features, as occurs in scRNA-seq data with random sampling and low detection rates, PCA structures become distorted [1]. In such cases, principal components become influenced by non-biological information such as sequencing depth (total UMI counts per cell) or the number of detected genes rather than genuine biological variation [1]. This represents a fundamental challenge for dimensionality reduction techniques that rely on PCA as a foundational step.

Mathematical Foundations

The curse of dimensionality in scRNA-seq data arises from the intrinsic properties of high-dimensional spaces. In datasets where the number of features (genes) vastly exceeds the number of observations (cells), the concept of distance becomes diluted as all cells appear nearly equidistant from one another [2]. This occurs because the accumulation of technical noise across thousands of genes causes distance errors to grow proportionally with dimensionality, eventually obscuring true biological differences between cells [1].

The technical noise in scRNA-seq data originates from multiple sources, including variations in copying and sequencing errors during data generation [1]. For UMI-based protocols, this noise can be defined as the difference between observed UMI counts and true RNA counts divided by their total counts [1]. As each cell typically expresses more than 10,000 genes, with noise affecting all expressed genes, the cumulative effect creates the statistical challenges characteristic of COD.

Table 1: Types of Curse of Dimensionality in scRNA-seq Data

Type	Key Characteristic	Impact on Analysis
Loss of Closeness (COD1)	Noise accumulation obscures true distances between cells	Impairs clustering and detailed classification of cells
Inconsistency of Statistics (COD2)	Statistical measures become unreliable	Causes false inference in PCA contribution rates and clustering metrics
Inconsistency of Principal Components (COD3)	PCA structures distorted by technical factors	Principal components reflect technical artifacts rather than biology

The Critical Role of Highly Variable Gene Selection

Theoretical Rationale for Feature Selection

Filtering highly variable genes (HVGs) before applying dimensionality reduction techniques like Principal Component Analysis (PCA) represents a crucial strategy for mitigating the curse of dimensionality in scRNA-seq data [2]. The fundamental premise is that not all genes are equally informative for distinguishing cell types or states—while the transcriptome may contain over 20,000 genes, only a subset exhibits sufficient variability across cells to contribute meaningfully to population structure [2]. Selecting this informative subset serves multiple critical functions in combating COD.

From a statistical perspective, HVG selection directly addresses the noise accumulation problem by reducing the number of dimensions in which random technical variation can obscure biological signals. By focusing analysis on genes demonstrating high cell-to-cell variation relative to technical noise, researchers effectively filter out genes that contribute primarily noise to distance calculations and statistical summaries [1]. This process enhances the signal-to-noise ratio in downstream analyses, making biological patterns more discernible. Additionally, by reducing the dimensionality from tens of thousands of genes to a more manageable subset of typically 1,000-5,000 HVGs, computational efficiency is dramatically improved without sacrificing biological resolution [2].

Integration with Downstream Analysis

The selection of highly variable genes serves as a critical preprocessing step that fundamentally enables effective application of PCA and other dimensionality reduction techniques [2]. When PCA is applied to full scRNA-seq datasets without HVG selection, the principal components often capture technical artifacts such as sequencing depth or mitochondrial percentage rather than biological variation [1] [3]. This occurs because the substantial technical noise in lowly expressed or non-informative genes dominates the total variation in the dataset.

In contrast, when PCA is applied to a curated set of highly variable genes, the resulting principal components more reliably capture biologically meaningful variation [2]. The top PCs typically represent transcriptional programs that distinguish cell types, activation states, or other biologically significant patterns. This focused approach maintains the advantages of PCA—including noise reduction and data compression—while minimizing the distorting effects of high-dimensional noise. The selected HVGs subsequently form the feature set for clustering and visualization techniques like t-SNE and UMAP, which also benefit from reduced dimensionality and enhanced signal-to-noise ratios [2].

Diagram 1: scRNA-seq workflow with curse of dimensionality challenges. HVG selection mitigates COD effects before dimensionality reduction.

Experimental Protocols for Mitigating Dimensionality Challenges

Protocol 1: Standardized Workflow for HVG Selection and PCA

This protocol provides a robust method for selecting highly variable genes and performing PCA to mitigate dimensionality challenges in scRNA-seq analysis, compatible with both Scanpy and Seurat frameworks.

Materials and Reagents:

• Processed scRNA-seq count matrix (UMI-based)
• Computational environment with Scanpy (Python) or Seurat (R) installed
• High-performance computing resources for large datasets

Step-by-Step Procedure:

• Data Preprocessing and Normalization

  • Begin with a quality-controlled UMI count matrix. Cells with excessive mitochondrial reads or low gene counts should be filtered out.
  • Normalize counts to account for sequencing depth variations using the method appropriate for your analysis framework:
    • Scanpy: Use sc.pp.normalize_total() function to normalize total counts per cell followed by sc.pp.log1p() for log transformation [2].
    • Seurat: Employ the NormalizeData() function with the "LogNormalize" method [3].

• Highly Variable Gene Selection

  • Identify genes with high cell-to-cell variation using the following approaches:
    • Scanpy: Apply sc.pp.highly_variable_genes() with the 'seurat' method to select 2,000-5,000 highly variable genes [2].
    • Seurat: Use the FindVariableFeatures() function with the "vst" method to identify the top 2,000 variable genes [3].
  • Subset the dataset to include only selected highly variable genes for downstream analysis.

• Principal Component Analysis

  • Scale the data to have zero mean and unit variance using sc.pp.scale() in Scanpy or ScaleData() in Seurat.
  • Perform PCA on the scaled data of highly variable genes using sc.tl.pca() in Scanpy or RunPCA() in Seurat.
  • Determine the number of significant PCs to retain for downstream analysis using elbow plots (percentage variance explained) or JackStraw analysis in Seurat.

• Downstream Applications

  • Use the selected PCs for graph-based clustering (sc.tl.louvain in Scanpy or FindClusters in Seurat).
  • Generate 2D visualizations using UMAP (sc.tl.umap in Scanpy or RunUMAP in Seurat) or t-SNE.
  • The reduced dimensional space now enables effective identification of cell populations and states.

Troubleshooting Tips:

• If PCA results appear driven by technical artifacts (e.g., sequencing depth), revisit normalization steps and consider alternative methods like SCTransform [3].
• If biological signal remains weak after HVG selection, adjust the number of variable genes selected or employ different selection methods.
• For large datasets exceeding 50,000 cells, consider approximate PCA methods for computational efficiency.

Protocol 2: Advanced Noise Reduction Using RECODE

For datasets where standard normalization and HVG selection prove insufficient, RECODE (Resolution of the Curse of Dimensionality) offers a specialized approach to directly address COD in scRNA-seq data with unique molecular identifiers [1].

Principles and Applications: RECODE operates on a fundamentally different principle from imputation methods, focusing specifically on resolving the statistical problems associated with high-dimensional noise rather than attempting to recover missing values [1]. The method is particularly valuable for analyzing subtle cell fate transitions, identifying rare cell populations, and working with lowly expressed genes where technical noise often masks biological signals.

Implementation Workflow:

• Data Requirements and Preprocessing

  • Input: UMI-based scRNA-seq count matrix without prior normalization
  • Filter out low-quality cells and genes with zero counts across all cells
  • No additional normalization is required before applying RECODE

• RECODE Processing

  • Apply the RECODE algorithm to the count matrix using available implementations
  • The method operates by decomposing the high-dimensional noise structure and resolving the distance impairments caused by noise accumulation
  • The output is a noise-reduced expression matrix preserving all original genes

• Downstream Analysis

  • Proceed with standard HVG selection and PCA on the RECODE-processed data
  • The resulting dimensional reductions demonstrate enhanced ability to resolve subtle cellular states and transitions
  • All genes, including lowly expressed ones, remain available for analysis

Advantages and Limitations:

• Advantages: Parameter-free, data-driven, deterministic algorithm; preserves all gene information; specifically addresses COD rather than masking it [1]
• Limitations: Primarily designed for UMI-based data; may not be suitable for non-UMI protocols; computational requirements may be higher than standard normalization

Table 2: Comparison of Dimensionality Reduction Methods for scRNA-seq Data

Method	Primary Function	Key Advantages	Limitations
PCA on HVGs	Linear dimensionality reduction	Computationally efficient, highly interpretable	Assumes linear relationships, may miss nonlinear structures
t-SNE	Nonlinear visualization	Excellent cluster separation in 2D/3D	Computational intensity, parameters sensitive, global structure loss
UMAP	Nonlinear dimensionality reduction	Preserves global structure, faster than t-SNE	Parameter sensitivity, potential over-interpretation of clusters
RECODE	COD resolution for UMI data	Parameter-free, preserves all genes, resolves COD	Designed for UMI data, computationally intensive
scVI	Deep generative modeling	Integrates normalization and batch correction	Complex implementation, requires substantial data

The Scientist's Toolkit: Essential Research Reagents and Computational Solutions

Table 3: Essential Tools for scRNA-seq Analysis and COD Mitigation

Tool/Resource	Function	Application Context
Cell Ranger	Processing raw sequencing data	Generates count matrices from 10x Genomics FASTQ files [4]
Seurat	Comprehensive scRNA-seq analysis in R	Data normalization, HVG selection, clustering, and visualization [4]
Scanpy	Scalable scRNA-seq analysis in Python	Preprocessing, PCA, clustering, UMAP visualization for large datasets [4] [2]
SCTransform	Regularized negative binomial regression	Normalization and variance stabilization addressing technical confounding [3]
BASiCS	Bayesian hierarchical modeling	Quantifies technical and biological noise, identifies highly variable genes [5]
Harmony	Batch effect correction	Integrates datasets across experiments while preserving biological variation [4]
Velocyto	RNA velocity analysis	Models cellular dynamics by quantifying spliced/unspliced transcripts [4]
CellBender	Ambient RNA removal	Uses deep learning to remove background noise from droplet-based data [4]

Diagram 2: Solution strategies addressing the curse of dimensionality in scRNA-seq data analysis.

The curse of dimensionality represents a fundamental challenge in single-cell RNA sequencing data analysis, with profound implications for interpreting cellular heterogeneity [1] [2]. Through the systematic implementation of highly variable gene selection prior to dimensionality reduction techniques like PCA, researchers can effectively mitigate the distorting effects of high-dimensional noise while preserving biological signal. The protocols and methodologies outlined herein provide a robust framework for addressing COD across diverse experimental contexts, from standard cell type identification to sophisticated analyses of subtle cellular transitions.

As single-cell technologies continue to evolve, generating increasingly complex and high-dimensional datasets, the strategic implementation of COD-aware analytical approaches will become ever more critical. By recognizing the curse of dimensionality not as a peripheral technical concern but as a central determinant of analytical success, researchers can design more effective computational workflows, extract more meaningful biological insights, and ultimately advance our understanding of cellular systems in health and disease.

In single-cell RNA sequencing (scRNA-seq) analysis, the process of feature selection serves as a critical gateway to meaningful biological discovery. The identification of highly variable genes (HVGs) represents a fundamental step that directly influences all subsequent analytical outcomes, from dimensionality reduction to cell type classification. This procedure navigates a delicate balance: eliminating technical noise inherent in sparse single-cell data while preserving the biological heterogeneity that reveals true cellular function and identity. The high dimensionality and sparsity of scRNA-seq data pose substantial challenges for downstream analysis, making feature selection essential for reducing dimensionality and enhancing interpretability [6] [7]. Current computational analyses often overlook the rich information encoded by variability within single-cell gene expression data by focusing exclusively on mean expression, potentially missing crucial biological insights embedded within this variability [8].

The strategic selection of HVGs has evolved beyond mere variance filtering to incorporate sophisticated relationships between expression patterns, dropout rates, and biological significance. As scRNA-seq technologies advance, enabling researchers to capture increasingly complex cellular landscapes, the methods for identifying informative genes must similarly evolve to distinguish meaningful biological signals from technical artifacts more effectively. This article explores contemporary methodologies for HVG selection, their impact on downstream analyses, and provides detailed protocols for implementation—all within the critical context of preparing data for principal component analysis in research environments.

Current Methodologies for HVG Selection

The landscape of HVG selection methods has diversified significantly, moving beyond simple variance-based approaches to incorporate more sophisticated statistical models and computational frameworks. These methods can generally be categorized into several distinct paradigms, each with unique strengths for addressing specific analytical challenges.

Statistical and Distributional Model-Based Approaches

Traditional statistical methods leverage predefined distributional assumptions to identify genes exhibiting greater variability than expected by technical noise alone. The widely-used Seurat package offers two principal approaches: VST (Variance Stabilizing Transformation), which computes variance through standardized expression level, and SCTransform, which utilizes Pearson residuals derived from a generalized linear model to compute variance [6] [7]. Similarly, M3Drop fits a Michaelis-Menten function to model the relationship between mean expression levels and dropout rates, then employs a t-test to evaluate the significance of each gene, while NBDrop calculates expected dropout rates through a negative binomial distribution [6]. These methods operate under specific statistical priors, which may not always align with the intricate structure of single-cell data, potentially leading to biases in feature selection [6].

Clustering and Graph-Based Methods

These approaches identify important genes by constructing gene-gene or cell-cell relationship networks and applying clustering or graph-theoretical techniques. FEAST assesses gene significance using F-statistics derived from consensus clusters, while HRG constructs a cell-cell similarity network to select genes exhibiting regional expression patterns [6]. GeneBasisR selects genes with the maximum distance between the true manifold and the manifold constructed using the selected gene panel iteratively [6]. Both CellBRF and DELVE rely on reconstructing the graph of cell neighborhoods using the k-nearest neighbors (k-NN) algorithm, differing primarily in their strategies for adjusting cell types used to assess gene importance [6]. These methods can capture complex relationships but may be compromised by the prevalent dropout noise in scRNA-seq data, resulting in inaccurate correlation estimates and suboptimal clustering outcomes [6].

Regression-Based and Adaptive Frameworks

More recent approaches leverage flexible regression techniques to model the relationship between gene expression characteristics without strong distributional assumptions. The newly developed GLP (genes identified through LOESS with positive ratio) algorithm utilizes optimized LOESS regression to precisely capture the relationship between gene average expression level and positive ratio while minimizing overfitting [6] [7]. This method integrates a two-step LOESS regression procedure that applies Tukey's biweight robust statistical method to identify outlier genes, then assigns these outliers zero weights to minimize their influence [6]. Similarly, spline-DV employs a nonparametric, model-free framework that uses three gene-level metrics—mean expression, coefficient of variation (CV), and dropout rate—to create a 3D model for estimating gene expression variability [8]. These approaches aim to more effectively distinguish biologically meaningful genes while mitigating the adverse effects of dropout noise and sparsity [6] [8].

Random Matrix Theory and Sparse PCA Approaches

Methods leveraging Random Matrix Theory (RMT) represent another innovative approach to handling high-dimensional single-cell data. RMT-guided sparse PCA uses a biwhitening method, inspired by the Sinkhorn-Knopp algorithm, that simultaneously stabilizes variance across genes and cells [9]. This enables the use of an RMT-based criterion to automatically select the sparsity level, rendering sparse PCA nearly parameter-free while retaining the interpretability of PCA with improved robustness [9]. This mathematically grounded approach systematically improves the reconstruction of the principal subspace and consistently outperforms PCA-, autoencoder-, and diffusion-based methods in cell-type classification tasks across diverse single-cell RNA-seq technologies [9].

Quantitative Comparison of HVG Methods

The performance of different HVG selection methods has been systematically evaluated across multiple studies, providing valuable insights for method selection. The following tables summarize key performance metrics and characteristics of major approaches.

Table 1: Performance Comparison of HVG Selection Methods Across Benchmarking Studies

Method	Algorithm Type	Key Features	Performance Advantages	Limitations
GLP [6] [7]	Optimized LOESS regression	Uses positive ratio vs. expression relationship; BIC-optimized bandwidth	Consistently outperforms 8 other methods across ARI, NMI, and silhouette metrics [6]	Computational intensity; requires minimum 3 cells expressing each gene
spline-DV [8]	Spline-based differential variability	3D model (mean, CV, dropout); identifies differential variability between conditions	Detects functional genes in obesity, fibrosis, cancer; validates well in simulations [8]	Complex implementation; requires specialized software
RMT-guided sPCA [9]	Random Matrix Theory + sparse PCA	Biwhitening normalization; automatic sparsity parameter selection	Outperforms PCA, autoencoders, diffusion methods in cell classification [9]	Limited to biwhitened data; mathematical complexity
Seurat VST [10]	Variance stabilizing transformation	Mean-variance relationship; standardized expression	Established benchmark; effective for general purpose integration [10]	Can miss biologically relevant variable genes
scSGS [11]	Stochastic gene silencing pattern analysis	Leverages natural transcriptional bursting; binarized expression	Reveals regulatory relationships; avoids survivorship bias of KO studies [11]	Requires sufficient cells for active/silenced classification

Table 2: Impact of Feature Selection on Downstream Analysis Performance

Analysis Task	Optimal Method	Key Metric	Performance Improvement	Reference
Data Integration	Batch-aware HVG selection	Batch correction + bio conservation	High-quality integrations with effective batch removal [10]	Significant improvement over random features
Query Mapping	2,000 HVGs (scanpy)	Mapping accuracy	Better preservation of biological variation [10]	30-50% improvement over suboptimal selection
Cell Type Classification	RMT-guided sPCA	Classification accuracy	Consistently outperforms standard PCA [9]	Robust across 7 scRNA-seq technologies
Trajectory Inference	GLP	Trajectory accuracy	Enhances downstream trajectory inference [6]	Improved pseudotemporal ordering
Differential Expression	spline-DV	Functional relevance	Identifies representative, functionally relevant genes [8]	Captures obesity, fibrosis, cancer mechanisms

Detailed Experimental Protocols

Protocol 1: GLP Implementation for HVG Selection

The GLP algorithm provides a robust feature selection approach that leverages optimized LOESS regression to capture the relationship between gene average expression level and positive ratio.

Materials and Reagents:

• scRNA-seq count matrix (genes × cells)
• Computational environment: R or Python with statistical libraries
• Hardware: Standard workstation (8+ GB RAM for datasets <10,000 cells)

Procedure:

• Data Preprocessing

  • Begin with a raw count matrix of dimensions g × c, where g represents genes and c represents cells.
  • Filter out genes detected in fewer than 3 cells to ensure reliable positive ratio calculation [6].
  • Do not apply normalization at this stage, as the method operates on raw counts.

• Calculate Expression Metrics

  • For each gene j, compute the average expression level (λj) using the formula: λj = (1/c) × ΣXij (where Xij is the count for gene j in cell i) [6].
  • Calculate the positive ratio (fj) for each gene using: fj = (1/c) × Σmin(1, Xij) [6].
  • Construct a two-dimensional distribution with positive ratio (f) as the independent variable and average expression (λ) as the dependent variable.

• Optimize LOESS Regression Parameters

  • Define the smoothing parameter (⍺) range from 0.01 to 0.1 with a step size of 0.01 [6].
  • For each ⍺ value, perform LOESS regression and calculate the Bayesian Information Criterion (BIC): BIC = c × ln(RSS/c) + k × ln(c) [6] where RSS is the residual sum of squares, c is the number of observations, and k is the degrees of freedom.
  • Select the ⍺ value that minimizes BIC to prevent overfitting while capturing data structure.

• Two-Step LOESS Regression with Outlier Handling

  • Perform initial LOESS regression using the optimized ⍺ parameter.
  • Apply Tukey's biweight robust statistical method to identify outlier genes [6].
  • Assign zero weights to identified outliers and reperform LOESS regression.
  • Select genes with average expression levels significantly higher than LOESS-predicted values based on positive ratios.

• Downstream Application

  • Use the selected HVGs for PCA input, clustering, or trajectory analysis.
  • The typical number of selected genes ranges from 1,000-5,000 depending on dataset complexity.

Figure 1: GLP Workflow for HVG Selection. This diagram illustrates the step-by-step process for implementing the GLP algorithm, from initial data filtering through optimized LOESS regression to final HVG selection.

Protocol 2: spline-DV for Differential Variability Analysis

The spline-DV method identifies genes exhibiting significantly different variability between experimental conditions, capturing biological insights beyond mean expression differences.

Materials and Reagents:

• scRNA-seq data from two experimental conditions
• R environment with scGEAToolbox [8]
• Graphics capabilities for 3D visualization (optional)

Procedure:

• Data Preparation and Metric Calculation

  • For each condition, preprocess data separately using standard normalization (CP10k + log1p) [12].
  • Calculate three gene-level metrics for each condition:
    • Mean expression (x-coordinate)
    • Coefficient of variation (y-coordinate)
    • Dropout rate (z-coordinate) [8].

• Spline Curve Fitting

  • Generate independent spline-fit curves for each condition within the 3D coordinate space.
  • Use the spline-fit algorithm from scGEAToolbox to model the relationship between expression statistics [8].

• Vector Deviation Calculation

  • For each gene in condition 1, identify its position in the 3D space based on observed metrics.
  • Calculate vector v₁ from the nearest point on the condition 1 spline curve to the gene's position.
  • Repeat for condition 2 to calculate vector v₂.
  • Compute the differential variability vector: dv = v₂ - v₁ [8].

• DV Scoring and Gene Ranking

  • Calculate the magnitude of the dv vector as the DV score for each gene.
  • Rank all genes based on their DV scores (highest to lowest magnitude).
  • Select top-ranked DV genes for functional analysis.

• Biological Validation

  • Perform functional enrichment analysis on top DV genes.
  • Validate findings through experimental follow-up or comparison with existing literature.

Figure 2: spline-DV Analysis Workflow. This diagram illustrates the process for identifying differentially variable genes between two experimental conditions using 3D metric analysis and vector calculations.

The Scientist's Toolkit: Essential Research Reagents and Computational Solutions

Table 3: Essential Research Reagents and Computational Tools for HVG Analysis

Category	Item/Software	Specification/Version	Function in HVG Analysis	Implementation Notes
Wet Lab Reagents	10x Genomics Single Cell Kit	3' Gene Expression v3	scRNA-seq library preparation	Enables high-quality input data for HVG detection
Cell Viability Reagents	Propidium Iodide/Annexin V	Flow cytometry grade	Exclusion of dead/dying cells	Reduces technical variability in expression
RNA Quality Control	Bioanalyzer RNA Integrity	RIN >8.0	Pre-seq quality assessment	Ensures high-quality input material
Computational Tools	Seurat R Package	v4.0+	Standard HVG selection	Provides VST, SCTransform methods [6] [10]
Programming Environment	R/Python	R 4.0+/Python 3.8+	Statistical computing	Essential for custom HVG implementations
Specialized Algorithms	scGEAToolbox	Custom installation	spline-DV implementation	Enables differential variability analysis [8]
Data Integration	scanpy	v1.8+	Python-based HVG selection	Batch-aware feature selection [10]
Visualization	ggplot2/Matplotlib	Current versions	Visualization of HVG results	Critical for quality assessment

Integration with PCA and Impact on Downstream Analysis

The selection of HVGs profoundly influences PCA outcomes and all subsequent analytical steps in single-cell workflows. When appropriately implemented, HVG filtering enhances PCA performance by focusing on biologically informative dimensions while reducing technical noise.

The critical importance of proper feature selection was demonstrated in a comprehensive benchmark study examining scRNA-seq data integration and querying [10]. This research revealed that feature selection methods significantly affect integration performance, with highly variable feature selection proving effective for producing high-quality integrations [10]. The study further provided guidance on the number of features selected, batch-aware feature selection, lineage-specific feature selection, and the interaction between feature selection and integration models [10].

For PCA preparation, the RMT-guided sparse PCA approach offers particular advantages. By employing a biwhitening method that simultaneously stabilizes variance across genes and cells, this technique enables automatic sparsity parameter selection using RMT-based criteria [9]. This method maintains PCA's interpretability while providing robust performance across diverse single-cell RNA-seq technologies, consistently outperforming standard PCA, autoencoder-, and diffusion-based methods in cell-type classification tasks [9].

The strategic selection of HVGs before PCA extends beyond technical improvements to biological discovery. The "variation-is-function" hypothesis proposes that cell-to-cell gene expression variability is key to population-level cellular functions [8]. This perspective is supported by empirical evidence showing that HVGs in homogeneous cellular populations participate in biological processes and provide molecular functions specific to their respective cell types [8]. Interestingly, most HVGs are not highly expressed, whereas highly expressed genes—such as housekeeping genes and marker genes—tend not to convey cell-type specific functions [8]. This understanding reinforces the importance of variability-based approaches for capturing functional relevance in single-cell studies.

When designing analytical workflows, researchers should align their HVG selection strategy with specific research goals. For general atlas-building and cell type identification, batch-aware HVG selection with 2,000-3,000 features typically optimizes the balance between biological signal preservation and computational efficiency [10]. For investigating specific biological processes or responses, specialized methods like spline-DV or GLP may yield more targeted insights [6] [8]. Regardless of the specific method chosen, the critical role of thoughtful HVG selection in shaping all downstream analytical outcomes cannot be overstated.

Feature selection, specifically the identification of highly variable genes (HVGs), is a critical and foundational step in the analysis of single-cell RNA sequencing (scRNA-seq) data. This process reduces the immense dimensionality of the data—where tens of thousands of genes are measured across thousands of cells—to a manageable set of the most informative features, thereby enhancing computational efficiency and highlighting biological signal over technical noise [7]. The core assumption underpinning all HVG selection methods is that genes exhibiting high cell-to-cell variation in expression are more likely to be associated with biologically meaningful heterogeneity, such as differences in cell type, state, or response to stimuli, rather than mere technical artifacts [7] [13]. The precise definition of "variability," however, differs across methodologies, leading to a variety of algorithms and approaches. This article explores the core assumptions and practical protocols of the dominant HVG selection methods, providing a structured guide for researchers in drug development and biomedical science to implement these techniques effectively within the context of a broader research workflow that includes downstream dimensionality reduction via Principal Component Analysis (PCA).

Core Assumptions and Methodologies in HVG Selection

The landscape of HVG selection methods can be broadly categorized by the statistical principles and assumptions they employ. The table below summarizes the core assumptions and techniques of several key approaches.

Table 1: Core Assumptions and Techniques of HVG Selection Methods

Method / Flavor	Core Assumption	Mathematical Foundation	Key Technique
Seurat (vst, mvp) [14] [13]	Genes outliers in mean-variance relationship are biologically relevant.	Models mean-variance trend using local regression (LOESS).	Identifies genes with high normalized dispersion (vst) or z-scores (mvp) from the fitted trend.
Scanpy (seurat) [15] [16]	Reproduces Seurat's mean-variance approach.	Normalizes dispersions for genes within bins of mean expression.	Selects genes with high normalized dispersion based on user-defined cutoffs.
Seurat v3 / Scanpy (seurat_v3) [15]	Variance of stabilized count (variance stabilizing transform) is informative.	Standardizes gene expression using a regularized standard deviation.	Ranks genes by variance of standardized values; selects top N genes.
GLP (Novel Method) [7]	Relationship between positive ratio (f) and mean expression (λ) is locally predictable; genes above this trend are important.	Optimized LOESS regression on (positive ratio, mean expression) space.	Uses Bayesian Information Criterion for bandwidth selection; identifies genes with λ significantly > predicted λ.

The foundational methods, such as those implemented in Seurat and Scanpy, primarily rely on the relationship between a gene's mean expression and its variation (variance or dispersion). They assume that technical noise follows a predictable mean-variance trend, and genes that deviate significantly from this trend—possessing more variation than expected—are driven by biological factors [14] [13]. For instance, the "vst" method in Seurat's FindVariableFeatures function directly models this relationship using a polynomial regression or LOESS curve [14].

A more recent and innovative approach, GLP (Genes identified through LOESS with Positive Ratio), challenges the sole reliance on variance. It posits that the positive ratio (the proportion of cells in which a gene is detected) provides a more robust estimator of a gene's importance in the context of scRNA-seq data's characteristic sparsity [7]. Its core assumption is that for a given positive ratio, genes with a mean expression level significantly higher than the local expected value carry more relevant biological information. This method employs a two-step LOESS regression, optimized using the Bayesian Information Criterion to avoid overfitting, to precisely capture this non-linear relationship and identify outlier genes [7].

Comparative Analysis of HVG Selection Methods

Choosing an appropriate HVG selection method requires a nuanced understanding of their performance characteristics, advantages, and limitations. The following table provides a comparative summary based on the evaluated search results.

Table 2: Comparative Analysis of HVG Selection Methods

Method	Key Advantages	Key Limitations	Best-Suited For
Seurat (vst) [14] [13]	Well-established, widely adopted, directly models mean-variance trend.	Assumptions may be violated in very sparse data; performance can be dataset-dependent.	Standard workflows on datasets with moderate sparsity.
Scanpy (seurat) [15] [16]	Integrates well with Python/scverse ecosystem, offers multiple flavors.	Similar limitations to the original Seurat method it mimics.	Scanpy-based analysis pipelines.
Seurat v3 (vst) [15]	Designed for stability, uses variance stabilization, less sensitive to extreme outliers.	Requires pre-selection of n_features.	Large or complex datasets where robust stability is required.
GLP [7]	Specifically designed for sparsity and dropout noise; uses positive ratio; automated bandwidth optimization.	Newer method with less extensive independent benchmarking.	Datasets with high dropout rates where traditional variance-based methods may struggle.
Batch-Aware HVG [15] [17]	Prevents selection of batch-specific genes; lightweight batch correction.	Does not replace full batch integration; tie-breaking strategy can influence final gene list.	Datasets with strong batch effects (e.g., multiple donors, sequencing runs).

Benchmarking studies highlight that the choice of HVG method can significantly impact downstream analysis outcomes, such as clustering accuracy. A large-scale benchmarking preprint notes that performance differences across entire scRNA-seq analysis pipelines are largely driven by the choice of HVGs and the PCA implementation [18]. Furthermore, the novel GLP method has been shown to consistently outperform several leading methods across metrics like Adjusted Rand Index (ARI), Normalized Mutual Information (NMI), and the silhouette coefficient on a variety of real datasets [7].

A critical consideration in modern datasets is handling batch effects. Both Seurat and Scanpy allow for batch-aware HVG selection. When a batch_key is provided, HVGs are selected within each batch separately and then merged. This process prioritizes genes that are variable across multiple batches, effectively acting as a lightweight batch correction method and avoiding the dominance of the final gene list by batch-specific genes [15] [17].

Experimental Protocols for HVG Selection

This section provides detailed, step-by-step protocols for identifying HVGs using two of the most common software ecosystems: Seurat in R and Scanpy in Python.

Protocol 1: HVG Selection with Seurat in R

The following protocol is adapted from the official Seurat guided clustering tutorial [13].

1. Data Preprocessing and Quality Control:

2. Data Normalization:

Note: As an alternative, users are encouraged to explore the SCTransform function, which replaces the need for NormalizeData, FindVariableFeatures, and ScaleData [13].

3. Identification of Highly Variable Features:

The FindVariableFeatures function supports other selection.method arguments, such as "dispersion" (mean.var.plot) [14].

Protocol 2: HVG Selection with Scanpy in Python

This protocol follows the standard Scanpy preprocessing and clustering workflow [16].

1. Data Preprocessing and Quality Control:

2. Data Normalization:

3. Identification of Highly Variable Genes:

The batch_key parameter is crucial for datasets with multiple samples or batches. It selects HVGs within each batch and merges the results, which helps avoid the selection of batch-specific genes and acts as a lightweight batch correction method [15] [16] [17].

Workflow Visualization

The following diagram illustrates the standard pre-processing workflow encompassing HVG selection, leading into dimensionality reduction, a critical step for downstream analysis like clustering.

Diagram 1: Standard scRNA-seq Pre-processing Workflow. The identification of Highly Variable Genes (HVG) is a pivotal step that occurs after normalization and before scaling and PCA.

The logical relationship between the core assumption of the GLP method and its computational procedure can be visualized as follows:

Diagram 2: GLP Method: From Assumption to Result. The GLP algorithm calculates the positive ratio and mean expression for all genes, then uses an optimized LOESS regression to identify outliers based on its core assumption.

The Scientist's Toolkit: Research Reagent Solutions

The following table details key computational tools and resources essential for implementing the HVG selection protocols described in this article.

Table 3: Essential Research Reagents and Computational Tools for HVG Analysis

Item Name	Function / Role in Experiment	Specific Example / Package
Single-Cell Analysis Suite (R)	Comprehensive toolkit for scRNA-seq analysis in R, including QC, normalization, HVG selection, and clustering.	Seurat (e.g., FindVariableFeatures function) [14] [13]
Single-Cell Analysis Suite (Python)	Comprehensive toolkit for scRNA-seq analysis in Python, equivalent to Seurat.	Scanpy (e.g., pp.highly_variable_genes function) [15] [16]
Normalization Reagent	Adjusts raw count data for sequencing depth differences between cells.	Seurat's LogNormalize [13] or Scanpy's normalize_total [16]
HVG Selection Algorithm	Statistical core for identifying genes with high biological variability.	"vst", "mvp" in Seurat [14]; "seurat", "seurat_v3" in Scanpy [15]; GLP [7]
Batch Covariate Data	Metadata column used to group cells for batch-aware HVG selection.	batch_key argument in Scanpy's highly_variable_genes [15] [17]
Visualization Tool	Generates plots to inspect HVG selection results and QC metrics.	Seurat's VariableFeaturePlot [13] or Scanpy's pl.highly_variable_genes [16]

The Critical Link Between HVG Selection and Downstream PCA Efficacy

In single-cell RNA sequencing (scRNA-seq) analysis, feature selection serves as a critical gateway that determines all subsequent analytical outcomes. The process of selecting highly variable genes (HVGs) before applying principal component analysis (PCA) is not merely a computational convenience but a fundamental step that dictates the biological signals captured in downstream analyses. While PCA provides a powerful framework for dimensionality reduction, its efficacy is entirely dependent on the input features selected during HVG screening. Without careful HVG selection, PCA may amplify technical noise over biological signal, leading to misinterpretation of cellular heterogeneity.

Recent studies have demonstrated that inappropriate feature selection can significantly compromise the detection of biologically relevant patterns, particularly when analyzing subtle cellular variations or complex phenotypic outcomes. Research shows that increased gene expression variability itself may be a key biological indicator in certain neurodevelopmental conditions, separate from changes in mean expression levels [19]. This highlights why the initial HVG selection parameters directly influence which aspects of biological heterogeneity become accessible for discovery in the reduced-dimensional space.

This Application Note establishes robust protocols for HVG selection optimized for downstream PCA efficacy, enabling researchers to preserve critical biological signals while minimizing technical artifacts in their scRNA-seq analyses.

Comparative Performance of HVG Selection Methods

The choice of HVG selection method significantly impacts PCA performance and subsequent clustering outcomes. Different algorithms prioritize distinct aspects of transcriptional heterogeneity, making method selection crucial for specific research applications.

Table 1: Comparison of HVG Selection Method Performance Characteristics

Method	Underlying Principle	Strengths	Limitations	Optimal Application Context
GLP [7]	LOESS regression between positive ratio and mean expression	Superior performance in ARI, NMI, and silhouette coefficients; robust to dropout noise	Computational intensity for large datasets	General purpose clustering where cell states are subtle
HRG [20]	Graph-based regional expression patterns	Excellent for detecting spatially restricted patterns; identifies rare cell types	Dependent on initial graph construction	Spatial transcriptomics; rare cell type identification
Standard HVG [21]	Variance stabilization based on mean-expression relationship	Computational efficiency; widely implemented	Marginal improvement over random selection for abundant cell types	Initial exploratory analysis of well-separated cell types
DE-Based [22]	Differential expression across predefined groups	Excellent for cell type identification	Requires prior cell type knowledge; misses novel states	Targeted studies with established markers
PCA-Based [22]	Gene contributions to principal components	Superior variation recovery; captures continuous transitions	Bias toward highly expressed genes	Recovering full transcriptional heterogeneity

Basic feature selection methods perform well at specific objectives but fail to address all analytical needs simultaneously. For instance, PCA-based selection from unscaled data excels at variation recovery, while differentially expressed genes score highest for cell type identification [22]. Surprisingly, for routine cell type identification where cell types are abundant and well-separated, even randomly chosen features can perform nearly as well as algorithmically chosen features when sufficient in number [21]. However, for more challenging tasks like identifying rare cell populations or subtle state differences, both the feature selection method and the number of features selected markedly influence clustering outcomes [21].

Table 2: Quantitative Benchmarking of HVG Methods Across Diverse Tissue Contexts

Method	Cell Type Recovery (ARI)	Rare Cell Detection	State Transition Resolution	Computational Efficiency
GLP	0.89 ± 0.07	Excellent	Superior	Medium
HRG	0.85 ± 0.09	Outstanding	Good	Low
Standard HVG	0.76 ± 0.12	Poor	Moderate	High
DE-Based	0.82 ± 0.08	Good	Poor	Medium
PCA-Based	0.80 ± 0.10	Moderate	Excellent	High

Recent innovations like the GLP method leverage optimized LOESS regression with positive ratio, demonstrating consistent outperformance across multiple benchmark criteria including adjusted rand index (ARI), normalized mutual information (NMI), and silhouette coefficients [7]. This approach identifies biologically informative genes by examining the relationship between positive ratio and average expression level, effectively mitigating the adverse effects of dropout noise and sparsity characteristic of scRNA-seq data.

Disentangling Cell Type and State Through Strategic Feature Selection

A critical challenge in scRNA-seq analysis involves distinguishing between stable cell type identities and transient cell states. Recent research indicates that strategic feature selection can disentangle these overlapping biological programs by prioritizing genes with distinct expression patterns.

When analyzing multi-condition experiments, feature scoring metrics can be specifically designed to quantify "typeness" versus "stateness" in gene expression patterns [23]. Type-specific scores (e.g., moderated F-statistics across clusters) capture genes with stable expression limited to specific subpopulations, while state-specific scores (e.g., percent variance explained by condition) identify genes exhibiting condition-specific responses within subpopulations [23].

Experimental evidence demonstrates that selecting features based on type-not-state prioritization produces embeddings that yield more comparable results between cluster-based and cluster-free differential testing approaches [23]. This separation of biological programs is particularly valuable when studying cellular responses to perturbations, disease progression, or developmental processes where both permanent and transient transcriptional changes coexist.

Figure 1: Strategic Feature Selection Workflow for Disentangling Cell Type and State Effects

Experimental Protocols and Procedures

Protocol 1: GLP-Based HVG Selection for Enhanced PCA

Purpose: To identify highly variable genes using optimized LOESS regression with positive ratio for superior downstream PCA performance.

Materials:

• Raw count matrix from scRNA-seq data
• Computational environment with R or Python installed
• GLP implementation (available from referenced publication)

Procedure:

• Data Preprocessing: Filter genes detected in fewer than 3 cells to minimize sparsity artifacts [7].
• Parameter Calculation: For each gene, compute:
  • Mean expression (λ) = ΣX_ij/c
  • Positive ratio (f) = Σmin(1,Xij)/c where Xij represents expression count of gene j in cell i, and c represents total cell count [7].
• Bandwidth Optimization: Determine optimal LOESS smoothing parameter (⍺) using Bayesian Information Criterion (BIC) across range 0.01-0.1 [7].
• Two-Step LOESS Regression:
  • Initial LOESS fit using Tukey's biweight robust method to identify outlier genes
  • Secondary LOESS fit with outlier weights set to zero
• Feature Selection: Identify genes with expression levels significantly above LOESS-predicted values based on positive ratio
• PCA Application: Perform principal component analysis using selected features

Technical Notes: The adaptive bandwidth selection prevents overfitting while capturing non-linear relationships between positive ratio and mean expression. The two-step regression with robust weighting minimizes influence of technical outliers [7].

Protocol 2: Type-Specific Feature Selection for Multi-Condition Experiments

Purpose: To select features that capture cell type identity while minimizing condition-specific state effects in multi-condition studies.

Materials:

• scRNA-seq data from multiple conditions
• Preliminary cluster assignments
• Condition metadata

Procedure:

• Compute Type Scores: Calculate moderated F-statistics (tF) testing expression differences across high-resolution clusters [23].
• Compute State Scores: Calculate percent variance explained by condition (sPVE) within each cluster [23].
• Rank Combination: Rank genes by difference between type score rank and state score rank (tF-sPVE strategy) [23].
• Feature Selection: Select top-ranked genes based on combined ranking.
• Embedding Generation: Perform PCA on selected features to generate type-informed embedding space.

Validation: Assess embedding quality using Local Inverse Simpson's Index (LISI) to confirm separation of cell types while maintaining condition mixing within types [23].

The Scientist's Toolkit: Essential Research Reagents and Computational Solutions

Table 3: Essential Research Reagents and Computational Tools for HVG Selection

Category	Item/Software	Specification/Function	Application Context
Experimental Platforms	10x Genomics Chromium	Single Cell 3' Library Kit	High-throughput scRNA-seq profiling
Reference Datasets	PBMC 10x Dataset	Peripheral blood mononuclear cells	Method benchmarking and validation
Analysis Packages	Seurat	HVG, SCTransform, PCA	General scRNA-seq analysis pipeline
Specialized Tools	Spapros	Combinatorial probe selection	Targeted spatial transcriptomics design
Benchmarking Suites	ScRNA-seq Evaluation Metrics	ARI, NMI, silhouette coefficients	Method performance quantification
Visualization Tools	UMAP	Non-linear dimensionality reduction	Visualization of PCA outputs

Analytical Framework for Method Selection and Evaluation

Establishing a systematic approach for HVG method selection and evaluation is essential for generating biologically meaningful results. The analytical framework should align method capabilities with specific research objectives.

Figure 2: Decision Framework for HVG Method Selection Based on Research Objectives

For studies focusing on subtle cellular states or developmental transitions, GLP-based selection provides superior resolution due to its robustness to technical noise [7]. When investigating rare cell populations or spatial expression patterns, HRG offers enhanced detection capability by identifying regionally distributed genes [20]. In targeted studies with established marker panels, DE-based selection ensures precise cell type identification, while PCA-based feature selection maximizes capture of continuous transcriptional variation [22].

Evaluation should incorporate multiple metrics assessing both technical performance and biological relevance. Cluster quality metrics (ARI, NMI) quantify agreement with known cell type labels, while neighborhood preservation scores assess maintenance of global data structure [22]. Additionally, biological validation through known marker expression and functional enrichment analysis ensures selected features capture biologically meaningful signals rather than technical artifacts.

The critical link between HVG selection and downstream PCA efficacy underscores the importance of deliberate methodological choices in scRNA-seq analysis. By aligning feature selection strategies with specific biological questions and implementing robust evaluation frameworks, researchers can significantly enhance the biological insights derived from their data. The protocols and analytical frameworks presented here provide a foundation for optimizing this critical step in the single-cell analysis workflow, ensuring that PCA operates on features that faithfully represent the underlying biological heterogeneity of interest.

From Theory to Practice: A Technical Walkthrough of Leading HVG Selection Methods

In single-cell RNA sequencing (scRNA-seq) analysis, the selection of highly variable genes (HVGs) represents a critical preprocessing step that profoundly influences all downstream analyses, including dimensional reduction, clustering, and cell type identification. Technical variation and the high dimensionality of scRNA-seq data necessitate effective strategies to isolate biological heterogeneity from confounding technical artifacts. This protocol examines two principal methodologies within the Seurat ecosystem for HVG selection: the variance stabilizing transformation (VST) algorithm implemented in FindVariableFeatures() and the regularized negative binomial regression approach underlying SCTransform(). We evaluate their performance within the context of a comprehensive workflow that includes data normalization, dimensionality reduction, and dataset integration, providing structured guidance for researchers making informed decisions about their analytical pipelines.

Methodological Comparison: VST vs. SCTransform

The following table summarizes the core characteristics, technical approaches, and optimal use cases for the two primary HVG selection methods.

Table 1: Comparative Analysis of HVG Selection Methods in Seurat

Feature	FindVariableFeatures (VST)	SCTransform
Method Name	Variance Stabilizing Transformation	Regularized Negative Binomial Regression
Underlying Model	Loess regression on mean-dispersion relationship [14]	Regularized Negative Binomial [24] [25]
Input Data	Raw counts from the counts slot [26]	Raw counts
Primary Output	Log-normalized data (RNA assay data slot)	Pearson residuals (SCT assay scale.data slot) [24] [25]
HVG Selection	Identifies 2,000 genes by default [13]	Identifies 3,000 genes by default [24]
Technical Covariates	Requires explicit regression (e.g., vars.to.regress in ScaleData) [13]	Integrated into the model (e.g., vars.to.regress in SCTransform) [24]
Computational Speed	Faster	Slower, but improved in v2 [27]
Downstream Workflow	NormalizeData() -> FindVariableFeatures() -> ScaleData() [13]	Replaces the above three commands in a single step [24]

Experimental Protocols

Protocol 1: Standard Log-Normalization with VST

This protocol utilizes the FindVariableFeatures() function with the VST method, which operates directly on raw counts to identify genes that are highly variable across cells [26]. It is a well-established and computationally efficient approach.

Procedure:

• Data Preprocessing: Filter cells based on quality control metrics (e.g., number of detected genes and mitochondrial percentage).

• Log-Normalization: Normalize the data using a global-scaling method.

• HVG Selection: Identify the top 2,000 highly variable genes using the VST method.

• Data Scaling: Scale the data to give equal weight to all HVGs in downstream analysis. Technical covariates like mitochondrial percentage can be regressed out at this stage.

Protocol 2: Normalization and HVG Selection with SCTransform

This protocol uses SCTransform(), which leverages a regularized negative binomial model to simultaneously normalize data, account for technical confounding, and identify variable features. This method often provides improved biological distinction and more robust downstream results [24] [25].

Procedure:

• Data Preprocessing: Perform the same initial QC filtering as in Protocol 1.
• SCTransform Normalization: Run the SCTransform function, which replaces NormalizeData, FindVariableFeatures, and ScaleData. Technical variables can be directly regressed out during this step.

• Output Handling: The results are stored in a new SCT assay. By default, subsequent analyses use this assay.
  • pbmc[["SCT"]]$scale.data: Contains the Pearson residuals used for PCA and other dimensional reductions [24].
  • pbmc[["SCT"]]$data: Contains log-normalized corrected counts suitable for visualization [24] [25].

Workflow Diagram: VST vs. SCTransform

The diagram below illustrates the key steps and differences between the two protocols.

Integration of Datasets Using Anchors

Following HVG selection and individual dataset preprocessing, integrating multiple samples or datasets is often essential to compare biological conditions and correct for batch effects. Seurat's anchor-based integration is a powerful method for this purpose [28].

Integration Method Selection

The choice of integration method depends on the data characteristics and computational constraints.

Table 2: Overview of Data Integration Methods in Seurat

Method	Underlying Technique	Key Features	Recommended Use Case
CCA (Canonical Correlation Analysis)	Finds shared correlation structures [29]	Robust to large expression shifts; can be slow	Integrating datasets with strong biological differences (e.g., across species) [30]
RPCA (Reciprocal PCA)	Projects datasets into mutual PCA space [30]	Faster, more conservative; reduces over-correction	Large datasets; same platform; many non-overlapping cell types [30]
Harmony	Iterative clustering and correction [29]	Runs on PCA embeddings; efficient	General-purpose batch correction
Scanorama	Mutual Nearest Neighbors (MNN) in Python [29]	Panoramic data integration	Integrating datasets with complex batch effects

Protocol: Integration using RPCA with SCTransform

This protocol details a fast and conservative integration workflow using the RPCA method, which is suitable for datasets where a substantial fraction of cells in one dataset may lack a matching type in the other [30].

Procedure:

• Split Object and Normalize: Split the combined Seurat object by batch or sample and apply SCTransform to each independently.

• Select Features and Run PCA: Identify integration features and run PCA on each dataset using these features.

• Find Anchors and Integrate: Identify integration anchors and perform the data integration.

• Integrated Analysis: Run the standard downstream workflow on the integrated data.

Integration Workflow Diagram

The following diagram outlines the key decision points and steps for a typical data integration pipeline in Seurat.

The Scientist's Toolkit: Essential Research Reagents

This table lists the key software tools and functions essential for implementing the workflows described in this protocol.

Table 3: Key Research Reagents and Computational Tools

Item Name	Function / Purpose	Application in Workflow
Seurat R Package	A comprehensive toolkit for single-cell genomics [13]	Primary environment for all data analysis, normalization, integration, and visualization.
FindVariableFeatures()	Identifies genes with high cell-to-cell variation [14]	HVG selection in the standard pre-processing workflow.
SCTransform()	Performs normalization, variance stabilization, and HVG selection using a regularized negative binomial model [24] [27]	Replaces the standard normalization, HVG selection, and scaling steps in a single command.
IntegrateData()	Combines multiple datasets into a single integrated object using identified anchors [28]	The final step in data integration, creating a new "integrated" assay for joint analysis.
glmGamPoi Package	Accelerates the estimation of parameters for negative binomial models [24] [27]	Significantly improves the speed of SCTransform() when installed and used.

The selection of highly variable genes is a foundational step that dictates the quality of all subsequent analysis in single-cell RNA sequencing. This protocol has detailed two robust methodologies for this task. The VST-based workflow offers a fast, transparent, and well-established approach. In contrast, SCTransform provides a more sophisticated, unified model that often yields superior biological distinction by simultaneously handling normalization, variance stabilization, and technical noise modeling.

For integration, the choice between CCA and RPCA should be guided by the dataset size and biological context. RPCA is recommended for larger datasets or when integrating samples where not all cell types are expected to overlap, as it is faster and more conservative, mitigating over-correction [30].

By following these detailed application notes, researchers can make informed decisions to optimize their scRNA-seq analysis pipeline, ensuring that the identification of highly variable genes provides a solid foundation for revealing meaningful biological insights.

Single-cell RNA sequencing (scRNA-seq) has revolutionized biological research by enabling the characterization of gene expression at unprecedented resolution. However, the raw data generated from scRNA-seq experiments contains technical variations that must be addressed before meaningful biological insights can be extracted. Scanpy, a scalable Python-based toolkit for analyzing single-cell gene expression data, provides comprehensive functionality for preprocessing, visualization, clustering, and trajectory inference [31]. Within this workflow, two critical steps that significantly impact downstream analyses are the selection of highly variable genes (HVGs) and the appropriate transformation of count data. These preprocessing steps are particularly crucial when preparing data for principal component analysis (PCA), as they directly influence which biological signals are preserved and amplified versus which technical noises are suppressed.

The fundamental challenge in scRNA-seq analysis lies in distinguishing biological heterogeneity from technical artifacts. Biological heterogeneity represents meaningful differences in gene expression patterns across cell types and states, while technical artifacts arise from variations in sequencing depth, capture efficiency, and other experimental factors. HVG selection addresses this challenge by identifying genes that exhibit more variability than would be expected due to technical noise alone, thus prioritizing genes likely to inform meaningful biological distinctions [15]. Similarly, logarithmic transformations help stabilize the mean-variance relationship in count data, preventing highly expressed genes from dominating downstream analyses simply due to their abundance rather than their biological information content.

This protocol focuses specifically on Scanpy's implementations of these critical preprocessing steps, with particular attention to their impact on PCA-based dimensional reduction. We provide detailed application notes for researchers, scientists, and drug development professionals working with single-cell transcriptomics data, framed within the context of optimizing feature selection before PCA to enhance biological discovery.

Theoretical Foundation: Highly Variable Genes and Logarithmic Transformations

The Biological and Statistical Rationale for HVG Selection

In scRNA-seq data analysis, the selection of highly variable genes serves as a crucial feature selection step that enhances signal-to-noise ratio by focusing computational resources on the most informative features. The underlying assumption is that genes exhibiting high variability across cells beyond what would be expected from technical noise are more likely to be associated with biologically meaningful processes, such as differentiation, activation, or response to stimuli. From a statistical perspective, this involves modeling the expected relationship between a gene's mean expression and its dispersion (variance-to-mean ratio), then identifying genes that significantly deviate from this expected relationship [15].

Scanpy provides multiple algorithmic approaches for HVG selection through its pp.highly_variable_genes() function, each with distinct statistical foundations [15]. The 'seurat' flavor, inspired by Satija et al. (2015), calculates normalized dispersions by scaling with the mean and standard deviation of dispersions for genes within expression bins. The 'cellranger' flavor, based on Zheng et al. (2017), employs a similar dispersion-based approach but uses different parameter defaults. Meanwhile, the 'seuratv3'/'seuratv3paper' flavor, implementing the method from Stuart et al. (2019), computes a normalized variance for each gene after standardizing the data with a regularized standard deviation. These methods all share the common goal of identifying genes whose variability exceeds technical expectations, though their specific implementations differ in how they model and normalize this relationship.

The Role of Logarithmic Transformation in scRNA-seq Analysis

Logarithmic transformation of count data serves multiple essential purposes in scRNA-seq analysis. First, it helps stabilize variance across the dynamic range of expression values, addressing the mean-variance relationship inherent to count data where highly expressed genes naturally exhibit greater variance. Second, it makes the data more approximately normally distributed, which improves the performance of many statistical methods used in downstream analyses. Third, it reduces the influence of extreme values, preventing a small number of highly expressed genes from dominating dimensional reduction techniques like PCA.

Scanpy implements logarithmic transformation through the pp.log1p() function, which computes (X = \log(X + 1)) [32]. The addition of 1 (pseudocount) prevents undefined values when dealing with zero counts. While this transformation is widely applied, it's important to note that it interacts with prior normalization steps in non-trivial ways. For instance, after normalization to total counts per cell, the relationship between cells is altered by the subsequent log transformation, as the mean of log-counts differs from the log-mean count [33]. This has implications for interpreting distances between cells in reduced-dimensional spaces.

Table 1: Common Data Transformations in Scanpy

Transformation	Function	Formula	Primary Purpose	Considerations
Log1p	scanpy.pp.log1p()	(X = \log(X + 1))	Variance stabilization	Standard approach; assumes base (e) unless specified
Normalize_total	scanpy.pp.normalize_total()	(X{\text{norm}} = X \cdot \frac{\text{targetsum}}{\sum X})	Depth normalization	Often followed by log1p; target_sum defaults to median
Pearson Residuals	scanpy.experimental.pp.highly_variable_genes()	(\frac{X - \mu}{\sqrt{\mu}})	Mean-variance stabilization	Alternative to log1p; does not require prior normalization

Integration of HVG Selection and Logarithmic Transformation in the Analytical Workflow

The sequential application of logarithmic transformation and HVG selection represents a critical preprocessing pipeline that directly influences the quality of subsequent dimensional reduction through PCA. When properly executed, this pipeline enhances biological signal by: (1) reducing technical artifacts through normalization and transformation; (2) focusing on genes most likely to represent biological heterogeneity; and (3) creating a feature space amenable to linear decomposition methods like PCA.

The order of operations is crucial. Typically, normalization (e.g., pp.normalize_total()) is applied first to account for differences in sequencing depth between cells. This is followed by logarithmic transformation (pp.log1p()) to stabilize variance. Finally, HVG selection (pp.highly_variable_genes()) identifies the most informative genes for downstream analysis [16]. However, emerging approaches like analytic Pearson residuals offer an alternative that combines normalization and variance stabilization in a single step [34].

The following diagram illustrates the standard preprocessing workflow in Scanpy and the key decision points researchers encounter:

Figure 1: Standard Scanpy Preprocessing Workflow for PCA. This diagram outlines the sequential steps in preprocessing scRNA-seq data, from raw counts to PCA-ready data, with Scanpy function calls indicated for key steps.

Scanpy Protocols for HVG Selection and Logarithmic Transformation

Quality Control and Data Preparation

Before applying HVG selection or transformations, proper quality control (QC) and data preparation are essential. The following protocol establishes a robust foundation for subsequent analysis:

This QC protocol removes low-quality cells and genes while preserving the raw count data for reference. The specific thresholds should be adjusted based on experimental conditions and cell types, but the provided values represent reasonable starting points for most human tissue scRNA-seq datasets [16] [35].

Normalization and Logarithmic Transformation Protocols

Following quality control, normalization and transformation create comparable expression values across cells:

This protocol uses median normalization rather than fixed-sum normalization (e.g., 10,000 reads per cell), as it is more robust to outliers [16]. The log1p transformation creates a more approximately normal distribution while preserving zero values. The normalized and transformed data is stored in adata.raw to preserve it before subsequent feature selection and scaling.

It's important to recognize that while this normalization and log transformation approach is widely used, it does alter the relationship between cells. As noted in the Scanpy documentation, the mean of log-counts after normalization is not generally the same as the log-mean count, which affects intercellular distances in downstream analyses [33].

Highly Variable Gene Selection Methods

Scanpy provides multiple flavors for HVG selection, each with distinct characteristics:

For datasets with multiple batches, setting the batch_key parameter enables batch-aware HVG selection, which identifies genes that are variable across batches rather than batch-specific genes [15]. This acts as a lightweight batch correction method before more sophisticated integration techniques.

Table 2: Comparison of HVG Selection Methods in Scanpy

Method	Flavor Parameter	Expected Input	Key Parameters	Advantages	Limitations
Seurat	'seurat'	Logarithmized data	min_mean, max_mean, min_disp	Well-established, interpretable	Sensitive to parameter choices
Cell Ranger	'cell_ranger'	Logarithmized data	n_top_genes	Simple interface, consistent number of genes	Less customizable
Seurat v3	'seurat_v3'	Count data	n_top_genes, batch_key	Batch integration, modern standard	Requires additional dependency (scikit-misc)
Pearson Residuals	'pearson_residuals' (experimental)	Count data	theta, clip	No prior normalization needed, analytically exact	Newer, less extensively validated

Alternative Protocol: Analytic Pearson Residuals

For researchers seeking an alternative to the standard normalization-log transformation-HVG selection pipeline, Scanpy's experimental Pearson residuals approach provides a unified solution:

This approach simultaneously normalizes, transforms, and selects features without the need for explicit normalization and log transformation steps [34]. Pearson residuals transform raw UMI counts into a representation that removes technical variation from differences in total counts between cells while stabilizing the mean-variance relationship across genes.

The following diagram compares these two alternative workflows for preparing data for PCA:

Figure 2: Comparison of Standard and Pearson Residuals Workflows. This diagram illustrates two alternative pathways for preprocessing scRNA-seq data before PCA, highlighting the simplified workflow enabled by Pearson residuals.

Research Reagent Solutions and Computational Tools

Table 3: Essential Research Reagents and Computational Tools for scRNA-seq Analysis

Item	Function/Purpose	Example/Format	Considerations
10X Genomics Kit	Single-cell RNA library preparation	3' Gene Expression, 5' Gene Expression, Multiome	Different chemistry versions affect gene detection
Parse Biosciences Evercode	Fixed RNA profiling with combinatorial barcoding	WT Mini, WT Mega	Enables massive scaling without specialized equipment
Cell Ranger	Process raw sequencing data to count matrices	Output: .mtx, .tsv, .csv files	10X-specific; alternative: kallisto, STARsolo
Scanpy	Comprehensive scRNA-seq analysis	Python package	Requires Python environment; part of scverse ecosystem
Harmony	Batch effect integration	Python package (harmonypy)	Accessed via scanpy.external as sce.tl.harmony_integrate
Scikit-misc	Additional statistical functions	Python package	Required for Seurat v3 flavor in highlyvariablegenes()

Application Notes and Troubleshooting

Optimizing Parameters for Specific Biological Contexts

The default parameters in Scanpy's HVG selection functions may not be optimal for all biological contexts. For example, datasets with particularly low or high sequencing depth may require adjustment of the min_mean and max_mean parameters to focus on appropriately expressed genes. Similarly, the n_top_genes parameter should be considered in the context of the overall experimental design - while 2000 is a reasonable default, larger datasets with more complex cellular heterogeneity may benefit from including more variable genes.

When working with specific tissue types, it can be valuable to verify that known cell-type marker genes are included in the selected HVGs. The following diagnostic approach helps assess the quality of HVG selection:

If key marker genes are not being selected as highly variable, consider adjusting the min_mean and max_mean parameters to ensure the expression levels of these genes fall within the considered range [34].

Addressing Common Challenges in Preprocessing

Batch Effects: When working with multiple samples or batches, HVG selection should be performed in a batch-aware manner to identify genes variable across batches rather than batch-specific genes. Setting the batch_key parameter in pp.highly_variable_genes() enables this functionality [15]. For severe batch effects, consider additional integration techniques such as Harmony, which can be accessed through Scanpy's external API [35].

Sparse Data: Extremely sparse datasets (e.g., with many dropouts) may benefit from the Seurat v3 flavor or Pearson residuals approach, as these methods were specifically designed to handle high dropout rates more effectively than dispersion-based methods.

Large Datasets: For datasets exceeding 1 million cells, computational efficiency becomes critical. The sc.pp.highly_variable_genes() function is generally efficient, but for extremely large datasets, consider the sc.experimental.pp.highly_variable_genes() with Pearson residuals, which can be computed in a single pass [34]. Additionally, ensure sufficient computational resources are available - at least 8 threads and 160 GB of RAM for million-cell datasets [35].

Validation of Preprocessing Quality

The quality of HVG selection and transformation should be validated through multiple approaches:

• Visual Inspection: Use sc.pl.highly_variable_genes(adata) to visualize the relationship between mean expression and dispersion, with selected HVGs highlighted. This helps identify whether the selection parameters are appropriately capturing the most variable genes.

• PCA Variance Exploration: Examine the variance ratio plot after PCA (sc.pl.pca_variance_ratio(adata)) to ensure that the selected principal components capture biological rather than technical variation. A sharp drop in explained variance after the first few PCs may indicate that technical artifacts still dominate the signal.

• Biological Validation: Cluster the data after complete processing (sc.tl.leiden()) and examine whether known cell types separate appropriately and express expected marker genes. Poor separation of known distinct cell types may indicate suboptimal HVG selection or transformation.

• Comparison to Ground Truth: When available, compare the results to established classifications or FACS sorting labels to quantify the preservation of biological signal.

By following these application notes and troubleshooting guidelines, researchers can optimize their preprocessing pipeline to maximize biological insights from their scRNA-seq data, particularly when preparing for PCA and subsequent dimensional reduction techniques.

The analysis of single-cell RNA sequencing (scRNA-seq) data presents substantial challenges due to its high dimensionality and sparsity. Effective feature selection to identify highly variable genes (HVGs) is a critical preprocessing step before principal components analysis (PCA) and other downstream analyses. This application note details two innovative computational frameworks—GLP's optimized LOESS regression and scGBM's Poisson bilinear model—that address fundamental limitations in current HVG selection and dimensionality reduction methodologies. By providing detailed protocols and performance comparisons, we equip researchers with advanced tools for enhancing biological discovery in transcriptomic studies.

GLP: Optimized LOESS Regression for Feature Selection

Theoretical Foundation

The GLP (genes identified through LOESS with positive ratio) method introduces a novel approach to feature selection by examining the relationship between gene average expression level (λ) and positive ratio (f). The positive ratio represents the proportion of cells in which a gene is detected, providing a more precise estimator of true population parameters than variance-based measures in sparse scRNA-seq data [6] [7].

The core mathematical assumption derives from Poisson distribution properties, where the relationship between average expression and positive ratio can be expressed as:

$$f = P(X ≥ 1) = 1 - e^{-\lambda}$$

For small values of λ, Taylor approximation yields $f ≈ λ$, though empirical evidence shows this relationship becomes nonlinear across the full gene expression range [6]. GLP leverages this deviation to identify biologically informative genes that demonstrate higher expression levels than expected given their positive ratios.

Algorithmic Implementation

The GLP algorithm implements a sophisticated two-step LOESS (locally estimated scatterplot smoothing) regression procedure with adaptive parameter optimization [6] [7]:

Table 1: Key Algorithmic Parameters in GLP Implementation

Parameter	Default Value	Description	Biological Rationale
Minimum cells	3	Filters genes detected in fewer cells	Ensures reliable positive ratio calculation
α range	0.01-0.1	Smoothing parameter range for LOESS	Balances local flexibility and overfitting
Step size	0.01	Increment for α optimization	Provides granular optimization within computational constraints
Polynomial degree	1 or 2	Local regression complexity	Captures local nonlinearities while maintaining stability

The algorithm employs Bayesian Information Criterion (BIC) for automatic bandwidth selection and incorporates Tukey's biweight robust statistical method to minimize outlier influence [6]. The complete workflow integrates these components into a cohesive feature selection pipeline, illustrated below:

Performance Benchmarks

GLP was rigorously evaluated against eight state-of-the-art feature selection methods across 20 diverse scRNA-seq datasets [6]. Performance was assessed using three established metrics: adjusted rand index (ARI) for clustering accuracy, normalized mutual information (NMI) for cluster similarity, and silhouette coefficient for cluster separation.

Table 2: Comparative Performance of GLP Against Leading Methods

Method	ARI	NMI	Silhouette Coefficient	Computational Efficiency
GLP	0.781 ± 0.042	0.812 ± 0.038	0.745 ± 0.051	Moderate
VST	0.692 ± 0.056	0.734 ± 0.049	0.681 ± 0.048	High
SCTransform	0.705 ± 0.051	0.748 ± 0.045	0.694 ± 0.052	Moderate
M3Drop	0.634 ± 0.063	0.682 ± 0.058	0.627 ± 0.061	High
FEAST	0.673 ± 0.055	0.715 ± 0.052	0.665 ± 0.056	Low
DELVE	0.718 ± 0.048	0.761 ± 0.044	0.703 ± 0.049	Moderate

GLP consistently outperformed competing methods across all evaluation metrics, demonstrating particular strength in preserving biologically relevant gene features for downstream clustering applications [6]. The method showed robust performance across diverse biological contexts, from developmental studies to disease models.

scGBM: Model-Based Dimensionality Reduction

Theoretical Framework

The scGBM (single-cell Poisson bilinear model) approach addresses fundamental limitations in traditional PCA applied to transformed count data. Conventional transformations like log(1+x) or Pearson residuals can induce spurious heterogeneity and mask true biological variability [36]. scGBM directly models UMI counts using a Poisson distribution, providing a more biologically grounded foundation for dimensionality reduction.

The model specification follows:

$$Y{ij} \sim \text{Poisson}(\mu{ij})$$ $$\log(\mu{ij}) = (UV^T){ij}$$

Where $U$ represents cell-level latent factors and $V$ represents gene-level loadings, analogous to traditional PCA but operating in the natural parameter space of the count distribution [36].

Computational Innovations

scGBM introduces three key advancements over existing model-based approaches:

• Fast estimation algorithm: Utilizes iteratively reweighted singular value decompositions to efficiently handle datasets with millions of cells [36]
• Uncertainty quantification: Provides calibrated inference for each cell's latent position
• Cluster cohesion index (CCI): Leverages embedding uncertainties to assess cluster stability and relationships

The integration of these components creates a comprehensive framework for dimensionality reduction that simultaneously addresses computational scalability and statistical rigor:

Empirical Validation

scGBM was evaluated against leading methods including scTransform and GLM-PCA using both simulated and real scRNA-seq datasets [36]. In a simulation with three cell types distinguished by single marker genes, scTransform failed to separate cell types in the principal component space, while scGBM successfully captured the biologically relevant structure.

Table 3: Performance Comparison in Single Marker Gene Simulation

Method	Cell Type Separation	Marker Gene Detection	Runtime (seconds)
scGBM	Clear separation of all three types	Accurate identification of markers	145.2
scTransform	Poor separation, mixed clusters	Masked by noise genes	98.7
GLM-PCA	Moderate separation	Partial identification	310.5
ZINB-WAVE	Good separation	Accurate identification	892.4

The computational efficiency of scGBM represents a significant advancement over other model-based approaches, enabling application to the increasingly large datasets generated by modern single-cell technologies [36].

Integrated Protocol for HVG Selection and Dimensionality Reduction

Comprehensive Experimental Workflow

This protocol integrates GLP and scGBM into a cohesive pipeline for analyzing scRNA-seq data, from raw counts to biological interpretation:

Step-by-Step GLP Implementation

Phase 1: Data Preparation and Positive Ratio Calculation

• Load the count matrix (genes × cells) into R
• Filter genes detected in fewer than 3 cells
• Calculate average expression (λ) for each gene: $λj = \frac{1}{c}∑{i=1}^c X_{ij}$
• Calculate positive ratio (f) for each gene: $fj = \frac{1}{c}∑{i=1}^c \min(1,X_{ij})$

Phase 2: Optimized LOESS Regression

• Define α range from 0.01 to 0.1 with step size 0.01
• Compute BIC for each α value: $BIC = c × \ln(RSS/c) + k × \ln(c)$
• Select α with minimum BIC value
• Perform first LOESS regression with Tukey's biweight function
• Assign zero weights to outlier genes
• Perform second LOESS regression with updated weights
• Select genes with expression levels significantly above the fitted curve

Phase 3: Downstream Analysis

• Use selected HVGs for scGBM dimensionality reduction
• Proceed with clustering and trajectory analysis

Step-by-Step scGBM Implementation

Phase 1: Model Initialization

• Format filtered count matrix with cells as rows and GLP-selected genes as columns
• Initialize latent factors using randomized SVD
• Set convergence tolerance (default: 1e-6) and maximum iterations (default: 100)

Phase 2: Model Fitting

• While (not converged and iterations < maximum):
  • Compute expected counts: $\mu{ij} = \exp((UV^T){ij})$
  • Compute weights: $w{ij} = \mu{ij}$
  • Compute working response: $z{ij} = (UV^T){ij} + (Y{ij} - \mu{ij})/\mu{ij}$
  • Update U and V via weighted SVD of z
  • Check convergence: $\|U^{(t)} - U^{(t-1)}\|F < \text{tolerance}$

Phase 3: Uncertainty Quantification and Cluster Validation

• Compute bootstrap confidence intervals for latent positions
• Calculate Cluster Cohesion Index (CCI) for each cluster
• Filter clusters with CCI below significance threshold
• Visualize stable clusters in low-dimensional space

The Scientist's Toolkit

Table 4: Essential Research Reagent Solutions for Implementation

Resource	Type	Function	Availability
GLP R Package	Software	Feature selection using optimized LOESS	GitHub Repository
scGBM R Package	Software	Poisson bilinear dimensionality reduction	GitHub Repository
Seurat	Software	Integration with standard scRNA-seq workflows	CRAN/Bioconductor
Scanpy	Software	Python-based alternative for full workflow	PyPI
10X Genomics Data	Reference data	Benchmarking and method validation	10X Genomics Website
BUSPA Suite	Quality control	Preprocessing and quality assessment	GitHub Repository

The integration of GLP's optimized LOESS regression for feature selection with scGBM's model-based dimensionality reduction represents a significant advancement in scRNA-seq analysis methodology. This combined approach addresses critical limitations in standard workflows, particularly the tendency of transformation-based PCA to introduce technical artifacts and mask biological signals.

GLP provides a robust foundation for HVG selection by leveraging the relationship between positive ratio and expression level, outperforming eight state-of-the-art methods across diverse datasets. The optimized LOESS regression with BIC-based bandwidth selection ensures adaptive fitting while minimizing overfitting, producing gene subsets that maximize preservation of biological information.

scGBM builds upon this foundation by directly modeling count data using a Poisson bilinear framework, avoiding the biases introduced by conventional transformations. The method's computational efficiency enables application to large-scale datasets, while its uncertainty quantification provides novel insights into cluster stability through the Cluster Cohesion Index.

For researchers investigating complex biological systems with single-cell resolution, this integrated pipeline offers enhanced capability to identify meaningful patterns, characterize cellular heterogeneity, and derive biologically valid conclusions from high-dimensional transcriptomic data.

Step-by-Step Code Implementation in R and Python Environments

High-dimensional genomic data, particularly from single-cell RNA sequencing (scRNA-seq) technologies, presents significant analytical challenges. The process of identifying highly variable genes (HVGs) prior to dimensionality reduction techniques like Principal Component Analysis (PCA) is a critical step in extracting biologically meaningful information from transcriptomic data. This protocol outlines standardized computational procedures for HVG selection and subsequent PCA implementation across both R and Python environments, enabling researchers to efficiently navigate the complexities of modern genomic datasets.

The core challenge in scRNA-seq analysis stems from characteristic data properties including high sparsity and technical artifacts like dropout noise, where expressed genes are mistakenly measured as zero. These properties complicate the distinction between technical variability and true biological heterogeneity. As noted in recent research, "sparsity arises from a complex interplay between true biological expression selectivity and technical dropout artifacts, making it difficult to accurately model using simple statistical priors" [7]. Therefore, robust feature selection methods that can effectively distinguish biologically meaningful genes while mitigating the adverse effects of dropout noise and sparsity are essential for meaningful downstream analysis.

Theoretical Background

Highly Variable Gene Selection Concepts

The fundamental assumption underlying HVG selection is that genes exhibiting high variation across cells beyond technical noise are most likely to represent biologically significant heterogeneity. Traditional HVG methods rely on measuring the relationship between gene expression mean and variance [7]. However, recent methodological advances propose more sophisticated approaches.

The newly developed GLP (genes identified through LOESS with positive ratio) method introduces a novel approach by examining the relationship between positive ratio (the proportion of cells where a gene is detected) and average expression level. This method hypothesizes that "genes with higher expression levels at a given positive ratio are likely to be more significant and better reflect true biological information" [7]. Instead of defining a dropout-rate, GLP utilizes the positive ratio of genes as a more straightforward yet precise estimator, integrating optimized LOESS regression with Bayesian Information Criterion for automatic bandwidth selection to prevent overfitting [7].

Principal Component Analysis Fundamentals

Principal Component Analysis is a dimensionality reduction technique that transforms high-dimensional data into a set of orthogonal components called principal components, which capture maximum variance in the data [37]. The principal components are constructed as linear combinations of the initial variables, with the first component accounting for the largest possible variance, the second component capturing the next highest variance while being uncorrelated with the first, and so on [37].

In biological contexts, PCA simplifies complex datasets while preserving critical patterns and trends. The eigenvalue associated with each eigenvector indicates the amount of variance carried by that principal component, allowing researchers to prioritize components that capture the most biologically relevant information [38]. When applied to genomic data, PCA performed on properly selected HVGs typically reveals more biologically meaningful structures than analysis conducted on the full gene set.

Computational Environment Setup

Research Reagent Solutions

Table 1: Essential Software Tools and Packages

Category	Tool/Package	Purpose/Function	Environment
Programming Languages	R (v4.0+)	Statistical computing and graphics	Base environment
	Python (v3.8+)	General-purpose programming	Base environment
R Packages	Seurat	Single-cell RNA-seq analysis	R
	Bioconductor	Genomic data analysis	R
	limma	Differential expression analysis	R
Python Packages	Scanpy	Single-cell analysis	Python
	Scikit-learn (sklearn)	Machine learning applications	Python
	NumPy	Scientific computing	Python
	Pandas	Data manipulation	Python
	SciPy	Advanced mathematical functions	Python
Analysis Methods	GLP	Feature selection using LOESS	Both
	PCA	Dimensionality reduction	Both
	SCTransform	Normalization and variance stabilization	R

Installation Commands

Experimental Protocols

The following diagram illustrates the complete analytical pipeline for HVG selection and PCA implementation:

Benchmarking Metrics for HVG Selection

Table 2: Evaluation Metrics for HVG Method Performance

Metric	Calculation	Interpretation	Optimal Value
Adjusted Rand Index (ARI)	Measures similarity between clustering results	Higher values indicate better cluster agreement with biological truth	Maximum 1.0
Normalized Mutual Information (NMI)	Information-theoretic measure of clustering similarity	Higher values indicate more biological information captured	Maximum 1.0
Silhouette Coefficient	Measures separation between clusters	Higher values indicate better-defined clusters	Maximum 1.0
Statistical Calibration	Agreement between reported and empirical p-values	Well-calibrated methods maintain false discovery rates	Close to nominal level
Computational Scalability	Execution time relative to dataset size	Methods should scale efficiently with increasing cells/genes	Minimal time increase

Recent benchmarking studies evaluating 14 computational methods across 96 spatial datasets using these six metrics revealed that "SPARK-X outperforms other benchmarked methods and Moran's I achieves a competitive performance" [39]. Similarly, evaluations of HVG selection methods for scRNA-seq data demonstrate that the GLP method "consistently outperforms eight state-of-the-art feature selection methods across three benchmark criteria" including ARI, NMI, and silhouette coefficient [7].

Code Implementation

Highly Variable Gene Selection in R

The following R code implements HVG selection using the GLP method:

Highly Variable Gene Selection in Python

The following Python code implements comparable HVG selection functionality:

Principal Component Analysis Implementation

After HVG selection, implement PCA using the following standardized approaches:

Results Interpretation and Downstream Analysis

Evaluating HVG Selection Effectiveness

The following diagram illustrates the process for evaluating HVG method performance:

Method Selection Guidelines

Table 3: HVG Method Selection Guide Based on Data Characteristics

Data Type	Recommended Method	Rationale	Key Parameters
Standard scRNA-seq	GLP	Robust to dropout noise, captures biological heterogeneity	alpha: 0.01-0.1
Spatial Transcriptomics	SPARK-X	Specifically designed for spatial patterns	Not specified in sources
Large datasets (>10,000 cells)	Moran's I	Computational efficiency while maintaining performance	Not specified in sources
Datasets with batch effects	SCTransform (Seurat)	Integrated normalization and variance stabilization	Not specified in sources
Rapid prototyping	Seurat VST	Fast, standardized implementation	nfeatures: 2000-3000

When evaluating results, consider that comprehensive benchmarking demonstrates the GLP method "consistently outperforms eight state-of-the-art feature selection methods across three benchmark criteria: adjusted rand index (ARI), normalized mutual information (NMI), and the silhouette coefficient" [7]. This performance advantage translates to more biologically meaningful PCA visualizations and improved downstream analysis including clustering and trajectory inference.

Troubleshooting and Optimization

Common challenges in HVG selection and PCA implementation include insufficient HVGs capturing biological variation, technical artifacts dominating principal components, and poor integration between analysis steps. To address these:

• HVG Selection Optimization: Adjust the number of selected HVGs based on dataset size. For larger datasets (>10,000 cells), consider increasing from the default 2000 genes to 3000-5000.

• Parameter Tuning: Utilize the built-in optimization procedures in the GLP method, which "integrates optimized LOESS regression with Bayesian Information Criterion for automatic bandwidth selection" [7].

• Visual Validation: Always visualize PCA results colored by technical covariates (batch, sequencing depth) to identify potential confounding factors.

• Method Combination: Consider combining multiple HVG selection approaches and comparing their concordance as a quality control measure.

The provided code implementations offer researchers standardized, reproducible frameworks for implementing these critical bioinformatics procedures, forming a solid foundation for robust single-cell genomic analysis.

Navigating Common Pitfalls and Advanced Scenarios in HVG Filtering

In single-cell RNA sequencing (scRNA-seq) analysis, feature selection is an indispensable step that significantly influences all downstream interpretations. Reducing the dimensionality from tens of thousands of genes to a manageable set of highly variable genes (HVGs) enhances computational efficiency, mitigates noise, and preserves fundamental biological signals for procedures like clustering and trajectory inference [40] [41]. The central challenge lies not merely in selecting variable genes but in determining the optimal number of HVGs to capture meaningful biological heterogeneity without introducing excessive noise. The performance of this step is context-dependent; while abundant, well-separated cell types can be identified with even random gene sets, discerning subtle cell states or rare populations demands careful, method-driven HVG selection [42]. This protocol provides a structured framework to navigate this critical optimization parameter, ensuring robust and reproducible single-cell analysis.

Key Considerations for Determining the Number of HVGs

There is no universal number of HVGs applicable to all datasets. The optimal quantity depends on several factors related to your data's characteristics and biological questions. The table below summarizes the primary considerations and their impacts on the selection strategy.

Table 1: Key Factors Influencing the Optimal Number of Highly Variable Genes

Factor	Impact on HVG Number	Practical Implication
Dataset Complexity	Increases required HVG number	Datasets with many distinct cell types or continuous differentiation trajectories require more features to capture the full spectrum of heterogeneity [42].
Subtlety of Biological Variation	Increases required HVG number	Identifying rare cell types or fine-grained cell states requires a more sensitive and specific gene set, which can be compromised by using random genes [42].
Sequencing Depth & Quality	Influences reliable HVG number	Noisier datasets or those with high technical variation may benefit from a more conservative selection to avoid highlighting technical artifacts.
Downstream Task	Dictates optimal HVG number	Simple visualization may require fewer genes, while detailed clustering or integration of multiple samples often performs better with more features [10].
Technical Performance	Varies with selection method	The choice of feature selection method itself influences the optimal number, as some algorithms are more precise or reproducible than others [41].

Current Methods for HVG Selection and Their Performance

Numerous algorithms have been developed for HVG selection, each with different strengths and weaknesses. A key finding from recent benchmarks is that the choice of method can be as important as the number of genes selected.

• Variance-Stabilizing Transformation (VST): Models the relationship between mean expression and variance using local polynomial regression, ranking genes by their standardized variance [40] [41].
• SCTransform: Uses regularized negative binomial regression to account for technical noise and selects genes based on Pearson residuals [41].
• Highly Variable Genes (HVG) in Scanpy: Bins genes by mean expression and ranks them by normalized dispersion (z-score of Fano factor) [42].
• GLP: A recently developed method that uses an optimized LOESS regression between a gene's average expression level and its "positive ratio" (the proportion of cells where the gene is detected) to identify biological signals while minimizing overfitting [7].
• HRG (Highly Regional Genes): A graph-based approach that identifies genes showing regional expression patterns in a cell-cell similarity network, mimicking manual feature selection from 2D visualizations [20].
• BigSur: A statistically principled method based on an analytical model of un-normalized scRNA-seq data, designed to improve the identification of subtle cell states and rare populations [42].

Comparative Performance and Reproducibility

Evaluations reveal significant differences in the performance and reproducibility of these methods. A comprehensive benchmark of over 20 feature selection methods for data integration and query mapping confirmed that HVG selection generally produces high-quality integrations, but the specific implementation matters [10]. Another study comparing nine methods found that SCHS demonstrated superior reproducibility in identifying HVGs across random subsamples of cells, though it exhibited a preference for highly expressed genes [41]. Methods like M3Drop showed lower reproducibility (50-70% overlap in selected genes across tests), whereas SCHS maintained consistently high reproducibility above 90% [41].

Table 2: Performance Comparison of Selected HVG Selection Methods

Method	Core Principle	Reproducibility	Notable Characteristics
SCHS	Detects genes by the spatial distribution of cells [41].	High (>90%) [41]	Prefers highly expressed genes; high reproducibility [41].
Seurat (VST, SCT, Disp)	Models mean-variance relationship or uses Pearson residuals [41].	Medium (50-70%) [41]	A quarter of selected genes are often lowly expressed [41].
GLP	LOESS regression on positive ratio vs. expression [7].	N/A	Designed to be robust to sparsity and dropout noise [7].
Scran	Uses variance of log-expression values [41].	Low to Medium (80-90%) [41]	Selects almost no lowly expressed genes [41].
M3Drop	Leverages dropout rates modeled with Michaelis-Menten function [41].	Low (50-70%) [41]	Can identify lowly expressed genes; lower distinguishing capability in clustering [41].

Recommended Protocols for HVG Number Optimization

A General Workflow for Robust HVG Selection

The following diagram outlines a systematic workflow for determining the number of HVGs, integrating robustness checks and biological validation.

Diagram 1: A systematic workflow for determining the number of HVGs.

Protocol 1: Benchmarking with the SIEVE Strategy

To enhance the robustness of any HVG method, the SIEVE (SIngle-cEll Variable gEnes) strategy can be employed. This approach uses multiple rounds of random sampling to identify a stable set of variable genes, substantially improving cell classification accuracy, especially for methods with lower inherent reproducibility [41].

• Input: A normalized scRNA-seq expression matrix (cells × genes).
• Resampling: Randomly sample (without replacement) 70% of the total cells to create a reference set. Repeat this process at least 50 times.
• HVG Application: In each resampling iteration, apply your chosen HVG selection method (e.g., VST, SCHS) to the reference set to select a pre-defined number of top genes (e.g., 2,000).
• Gene Frequency Calculation: For each gene, calculate the frequency at which it is selected as an HVG across all resampling iterations.
• Final Selection: Construct the final robust HVG set by including genes that appear in a high proportion (e.g., >80%) of the resampling runs. The final number of HVGs is the count of genes passing this frequency threshold.

Protocol 2: Empirical Optimization via Downstream Stability

This protocol uses the stability of downstream results to guide the selection of the HVG number.

• Generate a Range: Using your chosen method, generate ranked lists of HVGs. Create a series of gene sets of increasing size (e.g., 500, 1,000, 1,500, 2,000, 2,500, 3,000).
• Perform Downstream Analysis: For each gene set size, perform the core downstream analysis (e.g., clustering using the same resolution parameters, or data integration).
• Evaluate Stability: Calculate clustering metrics like Adjusted Rand Index (ARI) or Normalized Mutual Information (NMI) between consecutive gene set sizes. Alternatively, monitor the stability of known cell type markers and biological patterns in visualizations.
• Identify the Elbow Point: The optimal number is often found at the point where increasing the number of HVGs no longer leads to substantial improvements in cluster stability or biological coherence. This is the point of diminishing returns.

The Scientist's Toolkit: Essential Reagents & Computational Tools

Table 3: Key Computational Tools for HVG Selection and Analysis

Tool / Resource	Function	Implementation
Seurat	A comprehensive R toolkit for single-cell analysis. Includes VST, SCTransform, and dispersion-based methods for HVG selection [41].	R
Scanpy	A scalable Python toolkit for analyzing single-cell data. Includes the "pp.highlyvariablegenes" function for HVG selection [10] [42].	Python
SIEVE	A strategy, not a standalone tool, that can be implemented to improve the robustness of existing HVG selection methods via resampling [41].	R/Python
GLP	A feature selection algorithm based on optimized LOESS regression with positive ratio [7].	R
HRG	A graph-based feature selection method for identifying regionally distributed genes in a cell-cell network [20].	R/Python
BigSur	A statistically principled method for feature selection and network inference from unnormalized count data [42].	R

Determining the optimal number of highly variable genes is a critical, dataset-specific parameter in scRNA-seq analysis. There is no single magic number; the ideal value balances the capture of meaningful biological heterogeneity against the introduction of noise. As benchmarks show that method choice significantly impacts performance and reproducibility, researchers should adopt a systematic approach. This involves selecting a modern, robust HVG selection method and employing empirical validation—through resampling strategies like SIEVE or downstream stability checks—to identify the gene set size that best stabilizes the biological structures of interest. By following these structured protocols, researchers can ensure their feature selection step enhances, rather than hinders, their scientific discovery.

Addressing Batch Effects and Integrating Multi-Sample Data

In the era of large-scale genomic and single-cell RNA sequencing (scRNA-seq) studies, researchers frequently combine data from multiple samples, experiments, and technologies. This integration is essential for constructing comprehensive cell atlases and achieving sufficient statistical power in biomarker discovery [43]. However, technical variations between datasets—known as batch effects—can obscure biological signals and compromise downstream analyses. Effectively addressing these artifacts is a critical prerequisite for meaningful biological interpretation. The strategic filtering of highly variable genes (HVGs) before applying principal component analysis (PCA) serves as a foundational step in this process, enhancing the signal-to-noise ratio for more robust data integration [10].

Various computational strategies have been developed to correct for batch effects and integrate multi-sample data, ranging from classical statistical approaches to modern deep learning models. The table below summarizes the key methods, their underlying principles, and their applicability to different data scenarios.

Table 1: Data Integration Methods for Multi-Sample Omics Data

Method Name	Underlying Approach	Key Features	Typical Applications
ComBat / limma [43]	Linear model-based adjustment	Effective for known, linear batch effects; fast execution.	Bulk transcriptomics, proteomics data integration.
Batch-Effect Reduction Trees (BERT) [43]	Tree-based hierarchical integration of ComBat/limma	Handles incomplete omic profiles (missing values); considers covariates.	Large-scale integration of proteomics, transcriptomics, metabolomics.
Conditional VAE (cVAE) [44] [45]	Deep generative model (Variational Autoencoder)	Corrects non-linear batch effects; scalable to large datasets.	Single-cell atlases, integrating datasets with complex technical differences.
sysVI (VAMP + CYC) [44]	Enhanced cVAE with VampPrior and cycle-consistency	Improves integration for substantial batch effects (e.g., cross-species).	Integrating datasets from different biological systems (e.g., organoid vs. tissue).
Canonical Correlation Analysis (CCA) [45]	Correlation-based linear projection	Finds correlated components across datasets; interpretable.	Multi-omics integration on matched samples.
Integrative NMF (iNMF) [45]	Matrix factorization	Decomposes data into shared and dataset-specific factors.	Joint clustering of multi-omics data; single-cell data integration (LIGER).

Experimental Protocols

Standard Protocol for scRNA-seq Integration with HVG Filtering

This protocol outlines a standard workflow for integrating multiple scRNA-seq datasets, highlighting the critical role of feature selection prior to dimensionality reduction.

I. Preprocessing and Feature Selection

• Data Input: Start with a list of raw or normalized count matrices (cells x genes) from all batches to be integrated.
• Quality Control: Filter out low-quality cells based on metrics like mitochondrial gene percentage, total counts, and number of detected genes. Remove genes expressed in only a minimal number of cells (e.g., < 3 cells) [7].
• Normalization: Normalize the count data for each cell to correct for differences in sequencing depth. A common approach is to use library size normalization.
• Highly Variable Gene (HVG) Selection: This is a critical step for enhancing the signal in downstream PCA.
  • Method: Apply a robust feature selection algorithm. The GLP method is a recent advancement that uses optimized LOESS regression on the relationship between a gene's average expression level (λ) and its positive ratio (f), which is the proportion of cells where the gene is detected [7].
  • Calculation:
    • For each gene j, calculate:
      • Average expression: λ_j = (1/c) * Σ X_ij (sum over all cells i)
      • Positive ratio: f_j = (1/c) * Σ min(1, X_ij) (sum over all cells i)
  • Selection: Genes with expression levels significantly higher than the value predicted by the LOESS regression curve (given their positive ratio) are selected as HVGs [7]. This identifies genes whose variability is likely biologically meaningful.
  • Number of Genes: The number of HVGs to select can impact performance. Benchmarks suggest that using 2000-3000 HVGs is often effective, but this should be validated for specific datasets [10].

II. Data Integration and Downstream Analysis

• Dimensionality Reduction: Apply PCA to the scaled and normalized expression matrix containing only the selected HVGs. This projects the high-dimensional data into a lower-dimensional space that captures the main axes of biological variation.
• Batch Effect Correction: Choose an integration method from Table 1 (e.g., a cVAE-based method like sysVI for complex batches, or BERT for data with many missing values) and apply it using the PCA-reduced representation or the HVG expression matrix as input.
• Visualization and Clustering: Project the integrated data into 2D using UMAP or t-SNE for visualization. Perform clustering on the integrated latent space to identify cell populations.
• Evaluation: Assess integration quality using metrics like:
  • Batch ASW: Average silhouette width of batches; values closer to 0 indicate good mixing [43] [10].
  • iLISI: Integration local inverse Simpson's index; measures diversity of batches in local neighborhoods [44].
  • Bio-conservation: Use metrics like normalized mutual information (NMI) or cLISI to ensure biological variation (e.g., cell types) is preserved [44] [10].

Figure 1: Standard workflow for scRNA-seq data integration, highlighting the pivotal role of HVG selection before PCA.

Advanced Protocol for Complex System Integration

For challenging integration tasks with substantial batch effects—such as merging data from different species, sequencing technologies (e.g., single-cell vs. single-nuclei), or model systems (e.g., organoids vs. primary tissue)—standard methods may fail. The sysVI protocol is designed for these scenarios.

I. System-Specific Data Preparation

• Dataset Assembly: Collect scRNA-seq datasets from the different systems to be integrated (e.g., human and mouse pancreatic islets).
• Independent Preprocessing: Perform quality control and normalization on each dataset independently. It is crucial to use consistent parameters across datasets where possible.
• Confirm Substantial Batch Effects: Quantify the strength of the batch effect. Calculate the distance between cell types within the same system and compare it to the distance between systems. A significant difference confirms the need for a robust method like sysVI [44].

II. Integration with sysVI

• Model Configuration: Use a conditional Variational Autoencoder (cVAE) architecture extended with two key components:
  • VampPrior (VAMP): A multimodal prior for the latent space that helps preserve biological heterogeneity [44].
  • Cycle-Consistency Loss (CYC): A constraint that ensures a cell's representation can be mapped back to its original expression profile, improving the stability of integration [44].
• Training: Train the sysVI model on the concatenated datasets from all systems. The model will learn a joint latent representation where technical differences are minimized, but biological variance is retained.
• Latent Space Extraction: Use the trained model to project all cells into the integrated latent space for downstream analysis.

III. Evaluation of Complex Integration

• Assess Cell Type Alignment: Check if homologous cell types from different systems (e.g., human and mouse beta cells) co-localize in the integrated latent space.
• Check for Over-Correction: Ensure that the model does not "over-integrate" and mix biologically distinct cell types, a common failure mode of adversarial learning methods [44]. The use of VampPrior and cycle-consistency in sysVI is specifically designed to mitigate this.

Figure 2: The sysVI model architecture for integrating datasets with substantial biological and technical differences.

The Scientist's Toolkit

This section details key computational tools and reagents essential for implementing the protocols described above.

Table 2: Essential Reagents and Computational Tools for Data Integration

Item Name	Type	Function / Application	Key Characteristic
GLP Algorithm [7]	Feature Selection Method	Identifies Highly Variable Genes (HVGs) by modeling the gene average expression vs. positive ratio relationship with LOESS regression.	Robust to dropout noise; uses BIC for automatic parameter tuning.
BERT [43]	Data Integration Algorithm	Integrates incomplete omic profiles by building a binary tree of batch-effect correction steps.	Handles missing data without imputation; considers covariates.
sysVI [44]	Data Integration Algorithm	A cVAE-based method for integrating datasets with substantial batch effects (cross-system).	Uses VampPrior and cycle-consistency to preserve biological signals.
ComBat [43]	Data Integration Algorithm	Linear model-based adjustment for removing batch effects.	Established, fast method; forms the core of BERT's pairwise corrections.
VampPrior [44]	Model Component (in sysVI)	A multimodal prior in a VAE that helps the model capture diverse biological cell states.	Preovers loss of biological heterogeneity during integration.
scRNA-seq Data	Biological Data Input	The primary data type for protocols, measuring transcriptomes of individual cells.	Inherently sparse and noisy; requires specialized preprocessing [7].
Batch ASW / iLISI [44] [10]	Quality Metric	Quantitative scores to evaluate the success of batch effect removal and biological conservation.	Essential for objective benchmarking of integration performance.

Quantitative Data and Performance Metrics

The performance of different integration methods can be quantitatively evaluated using established metrics. The choice of feature selection method and the number of HVGs selected significantly impact these outcomes [10].

Table 3: Performance Comparison of Integration Workflows and Metrics

Metric Category	Metric Name	Description	Ideal Value	Notes on Feature Selection Impact
Batch Removal	Batch ASW [43] [10]	Average silhouette width of batches; measures batch mixing.	Closer to 0	HVG selection improves this by removing uninformative noise.
	iLISI [44] [10]	Integration Local Inverse Simpson’s Index; diversity of batches in a cell's neighborhood.	Closer to 1 (high diversity)	Sensitive to the biological signal preserved by feature selection.
Bio-conservation	NMI / bNMI [44] [10]	(Batch-balanced) Normalized Mutual Information; measures clustering similarity to cell type labels.	Closer to 1	Using too few HVGs can degrade this score by losing biological signal [10].
	cLISI [10]	Cell-type Local Inverse Simpson’s Index; purity of cell type labels in a cell's neighborhood.	Closer to 1	Depends on HVGs capturing true cell-type markers.
Data Retention	Numeric Value Retention [43]	Proportion of original data values retained after integration.	Closer to 1	Methods like BERT retain more values than dissection-based approaches.

Benchmarking studies reinforce that using HVGs for feature selection is effective for producing high-quality integrations [10]. The performance of the integration is not solely dependent on the integration algorithm itself but is also a function of the preparatory steps, with HVG filtering being a critical factor that sets the stage for successful batch effect correction and biological discovery.

Preserving Rare Cell Types and Subtle Biological Signals

The emergence of single-cell RNA sequencing (scRNA-seq) has fundamentally transformed biological research, enabling the exploration of cellular heterogeneity at unprecedented resolution [46]. A central challenge in this field is the accurate preservation of rare cell types and subtle biological signals during computational analysis, particularly when performing dimensionality reduction techniques like Principal Component Analysis (PCA). Rare cell populations, such as circulating tumor cells (CTCs), stem cells, and infected cells, often constitute less than 0.1% of a sample yet hold critical information about disease mechanisms, developmental processes, and therapeutic responses [47]. The high dimensionality, technical noise, and sparsity inherent to scRNA-seq data can easily obscure these fragile biological signals. This application note explores integrated experimental and computational strategies to safeguard these valuable signals throughout the single-cell analysis workflow, with particular emphasis on the crucial step of feature selection prior to PCA within the context of advanced thesis research.

The Critical Challenge of Rare Cells in Single-Cell Analysis

Rare cells are defined as populations occurring in frequencies of fewer than 1,000 cells per milliliter of sample and present unique challenges across both wet-lab and computational domains [47]. Their scarcity and fragility mean that conventional analytical approaches often lead to their loss or mischaracterization.

Experimental Obstacles in Rare Cell Isolation

From an experimental perspective, isolating rare cells requires exceptionally gentle and precise methods to maintain cell viability and integrity. Harsh mechanical or chemical treatments can easily compromise these fragile populations. Furthermore, rare cells frequently exhibit substantial heterogeneity and may not uniformly express canonical marker genes, making them difficult to enrich using standard affinity-based approaches [47]. Techniques like Buoyancy Activated Cell Sorting (BACS) have emerged as promising solutions, using microbubbles to gently separate target cells based on gravity rather than more stressful magnetic or column-based methods [47].

Computational Hurdles in Signal Preservation

Computationally, the high sparsity and dropout noise characteristic of scRNA-seq data pose significant challenges for distinguishing true biological variation from technical artifacts [7]. When performing dimensionality reduction via PCA, the dominant sources of variation from abundant cell types can easily overwhelm the subtle signals originating from rare populations. This problem is exacerbated when using all genes for PCA, as the high-dimensional noise can obscure biologically meaningful patterns. Consequently, judicious feature selection before PCA becomes paramount for preserving signals from rare cell types and enabling their subsequent discovery and characterization.

Feature Selection Methods: Preserving Biological Signals Before PCA

Feature selection serves as a critical preprocessing step that directly influences the performance of downstream PCA and the preservation of rare cell signals. The table below summarizes key feature selection methods and their impact on biological signal preservation.

Table 1: Feature Selection Methods for Preserving Rare Cell Signals

Method	Underlying Principle	Advantages for Rare Cells	Limitations
GLP (LOESS with Positive Ratio)	Models relationship between gene average expression and positive ratio using optimized LOESS regression [7]	Precisely captures true biological signal by leveraging positive ratio; minimizes overfitting; outperforms 8 other methods in benchmarks [7]	Requires sufficient cells for stable positive ratio estimation
Highly Variable Genes (HVG) Selection	Identifies genes with higher-than-expected variance across cells [10]	Standardized, widely adopted; effective for general population structure	May miss rare population-specific markers; sensitive to technical noise
RMT-guided Sparse PCA	Uses Random Matrix Theory to denoise covariance matrix estimation; guides sparsity parameter selection [9]	Provides mathematical foundation for noise/signal separation; improves cell type classification accuracy [9]	Computationally intensive; requires large cell numbers for RMT assumptions
Batch-Aware Feature Selection	Accounts for technical batch effects during feature selection [10]	Prevents batch-confounded genes from dominating PCA; preserves biologically relevant variation	Requires careful parameter tuning; may be implementation-specific

The performance of these feature selection methods directly impacts the quality of subsequent data integration and rare cell detection. Benchmarking studies have demonstrated that proper feature selection significantly improves metrics relevant to rare cell preservation, including:

• Cell-type Local Inverse Simpson's Index (cLISI): Measures cell-type mixing quality [10]
• Isolated Label Scores: Evaluate preservation of small populations [10]
• Mapping Accuracy: Assesses ability to correctly project query datasets containing rare cells [10]

Table 2: Benchmark Performance of Feature Selection Methods on Rare Cell Metrics

Method	cLISI Score	Isolated Label F1	Query Mapping Accuracy	Batch Correction
GLP	0.89	0.85	0.91	0.82
HVG (2000 genes)	0.76	0.72	0.84	0.79
RMT-sPCA	0.87	0.83	0.89	0.85
Random Features	0.45	0.38	0.52	0.41

Experimental Protocols for Rare Cell Analysis

Protocol 1: Gentle Rare Cell Isolation using Buoyancy Activated Cell Sorting (BACS)

Principle: Leverages microbubbles conjugated to specific antibodies that gently float target cells to the surface for aspiration-free collection [47].

Reagents Required:

• BACS microbubble cocktail (target-specific)
• Appropriate cell suspension sample
• Collection buffers compatible with downstream applications

Procedure:

• Sample Preparation: Prepare single-cell suspension using gentle dissociation protocols to preserve cell viability and surface epitopes.
• Incubation: Mix BACS microbubbles with cell suspension and incubate for 15-30 minutes with gentle agitation.
• Separation: Allow microbubble-cell complexes to float to surface by gravity (5-10 minutes).
• Collection: Carefully remove supernatant containing untouched cells of interest from the bottom of the tube.
• Wash: Perform gentle washing steps if required for downstream applications.
• Validation: Assess purity, viability, and recovery using appropriate methods (flow cytometry, microscopy).

Downstream Applications: Isolated cells are suitable for scRNA-seq, functional assays, and culture due to maintained viability and physiology [47].

Protocol 2: Integrated Computational Pipeline for Rare Cell Signal Preservation

Principle: Implements a feature selection and dimensionality reduction workflow optimized for rare cell preservation.

Software Requirements:

• Python/R environments with scRNA-seq analysis packages (Scanpy, Seurat)
• GLP implementation or alternative feature selection tools
• PCA and visualization libraries

Procedure:

• Quality Control: Filter low-quality cells and genes while preserving potential rare populations.
• Normalization: Apply appropriate normalization (e.g., SCTransform) to address technical variation.
• Feature Selection: Implement GLP feature selection algorithm:
  • Calculate average expression (λ) and positive ratio (f) for each gene [7]
  • Apply optimized LOESS regression to model λ ~ f relationship [7]
  • Select genes with expression significantly higher than expected based on positive ratio
• PCA on Selected Features: Perform dimensionality reduction using selected gene set.
• Rare Population Identification: Apply clustering algorithms sensitive to small populations (e.g., Leiden clustering with high resolution).
• Validation: Use marker gene expression, differential expression, and trajectory inference to validate rare populations.

Key Parameters:

• GLP bandwidth selection via Bayesian Information Criterion [7]
• PCA dimensions: 20-50 based on elbow plot of explained variance
• Clustering resolution: 2.0-5.0 for rare population detection

Visualization of Integrated Workflow

The following diagram illustrates the complete integrated workflow for preserving rare cell types from sample preparation through computational analysis:

Integrated Workflow for Rare Cell Analysis

Table 3: Key Research Reagent Solutions for Rare Cell Studies

Category	Product/Resource	Specific Function	Application Context
Cell Isolation	BACS Microbubble Kits	Gentle antibody-based cell separation using buoyancy principle [47]	Rare immune cell isolation; CTC enrichment
Cell Isolation	Microfluidic Chips (ChipShop)	Microchannel structures for cell focusing and imaging [48]	Label-free rare cell detection; imaging flow cytometry
scRNA-seq Protocols	10x Genomics Chromium	Droplet-based single-cell partitioning with UMIs [46]	High-throughput rare cell transcriptomics
scRNA-seq Protocols	Smart-Seq2	Full-length transcript coverage for low-input samples [46]	Deep characterization of rare populations
Computational Tools	GLP Feature Selection	Identifies biologically informative genes using positive ratio [7]	Preserving rare cell signals before PCA
Computational Tools	RMT-guided Sparse PCA	Denoised PCA implementation using Random Matrix Theory [9]	Improved cell type classification accuracy
Analysis Platforms	Seurat/Scanpy	Comprehensive scRNA-seq analysis ecosystems [10]	End-to-end rare cell analysis workflows

Troubleshooting and Optimization Strategies

Addressing Common Experimental Challenges

• Low Rare Cell Recovery: Optimize antibody concentrations during BACS; validate with spike-in controls; minimize processing steps.
• Poor Cell Viability: Use gentler dissociation protocols; reduce processing time; incorporate viability-preserving buffers.
• Insufficient Sequencing Depth: Increase read depth to 50,000-100,000 reads per cell for adequate rare population transcript capture.

Computational Parameter Optimization

• Feature Selection Quantity: Benchmark with 1,000-5,000 features to balance signal preservation and noise reduction [10].
• PCA Dimensions: Determine optimal dimensions using elbow plots and rare population stability metrics.
• Cluster Resolution: Use hierarchical or multi-resolution clustering approaches to ensure rare population detection without excessive fragmentation.

Preserving rare cell types and subtle biological signals requires an integrated approach spanning careful experimental design, gentle isolation methodologies, and computationally sophisticated feature selection strategies. The critical step of filtering highly variable genes before PCA directly influences our ability to detect and characterize these biologically significant but technically elusive populations. By implementing the protocols and recommendations outlined in this application note, researchers can significantly enhance their capability to uncover novel biology hidden within rare cellular compartments, ultimately advancing drug discovery, disease mechanism elucidation, and therapeutic development.

Mitigating Overfitting and Model-Specific Biases

In the analysis of high-dimensional biological data, such as single-cell RNA-sequencing (scRNA-seq), the initial step of feature selection—identifying highly variable genes (HVGs) for dimensionality reduction via Principal Component Analysis (PCA)—is critical. The performance and biological relevance of all subsequent analyses depend on this step. However, this process is vulnerable to two significant challenges: overfitting, where a model learns noise and technical artifacts instead of true biological signal, and model-specific biases, which are systematic errors introduced by the algorithms or data themselves. Left unaddressed, these issues can compromise the reliability of cell type classification, trajectory inference, and differential expression analysis. This document outlines the core concepts of these challenges and provides detailed application notes and protocols to mitigate them, ensuring robust and reproducible research outcomes.

Core Concepts and Challenges

The Problem of Overfitting in Feature Selection

In machine learning, overfitting occurs when a model learns the training data too well, including its noise and random fluctuations, resulting in poor performance on new, unseen data [49] [50]. In the context of HVG selection, this translates to selecting genes whose apparent variability is driven by technical artifacts (e.g., dropout events, amplification bias) rather than genuine biological heterogeneity [7].

• Consequences: An overfit feature selection model will lack generalizability. It may perform well on the initial dataset but fail to identify biologically meaningful genes in new data, leading to reduced predictive power and unreliable downstream analyses [49] [51].
• Causes: Primary causes include model complexity that is too high for the available data, insufficient training data (a small number of cells), and noisy data [49] [52] [50].

Understanding Model-Specific Biases

Model-specific bias refers to a systematic deviation in a model's predictions caused by flaws in the algorithm design, feature selection, or underlying data [53]. Unlike overfitting, which generally degrades overall performance, bias can lead to systematically unfair or inaccurate outcomes for specific subgroups within the data.

• Position Bias: Recently identified in Large Language Models, this is a tendency to overemphasize information at the beginning and end of an input sequence while neglecting the middle [54]. While demonstrated in transformers, analogous biases can exist in genomic pipelines based on how data is ordered or processed.
• Selection and Exclusion Bias: This occurs when the data used to train or develop a model does not accurately represent the target population [53] [55]. In scRNA-seq, this could manifest as a protocol that underrepresents certain, perhaps more fragile, cell types.
• Measurement Bias: Arises from systematic errors in data collection, leading to inaccurate or unrepresentative data [53]. In genomics, this could stem from batch effects or non-uniform sequencing depth across samples.

Table 1: Common Biases in Analytical Models and Their Potential Manifestations in Genomics

Bias Type	Definition	Potential Manifestation in Genomics
Position Bias [54]	Systematic overemphasis on early and late parts of an input sequence.	A model's performance varies based on the order of genes or cells in the input matrix.
Selection Bias [53]	Training data is not representative of the target population.	Certain cell types are systematically under-represented in the dataset due to experimental protocol.
Exclusion Bias [53]	Systematic exclusion of specific data segments during acquisition or processing.	Filtering steps inadvertently remove genes expressed in a rare but biologically important cell population.
Measurement Bias [53]	Systematic errors introduced during data collection.	Batch effects or inconsistent library preparation across samples skew gene expression measurements.

The following table summarizes key metrics and outcomes from recent studies that developed methods to address overfitting and bias in feature selection and dimensionality reduction.

Table 2: Quantitative Performance of Advanced Methods Mitigating Overfitting and Bias

Method (Citation)	Core Approach	Benchmark Criteria	Key Performance Outcome
GLP [7]	Optimized LOESS regression for HVG selection using the relationship between average expression and positive ratio, with BIC for bandwidth selection.	Adjusted Rand Index (ARI), Normalized Mutual Information (NMI), Silhouette Coefficient	Consistently outperformed 8 state-of-the-art feature selection methods across 20 scRNA-seq datasets.
RMT-guided Sparse PCA [9]	Uses Random Matrix Theory and biwhitening to guide parameter-free inference of sparse principal components.	Cell-type classification accuracy	Consistently outperformed PCA-, autoencoder-, and diffusion-based methods in cell-type classification tasks across 7 scRNA-seq technologies.
Theoretical Framework for Position Bias [54]	A graph-based framework to diagnose position bias in transformer models, linking it to attention masking.	Information retrieval accuracy	Identified a "lost-in-the-middle" U-shaped performance pattern, with high accuracy for data at sequence ends and low accuracy in the middle.

Detailed Experimental Protocols

Protocol 1: Robust Highly Variable Gene Selection with GLP

This protocol implements the GLP (LOESS with Positive Ratio) algorithm to minimize overfitting during HVG selection [7].

I. Research Reagent Solutions

• Computing Environment: R or Python with necessary statistical libraries.
• Input Data: A pre-processed scRNA-seq count matrix (cells x genes).
• Key Software Functions: loess regression function, Bayesian Information Criterion (BIC) calculation.

II. Step-by-Step Workflow

• Data Preprocessing:
  • Input the raw count matrix.
  • Filter out genes detected in fewer than 3 cells to reduce noise [7].

• Feature Metric Calculation:

  • For each gene ( j ), compute its average expression level ( \lambdaj ): ( \lambdaj = \frac{1}{c} \sum{i=1}^{c} X{ij} ) where ( c ) is the total number of cells and ( X_{ij} ) is the expression count of gene ( j ) in cell ( i ) [7].
  • For each gene ( j ), compute its positive ratio ( fj ): ( fj = \frac{1}{c} \sum{i=1}^{c} \min(1, X{ij}) ) This represents the fraction of cells where the gene is detected [7].

• Model Optimization:

  • Set the LOESS smoothing parameter ( \alpha ) to a range from 0.01 to 0.1 (in steps of 0.01). This range is empirically defined for mammalian gene counts [7].
  • For each ( \alpha ) value, fit a LOESS regression model predicting ( \lambda ) from ( f ).
  • Calculate the Bayesian Information Criterion (BIC) for each fitted model: ( \text{BIC} = c \times \ln\left(\frac{\text{RSS}}{c}\right) + k \times \ln(c) ) where RSS is the residual sum of squares, and ( k ) is the model degrees of freedom [7].
  • Select the ( \alpha ) value that yields the lowest BIC, thus avoiding overfitting by automatically determining the optimal bandwidth [7].

• Two-Step Robust Regression:

  • Using the optimal ( \alpha ), perform an initial LOESS fit.
  • Apply Tukey's biweight robust method to identify and assign zero weight to outlier genes that may disproportionately influence the model [7].
  • Perform a final LOESS regression using the updated weights.

• Gene Selection:

  • Select the final set of HVGs as those genes whose actual average expression level ( \lambdaj ) is significantly higher than the value predicted by the final LOESS model for their positive ratio ( fj ) [7].

Protocol 2: RMT-Guided Sparse PCA for Debiased Dimensionality Reduction

This protocol uses Random Matrix Theory (RMT) to guide Sparse PCA, denoising the principal components and mitigating bias from high-dimensional noise [9].

I. Research Reagent Solutions

• Computing Environment: Python with scientific computing libraries (NumPy, SciPy).
• Input Data: A pre-processed and possibly filtered scRNA-seq matrix.
• Key Algorithms: Sinkhorn-Knopp algorithm for biwhitening, any standard sparse PCA algorithm (e.g., from scikit-learn).

II. Step-by-Step Workflow

• Data Biwhitening:
  • Goal: Simultaneously stabilize variance across genes and cells to create a matrix where the theoretical spectral distribution is known.
  • Estimate two diagonal matrices ( C ) (for cells) and ( D ) (for genes) with positive entries such that the cell-wise and gene-wise variances of the transformed matrix ( Z = C X D ) are approximately 1 [9].
  • This step helps resolve the obstacle of estimating the underlying covariance matrices needed for RMT.

• Identify the Outlier Eigenspace:

  • Compute the sample covariance matrix ( S ) of the biwhitened data.
  • Use RMT to separate the noise spectrum from the signal (outlier) eigenvalues. The RMT mapping identifies outlier eigenvalues ( \lambda ) that lie outside the support of the theoretical noise distribution ( \rho_S ) [9].

• Guide Sparse PCA with RMT:

  • The core of the mitigation strategy is to use the RMT results to automatically select the sparsity parameter in a sparse PCA algorithm.
  • RMT provides a prediction for the angle between the true signal eigenvector and the outlier eigenspace [9].
  • The sparsity level is chosen so that the inferred sparse principal components form a subspace that matches this RMT-predicted angle as closely as possible. This renders the sparse PCA process nearly parameter-free and robust [9].

• Downstream Analysis:

  • Use the resulting sparse principal components for clustering (e.g., cell type classification) and visualization. The method has been shown to systematically improve the reconstruction of the principal subspace [9].

The Scientist's Toolkit

This section details essential computational reagents and resources for implementing the protocols described in this document.

Table 3: Essential Research Reagents and Computational Tools

Item Name	Function / Purpose	Example Use Case / Note
LOESS Regression	Fits a smooth curve to a scatterplot using local weighted regression.	Core of the GLP algorithm for modeling the non-linear relationship between gene expression and positive ratio [7].
Bayesian Information Criterion (BIC)	A criterion for model selection that penalizes model complexity.	Used in GLP to automatically select the optimal LOESS bandwidth, preventing overfitting [7].
Random Matrix Theory (RMT)	A field of mathematics that studies the properties of large random matrices.	Used to distinguish signal (outlier) eigenvalues from noise in high-dimensional covariance matrices [9].
Biwhitening Algorithm	A preprocessing technique that normalizes the row and column variances of a data matrix.	Prepares scRNA-seq data for RMT analysis by stabilizing variances and ensuring a known noise distribution [9].
Sparse PCA Algorithms	A variant of PCA that produces principal components with many zero loadings, enhancing interpretability.	Denoises the principal components; its sparsity parameter is set automatically via RMT guidance [9].
K-fold Cross-Validation	A resampling technique used to evaluate models on limited data by splitting it into K subsets.	A general-purpose method for detecting overfitting by assessing model performance on held-out data [49] [50].

Benchmarking HVG Methods: How to Quantify Performance and Biological Fidelity

In the analysis of high-dimensional genomic data, particularly single-cell RNA sequencing (scRNA-seq), feature selection serves as a critical preprocessing step that profoundly impacts all subsequent analyses. The process of identifying highly variable genes (HVGs) before applying dimensionality reduction techniques like Principal Component Analysis (PCA) presents researchers with a fundamental choice: which methodological approach will most accurately distinguish biological signal from technical noise? This decision is frequently complicated by the proliferation of computational methods, with over 400 available approaches for scRNA-seq analysis alone [56]. Establishing gold standards through rigorous benchmarking is therefore not merely an academic exercise but an essential practice for generating biologically valid and reproducible results.

Benchmarking studies provide the scientific community with neutral, evidence-based comparisons of method performance using well-characterized reference datasets [56]. In the context of HVG selection for PCA, such evaluations are particularly crucial because the chosen features directly influence the preservation of biological variation during dimensionality reduction. As recent research confirms, "feature selection methods affect the performance of scRNA-seq data integration and querying," reinforcing that this initial analytical step has cascading effects throughout the entire research pipeline [10]. This application note examines the complementary roles of real and simulated data in establishing benchmarking gold standards, with a specific focus on protocols for HVG selection prior to PCA in transcriptomic studies.

Theoretical Foundations: Real Data vs. Simulations

Benchmarking datasets generally fall into two categories: simulated (synthetic) and real (experimental) data, each with distinct advantages and limitations that must be strategically balanced in any comprehensive evaluation.

The Role of Simulated Data

Simulated data offers the fundamental advantage of known ground truth, enabling precise quantitative assessment of method performance. By introducing controlled biological signals into a realistic noise structure, simulations allow researchers to systematically evaluate how well HVG selection methods recover known differentially expressed genes or cell-type markers. Well-designed simulations incorporate empirical data properties to ensure biological relevance, examining characteristics such as dropout profiles and dispersion-mean relationships in scRNA-seq data [56]. The recently developed "exvar" R package exemplifies how simulated data can validate analytical pipelines for genetic variation analysis [57].

A key application of simulated data lies in testing specific methodological aspects under controlled conditions. For instance, simulations can systematically vary parameters like cell population abundances, effect sizes, or technical noise levels to determine how these factors influence HVG selection performance. Furthermore, simulations are invaluable for scalability testing, allowing researchers to generate datasets of virtually any size to evaluate computational efficiency and memory requirements [56].

The Role of Real Benchmarking Data

In contrast to simulated data, real benchmarking datasets capture the full complexity of biological systems and technical artifacts inherent to experimental protocols. These datasets often serve as community gold standards, enabling direct comparison across methods and laboratories. In supervised learning contexts, manually annotated real datasets provide reference labels—such as expertly curated cell type markers—that serve as provisional ground truth for evaluation [56].

Several strategies enhance the utility of real data for benchmarking HVG selection methods. Experimental designs with spiked-in synthetic RNA molecules at known concentrations create hybrid datasets that combine biological complexity with controlled elements of ground truth [56]. Similarly, fluorescence-activated cell sorting (FACS) of known cell populations prior to scRNA-seq provides validated cellular identities for evaluating whether HVG selection preserves biologically relevant variation [56]. Community resources such as the BioTuring Single-Cell Atlas exemplify how aggregated real datasets from diverse biological contexts can serve as benchmarking references [58].

Toward an Integrated Approach

The most robust benchmarking frameworks strategically integrate both simulated and real data approaches [59]. Simulations enable controlled stress-testing of methods under specific challenges, while real data ensures that performance evaluations reflect biologically meaningful scenarios. This mixed-methods approach mirrors the Clinical Scenario Evaluation framework, which systematically assesses analytical performance across diverse conditions [59].

Recent research emphasizes that the distinction between these approaches is becoming increasingly blurred, with real data informing realistic simulation models, and simulated data helping to interpret results from complex biological datasets [59]. For HVG selection prior to PCA, this integration is particularly valuable, as it allows researchers to simultaneously evaluate technical performance metrics (using simulations) and biological relevance (using real data).

Performance Comparison of Benchmarking Approaches

Table 1: Comparative analysis of benchmarking data types for evaluating HVG selection methods

Aspect	Simulated Data	Real Benchmarking Data
Ground Truth	Known and precisely controllable	Often incomplete or inferred; may use experimental proxies
Biological Complexity	Limited to modeled parameters	Captures full, often uncharacterized, biological complexity
Technical Artifacts	Must be explicitly incorporated	Inherently present from sample prep to sequencing
Scalability Testing	Virtually unlimited	Constrained by available samples and sequencing costs
Community Standards	Less established; models vary	Well-established gold standards (e.g., reference atlases)
Primary Strength	Controlled performance quantification	Biological relevance and methodological stress-testing
Common Applications	Method development, parameter sweeps, power analyses	Final performance validation, community challenges

Experimental Protocols

Comprehensive Benchmarking Workflow for HVG Selection

This protocol outlines a structured approach for benchmarking highly variable gene selection methods, integrating both simulated and real data evaluation strategies.

Stage 1: Experimental Design and Scope Definition

• Define Benchmarking Purpose: Clearly specify whether the study is a neutral method comparison or evaluation of a new method against established alternatives [56].
• Select Methods for Inclusion: For neutral benchmarks, include all available methods meeting predefined criteria (e.g., software availability, documentation). For method development benchmarks, compare against state-of-the-art and baseline approaches [56].
• Identify Evaluation Metrics: Select metrics that comprehensively assess different performance aspects. For HVG selection prior to PCA, include:
  • Biological conservation: Adjusted Rand Index (ARI), normalized mutual information (NMI), silhouette coefficients [7] [10]
  • Batch effect removal: Batch ASW, PCR batch [10]
  • Computational efficiency: Runtime, memory usage [58]
• Establish Baseline Comparisons: Include negative controls (e.g., random gene selection) and positive controls (e.g., established HVG selection methods) to contextualize performance [10].

Stage 2: Data Preparation and Curation

• Real Data Collection: Curate diverse datasets representing various biological contexts, technologies, and experimental designs. Recent benchmarks have utilized 20+ scRNA-seq datasets from different biological contexts to ensure comprehensive evaluation [7].
• Real Data Preprocessing: Apply consistent quality control thresholds across all datasets (e.g., minimum gene detection in 3 cells) [7].
• Simulated Data Generation: Develop simulation models that capture key characteristics of real data. For scRNA-seq benchmarking, ensure simulations accurately represent:
  • Mean-variance relationships [7]
  • Dropout characteristics and sparsity patterns [7]
  • Batch effects and technical variability [10]
• Ground Truth Definition: For real data, establish proxy ground truth through experimental validation (e.g., FACS sorting) or expert curation. For simulated data, precisely document introduced signals.

Stage 3: Method Execution and Evaluation

• Standardized Implementation: Execute all methods using consistent preprocessing steps and parameter settings appropriate for each approach.
• HVG Selection: Apply each method to identify highly variable genes across all benchmark datasets. The recently developed GLP method, for example, uses optimized LOESS regression between gene average expression and positive ratio to select biologically informative genes [7].
• Downstream Analysis: Apply PCA to the selected HVGs and proceed with standard analytical workflows (clustering, visualization) [9].
• Performance Quantification: Calculate predefined metrics comparing outcomes against ground truth or reference standards.
• Statistical Analysis: Perform appropriate statistical tests to determine significant performance differences, accounting for multiple comparisons.

Protocol for Evaluating HVG Selection in PCA Applications

This specialized protocol focuses specifically on benchmarking how HVG selection impacts PCA performance in scRNA-seq analysis.

Data Preparation

• Select Benchmark Datasets: Choose 3-5 publicly available scRNA-seq datasets with established cell type annotations. Ensure dataset diversity in:
  • Number of cells (1,000-50,000)
  • Number of cell types (5-20)
  • Sequencing technology (10X Genomics, Smart-seq2, etc.)
• Generate Simulated Data: Use tools like Splatter to simulate datasets with known cell populations and differential expression patterns. Incorporate realistic technical noise and batch effects.
• Preprocessing: Apply consistent normalization (e.g., logTPM) and quality control across all datasets.

Method Application

• Apply HVG Selection Methods: Execute a representative set of HVG selection approaches:
  • Statistical methods: Seurat's vst, SCTransform [10]
  • Dropout-based methods: M3Drop, NBDrop [7]
  • Novel approaches: GLP (LOESS-based) [7]
• Select Gene Subsets: For each method, select the top 2,000 highly variable genes for downstream analysis [10].
• Perform PCA: Apply PCA to the selected gene subsets using consistent parameters across all methods.
• Downstream Analysis: Proceed with clustering and visualization using the PCA embeddings.

Performance Evaluation

• Biological Conservation: Calculate ARI and NMI between known cell type labels and clustering results [7] [10].
• Batch Integration: For datasets with multiple batches, compute batch correction metrics (e.g., Batch ASW) on the PCA embeddings [10].
• Computational Efficiency: Record runtime and memory usage for each HVG selection method.
• Stability Assessment: Evaluate method robustness through subsampling analyses or jackknife resampling.

Essential Research Reagents and Tools

Table 2: Key computational tools for benchmarking HVG selection methods

Tool/Resource	Type	Primary Function	Relevance to HVG/PCA Benchmarking
GLP [7]	R Method	HVG selection using LOESS regression	Identifies biologically informative genes by relationship between positive ratio and expression
exvar [57]	R Package	Genetic variation analysis pipeline	Provides validated workflow for genetic data analysis benchmarking
BBrowserX [58]	Platform	scRNA-seq data analysis with reference atlas	Enables comparison against curated single-cell database
Seurat/Scanpy [58]	Toolkit	Standard scRNA-seq analysis	Provides established HVG selection implementations for baseline comparison
scIB [10]	Pipeline	Benchmarking integration methods	Offers metric collection and evaluation framework
Temporal GeneTerrain [60]	Visualization	Dynamic expression visualization	Enables interpretation of temporal patterns in selected HVGs

Advanced Applications and Future Directions

Temporal and Multi-Modal Benchmarking

As single-cell technologies evolve, benchmarking frameworks must adapt to address increasingly complex experimental designs. Temporal sequencing data introduces additional dimensions for evaluating HVG selection, where methods must identify genes with dynamic expression patterns across time points rather than simply high variance [61]. Approaches like Temporal GeneTerrain enable visualization of these dynamic patterns, providing new evaluation modalities beyond static clustering metrics [60].

Similarly, the integration of multi-modal data (transcriptome, epigenome, proteome) presents novel challenges for feature selection benchmarking. In these contexts, HVG selection must not only preserve cell-type resolution but also facilitate cross-modal integration. Recent benchmarks highlight that feature selection significantly impacts integration quality, with effects propagating through subsequent PCA and clustering steps [10].

Machine Learning-Enhanced Benchmarking

Artificial intelligence is transforming benchmarking methodologies through improved simulation fidelity and evaluation rigor. AI-based bioinformatics tools now demonstrate accuracy improvements up to 30% with simultaneous reductions in processing time [62]. These advances enable more realistic data generation for benchmarking and more sophisticated performance evaluation.

The emergence of large language models for sequence analysis presents particularly promising opportunities [62]. These models could generate synthetic benchmarking datasets with more realistic biological patterns or predict optimal HVG selection strategies for specific experimental contexts. Furthermore, AI-powered platforms like Nygen Insights offer automated cell annotation with confidence scores, creating new opportunities for standardized evaluation [58].

Community Standards and Reproducibility

The establishment of community-wide benchmarking standards remains essential for methodological progress. Initiatives such as the Open Problems in Single-Cell Analysis project provide standardized evaluation frameworks and metric implementations [10]. Similarly, the development of R packages like "exvar" that containerize analysis workflows using Docker enhances reproducibility across computing environments [57].

Future benchmarking efforts should prioritize documentation of failure modes in addition to success metrics, providing method developers with clear guidance for improvement. Furthermore, as genomic data volumes grow exponentially, benchmarking must increasingly address scalability and computational efficiency alongside statistical performance [62] [58].

Establishing gold standards for benchmarking HVG selection methods requires methodological rigor and strategic integration of complementary data types. Simulated data provides controlled performance quantification under known ground truth, while real benchmarking datasets capture the biological complexity and technical challenges of experimental data. The protocols outlined in this application note provide structured frameworks for comprehensive method evaluation, with particular attention to how HVG selection impacts subsequent PCA and downstream analyses.

As single-cell technologies continue to evolve, maintaining robust benchmarking practices will be essential for ensuring that methodological advances translate to biologically meaningful discoveries. By adopting standardized evaluation frameworks and balancing simulated with real data assessments, the research community can establish evidence-based best practices that enhance reproducibility and accelerate progress in genomic science.

In the analysis of high-dimensional biological data, such as single-cell RNA sequencing (scRNA-seq), clustering serves as a fundamental unsupervised learning technique for discovering distinct cell types and states. Evaluating the quality and robustness of these clusters is paramount for drawing accurate biological conclusions. This application note details four key clustering evaluation metrics—Adjusted Rand Index (ARI), Normalized Mutual Information (NMI), Silhouette Coefficient, and Lack of a dedicated LISI section—framed within the context of a broader thesis investigating the impact of filtering highly variable genes prior to Principal Component Analysis (PCA). These metrics provide complementary views on clustering validity, from alignment with ground truth to intrinsic data structure assessment, and are essential tools for researchers, scientists, and drug development professionals.

Metric Definitions and Theoretical Foundations

Adjusted Rand Index (ARI)

The Adjusted Rand Index (ARI) is a statistical measure that quantifies the similarity between two clusterings by comparing their pair-wise agreements, while correcting for chance. Its value ranges from -1 to 1, where 1 indicates perfect agreement, 0 indicates random labelling, and -1 indicates completely independent labellings [63]. It is computed using a contingency table and is defined as: [ \text{ARI}(P^, P) = \frac{ \sum_{i,j} \binom{N_{ij}}{2} - \frac{ \sum_i \binom{N_i}{2} \sum_j \binom{N_j}{2} }{ \binom{N}{2} } }{ \frac{1}{2} \left[ \sum_i \binom{N_i}{2} + \sum_j \binom{N_j}{2} \right] - \frac{ \sum_i \binom{N_i}{2} \sum_j \binom{N_j}{2} }{ \binom{N}{2} } } ] where ( N_{ij} ) is the number of data points of the true class ( C_j^ ) in partition ( P^* ) assigned to cluster ( Ci ) in partition ( P ), ( Ni ) is the number of points in cluster ( Ci ), and ( Nj ) is the number of points in class ( C_j^* ) [63].

Normalized Mutual Information (NMI)

Normalized Mutual Information (NMI) is an information-theoretic measure that assesses the agreement between two clusterings by quantifying the amount of information shared between them. It is based on the concepts of entropy and mutual information, normalized to a range of [0, 1], where 1 signifies perfect correlation [64]. Higher values indicate greater shared information between the cluster assignments and the ground truth labels.

Silhouette Coefficient

The Silhouette Coefficient is an intrinsic metric that evaluates clustering quality without requiring ground truth labels. It measures how similar a data point is to its own cluster compared to other clusters by analyzing within-cluster and between-cluster distances [64]. The coefficient ranges from -1 to 1, where values near 1 indicate dense, well-separated clusters, values near 0 suggest overlapping clusters, and negative values signify potential misassignment [64] [65]. The score for a sample is calculated as ( s = \frac{b - a}{\max(a, b)} ), where ( a ) is the mean intra-cluster distance and ( b ) is the mean distance to the nearest cluster [64].

LISI

While the search results do not contain a specific section detailing the theoretical foundations of the Local Inverse Simpson's Index (LISI), it is recognized in the field as a metric for assessing batch integration in single-cell analysis, with higher LISI values indicating better mixing of batches. This note acknowledges its importance, and protocols for its calculation can be found in specialized single-cell omics software documentation.

Performance Comparison and Benchmarking

Benchmarking studies across diverse datasets provide critical insights into the practical performance of clustering algorithms as evaluated by these metrics.

Table 1: Comparative Algorithm Performance (ARI) on High-Dimensional Datasets [65]

Algorithm	Dataset	Raw Data	PCA	t-SNE	UMAP
K-means	MNIST	0.367 ± 0.001	0.494 ± 0.002	0.721 ± 0.015	0.758 ± 0.012
K-means	Fashion-MNIST	0.341 ± 0.002	0.411 ± 0.003	0.587 ± 0.021	0.625 ± 0.019
K-means	UCI HAR	0.399 ± 0.006	0.586 ± 0.004	0.684 ± 0.009	0.712 ± 0.008
DBSCAN	MNIST	0.245 ± N/A	N/A	N/A	N/A

Table 2: scMINER Clustering Performance Benchmarking (ARI) on Ground-Truth Datasets [66]

Methodology	Average ARI	Key Finding
scMINER (MI-based)	0.84	Highest average ARI, superior accuracy and purity
Seurat / Scanpy (Graph-based)	Lower than 0.84	Inconsistent performance, worse on smaller datasets
SC3s (Consensus k-means)	Lower than 0.84	N/A
scVI / scDeepCluster (Deep Learning)	Lower than 0.84	N/A

A comprehensive comparative analysis on synthetic circular data also found that density-based and graph-based algorithms like DBSCAN and Spectral Clustering consistently outperform traditional techniques like K-means when evaluated with ARI and NMI on non-linear data structures [67]. These findings underscore that algorithm performance is highly dependent on data geometry, a critical consideration when analyzing the effects of gene filtering.

Experimental Protocols

Protocol 1: Evaluating Clustering with Ground Truth via ARI and NMI

This protocol is used when true cluster labels are available (e.g., known cell types) to benchmark clustering results.

• Input Preparation: True cluster labels (labels_true) and algorithm-predicted cluster labels (labels_pred).
• Software Environment: Python with scikit-learn library installed.
• Code Implementation:

• Interpretation: An ARI score close to 1 indicates near-perfect agreement with the ground truth, while a score of 0 suggests a random assignment. Similarly, an NMI score of 1 indicates perfect shared information [64] [63].

Protocol 2: Intrinsic Cluster Validation with Silhouette Coefficient

This protocol assesses cluster quality without ground truth by examining cluster cohesion and separation.

• Input Preparation: The feature matrix (data) used for clustering and the resulting cluster assignments (labels_pred).
• Software Environment: Python with scikit-learn.
• Code Implementation:

• Interpretation: Scores near 1 mean clusters are dense and well-separated. Scores around 0 indicate overlapping clusters, and negative scores imply that samples might have been assigned to the wrong clusters [64] [65].

Protocol 3: Integrated Workflow for Preprocessing and Evaluation

This workflow integrates the evaluation metrics into a pipeline that tests the core thesis hypothesis regarding highly variable gene filtering.

Diagram 1: HVG Filtering Evaluation Workflow (63 characters)

The Scientist's Toolkit: Research Reagent Solutions

Table 3: Essential Computational Tools for Clustering Evaluation

Tool / Resource	Function	Application Context
scikit-learn (Python)	Provides implementations of ARI, NMI, and Silhouette Coefficient.	General-purpose clustering evaluation and algorithm development [64] [65].
scMINER	An integrative, mutual information-based framework for single-cell data clustering and network inference.	Specifically designed for scRNA-seq data; demonstrated superior ARI in benchmarking [66].
Seurat / Scanpy	Comprehensive toolkits for single-cell analysis, including graph-based clustering methods.	Standard workflows in scRNA-seq analysis; performance can be benchmarked against scMINER [66].
UMAP	Non-linear dimensionality reduction technique.	Preprocessing step to improve clustering performance on complex datasets before evaluation [65].

The selection of appropriate evaluation metrics is critical for validating clustering results, especially in a specialized research context like filtering highly variable genes before PCA. ARI and NMI offer robust, externally validated measures of cluster accuracy against a known ground truth, making them ideal for benchmarking the biological plausibility of results. In contrast, the Silhouette Coefficient provides an intrinsic assessment of cluster compactness and separation, useful when true labels are unavailable.

Evidence suggests that mutual information-based frameworks like scMINER can achieve superior performance (ARI > 0.84) in single-cell clustering, highlighting the potential of non-linear dependency measures over traditional linear methods [66]. Furthermore, preprocessing with non-linear dimensionality reduction techniques like UMAP has been shown to consistently enhance clustering quality across various algorithms, as measured by these metrics [65]. This directly informs the thesis research: the choice of HVG filter and subsequent PCA can be rigorously evaluated by observing its impact on final ARI, NMI, and Silhouette scores, guiding the development of more accurate and robust analytical pipelines for drug development and clinical research.

The analysis of single-cell RNA sequencing (scRNA-seq) data presents significant challenges due to its high-dimensional and sparse nature. A critical preprocessing step in this analysis is the selection of highly variable genes (HVGs), which aims to identify genes that exhibit significant cell-to-cell variation beyond what would be expected by technical noise. This process is fundamentally important for filtering genes before principal component analysis (PCA), as it enhances the biological signal captured during dimensionality reduction, improves computational efficiency, and facilitates downstream analyses such as clustering and cell type identification [7] [68]. The core premise of HVG selection rests on the assumption that genes with high biological variability are more likely to be informative for distinguishing cell types and states.

This protocol provides a comprehensive comparative analysis of four prominent HVG selection methods: Seurat and Scanpy as established frameworks, and GLP (Genes identified through LOESS with Positive Ratio) and scGBM (single-cell Generalized Binary Model) as newer approaches. While Seurat and Scanpy represent widely adopted tools in the scRNA-seq analysis ecosystem, GLP and scGBM incorporate novel statistical frameworks designed to address specific limitations of existing methods. The performance of these methods has substantial implications for the entire analytical pipeline, as the choice of HVGs directly influences the quality of dimensionality reduction through PCA, which in turn affects clustering accuracy and biological interpretation [69] [7]. A recent comprehensive benchmark evaluating 59 computational methods for selecting marker genes in scRNA-seq data highlighted the critical importance of method selection, though interestingly found that simpler methods often demonstrate competitive performance [69].

Theoretical Foundations of HVG Methods

Statistical Principles and Algorithmic Approaches

The four methods compared in this analysis employ distinct statistical frameworks and algorithmic strategies for identifying highly variable genes, each with unique theoretical underpinnings and assumptions about the data structure.

Seurat implements multiple flavors of HVG selection, with its dispersion-based approach (flavor='seurat') being one of the most widely used. This method operates by binning genes based on their average expression and then calculating a normalized dispersion (variance-to-mean ratio) for each gene within its expression bin. Genes are ranked by their normalized dispersion, and those exceeding threshold values are selected as HVGs [70]. Seurat's more recent VST (Variance Stabilizing Transformation) flavor (flavor='seurat_v3') employs a robust variance stabilization technique based on a polynomial regression model, which regresses the relationship between log(variance) and log(mean) and identifies genes that significantly deviate from this trend [70] [71]. This approach is particularly effective for handling the mean-variance relationship characteristic of count-based scRNA-seq data.

Scanpy, as a Python-based framework, largely mirrors the statistical approaches implemented in Seurat, offering compatibility with multiple HVG selection flavors including 'seurat', 'cellranger', and 'seuratv3' [70] [16]. This intentional design facilitates comparative analyses and method benchmarking across programming environments. A key distinction lies in Scanpy's implementation within the Python ecosystem, which may offer computational advantages for certain analytical workflows and integration with other Python-based machine learning libraries.

GLP (Genes identified through LOESS with Positive Ratio) introduces a novel approach that fundamentally rethinks the relationship between gene expression characteristics and biological informativeness. Instead of relying on variance or dispersion metrics, GLP leverages the positive ratio (the proportion of cells where a gene is detected) in relation to its average expression level [7]. The method applies an optimized LOESS (Locally Estimated Scatterplot Smoothing) regression to model the non-linear relationship between these two variables. Genes with expression levels significantly higher than predicted by the LOESS curve given their positive ratio are selected as highly variable. GLP incorporates an adaptive bandwidth selection using the Bayesian Information Criterion (BIC) to prevent overfitting and employs a two-step regression procedure with Tukey's biweight robust statistical method to minimize the influence of outlier genes [7].

scGBM (single-cell Generalized Binary Model), while not detailed in the available search results, theoretically represents a class of methods that utilize generalized binary models to account for the excess zeros (dropouts) characteristic of scRNA-seq data. Such approaches typically model gene expression as a mixture of a binary component (representing detection vs. non-detection) and a conditional continuous component (representing expression level when detected).

Table 1: Theoretical Foundations of HVG Selection Methods

Method	Primary Statistical Approach	Key Assumptions	Distinguishing Features
Seurat	Normalized dispersion/VST	Mean-variance relationship follows definable pattern	Multiple algorithmic flavors; Extensive benchmarking
Scanpy	Compatible with Seurat approaches	Compatible with Seurat assumptions	Python implementation; Seamless integration with scverse ecosystem
GLP	LOESS regression on positive ratio	Positive ratio predicts expected expression for non-informative genes	Specifically addresses dropout noise; Adaptive bandwidth selection
scGBM	Generalized binary model	Expression can be modeled as mixture of binary and continuous components	Explicitly accounts for dropout events; Not benchmarked in available literature

Method-Specific Parameter Considerations

Each HVG selection method requires careful parameterization to optimize performance for specific dataset characteristics:

For Seurat's dispersion-based approach, critical parameters include n_bins (number of expression bins), mean.cutoff (threshold for mean expression), and dispersion.cutoff (threshold for normalized dispersion) [70]. The VST method requires specification of n_top_genes rather than explicit cutoffs, and uses span to control the degree of smoothing in the loess curve fit [70].

Scanpy maintains parameter compatibility with Seurat, with flavor='seurat' utilizing min_mean, max_mean, min_disp, and max_disp thresholds, while flavor='seurat_v3' uses n_top_genes and span parameters analogous to Seurat's implementation [70] [16].

GLP introduces specialized parameters including the alpha range for BIC-optimized bandwidth selection (typically 0.01 to 0.1) and robust regression parameters for outlier handling [7]. The method automatically filters genes detected in fewer than 3 cells to maintain reliability.

Performance Benchmarking and Comparative Analysis

Evaluation Metrics and Experimental Design

Rigorous benchmarking of HVG selection methods requires comprehensive evaluation using multiple complementary metrics that capture different aspects of performance. The most relevant metrics include:

• Adjusted Rand Index (ARI): Measures the similarity between clustering results and ground truth cell type labels, with higher values indicating better preservation of biological population structure [7].
• Normalized Mutual Information (NMI): Quantifies the mutual information between clustering results and reference labels, normalized by entropy [7].
• Silhouette Coefficient: Evaluates clustering quality based on the compactness of clusters and separation between clusters [7].
• Computational Efficiency: Assesses memory usage and computation time, particularly important for large-scale datasets.
• Downstream Analysis Performance: Evaluates how HVG selection impacts differential expression testing, trajectory inference, and other analytical tasks.

A recent comprehensive benchmark study evaluated 59 marker gene selection methods using 14 real scRNA-seq datasets and over 170 simulated datasets, providing valuable insights into the relative performance of different approaches [69]. While this study focused on marker gene selection rather than HVG selection specifically, the findings remain relevant as both tasks aim to identify biologically informative genes.

Comparative Performance Results

Evaluation studies have demonstrated distinct performance patterns across HVG selection methods:

GLP has shown particularly strong performance in comprehensive benchmarks, consistently outperforming eight state-of-the-art feature selection methods across three primary metrics (ARI, NMI, and silhouette coefficient) on 20 scRNA-seq datasets from diverse biological contexts [7]. The method's strength appears to derive from its robust handling of technical noise and dropout events, which are pervasive challenges in scRNA-seq data. Additionally, GLP has demonstrated capability to enhance downstream analysis performance, including differential gene expression analysis and trajectory inference [7].

The benchmark of 59 methods revealed that simpler statistical approaches, including the Wilcoxon rank-sum test and Student's t-test, often demonstrate competitive performance compared to more complex machine learning-based methods [69]. This finding suggests that the most computationally intensive approaches do not necessarily provide superior performance for HVG selection.

For Seurat and Scanpy, performance varies depending on the specific flavor employed. The seurat_v3/VST approach generally shows improved performance over earlier dispersion-based methods, particularly for datasets with strong mean-variance relationships [70] [71]. Both frameworks effectively handle batch effects when the batch_key parameter is specified, implementing a strategy that selects HVGs within each batch separately before merging results [70].

Table 2: Performance Comparison of HVG Selection Methods

Method	Clustering Accuracy (ARI)	Robustness to Dropouts	Computational Efficiency	Batch Effect Handling
Seurat	High with VST flavor	Moderate	High	Effective with batch_key parameter
Scanpy	High with seurat_v3 flavor	Moderate	High	Effective with batch_key parameter
GLP	Highest in benchmark studies	Excellent	Moderate	Not specifically evaluated
scGBM	Not available	Theoretical strength	Not available	Not available

Impact on Downstream PCA and Analysis

The selection of HVGs fundamentally shapes subsequent dimensionality reduction through PCA, as PCA operates exclusively on the gene subset identified as highly variable. The quality of HVG selection directly influences:

• Variance Explained: Informative HVG selection increases the biological signal captured in the first principal components.
• Cluster Separation: Appropriate HVGs enhance the separation of distinct cell types in PCA space.
• Batch Effect Mitigation: When batch-aware HVG selection is employed, technical artifacts are reduced in the PCA embedding.
• Trajectory Inference: For developmental datasets, biologically relevant HVGs preserve continuous processes in the reduced dimension space.

Evidence suggests that method choice significantly impacts these downstream results. For instance, the GLP method demonstrated improved preservation of biological information in downstream analyses compared to other feature selection methods [7]. Similarly, proper HVG selection in Seurat and Scanpy workflows has been shown to produce embeddings that effectively separate cell types and facilitate accurate clustering [16] [72].

Experimental Protocols and Implementation

Standardized Workflow for HVG Comparison

To conduct a rigorous comparative analysis of HVG methods, researchers should implement the following standardized workflow:

• Data Preprocessing:

  • Begin with a quality-controlled count matrix filtered for low-quality cells and genes.
  • Perform basic normalization (e.g., median normalization or CP10k) and log-transformation for methods requiring normalized input.
  • For Seurat's SCTransform or similar approaches, utilize raw counts as input.

• HVG Selection with Multiple Methods:

  • Apply each HVG method to the same preprocessed dataset using standardized parameters.
  • Select a consistent number of top genes (e.g., 2000-3000 HVGs) for fair comparison.
  • For batch-containing datasets, utilize batch-aware implementations where available.

• Downstream Processing:

  • Perform PCA on each HVG set using identical parameters (number of components, scaling options).
  • Construct neighbor graphs and generate UMAP/t-SNE embeddings for visualization.
  • Apply clustering algorithms (e.g., Leiden, Louvain) with fixed resolution parameters.

• Evaluation:

  • Calculate performance metrics (ARI, NMI, silhouette) using known cell type labels.
  • Assess computational requirements (memory, time) for each method.
  • Evaluate biological plausibility through marker gene expression patterns.

Method-Specific Implementation Protocols

Seurat Implementation:

Scanpy Implementation:

GLP Implementation:

The GLP method requires specialized implementation as described in the original publication [7]. The core algorithm involves:

• Calculating positive ratio and mean expression for each gene.
• Applying BIC to determine optimal LOESS bandwidth parameter.
• Performing two-step LOESS regression with outlier detection.
• Selecting genes with positive residuals exceeding a defined threshold.

The Scientist's Toolkit: Essential Research Reagents

Table 3: Essential Computational Tools for HVG Analysis

Tool/Resource	Function	Application Context
Seurat (R)	Comprehensive scRNA-seq analysis	End-to-end workflow with specialized HVG methods
Scanpy (Python)	Scalable scRNA-seq analysis	Python-based workflows with Seurat-compatible functions
scikit-misc (Python)	Additional statistical functions	Required for Scanpy's seurat_v3 flavor
DoubletDetection	Doublet identification	Quality control before HVG selection
SoupX	Ambient RNA correction	Data cleaning to improve HVG selection
scIB	Integration benchmarking	Evaluating HVG method performance

Visualization and Decision Support

HVG Selection Workflow Diagram

HVG Selection and Evaluation Workflow: This diagram illustrates the comprehensive process for comparing highly variable gene selection methods, from raw data processing through to performance evaluation.

Method Selection Decision Framework

HVG Method Selection Framework: This decision diagram provides guidance for selecting the most appropriate highly variable gene selection method based on dataset characteristics and research constraints.

Based on comprehensive benchmarking studies and methodological considerations, we provide the following recommendations for selecting HVG methods in different research contexts:

For most standard applications with datasets of moderate size (≤50,000 cells), GLP represents an excellent choice due to its superior performance in benchmark studies and robust handling of technical noise [7]. When GLP implementation is not feasible, Seurat's VST method (selection.method = "vst") or Scanpy's seurat_v3 flavor provide strong alternatives with extensive community validation [70] [16] [71].

In scenarios with pronounced batch effects, both Seurat and Scanpy offer effective batch-aware HVG selection through their batch_key parameters, which identifies genes that are variable across batches, thereby reducing technical artifacts in downstream analyses [70].

For large-scale datasets requiring computational efficiency, the standard Seurat and Scanpy implementations provide optimized performance while maintaining good biological fidelity. The Wilcoxon rank-sum test, identified as a top performer in the comprehensive marker gene benchmark, may also serve as a computationally efficient and effective alternative [69].

In methodological studies where benchmarking and comparison are primary goals, employing multiple HVG selection approaches provides valuable validation of result robustness. The complementary strengths of different statistical approaches can reveal consistent biological signals while highlighting method-specific biases.

Future methodological development should focus on integrating the strengths of these approaches—particularly combining GLP's innovative use of positive ratio relationships with the batch-correction capabilities of established frameworks. Additionally, continued benchmarking on diverse biological systems and experimental conditions will further refine our understanding of optimal HVG selection strategies for specific research contexts.

The critical importance of appropriate HVG selection for subsequent PCA and dimensionality reduction cannot be overstated, as this foundational step fundamentally shapes all downstream analyses and biological interpretations in single-cell RNA sequencing studies.

In single-cell RNA sequencing (scRNA-seq) analysis, feature selection through the identification of highly variable genes (HVGs) is a critical preprocessing step that profoundly influences all subsequent downstream results. This procedure filters the vast number of measured genes (often >20,000) down to a subset (typically 1,000-5,000) that exhibit high cell-to-cell variation, under the assumption that these genes most likely reflect biologically interesting heterogeneity rather than technical noise [13] [15]. The choice of HVG selection method is not merely a technical detail; it directly determines the biological signals that will be accessible for further investigation. This case study examines how different HVG selection strategies impact clustering outcomes and differential gene expression (DGE) analysis, providing evidence-based protocols for researchers to optimize their analytical workflows within the broader context of feature selection prior to dimensionality reduction techniques like Principal Component Analysis (PCA).

Impact of HVG Selection on Downstream Analyses

Effects on Clustering and Cell Type Identification

The specific set of HVGs selected can significantly alter the topology of the cell-cell similarity graph, thereby changing the resulting clusters. Different HVG selection methods exhibit distinct biases in the expression levels of genes they prioritize, which subsequently affects the resolution of cell types and states.

Table 1: Characteristics of HVG Selection Methods

Method	Underlying Principle	Expression Bias	Reproducibility	Clustering Performance
Seurat (VST) [13] [41]	Fits a loess curve to the log(variance) and log(mean) relationship	Prefers a mix of highly and lowly expressed genes	Moderate (50-70%)	Good separation for distinct cell types
SCTransform [41]	Regularized negative binomial regression on Pearson residuals	Prefers a mix of highly and lowly expressed genes	Moderate (50-70%)	Good separation for distinct cell types
SCHS [41]	Detects genes with non-random spatial distribution in expression space	Strong preference for highly expressed genes (>85%)	High (>90%)	Good for mature cell types, may miss fine structure in progenitors
M3Drop [41]	Models relationship between mean expression and dropout rates	Prefers a mix of highly and lowly expressed genes	Low (50-70%)	Lower distinguishing capability for similar cell types
Scran [41]	Fits a trend to the variance of log-expression against the mean	Almost no lowly expressed genes selected	Lower (80-90%)	Moderate performance
GLP [7]	Optimized LOESS regression between positive ratio and mean expression	Wider dynamic range of expression	High	Improves downstream clustering accuracy

A comprehensive evaluation of nine HVG methods using hematopoietic stem/progenitor cells (HSPCs) and mature blood cells revealed that methods selecting predominantly highly expressed genes (like SCHS and Scmap) can identify broad cell classes effectively but may miss subtler developmental transitions [41]. In contrast, methods that include lowly expressed genes with high biological variability (like M3Drop) sometimes fail to separate transcriptionally similar progenitor populations. The purity of cell clusters—defined as clusters containing cells of identical function and state—varied across methods, though most maintained above 90% purity in both HSPCs and mature blood cells [41].

Effects on Differential Gene Expression Analysis

HVG selection creates a fundamental constraint on downstream DGE analysis, as only selected genes are typically considered in subsequent steps. This can inadvertently lead to missing biologically important genes that don't pass variability thresholds.

Housekeeping genes with stable expression across cell types are often filtered out during HVG selection, which can negatively impact DGE analysis if these genes are of biological interest [73]. As one expert notes: "HVGs are usually used for clustering and visualization. For everything else use all genes, or those filtered for some minimal expression in x% of all cells or clusters" [73]. This highlights the importance of aligning feature selection strategy with analytical goals.

The reproducibility of differentially expressed genes (DEGs) identified in single-cell studies of neurodegenerative diseases is notably poor, particularly for complex disorders like Alzheimer's disease (AD) and schizophrenia (SCZ), where over 85% of DEGs detected in one dataset fail to reproduce in others [74]. This reproducibility challenge is partly attributable to the initial HVG selection, which determines the candidate gene pool for differential testing.

Table 2: Reproducibility of DEGs Across Neurodegenerative Disease Studies

Disease	Number of Studies	Reproducibility of DEGs	Predictive Power (AUC)
Alzheimer's (AD)	17	<0.1% of DEGs reproduced in >3 studies	0.68
Parkinson's (PD)	6	Moderate	0.77
Huntington's (HD)	4	Moderate	0.85
Schizophrenia (SCZ)	3	Poor	0.55
COVID-19	16	Moderate	0.75

A meta-analysis of 17 AD snRNA-seq studies found that genes initially identified as DEGs frequently showed inconsistent directionality across studies. For example, LINGO1—previously highlighted as a key oligodendrocyte DEG in AD—was not consistently upregulated across multiple datasets and was even downregulated in several studies [74]. This underscores how HVG selection and subsequent DEG identification can be sensitive to dataset-specific technical artifacts and biological variability.

Experimental Protocols for HVG Evaluation

Protocol 1: Comparing HVG Methods Using the SIEVE Framework

The SIEVE (SIngle-cEll Variable gEnes) framework provides a robust strategy for identifying reproducible HVGs through multiple rounds of random sampling, minimizing stochastic noise [41].

Workflow:

• Random Sampling: Repeatedly randomly sample 70% of cells from the complete dataset (minimum 50 iterations)
• HVG Selection: Apply HVG selection methods to each subsample
• Reproducibility Assessment: Calculate the proportion of overlapping genes identified across iterations
• Robust HVG Set: Select genes with high recurrence across sampling iterations for downstream analysis

Performance Evaluation:

• Calculate classification accuracy using random forest classifiers trained on reference sets and applied to query sets
• Assess cluster purity using Calinski-Harabasz index (higher values indicate better separation) and Davies-Bouldin index (lower values indicate better separation)
• Evaluate enrichment for known cell-type marker genes

Implementation of SIEVE significantly improves classification accuracy across multiple HVG selection methods and enhances the biological relevance of selected genes, as measured by enrichment for established cell cluster markers [41].

Protocol 2: Evaluating Impact on Differential Expression Analysis

To systematically assess how HVG selection affects differential expression results, we recommend the following protocol:

• Apply Multiple HVG Methods: Process the same dataset using at least 3-4 different HVG selection methods (e.g., Seurat VST, SCTransform, SCHS)
• Perform Clustering: For each HVG set, conduct standard dimensionality reduction (PCA) and clustering (Leiden algorithm) [16]
• Differential Expression Testing: For each clustering result, perform differential expression analysis between conditions within matched cell types
• Reproducibility Assessment: Compare the resulting DEG lists across HVG methods using:
  • Overlap coefficients (e.g., Jaccard index)
  • Correlation of effect sizes for common genes
  • Consistency in the top-ranking biological pathways

This approach helps identify robust biological findings that persist across HVG selection strategies versus method-specific artifacts.

Protocol 3: Direct Marker Selection with Festem

For studies specifically focused on cell type identification, Festem provides an alternative approach that directly selects cluster-informative marker genes rather than relying on surrogate criteria like high variability [75].

Procedure:

• Initial Clustering: Perform preliminary clustering with standard HVG selection
• Marker Screening: Apply Festem to distinguish genes with heterogeneous expression distributions across cells
• Iterative Refinement: Re-cluster using Festem-selected markers
• Validation: Compare resulting clusters with known cell type markers and assess false discovery rates

Festem demonstrates improved sensitivity for identifying rare cell populations that might be overlooked by conventional HVG-based approaches [75].

Workflow Visualization

Figure 1: HVG Method Comparison Workflow. This workflow illustrates the parallel application of multiple HVG selection methods to evaluate their impact on downstream clustering and differential expression analysis.

Figure 2: SIEVE Resampling Framework. The SIEVE method identifies robust HVGs through multiple rounds of random sampling, enhancing reproducibility and biological relevance of selected genes.

Research Reagent Solutions

Table 3: Essential Tools for HVG Analysis

Tool/Software	Function	Application Context
Seurat [13]	Comprehensive scRNA-seq analysis platform	HVG selection using VST, dispersion, or SCTransform methods; R environment
Scanpy [15] [16]	Scalable Python-based single-cell analysis	HVG selection compatible with Seurat methods; Python environment
SIEVE [41]	Resampling framework for robust HVG identification	Improving reproducibility of any base HVG method; available on GitHub
Festem [75]	Direct marker gene selection for clustering	Identifying cluster-informative genes beyond high variability
GLP [7]	LOESS-based feature selection using positive ratio	Alternative approach focusing on expression positivity patterns
SCHS [41]	Gene selection based on spatial expression distribution	Identifying genes with non-random expression patterns

Discussion and Recommendations

The choice of HVG selection method should be guided by the specific biological question and experimental context. For researchers focusing on well-separated cell types, most HVG methods perform adequately, though methods with higher reproducibility (like SCHS) may provide more consistent results [41]. For investigating subtle transitions along developmental trajectories, methods that include moderately variable genes (like GLP) may be preferable [7].

When designing studies where differential expression analysis is the primary goal, we recommend:

• Using full gene sets (with minimal expression filters) for differential expression testing rather than restricting to HVGs [73]
• Applying multiple HVG methods and comparing consistency of key findings
• Implementing resampling approaches like SIEVE to identify robust gene signatures
• Considering direct marker selection methods like Festem when cell type identification is the primary objective [75]

For meta-analyses across multiple datasets, methods that prioritize reproducibility (like SumRank) can substantially improve the identification of robust differentially expressed genes [74].

The integration of HVG selection with downstream PCA is not merely a technical preprocessing step but a fundamental analytical decision that shapes biological interpretation. As single-cell technologies advance to profile increasingly complex biological systems and clinical samples, the development and rigorous evaluation of feature selection methods remains crucial for extracting meaningful biological insights from high-dimensional transcriptomic data.

Conclusion

The strategic selection of highly variable genes is not merely a preliminary step but a foundational decision that profoundly influences all subsequent single-cell RNA-seq analyses. By understanding the core principles, adeptly applying and troubleshooting modern methods, and rigorously validating outcomes, researchers can ensure their PCA and downstream results are driven by robust biological signal rather than technical noise. Future directions point towards more adaptive, model-based algorithms that automatically quantify uncertainty and better preserve rare cell populations. Mastering HVG selection is therefore paramount for unlocking the full potential of scRNA-seq in discovering novel cell types, understanding disease mechanisms, and advancing personalized therapeutics.

References

• 1. Resolution of the curse of dimensionality in single-cell RNA ... [https://pmc.ncbi.nlm.nih.gov/articles/PMC9363502/]
• 2. 9. Dimensionality Reduction [https://www.sc-best-practices.org/preprocessing_visualization/dimensionality_reduction.html]
• 3. Normalization and variance stabilization of single-cell RNA ... [https://genomebiology.biomedcentral.com/articles/10.1186/s13059-019-1874-1]
• 4. Top 10 Bioinformatics Tools for scRNA-seq in 2025 [https://neovarsity.org/blogs/top-bioinformatics-tools-for-scrna-seq]
• 5. BASiCS workflow: a step-by-step analysis of expression ... [https://pmc.ncbi.nlm.nih.gov/articles/PMC11109695/]
• 6. Redefining the high variable genes by optimized LOESS ... [https://pmc.ncbi.nlm.nih.gov/articles/PMC12001687/]
• 7. Redefining the high variable genes by optimized LOESS ... [https://bmcbioinformatics.biomedcentral.com/articles/10.1186/s12859-025-06112-5]
• 8. Exploring cell-to-cell variability and functional insights ... [https://www.nature.com/articles/s41540-025-00507-z]
• 9. Random Matrix Theory-guided sparse PCA for single-cell ... [https://arxiv.org/html/2509.15429v1]
• 10. Feature selection methods affect the performance of ... [https://www.nature.com/articles/s41592-025-02624-3]
• 11. Gene function revealed at the moment of stochastic gene ... [https://pmc.ncbi.nlm.nih.gov/articles/PMC11743767/]
• 12. High-Fidelity Transcriptome Reconstruction of Degraded ... [https://www.mdpi.com/2079-7737/14/12/1652]
• 13. Seurat - Guided Clustering Tutorial [https://satijalab.org/seurat/articles/pbmc3k_tutorial.html]
• 14. Find variable features — FindVariableFeatures • Seurat [https://satijalab.org/seurat/reference/findvariablefeatures]
• 15. scanpy.pp.highly_variable_genes - Read the Docs [https://scanpy.readthedocs.io/en/stable/api/scanpy.pp.highly_variable_genes.html]
• 16. Preprocessing and clustering - Scanpy - Read the Docs [https://scanpy.readthedocs.io/en/stable/tutorials/basics/clustering.html]
• 17. Scanpy HVG batch_key [https://discourse.scverse.org/t/scanpy-hvg-batch-key/1158]
• 18. Benchmarking large-scale single-cell RNA-seq analysis [https://pubmed.ncbi.nlm.nih.gov/41279840/]
• 19. Single-cell analyses reveal increased gene expression ... [https://pmc.ncbi.nlm.nih.gov/articles/PMC12081229/]
• 20. graph-based gene selection for single-cell RNA-seq data [https://www.sciencedirect.com/science/article/abs/pii/S1673852722000339]
• 21. Statistically principled feature selection for single cell ... [https://bmcbioinformatics.biomedcentral.com/articles/10.1186/s12859-025-06240-y]
• 22. Probe set selection for targeted spatial transcriptomics [https://www.nature.com/articles/s41592-024-02496-z]
• 23. On feature selection to disentangle cell type and state ... [https://pmc.ncbi.nlm.nih.gov/articles/PMC12593937/]
• 24. Using sctransform in Seurat [https://satijalab.org/seurat/articles/sctransform_vignette.html]
• 25. SCTransform – simple and intuitive explanation [https://biostatsquid.com/sctransform-simple-explanation/]
• 26. Issue #5832 · satijalab/seurat - FindVariableFeatures VST [https://github.com/satijalab/seurat/issues/5832]
• 27. Introduction to SCTransform, v2 regularization • Seurat [https://satijalab.org/seurat/archive/v4.3/sctransform_v2_vignette]
• 28. Introduction to scRNA-seq integration • Seurat [https://satijalab.org/seurat/archive/v4.3/integration_introduction]
• 29. Data Integration - GitHub Pages [https://nbisweden.github.io/workshop-scRNAseq/labs/seurat/seurat_03_integration.html]
• 30. Fast integration using reciprocal PCA (RPCA) [https://satijalab.org/seurat/articles/integration_rpca.html]
• 31. Scanpy – Single-Cell Analysis in Python — scanpy [https://scanpy.readthedocs.io/]
• 32. scanpy.pp.log1p - Read the Docs [https://scanpy.readthedocs.io/en/1.11.x/generated/scanpy.pp.log1p.html]
• 33. Understanding normalization and log transformation #1364 [https://github.com/scverse/scanpy/issues/1364]
• 34. How to preprocess UMI count data with analytic Pearson ... [https://scanpy-tutorials.readthedocs.io/en/latest/tutorial_pearson_residuals.html]
• 35. How to analyze a 1 million cell data set using Scanpy and ... [https://support.parsebiosciences.com/hc/en-us/articles/7704577188500-How-to-analyze-a-1-million-cell-data-set-using-Scanpy-and-Harmony]
• 36. Model-based dimensionality reduction for single-cell RNA- ... [https://pmc.ncbi.nlm.nih.gov/articles/PMC10168202/]
• 37. Principal Component Analysis (PCA): Explained Step-by- ... [https://builtin.com/data-science/step-step-explanation-principal-component-analysis]
• 38. Principal Component Analysis: A Method for Determining ... [https://pmc.ncbi.nlm.nih.gov/articles/PMC4676806/]
• 39. Systematic benchmarking of computational methods to ... [https://pubmed.ncbi.nlm.nih.gov/40968359/]
• 40. Chapter 3 Feature selection | Basics of Single-Cell ... [https://bioconductor.org/books/3.14/OSCA.basic/feature-selection.html]
• 41. SIEVE: identifying robust single cell variable genes for single ... [https://pmc.ncbi.nlm.nih.gov/articles/PMC8974938/]
• 42. Statistically principled feature selection for single cell ... [https://pmc.ncbi.nlm.nih.gov/articles/PMC12490061/]
• 43. High performance data integration for large-scale analyses ... [https://www.nature.com/articles/s41467-025-62237-4]
• 44. Integrating single-cell RNA-seq datasets with substantial ... [https://pmc.ncbi.nlm.nih.gov/articles/PMC12577435/]
• 45. A technical review of multi-omics data integration methods [https://pmc.ncbi.nlm.nih.gov/articles/PMC12315550/]
• 46. Navigating single-cell RNA-sequencing: protocols, tools ... [https://genomicsinform.biomedcentral.com/articles/10.1186/s44342-025-00044-5]
• 47. Rare Cell Isolation: Methods for Rare Immune Cells [https://www.akadeum.com/types-of-cell-separation/rare-cell-isolation-methods/?srsltid=AfmBOopTCz7e1OatS1VCtpJu1iUPUgMeWT0DUkwE0JlWoZCz2c5HF4Bx]
• 48. Rare cell classification using label-free imaging flow cytometry ... [https://pubs.rsc.org/en/content/articlehtml/2025/lc/d5lc00634a]
• 49. How to Avoid Overfitting in Machine Learning? [https://www.geeksforgeeks.org/machine-learning/how-to-avoid-overfitting-in-machine-learning/]
• 50. Overfitting in Machine Learning Explained [https://aws.amazon.com/what-is/overfitting/]
• 51. Understanding Overfitting: Strategies and Solutions [https://www.lyzr.ai/glossaries/overfitting/]
• 52. Underfitting vs. Overfitting in Machine Learning [https://www.ayadata.ai/a-guide-to-overfitting-and-underfitting-in-machine-learning/]
• 53. Biases in Machine Learning Models: Understanding and ... [https://www.aryaxai.com/article/biases-in-machine-learning-models-understanding-and-overcoming-them]
• 54. Unpacking the bias of large language models [https://news.mit.edu/2025/unpacking-large-language-model-bias-0617]
• 55. Bias recognition and mitigation strategies in artificial ... [https://www.nature.com/articles/s41746-025-01503-7]
• 56. Essential guidelines for computational method benchmarking [https://genomebiology.biomedcentral.com/articles/10.1186/s13059-019-1738-8]
• 57. Efficient and easy gene expression and genetic variation ... [https://www.nature.com/articles/s41598-025-93067-5]
• 58. Powerful scRNA-seq Data Analysis Tools in 2025 [https://www.nygen.io/resources/blog/single-cell-scrna-seq-multi-omics-data-analysis-tools]
• 59. On the role of benchmarking data sets and simulations in ... [https://pubmed.ncbi.nlm.nih.gov/36810737/]
• 60. Temporal GeneTerrain: advancing precision medicine ... [https://www.frontiersin.org/journals/bioinformatics/articles/10.3389/fbinf.2025.1602850/full]
• 61. Computational workflow for investigating highly variable ... [https://www.sciencedirect.com/science/article/pii/S2666166723003544]
• 62. The 2025 Bioinformatics Toolkit: From Sample Prep to AI- ... [https://blog.omni-inc.com/blog/the-2025-bioinformatics-toolkit-from-sample-prep-to-ai-driven-analysis]
• 63. Adjusted Rand Index - an overview | ScienceDirect Topics [https://www.sciencedirect.com/topics/computer-science/adjusted-rand-index]
• 64. 7 Evaluation Metrics for Clustering Algorithms [https://towardsdatascience.com/7-evaluation-metrics-for-clustering-algorithms-bdc537ff54d2/]
• 65. Comparative Analysis of Clustering Techniques on High ... [https://arxiv.org/html/2503.23215v1]
• 66. scMINER: a mutual information-based framework for ... [https://pmc.ncbi.nlm.nih.gov/articles/PMC12062461/]
• 67. Comparative Analysis Of Clustering Algorithms On ... [https://papers.ssrn.com/sol3/papers.cfm?abstract_id=5086452]
• 68. Computational workflow for investigating highly variable ... [https://pmc.ncbi.nlm.nih.gov/articles/PMC10331590/]
• 69. A comparison of marker gene selection methods for single-cell ... [https://genomebiology.biomedcentral.com/articles/10.1186/s13059-024-03183-0]
• 70. scanpy.pp.highly_variable_genes — scanpy [https://scanpy.readthedocs.io/en/stable/generated/scanpy.pp.highly_variable_genes.html]
• 71. Seurat - Differential expression testing [https://satijalab.org/seurat/articles/de_vignette]
• 72. Preprocessing and clustering 3k PBMCs (legacy workflow) [https://scanpy-tutorials.readthedocs.io/en/latest/pbmc3k.html]
• 73. Will filter cell by highly variable genes ( HVG ) selection marks DGE for... [https://www.biostars.org/p/9556311/]
• 74. Improving reproducibility of differentially expressed genes ... [https://www.nature.com/articles/s41467-025-62579-z]
• 75. Article Directly selecting cell-type marker genes for single ... [https://www.sciencedirect.com/science/article/pii/S2667237524001735]