Restore archived Python bulk SNP counting and analysis compatibility

.github/workflows/ecosystem.yml

@@ -34,7 +34,7 @@ jobs:
           sudo apt-get install -y libhts-dev libbz2-dev liblzma-dev zlib1g-dev libssl-dev pkg-config libclang-dev samtools bcftools bedtools
           uv venv .venv
           source .venv/bin/activate
-          uv pip install "maturin==1.14.1" "planemo==0.75.44"
+          uv pip install "maturin==1.14.1" "planemo==0.75.44" "galaxy-tool-util==26.0.1" "galaxy-tool-util-models==26.0.1"
           maturin develop --release
       - name: Lint Galaxy wrappers
         run: |

CHANGELOG.md

@@ -4,6 +4,20 @@ All notable changes to WASP2 will be documented in this file.
 
 ## [Unreleased]
 
+### Fixed
+- Bulk SNP counting preserves archived Python input-order priority across paired
+  mates, with read names consumed before checking aligned bases at each SNP.
+- Bulk region outputs remove one pseudocount per retained SNP and allele, restoring
+  the archived Python raw-count aggregation without changing the fitted statistics.
+- Single-cell analysis input preparation and differential comparisons now follow
+  the cells-by-variants AnnData layout throughout, including barcode groups and genotypes.
+- Differential pseudobulk sums widen to UInt64 rather than wrapping at 65,535.
+
+### Tests
+- Frozen outputs from archived Python commit `41354fc` cover bulk, single-cell,
+  and differential analysis. Linear-reference array-shape adaptations are explicit;
+  the archive's ignored phased-linear flag is reported as a reference limitation.
+
 ## [1.5.0] - 2026-07-14
 
 ### Added

environment.yml

@@ -13,7 +13,7 @@ dependencies:
   # Data processing (aligned with pyproject.toml)
   - numpy>=1.21,<2.0  # <2.0 for ABI compatibility with pysam, pyarrow
   - pandas>=1.5,<3.0
-  - polars>=0.19
+  - polars>=0.19,<2.0
   - scipy>=1.10
 
   # Bioinformatics (bioconda)

pixi.toml

@@ -27,7 +27,7 @@ bwa = ">=0.7"
 # Data processing
 numpy = ">=1.21,<2.0"
 pandas = ">=1.5,<3.0"
-polars = ">=0.19"
+polars = ">=0.19,<2.0"
 scipy = ">=1.10"
 pysam = ">=0.21"
 pybedtools = ">=0.9"

pyproject.toml

