Function subset_genes

src/memory/processing/hvg/mod.rs:404–466

Select highly variable genes based on dispersion scores and expression filters. Final step of HVG detection that applies thresholds to identify the most biologically relevant genes. Supports both top-N selection and threshold-based selection. ## Parameters `log_means` - Log-transformed mean expression values `dispersion_norm` - Normalized dispersion scores `n_top_genes` - Optional: select this m

(
    log_means: &[f64], // These are already log-transformed
    dispersion_norm: &[f64],
    n_top_genes: Option<usize>,
    min_mean: f64,
    max_mean: f64,
    min_dispersion: f64,
)

/// - Otherwise: Apply all threshold filters (mean and dispersion)
/// - Genes with NaN dispersions are handled appropriately for each mode
fn subset_genes(
    log_means: &[f64], // These are already log-transformed
    dispersion_norm: &[f64],
    n_top_genes: Option<usize>,
    min_mean: f64,
    max_mean: f64,
    min_dispersion: f64,
) -> anyhow::Result<Vec<bool>> {
    let mut highly_variable = vec![false; log_means.len()];

    if let Some(n_top) = n_top_genes {
        // Python's approach for n_top_genes:
        // 1. First, remove NaN values to compute threshold
        let non_nan_dispersions: Vec<f64> = dispersion_norm
            .iter()
            .filter(|&&d| !d.is_nan())
            .copied()
            .collect();

        if non_nan_dispersions.is_empty() {
            return Ok(highly_variable);
        }

        // Find the nth highest value
        let n_to_select = n_top.min(non_nan_dispersions.len());
        let mut sorted_dispersions = non_nan_dispersions.clone();
        sorted_dispersions.sort_by(|a, b| b.partial_cmp(a).unwrap_or(std::cmp::Ordering::Equal));

        let threshold = sorted_dispersions[n_to_select - 1];

        // 2. Now apply threshold to nan_to_num version (NaN → -inf)
        for i in 0..dispersion_norm.len() {
            let disp_value = if dispersion_norm[i].is_nan() {
                f64::NEG_INFINITY // Python uses -inf for NaN in final selection
            } else {
                dispersion_norm[i]
            };

            if disp_value >= threshold {
                highly_variable[i] = true;
            }
        }
    } else {
        // Original cutoff-based selection
        // Python applies nan_to_num (NaN → 0) before checking bounds
        let clean_dispersions: Vec<f64> = dispersion_norm
            .iter()
            .map(|&d| if d.is_nan() { 0.0 } else { d })
            .collect();

        // Apply mean filters
        let valid_by_mean: Vec<bool> = log_means
            .iter()
            .map(|&log_mean| log_mean > min_mean && log_mean < max_mean)
            .collect();

        for i in 0..log_means.len() {
            highly_variable[i] = valid_by_mean[i] && clean_dispersions[i] > min_dispersion;