From e595d805b6a870a80acb01574e7e432ed6bd1a1a Mon Sep 17 00:00:00 2001 From: Lucas Arnoldt Date: Tue, 25 Aug 2026 21:54:31 +0100 Subject: [PATCH 1/3] fixes, more tests --- src/cellink/_core/donordata.py | 19 ++-- src/cellink/io/_export.py | 7 ++ src/cellink/io/_pgen.py | 65 +++++++---- src/cellink/tl/external/_ld.py | 2 +- src/cellink/tl/external/_seismic_torch.py | 6 +- src/cellink/tl/external/_sldsc_utils.py | 2 +- tests/test_categorical_var_dtype.py | 122 ++++++++++++++++++++ tests/test_coloc.py | 130 ++++++++++++++++++++++ tests/test_donor_data.py | 12 ++ tests/test_io.py | 25 +++++ tests/test_scooby_resolve_snp.py | 66 ++++++++++- tests/test_seismic_torch.py | 122 ++++++++++++++++++++ tests/test_sldsc_utils.py | 107 ++++++++++++++++++ tests/test_tl_external.py | 68 +++++++++++ tests/test_tl_external_pc_ld_gsmap.py | 121 ++++++++++++++++++++ 15 files changed, 834 insertions(+), 40 deletions(-) create mode 100644 tests/test_categorical_var_dtype.py create mode 100644 tests/test_coloc.py create mode 100644 tests/test_seismic_torch.py create mode 100644 tests/test_sldsc_utils.py create mode 100644 tests/test_tl_external_pc_ld_gsmap.py diff --git a/src/cellink/_core/donordata.py b/src/cellink/_core/donordata.py index 75f6177..55c77c9 100644 --- a/src/cellink/_core/donordata.py +++ b/src/cellink/_core/donordata.py @@ -1,5 +1,6 @@ from __future__ import annotations +import copy as _copy import logging from collections.abc import Callable @@ -167,11 +168,13 @@ def _match_donors(self, G: AnnData | MuData, C: AnnData | MuData) -> None: self._G = G def copy(self) -> DonorData: - if self._G.is_view: - self._G = self._G.copy() - if self._C.is_view: - self._C = self._C.copy() - return self + new = DonorData.__new__(DonorData) + new._var_dims_to_sync = list(self._var_dims_to_sync) + new.donor_id = self.donor_id + new._G = self._G.copy() + new._C = self._C.copy() + new.uns = _copy.deepcopy(self.uns) + return new def _write_dd(self, f: h5py.File, zarr_path: str | None = None, x_chunks=None): is_zarr = isinstance(f, zarr.Group) @@ -312,10 +315,8 @@ def sel( C_obs: slice = slice(None), C_var: slice = slice(None), ): - _G = self.G[G_obs] - _G = _G[:, G_var] - _C = self.C[C_obs] - _C = _C[:, C_var] + _G = self.G[G_obs, G_var] + _C = self.C[C_obs, C_var] _G = self._sync_var_dims(_G, _C) return DonorData(G=_G, C=_C, donor_id=self.donor_id, var_dims_to_sync=self._var_dims_to_sync) diff --git a/src/cellink/io/_export.py b/src/cellink/io/_export.py index 621a38e..1e95caa 100644 --- a/src/cellink/io/_export.py +++ b/src/cellink/io/_export.py @@ -2,6 +2,7 @@ import sys import numpy as np +import pandas as pd import xarray as xr from anndata import AnnData from pandas_plink import write_plink1_bin @@ -72,6 +73,12 @@ def to_plink( if not output_prefix.endswith(".bed"): output_prefix += ".bed" + categorical_cols = [c for c in (chrom, a0, a1) if isinstance(gdata.var[c].dtype, pd.CategoricalDtype)] + if categorical_cols: + gdata = gdata.copy() + for c in categorical_cols: + gdata.var[c] = gdata.var[c].astype(str) + xarr = xr.DataArray( gdata.X.astype("float32", copy=False), dims=("sample", "variant"), diff --git a/src/cellink/io/_pgen.py b/src/cellink/io/_pgen.py index 28e6f95..8aa974f 100644 --- a/src/cellink/io/_pgen.py +++ b/src/cellink/io/_pgen.py @@ -29,16 +29,24 @@ def _is_matrix_elem(elem_name: str) -> bool: return elem_name.endswith("/X") or elem_name.rsplit("/", 1)[0].endswith("/layers") -def lazy_anndata_zarr_callback(func, elem_name: str, elem, iospec): + + +def lazy_anndata_zarr_callback(func, elem_name: str, elem, iospec, backend: Literal["dask", "zarr"] = "dask"): """``read_dispatched`` callback that reconstructs an AnnData (at any nesting depth, e.g. as ``G``/``C`` inside a larger DonorData zarr store) - while keeping a dense ``X``/``layers`` entry Dask-backed instead of + while keeping a dense ``X``/``layers`` entry lazily backed instead of materializing it. + + ``backend="dask"`` (default, unchanged behavior) wraps the dense array in + a Dask array. ``backend="zarr"`` instead returns the raw Zarr array + directly. """ if iospec.encoding_type == "anndata" or elem_name.endswith("/"): return ad.AnnData( **{ - k: read_dispatched(v, lazy_anndata_zarr_callback) + k: read_dispatched( + v, lambda f, n, e, iospec: lazy_anndata_zarr_callback(f, n, e, iospec, backend=backend) + ) for k, v in dict(elem).items() if not k.startswith("raw.") } @@ -51,22 +59,23 @@ def lazy_anndata_zarr_callback(func, elem_name: str, elem, iospec): ): return read_elem(elem) elif _is_matrix_elem(elem_name) and iospec.encoding_type == "array": - return da.from_zarr(elem) + return elem if backend == "zarr" else da.from_zarr(elem) else: return func(elem) -def read_pgen_zarr(store: str | Path) -> ad.AnnData: +def read_pgen_zarr(store: str | Path, backend: Literal["dask", "zarr"] = "dask") -> ad.AnnData: """ Lazily read an AnnData Zarr v3 store written by `stream_pgen_to_zarr`. This function reconstructs an :class:`anndata.AnnData` object from a Zarr - store while keeping the primary data matrix (`X`) backed by Dask arrays. - It is designed for large genotype matrices that cannot be loaded fully - into memory. + store while keeping the primary data matrix (`X`) lazily backed rather + than materializing it. It is designed for large genotype matrices that + cannot be loaded fully into memory. The reader preserves: - - Dense X stored as a Zarr array (returned as a Dask-backed array) + - Dense X stored as a Zarr array (backed by Dask, or the raw Zarr + array directly; see `backend`) - Sparse matrices (CSR/CSC) - DataFrames (obs, var) - Awkward arrays @@ -77,12 +86,17 @@ def read_pgen_zarr(store: str | Path) -> ad.AnnData: store : str or pathlib.Path Path to a Zarr directory created by `stream_pgen_to_zarr` or a compatible AnnData Zarr v3 store. + backend : {"dask", "zarr"} + How to back a dense `X`/`layers` entry. ``"dask"`` (default) wraps it in a Dask array, useful if + you need Dask's own lazy-graph chaining on X. ``"zarr"`` returns the + raw Zarr array directly instead. Returns ------- anndata.AnnData AnnData object with: - - `X` as a Dask-backed array (for dense storage) + - `X` as a Dask-backed or raw-Zarr-backed array (for dense storage), + per `backend` - `obs` and `var` as pandas DataFrames - empty container groups (`uns`, `obsm`, `varm`, `layers`, etc.) if present in the store @@ -90,9 +104,10 @@ def read_pgen_zarr(store: str | Path) -> ad.AnnData: Notes ----- - The returned object is **lazy** when X is dense. Computation is triggered - only when `.compute()` or in-memory materialization is requested. + only when `.compute()` (Dask backend) or ordinary indexing (Zarr + backend) is requested. - For sparse X written via `stream_pgen_to_zarr(..., sparse=True)`, - the matrix is loaded as a SciPy sparse matrix. + the matrix is loaded as a SciPy sparse matrix regardless of `backend`. - This function relies on AnnData's experimental dispatched I/O API. Examples @@ -101,12 +116,16 @@ def read_pgen_zarr(store: str | Path) -> ad.AnnData: >>> adata = cellink.io.read_pgen_zarr("genotypes.zarr") >>> adata AnnData object with n_obs x n_vars = ... + >>> fast = cellink.io.read_pgen_zarr("genotypes.zarr", backend="zarr") + >>> fast.X[donor_idx, :] # direct Zarr read, no Dask task overhead >>> # Trigger computation >>> X = adata.X.compute() """ f = zarr.open(str(store), mode="r") - return read_dispatched(f, callback=lazy_anndata_zarr_callback) + return read_dispatched( + f, callback=lambda func, name, elem, iospec: lazy_anndata_zarr_callback(func, name, elem, iospec, backend=backend) + ) def _read_pvar(pvar_file: Path) -> pd.DataFrame: @@ -141,12 +160,16 @@ def _read_pvar(pvar_file: Path) -> pd.DataFrame: pv = pv.rename(columns={k: v for k, v in rename_map.items() if k in pv.columns}) if VAnn.index in pv.columns: pv[VAnn.index] = pv[VAnn.index].astype(str) + and + if VAnn.chrom in pv.columns: - pv[VAnn.chrom] = pv[VAnn.chrom].astype(str) + pv[VAnn.chrom] = pv[VAnn.chrom].astype(str).astype("category") + if "FILTER" in pv.columns: + pv["FILTER"] = pv["FILTER"].astype(str).astype("category") if VAnn.a0 in pv.columns: - pv[VAnn.a0] = pv[VAnn.a0].astype(str) + pv[VAnn.a0] = pv[VAnn.a0].astype(str).astype("category") if VAnn.a1 in pv.columns: - pv[VAnn.a1] = pv[VAnn.a1].astype(str) + pv[VAnn.a1] = pv[VAnn.a1].astype(str).astype("category") return pv @@ -159,7 +182,7 @@ def stream_pgen_to_zarr( chunk_samples: int = 4096, chunk_variants: int = 2048, memory_limit_gb: float = 10.0, - compressor: str = "zstd", + compressor: str | None = "zstd", compression_level: int = 7, sparse: bool = False, sparse_format: Literal["csc", "csr"] = "csc", @@ -284,11 +307,7 @@ def _base(p: str) -> str: pvar.index = pvar.index.astype(str) output_path = Path(output_path) - blosc = BloscCodec( - cname=compressor, - clevel=compression_level, - shuffle=BloscShuffle.bitshuffle, - ) + codecs = () if compressor is None else (BloscCodec(cname=compressor, clevel=compression_level, shuffle=BloscShuffle.bitshuffle),) if sparse: if sparse_format not in ("csr", "csc"): @@ -356,7 +375,7 @@ def _base(p: str) -> str: shape=(n_samples, n_variants_total), chunks=(chunk_samples, chunk_variants), dtype="i1", - compressors=(blosc,), + compressors=codecs, ) Xz.attrs["encoding-type"] = "array" Xz.attrs["encoding-version"] = "0.2.0" diff --git a/src/cellink/tl/external/_ld.py b/src/cellink/tl/external/_ld.py index ac12412..a4aa5ec 100644 --- a/src/cellink/tl/external/_ld.py +++ b/src/cellink/tl/external/_ld.py @@ -54,7 +54,7 @@ def calculate_ld( plink_export_kwargs = {} if run and shutil.which("plink") is None: - raise ImportError("plink is required for `calculate_pcs`. Please install it.") + raise ImportError("plink is required for `calculate_ld`. Please install it.") if out is None: out = f"{prefix}_ld" diff --git a/src/cellink/tl/external/_seismic_torch.py b/src/cellink/tl/external/_seismic_torch.py index 03e3212..103401b 100644 --- a/src/cellink/tl/external/_seismic_torch.py +++ b/src/cellink/tl/external/_seismic_torch.py @@ -189,11 +189,7 @@ def forward(self, G: torch.Tensor, verbose: bool = False, return_all: bool = Fal if return_all: pval_two_sided = torch.tensor(st.chi2(1).sf(self.lrt.cpu().data.numpy()), device=self.F.device) pval_one_sided = torch.where(self.beta_g > 0, pval_two_sided / 2.0, 1.0 - (pval_two_sided / 2.0)) - z = np.sign(self.beta_g.cpu().data.numpy()) * np.sqrt( - st.chi2.ppf(1.0 - pval_two_sided.cpu().data.numpy(), df=1) - ) - z = torch.tensor(z, device=self.F.device) - ste = self.beta_g / z + ste = torch.sqrt(self.s2 * n[:, None]) return nll, pval_one_sided, self.beta_g, ste return nll diff --git a/src/cellink/tl/external/_sldsc_utils.py b/src/cellink/tl/external/_sldsc_utils.py index b788a04..1b5afaf 100644 --- a/src/cellink/tl/external/_sldsc_utils.py +++ b/src/cellink/tl/external/_sldsc_utils.py @@ -549,7 +549,7 @@ def _pick_var_col(adata: AnnData, candidates: list[str], default: str | None) -> def _normalize_chromosome(chr_series: pd.Series) -> pd.Series: """Normalize chromosome labels to standard format.""" - normalized = chr_series.astype(str).str.replace("^chr", "", regex=True).str.upper() + normalized = chr_series.astype(str).str.replace("^chr", "", regex=True, case=False).str.upper() return normalized.str.extract(r"^([0-9XYM]+)", expand=False) diff --git a/tests/test_categorical_var_dtype.py b/tests/test_categorical_var_dtype.py new file mode 100644 index 0000000..bc0433b --- /dev/null +++ b/tests/test_categorical_var_dtype.py @@ -0,0 +1,122 @@ +import numpy as np +import pandas as pd +import pytest + +from cellink._core.dummy_data import sim_gdata +from cellink.io._export import write_variants_to_vcf +from cellink.tl._subset_region import subset_genomic_region + + +def _str_and_categorical_gdata(): + """Two AnnDatas with byte-identical underlying data, differing only in + whether chrom/a0/a1 are plain string or category dtype -- sim_gdata() + draws random alleles internally, so calling it twice independently (as + an earlier version of this fixture did) compares two different random + genotypes, not the same data under two dtypes. + """ + gdata_str = sim_gdata(n_donors=20, n_snps=30) + gdata_str.var["chrom"] = gdata_str.var["chrom"].astype(str) + gdata_cat = gdata_str.copy() + gdata_cat.var["chrom"] = gdata_cat.var["chrom"].astype("category") + gdata_cat.var["a0"] = gdata_cat.var["a0"].astype("category") + gdata_cat.var["a1"] = gdata_cat.var["a1"].astype("category") + return gdata_str, gdata_cat + + +def test_subset_genomic_region_matches_plain_string_dtype(): + gdata_str, gdata_cat = _str_and_categorical_gdata() + + start, end = int(gdata_str.var["pos"].min()), int(gdata_str.var["pos"].max()) + 1 + sub_str = subset_genomic_region(gdata_str, chrom="1", start=start, end=end) + sub_cat = subset_genomic_region(gdata_cat, chrom="1", start=start, end=end) + + assert sub_str.shape == sub_cat.shape + assert list(sub_str.var.index) == list(sub_cat.var.index) + + +def test_np_unique_and_equality_on_categorical_chrom(): + _, gdata = _str_and_categorical_gdata() + uniq = np.unique(gdata.var["chrom"]) + assert list(uniq) == ["1"] + mask = gdata.var["chrom"] == "1" + assert mask.all() + assert type(gdata.var["chrom"].iloc[0]) is str + + +def test_write_variants_to_vcf_identical_with_categorical(tmp_path): + gdata_str, gdata_cat = _str_and_categorical_gdata() + + out_str, out_cat = tmp_path / "str.vcf", tmp_path / "cat.vcf" + write_variants_to_vcf(gdata_str, out_file=str(out_str)) + write_variants_to_vcf(gdata_cat, out_file=str(out_cat)) + assert out_str.read_text() == out_cat.read_text() + + +def test_to_plink_roundtrip_identical_with_categorical(tmp_path): + bed_reader = pytest.importorskip("bed_reader") + from cellink.io._export import to_plink + + gdata_str, gdata_cat = _str_and_categorical_gdata() + gdata_str.obs["donor_id"] = gdata_str.obs.index + gdata_str.obs["sex"] = 0 + gdata_cat.obs["donor_id"] = gdata_cat.obs.index + gdata_cat.obs["sex"] = 0 + + prefix_str, prefix_cat = str(tmp_path / "str"), str(tmp_path / "cat") + to_plink(gdata_str, output_prefix=prefix_str) + to_plink(gdata_cat, output_prefix=prefix_cat) + + b_str = bed_reader.open_bed(prefix_str + ".bed") + b_cat = bed_reader.open_bed(prefix_cat + ".bed") + np.testing.assert_array_equal(b_str.read(), b_cat.read()) + assert list(b_str.chromosome) == list(b_cat.chromosome) + assert list(b_str.allele_1) == list(b_cat.allele_1) + assert list(b_str.allele_2) == list(b_cat.allele_2) + + +def test_tensorqtl_input_generator_cis_matches_plain_string_dtype(): + """The real integration risk flagged for this change: cellink's + run_tensorqtl(use_python_api=True) hands variant_df straight to + tensorqtl's own genotypeio.InputGeneratorCis, which does + variant_df['chrom'].unique() / .groupby('chrom') / membership checks -- + code cellink does not control. Verify categorical dtype produces + byte-identical cis-window results there directly, not just in cellink's + own call sites. + """ + tensorqtl_genotypeio = pytest.importorskip("tensorqtl.genotypeio") + + rng = np.random.default_rng(0) + n_var, n_genes, n_samples = 60, 6, 12 + chrom_str = np.array(["1"] * 30 + ["2"] * 30) + pos = np.concatenate([np.sort(rng.choice(1_000_000, 30, replace=False)), + np.sort(rng.choice(1_000_000, 30, replace=False))]) + variant_df_str = pd.DataFrame({"chrom": chrom_str, "pos": pos}, + index=[f"var{i}" for i in range(n_var)]) + variant_df_cat = variant_df_str.copy() + variant_df_cat["chrom"] = variant_df_cat["chrom"].astype("category") + + phenotype_pos_df = pd.DataFrame({ + "chr": np.array(["1"] * 3 + ["2"] * 3), + "start": np.sort(rng.choice(1_000_000, n_genes, replace=False)), + }, index=[f"gene{i}" for i in range(n_genes)]) + phenotype_pos_df["end"] = phenotype_pos_df["start"] + 1000 + + genotype_df = pd.DataFrame(rng.integers(0, 3, size=(n_var, n_samples)), + index=variant_df_str.index, columns=[f"s{i}" for i in range(n_samples)]) + phenotype_df = pd.DataFrame(rng.normal(size=(n_genes, n_samples)), + index=phenotype_pos_df.index, columns=[f"s{i}" for i in range(n_samples)]) + + def cis_ranges_for(variant_df): + gen = tensorqtl_genotypeio.InputGeneratorCis( + genotype_df, variant_df, phenotype_df, phenotype_pos_df, window=1_000_000 + ) + return dict(gen.cis_ranges), gen.chrs, gen.phenotype_df.index.tolist() + + ranges_str, chrs_str, kept_str = cis_ranges_for(variant_df_str) + ranges_cat, chrs_cat, kept_cat = cis_ranges_for(variant_df_cat) + + assert chrs_str == chrs_cat + assert kept_str == kept_cat + assert set(ranges_str) == set(ranges_cat) + for k in ranges_str: + np.testing.assert_array_equal(ranges_str[k], ranges_cat[k]) diff --git a/tests/test_coloc.py b/tests/test_coloc.py new file mode 100644 index 0000000..70e2096 --- /dev/null +++ b/tests/test_coloc.py @@ -0,0 +1,130 @@ +from __future__ import annotations + +import numpy as np +import pytest + +from cellink.tl import coloc_abf, coloc_susie +from cellink.tl._coloc import DEFAULT_PRIOR_VAR, _combine_log_abf, _log_abf + +rpy2 = pytest.importorskip("rpy2", reason="rpy2 not installed") +robjects = pytest.importorskip("rpy2.robjects", reason="rpy2 not installed") +from rpy2.robjects.packages import PackageNotInstalledError, importr # noqa: E402 + +try: + coloc_r = importr("coloc") +except PackageNotInstalledError: + pytest.skip("R package 'coloc' not installed", allow_module_level=True) + + +def _r_log_abf(beta: np.ndarray, se: np.ndarray, prior_var: float) -> np.ndarray: + """R coloc's own `approx.bf.estimates(z, V, type="quant", sdY=sqrt(prior_var)/0.15)`, + called with sdY chosen so that R's internal `sd.prior = 0.15 * sdY` equals + `sqrt(prior_var)` exactly, matching cellink's `prior_var` parameterization.""" + z = np.asarray(beta) / np.asarray(se) + v = np.asarray(se) ** 2 + sdy = float(np.sqrt(prior_var) / 0.15) + out = coloc_r.approx_bf_estimates( + z=robjects.FloatVector(z), V=robjects.FloatVector(v), type="quant", sdY=robjects.FloatVector([sdy]) + ) + labf_col = list(out.names).index("lABF") + return np.asarray(out[labf_col]) + + +def _r_combine_abf(l1: np.ndarray, l2: np.ndarray, p1: float, p2: float, p12: float) -> dict[str, float]: + out = coloc_r.combine_abf( + robjects.FloatVector(np.asarray(l1)), + robjects.FloatVector(np.asarray(l2)), + robjects.FloatVector([p1]), + robjects.FloatVector([p2]), + robjects.FloatVector([p12]), + quiet=True, + ) + # R names its output PP.H0.abf..PP.H4.abf; remap to PP0..PP4 to match cellink's keys + r_to_py_key = {f"PP.H{i}.abf": f"PP{i}" for i in range(5)} + return {r_to_py_key[name]: float(out.rx2(name)[0]) for name in out.names} + + +@pytest.fixture +def rng(): + return np.random.default_rng(0) + + +def test_log_abf_matches_r_coloc(rng): + n = 200 + beta = rng.normal(0, 0.03, n) + se = rng.uniform(0.01, 0.06, n) + py_labf = _log_abf(beta, se, DEFAULT_PRIOR_VAR) + r_labf = _r_log_abf(beta, se, DEFAULT_PRIOR_VAR) + np.testing.assert_allclose(py_labf, r_labf, rtol=1e-10, atol=1e-12) + + +@pytest.mark.parametrize("p1,p2,p12", [(1e-4, 1e-4, 1e-5), (5e-5, 2e-4, 1e-6)]) +def test_combine_log_abf_matches_r_coloc(rng, p1, p2, p12): + n = 150 + l1 = rng.normal(0, 3, n) + l2 = rng.normal(0, 3, n) + # give a handful of SNPs a real shared/distinct signal so all 5 hypotheses + # get non-negligible mass, not just a diffuse null + l1[10] += 15 + l2[10] += 12 + l1[50] += 20 + l2[90] += 18 + + py_pp = _combine_log_abf(l1, l2, p1, p2, p12) + r_pp = _r_combine_abf(l1, l2, p1, p2, p12) + + assert set(py_pp) == {"PP0", "PP1", "PP2", "PP3", "PP4"} + for key in py_pp: + assert py_pp[key] == pytest.approx(r_pp[key], abs=1e-8), f"{key}: py={py_pp[key]!r} r={r_pp[key]!r}" + assert sum(py_pp.values()) == pytest.approx(1.0, abs=1e-10) + + +def test_coloc_abf_end_to_end_matches_r_coloc(rng): + """Full coloc_abf pipeline (beta/se -> log-ABF -> H0-H4), not just its + two halves tested in isolation above.""" + n = 100 + beta1 = rng.normal(0, 0.02, n) + se1 = rng.uniform(0.015, 0.05, n) + beta2 = rng.normal(0, 0.02, n) + se2 = rng.uniform(0.015, 0.05, n) + # SNP 7: a real shared causal signal, strong in both studies + beta1[7], se1[7] = 0.6, 0.04 + beta2[7], se2[7] = 0.5, 0.05 + + py_pp = coloc_abf(beta1, se1, beta2, se2) + + l1 = _r_log_abf(beta1, se1, DEFAULT_PRIOR_VAR) + l2 = _r_log_abf(beta2, se2, DEFAULT_PRIOR_VAR) + r_pp = _r_combine_abf(l1, l2, p1=1e-4, p2=1e-4, p12=1e-5) + + for key in py_pp: + assert py_pp[key] == pytest.approx(r_pp[key], abs=1e-6) + assert py_pp["PP4"] > 0.9 # sanity: the planted shared signal should dominate + + +def test_coloc_susie_pairwise_combination_matches_r_coloc(rng): + """coloc_susie delegates each (signal, signal) pair to the same + _combine_log_abf core coloc_abf uses; check that per-pair result + against R directly rather than trusting the delegation blindly.""" + n = 40 + lbf1 = rng.normal(0, 1, (2, n)) + lbf2 = rng.normal(0, 1, (2, n)) + lbf1[0, 15] = 25.0 # effect 0 in study 1 points at SNP 15 + lbf2[0, 15] = 22.0 # effect 0 in study 2 also points at SNP 15: real colocalization + lbf1[1, 30] = 20.0 # effect 1 in study 1 points at SNP 30, no match in study 2 + + res = coloc_susie(lbf1, lbf2, cs1_index=[0, 1], cs2_index=[0]) + + for _, row in res.iterrows(): + l1 = lbf1[int(row["idx1"])] + l2 = lbf2[int(row["idx2"])] + r_pp = _r_combine_abf(l1, l2, p1=1e-4, p2=1e-4, p12=1e-5) + for key in ["PP0", "PP1", "PP2", "PP3", "PP4"]: + assert row[key] == pytest.approx(r_pp[key], abs=1e-6), f"idx1={row['idx1']} idx2={row['idx2']} {key}" + + # the matched pair (0, 0) should show real colocalization; the unmatched + # pair (1, 0) should not + hit00 = res[(res["idx1"] == 0) & (res["idx2"] == 0)].iloc[0] + hit10 = res[(res["idx1"] == 1) & (res["idx2"] == 0)].iloc[0] + assert hit00["PP4"] > 0.9 + assert hit10["PP4"] < hit00["PP4"] diff --git a/tests/test_donor_data.py b/tests/test_donor_data.py index 56874ca..8bdeae7 100644 --- a/tests/test_donor_data.py +++ b/tests/test_donor_data.py @@ -67,6 +67,18 @@ def test_slice_all(adata, gdata): assert dd.shape == (1, 1, 1, 1) +def test_copy_returns_independent_object(adata, gdata): + dd = DonorData(G=gdata, C=adata) + dd_copy = dd.copy() + assert dd_copy is not dd + assert dd_copy.G is not dd.G + assert dd_copy.C is not dd.C + dd_copy.uns["marker"] = "mutated_via_copy" + dd_copy.G.obs["new_col"] = 0 + assert "marker" not in dd.uns + assert "new_col" not in dd.G.obs + + def test_donordata_aggregate(adata, gdata, dummy_covariates): previous_adata_shape = adata.shape diff --git a/tests/test_io.py b/tests/test_io.py index d39c0a2..1c46d81 100644 --- a/tests/test_io.py +++ b/tests/test_io.py @@ -107,6 +107,31 @@ def test_stream_pgen_to_zarr_roundtrip(tmp_path): np.testing.assert_array_equal(X, expected) +@pytest.mark.slow +def test_stream_pgen_to_zarr_no_compression(tmp_path): + """compressor=None writes a plain, uncompressed Zarr array -- the fair + comparison point against PGEN's own storage_mode=0x10 (fixed 2-bit, no + compression, no difference lists): both stores should decode to the same + genotypes, and skipping BloscCodec entirely (not clevel=0, which still + pays codec framing overhead) should not change any value. + """ + pytest.importorskip("pgenlib") + from cellink.io import read_pgen_zarr, stream_pgen_to_zarr + + pgen_file = DATA / "simulated_genotype_calls.pgen" + out_path = tmp_path / "pgen_nocompression.zarr" + + stream_pgen_to_zarr(str(pgen_file), str(out_path), chunk_samples=50, chunk_variants=200, compressor=None) + reloaded = read_pgen_zarr(str(out_path), backend="zarr") + assert reloaded.X.compressors == () + X = np.asarray(reloaded.X[:]) + + out_path_compressed = tmp_path / "pgen_compressed.zarr" + stream_pgen_to_zarr(str(pgen_file), str(out_path_compressed), chunk_samples=50, chunk_variants=200) + X_compressed = np.asarray(read_pgen_zarr(str(out_path_compressed), backend="zarr").X[:]) + np.testing.assert_array_equal(X, X_compressed) + + @pytest.mark.slow def test_stream_pgen_to_zarr_sparse(tmp_path): pytest.importorskip("pgenlib") diff --git a/tests/test_scooby_resolve_snp.py b/tests/test_scooby_resolve_snp.py index b181f32..5ecbfe0 100644 --- a/tests/test_scooby_resolve_snp.py +++ b/tests/test_scooby_resolve_snp.py @@ -1,8 +1,58 @@ from __future__ import annotations +import sys +import types +from dataclasses import dataclass +from typing import Literal + +import numpy as np import pytest -pytest.importorskip("embpy", reason="resolve_snp_and_exon_bins needs embpy, install with `pip install cellink[embpy]`") +try: + import embpy.tl.genomics as _real_genomics # noqa: F401 + + _HAS_REAL_EMBPY = True +except ImportError: + _HAS_REAL_EMBPY = False + + +if not _HAS_REAL_EMBPY: + @dataclass + class _SNPContext: + position: int + ref_allele: str + alt_alleles: list + context_window: int = 512 + chrom: str = "" + strand: Literal["+", "-"] = "+" + variant_id: str = "" + + def __post_init__(self) -> None: + self.ref_allele = self.ref_allele.upper() + self.alt_alleles = [a.upper() for a in self.alt_alleles] + + def _genomic_to_bin_indices(intervals, window_start, bin_size, profile_offset_bp=0, num_bins=None): + bins = set() + for start, end in intervals: + rel_start = start - window_start - profile_offset_bp + rel_end = end - window_start - profile_offset_bp + b0 = int(rel_start // bin_size) + b1 = -(-int(rel_end) // bin_size) + for b in range(max(b0, 0), b1): + if num_bins is None or 0 <= b < num_bins: + bins.add(b) + return np.array(sorted(bins), dtype=int) + + _fake_embpy = types.ModuleType("embpy") + _fake_embpy_tl = types.ModuleType("embpy.tl") + _fake_embpy_genomics = types.ModuleType("embpy.tl.genomics") + _fake_embpy_genomics.SNPContext = _SNPContext + _fake_embpy_genomics.genomic_to_bin_indices = _genomic_to_bin_indices + _fake_embpy_tl.genomics = _fake_embpy_genomics + _fake_embpy.tl = _fake_embpy_tl + sys.modules.setdefault("embpy", _fake_embpy) + sys.modules.setdefault("embpy.tl", _fake_embpy_tl) + sys.modules.setdefault("embpy.tl.genomics", _fake_embpy_genomics) from cellink.tl.external import resolve_snp_and_exon_bins # noqa: E402 @@ -54,3 +104,17 @@ def test_skips_when_no_exon_bins_overlap(): spurious zero effect.""" result = resolve_snp_and_exon_bins(a0="A", a1="G", exon_intervals=[], **COMMON_KWARGS) assert result is None + + +def test_bin_indices_clipped_to_num_bins(): + """Regression test for the unclipped-bin-index bug: exons far outside the model's + cropped profile region must be dropped, not returned as out-of-range indices that + would crash downstream indexing.""" + result = resolve_snp_and_exon_bins( + a0="A", a1="G", exon_intervals=[(POS - 10, POS + 10), (POS + 100_000, POS + 100_020)], + **{**COMMON_KWARGS, "num_bins": 150}, + ) + assert result is not None + _, bins = result + assert len(bins) == 20 # only the near exon's bins survive + assert bins.max() < 150 diff --git a/tests/test_seismic_torch.py b/tests/test_seismic_torch.py new file mode 100644 index 0000000..66ee409 --- /dev/null +++ b/tests/test_seismic_torch.py @@ -0,0 +1,122 @@ +from __future__ import annotations + +import numpy as np +import statsmodels.api as sm +import torch +from anndata import AnnData +from scipy.sparse import csr_matrix + +from cellink.tl.external._seismic_torch import RegressionNLL, SparseScore, _adata_to_sparse_csr_tensor + + +def test_adata_to_sparse_csr_tensor_from_sparse_layer(): + X = csr_matrix(np.array([[1.0, 0.0, 3.0], [0.0, 2.0, 0.0]], dtype=np.float32)) + adata = AnnData(X) + t = _adata_to_sparse_csr_tensor(adata, layer=None) + assert t.layout == torch.sparse_csr + dense = t.to_dense().numpy() + np.testing.assert_allclose(dense, X.toarray()) + + +def test_adata_to_sparse_csr_tensor_from_dense_falls_back_to_dense_tensor(): + X = np.array([[1.0, 2.0], [3.0, 4.0]], dtype=np.float32) + adata = AnnData(X) + t = _adata_to_sparse_csr_tensor(adata, layer=None) + assert not t.is_sparse + np.testing.assert_allclose(t.numpy(), X) + + +def test_regression_nll_matches_statsmodels_ols(): + """Cross-validate RegressionNLL's closed-form likelihood-ratio test against + an independent OLS fit + LRT computed by hand from statsmodels, rather than + trusting the closed-form linear-algebra implementation on its own.""" + rng = np.random.default_rng(0) + n = 200 + z = rng.normal(size=n).astype(np.float64) + g = 0.6 * z + rng.normal(scale=0.5, size=n) # correlated covariate + g = g.astype(np.float64) + + model = RegressionNLL(torch.tensor(z)) + nll, pval, beta, se = model.forward(torch.tensor(g[:, None]), return_all=True) + + # independent reference: full model vs. intercept-only null, via statsmodels + X_full = sm.add_constant(g) + fit_full = sm.OLS(z, X_full).fit() + fit_null = sm.OLS(z, np.ones((n, 1))).fit() + ref_beta = fit_full.params[1] + ref_se = fit_full.bse[1] + + np.testing.assert_allclose(beta.item(), ref_beta, rtol=1e-4) + # small, expected gap vs. statsmodels' unbiased-variance convention, not a bug + np.testing.assert_allclose(se.item(), ref_se, rtol=5e-3) + # RegressionNLL's one-sided p-value, doubled, should recover statsmodels' own + # two-sided p-value for whichever tail the observed effect sign falls on + ref_p_two_sided = fit_full.pvalues[1] + doubled = min(2 * pval.item(), 1.0) + np.testing.assert_allclose(doubled, ref_p_two_sided, rtol=1e-2, atol=1e-6) + + +def test_regression_nll_multiple_columns_independent_of_each_other(): + """Each column of G is regressed independently; a strongly-associated + column and a pure-noise column in the same call must not contaminate + each other's p-value.""" + rng = np.random.default_rng(1) + n = 300 + z = rng.normal(size=n) + strong = 0.8 * z + rng.normal(scale=0.3, size=n) + noise = rng.normal(size=n) + G = np.stack([strong, noise], axis=1) + + model = RegressionNLL(torch.tensor(z)) + _, pval, beta, _ = model.forward(torch.tensor(G), return_all=True) + + assert pval[0, 0].item() < 0.001 # strong signal, should be highly significant + assert pval[1, 0].item() > 0.05 # pure noise, should not be + + +def test_sparse_score_perfectly_specific_gene_scores_near_one_in_its_own_celltype(): + """A gene expressed only in cluster A (zero elsewhere) should score near + 1.0 for cluster A and near 0 for cluster B under seismic's specificity + z-test, the basic sanity check that the sparse closed-form matches the + algorithm's own definition.""" + # 6 cells: 3 in cluster A (gene highly expressed), 3 in cluster B (gene ~0) + E = torch.tensor( + [ + [5.0, 0.0], + [6.0, 0.1], + [4.0, 0.0], + [0.0, 5.0], + [0.1, 6.0], + [0.0, 4.0], + ] + ) + masks = torch.tensor( + [ + [1.0, 0.0], + [1.0, 0.0], + [1.0, 0.0], + [0.0, 1.0], + [0.0, 1.0], + [0.0, 1.0], + ] + ) + scorer = SparseScore(E) + s = scorer.forward(masks) + assert s.shape == (2, 2) # [genes, celltypes] + # gene 0 (col 0 of E) is high in cluster A (mask col 0) -> should score highest there + assert s[0, 0] > s[0, 1] + # gene 1 is high in cluster B -> should score highest there + assert s[1, 1] > s[1, 0] + + +def test_sparse_score_uniform_expression_gives_low_specificity_everywhere(): + """A gene expressed identically across both clusters has no real + specificity signal; its score should not spike in either cluster.""" + E = torch.tensor([[5.0], [5.0], [5.0], [5.0], [5.0], [5.0]]) + masks = torch.tensor( + [[1.0, 0.0], [1.0, 0.0], [1.0, 0.0], [0.0, 1.0], [0.0, 1.0], [0.0, 1.0]] + ) + scorer = SparseScore(E) + s = scorer.forward(masks) + # both cluster scores should be small and roughly comparable (no specificity) + assert abs(s[0, 0].item() - s[0, 1].item()) < 0.3 diff --git a/tests/test_sldsc_utils.py b/tests/test_sldsc_utils.py new file mode 100644 index 0000000..1ddf9d6 --- /dev/null +++ b/tests/test_sldsc_utils.py @@ -0,0 +1,107 @@ +from __future__ import annotations + +import numpy as np +import pandas as pd +import pytest +from anndata import AnnData + +from cellink.tl.external._sldsc_utils import ( + _compute_celltype_means, + _compute_specificity, + _normalize_chromosome, + _pick_var_col, + _safe_filename, +) + + +@pytest.mark.parametrize( + "raw,expected", + [ + ("chr1", "1"), + ("CHR1", "1"), + ("1", "1"), + ("chrX", "X"), + ("chrY", "Y"), + ("chrM", "M"), + ("22", "22"), + ], +) +def test_normalize_chromosome_single(raw, expected): + result = _normalize_chromosome(pd.Series([raw])) + assert result.iloc[0] == expected + + +def test_normalize_chromosome_mixed_series(): + result = _normalize_chromosome(pd.Series(["chr1", "CHR2", "chrX", "10"])) + assert list(result) == ["1", "2", "X", "10"] + + +def test_pick_var_col_prefers_default_when_present(): + adata = AnnData(np.zeros((2, 2)), var=pd.DataFrame(index=["g1", "g2"], data={"symbol": ["A", "B"], "gene_name": ["A", "B"]})) + assert _pick_var_col(adata, candidates=["gene_name"], default="symbol") == "symbol" + + +def test_pick_var_col_falls_back_to_candidates(): + adata = AnnData(np.zeros((2, 2)), var=pd.DataFrame(index=["g1", "g2"], data={"gene_name": ["A", "B"]})) + assert _pick_var_col(adata, candidates=["gene_name", "symbol"], default="missing_col") == "gene_name" + + +def test_pick_var_col_returns_default_when_nothing_found(): + adata = AnnData(np.zeros((2, 2)), var=pd.DataFrame(index=["g1", "g2"])) + assert _pick_var_col(adata, candidates=["gene_name"], default="fallback") == "fallback" + + +def test_compute_celltype_means_basic(): + X = np.array([[1.0, 3.0], [3.0, 1.0], [10.0, 20.0]]) + adata = AnnData(X, obs=pd.DataFrame({"celltype": ["A", "A", "B"]}), var=pd.DataFrame(index=["g1", "g2"])) + means = _compute_celltype_means(adata, "celltype") + assert list(means.columns) == ["A", "B"] + np.testing.assert_allclose(means.loc["g1", "A"], 2.0) + np.testing.assert_allclose(means.loc["g2", "A"], 2.0) + np.testing.assert_allclose(means.loc["g1", "B"], 10.0) + + +def test_compute_celltype_means_empty_celltype_is_nan_not_zero(): + """Regression check: a cell type with 0 cells (e.g. after upstream filtering + drops every cell of that type) must be reported as NaN, not silently + scored as 0 expression, which would corrupt the specificity computation.""" + X = np.array([[1.0], [3.0]]) + adata = AnnData( + X, + obs=pd.DataFrame({"celltype": pd.Categorical(["A", "A"], categories=["A", "B"])}), + var=pd.DataFrame(index=["g1"]), + ) + means = _compute_celltype_means(adata, "celltype") + assert means.loc["g1", "A"] == 2.0 + assert np.isnan(means.loc["g1", "B"]) + + +def test_compute_specificity_sums_to_one_across_celltypes(): + mean_expr = pd.DataFrame({"A": [2.0, 0.0], "B": [8.0, 0.0]}, index=["g1", "g2"]) + spec = _compute_specificity(mean_expr) + np.testing.assert_allclose(spec.loc["g1", "A"], 0.2) + np.testing.assert_allclose(spec.loc["g1", "B"], 0.8) + np.testing.assert_allclose(spec.loc["g1"].sum(), 1.0) + + +def test_compute_specificity_zero_total_expression_is_zero_not_nan(): + """A gene expressed nowhere has no well-defined specificity; the + implementation must return 0, not propagate NaN into downstream LD-score + annotation files (which would silently corrupt them).""" + mean_expr = pd.DataFrame({"A": [0.0], "B": [0.0]}, index=["g1"]) + spec = _compute_specificity(mean_expr) + assert spec.loc["g1", "A"] == 0.0 + assert spec.loc["g1", "B"] == 0.0 + + +@pytest.mark.parametrize( + "raw,expected", + [ + ("CD4+ T cell", "CD4+_T_cell"), + ("Natural Killer (NK)", "Natural_Killer__NK_"), + (" leading/trailing ", "leading_trailing"), + ("already_safe-1.0", "already_safe-1.0"), + ], +) +def test_safe_filename(raw, expected): + assert _safe_filename(raw) == expected diff --git a/tests/test_tl_external.py b/tests/test_tl_external.py index 47654a9..a9b203b 100644 --- a/tests/test_tl_external.py +++ b/tests/test_tl_external.py @@ -1,7 +1,9 @@ import numpy as np import pandas as pd +import pytest from cellink.tl.external import JointNMFWrapper, compute_escore, scores_to_covar, scores_to_gmt +from cellink.tl.external._tensorqtl import build_known_cis_eqtls_from_tensorqtl def test_scores_to_gmt(tmp_path): @@ -97,3 +99,69 @@ def test_joint_nmf_wrapper(): assert wrapper.Hd.shape == (3, 12) assert (wrapper.Wh >= 0).all() assert (wrapper.Hh >= 0).all() + + +def _write_tensorqtl_nominal_parquet(path, rows): + pd.DataFrame(rows, columns=["gene", "variant_id", "pval"]).to_parquet(path) + return path + + +def test_build_known_cis_eqtls_picks_lowest_pval_per_gene(tmp_path): + parquet = _write_tensorqtl_nominal_parquet( + tmp_path / "nominal.parquet", + [ + ("ENSG1", "1:100:A:G", 0.5), + ("ENSG1", "1:200:A:G", 0.01), # lowest p for ENSG1 + ("ENSG2", "2:100:A:G", 0.2), + ("ENSG2", "2:200:A:G", 0.9), + ], + ) + known = build_known_cis_eqtls_from_tensorqtl(str(parquet), gene_names=["ENSG1", "ENSG2"]) + + # only the winning (lowest-pval) SNP per gene becomes a row at all: this is a + # sparse "selected loci" matrix, not a full candidate universe with explicit 0s + assert known.shape == (2, 2) + assert known.loc["1:200:A:G", "ENSG1"] == 1 + assert "1:100:A:G" not in known.index # the higher-pval candidate never appears + assert known.loc["2:100:A:G", "ENSG2"] == 1 + assert known.values.sum() == 2 # exactly one selected SNP per gene + + +def test_build_known_cis_eqtls_respects_pval_threshold(tmp_path): + parquet = _write_tensorqtl_nominal_parquet( + tmp_path / "nominal.parquet", + [ + ("ENSG1", "1:100:A:G", 0.5), + ("ENSG2", "2:100:A:G", 0.2), + ], + ) + # ENSG1's only variant fails the threshold and must not appear at all + known = build_known_cis_eqtls_from_tensorqtl( + str(parquet), gene_names=["ENSG1", "ENSG2"], pval_threshold=0.3, + ) + assert "ENSG1" not in known.columns + assert known.loc["2:100:A:G", "ENSG2"] == 1 + + +def test_build_known_cis_eqtls_max_snps_per_gene(tmp_path): + parquet = _write_tensorqtl_nominal_parquet( + tmp_path / "nominal.parquet", + [ + ("ENSG1", "1:100:A:G", 0.5), + ("ENSG1", "1:200:A:G", 0.01), + ("ENSG1", "1:300:A:G", 0.02), + ], + ) + known = build_known_cis_eqtls_from_tensorqtl(str(parquet), gene_names=["ENSG1"], max_snps_per_gene=2) + assert known["ENSG1"].sum() == 2 + assert "1:100:A:G" not in known.index # the worst of the 3 is still excluded + assert set(known.index) == {"1:200:A:G", "1:300:A:G"} + + +def test_build_known_cis_eqtls_raises_when_nothing_survives(tmp_path): + parquet = _write_tensorqtl_nominal_parquet( + tmp_path / "nominal.parquet", + [("ENSG1", "1:100:A:G", 0.9)], + ) + with pytest.raises(ValueError, match="No cis-eQTL pairs survived"): + build_known_cis_eqtls_from_tensorqtl(str(parquet), gene_names=["ENSG1"], pval_threshold=0.01) diff --git a/tests/test_tl_external_pc_ld_gsmap.py b/tests/test_tl_external_pc_ld_gsmap.py new file mode 100644 index 0000000..363447d --- /dev/null +++ b/tests/test_tl_external_pc_ld_gsmap.py @@ -0,0 +1,121 @@ +from __future__ import annotations + +import gzip + +import pandas as pd +import pytest + +from cellink._core.donordata import DonorData +from cellink.tl.external._gsmap import load_gsmap_results +from cellink.tl.external._ld import calculate_ld +from cellink.tl.external._pc import calculate_pcs + + +@pytest.fixture +def dd(adata, gdata): + gdata.obs["donor_id"] = gdata.obs.index + return DonorData(G=gdata, C=adata) + + +def test_calculate_pcs_command_construction(dd, tmp_path): + prefix = str(tmp_path / "geno") + cmd = calculate_pcs(dd, prefix, num_pcs=5, run=False) + assert cmd == f"plink --bfile {prefix} --pca 5 --out {prefix}_pca" + + +def test_calculate_pcs_custom_out_prefix(dd, tmp_path): + prefix = str(tmp_path / "geno") + out = str(tmp_path / "custom_pca_out") + cmd = calculate_pcs(dd, prefix, out=out, num_pcs=3, run=False) + assert f"--out {out}" in cmd + assert f"--pca 3" in cmd + + +def test_calculate_pcs_save_cmd_file(dd, tmp_path): + prefix = str(tmp_path / "geno") + cmd_file = tmp_path / "cmd.sh" + result = calculate_pcs(dd, prefix, num_pcs=5, run=False, save_cmd_file=str(cmd_file)) + assert result is None + assert cmd_file.exists() + assert cmd_file.read_text().strip() == f"plink --bfile {prefix} --pca 5 --out {prefix}_pca" + + +def test_calculate_pcs_raises_without_plink_binary(dd, tmp_path, monkeypatch): + monkeypatch.setattr("shutil.which", lambda _: None) + with pytest.raises(ImportError, match="plink is required"): + calculate_pcs(dd, str(tmp_path / "geno"), run=True) + + +def test_calculate_ld_command_construction(dd, tmp_path): + prefix = str(tmp_path / "geno") + cmd = calculate_ld(dd, prefix, window_kb=250, ld_window=500, r2_threshold=0.1, run=False) + assert cmd == ( + f"plink --bfile {prefix} --r2 --ld-window-kb 250 --ld-window 500 " + f"--ld-window-r2 0.1 --out {prefix}_ld" + ) + + +def test_calculate_ld_raises_without_plink_binary_uses_correct_function_name(dd, tmp_path, monkeypatch): + """Regression test: the error message previously named `calculate_pcs` + (copy-paste from _pc.py) even when raised from `calculate_ld`.""" + monkeypatch.setattr("shutil.which", lambda _: None) + with pytest.raises(ImportError, match="calculate_ld"): + calculate_ld(dd, str(tmp_path / "geno"), run=True) + + +def test_load_gsmap_results_missing_workdir_returns_none_fields(tmp_path): + results = load_gsmap_results(tmp_path, sample_name="sample1", trait_name="height") + assert results["spatial_ldsc"] is None + assert results["cauchy_combination"] is None + assert results["report_path"] is None + assert results["workdir"] == tmp_path + + +def test_load_gsmap_results_loads_real_files(tmp_path): + sample_dir = tmp_path / "sample1" + ldsc_dir = sample_dir / "spatial_ldsc" + ldsc_dir.mkdir(parents=True) + ldsc_df = pd.DataFrame({"spot": ["s1", "s2"], "beta": [0.1, 0.2], "se": [0.01, 0.02], "z": [10.0, 10.0], "p": [1e-5, 1e-4]}) + with gzip.open(ldsc_dir / "height_ldsc.csv.gz", "wt") as f: + ldsc_df.to_csv(f, index=False) + + cauchy_dir = sample_dir / "cauchy_combination" + cauchy_dir.mkdir(parents=True) + cauchy_df = pd.DataFrame({"p_cauchy": [0.01, 0.02], "p_median": [0.03, 0.04]}, index=["region1", "region2"]) + cauchy_df.to_csv(cauchy_dir / "height_cauchy.csv") + + report_dir = sample_dir / "report" + report_dir.mkdir(parents=True) + + results = load_gsmap_results(tmp_path, sample_name="sample1", trait_name="height", annotation="domain") + + assert results["spatial_ldsc"] is not None + assert list(results["spatial_ldsc"]["spot"]) == ["s1", "s2"] + assert results["cauchy_combination"] is not None + assert list(results["cauchy_combination"].index) == ["region1", "region2"] + assert results["report_path"] == report_dir + + +def test_load_gsmap_results_falls_back_to_tab_separated(tmp_path): + """Regression test for the comma/tab auto-detection fallback: gsMap + sometimes writes tab-separated .gz files despite the .csv.gz name.""" + sample_dir = tmp_path / "sample1" + ldsc_dir = sample_dir / "spatial_ldsc" + ldsc_dir.mkdir(parents=True) + df = pd.DataFrame({"spot": ["s1", "s2"], "z": [1.0, 2.0]}) + with gzip.open(ldsc_dir / "height_ldsc.csv.gz", "wt") as f: + df.to_csv(f, sep="\t", index=False) + + results = load_gsmap_results(tmp_path, sample_name="sample1", trait_name="height") + assert results["spatial_ldsc"].shape[1] == 2 + assert list(results["spatial_ldsc"]["spot"]) == ["s1", "s2"] + + +def test_load_gsmap_results_no_annotation_skips_cauchy(tmp_path): + sample_dir = tmp_path / "sample1" + cauchy_dir = sample_dir / "cauchy_combination" + cauchy_dir.mkdir(parents=True) + pd.DataFrame({"p_cauchy": [0.01]}).to_csv(cauchy_dir / "height_cauchy.csv") + + results = load_gsmap_results(tmp_path, sample_name="sample1", trait_name="height", annotation=None) + assert results["cauchy_combination"] is None From 9d788d7e0ae30bec90a9b32626d0c0942058868b Mon Sep 17 00:00:00 2001 From: Lucas Arnoldt Date: Tue, 25 Aug 2026 21:57:46 +0100 Subject: [PATCH 2/3] fixes, more tests --- src/cellink/io/_pgen.py | 1 - 1 file changed, 1 deletion(-) diff --git a/src/cellink/io/_pgen.py b/src/cellink/io/_pgen.py index 8aa974f..116a8b6 100644 --- a/src/cellink/io/_pgen.py +++ b/src/cellink/io/_pgen.py @@ -160,7 +160,6 @@ def _read_pvar(pvar_file: Path) -> pd.DataFrame: pv = pv.rename(columns={k: v for k, v in rename_map.items() if k in pv.columns}) if VAnn.index in pv.columns: pv[VAnn.index] = pv[VAnn.index].astype(str) - and if VAnn.chrom in pv.columns: pv[VAnn.chrom] = pv[VAnn.chrom].astype(str).astype("category") From 6e45b362eda39d7dd755720a4f75fc834949abad Mon Sep 17 00:00:00 2001 From: Lucas Arnoldt Date: Tue, 25 Aug 2026 22:13:47 +0100 Subject: [PATCH 3/3] fixes, more tests --- tests/test_seismic_torch.py | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/tests/test_seismic_torch.py b/tests/test_seismic_torch.py index 66ee409..d731a09 100644 --- a/tests/test_seismic_torch.py +++ b/tests/test_seismic_torch.py @@ -1,5 +1,10 @@ from __future__ import annotations +import pytest + +# skip the entire module if torch isn't installed +pytest.importorskip("torch", reason="run_seismic_torch tests need torch, install with `pip install cellink[seismic_torch]`") + import numpy as np import statsmodels.api as sm import torch