@@ -46,7 +46,7 @@ classifiers = [
 dependencies = [
     "numpy>=1.21.0",
     "pandas>=1.5.0,<3.0.0",
-    "polars>=0.19.0",
+    "polars>=0.19.0,<2.0",
     "scipy>=1.10.0",
     "pysam>=0.21.0",
     "pybedtools>=0.9.0",

requirements.txt

@@ -6,7 +6,7 @@
 # Data processing
 numpy>=1.21.0
 pandas>=1.5.0,<3.0.0
-polars>=0.19.0
+polars>=0.19.0,<2.0
 scipy>=1.10.0
 
 # Bioinformatics

rust/src/analysis.rs

@@ -508,9 +508,9 @@ fn golden_section_search<F>(f: F, a: f64, mut b: f64, tol: f64) -> Result<f64>
 where
     F: Fn(f64) -> f64,
 {
-    const PHI: f64 = 1.618033988749895; // Golden ratio
-    let inv_phi = 1.0 / PHI;
-    let inv_phi2 = 1.0 / (PHI * PHI);
+    let phi = (1.0 + 5.0_f64.sqrt()) / 2.0;
+    let inv_phi = 1.0 / phi;
+    let inv_phi2 = 1.0 / (phi * phi);
 
     let mut a = a;
     let mut h = b - a;
@@ -1211,18 +1211,20 @@ pub fn analyze_imbalance(
 
     // Remove pseudocounts from results
     for result in results.iter_mut() {
-        if result.ref_count < config.pseudocount
-            || result.alt_count < config.pseudocount
-            || result.n < 2 * config.pseudocount
+        // Each retained SNP contributes one pseudocount per allele.
+        let region_pseudocount = config.pseudocount.saturating_mul(result.snp_count as u32);
+        if result.ref_count < region_pseudocount
+            || result.alt_count < region_pseudocount
+            || result.n < 2 * region_pseudocount
         {
             eprintln!(
                 "[WARN] Counts smaller than pseudocount for region {}: ref={}, alt={}, n={}, pc={}",
-                result.region, result.ref_count, result.alt_count, result.n, config.pseudocount
+                result.region, result.ref_count, result.alt_count, result.n, region_pseudocount
             );
         }
-        result.ref_count = result.ref_count.saturating_sub(config.pseudocount);
-        result.alt_count = result.alt_count.saturating_sub(config.pseudocount);
-        result.n = result.n.saturating_sub(2 * config.pseudocount);
+        result.ref_count = result.ref_count.saturating_sub(region_pseudocount);
+        result.alt_count = result.alt_count.saturating_sub(region_pseudocount);
+        result.n = result.n.saturating_sub(2 * region_pseudocount);
     }
 
     Ok(results)
@@ -1232,6 +1234,40 @@ pub fn analyze_imbalance(
 mod tests {
     use super::*;
 
+    #[test]
+    fn test_region_counts_remove_pseudocounts_per_retained_snp() {
+        let variant = |pos, ref_count, alt_count| VariantCounts {
+            chrom: "chr1".to_string(),
+            pos,
+            ref_count,
+            alt_count,
+            region: "peak".to_string(),
+            sample: None,
+            gt: None,
+        };
+        for pseudocount in [0, 1, 2] {
+            // The duplicate must not contribute counts or extra pseudocounts.
+            let variants = vec![
+                variant(101, 18, 2),
+                variant(201, 3, 17),
+                variant(101, 18, 2),
+            ];
+            let config = AnalysisConfig {
+                min_count: 10,
+                pseudocount,
+                method: AnalysisMethod::Single,
+                phased: false,
+            };
+            let results = analyze_imbalance(variants, &config).unwrap();
+            assert_eq!(results.len(), 1);
+            let result = &results[0];
+            assert_eq!(result.snp_count, 2);
+            assert_eq!(result.ref_count, 21);
+            assert_eq!(result.alt_count, 19);
+            assert_eq!(result.n, 40);
+        }
+    }
+
     #[test]
     fn test_opt_prob() {
         // Test beta-binomial likelihood calculation

rust/src/bam_counter.rs

@@ -5,7 +5,7 @@ use rust_htslib::{bam, bam::ext::BamRecordExtensions, bam::Read as BamRead};
 use rustc_hash::{FxHashMap, FxHashSet};
 use std::path::Path;
 
-/// BAM allele counter using rust-htslib with batched fetching
+/// BAM allele counter preserving archived Python SNP encounter order.
 #[pyclass]
 pub struct BamCounter {
     bam_path: String,
@@ -43,7 +43,7 @@ impl BamCounter {
         Ok(BamCounter { bam_path })
     }
 
-    /// Count alleles at variant positions (SNPs and INDELs) using batched fetching
+    /// Count alleles at variant positions (SNPs and INDELs).
     ///
     /// Args:
     ///     regions: List of (chrom, pos, ref, alt) tuples - supports both SNPs and INDELs
@@ -172,6 +172,9 @@ impl BamCounter {
         min_qual: u8,
         debug_sites: &FxHashMap<(String, u32), usize>,
     ) -> PyResult<FxHashMap<usize, (u32, u32, u32)>> {
+        if regions.iter().all(|(_, region)| region.is_snp()) {
+            return self.process_chromosome_snps(chrom, regions, min_qual, debug_sites);
+        }
         let mut bam = bam::IndexedReader::from_path(&self.bam_path).map_err(|e| {
             PyErr::new::<pyo3::exceptions::PyIOError, _>(format!("Failed to open BAM: {}", e))
         })?;
@@ -374,6 +377,65 @@ impl BamCounter {
         Ok(counts)
     }
 
+    /// Match archive count_snp_alleles: SNP order takes priority over BAM order.
+    /// The seen-name set includes mates and persists across all SNPs on a chromosome.
+    fn process_chromosome_snps(
+        &self,
+        chrom: &str,
+        regions: &[(usize, Region)],
+        min_qual: u8,
+        debug_sites: &FxHashMap<(String, u32), usize>,
+    ) -> PyResult<FxHashMap<usize, (u32, u32, u32)>> {
+        let mut bam = bam::IndexedReader::from_path(&self.bam_path).map_err(|e| {
+            PyErr::new::<pyo3::exceptions::PyIOError, _>(format!("Failed to open BAM: {e}"))
+        })?;
+        let mut seen_reads: FxHashSet<Vec<u8>> = FxHashSet::default();
+        let mut counts = FxHashMap::default();
+        for (idx, region) in regions {
+            let target = region.pos.checked_sub(1).ok_or_else(|| {
+                PyErr::new::<pyo3::exceptions::PyValueError, _>("SNP positions must be one-based")
+            })?;
+            if let Err(e) = bam.fetch((chrom, i64::from(target), i64::from(region.pos))) {
+                eprintln!(
+                    "[WARN] Failed to fetch {chrom}:{}: {e}. Skipping remaining SNPs.",
+                    region.pos
+                );
+                break;
+            }
+            let mut site_counts = (0u32, 0u32, 0u32);
+            for record in bam.records() {
+                let record = record.map_err(|e| {
+                    PyErr::new::<pyo3::exceptions::PyIOError, _>(format!("Failed to read BAM: {e}"))
+                })?;
+                // Python consumes the name before looking for an aligned base,
+                // even when the fetched record has a deletion or reference skip here.
+                if !seen_reads.insert(record.qname().to_vec()) {
+                    continue;
+                }
+                if let Some(base) = get_base_at_position(&record, target, min_qual) {
+                    if base as u8 == region.ref_allele.as_bytes()[0] {
+                        site_counts.0 += 1;
+                    } else if base as u8 == region.alt_allele.as_bytes()[0] {
+                        site_counts.1 += 1;
+                    } else {
+                        site_counts.2 += 1;
+                    }
+                    if let Some(limit) = debug_sites.get(&(chrom.to_string(), region.pos)) {
+                        if site_counts.0 + site_counts.1 + site_counts.2 <= *limit as u32 {
+                            eprintln!(
+                                "[DEBUG VARIANT] {chrom}:{} read={} -> idx={idx} base={base}",
+                                region.pos,
+                                String::from_utf8_lossy(record.qname())
+                            );
+                        }
+                    }
+                }
+            }
+            counts.insert(*idx, site_counts);
+        }
+        Ok(counts)
+    }
+
     /// Group regions by chromosome while preserving encounter order
     fn group_regions_by_chrom(&self, regions: &[Region]) -> Vec<(String, Vec<(usize, Region)>)> {
         let mut grouped: Vec<Vec<(usize, Region)>> = Vec::new();
@@ -397,7 +459,6 @@ impl BamCounter {
 
 /// Get base at genomic position, accounting for CIGAR operations
 /// Matches WASP2 behavior: NO quality filtering by default
-#[allow(dead_code)]
 fn get_base_at_position(
     record: &bam::Record,
     target_pos: u32, // 0-based genomic position
@@ -423,14 +484,7 @@ fn get_base_at_position(
             }
 
             // Get the base (using array indexing)
-            let base = match seq[qpos as usize] {
-                b'A' => 'A',
-                b'C' => 'C',
-                b'G' => 'G',
-                b'T' => 'T',
-                b'N' => 'N',
-                _ => return None,
-            };
+            let base = char::from(seq[qpos as usize]);
             return Some(base);
         }
     }