From ce5e507d3d8bd513debc833be187d17df07fc570 Mon Sep 17 00:00:00 2001 From: Jeff Jaureguy <67065808+Jaureguy760@users.noreply.github.com> Date: Tue, 13 Oct 2026 12:14:29 -0700 Subject: [PATCH 1/3] fix: restore archived Python count and single-cell compatibility --- CHANGELOG.md | 12 + rust/src/analysis.rs | 50 +- scripts/compare_archive_analysis.py | 329 +++ src/analysis/as_analysis_sc.py | 8 +- src/analysis/compare_ai.py | 69 +- src/analysis/run_analysis_sc.py | 35 +- .../fixtures/archive_python_reference.json | 2630 +++++++++++++++++ tests/analysis/test_archive_parity.py | 173 ++ 8 files changed, 3243 insertions(+), 63 deletions(-) create mode 100644 scripts/compare_archive_analysis.py create mode 100644 tests/analysis/fixtures/archive_python_reference.json create mode 100644 tests/analysis/test_archive_parity.py diff --git a/CHANGELOG.md b/CHANGELOG.md index c950273..a0295dc 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -4,6 +4,18 @@ All notable changes to WASP2 will be documented in this file. ## [Unreleased] +### Fixed +- 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 diff --git a/rust/src/analysis.rs b/rust/src/analysis.rs index 73f6d60..8b639da 100644 --- a/rust/src/analysis.rs +++ b/rust/src/analysis.rs @@ -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 diff --git a/scripts/compare_archive_analysis.py b/scripts/compare_archive_analysis.py new file mode 100644 index 0000000..9059852 --- /dev/null +++ b/scripts/compare_archive_analysis.py @@ -0,0 +1,329 @@ +#!/usr/bin/env python3 +"""Compare analysis with Aaron Ho's pinned, pre-migration Python implementation. + +Run from a full repository clone with an installed Rust extension: + python scripts/compare_archive_analysis.py --output comparison.json +The archived code is read from git without altering it. AnnData is transposed +at the boundary to compare the historical SNP-by-cell layout with today's +cell-by-SNP layout. This is a compatibility check, not statistical calibration. +""" + +from __future__ import annotations + +import argparse +import contextlib +import io +import json +import subprocess +import sys +import tempfile +import types +import warnings +from pathlib import Path + +import numpy as np +import pandas as pd +from anndata import AnnData +from scipy.sparse import csr_matrix + +ROOT = Path(__file__).resolve().parents[1] +ARCHIVE_REF = "41354fc" +STAT_COLUMNS = ["null_ll", "alt_ll", "mu", "lrt", "pval", "fdr_pval"] +COUNT_COLUMNS = ["ref_count", "alt_count", "N", "snp_count", "num_snps"] + + +def load_archive(ref: str = ARCHIVE_REF) -> dict[str, types.ModuleType]: + """Load unmodified historical modules, restoring their legacy import names.""" + names = ["as_analysis", "as_analysis_sc", "compare_ai"] + saved = {name: sys.modules.get(name) for name in names} + modules = {} + try: + for name in names: + source = subprocess.check_output( + ["git", "show", f"{ref}:src/analysis/{name}.py"], cwd=ROOT, text=True + ) + module = types.ModuleType(f"archive_{name}") + modules[name] = module + sys.modules[name] = module + exec(compile(source, f"{ref}:src/analysis/{name}.py", "exec"), module.__dict__) + finally: + for name, module in saved.items(): + if module is None: + sys.modules.pop(name, None) + else: + sys.modules[name] = module + return modules + + +def adapt_linear_reference(archive: dict[str, types.ModuleType]) -> None: + """Explicitly repair archived array shapes, without changing its likelihood. + + The archive passes length-one arrays where SciPy requires scalar objectives, + and passes the full dispersion vector to the first/subsequent SNP slices. + This optional reference-only adapter aligns those shapes. It deliberately + retains the archived optimizer, first-SNP anchor, chi-square p-values, and + linear model's failure to forward its phased argument. + """ + module = archive["as_analysis"] + original_prob = module.opt_prob + + def scalar_prob(*args, **kwargs): + result = original_prob(*args, **kwargs) + return float(np.asarray(result).item()) if np.asarray(result).size == 1 else result + + def aligned_dp(prob, disp, first_ref, first_n, phase_ref, phase_n): + dispersions = np.asarray(disp) + per_snp = dispersions.ndim and len(dispersions) > 1 + first_disp = dispersions[0] if per_snp else disp + phase_disp = dispersions[1:] if per_snp else disp + first_ll = scalar_prob(prob, first_disp, first_ref[0], first_n[0]) + p1 = np.atleast_1d(scalar_prob(prob, phase_disp, phase_ref, phase_n, log=False)) + p2 = np.atleast_1d(scalar_prob(1 - prob, phase_disp, phase_ref, phase_n, log=False)) + return float(first_ll - np.log(np.prod(0.5 * p1 + 0.5 * p2))) + + module.opt_prob = scalar_prob + module.opt_unphased_dp = aligned_dp + + +def bulk_counts() -> pd.DataFrame: + pairs = [ + (18, 2), + (3, 17), + (22, 18), + (25, 15), + (29, 11), + (10, 30), + (10, 10), + (8, 12), + (35, 5), + (6, 34), + (14, 6), + (7, 13), + ] + return pd.DataFrame( + { + "chrom": ["chr1"] * len(pairs), + "pos": np.arange(1, len(pairs) + 1) * 100, + "ref": ["A"] * len(pairs), + "alt": ["G"] * len(pairs), + "GT": ["0|1", "1|0"] * (len(pairs) // 2), + "region": [f"peak{i // 2}" for i in range(len(pairs))], + "gene": [f"gene{i // 3}" for i in range(len(pairs))], + "ref_count": [p[0] for p in pairs], + "alt_count": [p[1] for p in pairs], + "other_count": [0] * len(pairs), + } + ) + + +def single_cell_counts() -> AnnData: + """Historical layout; includes both single- and multi-SNP features.""" + ref = np.array( + [[9, 8, 2, 3], [2, 3, 8, 9], [6, 5, 5, 4], [8, 7, 7, 8], [3, 2, 2, 3], [5, 4, 6, 5]], + dtype=np.uint16, + ) + alt = 10 - ref + obs = pd.DataFrame( + { + "index": np.arange(6), + "sample1": ["0|1", "1|0"] * 3, + "ref_count": ref.sum(axis=1), + "alt_count": alt.sum(axis=1), + }, + index=[f"snp{i}" for i in range(6)], + ) + var = pd.DataFrame({"group": ["A", "A", "B", "B"]}, index=[f"cell{i}" for i in range(4)]) + adata = AnnData(X=csr_matrix(ref + alt), obs=obs, var=var) + adata.layers["ref"] = csr_matrix(ref) + adata.layers["alt"] = csr_matrix(alt) + adata.uns["samples"] = ["sample1"] + adata.uns["feature"] = pd.DataFrame( + {"index": np.arange(6), "region": ["p0", "p0", "p1", "p2", "p2", "p2"]} + ) + return adata + + +def compare_frames(expected: pd.DataFrame, observed: pd.DataFrame, key: str = "region") -> dict: + expected = expected.sort_values(key).reset_index(drop=True) + observed = observed.sort_values(key).reset_index(drop=True) + if expected[key].tolist() != observed[key].tolist(): + return {"status": "FAIL", "reason": "different tested identities"} + deltas = {} + failures = [] + for column in [*COUNT_COLUMNS, *STAT_COLUMNS, "mu1", "mu2", "combined_mu"]: + if column not in expected: + continue + if column not in observed: + failures.append(f"missing {column}") + continue + left = expected[column].to_numpy(dtype=float) + right = observed[column].to_numpy(dtype=float) + deltas[column] = float(np.max(np.abs(left - right))) if len(left) else 0.0 + if column in COUNT_COLUMNS: + equal = np.array_equal(left, right) + else: + equal = ( + np.all(np.isfinite(left)) + and np.all(np.isfinite(right)) + and np.allclose(left, right, rtol=5e-5, atol=5e-5) + ) + if not equal: + failures.append(column) + for alpha in [0.05, 0.1]: + if not np.array_equal(expected["fdr_pval"] < alpha, observed["fdr_pval"] < alpha): + failures.append(f"significance at q={alpha}") + return { + "status": "FAIL" if failures else "PASS", + "failures": failures, + "max_absolute_differences": deltas, + } + + +def run_comparison(archive: dict[str, types.ModuleType]) -> dict: + import wasp2_rust + from analysis.as_analysis_sc import get_imbalance_sc + from analysis.compare_ai import get_compared_imbalance, get_compared_imbalance_diff_snps + from analysis.run_analysis import run_ai_analysis + from analysis.run_analysis_sc import process_adata_inputs + + cases = {} + with tempfile.TemporaryDirectory(prefix="wasp2-parity-") as tmp: + temp = Path(tmp) + for model in ["single", "linear"]: + for unit in ["snv", "atac", "rna"]: + for phased in [False, True]: + for pc in [0, 1]: + name = f"bulk/{model}/{unit}/phased={phased}/pc={pc}" + region = {"snv": None, "atac": "region", "rna": "gene"}[unit] + counts = bulk_counts().drop( + columns=["region", "gene"] + if unit == "snv" + else ["gene"] + if unit == "atac" + else ["region"] + ) + counts.to_csv(temp / "counts.tsv", sep="\t", index=False) + try: + expected = archive["as_analysis"].get_imbalance( + counts.copy(), + method=model, + phased=phased, + region_col=region, + pseudocount=pc, + ) + except Exception as exc: + cases[name] = {"status": "REFERENCE_ERROR", "reason": repr(exc)} + continue + try: + run_ai_analysis( + temp / "counts.tsv", + model=model, + phased=phased, + region_col=region, + pseudocount=pc, + out_file=str(temp / "results.tsv"), + ) + actual = pd.read_csv(temp / "results.tsv", sep="\t") + if region != "region": + expected = expected.rename(columns={region or "variant": "region"}) + cases[name] = compare_frames(expected, actual) + except Exception as exc: + cases[name] = {"status": "ERROR", "reason": repr(exc)} + for phased in [False, True]: + for pc in [0, 1]: + old = single_cell_counts() + new = old.T.copy() + options = { + "sample": "sample1", + "phased": phased, + "pseudocount": pc, + "min_count": 10, + } + for label, reference, current in [ + ("single-cell", archive["as_analysis_sc"].get_imbalance_sc, get_imbalance_sc), + ( + "differential/shared", + archive["compare_ai"].get_compared_imbalance, + get_compared_imbalance, + ), + ( + "differential/different", + archive["compare_ai"].get_compared_imbalance_diff_snps, + get_compared_imbalance_diff_snps, + ), + ]: + name = f"{label}/phased={phased}/pc={pc}" + try: + expected = reference(old.copy(), **options) + except Exception as exc: + cases[name] = {"status": "REFERENCE_ERROR", "reason": repr(exc)} + continue + try: + prepared = process_adata_inputs(new.copy(), sample="sample1", phased=phased) + actual = current(prepared.adata, **options) + details = { + str(key): compare_frames(frame, actual[key]) + for key, frame in expected.items() + } + cases[name] = { + "status": "PASS" + if expected.keys() == actual.keys() + and all(d["status"] == "PASS" for d in details.values()) + else "FAIL", + "groups": details, + } + except Exception as exc: + cases[name] = {"status": "ERROR", "reason": repr(exc)} + return { + "archive_commit": subprocess.check_output( + ["git", "rev-parse", ARCHIVE_REF], cwd=ROOT, text=True + ).strip(), + "current_commit": subprocess.check_output( + ["git", "rev-parse", "HEAD"], cwd=ROOT, text=True + ).strip(), + "extension": wasp2_rust.__file__, + "cases": cases, + } + + +def main() -> None: + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument("--output", type=Path, required=True) + parser.add_argument( + "--adapt-linear-reference", + action="store_true", + help="Explicitly repair archived scalar/array shapes; report this adaptation in JSON.", + ) + args = parser.parse_args() + # Import the installed extension before src's stub package can shadow it. + import wasp2_rust # noqa: F401 + + sys.path.insert(0, str(ROOT / "src")) + with warnings.catch_warnings(), contextlib.redirect_stdout(io.StringIO()): + warnings.simplefilter("ignore", FutureWarning) + archive = load_archive() + if args.adapt_linear_reference: + adapt_linear_reference(archive) + report = run_comparison(archive) + report["reference_adaptations"] = ( + ["Scalarize one-element objectives and align per-SNP dispersion arrays"] + if args.adapt_linear_reference + else [] + ) + report["known_reference_limitations"] = [ + "Archived linear_model does not forward phased to parse_opt; phased multi-SNP " + "linear comparisons therefore test different declared behavior." + ] + args.output.write_text(json.dumps(report, indent=2) + "\n") + summary = { + status: sum(c["status"] == status for c in report["cases"].values()) + for status in ["PASS", "FAIL", "ERROR", "REFERENCE_ERROR"] + } + print(json.dumps(summary)) + for name, result in report["cases"].items(): + if result["status"] != "PASS": + print(name, result.get("failures", result.get("reason", result.get("groups")))) + + +if __name__ == "__main__": + main() diff --git a/src/analysis/as_analysis_sc.py b/src/analysis/as_analysis_sc.py index faba877..e42ede4 100644 --- a/src/analysis/as_analysis_sc.py +++ b/src/analysis/as_analysis_sc.py @@ -121,8 +121,8 @@ def get_imbalance_sc( snp_cutoff = 2 * pseudocount # Pseudobulk per-SNP sums reduce over cells (axis=0) now that cells are rows. - ref_counts = adata.layers["ref"].sum(axis=0).T.A1 + pseudocount - alt_counts = adata.layers["alt"].sum(axis=0).T.A1 + pseudocount + ref_counts = adata.layers["ref"].sum(axis=0, dtype=np.uint64).T.A1 + pseudocount + alt_counts = adata.layers["alt"].sum(axis=0, dtype=np.uint64).T.A1 + pseudocount n_counts = ref_counts + alt_counts # Calculate dispersion across dataset @@ -152,8 +152,8 @@ def opt_disp(rho: float, ref_data: NDArray[np.uint16], n_data: NDArray[np.uint16 adata_sub = adata[adata.obs["group"] == group_name, :] # Create count data per group (per-SNP sums reduce over cells, axis=0) - ref_counts_group = adata_sub.layers["ref"].sum(axis=0).T.A1 + pseudocount - alt_counts_group = adata_sub.layers["alt"].sum(axis=0).T.A1 + pseudocount + ref_counts_group = adata_sub.layers["ref"].sum(axis=0, dtype=np.uint64).T.A1 + pseudocount + alt_counts_group = adata_sub.layers["alt"].sum(axis=0, dtype=np.uint64).T.A1 + pseudocount n_counts_group = ref_counts_group + alt_counts_group nonzero_idx = np.where(n_counts_group > snp_cutoff) # Get indices where counts were found diff --git a/src/analysis/compare_ai.py b/src/analysis/compare_ai.py index 60822c5..abc7b0a 100644 --- a/src/analysis/compare_ai.py +++ b/src/analysis/compare_ai.py @@ -15,6 +15,7 @@ # Local imports from .as_analysis import clamp_rho, opt_phased_new, opt_prob, opt_unphased_dp +from .as_analysis_sc import _assert_obs_cells_var_snps logger = logging.getLogger(__name__) @@ -100,12 +101,14 @@ def get_compared_imbalance( :param groups: List of groups to compare (if None, compare all) :return: Dict mapping (group1, group2) tuples to comparison DataFrames """ + _assert_obs_cells_var_snps(adata) + # Failsafe in case preparse somehow misses these if sample is None: phased = False if groups is None: - groups = list(adata.var["group"].dropna().unique()) + groups = list(adata.obs["group"].dropna().unique()) logger.info("Comparing all combinations of available groups") elif len(groups) == 1: raise ValueError("Please provide 2 or more groups to compare.") @@ -113,16 +116,16 @@ def get_compared_imbalance( # Process initial minimums for whole data dispersion snp_cutoff: int = min_count + (2 * pseudocount) - ref_counts: NDArray[np.uint16] = ( - adata.layers["ref"].sum(axis=1, dtype=np.uint16).T.A1 + pseudocount + ref_counts: NDArray[np.uint64] = ( + adata.layers["ref"].sum(axis=0, dtype=np.uint64).T.A1 + pseudocount ) - alt_counts: NDArray[np.uint16] = ( - adata.layers["alt"].sum(axis=1, dtype=np.uint16).T.A1 + pseudocount + alt_counts: NDArray[np.uint64] = ( + adata.layers["alt"].sum(axis=0, dtype=np.uint64).T.A1 + pseudocount ) - n_counts: NDArray[np.uint16] = ref_counts + alt_counts + n_counts: NDArray[np.uint64] = ref_counts + alt_counts # Calculate dispersion across dataset - def opt_disp(rho: float, ref_data: NDArray[np.uint16], n_data: NDArray[np.uint16]) -> float: + def opt_disp(rho: float, ref_data: NDArray[np.uint64], n_data: NDArray[np.uint64]) -> float: rho_safe = float(clamp_rho(rho)) # Prevent division by zero (Issue #228) return float( -np.sum( @@ -145,7 +148,7 @@ def opt_disp(rho: float, ref_data: NDArray[np.uint16], n_data: NDArray[np.uint16 if phased: gt_array: NDArray[np.uint8] | None = ( - adata.obs[sample].str.split("|", n=1).str[0].to_numpy(dtype=np.uint8) + adata.var[sample].str.split("|", n=1).str[0].to_numpy(dtype=np.uint8) ) else: gt_array = None @@ -156,10 +159,10 @@ def opt_disp(rho: float, ref_data: NDArray[np.uint16], n_data: NDArray[np.uint16 for group_name in groups: # Subset by group - adata_sub = adata[:, adata.var["group"] == group_name] + adata_sub = adata[adata.obs["group"] == group_name, :] - ref_counts_group = adata_sub.layers["ref"].sum(axis=1, dtype=np.uint16).T.A1 + pseudocount - alt_counts_group = adata_sub.layers["alt"].sum(axis=1, dtype=np.uint16).T.A1 + pseudocount + ref_counts_group = adata_sub.layers["ref"].sum(axis=0, dtype=np.uint64).T.A1 + pseudocount + alt_counts_group = adata_sub.layers["alt"].sum(axis=0, dtype=np.uint64).T.A1 + pseudocount n_counts_group = ref_counts_group + alt_counts_group nonzero_idx = np.where( @@ -247,10 +250,10 @@ def opt_disp(rho: float, ref_data: NDArray[np.uint16], n_data: NDArray[np.uint16 def compare_imbalance_between_groups( disp: float, - ref_counts1: NDArray[np.uint16], - n_counts1: NDArray[np.uint16], - ref_counts2: NDArray[np.uint16], - n_counts2: NDArray[np.uint16], + ref_counts1: NDArray[np.uint64], + n_counts1: NDArray[np.uint64], + ref_counts2: NDArray[np.uint64], + n_counts2: NDArray[np.uint64], region_snp_dict: dict[str, tuple[int, ...]], gt_array: NDArray[np.uint8] | None = None, ) -> pd.DataFrame: @@ -402,12 +405,14 @@ def get_compared_imbalance_diff_snps( :param groups: List of groups to compare (if None, compare all) :return: Dict mapping (group1, group2) tuples to comparison DataFrames """ + _assert_obs_cells_var_snps(adata) + # Failsafe in case preparse somehow misses these if sample is None: phased = False if groups is None: - groups = list(adata.var["group"].dropna().unique()) + groups = list(adata.obs["group"].dropna().unique()) logger.info("Comparing all combinations of available groups") elif len(groups) == 1: raise ValueError("Please provide 2 or more groups to compare.") @@ -415,25 +420,25 @@ def get_compared_imbalance_diff_snps( # Process initial minimums for whole data dispersion cutoff: int = min_count + (2 * pseudocount) - ref_counts: NDArray[np.uint16] = ( - adata.layers["ref"].sum(axis=1, dtype=np.uint16).T.A1 + pseudocount + ref_counts: NDArray[np.uint64] = ( + adata.layers["ref"].sum(axis=0, dtype=np.uint64).T.A1 + pseudocount ) - alt_counts: NDArray[np.uint16] = ( - adata.layers["alt"].sum(axis=1, dtype=np.uint16).T.A1 + pseudocount + alt_counts: NDArray[np.uint64] = ( + adata.layers["alt"].sum(axis=0, dtype=np.uint64).T.A1 + pseudocount ) - n_counts: NDArray[np.uint16] = ref_counts + alt_counts + n_counts: NDArray[np.uint64] = ref_counts + alt_counts min_idx: tuple[NDArray[np.intp], ...] = np.where( n_counts >= cutoff ) # Get indices for min_count - ref_counts_filt: NDArray[np.uint16] - n_counts_filt: NDArray[np.uint16] + ref_counts_filt: NDArray[np.uint64] + n_counts_filt: NDArray[np.uint64] ref_counts_filt, n_counts_filt = ref_counts[min_idx], n_counts[min_idx] # Calculate dispersion across dataset def opt_disp_filt( - rho: float, ref_data: NDArray[np.uint16], n_data: NDArray[np.uint16] + rho: float, ref_data: NDArray[np.uint64], n_data: NDArray[np.uint64] ) -> float: rho_safe = float(clamp_rho(rho)) # Prevent division by zero (Issue #228) return float( @@ -466,10 +471,10 @@ def opt_disp_filt( for group_name in groups: # Subset by group - adata_sub = adata[:, adata.var["group"] == group_name] + adata_sub = adata[adata.obs["group"] == group_name, :] - ref_counts_group = adata_sub.layers["ref"].sum(axis=1, dtype=np.uint16).T.A1 + pseudocount - alt_counts_group = adata_sub.layers["alt"].sum(axis=1, dtype=np.uint16).T.A1 + pseudocount + ref_counts_group = adata_sub.layers["ref"].sum(axis=0, dtype=np.uint64).T.A1 + pseudocount + alt_counts_group = adata_sub.layers["alt"].sum(axis=0, dtype=np.uint64).T.A1 + pseudocount n_counts_group = ref_counts_group + alt_counts_group min_idx_group = np.where(n_counts_group >= cutoff) # Get indices for min_count @@ -481,7 +486,7 @@ def opt_disp_filt( if phased: phase_array = ( - adata.obs.iloc[min_idx_group][sample] + adata.var.iloc[min_idx_group][sample] .str.split("|", n=1) .str[0] .to_numpy(dtype=np.uint8) @@ -526,12 +531,12 @@ def opt_disp_filt( def compare_imbalance_between_groups_diff_snps( disp: float, - ref_counts1: NDArray[np.uint16], - n_counts1: NDArray[np.uint16], + ref_counts1: NDArray[np.uint64], + n_counts1: NDArray[np.uint64], phase_array1: NDArray[np.uint8] | None, region_snp_dict1: dict[str, tuple[int, ...]], - ref_counts2: NDArray[np.uint16], - n_counts2: NDArray[np.uint16], + ref_counts2: NDArray[np.uint64], + n_counts2: NDArray[np.uint64], phase_array2: NDArray[np.uint8] | None, region_snp_dict2: dict[str, tuple[int, ...]], ) -> pd.DataFrame: diff --git a/src/analysis/run_analysis_sc.py b/src/analysis/run_analysis_sc.py index 3bf98f7..5a24212 100644 --- a/src/analysis/run_analysis_sc.py +++ b/src/analysis/run_analysis_sc.py @@ -8,7 +8,7 @@ from anndata import AnnData # local imports -from .as_analysis_sc import adata_count_qc, get_imbalance_sc +from .as_analysis_sc import _assert_obs_cells_var_snps, adata_count_qc, get_imbalance_sc # Class that stores relevant data @@ -82,7 +82,7 @@ def update_data(self, data: NamedTuple) -> None: # Define namedtuple for adata inputs class AdataInputs(NamedTuple): adata: AnnData - sample: str + sample: str | None groups: list[str] phased: bool @@ -103,6 +103,8 @@ def process_adata_inputs( groups = ai_files.groups if isinstance(ai_files.groups, list) else None phased = ai_files.phased + _assert_obs_cells_var_snps(adata) + # Check genotype and phasing input if "samples" not in adata.uns_keys(): if sample is not None: @@ -130,39 +132,33 @@ def process_adata_inputs( f"Please input a sample from list: {sample_list}" ) else: - if not any(i in ["1|0", "0|1", "1/0", "0/1"] for i in adata.obs[sample].unique()): + if not any(i in ["1|0", "0|1", "1/0", "0/1"] for i in adata.var[sample].unique()): raise ValueError(f"Heterozygous genotypes not found for sample: {sample}.") # Copy (not view) avoids implicit-modification warning downstream. - adata = adata[adata.obs[sample].isin(["1|0", "0|1", "1/0", "0/1"])].copy() - adata.obs = adata.obs.reset_index( - drop=True - ) # Have to reset index every time i subset adata + adata = adata[:, adata.var[sample].isin(["1|0", "0|1", "1/0", "0/1"])].copy() + adata.var = adata.var.reset_index(drop=True) # Have to reindex the regions after filtering GT's if "feature" in adata.uns_keys(): - # idx_df = adata.obs[["index"]].reset_index( - # drop=True).copy().reset_index(names="filt_index") - adata.uns["feature"] = ( adata.uns["feature"] - .merge(adata.obs[["index"]].reset_index(names="filt_index"), on="index")[ + .merge(adata.var[["index"]].reset_index(names="filt_index"), on="index")[ ["region", "filt_index"] ] .rename(columns={"filt_index": "index"}) ) - # Need to update adata.obs index col as well - adata.obs["index"] = adata.obs.index + adata.var["index"] = adata.var.index # Check phasing if True or None if phased is not False: - if {"0|1", "1|0"}.issubset(set(adata.obs[sample].unique())): + if {"0|1", "1|0"}.issubset(set(adata.var[sample].unique())): phased = True else: phased = False warning_msg = ( - f"Phased model selected for unphased genotypes ({adata.obs[sample].unique()}). " + f"Phased model selected for unphased genotypes ({adata.var[sample].unique()}). " "Switching to unphased model" ) warnings.warn(warning_msg, stacklevel=2) @@ -172,13 +168,13 @@ def process_adata_inputs( map_df = pd.read_csv( bc_map, sep="\t", header=None, names=["group"], index_col=0, dtype="category" ) - adata.var = adata.var.join(map_df, how="left") + adata.obs = adata.obs.drop(columns="group", errors="ignore").join(map_df, how="left") # No existing groups or mapping provided - if "group" not in adata.var_keys(): + if "group" not in adata.obs_keys(): raise KeyError("groups not found in dataset, please provide a barcode mapping") elif groups is not None: - valid_groups = list(adata.var["group"].dropna().unique()) + valid_groups = list(adata.obs["group"].dropna().unique()) new_groups = [i for i in groups if i in valid_groups] if len(new_groups) == 0: @@ -190,10 +186,9 @@ def process_adata_inputs( else: groups = new_groups else: - groups = list(adata.var["group"].dropna().unique()) + groups = list(adata.obs["group"].dropna().unique()) # Ensure all required values are set (type narrowing for mypy) - assert sample is not None, "sample must be set by this point" assert groups is not None, "groups must be set by this point" assert phased is not None, "phased must be set by this point" diff --git a/tests/analysis/fixtures/archive_python_reference.json b/tests/analysis/fixtures/archive_python_reference.json new file mode 100644 index 0000000..47efd5b --- /dev/null +++ b/tests/analysis/fixtures/archive_python_reference.json @@ -0,0 +1,2630 @@ +{ + "archive_ref": "41354fc", + "description": "Outputs from unmodified archived Python; synthetic count fixtures from compare_archive_analysis.py. No expected values are derived from Rust.", + "bulk": { + "snv/False/0": [ + { + "region": "chr1_100", + "ref_count": 18, + "alt_count": 2, + "N": 20, + "snp_count": 1, + "null_ll": -3.4070034111360483, + "alt_ll": -2.1513719236806397, + "mu": 0.8212618672179456, + "lrt": 2.5112629749108173, + "pval": 0.11303530907818209, + "fdr_pval": 0.5791591824660205 + }, + { + "region": "chr1_200", + "ref_count": 3, + "alt_count": 17, + "N": 20, + "snp_count": 1, + "null_ll": -3.1872407677079115, + "alt_ll": -2.34015911092653, + "mu": 0.2230666188870232, + "lrt": 1.694163313562763, + "pval": 0.19305306082200682, + "fdr_pval": 0.5791591824660205 + }, + { + "region": "chr1_300", + "ref_count": 22, + "alt_count": 18, + "N": 40, + "snp_count": 1, + "null_ll": -3.3583982378979256, + "alt_ll": -3.3444580108274016, + "mu": 0.5385183061467373, + "lrt": 0.027880454141048006, + "pval": 0.8673901828057538, + "fdr_pval": 0.9462438357880949 + }, + { + "region": "chr1_400", + "ref_count": 25, + "alt_count": 15, + "N": 40, + "snp_count": 1, + "null_ll": -3.4043794875799938, + "alt_ll": -3.3149450393612003, + "mu": 0.5964567253341145, + "lrt": 0.17886889643758686, + "pval": 0.6723471003398727, + "fdr_pval": 0.8934674479614134 + }, + { + "region": "chr1_500", + "ref_count": 29, + "alt_count": 11, + "N": 40, + "snp_count": 1, + "null_ll": -3.5394151282618713, + "alt_ll": -3.2267953963654117, + "mu": 0.6745400281019851, + "lrt": 0.6252394637929193, + "pval": 0.42910690592944667, + "fdr_pval": 0.8582138118588933 + }, + { + "region": "chr1_600", + "ref_count": 10, + "alt_count": 30, + "N": 40, + "snp_count": 1, + "null_ll": -3.5900754173386877, + "alt_ll": -3.1931065599727684, + "mu": 0.30565997614099266, + "lrt": 0.7939377147318387, + "pval": 0.3729120966279318, + "fdr_pval": 0.8582138118588933 + }, + { + "region": "chr1_700", + "ref_count": 10, + "alt_count": 10, + "N": 20, + "snp_count": 1, + "null_ll": -2.700399614229547, + "alt_ll": -2.700399614229547, + "mu": 0.49999999999999917, + "lrt": -0.0, + "pval": 1.0, + "fdr_pval": 1.0 + }, + { + "region": "chr1_800", + "ref_count": 8, + "alt_count": 12, + "N": 20, + "snp_count": 1, + "null_ll": -2.7318744134268096, + "alt_ll": -2.6787932086415145, + "mu": 0.4227218600053489, + "lrt": 0.1061624095705902, + "pval": 0.7445562066345112, + "fdr_pval": 0.8934674479614134 + }, + { + "region": "chr1_900", + "ref_count": 35, + "alt_count": 5, + "N": 40, + "snp_count": 1, + "null_ll": -4.017183097998134, + "alt_ll": -2.8920926518336545, + "mu": 0.7974670453615643, + "lrt": 2.2501808923289595, + "pval": 0.13359878439826106, + "fdr_pval": 0.5791591824660205 + }, + { + "region": "chr1_1000", + "ref_count": 6, + "alt_count": 34, + "N": 40, + "snp_count": 1, + "null_ll": -3.899242760709291, + "alt_ll": -2.978694573661204, + "mu": 0.224039033024206, + "lrt": 1.8410963740961748, + "pval": 0.1748224998559644, + "fdr_pval": 0.5791591824660205 + }, + { + "region": "chr1_1100", + "ref_count": 14, + "alt_count": 6, + "N": 20, + "snp_count": 1, + "null_ll": -2.8332364714183322, + "alt_ll": -2.607895165762728, + "mu": 0.6552542282691142, + "lrt": 0.45068261131120835, + "pval": 0.5020109712836058, + "fdr_pval": 0.8605902364861814 + }, + { + "region": "chr1_1200", + "ref_count": 7, + "alt_count": 13, + "N": 20, + "snp_count": 1, + "null_ll": -2.772776107320106, + "alt_ll": -2.650434300377675, + "mu": 0.38387824581387486, + "lrt": 0.24468361388486226, + "pval": 0.6208436393900352, + "fdr_pval": 0.8934674479614134 + } + ], + "snv/False/1": [ + { + "region": "chr1_100", + "ref_count": 18, + "alt_count": 2, + "N": 20, + "snp_count": 1, + "null_ll": -3.44965243216541, + "alt_ll": -2.3048139074687706, + "mu": 0.7989369847225031, + "lrt": 2.289677049393279, + "pval": 0.13023701945544647, + "fdr_pval": 0.6147556441484411 + }, + { + "region": "chr1_200", + "ref_count": 3, + "alt_count": 17, + "N": 20, + "snp_count": 1, + "null_ll": -3.232319597978043, + "alt_ll": -2.4288356545879752, + "mu": 0.2403029858365358, + "lrt": 1.6069678867801356, + "pval": 0.20491854804948037, + "fdr_pval": 0.6147556441484411 + }, + { + "region": "chr1_300", + "ref_count": 22, + "alt_count": 18, + "N": 40, + "snp_count": 1, + "null_ll": -3.31439686754961, + "alt_ll": -3.2987916473063317, + "mu": 0.538295630876192, + "lrt": 0.0312104404865563, + "pval": 0.8597717160631849, + "fdr_pval": 0.937932781159838 + }, + { + "region": "chr1_400", + "ref_count": 25, + "alt_count": 15, + "N": 40, + "snp_count": 1, + "null_ll": -3.3723930751946427, + "alt_ll": -3.272580559403795, + "mu": 0.5958403200402597, + "lrt": 0.19962503158169564, + "pval": 0.6550236800869007, + "fdr_pval": 0.8733649067825342 + }, + { + "region": "chr1_500", + "ref_count": 29, + "alt_count": 11, + "N": 40, + "snp_count": 1, + "null_ll": -3.5408626649084582, + "alt_ll": -3.1952607748953668, + "mu": 0.6730900224097598, + "lrt": 0.691203780026183, + "pval": 0.40575523209466546, + "fdr_pval": 0.8115104641893309 + }, + { + "region": "chr1_600", + "ref_count": 10, + "alt_count": 30, + "N": 40, + "snp_count": 1, + "null_ll": -3.603307603483897, + "alt_ll": -3.1661170855428624, + "mu": 0.30742429956570055, + "lrt": 0.8743810358820694, + "pval": 0.3497453011477637, + "fdr_pval": 0.8115104641893309 + }, + { + "region": "chr1_700", + "ref_count": 10, + "alt_count": 10, + "N": 20, + "snp_count": 1, + "null_ll": -2.7016027910294853, + "alt_ll": -2.7016027910294853, + "mu": 0.5, + "lrt": -0.0, + "pval": 1.0, + "fdr_pval": 1.0 + }, + { + "region": "chr1_800", + "ref_count": 8, + "alt_count": 12, + "N": 20, + "snp_count": 1, + "null_ll": -2.737729196484674, + "alt_ll": -2.684067065544477, + "mu": 0.42674217074194676, + "lrt": 0.10732426188039401, + "pval": 0.7432112344153309, + "fdr_pval": 0.891853481298397 + }, + { + "region": "chr1_900", + "ref_count": 35, + "alt_count": 5, + "N": 40, + "snp_count": 1, + "null_ll": -4.110468505481899, + "alt_ll": -2.9173309425077507, + "mu": 0.7924860660266098, + "lrt": 2.3862751259482966, + "pval": 0.12240497843235777, + "fdr_pval": 0.6147556441484411 + }, + { + "region": "chr1_1000", + "ref_count": 6, + "alt_count": 34, + "N": 40, + "snp_count": 1, + "null_ll": -3.974184081191121, + "alt_ll": -2.9865857759309193, + "mu": 0.22801198005748452, + "lrt": 1.9751966105204035, + "pval": 0.15989734974073055, + "fdr_pval": 0.6147556441484411 + }, + { + "region": "chr1_1100", + "ref_count": 14, + "alt_count": 6, + "N": 20, + "snp_count": 1, + "null_ll": -2.852486330351612, + "alt_ll": -2.62749785183385, + "mu": 0.6468920377476061, + "lrt": 0.4499769570355241, + "pval": 0.5023458972991031, + "fdr_pval": 0.8611643953698909 + }, + { + "region": "chr1_1200", + "ref_count": 7, + "alt_count": 13, + "N": 20, + "snp_count": 1, + "null_ll": -2.784331711827282, + "alt_ll": -2.661257951381625, + "mu": 0.39000229918098056, + "lrt": 0.24614752089131375, + "pval": 0.6198008856418893, + "fdr_pval": 0.8733649067825342 + } + ], + "snv/True/0": [ + { + "region": "chr1_100", + "ref_count": 18, + "alt_count": 2, + "N": 20, + "snp_count": 1, + "null_ll": -3.4070034111360483, + "alt_ll": -2.1513719236806397, + "mu": 0.8212618672179456, + "lrt": 2.5112629749108173, + "pval": 0.11303530907818209, + "fdr_pval": 0.5791591824660205 + }, + { + "region": "chr1_200", + "ref_count": 3, + "alt_count": 17, + "N": 20, + "snp_count": 1, + "null_ll": -3.1872407677079115, + "alt_ll": -2.34015911092653, + "mu": 0.2230666188870232, + "lrt": 1.694163313562763, + "pval": 0.19305306082200682, + "fdr_pval": 0.5791591824660205 + }, + { + "region": "chr1_300", + "ref_count": 22, + "alt_count": 18, + "N": 40, + "snp_count": 1, + "null_ll": -3.3583982378979256, + "alt_ll": -3.3444580108274016, + "mu": 0.5385183061467373, + "lrt": 0.027880454141048006, + "pval": 0.8673901828057538, + "fdr_pval": 0.9462438357880949 + }, + { + "region": "chr1_400", + "ref_count": 25, + "alt_count": 15, + "N": 40, + "snp_count": 1, + "null_ll": -3.4043794875799938, + "alt_ll": -3.3149450393612003, + "mu": 0.5964567253341145, + "lrt": 0.17886889643758686, + "pval": 0.6723471003398727, + "fdr_pval": 0.8934674479614134 + }, + { + "region": "chr1_500", + "ref_count": 29, + "alt_count": 11, + "N": 40, + "snp_count": 1, + "null_ll": -3.5394151282618713, + "alt_ll": -3.2267953963654117, + "mu": 0.6745400281019851, + "lrt": 0.6252394637929193, + "pval": 0.42910690592944667, + "fdr_pval": 0.8582138118588933 + }, + { + "region": "chr1_600", + "ref_count": 10, + "alt_count": 30, + "N": 40, + "snp_count": 1, + "null_ll": -3.5900754173386877, + "alt_ll": -3.1931065599727684, + "mu": 0.30565997614099266, + "lrt": 0.7939377147318387, + "pval": 0.3729120966279318, + "fdr_pval": 0.8582138118588933 + }, + { + "region": "chr1_700", + "ref_count": 10, + "alt_count": 10, + "N": 20, + "snp_count": 1, + "null_ll": -2.700399614229547, + "alt_ll": -2.700399614229547, + "mu": 0.49999999999999917, + "lrt": -0.0, + "pval": 1.0, + "fdr_pval": 1.0 + }, + { + "region": "chr1_800", + "ref_count": 8, + "alt_count": 12, + "N": 20, + "snp_count": 1, + "null_ll": -2.7318744134268096, + "alt_ll": -2.6787932086415145, + "mu": 0.4227218600053489, + "lrt": 0.1061624095705902, + "pval": 0.7445562066345112, + "fdr_pval": 0.8934674479614134 + }, + { + "region": "chr1_900", + "ref_count": 35, + "alt_count": 5, + "N": 40, + "snp_count": 1, + "null_ll": -4.017183097998134, + "alt_ll": -2.8920926518336545, + "mu": 0.7974670453615643, + "lrt": 2.2501808923289595, + "pval": 0.13359878439826106, + "fdr_pval": 0.5791591824660205 + }, + { + "region": "chr1_1000", + "ref_count": 6, + "alt_count": 34, + "N": 40, + "snp_count": 1, + "null_ll": -3.899242760709291, + "alt_ll": -2.978694573661204, + "mu": 0.224039033024206, + "lrt": 1.8410963740961748, + "pval": 0.1748224998559644, + "fdr_pval": 0.5791591824660205 + }, + { + "region": "chr1_1100", + "ref_count": 14, + "alt_count": 6, + "N": 20, + "snp_count": 1, + "null_ll": -2.8332364714183322, + "alt_ll": -2.607895165762728, + "mu": 0.6552542282691142, + "lrt": 0.45068261131120835, + "pval": 0.5020109712836058, + "fdr_pval": 0.8605902364861814 + }, + { + "region": "chr1_1200", + "ref_count": 7, + "alt_count": 13, + "N": 20, + "snp_count": 1, + "null_ll": -2.772776107320106, + "alt_ll": -2.650434300377675, + "mu": 0.38387824581387486, + "lrt": 0.24468361388486226, + "pval": 0.6208436393900352, + "fdr_pval": 0.8934674479614134 + } + ], + "snv/True/1": [ + { + "region": "chr1_100", + "ref_count": 18, + "alt_count": 2, + "N": 20, + "snp_count": 1, + "null_ll": -3.44965243216541, + "alt_ll": -2.3048139074687706, + "mu": 0.7989369847225031, + "lrt": 2.289677049393279, + "pval": 0.13023701945544647, + "fdr_pval": 0.6147556441484411 + }, + { + "region": "chr1_200", + "ref_count": 3, + "alt_count": 17, + "N": 20, + "snp_count": 1, + "null_ll": -3.232319597978043, + "alt_ll": -2.4288356545879752, + "mu": 0.2403029858365358, + "lrt": 1.6069678867801356, + "pval": 0.20491854804948037, + "fdr_pval": 0.6147556441484411 + }, + { + "region": "chr1_300", + "ref_count": 22, + "alt_count": 18, + "N": 40, + "snp_count": 1, + "null_ll": -3.31439686754961, + "alt_ll": -3.2987916473063317, + "mu": 0.538295630876192, + "lrt": 0.0312104404865563, + "pval": 0.8597717160631849, + "fdr_pval": 0.937932781159838 + }, + { + "region": "chr1_400", + "ref_count": 25, + "alt_count": 15, + "N": 40, + "snp_count": 1, + "null_ll": -3.3723930751946427, + "alt_ll": -3.272580559403795, + "mu": 0.5958403200402597, + "lrt": 0.19962503158169564, + "pval": 0.6550236800869007, + "fdr_pval": 0.8733649067825342 + }, + { + "region": "chr1_500", + "ref_count": 29, + "alt_count": 11, + "N": 40, + "snp_count": 1, + "null_ll": -3.5408626649084582, + "alt_ll": -3.1952607748953668, + "mu": 0.6730900224097598, + "lrt": 0.691203780026183, + "pval": 0.40575523209466546, + "fdr_pval": 0.8115104641893309 + }, + { + "region": "chr1_600", + "ref_count": 10, + "alt_count": 30, + "N": 40, + "snp_count": 1, + "null_ll": -3.603307603483897, + "alt_ll": -3.1661170855428624, + "mu": 0.30742429956570055, + "lrt": 0.8743810358820694, + "pval": 0.3497453011477637, + "fdr_pval": 0.8115104641893309 + }, + { + "region": "chr1_700", + "ref_count": 10, + "alt_count": 10, + "N": 20, + "snp_count": 1, + "null_ll": -2.7016027910294853, + "alt_ll": -2.7016027910294853, + "mu": 0.5, + "lrt": -0.0, + "pval": 1.0, + "fdr_pval": 1.0 + }, + { + "region": "chr1_800", + "ref_count": 8, + "alt_count": 12, + "N": 20, + "snp_count": 1, + "null_ll": -2.737729196484674, + "alt_ll": -2.684067065544477, + "mu": 0.42674217074194676, + "lrt": 0.10732426188039401, + "pval": 0.7432112344153309, + "fdr_pval": 0.891853481298397 + }, + { + "region": "chr1_900", + "ref_count": 35, + "alt_count": 5, + "N": 40, + "snp_count": 1, + "null_ll": -4.110468505481899, + "alt_ll": -2.9173309425077507, + "mu": 0.7924860660266098, + "lrt": 2.3862751259482966, + "pval": 0.12240497843235777, + "fdr_pval": 0.6147556441484411 + }, + { + "region": "chr1_1000", + "ref_count": 6, + "alt_count": 34, + "N": 40, + "snp_count": 1, + "null_ll": -3.974184081191121, + "alt_ll": -2.9865857759309193, + "mu": 0.22801198005748452, + "lrt": 1.9751966105204035, + "pval": 0.15989734974073055, + "fdr_pval": 0.6147556441484411 + }, + { + "region": "chr1_1100", + "ref_count": 14, + "alt_count": 6, + "N": 20, + "snp_count": 1, + "null_ll": -2.852486330351612, + "alt_ll": -2.62749785183385, + "mu": 0.6468920377476061, + "lrt": 0.4499769570355241, + "pval": 0.5023458972991031, + "fdr_pval": 0.8611643953698909 + }, + { + "region": "chr1_1200", + "ref_count": 7, + "alt_count": 13, + "N": 20, + "snp_count": 1, + "null_ll": -2.784331711827282, + "alt_ll": -2.661257951381625, + "mu": 0.39000229918098056, + "lrt": 0.24614752089131375, + "pval": 0.6198008856418893, + "fdr_pval": 0.8733649067825342 + } + ], + "atac/False/0": [ + { + "region": "peak0", + "ref_count": 21, + "alt_count": 19, + "N": 40, + "snp_count": 2, + "null_ll": -6.59424417884396, + "alt_ll": -5.165810742938646, + "mu": 0.7943947354422328, + "lrt": 2.8568668718106274, + "pval": 0.0909845596630293, + "fdr_pval": 0.28991576343191233 + }, + { + "region": "peak1", + "ref_count": 47, + "alt_count": 33, + "N": 80, + "snp_count": 2, + "null_ll": -6.762777725477919, + "alt_ll": -6.755094501416779, + "mu": 0.5212521049143231, + "lrt": 0.015366448122280829, + "pval": 0.9013457755508965, + "fdr_pval": 0.9999999524423742 + }, + { + "region": "peak2", + "ref_count": 39, + "alt_count": 41, + "N": 80, + "snp_count": 2, + "null_ll": -7.129490545600559, + "alt_ll": -6.871329403787703, + "mu": 0.6405801423804764, + "lrt": 0.516322283625712, + "pval": 0.472414721379269, + "fdr_pval": 0.9144978736886265 + }, + { + "region": "peak3", + "ref_count": 18, + "alt_count": 22, + "N": 40, + "snp_count": 2, + "null_ll": -5.432274027656357, + "alt_ll": -5.432274027656355, + "mu": 0.499999999999999, + "lrt": 3.552713678800501e-15, + "pval": 0.9999999524423742, + "fdr_pval": 0.9999999524423742 + }, + { + "region": "peak4", + "ref_count": 41, + "alt_count": 39, + "N": 80, + "snp_count": 2, + "null_ll": -7.916425858707425, + "alt_ll": -6.536345991339368, + "mu": 0.7814279138526817, + "lrt": 2.7601597347361153, + "pval": 0.09663858781063744, + "fdr_pval": 0.28991576343191233 + }, + { + "region": "peak5", + "ref_count": 21, + "alt_count": 19, + "N": 40, + "snp_count": 2, + "null_ll": -5.606012578738438, + "alt_ll": -5.475681238145786, + "mu": 0.5913652469015125, + "lrt": 0.26066268118530544, + "pval": 0.6096652491257509, + "fdr_pval": 0.9144978736886265 + } + ], + "atac/False/1": [ + { + "region": "peak0", + "ref_count": 21, + "alt_count": 19, + "N": 40, + "snp_count": 2, + "null_ll": -6.681972030143453, + "alt_ll": -5.398518037035643, + "mu": 0.772913502606682, + "lrt": 2.56690798621562, + "pval": 0.10912083165044163, + "fdr_pval": 0.3273624949513249 + }, + { + "region": "peak1", + "ref_count": 47, + "alt_count": 33, + "N": 80, + "snp_count": 2, + "null_ll": -6.6867899427442525, + "alt_ll": -6.678089819457147, + "mu": 0.521371033561547, + "lrt": 0.017400246574210243, + "pval": 0.8950554834685008, + "fdr_pval": 1.0 + }, + { + "region": "peak2", + "ref_count": 39, + "alt_count": 41, + "N": 80, + "snp_count": 2, + "null_ll": -7.144170268392355, + "alt_ll": -6.844449963755999, + "mu": 0.6442901453856908, + "lrt": 0.5994406092727136, + "pval": 0.4387915386081115, + "fdr_pval": 0.877583077216223 + }, + { + "region": "peak3", + "ref_count": 18, + "alt_count": 22, + "N": 40, + "snp_count": 2, + "null_ll": -5.4393319875141595, + "alt_ll": -5.43933198751416, + "mu": 0.5000000000000001, + "lrt": -1.7763568394002505e-15, + "pval": 1.0, + "fdr_pval": 1.0 + }, + { + "region": "peak4", + "ref_count": 41, + "alt_count": 39, + "N": 80, + "snp_count": 2, + "null_ll": -8.08465258667302, + "alt_ll": -6.577092239063196, + "mu": 0.7781676485131456, + "lrt": 3.0151206952196468, + "pval": 0.08249131410228634, + "fdr_pval": 0.3273624949513249 + }, + { + "region": "peak5", + "ref_count": 21, + "alt_count": 19, + "N": 40, + "snp_count": 2, + "null_ll": -5.636818042178894, + "alt_ll": -5.506932579076578, + "mu": 0.5860473268196926, + "lrt": 0.2597709262046326, + "pval": 0.610277575224037, + "fdr_pval": 0.9154163628360555 + } + ], + "atac/True/0": [ + { + "region": "peak0", + "ref_count": 21, + "alt_count": 19, + "N": 40, + "snp_count": 2, + "null_ll": -6.59424417884396, + "alt_ll": -4.5090968785153365, + "mu": 0.8003311941197009, + "lrt": 4.170294600657247, + "pval": 0.041138645929138976, + "fdr_pval": 0.12994807079238024 + }, + { + "region": "peak1", + "ref_count": 47, + "alt_count": 33, + "N": 80, + "snp_count": 2, + "null_ll": -6.762777725477919, + "alt_ll": -6.746218242590604, + "mu": 0.4702584390542891, + "lrt": 0.03311896577463003, + "pval": 0.855593580045347, + "fdr_pval": 0.855593580045347 + }, + { + "region": "peak2", + "ref_count": 39, + "alt_count": 41, + "N": 80, + "snp_count": 2, + "null_ll": -7.129490545600559, + "alt_ll": -6.422209788150893, + "mu": 0.6845625456249463, + "lrt": 1.414561514899333, + "pval": 0.2343006293423302, + "fdr_pval": 0.4686012586846604 + }, + { + "region": "peak3", + "ref_count": 18, + "alt_count": 22, + "N": 40, + "snp_count": 2, + "null_ll": -5.432274027656357, + "alt_ll": -5.405661547184883, + "mu": 0.538945255942631, + "lrt": 0.053224960942948485, + "pval": 0.8175437128606492, + "fdr_pval": 0.855593580045347 + }, + { + "region": "peak4", + "ref_count": 41, + "alt_count": 39, + "N": 80, + "snp_count": 2, + "null_ll": -7.916425858707425, + "alt_ll": -5.874912929583036, + "mu": 0.787014880525805, + "lrt": 4.08302585824878, + "pval": 0.043316023597460084, + "fdr_pval": 0.12994807079238024 + }, + { + "region": "peak5", + "ref_count": 21, + "alt_count": 19, + "N": 40, + "snp_count": 2, + "null_ll": -5.606012578738438, + "alt_ll": -5.265877171943071, + "mu": 0.6359861491008233, + "lrt": 0.6802708135907345, + "pval": 0.40949348604417024, + "fdr_pval": 0.6142402290662554 + } + ], + "atac/True/1": [ + { + "region": "peak0", + "ref_count": 21, + "alt_count": 19, + "N": 40, + "snp_count": 2, + "null_ll": -6.681972030143453, + "alt_ll": -4.747055894259569, + "mu": 0.7801213818657149, + "lrt": 3.869832271767768, + "pval": 0.0491614341908432, + "fdr_pval": 0.1474843025725296 + }, + { + "region": "peak1", + "ref_count": 47, + "alt_count": 33, + "N": 80, + "snp_count": 2, + "null_ll": -6.6867899427442525, + "alt_ll": -6.668366943333292, + "mu": 0.4705220111943598, + "lrt": 0.036845998821920745, + "pval": 0.8477789302581967, + "fdr_pval": 0.8477789302581967 + }, + { + "region": "peak2", + "ref_count": 39, + "alt_count": 41, + "N": 80, + "snp_count": 2, + "null_ll": -7.144170268392355, + "alt_ll": -6.363867286815808, + "mu": 0.6829426576911672, + "lrt": 1.560605963153094, + "pval": 0.21157651767291047, + "fdr_pval": 0.42315303534582094 + }, + { + "region": "peak3", + "ref_count": 18, + "alt_count": 22, + "N": 40, + "snp_count": 2, + "null_ll": -5.4393319875141595, + "alt_ll": -5.41244260857656, + "mu": 0.5368716864477697, + "lrt": 0.05377875787519848, + "pval": 0.8166137622160525, + "fdr_pval": 0.8477789302581967 + }, + { + "region": "peak4", + "ref_count": 41, + "alt_count": 39, + "N": 80, + "snp_count": 2, + "null_ll": -8.08465258667302, + "alt_ll": -5.907917118477915, + "mu": 0.7824879448667197, + "lrt": 4.3534709363902095, + "pval": 0.03693360053204512, + "fdr_pval": 0.1474843025725296 + }, + { + "region": "peak5", + "ref_count": 21, + "alt_count": 19, + "N": 40, + "snp_count": 2, + "null_ll": -5.636818042178894, + "alt_ll": -5.296181783547893, + "mu": 0.6286760026076096, + "lrt": 0.681272517262002, + "pval": 0.40914888153659645, + "fdr_pval": 0.6137233223048947 + } + ], + "rna/False/0": [ + { + "region": "gene0", + "ref_count": 43, + "alt_count": 37, + "N": 80, + "snp_count": 3, + "null_ll": -9.952642416741885, + "alt_ll": -9.13170055355274, + "mu": 0.7046744212826972, + "lrt": 1.641883726378289, + "pval": 0.2000672092033727, + "fdr_pval": 0.6976199949899382 + }, + { + "region": "gene1", + "ref_count": 64, + "alt_count": 56, + "N": 120, + "snp_count": 3, + "null_ll": -10.533870033180552, + "alt_ll": -10.475158348133196, + "mu": 0.5629361767249976, + "lrt": 0.11742337009471271, + "pval": 0.7318460303324834, + "fdr_pval": 0.9757947071099778 + }, + { + "region": "gene2", + "ref_count": 53, + "alt_count": 27, + "N": 80, + "snp_count": 3, + "null_ll": -9.44945712565449, + "alt_ll": -9.4494571256545, + "mu": 0.49999999999999883, + "lrt": -2.1316282072803006e-14, + "pval": 1.0, + "fdr_pval": 1.0 + }, + { + "region": "gene3", + "ref_count": 27, + "alt_count": 53, + "N": 80, + "snp_count": 3, + "null_ll": -9.505255339447729, + "alt_ll": -9.066364531324112, + "mu": 0.35947109819317846, + "lrt": 0.8777816162472334, + "pval": 0.3488099974949691, + "fdr_pval": 0.6976199949899382 + } + ], + "rna/False/1": [ + { + "region": "gene0", + "ref_count": 43, + "alt_count": 37, + "N": 80, + "snp_count": 3, + "null_ll": -9.996368897693063, + "alt_ll": -9.273712969605496, + "mu": 0.6836629820181194, + "lrt": 1.4453118561751346, + "pval": 0.2292816970922187, + "fdr_pval": 0.6666380765136606 + }, + { + "region": "gene1", + "ref_count": 64, + "alt_count": 56, + "N": 120, + "snp_count": 3, + "null_ll": -10.516563343586999, + "alt_ll": -10.44473241543556, + "mu": 0.5679960287587807, + "lrt": 0.1436618563028773, + "pval": 0.7046674324063067, + "fdr_pval": 0.9395565765417422 + }, + { + "region": "gene2", + "ref_count": 53, + "alt_count": 27, + "N": 80, + "snp_count": 3, + "null_ll": -9.54980049299606, + "alt_ll": -9.549800492996058, + "mu": 0.5000000000000002, + "lrt": 3.552713678800501e-15, + "pval": 0.9999999524423742, + "fdr_pval": 0.9999999524423742 + }, + { + "region": "gene3", + "ref_count": 27, + "alt_count": 53, + "N": 80, + "snp_count": 3, + "null_ll": -9.611002123370016, + "alt_ll": -9.143022204066911, + "mu": 0.36311633131200055, + "lrt": 0.9359598386062089, + "pval": 0.3333190382568303, + "fdr_pval": 0.6666380765136606 + } + ], + "rna/True/0": [ + { + "region": "gene0", + "ref_count": 43, + "alt_count": 37, + "N": 80, + "snp_count": 3, + "null_ll": -9.952642416741885, + "alt_ll": -8.381118647873514, + "mu": 0.7280834106829174, + "lrt": 3.1430475377367415, + "pval": 0.07625120901200119, + "fdr_pval": 0.27858243792718723 + }, + { + "region": "gene1", + "ref_count": 64, + "alt_count": 56, + "N": 120, + "snp_count": 3, + "null_ll": -10.533870033180552, + "alt_ll": -10.261142539881286, + "mu": 0.40233617953597134, + "lrt": 0.5454549865985321, + "pval": 0.4601807540346815, + "fdr_pval": 0.4601807540346815 + }, + { + "region": "gene2", + "ref_count": 53, + "alt_count": 27, + "N": 80, + "snp_count": 3, + "null_ll": -9.44945712565449, + "alt_ll": -8.844412633359365, + "mu": 0.6463826801555765, + "lrt": 1.21008898459025, + "pval": 0.27131449943905395, + "fdr_pval": 0.36175266591873856 + }, + { + "region": "gene3", + "ref_count": 27, + "alt_count": 53, + "N": 80, + "snp_count": 3, + "null_ll": -9.505255339447729, + "alt_ll": -8.412369407891662, + "mu": 0.30871553168813537, + "lrt": 2.1857718631121337, + "pval": 0.13929121896359362, + "fdr_pval": 0.27858243792718723 + } + ], + "rna/True/1": [ + { + "region": "gene0", + "ref_count": 43, + "alt_count": 37, + "N": 80, + "snp_count": 3, + "null_ll": -9.996368897693063, + "alt_ll": -8.533233552783928, + "mu": 0.7096022032957561, + "lrt": 2.9262706898182707, + "pval": 0.08714858999658096, + "fdr_pval": 0.2633404423264745 + }, + { + "region": "gene1", + "ref_count": 64, + "alt_count": 56, + "N": 120, + "snp_count": 3, + "null_ll": -10.516563343586999, + "alt_ll": -10.217770980958651, + "mu": 0.40381770662443445, + "lrt": 0.597584725256695, + "pval": 0.4395010494688114, + "fdr_pval": 0.4395010494688114 + }, + { + "region": "gene2", + "ref_count": 53, + "alt_count": 27, + "N": 80, + "snp_count": 3, + "null_ll": -9.54980049299606, + "alt_ll": -8.916531377779945, + "mu": 0.6414319846009774, + "lrt": 1.2665382304322286, + "pval": 0.2604171114729931, + "fdr_pval": 0.3472228152973241 + }, + { + "region": "gene3", + "ref_count": 27, + "alt_count": 53, + "N": 80, + "snp_count": 3, + "null_ll": -9.611002123370016, + "alt_ll": -8.474651467232219, + "mu": 0.31514336109738306, + "lrt": 2.2727013122755935, + "pval": 0.13167022116323726, + "fdr_pval": 0.2633404423264745 + } + ] + }, + "single_cell": { + "single/False/0": { + "A": [ + { + "region": "p0", + "num_snps": 2, + "mu": 0.7825601441901614, + "null_ll": -7.424952555448208, + "alt_ll": -4.696477217962663, + "pval": 0.019490631620863782, + "fdr_pval": 0.05847189486259134 + }, + { + "region": "p1", + "num_snps": 1, + "mu": 0.55, + "null_ll": -2.2092614298154665, + "alt_ll": -2.1683943841147517, + "pval": 0.774961034462935, + "fdr_pval": 0.774961034462935 + }, + { + "region": "p2", + "num_snps": 3, + "mu": 0.6569575554611818, + "null_ll": -8.50999611779697, + "alt_ll": -7.618515033354983, + "pval": 0.1817860094595569, + "fdr_pval": 0.2726790141893353 + } + ], + "B": [ + { + "region": "p0", + "num_snps": 2, + "mu": 0.2166705838515093, + "null_ll": -7.424952555448208, + "alt_ll": -4.701281319069672, + "pval": 0.019598117265708728, + "fdr_pval": 0.05879435179712618 + }, + { + "region": "p1", + "num_snps": 1, + "mu": 0.45, + "null_ll": -2.2092614298154665, + "alt_ll": -2.1683943841147535, + "pval": 0.7749610344629396, + "fdr_pval": 0.7749610344629396 + }, + { + "region": "p2", + "num_snps": 3, + "mu": 0.6569575554611818, + "null_ll": -8.50999611779697, + "alt_ll": -7.618515033354983, + "pval": 0.1817860094595569, + "fdr_pval": 0.2726790141893353 + } + ] + }, + "shared/False/0": { + "('A', 'B')": [ + { + "region": "p0", + "num_snps": 2, + "combined_mu": 0.6787471745078305, + "mu1": 0.7825601441901614, + "mu2": 0.2166705838515093, + "null_ll": -13.707277129349903, + "alt_ll": -9.397758537032335, + "pval": 0.0033266764056388686, + "fdr_pval": 0.009980029216916606 + }, + { + "region": "p1", + "num_snps": 1, + "combined_mu": 0.4999999999999996, + "mu1": 0.5467395941580712, + "mu2": 0.4532604058413437, + "null_ll": -4.418524248065584, + "alt_ll": -4.336386031479414, + "pval": 0.6852494418658779, + "fdr_pval": 1.0 + }, + { + "region": "p2", + "num_snps": 3, + "combined_mu": 0.6569575554611818, + "mu1": 0.6569575554611818, + "mu2": 0.6569575554611818, + "null_ll": -15.237030066709966, + "alt_ll": -15.237030066709966, + "pval": 1.0, + "fdr_pval": 1.0 + } + ] + }, + "different/False/0": { + "('A', 'B')": [ + { + "region": "p0", + "num_snps_group1": 2, + "num_snps_group2": 2, + "combined_mu": 0.6787471745078305, + "mu1": 0.7825601441901614, + "mu2": 0.2166705838515093, + "null_ll": -13.707277129349903, + "alt_ll": -9.397758537032335, + "pval": 0.0033266764056388686, + "fdr_pval": 0.009980029216916606 + }, + { + "region": "p1", + "num_snps_group1": 1, + "num_snps_group2": 1, + "combined_mu": 0.4999999999999996, + "mu1": 0.5467395941580712, + "mu2": 0.4532604058413437, + "null_ll": -4.418524248065584, + "alt_ll": -4.336386031479414, + "pval": 0.6852494418658779, + "fdr_pval": 1.0 + }, + { + "region": "p2", + "num_snps_group1": 3, + "num_snps_group2": 3, + "combined_mu": 0.6569575554611818, + "mu1": 0.6569575554611818, + "mu2": 0.6569575554611818, + "null_ll": -15.237030066709966, + "alt_ll": -15.237030066709966, + "pval": 1.0, + "fdr_pval": 1.0 + } + ] + }, + "single/False/1": { + "A": [ + { + "region": "p0", + "num_snps": 2, + "mu": 0.7576250662233635, + "null_ll": -7.243552812325376, + "alt_ll": -4.846331599764509, + "pval": 0.028551696581524982, + "fdr_pval": 0.08565508974457495 + }, + { + "region": "p1", + "num_snps": 1, + "mu": 0.5454545454545454, + "null_ll": -2.243298749477013, + "alt_ll": -2.2052443378160547, + "pval": 0.7826415596801843, + "fdr_pval": 0.7826415596801843 + }, + { + "region": "p2", + "num_snps": 3, + "mu": 0.6402051023191387, + "null_ll": -8.479890671235403, + "alt_ll": -7.706996371536085, + "pval": 0.21375823610460162, + "fdr_pval": 0.3206373541569024 + } + ], + "B": [ + { + "region": "p0", + "num_snps": 2, + "mu": 0.24125059985110295, + "null_ll": -7.243552812325376, + "alt_ll": -4.853517508661781, + "pval": 0.028790940177669772, + "fdr_pval": 0.08637282053300932 + }, + { + "region": "p1", + "num_snps": 1, + "mu": 0.45454545454545453, + "null_ll": -2.243298749477013, + "alt_ll": -2.2052443378160724, + "pval": 0.7826415596802337, + "fdr_pval": 0.7826415596802337 + }, + { + "region": "p2", + "num_snps": 3, + "mu": 0.6402051023191387, + "null_ll": -8.479890671235403, + "alt_ll": -7.706996371536085, + "pval": 0.21375823610460162, + "fdr_pval": 0.3206373541569024 + } + ] + }, + "shared/False/1": { + "('A', 'B')": [ + { + "region": "p0", + "num_snps": 2, + "combined_mu": 0.6568321069489483, + "mu1": 0.7576250662233635, + "mu2": 0.24125059985110295, + "null_ll": -13.597985992310939, + "alt_ll": -9.69984910842629, + "pval": 0.005235408951154251, + "fdr_pval": 0.015706226853462756 + }, + { + "region": "p1", + "num_snps": 1, + "combined_mu": 0.5000000000000001, + "mu1": 0.5428342230858644, + "mu2": 0.4571657769149307, + "null_ll": -4.486598462325851, + "alt_ll": -4.410200982537551, + "pval": 0.6958785766918056, + "fdr_pval": 1.0 + }, + { + "region": "p2", + "num_snps": 3, + "combined_mu": 0.6402051023191387, + "mu1": 0.6402051023191387, + "mu2": 0.6402051023191387, + "null_ll": -15.41399274307217, + "alt_ll": -15.41399274307217, + "pval": 1.0, + "fdr_pval": 1.0 + } + ] + }, + "different/False/1": { + "('A', 'B')": [ + { + "region": "p0", + "num_snps_group1": 2, + "num_snps_group2": 2, + "combined_mu": 0.6568321069489483, + "mu1": 0.7576250662233635, + "mu2": 0.24125059985110295, + "null_ll": -13.597985992310939, + "alt_ll": -9.69984910842629, + "pval": 0.005235408951154251, + "fdr_pval": 0.015706226853462756 + }, + { + "region": "p1", + "num_snps_group1": 1, + "num_snps_group2": 1, + "combined_mu": 0.5000000000000001, + "mu1": 0.5428342230858644, + "mu2": 0.4571657769149307, + "null_ll": -4.486598462325851, + "alt_ll": -4.410200982537551, + "pval": 0.6958785766918056, + "fdr_pval": 1.0 + }, + { + "region": "p2", + "num_snps_group1": 3, + "num_snps_group2": 3, + "combined_mu": 0.6402051023191387, + "mu1": 0.6402051023191387, + "mu2": 0.6402051023191387, + "null_ll": -15.41399274307217, + "alt_ll": -15.41399274307217, + "pval": 1.0, + "fdr_pval": 1.0 + } + ] + }, + "single/True/0": { + "A": [ + { + "region": "p0", + "num_snps": 2, + "mu": 0.7834514219321285, + "null_ll": -7.424952555448208, + "alt_ll": -4.0086271810134875, + "pval": 0.00895063353393758, + "fdr_pval": 0.026851900601812737 + }, + { + "region": "p1", + "num_snps": 1, + "mu": 0.55, + "null_ll": -2.2092614298154665, + "alt_ll": -2.1683943841147517, + "pval": 0.774961034462935, + "fdr_pval": 0.774961034462935 + }, + { + "region": "p2", + "num_snps": 3, + "mu": 0.6466635007634575, + "null_ll": -8.50999611779697, + "alt_ll": -7.289751714388885, + "pval": 0.11823937561912319, + "fdr_pval": 0.1773590634286848 + } + ], + "B": [ + { + "region": "p0", + "num_snps": 2, + "mu": 0.216548578068063, + "null_ll": -7.424952555448208, + "alt_ll": -4.008627181013487, + "pval": 0.00895063353393757, + "fdr_pval": 0.026851900601812712 + }, + { + "region": "p1", + "num_snps": 1, + "mu": 0.45, + "null_ll": -2.2092614298154665, + "alt_ll": -2.1683943841147535, + "pval": 0.7749610344629396, + "fdr_pval": 0.7749610344629396 + }, + { + "region": "p2", + "num_snps": 3, + "mu": 0.6752593726355599, + "null_ll": -8.50999611779697, + "alt_ll": -6.720902075037822, + "pval": 0.05854297334097624, + "fdr_pval": 0.08781446001146435 + } + ] + }, + "shared/True/0": { + "('A', 'B')": [ + { + "region": "p0", + "num_snps": 2, + "combined_mu": 0.49999999999999994, + "mu1": 0.7834514219321285, + "mu2": 0.216548578068063, + "null_ll": -14.849905110896408, + "alt_ll": -8.017254362026975, + "pval": 0.00021845386157032565, + "fdr_pval": 0.0006553615847109769 + }, + { + "region": "p1", + "num_snps": 1, + "combined_mu": 0.4999999999999996, + "mu1": 0.5467395941580712, + "mu2": 0.4532604058413437, + "null_ll": -4.418524248065584, + "alt_ll": -4.336386031479414, + "pval": 0.6852494418658779, + "fdr_pval": 0.8210187066840637 + }, + { + "region": "p2", + "num_snps": 3, + "combined_mu": 0.6611479228705941, + "mu1": 0.6466635007634575, + "mu2": 0.6752593726355599, + "null_ll": -14.036244910588938, + "alt_ll": -14.010653789426707, + "pval": 0.8210187066840637, + "fdr_pval": 0.8210187066840637 + } + ] + }, + "different/True/0": { + "('A', 'B')": [ + { + "region": "p0", + "num_snps_group1": 2, + "num_snps_group2": 2, + "combined_mu": 0.49999999999999994, + "mu1": 0.7834514219321285, + "mu2": 0.216548578068063, + "null_ll": -14.849905110896408, + "alt_ll": -8.017254362026975, + "pval": 0.00021845386157032565, + "fdr_pval": 0.0006553615847109769 + }, + { + "region": "p1", + "num_snps_group1": 1, + "num_snps_group2": 1, + "combined_mu": 0.4999999999999996, + "mu1": 0.5467395941580712, + "mu2": 0.4532604058413437, + "null_ll": -4.418524248065584, + "alt_ll": -4.336386031479414, + "pval": 0.6852494418658779, + "fdr_pval": 0.8210187066843629 + }, + { + "region": "p2", + "num_snps_group1": 3, + "num_snps_group2": 3, + "combined_mu": 0.3388520771352042, + "mu1": 0.35333649923869614, + "mu2": 0.3247406273634475, + "null_ll": -14.036244910588955, + "alt_ll": -14.010653789426811, + "pval": 0.8210187066843629, + "fdr_pval": 0.8210187066843629 + } + ] + }, + "single/True/1": { + "A": [ + { + "region": "p0", + "num_snps": 2, + "mu": 0.758982457014652, + "null_ll": -7.243552812325376, + "alt_ll": -4.161349689458703, + "pval": 0.013034611692363305, + "fdr_pval": 0.039103835077089914 + }, + { + "region": "p1", + "num_snps": 1, + "mu": 0.5454545454545454, + "null_ll": -2.243298749477013, + "alt_ll": -2.2052443378160547, + "pval": 0.7826415596801843, + "fdr_pval": 0.7826415596801843 + }, + { + "region": "p2", + "num_snps": 3, + "mu": 0.6332105568137278, + "null_ll": -8.479890671235403, + "alt_ll": -7.366466615233957, + "pval": 0.13563027815445183, + "fdr_pval": 0.20344541723167775 + } + ], + "B": [ + { + "region": "p0", + "num_snps": 2, + "mu": 0.2410175429961819, + "null_ll": -7.243552812325376, + "alt_ll": -4.161349689458719, + "pval": 0.013034611692363534, + "fdr_pval": 0.0391038350770906 + }, + { + "region": "p1", + "num_snps": 1, + "mu": 0.45454545454545453, + "null_ll": -2.243298749477013, + "alt_ll": -2.2052443378160724, + "pval": 0.7826415596802337, + "fdr_pval": 0.7826415596802337 + }, + { + "region": "p2", + "num_snps": 3, + "mu": 0.659871232807733, + "null_ll": -8.479890671235403, + "alt_ll": -6.841264146507394, + "pval": 0.07024641070689956, + "fdr_pval": 0.10536961606034934 + } + ] + }, + "shared/True/1": { + "('A', 'B')": [ + { + "region": "p0", + "num_snps": 2, + "combined_mu": 0.5000000000000002, + "mu1": 0.758982457014652, + "mu2": 0.2410175429961819, + "null_ll": -14.487105624650727, + "alt_ll": -8.322699378917422, + "pval": 0.0004460200551571753, + "fdr_pval": 0.0013380601654715259 + }, + { + "region": "p1", + "num_snps": 1, + "combined_mu": 0.5000000000000001, + "mu1": 0.5428342230858644, + "mu2": 0.4571657769149307, + "null_ll": -4.486598462325851, + "alt_ll": -4.410200982537551, + "pval": 0.6958785766918056, + "fdr_pval": 0.8260474537950093 + }, + { + "region": "p2", + "num_snps": 3, + "combined_mu": 0.646682000767002, + "mu1": 0.6332105568137278, + "mu2": 0.659871232807733, + "null_ll": -14.231880987079016, + "alt_ll": -14.20773076174135, + "pval": 0.8260474537950093, + "fdr_pval": 0.8260474537950093 + } + ] + }, + "different/True/1": { + "('A', 'B')": [ + { + "region": "p0", + "num_snps_group1": 2, + "num_snps_group2": 2, + "combined_mu": 0.5000000000000002, + "mu1": 0.758982457014652, + "mu2": 0.2410175429961819, + "null_ll": -14.487105624650727, + "alt_ll": -8.322699378917422, + "pval": 0.0004460200551571753, + "fdr_pval": 0.0013380601654715259 + }, + { + "region": "p1", + "num_snps_group1": 1, + "num_snps_group2": 1, + "combined_mu": 0.5000000000000001, + "mu1": 0.5428342230858644, + "mu2": 0.4571657769149307, + "null_ll": -4.486598462325851, + "alt_ll": -4.410200982537551, + "pval": 0.6958785766918056, + "fdr_pval": 0.8260474537934291 + }, + { + "region": "p2", + "num_snps_group1": 3, + "num_snps_group2": 3, + "combined_mu": 0.3533179992352407, + "mu1": 0.36678943923436186, + "mu2": 0.340128767195461, + "null_ll": -14.231880987079013, + "alt_ll": -14.207730761740901, + "pval": 0.8260474537934291, + "fdr_pval": 0.8260474537934291 + } + ] + } + }, + "linear_reference_adaptations": [ + "Scalarize one-element objectives", + "Align per-SNP dispersion arrays; preserve archived first-SNP likelihood and chi-square" + ], + "linear": { + "snv/False/0": [ + { + "region": "chr1_100", + "ref_count": 18, + "alt_count": 2, + "N": 20, + "snp_count": 1, + "null_ll": -3.4961229205272923, + "alt_ll": -2.109317536064034, + "mu": 0.8275747474065025, + "lrt": 2.773610768926517, + "pval": 0.09582977056890955, + "fdr_pval": 0.580267974343758 + }, + { + "region": "chr1_200", + "ref_count": 3, + "alt_count": 17, + "N": 20, + "snp_count": 1, + "null_ll": -3.234950705529754, + "alt_ll": -2.29533327249182, + "mu": 0.2170612477601647, + "lrt": 1.8792348660758682, + "pval": 0.17042113568024736, + "fdr_pval": 0.580267974343758 + }, + { + "region": "chr1_300", + "ref_count": 22, + "alt_count": 18, + "N": 40, + "snp_count": 1, + "null_ll": -3.396439248411714, + "alt_ll": -3.383683487896711, + "mu": 0.5377508602789496, + "lrt": 0.025511521030005646, + "pval": 0.873099074231311, + "fdr_pval": 0.9524717173432482 + }, + { + "region": "chr1_400", + "ref_count": 25, + "alt_count": 15, + "N": 40, + "snp_count": 1, + "null_ll": -3.435773013757093, + "alt_ll": -3.3539133805262407, + "mu": 0.5945585857280776, + "lrt": 0.16371926646170465, + "pval": 0.6857551089025606, + "fdr_pval": 0.8766593045255179 + }, + { + "region": "chr1_500", + "ref_count": 29, + "alt_count": 11, + "N": 40, + "snp_count": 1, + "null_ll": -3.55136480569324, + "alt_ll": -3.264962640283396, + "mu": 0.6712351175796699, + "lrt": 0.5728043308196886, + "pval": 0.44914669701506194, + "fdr_pval": 0.8198962466650255 + }, + { + "region": "chr1_600", + "ref_count": 10, + "alt_count": 30, + "N": 40, + "snp_count": 1, + "null_ll": -3.594763043107044, + "alt_ll": -3.2309555294775283, + "mu": 0.30928566206817376, + "lrt": 0.7276150272590316, + "pval": 0.3936573449933354, + "fdr_pval": 0.8198962466650255 + }, + { + "region": "chr1_700", + "ref_count": 10, + "alt_count": 10, + "N": 20, + "snp_count": 1, + "null_ll": -2.6510455760508753, + "alt_ll": -2.651045576050875, + "mu": 0.4999999999999988, + "lrt": 8.881784197001252e-16, + "pval": 0.9999999762211871, + "fdr_pval": 0.9999999762211871 + }, + { + "region": "chr1_800", + "ref_count": 8, + "alt_count": 12, + "N": 20, + "snp_count": 1, + "null_ll": -2.6889798571512484, + "alt_ll": -2.629675977852769, + "mu": 0.4207515946822724, + "lrt": 0.11860775859695849, + "pval": 0.7305494204379316, + "fdr_pval": 0.8766593045255179 + }, + { + "region": "chr1_900", + "ref_count": 35, + "alt_count": 5, + "N": 40, + "snp_count": 1, + "null_ll": -3.961516793056401, + "alt_ll": -2.9268713261072543, + "mu": 0.7926379229113844, + "lrt": 2.0692909338982934, + "pval": 0.15029138166082384, + "fdr_pval": 0.580267974343758 + }, + { + "region": "chr1_1000", + "ref_count": 6, + "alt_count": 34, + "N": 40, + "snp_count": 1, + "null_ll": -3.8600674080728443, + "alt_ll": -3.014390711076592, + "mu": 0.2287050839496923, + "lrt": 1.6913533939925047, + "pval": 0.193422658114586, + "fdr_pval": 0.580267974343758 + }, + { + "region": "chr1_1100", + "ref_count": 14, + "alt_count": 6, + "N": 20, + "snp_count": 1, + "null_ll": -2.810988420177069, + "alt_ll": -2.559588885453694, + "mu": 0.6590894791618173, + "lrt": 0.5027990694467501, + "pval": 0.47827281055459825, + "fdr_pval": 0.8198962466650255 + }, + { + "region": "chr1_1200", + "ref_count": 7, + "alt_count": 13, + "N": 20, + "snp_count": 1, + "null_ll": -2.738242057227434, + "alt_ll": -2.6016352880739664, + "mu": 0.3809538526189891, + "lrt": 0.2732135383069352, + "pval": 0.6011842044980265, + "fdr_pval": 0.8766593045255179 + } + ], + "snv/False/1": [ + { + "region": "chr1_100", + "ref_count": 18, + "alt_count": 2, + "N": 20, + "snp_count": 1, + "null_ll": -3.594636875735214, + "alt_ll": -2.2247235605535174, + "mu": 0.8087630535565404, + "lrt": 2.7398266303633934, + "pval": 0.09787538499393948, + "fdr_pval": 0.5729834993939965 + }, + { + "region": "chr1_200", + "ref_count": 3, + "alt_count": 17, + "N": 20, + "snp_count": 1, + "null_ll": -3.3121770622913194, + "alt_ll": -2.3460171799505356, + "mu": 0.23125809723699667, + "lrt": 1.9323197646815675, + "pval": 0.16450487428652077, + "fdr_pval": 0.5729834993939965 + }, + { + "region": "chr1_300", + "ref_count": 22, + "alt_count": 18, + "N": 40, + "snp_count": 1, + "null_ll": -3.379334764082913, + "alt_ll": -3.365902024229376, + "mu": 0.5371420574032392, + "lrt": 0.02686547970707398, + "pval": 0.869804416325196, + "fdr_pval": 0.9488775450820319 + }, + { + "region": "chr1_400", + "ref_count": 25, + "alt_count": 15, + "N": 40, + "snp_count": 1, + "null_ll": -3.4253021326233513, + "alt_ll": -3.3393471707302593, + "mu": 0.5929817351834095, + "lrt": 0.17190992378618386, + "pval": 0.6784199463914897, + "fdr_pval": 0.8617687681803181 + }, + { + "region": "chr1_500", + "ref_count": 29, + "alt_count": 11, + "N": 40, + "snp_count": 1, + "null_ll": -3.558981844278041, + "alt_ll": -3.260962606869964, + "mu": 0.6680833775015869, + "lrt": 0.5960384748161536, + "pval": 0.4400935318227224, + "fdr_pval": 0.7891104746897067 + }, + { + "region": "chr1_600", + "ref_count": 10, + "alt_count": 30, + "N": 40, + "snp_count": 1, + "null_ll": -3.6085931085126366, + "alt_ll": -3.231398683524503, + "mu": 0.31292874163871576, + "lrt": 0.7543888499762668, + "pval": 0.3850902424105198, + "fdr_pval": 0.7891104746897067 + }, + { + "region": "chr1_700", + "ref_count": 10, + "alt_count": 10, + "N": 20, + "snp_count": 1, + "null_ll": -2.613772095346134, + "alt_ll": -2.613772095346133, + "mu": 0.5000000000000013, + "lrt": 1.7763568394002505e-15, + "pval": 0.9999999663716803, + "fdr_pval": 0.9999999663716803 + }, + { + "region": "chr1_800", + "ref_count": 8, + "alt_count": 12, + "N": 20, + "snp_count": 1, + "null_ll": -2.6616558184050967, + "alt_ll": -2.596515209229925, + "mu": 0.4239080527959045, + "lrt": 0.13028121835034323, + "pval": 0.7181406401502651, + "fdr_pval": 0.8617687681803181 + }, + { + "region": "chr1_900", + "ref_count": 35, + "alt_count": 5, + "N": 40, + "snp_count": 1, + "null_ll": -4.012995914768656, + "alt_ll": -2.9786613774158957, + "mu": 0.7849799337518385, + "lrt": 2.06866907470552, + "pval": 0.15035268027878387, + "fdr_pval": 0.5729834993939965 + }, + { + "region": "chr1_1000", + "ref_count": 6, + "alt_count": 34, + "N": 40, + "snp_count": 1, + "null_ll": -3.9040396471475494, + "alt_ll": -3.049074827895322, + "mu": 0.23521267826123604, + "lrt": 1.7099296385044545, + "pval": 0.19099449979799885, + "fdr_pval": 0.5729834993939965 + }, + { + "region": "chr1_1100", + "ref_count": 14, + "alt_count": 6, + "N": 20, + "snp_count": 1, + "null_ll": -2.813451792008641, + "alt_ll": -2.540886808994198, + "mu": 0.6524545528674803, + "lrt": 0.5451299660288864, + "pval": 0.46031444356899565, + "fdr_pval": 0.7891104746897067 + }, + { + "region": "chr1_1200", + "ref_count": 7, + "alt_count": 13, + "N": 20, + "snp_count": 1, + "null_ll": -2.7233572320579715, + "alt_ll": -2.5740778587163065, + "mu": 0.3857828701149865, + "lrt": 0.29855874668333016, + "pval": 0.5847873709891638, + "fdr_pval": 0.8617687681803181 + } + ], + "snv/True/0": [ + { + "region": "chr1_100", + "ref_count": 18, + "alt_count": 2, + "N": 20, + "snp_count": 1, + "null_ll": -3.4961229205272923, + "alt_ll": -2.109317536064034, + "mu": 0.8275747474065025, + "lrt": 2.773610768926517, + "pval": 0.09582977056890955, + "fdr_pval": 0.580267974343758 + }, + { + "region": "chr1_200", + "ref_count": 3, + "alt_count": 17, + "N": 20, + "snp_count": 1, + "null_ll": -3.234950705529754, + "alt_ll": -2.29533327249182, + "mu": 0.2170612477601647, + "lrt": 1.8792348660758682, + "pval": 0.17042113568024736, + "fdr_pval": 0.580267974343758 + }, + { + "region": "chr1_300", + "ref_count": 22, + "alt_count": 18, + "N": 40, + "snp_count": 1, + "null_ll": -3.396439248411714, + "alt_ll": -3.383683487896711, + "mu": 0.5377508602789496, + "lrt": 0.025511521030005646, + "pval": 0.873099074231311, + "fdr_pval": 0.9524717173432482 + }, + { + "region": "chr1_400", + "ref_count": 25, + "alt_count": 15, + "N": 40, + "snp_count": 1, + "null_ll": -3.435773013757093, + "alt_ll": -3.3539133805262407, + "mu": 0.5945585857280776, + "lrt": 0.16371926646170465, + "pval": 0.6857551089025606, + "fdr_pval": 0.8766593045255179 + }, + { + "region": "chr1_500", + "ref_count": 29, + "alt_count": 11, + "N": 40, + "snp_count": 1, + "null_ll": -3.55136480569324, + "alt_ll": -3.264962640283396, + "mu": 0.6712351175796699, + "lrt": 0.5728043308196886, + "pval": 0.44914669701506194, + "fdr_pval": 0.8198962466650255 + }, + { + "region": "chr1_600", + "ref_count": 10, + "alt_count": 30, + "N": 40, + "snp_count": 1, + "null_ll": -3.594763043107044, + "alt_ll": -3.2309555294775283, + "mu": 0.30928566206817376, + "lrt": 0.7276150272590316, + "pval": 0.3936573449933354, + "fdr_pval": 0.8198962466650255 + }, + { + "region": "chr1_700", + "ref_count": 10, + "alt_count": 10, + "N": 20, + "snp_count": 1, + "null_ll": -2.6510455760508753, + "alt_ll": -2.651045576050875, + "mu": 0.4999999999999988, + "lrt": 8.881784197001252e-16, + "pval": 0.9999999762211871, + "fdr_pval": 0.9999999762211871 + }, + { + "region": "chr1_800", + "ref_count": 8, + "alt_count": 12, + "N": 20, + "snp_count": 1, + "null_ll": -2.6889798571512484, + "alt_ll": -2.629675977852769, + "mu": 0.4207515946822724, + "lrt": 0.11860775859695849, + "pval": 0.7305494204379316, + "fdr_pval": 0.8766593045255179 + }, + { + "region": "chr1_900", + "ref_count": 35, + "alt_count": 5, + "N": 40, + "snp_count": 1, + "null_ll": -3.961516793056401, + "alt_ll": -2.9268713261072543, + "mu": 0.7926379229113844, + "lrt": 2.0692909338982934, + "pval": 0.15029138166082384, + "fdr_pval": 0.580267974343758 + }, + { + "region": "chr1_1000", + "ref_count": 6, + "alt_count": 34, + "N": 40, + "snp_count": 1, + "null_ll": -3.8600674080728443, + "alt_ll": -3.014390711076592, + "mu": 0.2287050839496923, + "lrt": 1.6913533939925047, + "pval": 0.193422658114586, + "fdr_pval": 0.580267974343758 + }, + { + "region": "chr1_1100", + "ref_count": 14, + "alt_count": 6, + "N": 20, + "snp_count": 1, + "null_ll": -2.810988420177069, + "alt_ll": -2.559588885453694, + "mu": 0.6590894791618173, + "lrt": 0.5027990694467501, + "pval": 0.47827281055459825, + "fdr_pval": 0.8198962466650255 + }, + { + "region": "chr1_1200", + "ref_count": 7, + "alt_count": 13, + "N": 20, + "snp_count": 1, + "null_ll": -2.738242057227434, + "alt_ll": -2.6016352880739664, + "mu": 0.3809538526189891, + "lrt": 0.2732135383069352, + "pval": 0.6011842044980265, + "fdr_pval": 0.8766593045255179 + } + ], + "snv/True/1": [ + { + "region": "chr1_100", + "ref_count": 18, + "alt_count": 2, + "N": 20, + "snp_count": 1, + "null_ll": -3.594636875735214, + "alt_ll": -2.2247235605535174, + "mu": 0.8087630535565404, + "lrt": 2.7398266303633934, + "pval": 0.09787538499393948, + "fdr_pval": 0.5729834993939965 + }, + { + "region": "chr1_200", + "ref_count": 3, + "alt_count": 17, + "N": 20, + "snp_count": 1, + "null_ll": -3.3121770622913194, + "alt_ll": -2.3460171799505356, + "mu": 0.23125809723699667, + "lrt": 1.9323197646815675, + "pval": 0.16450487428652077, + "fdr_pval": 0.5729834993939965 + }, + { + "region": "chr1_300", + "ref_count": 22, + "alt_count": 18, + "N": 40, + "snp_count": 1, + "null_ll": -3.379334764082913, + "alt_ll": -3.365902024229376, + "mu": 0.5371420574032392, + "lrt": 0.02686547970707398, + "pval": 0.869804416325196, + "fdr_pval": 0.9488775450820319 + }, + { + "region": "chr1_400", + "ref_count": 25, + "alt_count": 15, + "N": 40, + "snp_count": 1, + "null_ll": -3.4253021326233513, + "alt_ll": -3.3393471707302593, + "mu": 0.5929817351834095, + "lrt": 0.17190992378618386, + "pval": 0.6784199463914897, + "fdr_pval": 0.8617687681803181 + }, + { + "region": "chr1_500", + "ref_count": 29, + "alt_count": 11, + "N": 40, + "snp_count": 1, + "null_ll": -3.558981844278041, + "alt_ll": -3.260962606869964, + "mu": 0.6680833775015869, + "lrt": 0.5960384748161536, + "pval": 0.4400935318227224, + "fdr_pval": 0.7891104746897067 + }, + { + "region": "chr1_600", + "ref_count": 10, + "alt_count": 30, + "N": 40, + "snp_count": 1, + "null_ll": -3.6085931085126366, + "alt_ll": -3.231398683524503, + "mu": 0.31292874163871576, + "lrt": 0.7543888499762668, + "pval": 0.3850902424105198, + "fdr_pval": 0.7891104746897067 + }, + { + "region": "chr1_700", + "ref_count": 10, + "alt_count": 10, + "N": 20, + "snp_count": 1, + "null_ll": -2.613772095346134, + "alt_ll": -2.613772095346133, + "mu": 0.5000000000000013, + "lrt": 1.7763568394002505e-15, + "pval": 0.9999999663716803, + "fdr_pval": 0.9999999663716803 + }, + { + "region": "chr1_800", + "ref_count": 8, + "alt_count": 12, + "N": 20, + "snp_count": 1, + "null_ll": -2.6616558184050967, + "alt_ll": -2.596515209229925, + "mu": 0.4239080527959045, + "lrt": 0.13028121835034323, + "pval": 0.7181406401502651, + "fdr_pval": 0.8617687681803181 + }, + { + "region": "chr1_900", + "ref_count": 35, + "alt_count": 5, + "N": 40, + "snp_count": 1, + "null_ll": -4.012995914768656, + "alt_ll": -2.9786613774158957, + "mu": 0.7849799337518385, + "lrt": 2.06866907470552, + "pval": 0.15035268027878387, + "fdr_pval": 0.5729834993939965 + }, + { + "region": "chr1_1000", + "ref_count": 6, + "alt_count": 34, + "N": 40, + "snp_count": 1, + "null_ll": -3.9040396471475494, + "alt_ll": -3.049074827895322, + "mu": 0.23521267826123604, + "lrt": 1.7099296385044545, + "pval": 0.19099449979799885, + "fdr_pval": 0.5729834993939965 + }, + { + "region": "chr1_1100", + "ref_count": 14, + "alt_count": 6, + "N": 20, + "snp_count": 1, + "null_ll": -2.813451792008641, + "alt_ll": -2.540886808994198, + "mu": 0.6524545528674803, + "lrt": 0.5451299660288864, + "pval": 0.46031444356899565, + "fdr_pval": 0.7891104746897067 + }, + { + "region": "chr1_1200", + "ref_count": 7, + "alt_count": 13, + "N": 20, + "snp_count": 1, + "null_ll": -2.7233572320579715, + "alt_ll": -2.5740778587163065, + "mu": 0.3857828701149865, + "lrt": 0.29855874668333016, + "pval": 0.5847873709891638, + "fdr_pval": 0.8617687681803181 + } + ], + "atac/False/0": [ + { + "region": "peak0", + "ref_count": 21, + "alt_count": 19, + "N": 40, + "snp_count": 2, + "null_ll": -6.7310736260570465, + "alt_ll": -5.091513640623601, + "mu": 0.8023617349407871, + "lrt": 3.2791199708668906, + "pval": 0.07016654319656014, + "fdr_pval": 0.352477499453001 + }, + { + "region": "peak1", + "ref_count": 47, + "alt_count": 33, + "N": 80, + "snp_count": 2, + "null_ll": -6.832212262168807, + "alt_ll": -6.825241422784275, + "mu": 0.5206521913613598, + "lrt": 0.01394167876906316, + "pval": 0.9060083220812054, + "fdr_pval": 1.0 + }, + { + "region": "peak2", + "ref_count": 39, + "alt_count": 41, + "N": 80, + "snp_count": 2, + "null_ll": -7.146127848800284, + "alt_ll": -6.919578893370216, + "mu": 0.6334729072333185, + "lrt": 0.4530979108601372, + "pval": 0.5008674654547305, + "fdr_pval": 0.8801202029620296 + }, + { + "region": "peak3", + "ref_count": 18, + "alt_count": 22, + "N": 40, + "snp_count": 2, + "null_ll": -5.340025433202124, + "alt_ll": -5.340025433202125, + "mu": 0.4999999999999987, + "lrt": -1.7763568394002505e-15, + "pval": 1.0, + "fdr_pval": 1.0 + }, + { + "region": "peak4", + "ref_count": 41, + "alt_count": 39, + "N": 80, + "snp_count": 2, + "null_ll": -7.821584201129245, + "alt_ll": -6.59636800922164, + "mu": 0.774875974762779, + "lrt": 2.45043238381521, + "pval": 0.117492499817667, + "fdr_pval": 0.352477499453001 + }, + { + "region": "peak5", + "ref_count": 21, + "alt_count": 19, + "N": 40, + "snp_count": 2, + "null_ll": -5.549230477404503, + "alt_ll": -5.401503722915232, + "mu": 0.595083388542957, + "lrt": 0.29545350897854306, + "pval": 0.5867468019746864, + "fdr_pval": 0.8801202029620296 + } + ], + "atac/False/1": [ + { + "region": "peak0", + "ref_count": 21, + "alt_count": 19, + "N": 40, + "snp_count": 2, + "null_ll": -6.906813938026533, + "alt_ll": -5.257960673034265, + "mu": 0.7857735703258241, + "lrt": 3.297706529984536, + "pval": 0.06937668532440518, + "fdr_pval": 0.34937904715898405 + }, + { + "region": "peak1", + "ref_count": 47, + "alt_count": 33, + "N": 80, + "snp_count": 2, + "null_ll": -6.804636896706264, + "alt_ll": -6.797264153569046, + "mu": 0.5204058828114936, + "lrt": 0.01474548627443717, + "pval": 0.9033496689512466, + "fdr_pval": 0.9999999327433605 + }, + { + "region": "peak2", + "ref_count": 39, + "alt_count": 41, + "N": 80, + "snp_count": 2, + "null_ll": -7.167574952790678, + "alt_ll": -6.928274455404915, + "mu": 0.6324488278932066, + "lrt": 0.47860099477152573, + "pval": 0.4890566927689526, + "fdr_pval": 0.8538270915863685 + }, + { + "region": "peak3", + "ref_count": 18, + "alt_count": 22, + "N": 40, + "snp_count": 2, + "null_ll": -5.275427913751231, + "alt_ll": -5.275427913751227, + "mu": 0.5000000000000011, + "lrt": 7.105427357601002e-15, + "pval": 0.9999999327433605, + "fdr_pval": 0.9999999327433605 + }, + { + "region": "peak4", + "ref_count": 41, + "alt_count": 39, + "N": 80, + "snp_count": 2, + "null_ll": -7.9170355619162045, + "alt_ll": -6.68488628548614, + "mu": 0.7680124049311785, + "lrt": 2.4642985528601287, + "pval": 0.116459682386328, + "fdr_pval": 0.34937904715898405 + }, + { + "region": "peak5", + "ref_count": 21, + "alt_count": 19, + "N": 40, + "snp_count": 2, + "null_ll": -5.536809024066613, + "alt_ll": -5.3748127247284865, + "mu": 0.5918554050710805, + "lrt": 0.3239925986762522, + "pval": 0.569218061057579, + "fdr_pval": 0.8538270915863685 + } + ], + "rna/False/0": [ + { + "region": "gene0", + "ref_count": 43, + "alt_count": 37, + "N": 80, + "snp_count": 3, + "null_ll": -10.12751287446876, + "alt_ll": -9.131296791027125, + "mu": 0.7198438329549017, + "lrt": 1.992432166883269, + "pval": 0.15808680710220702, + "fdr_pval": 0.6323472284088281 + }, + { + "region": "gene1", + "ref_count": 64, + "alt_count": 56, + "N": 120, + "snp_count": 3, + "null_ll": -10.581900862557378, + "alt_ll": -10.531979947390589, + "mu": 0.5575724445131891, + "lrt": 0.09984183033357752, + "pval": 0.752019526430494, + "fdr_pval": 1.0 + }, + { + "region": "gene2", + "ref_count": 53, + "alt_count": 27, + "N": 80, + "snp_count": 3, + "null_ll": -9.301542226258524, + "alt_ll": -9.301542226258535, + "mu": 0.4999999999999985, + "lrt": -2.1316282072803006e-14, + "pval": 1.0, + "fdr_pval": 1.0 + }, + { + "region": "gene3", + "ref_count": 27, + "alt_count": 53, + "N": 80, + "snp_count": 3, + "null_ll": -9.409297885477347, + "alt_ll": -9.015436446390648, + "mu": 0.3653537784030239, + "lrt": 0.7877228781733976, + "pval": 0.374789566454941, + "fdr_pval": 0.749579132909882 + } + ], + "rna/False/1": [ + { + "region": "gene0", + "ref_count": 43, + "alt_count": 37, + "N": 80, + "snp_count": 3, + "null_ll": -10.286148702109447, + "alt_ll": -9.256236250752744, + "mu": 0.7094500879382742, + "lrt": 2.0598249027134052, + "pval": 0.15122754685170467, + "fdr_pval": 0.6049101874068187 + }, + { + "region": "gene1", + "ref_count": 64, + "alt_count": 56, + "N": 120, + "snp_count": 3, + "null_ll": -10.59287708541403, + "alt_ll": -10.539030769648871, + "mu": 0.558048096412444, + "lrt": 0.10769263153031616, + "pval": 0.7427864900537504, + "fdr_pval": 0.9903819867383339 + }, + { + "region": "gene2", + "ref_count": 53, + "alt_count": 27, + "N": 80, + "snp_count": 3, + "null_ll": -9.288423828519885, + "alt_ll": -9.288423828519882, + "mu": 0.5000000000000016, + "lrt": 7.105427357601002e-15, + "pval": 0.9999999327433605, + "fdr_pval": 0.9999999327433605 + }, + { + "region": "gene3", + "ref_count": 27, + "alt_count": 53, + "N": 80, + "snp_count": 3, + "null_ll": -9.440848671214162, + "alt_ll": -9.054190555429845, + "mu": 0.3733528758302787, + "lrt": 0.7733162315686357, + "pval": 0.3791931006762539, + "fdr_pval": 0.7583862013525078 + } + ] + } +} diff --git a/tests/analysis/test_archive_parity.py b/tests/analysis/test_archive_parity.py new file mode 100644 index 0000000..0c5bcd5 --- /dev/null +++ b/tests/analysis/test_archive_parity.py @@ -0,0 +1,173 @@ +"""Production-code parity with outputs frozen from unmodified Python 41354fc. + +The fixture is independent of Rust and does not need archived git objects in CI. +These tests cover statistical compatibility on synthetic counts, not calibration. +""" + +import json +from pathlib import Path + +import numpy as np +import pandas as pd +import pytest +from scripts.compare_archive_analysis import bulk_counts, compare_frames, single_cell_counts +from typer.testing import CliRunner + +from analysis.as_analysis_sc import get_imbalance_sc +from analysis.compare_ai import get_compared_imbalance, get_compared_imbalance_diff_snps +from analysis.run_analysis import run_ai_analysis +from analysis.run_analysis_sc import process_adata_inputs + +REFERENCE = json.loads( + (Path(__file__).parent / "fixtures/archive_python_reference.json").read_text() +) + + +@pytest.mark.rust +@pytest.mark.parametrize("unit", ["snv", "atac", "rna"]) +@pytest.mark.parametrize("phased", [False, True]) +@pytest.mark.parametrize("pc", [0, 1]) +def test_bulk_archive_parity(tmp_path, unit, phased, pc): + region = {"snv": None, "atac": "region", "rna": "gene"}[unit] + counts = bulk_counts().drop( + columns=["region", "gene"] if unit == "snv" else ["gene"] if unit == "atac" else ["region"] + ) + count_file, output = tmp_path / "counts.tsv", tmp_path / "results.tsv" + counts.to_csv(count_file, sep="\t", index=False) + run_ai_analysis( + count_file, + model="single", + phased=phased, + region_col=region, + pseudocount=pc, + out_file=str(output), + ) + expected = pd.DataFrame(REFERENCE["bulk"][f"{unit}/{phased}/{pc}"]) + actual = pd.read_csv(output, sep="\t") + comparison = compare_frames(expected, actual) + assert comparison["status"] == "PASS", comparison + + +@pytest.mark.rust +@pytest.mark.parametrize( + "unit,phased", [("snv", False), ("snv", True), ("atac", False), ("rna", False)] +) +@pytest.mark.parametrize("pc", [0, 1]) +def test_linear_shape_adapted_reference(tmp_path, unit, phased, pc): + # The fixture records the reference-only adaptations. Phased multi-SNP + # linear results are excluded because the archive silently ignored phasing. + region = {"snv": None, "atac": "region", "rna": "gene"}[unit] + counts = bulk_counts().drop( + columns=["region", "gene"] if unit == "snv" else ["gene"] if unit == "atac" else ["region"] + ) + count_file, output = tmp_path / "counts.tsv", tmp_path / "results.tsv" + counts.to_csv(count_file, sep="\t", index=False) + run_ai_analysis( + count_file, + model="linear", + phased=phased, + region_col=region, + pseudocount=pc, + out_file=str(output), + ) + expected = pd.DataFrame(REFERENCE["linear"][f"{unit}/{phased}/{pc}"]) + comparison = compare_frames(expected, pd.read_csv(output, sep="\t")) + assert comparison["status"] == "PASS", comparison + + +@pytest.mark.parametrize("mode", ["single", "shared", "different"]) +@pytest.mark.parametrize("phased", [False, True]) +@pytest.mark.parametrize("pc", [0, 1]) +def test_single_cell_archive_parity(mode, phased, pc): + current = single_cell_counts().T.copy() + prepared = process_adata_inputs(current, sample="sample1", phased=phased) + function = { + "single": get_imbalance_sc, + "shared": get_compared_imbalance, + "different": get_compared_imbalance_diff_snps, + }[mode] + actual = function( + prepared.adata, + sample=prepared.sample, + phased=prepared.phased, + groups=prepared.groups, + pseudocount=pc, + min_count=10, + ) + expected = REFERENCE["single_cell"][f"{mode}/{phased}/{pc}"] + assert set(expected) == {str(k) for k in actual} + for key, frame in actual.items(): + comparison = compare_frames(pd.DataFrame(expected[str(key)]), frame) + assert comparison["status"] == "PASS", comparison + + +def test_barcode_mapping_and_genotype_filter_preserve_feature_indices(tmp_path): + current = single_cell_counts().T.copy() + current.var.loc[current.var.index[1], "sample1"] = "0/0" + mapping = tmp_path / "groups.tsv" + mapping.write_text("cell0\tnewA\ncell1\tnewA\ncell2\tnewB\ncell3\tnewB\n") + prepared = process_adata_inputs(current, bc_map=mapping, sample="sample1", phased=False) + assert prepared.adata.shape == (4, 5) + assert prepared.adata.obs["group"].tolist() == ["newA", "newA", "newB", "newB"] + assert prepared.adata.var["index"].tolist() == list(range(5)) + assert prepared.adata.uns["feature"]["index"].tolist() == list(range(5)) + assert prepared.adata.uns["feature"]["region"].tolist() == ["p0", "p1", "p2", "p2", "p2"] + + +def test_unphased_counts_without_genotype_metadata(): + current = single_cell_counts().T.copy() + del current.uns["samples"] + current.var = current.var.drop(columns="sample1") + prepared = process_adata_inputs(current) + assert prepared.sample is None + assert prepared.phased is False + assert set(get_imbalance_sc(prepared.adata, groups=prepared.groups)) == {"A", "B"} + + +@pytest.mark.parametrize("mode", ["shared", "different"]) +def test_differential_sums_do_not_wrap_at_uint16(monkeypatch, mode): + import analysis.compare_ai as comparison + + current = single_cell_counts().T.copy() + current.layers["ref"] = current.layers["ref"].astype(np.uint32) * 10000 + current.layers["alt"] = current.layers["alt"].astype(np.uint32) * 10000 + captured = [] + if mode == "shared": + + def capture(disp, ref1, n1, ref2, n2, regions, gt): + captured.append((ref1.copy(), n1.copy(), ref2.copy(), n2.copy())) + return pd.DataFrame({"region": ["p0"]}) + + monkeypatch.setattr(comparison, "compare_imbalance_between_groups", capture) + function = comparison.get_compared_imbalance + else: + + def capture(disp, ref1, n1, gt1, regions1, ref2, n2, gt2, regions2): + captured.append((ref1.copy(), n1.copy(), ref2.copy(), n2.copy())) + return pd.DataFrame({"region": ["p0"]}) + + monkeypatch.setattr(comparison, "compare_imbalance_between_groups_diff_snps", capture) + function = comparison.get_compared_imbalance_diff_snps + function(current, sample="sample1", phased=False, pseudocount=0) + assert len(captured) == 1 + ref1, n1, ref2, n2 = captured[0] + np.testing.assert_array_equal(ref1, np.asarray(current.layers["ref"][:2].sum(axis=0)).ravel()) + np.testing.assert_array_equal(ref2, np.asarray(current.layers["ref"][2:].sum(axis=0)).ravel()) + np.testing.assert_array_equal(n1, np.full(6, 200000)) + np.testing.assert_array_equal(n2, np.full(6, 200000)) + + +@pytest.mark.parametrize("command", ["find-imbalance-sc", "compare-imbalance"]) +def test_single_cell_cli_reads_current_count_layout(tmp_path, command): + from analysis.__main__ import app + + counts, groups, output = tmp_path / "counts.h5ad", tmp_path / "groups.tsv", tmp_path / "ai.tsv" + single_cell_counts().T.copy().write_h5ad(counts) + groups.write_text("cell0\tA\ncell1\tA\ncell2\tB\ncell3\tB\n") + result = CliRunner().invoke( + app, [command, str(counts), str(groups), "--sample", "sample1", "--out_file", str(output)] + ) + assert result.exit_code == 0, result.output + repr(result.exception) + files = sorted(tmp_path.glob("ai_*.tsv")) + assert len(files) == (2 if command == "find-imbalance-sc" else 1) + assert all(not pd.read_csv(file, sep="\t").empty for file in files) From e23a447bb0729108deb8bc242fceba9c8c62b4cd Mon Sep 17 00:00:00 2001 From: Jeff Jaureguy <67065808+Jaureguy760@users.noreply.github.com> Date: Tue, 13 Oct 2026 12:24:08 -0700 Subject: [PATCH 2/3] fix: constrain compatible runtime and CI dependency versions --- .github/workflows/ecosystem.yml | 2 +- environment.yml | 2 +- pixi.toml | 2 +- pyproject.toml | 2 +- requirements.txt | 2 +- rust/src/analysis.rs | 6 +++--- 6 files changed, 8 insertions(+), 8 deletions(-) diff --git a/.github/workflows/ecosystem.yml b/.github/workflows/ecosystem.yml index 55937c6..3855f20 100644 --- a/.github/workflows/ecosystem.yml +++ b/.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: | diff --git a/environment.yml b/environment.yml index bf52ff5..9e358df 100644 --- a/environment.yml +++ b/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) diff --git a/pixi.toml b/pixi.toml index 97265a5..a14060a 100644 --- a/pixi.toml +++ b/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" diff --git a/pyproject.toml b/pyproject.toml index ba4765e..1a520d6 100644 --- a/pyproject.toml +++ b/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", diff --git a/requirements.txt b/requirements.txt index 8a379a4..c098ce1 100644 --- a/requirements.txt +++ b/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 diff --git a/rust/src/analysis.rs b/rust/src/analysis.rs index 8b639da..822ba51 100644 --- a/rust/src/analysis.rs +++ b/rust/src/analysis.rs @@ -508,9 +508,9 @@ fn golden_section_search(f: F, a: f64, mut b: f64, tol: f64) -> Result 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; From b29b472752feb957c25f90dd7a0194fcec8351ba Mon Sep 17 00:00:00 2001 From: Jeff Jaureguy <67065808+Jaureguy760@users.noreply.github.com> Date: Tue, 13 Oct 2026 13:40:45 -0700 Subject: [PATCH 3/3] fix: match archived Python SNP counting priority --- CHANGELOG.md | 2 + rust/src/bam_counter.rs | 76 ++++++- tests/counting/fixtures/README.md | 14 ++ .../fixtures/archive_count_alleles.py.txt | 124 +++++++++++ tests/test_rust_python_counting_parity.py | 198 ++++++++++-------- 5 files changed, 311 insertions(+), 103 deletions(-) create mode 100644 tests/counting/fixtures/README.md create mode 100644 tests/counting/fixtures/archive_count_alleles.py.txt diff --git a/CHANGELOG.md b/CHANGELOG.md index a0295dc..0c3bfd1 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -5,6 +5,8 @@ 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 diff --git a/rust/src/bam_counter.rs b/rust/src/bam_counter.rs index 7cdaa44..39796c7 100644 --- a/rust/src/bam_counter.rs +++ b/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> { + 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::(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> { + let mut bam = bam::IndexedReader::from_path(&self.bam_path).map_err(|e| { + PyErr::new::(format!("Failed to open BAM: {e}")) + })?; + let mut seen_reads: FxHashSet> = FxHashSet::default(); + let mut counts = FxHashMap::default(); + for (idx, region) in regions { + let target = region.pos.checked_sub(1).ok_or_else(|| { + PyErr::new::("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::(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::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); } } diff --git a/tests/counting/fixtures/README.md b/tests/counting/fixtures/README.md new file mode 100644 index 0000000..3e08038 --- /dev/null +++ b/tests/counting/fixtures/README.md @@ -0,0 +1,14 @@ +`archive_count_alleles.py.txt` is a verbatim copy of +`src/counting/count_alleles.py` from WASP2 commit +`41354fc087404ba1fa99542803ee9cb9b7e0e6bf` (`archive/master-pre-v1.3.0`). + +The counting parity tests compile only `count_snp_alleles` from this source. +They execute its original SNP-ordered BAM fetches and chromosome-wide read-name +set without translating the algorithm or changing expected counts. + +The Rust bulk biallelic SNP path now uses the same per-SNP fetch order. This can +require more index seeks and repeated BAM decoding than a chromosome-span scan; +large-cohort performance has not been measured here. Non-SNP counting keeps its +existing implementation. Count storage remains UInt32 to avoid legacy UInt16 +overflow, and an explicitly requested base-quality threshold remains available; +archive comparisons use the default `min_qual=0`. diff --git a/tests/counting/fixtures/archive_count_alleles.py.txt b/tests/counting/fixtures/archive_count_alleles.py.txt new file mode 100644 index 0000000..8df6c7b --- /dev/null +++ b/tests/counting/fixtures/archive_count_alleles.py.txt @@ -0,0 +1,124 @@ +import timeit +from pathlib import Path +from bisect import bisect_left + +import polars as pl + +from pysam.libcalignmentfile import AlignmentFile + +# Helper that does binary search +def find_read_aln_pos(read, pos): + + aln_list = read.get_aligned_pairs(True) + + i = bisect_left(aln_list, pos, key=lambda x: x[1]) + + if i != len(aln_list) and aln_list[i][1] == pos: + return aln_list[i][0] + else: + return None + + +def make_count_df(bam_file, df): + """ + Make DF containing all intersections and allele counts + + :param str bam_file: Path to BAM file + :param DataFrame df: Dataframe of intersections, output from + parse_(intersect/gene)_df() + :return DataFrame: DataFrame of counts + """ + count_list = [] + + chrom_list = df.get_column("chrom").unique( + maintain_order=True) + + total_start = timeit.default_timer() + + with AlignmentFile(bam_file, "rb") as bam: + + for chrom in chrom_list: + chrom_df = df.filter(pl.col("chrom") == chrom) + + snp_list = chrom_df.select( + ["pos", "ref", "alt"]).unique( + subset=["pos"], maintain_order=True).iter_rows() + + start = timeit.default_timer() + + try: + count_list.extend(count_snp_alleles(bam, chrom, snp_list)) + except ValueError: + print(f"Skipping {chrom}: Contig not found\n") + else: + print(f"{chrom}: Counted {chrom_df.height} SNP's in {timeit.default_timer() - start:.2f} seconds!") + + + total_end = timeit.default_timer() + print(f"Counted all SNP's in {total_end - total_start:.2f} seconds!") + + # Previously used str as chrom instead of cat + chrom_enum = pl.Enum(df.get_column("chrom").cat.get_categories()) + + count_df = pl.DataFrame( + count_list, + schema={"chrom": chrom_enum, + "pos": pl.UInt32, + "ref_count": pl.UInt16, + "alt_count": pl.UInt16, + "other_count": pl.UInt16 + } + ) + + # possibly find better solution + df = df.with_columns([pl.col("chrom").cast(chrom_enum)] + ).join(count_df, on=["chrom", "pos"], how="left") + + # df = df.join(count_df, on=["chrom", "pos"], how="left") + + return df + + +def count_snp_alleles(bam, chrom, snp_list): + """ + Helper function called by... + make_count_df() + """ + + read_set = set() + allele_counts = [] + + for pos, ref, alt in snp_list: + + # read_set = set() + ref_count, alt_count, other_count = 0, 0, 0 + + # Got make sure read is not double counted + for read in bam.fetch(chrom, pos-1, pos): + + # If already counted allele + if read.query_name in read_set: + continue + + read_set.add(read.query_name) + + seq = read.query_sequence + + for qpos, refpos in read.get_aligned_pairs(True): + + # TODO Update with binary search + if refpos == pos-1: + + if seq[qpos] == ref: + ref_count+=1 + elif seq[qpos] == alt: + alt_count+=1 + else: + other_count+=1 + + # Found no longer need to loop + break + + allele_counts.append((chrom, pos, ref_count, alt_count, other_count)) + + return allele_counts \ No newline at end of file diff --git a/tests/test_rust_python_counting_parity.py b/tests/test_rust_python_counting_parity.py index b683f84..c7a7132 100644 --- a/tests/test_rust_python_counting_parity.py +++ b/tests/test_rust_python_counting_parity.py @@ -1,21 +1,11 @@ -""" -Rust vs Python parity test for allele counting. - -Reimplements the Rust BamCounter algorithm in pure Python + pysam and -compares allele counts at every variant position. The Python reference -matches the Rust semantics exactly: - -- Single BAM fetch per chromosome spanning all variant positions -- Each read assigned to the earliest-encounter-index SNP it overlaps -- Seen-reads set accumulates across positions (each read counted once) -- BAM flag filtering: only unmapped skipped (canonical WASP contract — - caller is responsible for BAM pre-filtering) +"""Bulk SNP counting parity against archived Python commit +41354fc087404ba1fa99542803ee9cb9b7e0e6bf. -This catches any numerical differences between the Rust and Python -implementations that golden-file tests (which compare Rust output to -Rust-generated baselines) would miss. +The fixture is the verbatim archived src/counting/count_alleles.py. Only its +count_snp_alleles function is loaded, avoiding unrelated legacy imports. """ +import ast from pathlib import Path import pysam @@ -27,84 +17,23 @@ # --------------------------------------------------------------------------- -# Python counting algorithm matching Rust BamCounter semantics +# Unmodified archived Python counting reference # --------------------------------------------------------------------------- def python_count_snp_alleles(bam_path: str, chrom: str, snp_list: list[tuple]): - """Count ref/alt/other alleles using pure Python + pysam. - - Reimplements the Rust BamCounter algorithm (bam_counter.rs:228-370) in - Python to enable exact parity comparison. Key semantics matched: - - 1. Single fetch for the entire chromosome span (min_pos to max_pos). - 2. Each read is assigned to the earliest-encounter-index SNP it has an - aligned base at. A read is counted at exactly one position. - 3. Seen-reads set accumulates across all positions (paired-end mates - and multi-overlap reads are counted only once per chromosome). - 4. BAM flag filtering: skip only unmapped (canonical WASP contract). - - Args: - bam_path: Path to BAM file. - chrom: Chromosome name. - snp_list: List of (pos_1based, ref, alt) tuples, in encounter order. - - Returns: - List of (chrom, pos_1based, ref_count, alt_count, other_count) tuples. - """ - if not snp_list: - return [] - - # Build position -> list of (encounter_index, ref, alt) - pos_map: dict[int, list[tuple]] = {} - for enc_idx, (pos, ref, alt) in enumerate(snp_list): - pos_map.setdefault(pos, []).append((enc_idx, ref, alt)) - - # Initialize counts per encounter index - counts: dict[int, list[int]] = {i: [0, 0, 0] for i in range(len(snp_list))} - - min_pos = min(pos for pos, _, _ in snp_list) - max_pos = max(pos for pos, _, _ in snp_list) - - seen_reads: set[str] = set() - + """Run the unmodified archived function, not a translation of Rust.""" + source = (ROOT / "tests/counting/fixtures/archive_count_alleles.py.txt").read_text() + tree = ast.parse(source) + function = next( + node + for node in tree.body + if isinstance(node, ast.FunctionDef) and node.name == "count_snp_alleles" + ) + namespace = {} + exec( + compile(ast.Module(body=[function], type_ignores=[]), "archived_python", "exec"), namespace + ) with pysam.AlignmentFile(bam_path, "rb") as bam: - for read in bam.fetch(chrom, max(0, min_pos - 1), max_pos + 1): - if read.is_unmapped: - continue - - qname = read.query_name - if qname in seen_reads: - continue - - seq = read.query_sequence - if seq is None: - continue - - # Find earliest-encounter-index SNP this read has an aligned base at - best = None # (encounter_idx, ref, alt, qpos) - for qpos, refpos in read.get_aligned_pairs(True): - if qpos is None or refpos is None: - continue - pos1 = refpos + 1 - if pos1 in pos_map: - for enc_idx, ref, alt in pos_map[pos1]: - if best is None or enc_idx < best[0]: - best = (enc_idx, ref, alt, qpos) - - if best is not None: - enc_idx, ref, alt, qpos = best - base = seq[qpos] - if base == ref: - counts[enc_idx][0] += 1 - elif base == alt: - counts[enc_idx][1] += 1 - else: - counts[enc_idx][2] += 1 - seen_reads.add(qname) - - return [ - (chrom, snp_list[i][0], counts[i][0], counts[i][1], counts[i][2]) - for i in range(len(snp_list)) - ] + return namespace["count_snp_alleles"](bam, chrom, snp_list) def parse_het_variants_from_vcf(vcf_path: str, sample: str | None = None): @@ -217,8 +146,7 @@ def sanity_data(self): def test_counting_parity_chr21(self, sanity_data): """Rust and Python must produce identical counts on real chr21 data. - This is the definitive parity test — 33k+ variants, real sequencing - data from HG00731. Any systematic bias would surface here. + Compare the available HG00731 data to the archived algorithm. """ import wasp2_rust @@ -270,3 +198,89 @@ def test_counting_parity_chr21(self, sanity_data): + "\n".join(mismatches[:20]) + (f"\n... and {len(mismatches) - 20} more" if len(mismatches) > 20 else "") ) + + +def _write_counting_bam(tmp_path, scenario): + bam_path = tmp_path / "archive_priority.bam" + # Records are coordinate sorted. The same name on chr2 must remain eligible. + records = [("chr1", 100, "shared", "A" * 20, [(0, 20)], 65)] + if scenario == "mates": + records.append(("chr1", 200, "shared", "G" * 20, [(0, 20)], 129)) + elif scenario == "spanning": + records = [("chr1", 100, "shared", "A" * 120, [(0, 120)], 0)] + elif scenario in {"deletion", "skip"}: + gap = 2 if scenario == "deletion" else 3 + records = [("chr1", 99, "shared", "A" * 10, [(0, 1), (gap, 100), (0, 9)], 0)] + elif scenario == "repeated": + records.append(("chr1", 100, "shared", "G" * 20, [(0, 20)], 0)) + records.append(("chr1", 200, "shared", "G" * 20, [(0, 20)], 0)) + elif scenario == "other": + records = [("chr1", 100, "shared", "N" + "A" * 19, [(0, 20)], 0)] + records.append(("chr2", 100, "shared", "G" * 20, [(0, 20)], 0)) + with pysam.AlignmentFile( + str(bam_path), + "wb", + header={ + "HD": {"VN": "1.6", "SO": "coordinate"}, + "SQ": [{"SN": "chr1", "LN": 1000}, {"SN": "chr2", "LN": 1000}], + }, + ) as bam: + for chrom, start, name, sequence, cigar, flag in records: + read = pysam.AlignedSegment() + read.query_name = name + read.query_sequence = sequence + read.flag = flag + read.reference_id = 0 if chrom == "chr1" else 1 + read.reference_start = start + read.mapping_quality = 60 + read.cigar = cigar + read.query_qualities = pysam.qualitystring_to_array("I" * len(sequence)) + if flag & 1: + read.next_reference_id = 0 + read.next_reference_start = 200 if flag == 65 else 100 + bam.write(read) + pysam.index(str(bam_path)) + return str(bam_path) + + +@pytest.mark.parametrize("scenario", ["mates", "spanning", "deletion", "skip", "repeated", "other"]) +@pytest.mark.parametrize("positions", [[101, 201], [201, 101]]) +@pytest.mark.parametrize("threads", [1, 2]) +def test_archived_snp_priority(tmp_path, scenario, positions, threads): + import wasp2_rust + + bam_path = _write_counting_bam(tmp_path, scenario) + snps = [(pos, "A", "G") for pos in positions] + expected = [row[2:] for row in python_count_snp_alleles(bam_path, "chr1", snps)] + expected += [row[2:] for row in python_count_snp_alleles(bam_path, "chr2", [(101, "A", "G")])] + regions = [("chr1", *snp) for snp in snps] + [("chr2", 101, "A", "G")] + observed = wasp2_rust.BamCounter(bam_path).count_alleles(regions, min_qual=0, threads=threads) + assert observed == expected + assert observed[-1] == (0, 1, 0), "Read-name deduplication must reset per chromosome" + if scenario in {"deletion", "skip"} and positions == [101, 201]: + assert observed[:2] == [(0, 0, 0), (0, 0, 0)] + + +def test_archived_counts_are_reused_for_overlapping_regions(tmp_path): + import polars as pl + + from counting.count_alleles import make_count_df + + bam_path = _write_counting_bam(tmp_path, "spanning") + variants = pl.DataFrame( + { + "chrom": ["chr1"] * 4, + "pos": [201, 101, 201, 101], + "ref": ["A"] * 4, + "alt": ["G"] * 4, + "region": ["peak_a", "peak_a", "peak_b", "peak_b"], + } + ).with_columns(pl.col("chrom").cast(pl.Categorical)) + expected = { + row[1]: row[2:] + for row in python_count_snp_alleles(bam_path, "chr1", [(201, "A", "G"), (101, "A", "G")]) + } + observed = make_count_df(bam_path, variants) + assert observed.height == 4 + for row in observed.iter_rows(named=True): + assert (row["ref_count"], row["alt_count"], row["other_count"]) == expected[row["pos"]]